### Problem implementing Schwarz Method for Laplace's equation

187
views
0
8 months ago by
We are trying to couple two domains (on which we solve Laplace's equation) using a Schwarz domain decomposition.

In our Problem setting one (inside) domain will be embedded in another (outside) domain (imagine a biological cell in a medium).

The mesh:
https://bitbucket.org/CRBFE/fexp/src/11410f1f716f80965a020a26e0ff08a97ba2a6d8/screen/mesh.png?at=WORK&fileviewer=file-view-default

To get a feeling for the method we started with a very small example, a square domain (10 um) with Dirichlet b.c. on the top and bottom (+5e-3 and -5e-3) and homogeneous Neumann b.c. on both sides.

To couple the inner (circular) domain to the outer domain we used identity to couple u and gradient(u) in a Dirichlet-Neumann scheme. (later this will be a function of the gradient and the drop in u across the inner interface)

We would have prefered to implement a multiplicative Schwarz scheme, but we were unable to compose a single bilinear form with functions "living" on two separates subspaces (or sub-meshes).

The next problem while implementing the Schwarz method was, that, for some reason, the points on the inner boundary between both domains are not considered parts of the domains:

*** Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.

Unfortunately the method we get does not converge.

On the square domain, with the b.c. described above one should get a homogeous gradient (like the field between two plates of a plate capacitor).

The way it diverges from the correct solution is quite particular, the field in the inner domain seems to turn a little bit counter-clockwise in each iteration.

After about 200 iterations it goes to a steady state, where the field in the inner domain is off bei 90° and far away from homogeneous.

The code is here:

https://bitbucket.org/CRBFE/fexp/src/11410f1f716f?at=WORK

https://bitbucket.org/CRBFE/fexp/src/11410f1f716f80965a020a26e0ff08a97ba2a6d8/test_contact_impedance-2mesh.py?at=WORK&fileviewer=file-view-default

from dolfin import *
import numpy as np

# PARAMETERS
width = 10*1E-6
height = 10*1E-6

epsrCyt = 50 # [1] relative permitivity cytoplasm
epsrMed = 80 # [1] membrane
epsrMem = 9.0 # [1] medium

sigmaCyt = 0.53 # [S/m] cytoplasm
sigmaMem = 5*1E-7 # [S/m] membrane
sigmaMed = 1.2 # [S/m] medium

memThickness = 5*1E-9 # [m] membrane thickness

VUpper = 0.5*1000*height # E = 1V/mm = 1000V/m
VLower = -0.5*1000*height #

# MESHES
mesh = Mesh('./geometry/membrane-potential.xml')
subdomains = MeshFunction("size_t", mesh, "./geometry/membrane-potential_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "./geometry/membrane-potential_facet_region.xml")

meshOutside = SubMesh(mesh, subdomains, 0)
meshInside = SubMesh(mesh, subdomains, 1)

# plot(meshInside, "Mesh", interactive=True)

# FUNCTION SPACES
VOutside = FunctionSpace(meshOutside, "CG", 1) # medium == external domain
VInside = FunctionSpace(meshInside, "CG", 1) # cell == internal domain

W = VectorFunctionSpace(mesh, "DG", 0)
WOutside = VectorFunctionSpace(meshOutside, "DG", 0)
WInside = VectorFunctionSpace(meshInside, "DG", 0)

# Define test and trial functions
uOutside = TrialFunction(VOutside)
vOutside = TestFunction(VOutside)

uInside = TrialFunction(VInside)
vInside = TestFunction(VInside)

# boundaries
# left = CompiledSubDomain("near(x[0], side) && on_boundary", side = -0.5*width)
# right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.5*width)
circOutside = CompiledSubDomain("(x[0] < 0.48*width && x[0] > -0.48*width && x[1] < 0.48*height && x[1] > -0.48*height) && on_boundary", width=width, height=height)
circInside = CompiledSubDomain("on_boundary")

# WEAK FORMS
fOutside = Constant(0)
fInside = Constant(0)

# for connecting inside and outside domain
inbound = FacetFunction('size_t', meshOutside, 0)
inbound.set_all(0)
circOutside.mark(inbound, 1)
# plot(inbound, "Mesh", interactive=True)
dsInsideOutside = Measure('ds', domain=meshOutside, subdomain_data=inbound)
normalOutside = FacetNormal(meshOutside)

upper = CompiledSubDomain("near(x[1], side) && on_boundary", side = 0.5*height)
lower = CompiledSubDomain("near(x[1], side) && on_boundary", side = -0.5*height)
# upper.mark(inbound, 2)
# lower.mark(inbound, 3)
# plot(inbound, "Mesh", interactive=True)

EOutside = Function(WOutside) # electric field on and outside of membrane
solInside = Function(VInside) # potential in the cell
solOutside = Function(VOutside) # potential in the medium
u2OnBoundary = Function(VInside)

# solInside.interpolate(Constant(0))

set_log_level(WARNING)
iteration = 0
while iteration < 50:
iteration += 1
print "Iter " + str(iteration) + ":"

# solve medium domain
EInside.set_allow_extrapolation(True)
# EOutside.interpolate(EInside)
# plot(project(EInsideProjected, WOutside))

LOutside = fOutside*vOutside*dx + dot(EInside, normalOutside)*vOutside*dsInsideOutside(1)
bc1a = DirichletBC(VOutside, Constant(VUpper), upper)
bc1b = DirichletBC(VOutside, Constant(VLower), lower)
problemOutside = LinearVariationalProblem(aOutside, LOutside, solOutside, [bc1a, bc1b])
solverOutside = LinearVariationalSolver(problemOutside)

prm = solverOutside.parameters.krylov_solver # short form
prm.absolute_tolerance = 1E-12
prm.relative_tolerance = 1E-7

solverOutside.solve()

# solve cell domain
LInside = fInside*vInside*dx

solOutside.set_allow_extrapolation(True)

# u2OnBoundary.interpolate(solOutside)
bc2a = DirichletBC(VInside, solOutside, circInside)
problemInside = LinearVariationalProblem(aInside, LInside, solInside, [bc2a])
solverInside = LinearVariationalSolver(problemInside)
prm = solverInside.parameters.krylov_solver # short form
prm.absolute_tolerance = 1E-16
prm.relative_tolerance = 1E-7

solverInside.solve()