Coupled non linear eqautions on multiple subdomains


74
views
0
12 weeks ago by
Moritz  
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)

F_am = sigma*inner(grad(u[0]), grad(v[0]))*(dx(8)+dx(9)+dx(10)+dx(11))- i_BV(u)*v[0]*(dx(10)+dx(11))

F_electrolyte = kappa*inner(grad(u[1]), grad(v[1]))*(dx(10)+dx(11)+dx(12))+ i_BV(u)*v[1]*(dx(10)+dx(11))

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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »