Hi everyone,
Am solving the following time dependent nonlinear coupled pdes using the automatic derivative for the jacobian, but after iteration zero, the rest is -nan and am not sure how to fix it. Below is complete code. Thanks to anyone that helps.
from dolfin import *
class PeriodicBoundary(SubDomain):
# left boundary is "target domain"
def inside(self, x, on_boundary):
return bool(x[0]< DOLFIN_EPS and x[0]> - DOLFIN_EPS and on_boundary)
# Map right boundary to left boundary
def map(self, x, y):
y[0] = x[0] - 1.0
y[1] = x[1]
pbc = PeriodicBoundary()
mesh = RectangleMesh(Point(-1,-1),Point(1,1),20,20) #Domain[-1,1]*[-1,1]
V1= VectorFunctionSpace(mesh, 'CG', degree=2, constrained_domain = pbc)
V2=VectorFunctionSpace(mesh, 'CG', degree=2, constrained_domain= pbc)
Q1 = FunctionSpace(mesh, 'CG', degree=1, constrained_domain = pbc )
Q2 = FunctionSpace(mesh, 'CG', degree=1, constrained_domain = pbc)
VQ = MixedFunctionSpace([V1,V2, Q1, Q2])
# Define trial and test functions
(U,u, phi,p ) = TrialFunctions(VQ)
(m,h,q,r) = TestFunctions(VQ)
H= Function(VQ)
# Dirichlet boundary condition
# Sub domain for Dirichlet boundary condition
class Top(SubDomain):
def inside(self, x, on_boundary):
return (x[1]>1 - DOLFIN_EPS and on_boundary)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return (x[1]< -1 + DOLFIN_EPS and on_boundary)
UB= Constant((-1,0))
UT= Constant((1,0))
ub= Constant((-1,0))
ut= Constant((1,0))
top= Top()
bott = Bottom()
bcUB= DirichletBC(VQ.sub(0), UB ,bott)
bcUT= DirichletBC(VQ.sub(0), UB ,top)
bcub= DirichletBC(VQ.sub(1), ub ,bott )
bcut= DirichletBC(VQ.sub(1), ub ,top )
bcs = [bcUB,bcUT, bcub,bcut]
#constants
c = Constant (2/3)
xi0= Constant(1)
n = Constant(2)
alpha = Constant(25)
# initial conditions
U_int = Expression(["x[1]", "x[1]"]) # U=u=y
u_int = Expression(["x[1]", "x[1]"])
phi_int= Constant(0.01)
p_int = Constant(0.0)
U00 = interpolate(U_int, V1) # REPLACE H BY V1. SAME FOR THE REST
u00= interpolate(u_int,V2)
phi00=interpolate(phi_int,Q1)
p00 = interpolate(p_int, Q2)
U0 = U00
u0 = u00
phi0 = phi00
p0 = p00
k_phi = (phi/phi0)**n
eta = exp(-alpha*(phi-phi0))
xi = xi0* eta
dt=0.01
f1 = phi*q*dx - phi0*q*dx - dt*div((1-phi)*U)*q*dx #=0
f2 = inner(phi*(u-U), m)*dx + k_phi * p*div(m)*dx # =0
f3 = p*div(h)*dx + eta*inner(grad(U) + grad(U).T, grad(h))*dx - (xi*div(U)) *div(h)*dx + c*(eta*div(U))*div(h)*dx # =0
f4 = div(phi*u + (1-phi)*U)*r*dx #=0
F = f1 + f2 + f3 + f4
T_total = 1
t=dt
while (t<=T_total):
F1= action(F,H)
# Automatic computation of the Jacobian using Gateaux derivative#
Jac = derivative(F1,H)
# Automatic Defining the nonlinear problem with boundary conditions
problem = NonlinearVariationalProblem(F1,H,bcs, Jac)
# Creating automatic solver
solver = NonlinearVariationalSolver(problem)
solver.solve()
#ALL_initial(0,0).assign(H)
(US, uS, phiS, ps) = split(H)
assigner = FunctionAssigner(Q1, VQ.sub(2))
assigner.assign(phi0, phiS)
t+=dt
The error is as below
Newton iteration 0: r (abs) = 1.265e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 2: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
................................................................
............................................................
Newton iteration 50: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Traceback (most recent call last):
File "mat.py", line 125, in <module>
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics@fenicsproject.org
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** DOLFIN version: 1.6.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
