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.1751x[0]*x[0]+0.1741)’, degree =2 ), facets, 4)] # Emperical BC applied at the upper and lower edges.
=============================================================================
Define source terms
f_1 = project(-kaDH.dx(0)+2kaFH.dx(0)-kaGH.dx(1),q)
f_2 = 0.5E2(D-1)(x1.dx(0)f+x2.dx(0)g+x3.dx(0)h)
f_3 = project(0.5E2(D-1)(x1.dx(1)f+x2.dx(1)g+x3.dx(1)h),q)
f_4 = project((Db11-2Fb12+Gb22)pow(a,-1.5),q)
f_5 = project((DG-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((b22b11-b12b12)*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
FF = -(aDv_1.dx(0)H.dx(0)+v_1a.dx(0)DH.dx(0)+v_1aD.dx(0)H.dx(0))kdx
+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_1D(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_12F(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_1G(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(4kpow(a,2)pow(H,3)+2(k1-k)Hpow(a,2)K)dx
-2v_1pow(a,2)H(kHH+k1K+E1/8(G-1)(G-1)+E2/8(D-1)(D-1))dx
-0.5v_1E1*(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.5v_1E2(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_1ppow(a,2)dx
-v_2(E10.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(E10.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,