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

187

views

0

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)

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.