Coupled Equations Over Multiple Domains

The problem I am tring to simulate involves multiple domains. It consists of 3 layers stacked on one another. One set of equations must be solved over all three layers, and another set of equations only applied to the middle layer. In that middle layer the equations are coupled, and must therefore be solved simultaneously.

In order to figure out how to do this, I have started with a much simpler problem. 2 diffusion equations that are not even coupled:

  • domain 1, D2, is unit square
  • domain 2, D2, is the left half of D1
  • basic diffusion equation of species “a” over D1
    • del(a) = 1
    • a = a0 @ x=0
    • 0 gradient on all other boundaries
  • basic diffusion equation of species “b” over only D2
    • del(b) = 1
    • b = b0 @ x=0
    • 0 gradient on all other boundaries

I have included my code below. It does not work.

Could someone please let me know what I am doing wrong with my current code, or if this is even the correct approach in the first place.

Resulting garbage output

In[1]: a.vector()[:]
Out[1]: array([nan, nan, nan, ...,  1.,  1.,  1.])
In[2]: b.vector()[:]
Out[2]: array([nan, nan, nan, ...,  1.,  1.,  1.])

Code

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt


# Define Boundary equations
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

# Define mesh
resolution = 64

full_domain = Rectangle(Point(0, 0), Point(1, 1))
b_domain = Rectangle(Point(0, 0), Point(0.5, 1))

full_domain.set_subdomain(1, b_domain)

mesh = generate_mesh(full_domain, resolution)

# Initialize boundary instance(s)
left_boundary = LeftBoundary()

# Define domains
domains = MeshFunction("size_t", mesh, 2, mesh.domains())

# Define boundaries on main mesh
boundaries = MeshFunction("size_t", mesh, 1)
boundaries.set_all(0)
left_boundary.mark(boundaries, 1)

# Plot subdomains to make sure that's working properly
plot(domains)
plt.show()

# Generate dx and ds for all meshes and submeshes
dx = Measure('dx', domain=mesh, subdomain_data=domains)

# Define input Data
a0 = 1
b0 = 1
f_a = Constant(1)
f_b = Constant(1)

# Define function space and basis functions
degree: int = 1
fe = FiniteElement("CG", mesh.ufl_cell(), degree=degree)
me = MixedElement([fe, fe])
V = FunctionSpace(mesh, me)

# Define trial and test functions
a, b = TrialFunction(V)
v_a, v_b = TestFunction(V)

# Define Dirichlet Boundary conditions
bcs = [DirichletBC(V.sub(0), a0, boundaries, 1),
       DirichletBC(V.sub(1), b0, boundaries, 1)]

# Define variational form
F = inner(grad(a), grad(v_a))*dx(0) + f_a*v_a*dx(0) + \
    inner(grad(b), grad(v_b))*dx(1) + f_b*v_b*dx(1)

# Seperate LHS and RHS of F
A, L = lhs(F), rhs(F)

# Solve the problem
u = Function(V)
solve(A == L, u, bcs)
a, b = u.split()

#Graph the results

plt.figure()

im = plot(a)
ax = plt.gca()
cbar = ax.figure.colorbar(im, ax=ax)
cbar.ax.set_ylabel("", rotation=-90, va="bottom")
plt.title("a")
plt.show()

plt.figure()

im = plot(b)
ax = plt.gca()
cbar = ax.figure.colorbar(im, ax=ax)
cbar.ax.set_ylabel("", rotation=-90, va="bottom")
plt.title("b")
plt.show()
1 Like

Solution provided by @hernan in the discord group. I am placing it here so others can see it.

The solution to this is to use the mixed dimensional branch. I compiled it from source, but a docker image is also available:

docker pull ceciledc/fenics_mixed_dimensional

I was able to implement the problem I described above. Here is the code for it:

from dolfin import Constant, between, File, Function, inner, Measure, MixedFunctionSpace, \
    MeshFunction, nabla_grad, near, SubDomain, RectangleMesh, Point, MeshView, FunctionSpace, DirichletBC,\
    TestFunctions, TrialFunctions, solve

# Width of total domain
width = 1.0
# Height of total domain
height = 1.0
# Total resolution
resolution = 16
# Forcing function
f = Constant(10)
# Boundary values for the two equations
a_0 = 1
b_0 = 1


# Define SubDomain
class LeftDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, height/2))


# Define Boundary
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)


# Mesh for the entire domain
mesh_full = RectangleMesh(Point(0, 0),
                          Point(width, height),
                          resolution,
                          resolution)

# Define domains
domains = MeshFunction("size_t", mesh_full, 2)
domains.set_all(0)
left_domain = LeftDomain()
left_domain.mark(domains, 1)

# Define boundaries
boundaries_full = MeshFunction("size_t", mesh_full, 1)
boundaries_full.set_all(0)
left_boundary = LeftBoundary()
left_boundary.mark(boundaries_full, 1)

# Create mesh for the left domain
mesh_left = MeshView.create(domains, 1)

# Define boundaries for left domain
boundaries_left = MeshFunction("size_t", mesh_left, 1)
boundaries_left.set_all(0)
left_boundary_2 = LeftBoundary()
left_boundary_2.mark(boundaries_left, 1)

# Create the function spaces
W_a = FunctionSpace(mesh_full, "Lagrange", 1)
W_b = FunctionSpace(mesh_left, "Lagrange", 1)
W = MixedFunctionSpace(W_a, W_b)

# Define boundary conditions
bc_a = DirichletBC(W.sub_space(0), a_0, boundaries_full, 1)
bc_b = DirichletBC(W.sub_space(1), b_0, boundaries_left, 1)
bcs = [bc_a, bc_b]

# Create all of the weak form variables
u = Function(W)
ua, ub = TrialFunctions(W)
va, vb = TestFunctions(W)

# Need to manually specify these.
dxa = Measure("dx", domain=W.sub_space(0).mesh())
dxb = Measure("dx", domain=W.sub_space(1).mesh())

# The weak formulation
a = inner(nabla_grad(ua), nabla_grad(va)) * dxa + inner(nabla_grad(ub), nabla_grad(vb)) * dxb
L = f * va * dxa + f * vb * dxb

solve(a == L, u, bcs)

File("sd/ua.pvd") << u.sub(0)
File("sd/ub.pvd") << u.sub(1)

Additionally, you can couple the two equations D2. Though coupling so far appears to be very finicky, especially when trying to run with MPI.

1 Like