Coupled Equations Over Multiple Domains

variational formulation
1

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()
2 Likes
2

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.

3 Likes
3

Dear @Alex and @cdaversin,

I am using docker containr with the image ceciledc/fenics_mixed_dimensional:latest as suggested by user cdaversin for solving a mixed dimensional problem.

Thank you for this code. I am working with mixed dimensional examples and running this code. While the code is working perfectly well with no errors, I wanted to touch base on validation with a known solution.

For that I chose, f = 1.0, and the boundary values on the left boundary a0, b0 =1. When I do the analytical solution for this system, the solutions of a and b as a function of x-direction coordinate takes the form:

a = 0.5(x^2) - x + 1 for 0<x<1 (Defined in the entire square domain)
b = 0.5(x^2) - 0.5x + 1 for 0<x<0.5 (Defined in the left half of the square domain)

Based on the analytical solution, the value of a should vary from a= 1 at x=0 to a = 0.5 at x=1. However, when I run the code, I do not get this match instead the value from the simulation changes from a = 1 to a=1.5 as x goes from 0 to 1. I am attaching the simulation results as well as how the analytical solution was derived here. Please let me know your thoughts on what I might be doing wrong.

ua
ua3159×1772 47.4 KB

20200409_215102
20200409_2151022448×3264 1.92 MB

Code:

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(1.0)
# 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)

## Export last solution
out_sub0 = File("leftRightAlex_ua.pvd")
out_sub1 = File("leftRightAlex_ub.pvd")
# out_exact = File("poisson-lm-2D-exact.pvd")
out_sub0 << u.sub(0)
out_sub1 << u.sub(1)