Mixed dimensional Diffusion Problem

[Ekrem_Ekici]

Hi there,

I would like to model 2D diffusion problem with spatially changing boundary condition. That spatial change is the result of 1D diffusion problem solution on that boundary. When I have a certain values over that 1D boundary, I would like to update the markers to make boundary smaller after every timestep (The problem is time-dependent) .

I would be grateful if you can guide me what approach would work best as of dolfinx v0.9.0?

Should I follow mixed function spaces with variational approach?

I checked Dokken’s great tutorial about Coupling PDEs of multiple dimensions, but need further guidance about how I can tackle the diffusion problem.

Thank you!

[dokken]

Hi Ekrem,

I would probably solve this as direct, coupled problem using the mixed-dimensional approach.

@Ekrem_Ekici

I checked Dokken’s great tutorial about Coupling PDEs of multiple dimensions, but need further guidance about how I can tackle the diffusion problem.

It would be similar to the above problem, but add a ufl.conditional that would turn on or off the boundary condition spatially.

Could you write down the mathematical formulation of your two problems, and how they relate through a boundary condition?
Is it a Dirichlet condition, or a Neumann/Robin kind?

[Ekrem_Ekici]

Dear Dokken,

Thank you dor your quick response, please find the schematic of the problem below;

I would like to solve:

\frac{\partial \phi_1}{\partial t} = D_1\nabla \cdot \nabla \phi_1 \quad \text{in $\Omega_1$},

and

\frac{\partial \phi_2}{\partial t} = D_2\nabla \cdot \nabla \phi_2 + s\phi_2 (1-\phi_2/\phi_0) \quad \text{in $\Omega_2$},

at the same time. Both domains have symmetry (Neumann) BC on the left. Blue domain has a BC of \phi_2(right, t)=1.0 on the right end and zero elsewhere. As time progress, blue domain fills up with 1.0s and the cells in \Omega_2 with value 1.0 should become Neumann BC for \Omega_1.

What I made so far is to manually shrink the blue domain, as 2D MWE for full domain shows below:

from dolfinx.fem import functionspace, Constant, Function, form, assemble_matrix
from dolfinx.fem.petsc import assemble_matrix, apply_lifting, set_bc, assemble_vector, create_vector
from petsc4py import PETSc
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI
from dolfinx import default_scalar_type
from ufl import TrialFunction, TestFunction, Measure, dot, grad
from dolfinx.io import XDMFFile
from dolfinx import fem
import numpy as np
import sys
import os
import gmsh

projectRoot = os.path.join(os.getcwd())
resultsDir = os.path.join(projectRoot,'ResultsDir')
screenshotsDir = os.path.join(resultsDir,'screenshots')
test_name = "transient_2d.xdmf"
xdmf_filename = "/"+test_name 

def twoDimHydro(lx, ly, l_drug, drug_start_x, lc=0.005):
    path = os.path.dirname(os.path.abspath(__file__))
    mesh_dir = "/MeshDir"
    mesh_name = "/twodim"

    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add(__name__)
    
    '''
             p_drug_start    p_drug_end
    1-------------|--------------|-------------------2
    |     l1           l_drug          l_2          |
    |                                               |
    |                                               |
    | l_4                                           | l_3
    |                                               |
    |-----------------------------------------------|
    4                   l_wound                      3
    '''
    
    p_1 = gmsh.model.occ.addPoint(0, ly, 0)
    p_2 = gmsh.model.occ.addPoint(lx, ly, 0)
    p_drug_start = gmsh.model.occ.addPoint(drug_start_x, ly, 0)
    p_drug_end =   gmsh.model.occ.addPoint(drug_start_x+l_drug, ly, 0)
    p_3 = gmsh.model.occ.addPoint(lx, 0, 0, lc/5)
    p_4 = gmsh.model.occ.addPoint(0, 0, 0, lc/5)

    l_1 = gmsh.model.occ.addLine(p_1, p_drug_start)
    l_drug = gmsh.model.occ.addLine(p_drug_start, p_drug_end)
    l_2 = gmsh.model.occ.addLine(p_drug_end, p_2)
    l_3 = gmsh.model.occ.addLine(p_2, p_3)
    l_wound = gmsh.model.occ.addLine(p_3, p_4)
    l_4 = gmsh.model.occ.addLine(p_4, p_1)

    cl_1 = gmsh.model.occ.addCurveLoop([l_1,l_drug,l_2, l_3, l_wound, l_4])
    s_1 = gmsh.model.occ.addPlaneSurface([cl_1])
    
    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(1, [l_1,l_2,l_3,l_4], tag=1)  # all other tag
    gmsh.model.addPhysicalGroup(1, [l_drug], tag=2)  # drug tag
    gmsh.model.addPhysicalGroup(1, [l_wound], tag=3)  # wound tag
    gmsh.model.addPhysicalGroup(2, [1], tag=1)  # hydrogel tag

    gmsh.option.setNumber("Mesh.MeshSizeMax", lc)
    gmsh.option.setNumber("Mesh.Optimize", 1)

    gmsh.model.mesh.generate(2)

    if "-nopopup" not in sys.argv:
        gmsh.fltk.run()

    return gmsh.model

lx=1
ly=1/2
l_drug=0.05
drug_start_x=lx/2-l_drug/2

gmsh_model = twoDimHydro(lx, ly, l_drug, drug_start_x, lc=0.01)

mesh, subdomains, facet_tags = model_to_mesh(gmsh_model, MPI.COMM_WORLD, rank=0, gdim=2)

drug_concentration = +1.0
wound_concentration = 0.5
other_concentration = 0.0
other_flux = 0.0

degree = 1

V = functionspace(mesh, ("Lagrange", degree))
u = TrialFunction(V)

# Static/tagged BCs (use your existing tags 1 and 2)
u_other = Function(V)
u_other.x.array[:] = other_concentration
u_drug = Function(V)
u_drug.x.array[:] = drug_concentration

bc_other = fem.dirichletbc(
    u_other,
    fem.locate_dofs_topological(V, mesh.topology.dim-1, facet_tags.find(1))
)
bc_drug = fem.dirichletbc(
    u_drug,
    fem.locate_dofs_topological(V, mesh.topology.dim-1, facet_tags.find(2))
)

def build_wound_bc(V, l_wound, value, lx, x_center=None, y_bottom=0.0, atol=1e-12):
    """
    Return a Dirichlet BC for the wound that shrinks symmetrically toward the middle.

    Args:
        V: Function space.
        l_wound: Remaining wound length along the bottom edge.
        value: Dirichlet value to impose on the wound segment.
        lx: Total domain width in x.
        x_center: Center of the wound segment (defaults to lx/2).
        y_bottom: y-coordinate of the bottom boundary.
        atol: tolerance for y≈y_bottom check.

    Behavior:
        - If l_wound <= 0 -> returns None (no wound left).
        - If l_wound > lx -> clamps to lx.
        - Selects bottom boundary dofs with x in [x_left, x_right], where
          x_left = x_center - l_wound/2 and x_right = x_center + l_wound/2.
    """
    if x_center is None:
        x_center = 0.5 * lx

    # Clamp length and short-circuit if nothing left
    L = float(np.clip(l_wound, 0.0, lx))
    if L <= 0.0:
        return None

    x_left  = x_center - 0.5 * L
    x_right = x_center + 0.5 * L

    # Pick boundary dofs on bottom (y≈y_bottom) & within [x_left, x_right]
    def on_wound(x):
        y_ok = np.isclose(x[1], y_bottom, atol=atol)
        x_ok = np.logical_and(x[0] >= x_left, x[0] <= x_right)
        return np.logical_and(y_ok, x_ok)

    dofs = fem.locate_dofs_geometrical(V, on_wound)
    if dofs.size == 0:
        return None  # no wound region (e.g., very small L wrt mesh)

    # Constant Dirichlet value (fast + robust):
    return fem.dirichletbc(PETSc.ScalarType(value), dofs, V)

l_wound = lx

bc_wound = build_wound_bc(V, l_wound, wound_concentration, lx)

xdmf = XDMFFile(V.mesh.comm, resultsDir+xdmf_filename, "w")
xdmf.write_mesh(V.mesh)

# Define initial condition
u_n = Function(V)
u_n.x.array[:] = 0.5
u_n.name = "u_n"

# Define temporal parameters
t = 0  # Start time
T = 1200  # Final time
num_steps = int(T/2)
dt = T / num_steps  # time step size

D = default_scalar_type(2.6e-5)

u, v = TrialFunction(V), TestFunction(V)
dx = Measure("dx", domain=V.mesh, subdomain_data=subdomains)
uh = Function(V)
uh.name = "phi"
uh.interpolate(u_n)
        
xdmf.write_function(uh, t)

dt = T / num_steps  # time step size

f = Constant(V.mesh, PETSc.ScalarType(0))

a = u * v * dx + dt * D * dot(grad(u), grad(v)) * dx
L = (u_n + dt * f) * v * dx

bcs = [bc_other, bc_drug] + ([bc_wound] if bc_wound is not None else [])

bilinear_form = form(a)
linear_form = form(L)

A = assemble_matrix(bilinear_form, bcs=bcs)
A.assemble()
b = create_vector(linear_form)

# Solve
print("\n- Transient problem is assembled.")
print("\n- Starting the solution..")

solver = PETSc.KSP().create(V.mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

for i in range(num_steps):
    t += dt

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, linear_form)

    # e.g., shrink every 30 steps (just an example):
    if (i+1) % 1 == 0:
        l_wound = max(0.0, l_wound - 0.01/6)
        bc_wound = build_wound_bc(V, l_wound, wound_concentration, lx)

    # Build list of BCs for this step
    bcs = [bc_other, bc_drug] + ([bc_wound] if bc_wound is not None else [])

    # Apply Dirichlet boundary condition to the vector
    apply_lifting(b, [bilinear_form], [bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, bcs)

    # Solve linear problem
    solver.solve(b, uh.x.petsc_vec)
    uh.x.scatter_forward()

    print("Iteration: ", i, "Time:", t)

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array

    # Write solution to file
    xdmf.write_function(uh, t)

xdmf.close()

solution = uh

gives;

at t=32 and

at t=340 and

at the end, t=1200.

What I want is to do that shrinking with the solution of 1D diffusion.

What would you recommend me to handle this?

[dokken]

Switching from Dirichlet to neumann makes it a bit painful.
However, if you start by making a splitting scheme, where you:

1. Solve the PDE on the Omega_2 for phi_2, given an initial phi_0
2. Check if avg phi_2=1 in any cell,if so mark cell in submesh and send marker back to parent.
3. Update marker function for dirichletbc in parent mesh to take the new marking into account.
4. Solve pde in Omega_1.

The more advanced approach, which I guess requires a is_neumann = ufl.conditional(ufl.eq(phi_0, 1), 1, 0)
would be to implement Dirichlet weakly with a (1-is_neumann)*nitsche_formulation + is_neumann*neumann_condition

[Ekrem_Ekici]

Thank you Dokken,
By saying splitting scheme, do you mean solving these diffusion problems with separate solvers on the same function space?
From my understanding, I need to generate a submesh from 2D mesh for the 1D diffusion, and I mark the cells as you mentioned and pass it to the parent mesh to impose new BCs.

I would be grateful if you guide me further for more advanced approach.

Also, would the solutions from these 2 approaches differ at the end? If not, I can proceed with the first approach, which seems easy to implement?

[dokken]

@Ekrem_Ekici

By saying splitting scheme, do you mean solving these diffusion problems with separate solvers on the same function space?

A splitting scheme would involve solving the 1D problem first, then the 2D problem with the 1D solution, and iterate until convergence. Such a method might seem simple and great at first, but it is hard to prove and/or guarantee convergence of such a scheme. Especially since the 1D equation is non-linear.

Thus, I would advice going for the fully coupled method, where you first:

1. Define your 2D mesh
2. Define the 1D submesh
3. Create a mixed variational form (using Nitsche for dirichlet bc as described above)
4. Solve the mixed problem.

[Ekrem_Ekici]

Thank you Dokken,
For step 3, is your Coupling PDEs of multiple domains good start to formulate mixed variational form?

I should also mention that, I am solving 2 different diffusion problems. They represent different physics. All I need from 1D diffusion is to get cells with 1.0s. Then non(1.0) cells become BC for 2D problem.
So in that context, do I still need fully coupled method?

[dokken]

I didn’t realize that there was no coupling, except to determine the boundary condition (I mixed phi_0 with phi_1).

In that case, a splitting scheme might be working nicely.
I would still enforce the Dirichlet weakly (using nitsche) as a starting point, as it is then easy to make the condition move in time without having a major refactor of DOLFINx code.

A good starting point is in between the demo you mention and: Add a demo showing how submeshes can be used to represent boundary data by jpdean · Pull Request #3711 · FEniCS/dolfinx · GitHub

To do it with a strong condition, I think the easiest is to:

1. integrate Neumann over full boundary
2. Modify bcs, which means re-assembling both a and apply correct lifting
3. Send new matrix into KsP and solve.

[Ekrem_Ekici]

Thanks again. Also, if I model the changing boundary to be Neumann or Robin, what makes it easier compared to Dirichlet? Do I need to generate weak forms again at each time step because of the change in ds?

[dokken]

You don’t need to generate the weak forms again at each time step, but you would need to modify the list of boundary conditions on each time step.

Here is an illustration using an indicator function:

# # Ilustration of varying boundary conditions through physical parameters
# Author: Jørgen S. Dokken <dokken@simula.no>
# SPDX-License-Identifier: MIT

from mpi4py import MPI
import dolfinx.fem.petsc
import ufl
import numpy as np

comm = MPI.COMM_WORLD
N = 40
mesh = dolfinx.mesh.create_unit_square(comm, N, N)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
uh = dolfinx.fem.Function(V, name="u")

mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
facet_mesh, entity_map = dolfinx.mesh.create_submesh(
    mesh, mesh.topology.dim - 1, exterior_facets
)[:2]

Q = dolfinx.fem.functionspace(facet_mesh, ("Lagrange", 1))
# Create an initial split of the boundary
indicator = dolfinx.fem.Function(Q, name="indicator")
indicator.interpolate(lambda x: x[0] + 0.6)

is_dirichlet = ufl.conditional(ufl.ge(indicator, 1.0), 1, 0)

# Define an analytical solution to enforce bcs
uD = dolfinx.fem.Function(V)
x = ufl.SpatialCoordinate(mesh)
u_ex = 1 + x[0] ** 2 + 2 * x[1] ** 2
uD.interpolate(dolfinx.fem.Expression(u_ex, V.element.interpolation_points))
f = -ufl.div(ufl.grad(u_ex))
n = ufl.FacetNormal(mesh)
g = -ufl.dot(ufl.grad(u_ex), n)

dx = ufl.dx(domain=mesh)
ds = ufl.ds(domain=mesh)
F = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
F -= f * v * dx
F += (1 - is_dirichlet) * g * v * ds
h = 2 * ufl.Circumradius(mesh)
alpha = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(10))
F += -is_dirichlet * ufl.inner(n, ufl.grad(u)) * v * ds
F += (
    -is_dirichlet * ufl.inner(n, ufl.grad(v)) * (u - uD) * ds
    + is_dirichlet * alpha / h * ufl.inner(u - uD, v) * ds
)

a, L = ufl.system(F)

problem = dolfinx.fem.petsc.LinearProblem(
    a,
    L,
    bcs=[],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
    },
    u=uh,
    entity_maps=[entity_map],
    petsc_options_prefix="mwe_",
)
problem.solve()

vtx_u = dolfinx.io.VTXWriter(comm, "u.bp", [uh])
vtx_u.write(0.0)

is_dirichlet_func = dolfinx.fem.Function(Q, name="is_dirichlet")
is_dirichlet_expr = dolfinx.fem.Expression(is_dirichlet, Q.element.interpolation_points)
is_dirichlet_func.interpolate(is_dirichlet_expr)
vtx_id = dolfinx.io.VTXWriter(comm, "id.bp", [indicator, is_dirichlet_func])
vtx_id.write(0.0)

L2_u = dolfinx.fem.form(ufl.inner(uh - u_ex, uh - u_ex) * dx)
L2_local = dolfinx.fem.assemble_scalar(L2_u)
L2_u_global = np.sqrt(comm.allreduce(L2_local, op=MPI.SUM))
print(f"L2 error: {L2_u_global:.6e}")

is_dirichlet_area = dolfinx.fem.form(is_dirichlet * ds, entity_maps=[entity_map])
print(
    "Dirichlet area",
    comm.allreduce(dolfinx.fem.assemble_scalar(is_dirichlet_area), op=MPI.SUM),
)

indicator.interpolate(lambda x: x[0] + 0.1)
print(
    "Dirichlet area post change",
    comm.allreduce(dolfinx.fem.assemble_scalar(is_dirichlet_area), op=MPI.SUM),
)

problem.solve()
L2_local = dolfinx.fem.assemble_scalar(L2_u)
L2_u_global = np.sqrt(comm.allreduce(L2_local, op=MPI.SUM))
print(f"L2 error: {L2_u_global:.6e}")
is_dirichlet_func.interpolate(is_dirichlet_expr)

vtx_u.write(1.0)
vtx_id.write(1.0)
vtx_u.close()
vtx_id.close()

First time solving:

After update

Of course, if we now use a g that is different than the analytical g, we will clearly see the effect on the solution. For instance, using: g = ufl.sin(x[0]) * ufl.cos(x[1] * ufl.pi):