### Coupled non linear eqautions on multiple subdomains

74
views
0
12 weeks ago by
Hello,
I am tying to solve a set of cpupled equations on a mesh with several subdomains. The states are only defined on some subdomain (the mesh is imported from gmsh and I will use its indecis).
The equations are:
$-\nabla (\sigma \nabla \Phi_1) - i_{bv}(\Phi_1-\Phi_2)=0 \text{ on subdomains 8,9,10,11} \\ -\nabla (\kappa \nabla \Phi_2 )+ i_{bv}(\Phi_1-\Phi_2)=0 \text{ on subdomains 10,11,12} \\ i_{bv}(\Delta\phi)=i_0(exp(\alpha \Delta\phi ) -exp(-(1-\alpha)\Delta\phi))$
The subdomain pairs   (8,10) and (9,11) are seperated by subdomain 12.
I can only impose Dirichlet COnditions for $\Phi_1$ on facets of 8 and 9.
The Newton solver exits after 0 iterations:
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-10)
Newton solver finished in 0 iterations and 0 linear solver iterations.
​

Edit(fixed a typo, new output) :

Newton iteration 0: r (abs) = 4.466e-12 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-10)
Newton solver finished in 0 iterations and 0 linear solver iterations.

The "result" is all zero, so it seems to ignore the boundary condition.

Since I am quite new to Fenics, I am not sure where to start debugging.
My guess is, that I have to much DoF since the states are defined on the whole mesh. But I didnt find a way to remove them.
Here is my code:

from fenics import *
mesh=Mesh("../test_mesh/test_mesh.xml")
cell_domains=MeshFunction('size_t',mesh,"../test_mesh/test_mesh_physical_region.xml")
face_domain=MeshFunction('size_t',mesh,"../test_mesh/test_mesh_facet_region.xml")
vtkfile_fd = File('fd.pvd')
vtkfile_fd<<face_domain;
vtkfile_cd = File('cd.pvd')
vtkfile_cd<<cell_domains;
#dx = Measure('dx', domain=mesh, subdomain_data=cell_domains)
cell = tetrahedron

scalar = FiniteElement("DG", cell, 1)

TH = MixedElement([scalar,scalar])
V=FunctionSpace(mesh,TH)
u=Function(V);

v=TestFunction(V)

u_n=Coefficient(V)
i0=Coefficient(scalar);
alpha=Coefficient(scalar);
sigma=Coefficient(scalar);
kappa=Coefficient(scalar);

def i_BV(x):
return i0*(exp(alpha*(x[0]-x[1]-0.5))-exp(-(1.-alpha)*(x[0]-x[1]-0.5)))

class ExchangeCurrentDensity(Expression):
def __init__(self,domains, j0_pos, j0_neg, **kwargs):
self.domains=domains
self.j0_pos=j0_pos
self.j0_neg=j0_neg

def eval_cell(self,value,x,cell):
if self.domains[cell.index]==10:
value[0]=self.j0_pos
elif self.domains[cell.index]==11:
value[0]=self.j0_neg
else:
value[0]=0.

class Sigma(Expression):
def __init__(self,domains, sigma_pos, sigma_neg, **kwargs):
self.domains=domains
self.sigma_pos=sigma_pos
self.sigma_neg=sigma_neg

def eval_cell(self,value,x,cell):
if (self.domains[cell.index]==8):
value[0]=self.sigma_pos
elif (self.domains[cell.index]==9):
value[0]=self.sigma_neg
elif (self.domains[cell.index]==10):
value[0]=self.sigma_pos
elif (self.domains[cell.index]==11):
value[0]=self.sigma_neg
elif (self.domains[cell.index]==12):
value[0]=0
else:
value[0]=0.

class Kappa(Expression):
def __init__(self,domains,sigma_electrolyte, **kwargs):
self.domains=domains
self.sigma_electrolyte=sigma_electrolyte

def eval_cell(self,value,x,cell):
if (self.domains[cell.index]==10):
value[0]=self.sigma_electrolyte
elif (self.domains[cell.index]==11):
value[0]=self.sigma_electrolyte
elif (self.domains[cell.index]==12):
value[0]=self.sigma_electrolyte
else:
value[0]=0.

i0=ExchangeCurrentDensity(dx,1e-3,-1e-3,degree=1)
alpha=Constant(0.5)
sigma=Sigma(cell_domains,1e5,1e6,degree=0);
kappa=Kappa(cell_domains,1e2,degree=0);

bc_neg=DirichletBC(V.sub(0), Constant(0.),face_domain,1);
bc_pos=DirichletBC(V.sub(0), Constant(1e3),face_domain,2);
bcs=[bc_pos,bc_neg]

u=interpolate(Constant((1.,0.)),V)

solve((F_am+F_electrolyte) == 0, u,bcs,
solver_parameters={"newton_solver":{"relative_tolerance":1e-10}})

vtkfile_u=File('u.pvd')
vtkfile_u<<(interpolate(u,V))

​

I got it to "iterate" by using

u=interpolate(init,V)
problem=NonlinearVariationalProblem(F,u,bcs,J)
solver=NonlinearVariationalSolver(problem)


and defining non zero initial conditions.
But the iterations still does nothing, having the same error and exiting with maximum iterations.

Community: FEniCS Project