Problem implementing Schwarz Method for Laplace's equation


187
views
0
8 months ago by
CRB  
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.

We had to enable extrapolation.

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.

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

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

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
radius = 3*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 = project(grad(solInside), WInside)
EInside.set_allow_extrapolation(True)
# EOutside.interpolate(EInside)
# plot(project(EInsideProjected, WOutside))

aOutside = inner(grad(uOutside), grad(vOutside))*dx
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
aInside = inner(grad(uInside), grad(vInside))*dx
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()

EOutside = project(grad(solOutside), WOutside)

plot(solOutside, mode="color")
plot(solInside, mode="color")
plot(EOutside)
plot(EInside, interactive=True)
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »