### Coupled System with disjoint domains

114
views
0
12 weeks ago by
This is an alternate appraoch to my problem here: https://www.allanswered.com/post/zxkvx/coupled-non-linear-eqautions-on-multiple-subdomains/
Here I am trying to use submeshes to restrict the functionspaces to the relevant domains and the solve them on each SubMesh. But I dont get to solve even one.
The output is:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.383e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
File "ChargeConservationSubMesh.py", line 166, 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
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
*** Where:   This error was encountered inside /build/dolfin-C0sq7Z/dolfin-2017.2.0.post0/dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------
​

Here my adapted code and msh file:
File attached: test_mesh.msh (3.45 MB)
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")

scalar = FiniteElement("CG",mesh.ufl_cell(), 1)

positive_electrode = MeshFunction("size_t", mesh,mesh.topology().dim())
positive_electrode.array()[cell_domains.array()==8]=1
positive_electrode.array()[cell_domains.array()==10]=1

negative_electrode = MeshFunction("size_t", mesh,mesh.topology().dim())
negative_electrode.array()[cell_domains.array()==9]=1
negative_electrode.array()[cell_domains.array()==11]=1

electrolyte_electrode = MeshFunction("size_t", mesh,mesh.topology().dim())
electrolyte_electrode.array()[cell_domains.array()==10]=1
electrolyte_electrode.array()[cell_domains.array()==11]=1
electrolyte_electrode.array()[cell_domains.array()==12]=1

mesh_positive=SubMesh(mesh,positive_electrode,1)
mesh_negative=SubMesh(mesh,negative_electrode,1)
mesh_electrolyte=SubMesh(mesh,electrolyte_electrode,1)

faces_pos=MeshFunction('size_t',mesh_positive,"../test_mesh/test_mesh_facet_region.xml")
faces_neg=MeshFunction('size_t',mesh_negative,"../test_mesh/test_mesh_facet_region.xml")
dx_pos = Measure('dx', domain=mesh_positive)
ds_pos = Measure('ds', domain=mesh_positive)
dx_neg = Measure('dx', domain=mesh_negative)
dx_elec = Measure('dx', domain=mesh_electrolyte)

V_positive=FunctionSpace(mesh_positive,scalar)
V_negative=FunctionSpace(mesh_negative,scalar)
V_electrolyte=FunctionSpace(mesh_electrolyte,scalar)

u_positive=Function(V_positive);
v_positive=TestFunction(V_positive)

u_negative=Function(V_negative);
v_negative=TestFunction(V_negative)

u_electrolyte=Function(V_electrolyte);
v_electrolyte=TestFunction(V_electrolyte)

0=Coefficient(scalar);
u0=Coefficient(scalar);
alpha=Coefficient(scalar);
sigma_positive=Coefficient(scalar);
kappa=Coefficient(V_electrolyte);

def i_BV(x_electrode,x_electrolyte):
return i0*(exp(alpha*(x_electrode-x_electrolyte-u0))-exp(-(1.-alpha)*(x_electrode-x_electrolyte-u0)))

class EquilibirumPotential(Expression):
def __init__(self,domains, u0_pos, u0_neg, **kwargs):
self.domains=domains
self.u0_pos=u0_pos
self.u0_neg=u0_neg

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

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(cell_domains,1e0,1e0,degree=0)
u0=EquilibirumPotential(cell_domains,1.75,-0.25,degree=0)
alpha=Constant(0.5)
sigma=Sigma(cell_domains,1e6,1e6,degree=0);
kappa=Kappa(cell_domains,1e2,degree=0);

bc_pos=DirichletBC(V_positive, Constant(2),faces_pos,2);
bc_neg=DirichletBC(V_negative, Constant(0),faces_neg,1);
bcs=[bc_pos]

u_positive=interpolate(Constant((1.9)),V_positive)
u_negative=interpolate(Constant(0.0),V_negative)
u_electrolyte=interpolate(Constant(0.25),V_electrolyte)
J=derivative(F_pos,u_positive)
problem=NonlinearVariationalProblem(F_pos , u_positive,bcs,J)
solver=NonlinearVariationalSolver(problem)
prm=solver.parameters

solver.solve()

​
Community: FEniCS Project
1
You might try looking into multiphenics:

https://github.com/mathLab/multiphenics

I haven't used multiphenics myself, but taking a closer look at it has been on my to-do list for a while.
written 12 weeks ago by David Kamensky
Hi,
I got it to work with MultiPhenics for my test problem. I have to see if it works with a more complex setup.
(FYI Tutorial 3 helped me here)
written 12 weeks ago by Moritz

1
12 weeks ago by
Hi Moritz Huck,

you may like the implementation under:
https://github.com/afqueiruga/EMSI-2018

Have a look at the python files in the src directory. SubMesh is working good in one cpu. If you want to use parallel computing, cbc_post submesh is useful, or you can also use the submesh function as used in the code published in github.

I believe that the problem in your implementation is mainly the use of different measures, namely dx_pos and dx_neg. I use the same measure and marking to integrate only on a part of the domain.

Best, Emek
Thank you for your sources, I will look through them.

I dont understand what you mean by different measures. As far as I can see you are also using different measures for the submesh and the complete mesh (dV and dv).
If I use a dx defined over the whole mesh I get a bunch of compilation errors.
written 12 weeks ago by Moritz
You are right, you use different submeshes for each measure. Another look into your code, I see that you use dx_pos(10) without giving the list into the measure. I believe you need to do something like that

cells_pos = MeshFunction('size_t', mesh_positive, 2/3, 0) #depending on 2D or 3D, default to set all 0
# here a loop to match with other list and write into cells_pos the marking number