# Fenics Nan for solving PDEs systems

Hi, Fenics staff and users, I had problem on solving a highly difficult PDEs system, the results are Nan. In the PDEs system, the core unknowns are x1 x2 x3 with first three equations (total 21 unknows ), while inside the PDEs system there are parameters(that is why it has 21 unknowns) expressed by x1, x2 x3, which make this PDEs system complicated. I used mixed method and weak form to solve such problem. Could any suggestions pointed out for this simulation ? Thanks for your attention !
The code are shown below
L = 1
H = 1
nx = 20
ny = 20
mesh = RectangleMesh(Point(-1, -1), Point(L, H), nx, ny)
sag = mesh.coordinates()

# Initializing material parameters

k = Constant(0.82)
k1 = Constant(-0.82)
E1 = Constant(0.2)
E2 = Constant(0.2)
p = Constant(0.5)

# Defining test and trial functions in splitted form

element = VectorElement(‘P’, triangle, 1, dim=21)
v = FunctionSpace(mesh, element)
q = FunctionSpace(mesh, “DG”, 0)

plot(mesh)
plt.show()

v_1, v_2, v_3, v_4, v_5, v_6, v_7, v_8, v_9, v_10, v_11, v_12, v_13, v_14, v_15, v_16, v_17, v_18, v_25, v_26, v_27 = TestFunctions(v)

u = Function(v)

x1, x2, x3, H, a, D, F, K, G, Y, Z, b, c, d, e, f, g, h, b11, b22, b12 = split(u)

x, y = SpatialCoordinate(mesh)

# Defining boundary conditions

class left(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return abs(x[0]+1) < tol

class right(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return abs(x[0]-1.0) < tol

class top(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return abs(x[1]-1.0) < tol

class bottom(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return abs(x[1]+1.0) < tol

class middleX(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return abs(x[0]) < tol

class middleY(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return abs(x[1]) < tol

facets = MeshFunction(“size_t”, mesh, 1)
facets.set_all(0)
left().mark(facets, 1)
right().mark(facets, 2)
top().mark(facets, 3)
bottom().mark(facets, 4)
middleX().mark(facets, 5)
middleY().mark(facets, 6)
ds = Measure(“ds”, subdomain_data=facets)

def boundary(x, on_boundary):
return on_boundary

## Boundary Conditions

u_R = Expression(‘x[1]’, degree=1)
u_R1 = Expression(‘x[0]’, degree=1)
bc = [
DirichletBC(v.sub(2),Constant(0),facets,1),
DirichletBC(v.sub(2),Constant(0),facets,2),
DirichletBC(v.sub(2),Constant(0),facets,3),
DirichletBC(v.sub(2),Constant(0),facets,4),
DirichletBC(v.sub(0),u_R1,facets,1),
DirichletBC(v.sub(0),u_R1,facets,2),
DirichletBC(v.sub(1),u_R,facets,3),
DirichletBC(v.sub(1),u_R,facets,4),
DirichletBC(v.sub(0),u_R1, facets, 5),
DirichletBC(v.sub(1),u_R, facets, 6),
##
DirichletBC(v.sub(3),Constant(0),facets,1),
DirichletBC(v.sub(3),Constant(0),facets,2),
DirichletBC(v.sub(3),Constant(0),facets,3),
DirichletBC(v.sub(3),Constant(0),facets,4),
##
]
# DirichletBC(V.sub(1), Expression(‘(1-(-0.1751*x[0]x[0]+0.1741))‘, degree =2 ), facets, 3), # Emperical BC applied at the upper and lower edges.
# DirichletBC(V.sub(1), Expression(’(-0.1751
x[0]*x[0]+0.1741)’, degree =2 ), facets, 4)] # Emperical BC applied at the upper and lower edges.

# Define source terms

f_2 = 0.5E2(D-1)(x1.dx(0)f+x2.dx(0)g+x3.dx(0)h)
f_3 = project(0.5
E2
(D-1)
(x1.dx(1)f+x2.dx(1)g+x3.dx(1)h),q)
f_4 = project((D
b11-2
F
b12+G
b22)pow(a,-1.5),q)
f_5 = project((D
G-F*F),q)
f_6 = project((x1.dx(1)*x1.dx(1)+x2.dx(1)*x2.dx(1)+x3.dx(1)*x3.dx(1)),q)
f_7 = project((x1.dx(0)*x1.dx(1)+x2.dx(0)x2.dx(1)+x3.dx(0)x3.dx(1)),q)
f_8 = project((b22
b11-b12
b12)*pow(a,-2),q)
f_9 = project((x1.dx(0)*x1.dx(0)+x2.dx(0)*x2.dx(0)+x3.dx(0)x3.dx(0)),q)
f_10 = project(x1.dx(0),q)
f_11 = project(x2.dx(0),q)
f_12 = project(x3.dx(0),q)
f_13 = project(x1.dx(0),q)
f_14 = project(x2.dx(0),q)
f_15 = project(x3.dx(0),q)
f_16 = project(x1.dx(1),q)
f_17 = project(x2.dx(1),q)
f_18 = project(x3.dx(1),q)
f_25 = project((Y
(x2.dx(0)*x3.dx(1)-x2.dx(1)x3.dx(0))+Z(x1.dx(1)*x3.dx(0)-x1.dx(0)x3.dx(1))+b(x1.dx(0)*x2.dx(1)-x2.dx(0)x1.dx(1))),q)
f_26 = project((f
(x2.dx(0)*x3.dx(1)-x2.dx(1)x3.dx(0))+g(x1.dx(1)*x3.dx(0)-x1.dx(0)x3.dx(1))+h(x1.dx(0)*x2.dx(1)-x2.dx(0)x1.dx(1))),q)
f_27 = project((c
(x2.dx(0)*x3.dx(1)-x2.dx(1)x3.dx(0))+d(x1.dx(1)*x3.dx(0)-x1.dx(0)x3.dx(1))+e(x1.dx(0)*x2.dx(1)-x2.dx(0)*x1.dx(1))),q)

# Define variational problem

+2
(v_1.dx(1)aFH.dx(0)+v_1a.dx(1)FH.dx(0)+v_1aF.dx(1)H.dx(0))kdx-(v_1.dx(1)aGH.dx(1)+v_1a.dx(1)GH.dx(1)+v_1aG.dx(1)H.dx(1))kdx
-v_1
D
(H.dx(0)(Y(x1.dx(0)D-x1.dx(1)F)+Z(x2.dx(0)D-x2.dx(1)F)+b(x3.dx(0)D-x3.dx(1)F))+H.dx(1)(Y(-x1.dx(0)F+x1.dx(1)G)+Z(-x2.dx(0)F+x2.dx(1)G)+b(-x3.dx(0)F+x3.dx(1)G)))kdx
+v_1
2
F
(H.dx(0)
(c
(x1.dx(0)D-x1.dx(1)F)+d(x2.dx(0)D-x2.dx(1)F)+e(x3.dx(0)D-x3.dx(1)F))+H.dx(1)(c(-x1.dx(0)F+x1.dx(1)G)+d(-x2.dx(0)F+x2.dx(1)G)+e(-x3.dx(0)F+x3.dx(1)G)))kdx
-v_1
G
(H.dx(0)
(f
(x1.dx(0)D-x1.dx(1)F)+g(x2.dx(0)D-x2.dx(1)F)+h(x3.dx(0)D-x3.dx(1)F))+H.dx(1)(f(-x1.dx(0)F+x1.dx(1)G)+g(-x2.dx(0)F+x2.dx(1)G)+h(-x3.dx(0)F+x3.dx(1)G)))kdx
+v_1
(4
k
pow(a,2)pow(H,3)+2(k1-k)Hpow(a,2)K)dx
-2
v_1
pow(a,2)H(k
H
H+k1
K+E1/8
(G-1)
(G-1)+E2/8
(D-1)(D-1))dx
-0.5
v_1
E1*(G-1)(Y(x2.dx(0)x3.dx(1)-x2.dx(1)x3.dx(0))+Z(x1.dx(1)x3.dx(0)-x1.dx(0)x3.dx(1))+b(x1.dx(0)x2.dx(1)-x2.dx(0)x1.dx(1)))pow(a,1.5)dx
-0.5
v_1
E2
(D-1)
(f
(x2.dx(0)x3.dx(1)-x2.dx(1)x3.dx(0))+g(x1.dx(1)x3.dx(0)-x1.dx(0)x3.dx(1))+h(x1.dx(0)x2.dx(1)-x2.dx(0)x1.dx(1)))pow(a,1.5)dx
-v_1
p
pow(a,2)dx
-v_2
(E1
0.5
(G-1)
(x1.dx(0)Y+x2.dx(0)Z+x3.dx(0)b)+E20.5(D-1)(x1.dx(1)c+x2.dx(1)d+x3.dx(1)e)+0.5E1G.dx(0)G+0.5E1(G-1)
(x1.dx(0)Y+x2.dx(0)Z+x3.dx(0)b)+0.5E2D.dx(1)F)dx
-v_3
(E1
0.5
(G-1)
(x1.dx(0)c+x2.dx(0)d+x3.dx(0)e)+E20.5(D-1)(x1.dx(1)f+x2.dx(1)g+x3.dx(1)h)+0.5E1G.dx(0)F+0.5E1(G-1)(x1.dx(1)Y+x2.dx(1)Z+x3.dx(1)b)+0.5E2D.dx(1)D)dx+(v_4H)dx+(v_5a)dx+(v_6D)dx+(v_7F)dx+(v_8K)dx+(v_9G)dx+(v_10Y+v_10.dx(0)x1.dx(0))dx+(v_11Z+v_11.dx(0)x2.dx(0))dx+(v_12b+v_12.dx(0)x3.dx(0))dx+(v_13c+v_13.dx(1)x1.dx(0))dx+(v_14d+v_14.dx(1)x2.dx(0))dx+(v_15e+v_15.dx(1)x3.dx(0))dx+(v_16f+v_16.dx(1)x1.dx(1))dx+(v_17g+v_17.dx(1)x2.dx(1))dx+(v_18h+v_18.dx(1)x3.dx(1))dx+(v_25b11)dx+(v_26b22)dx+(v_27b12)dx-f_01v_1ds(1)-f_1v_1ds- f_2v_2dx -f_3v_3dx- f_4v_4dx - f_5v_5dx - f_6v_6dx-f_7v_7dx -f_8v_8dx -f_9v_9dx -f_10v_10ds-f_11v_11ds-f_12v_12ds-f_13v_13ds-f_14v_14ds-f_15v_15ds-f_16v_16ds-f_17v_17ds-f_18v_18ds-f_25v_25dx-f_26v_26dx-f_27v_27*dx

##PICARD Iteration
u_ = Expression((‘0.01’,‘0.01’,‘0.01’,‘0.01’,‘1.0’,‘1.0’,‘0.01’,‘0.01’,‘1.0’,‘0.01’,‘0.01’,‘0.01’,‘0.01’,‘0.01’,‘0.01’,‘0.01’, ‘0.01’,‘0.01’,‘1’,‘1’,‘1’), degree =1)
u_k= interpolate(u_,v)
eps = 1.0
tol = 1.0E-2
iter = 0
maxiter = 1

while eps > tol and iter < maxiter:
iter +=1
solve(FF == 0, u, bc)
diff = u.vector().get_local() - u_k.vector().get_local()
eps = np.linalg.norm(diff, ord=np.Inf)
print ('iter= ', iter )
print ('eps= ', eps)
u_k.assign(u)

Best,

Your code is too long for me to parse, but consider the tips and ticks here: Default absolute tolerance and relative tolerance - #4 by nate

Also consider walking before you can run. By this I mean, start with a small simple system making sure that works; and, gradually add complexity until you reach your final full model.

Thanks for your attention and review on my problem !