# Coupling the equations defined on two different subdomains

I have a question about solving problem defined on two subdomains.
Namely, I have a system of equations consisting of a Poisson equation and several drift-diffusion equations (number varies depending on a model, but its always more than 4). Poisson equation is defined on the whole domain \Omega:
\varepsilon \nabla^2 \varphi = -\rho(n_i) in\:\: \Omega,
where \varepsilon= \begin{cases} 1, & in \:\: \Omega_1\\ c, & \text{otherwise} \end{cases}
while drift-diffusion equations are defined only on subdomain \Omega_1:
\frac{\partial n_i}{\partial t} + \nabla \cdot \mathbf{\Gamma_{\textit{i}}}(\varphi, n_i)=S_{\textit{i}}(\varphi, n_i).
For Poisson equation, Dirichlet boundary condition, along with Robin boundary condition
-\mathbf{n} \cdot \mathbf{\nabla \varphi} = \sigma(\varphi, n_i) at the interface between two subdomains, is used. For drift-diffusion equations, Robin boundary condition \mathbf{n} \cdot \mathbf{\Gamma} = f(\varphi, n_i) is used at all the boundaries. The problem is illustrated in a figure below.

I tried different approaches to solve this problem, like solving Poisson equation on one and drift-diffusion equations on the other mesh and then transferring the values between functions defined on different meshes by using transfer_matrix. This approach works well, however it requires pretty small time step. By using a commercial software, when all the equations are coupled, time step requirement is much less restrictive.

As the other option, I tried to use mixed-dimensional branch (ceciledc/fenics_mixed_dimensional:latest), however, I had a problem since it reports that mshr module is missing (and it is required for my use due to geometry). I tried to install it, but it still reports that mshr is missing.

Other tries to solve this problem were not successful. So my question is, how to couple these equations, if they are defined on two different subdomains and is this even possible with FEniCS?

Hello,

Regarding the use of the mixed-dimensional branch, I would suggest building your mesh using an external software such as Gmsh. This post gives more details on how to convert a Gmsh mesh using meshio and how to import it into FEniCS using xdmf format.

Hi,

Thank you for a suggestion. I already use gmsh, but for small tests I usually use mshr. However, if it is not available, it won’t be a problem.

Concerning the other question, is it possible to solve the problem, where one equation is defined over a whole domain, while the others are defined only on the part of it, using the mixed-dimensional branch. Also, is there any example available for a nonlinear problem defined on submeshes using mixed-dimensional branch?

Yes, it’s possible to create submeshes corresponding to the part of the domain your equations belong to, and to couple them in the mixed formulation.
As a (very basic) non-linear example, you can take a look to this demo

1 Like

Thank you for the example, I checked it out and it seems that it can be easily adapted to my problem. I have couple more questions concerning the surface integrals and custom solvers in a mixed-dimensional branch.

Concerning solvers, is it possible to use a custom made nonlinear solvers, such as one given here and is it possible to use PETSc SNES. I noticed a comment in the demo file stating that SNES is not available yet.

How is marking of the boundaries done in the mixed-dimensional branch, is it the same as in the old version, by using MeshFunction()?

Concerning solvers, there might still be some issues with PETSc SNES. I recommend using a NewtonSolver instead. Also you could try to define a customized NewtonSolver as in the post you mention (I haven’t tried doing that).

Regarding the marking of the boundaries, you can use MeshFunction as usual since the submeshes you build using MeshView are standard Mesh objects (with some additional mappings).

Thank you very much for your help. I think that mixed-dimensional branch solves my problem.

I am sorry for bumping the thread, however I have a problem with a solver. Namely, after implementing changes (required by mixed-dimensional branch) to my original code, the program crashes with an error stating: Unable to successfully call PETSc function 'KSPSolve. Since the original code relies on external modules, I have tried to build simplest possible MWE that can illustrate the problem with the solver:

from dolfin import *

r = Expression('x[0]', degree=1) # coordinate x redefinition to r required for cylindrical coordinates

parameters['allow_extrapolation'] = True

elementary_charge = 1.6021766208e-19 #[C]
epsilon_0 = 8.854187817e-12 #[F/m]
def mu(x):
return 2.3987*x**(-0.26)

def Diffusion_coefficient(x):
return 4.3628e-3*x**(0.22)

def alpha(x):
return (1.1944e6+4.3666e26 * x**(-3))*exp(-2.73e7/x)-340.75

dt = Expression("time_step", time_step = 5e-12, degree = 0) #time step [s]
t = 0.0

mesh = RectangleMesh(Point(0.0, 0.0),Point(1.25e-2, 1.35e-2), 15, 15, "crossed")

normal = FacetNormal(mesh)

marker = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for c in cells(mesh):
marker[c] = c.midpoint().y() <= 1.25e-2

submesh1 = MeshView.create(marker, 1)
submesh2 = MeshView.create(marker, 0)

V_dd = FunctionSpace(submesh1, "Lagrange", 1)   #for drift-diffusion equation
V_poisson = FunctionSpace(mesh, "Lagrange", 1)   #for Poisson's equation
ME = MixedFunctionSpace(V_dd, V_dd, V_poisson)

v = TestFunctions(ME)
u = Function(ME)
u_old = Function(ME)

u1 = u.sub(0)
u2 = u.sub(1)
Phi = u.sub(2)

u1_old = u_old.sub(0)
u2_old = u_old.sub(1)
Phi_old = u_old.sub(2)

dx1 = Measure("dx", domain = submesh1)
dx2 = Measure("dx", domain = mesh)
ds_dd = Measure('ds', domain = submesh1)
normal_dd = FacetNormal(submesh1)

u1_old = interpolate(Expression('std::log(1e13)', degree = 1), V_dd)
u2_old = interpolate(Expression('std::log(1e13+5e18*exp(-(pow(x[0], 2)+pow(x[1]-1e-2, 2))/pow(0.4e-3, 2)))', degree = 1), V_dd)

Phi_cathode = Constant(0.0)   #[V]
Phi_anode = Constant(18750.0) #[V]

def Potential_Cathode(x, on_boundary):
if near(x[1], 0.0) and on_boundary:
return True

def Potential_Anode(x, on_boundary):
if near(x[1], 1.35e-2) and on_boundary:
return True

Potential_Cathode_bc = DirichletBC(V_poisson, Phi_cathode, Potential_Cathode)
Potential_Anode_bc = DirichletBC(V_poisson, Phi_anode, Potential_Anode)

bcs_potential = [Potential_Cathode_bc, Potential_Anode_bc]

epsilon = Expression('x[1] - 1.25e-2 <= tol ? eps1*eps0 : eps2*eps0', eps1 = 1.0, eps2 = 9.0, eps0 = 8.854187817e-12, tol = DOLFIN_EPS, degree = 1)

rho_p = project((exp(u2_old)-exp(u1_old)) * elementary_charge, V_dd)
rho =  Expression('x[1] - 1.25e-2 <= tol ? rho1 : rho2', rho1 = rho_p, rho2 = 0.0, tol = DOLFIN_EPS, degree = 1)

solve(F_potential == 0, Phi, bcs_potential)   #solving Poisson equation to get the field required for transport and rate coefficients. This works.

E_mag = sqrt(inner(E, E))

F1 = 2*pi*exp(u1)*((u1 - u1_old)/dt)*v[0]*r*dx1\
- 2*pi*alpha(E_mag)*mu(E_mag)*E_mag*exp(u1)*v[0]*r*dx1 + 2*pi*dot(-mu(E_mag) * E, normal_dd)*exp(u1)*v[0]*r*ds_dd

F2 = 2*pi*exp(u2)*((u2 - u2_old)/dt)*v[1]*r*dx1\

F = F1 + F2 + F_potential

solve(F == 0, u, bcs_potential, solver_parameters={"nonlinear_solver": "newton"})   #This one reports the error


On a computing server with Xeon processor and docker image for the mixed-dimensional branch, this code gives the following error:

Traceback (most recent call last):
File "Benchmark.py", line 95, in <module>
solve(F == 0, u, bcs_potential)
File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 233, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 298, in _solve_varproblem
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 'KSPSolve'.
*** Reason:  PETSc error code is: 56 (No support for this operation for this object type).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  39253449ee544fe04eb8bf387eab115e094d325e
*** -------------------------------------------------------------------------


I tried to specify solver parameters such as solver_parameters={"nonlinear_solver": "newton", "newton_solver": {"linear_solver": "lu"}}, but the problem remains. Any suggestion how to solve this problem is appreciated.

Hi,
I have push minor fixes in the mixed-dimensional branch that might fix your issue (the Docker container has been updated accordingly, tag latest).

Also, mshr should now be available in the Docker container.

Thank you for the changes, MWE now works. However, when I try to solve my original problem, I get the following PETSc 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 /tmp/dolfin/dolfin/la/PETScMatrix.cpp.


If I try to use iterative solver, I get a different error:

*** Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 73 (Object is in wrong state).
*** Where: This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.

Can the problem be solved by using SNES, like described here?

Can you give more details on what is different between the previous MWE and the original problem ? The geometry I suppose, is there a difference in term of variational formulation ?
The error you get seems to happens during the assembly, and if so I don’t think using SNES would make any difference if the system is not correctly assembled.

Yes, you are right, the geometry is different. In the original problem, the geometry is as in the first post. In term of variational formulation, the code is basically the same. The only thing that I can think of that could cause the error is a mistake in defining surface integrals in a weak form.

Is it correct to define them like this?

F = ... - 2*pi*2.0/(1.0+ref_coeff[0][i])*gamma[0]*Max(dot(Flux, normal_plasma), 0)*v[i]*r*ds_plasma(0)\
- 2*pi*2.0/(1.0+ref_coeff[1][i])*gamma[1]*Max(dot(Flux, normal_plasma), 0)*v[i]*r*ds_plasma(1)


where:

submesh_plasma = MeshView.create(marker, 0)
normal_plasma = FacetNormal(submesh_plasma)


and Max() is custom made function for finding the maximum.

Could this be a cause for the error during the assembly?

I cannot see any clear reason why that could cause an error, at least from the code samples you provided. Difficult to tell more without having the full code.
I would suggest trying to isolate the term(s) in the variational form that cause the error (are the surface integrals you gave the only terms that are different from your MWE ?)

Thank you for the suggestion and help with this problem. I will try to identify which part of the variational form causes the error, by isolating terms. Even though the code is same as in MWE, it could have been possible that I made a mistake when adapting old code for mixed-dimensional branch.

Hello,

I’m trying to solve a poisson equation with non-linear boundary conditions in one subdomain of a domain.

First, I tried solving the poisson equation in the full domain to check if it is working. The code works fine.

But when i try to solve it only in a subdomain, i’m not able to get it to work. I’m trying to use a submesh for the subdomain. There is even no change in the geometry of the domains in the first and second case. I’m making some mistake while integrating the poisson solver with the submesh feature.

I have attached the working code in a single domain first, and then the code to solve in in only one subdomain (not working).

## Working (Single domain)

from future import print_function
from dolfin import *
import matplotlib.pyplot as plt
from fenics import *
from mshr import *
import numpy as np
from ufl import nabla_div
from ufl import sinh
import random

domain = Rectangle(Point(0,0), Point(100e-6, 200e-6))
mesh = generate_mesh(domain,100)
V = FunctionSpace(mesh, ‘P’, 1)

f = Constant(0.0)

i0 = 12.6
r = Constant(2*i0)

iApp = -100
g = iApp

F = 96485.33
R = 8.314
T = 300
C = Constant(F/(2RT))

kappa = 1

class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[1]) < tol

class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[1] - 200e-6) < tol

bottom_boundary = BottomBoundary()
top_boundary = TopBoundary()

boundary_parts = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
bottom_boundary.mark(boundary_parts, 1)
top_boundary.mark(boundary_parts, 2)
ds = Measure(“ds”)[boundary_parts]

u = TrialFunction(V)
v = TestFunction(V)

u = Function(V)

bcs = []

J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

import matplotlib.pyplot as plt
p = plot(u, mode=‘color’)
plt.colorbar§
plt.title(r"\Phi",fontsize=26)
plt.show()

Solving on subdomain (not working)

from future import print_function
from dolfin import *
from fenics import *
from mshr import *
import numpy as np
from ufl import nabla_div
from ufl import sinh
import random

domain = Rectangle(Point(0,0), Point(100e-6, 400e-6))
mesh = generate_mesh(domain,100)

f = Constant(0.0)

i0 = 12.6
r = Constant(2*i0)

iApp = -100
g = iApp

F = 96485.33
R = 8.314
T = 300
C = Constant(F/(2RT))

kappa = 1

class OmegaMiddle(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 200e-6
omegaM = OmegaMiddle()

# Initialize mesh function for interior domains

subdomains = MeshFunction(‘size_t’, mesh, 2)
subdomains.set_all(0)
omegaM.mark(subdomains, 1)

dx = Measure(‘dx’)[subdomains]

submesh_omegaL = SubMesh(mesh, subdomains, 1)

# define boundaries for the subdomain

class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[1]) < tol
Bottom_boundary = BottomBoundary()

class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[1]-200e-6) < tol
Top_boundary = TopBoundary()

omegaL_boundaries = MeshFunction(‘size_t’, submesh_omegaL,1)
omegaL_boundaries.set_all(0)
Bottom_boundary.mark(omegaL_boundaries, 1)
Top_boundary.mark(omegaL_boundaries, 2)

ds_omegaL = Measure(‘ds’)[omegaL_boundaries]

mydomains = MeshFunction(‘size_t’, submesh_omegaL,2)
mydomains.set_all(0)
dx_subdomain = Measure(‘dx’)[mydomains]

V = FunctionSpace(submesh_omegaL, ‘P’, 1)
u = TrialFunction(V)
v = TestFunction(V)

u = Function(V)

bcs = []

J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

import matplotlib.pyplot as plt
p = plot(u, mode=‘color’)
plt.colorbar§
plt.show()

Please encapsulate your code with ` and format it such that it is runnable by copy paste (indentation, correct use of " , etc), as it makes debugging this easier for other people.

1 Like

@Sabarish_Gopi as you made a new topic for this question, please remove the post above.

For those who might face such an issue anytime in future, please check the link below:

1 Like