# Diffusion coupled with boundary equation (mixed-dimensional)

**URL:** <https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580>\
**Category:** variational formulation\
**Created:** [April 14, 2021, 2:47pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580 "2021-04-14T14:47:33Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![cm324](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cm324/32/2541_2.png) [@cm324](https://fenicsproject.discourse.group/u/cm324)\
**Post date:** [April 14, 2021, 2:47pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/1 "2021-04-14T14:47:33Z")

</div>

I am trying to solve the 2D Diffusion equation with a 1D differential equation governing the flux through the bottom wall. To solve this I am using the cecile/mixed-dimensional branch installed through docker (tag 27-03-20) .

**Problem:**

> In 2D medium:
> 
> \frac{\partial c}{\partial t} = \frac{\partial^2 c}{\partial x^2} + \frac{\partial^2 c}{\partial y^2}
> 
> Bottom Wall:
> 
> \frac{\partial b}{\partial t} = c\_{y=0} (1-b) - b
> 
> \Gamma\_{bottom} : \frac{\partial b}{\partial t} = -\frac{\partial c}{\partial y}
> 
> Other Walls:
> 
> \Gamma\_{walls} : \mathbf{n}\cdot \nabla c = 0

**Variational Formulation:**

> F = \int\_{\Omega} \left[\frac{1}{\Delta t}\left(c-c^n \right) v\_c +\frac{\partial c}{\partial x}\frac{\partial v\_c}{\partial x} + \frac{\partial c}{\partial y}\frac{\partial v\_c}{\partial y} \right] d\Omega - \int\_{\Gamma\_s} \frac{\partial c}{\partial y} v\_c ~d\Gamma = 0
> 
> F \text{ +=} \int\_{\Gamma\_s} \frac{1}{\Delta t}(b-b^n) w - g(c,b) w ~d\Gamma = 0
> 
> where,
> 
> \frac{\partial c}{\partial y} = -g(c,b) \qquad , \qquad g(c,b) = c(1-b)-b

**MWE:**

```python
from dolfin import *

mesh = UnitSquareMesh(10, 10)

# Domains and Boundaries
domains = MeshFunction('size_t', mesh, mesh.topology().dim()) ; domains.set_all(0)
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1) ; boundaries.set_all(0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0, DOLFIN_EPS) and on_boundary

Bottom = Bottom() ; Bottom.mark(domains,1) ; Bottom.mark(boundaries,1)

# SubMesh
marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
Bottom.mark(marker,1)
SubMesh = MeshView.create(marker,1)

domains_submesh = MeshFunction('size_t', SubMesh, SubMesh.topology().dim(), 0)

# Define the mixed function space
Q1 = FunctionSpace(mesh, "Lagrange", 1)
Q2 = FunctionSpace(SubMesh, "Lagrange", 1)
W = MixedFunctionSpace(Q1,Q2)

# Measures
dV = Measure("dx", domain=W.sub_space(0).mesh(), subdomain_data = domains)
dL = Measure("dx", domain=W.sub_space(1).mesh(), subdomain_data = domains_submesh)

# Initial Conditions
c_n = interpolate(Constant(1), W.sub_space(0))
b_n = interpolate(Constant(0), W.sub_space(1))

# Define test functions
(v_c , v_b) = TestFunctions(W)

# Functions at Current time
c = Function(W.sub_space(0))
b = Function(W.sub_space(1))

# Constants
t = 0 ; dt = 0.001
Δt = Constant(dt)

def g(c,b): return c*(1-b)-b

# Variational Form
F = (c-c_n)/Δt*v_c*dV + dot(grad(c),grad(v_c))*dV + g(c,b)*v_c*dL

F += (b-b_n)/Δt*v_b*dL - g(c,b)*v_b*dL

sol = Function(W)

# solve(F==0, sol, bcs=[], solver_parameters= {"nonlinear_solver": "newton",
# "newton_solver": {"linear_solver": "mumps"}})

# Manually Define Problem

## Calculate Jacobian Matrix
# from ufl.algorithms import expand_derivatives
F_blocks = extract_blocks(F)
Jacobian = []
for Fi in F_blocks:
    for ui in sol.split():
        d = derivative(Fi, ui)
        # d = expand_derivatives(d) # results in Jacobian =[form[],form[],form[],form[]]
        Jacobian.append(d)

## Define Problem and Solver
problem = MixedNonlinearVariationalProblem(F_blocks, sol.split(), bcs=[], J=Jacobian)
solver = MixedNonlinearVariationalSolver(problem)

## Solver Parameters
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['linear_solver'] = 'mumps'
prm['newton_solver']['convergence_criterion'] = 'residual'

solver.solve()

```

**Error:**

```bash
fenics@73585bea6869:~/shared$ python3 MWE.py 
[problem] create list of residual forms OK
Calling FFC just-in-time (JIT) compiler, this may take some time.
[problem] create list of jacobian forms OK, J_list size = 4
[MixedNonlinearVariationalProblem]::check_forms - To be implemented
Solving mixed nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.559e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
  File "MWE.py", line 82, in <module>
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. ...
*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------

```

**Additional:**  
The reason I am using `mumps` is because of the following error, which arises for the default linear\_solver:

```bash
UMFPACK V5.7.1 (Oct 10, 2014): ERROR: required argument(s) missing

Traceback (most recent call last):
  File "MWE.py", line 81, in <module>
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. ...
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'KSPSolve'.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** -------------------------------------------------------------------------

```

I’d be very grateful for any help. Thanks in advance!

---

<div class="post-metadata">

**Author:** ![cm324](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cm324/32/2541_2.png) [@cm324](https://fenicsproject.discourse.group/u/cm324)\
**Post date:** [April 15, 2021, 5:03pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/2 "2021-04-15T17:03:01Z")

</div>

I have since updated to the latest Docker container (built 15-04-2021). The MWE now throws the following error:

```bash
fenics@9b3d9606959f:~/shared$ python3 MWE.py 
[problem] create list of residual forms OK
[problem] create list of jacobian forms OK, J_list size = 4
[MixedNonlinearVariationalProblem]::check_forms - To be implemented
Solving mixed nonlinear variational problem.
  Newton iteration 0: r (abs) = 9.501e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
  File "MWE.py", line 80, in <module>
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an 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: 53a31ce2dc704214325c989d7dfc75b5e9a9d111
*** -------------------------------------------------------------------------

```

- See related post: [mixednonlinearvariatinalsolver-problem-with-kspsolve](https://fenicsproject.discourse.group/t/mixednonlinearvariatinalsolver-problem-with-kspsolve/2261)

---

<div class="post-metadata">

**Author:** ![cm324](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cm324/32/2541_2.png) [@cm324](https://fenicsproject.discourse.group/u/cm324)\
**Post date:** [April 26, 2021, 5:00pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/3 "2021-04-26T17:00:36Z")

</div>

@cdaversin do you see any correlation between my MWE and the one in the post I linked? I’m just wondering if there is anything I can do on my end to solve this KSPSolve problem. I understand that you are not 24/7 support; I appreciate your time and efforts very much!

---

<div class="post-metadata">

**Author:** ![cm324](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cm324/32/2541_2.png) [@cm324](https://fenicsproject.discourse.group/u/cm324)\
**Post date:** [May 15, 2021, 1:07pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/4 "2021-05-15T13:07:27Z")

</div>

Sorry to bump this, but i’m really stuck here. Maybe there is an alternative way to solve this type of problem, at least until the issue gets resolved. I know I could let b be a 2d function b(x,y), but that method doesn’t seem ideal. Any suggestions?

---

<div class="post-metadata">

**Author:** ![cdaversin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cdaversin/32/44_2.png) [@cdaversin](https://fenicsproject.discourse.group/u/cdaversin)\
**Post date:** [May 20, 2021, 12:44pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/5 "2021-05-20T12:44:14Z")

</div>

Hello, and sorry for my late reply.  
I do see a correlation between your MWE and the post you linked. The change I made regarding the other post didn’t apply to mumps, and it does now. I have pushed a new Docker image including this update.  
It seems that the problem in your MWE doesn’t converge (I haven’t looked at that in details), but it doesn’t throw the error you mentioned anymore.

---

<div class="post-metadata">

**Author:** ![cm324](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cm324/32/2541_2.png) [@cm324](https://fenicsproject.discourse.group/u/cm324)\
**Post date:** [May 20, 2021, 3:52pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/6 "2021-05-20T15:52:12Z")

</div>

Thank you so much @cdaversin for the ongoing support, I am very thankful for your help! About the convergence issue… I did experiment with another package called [multiphenics](https://github.com/mathLab/multiphenics), which you probably know of, and the MWE appears to work fine. The major difference is that multiphenics is using a SNES solver rather than Krylov. Does that mean I just need to use a different solver with the mixed-dimensional branch?

I have attached the working example below in hopes that it will make clear if the convergence issue is due to my variational form, boundary conditions (or lack of), or something else.

```python
# Mulitphenics Example (Working)
from dolfin import *
import multiphenics as mp
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(15, 15)

# Domains and Boundaries
domains = MeshFunction('size_t', mesh, mesh.topology().dim()) ; domains.set_all(0)
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1) ; boundaries.set_all(0)

class OnBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

class OnInterface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0, DOLFIN_EPS) and on_boundary

OnBoundary = OnBoundary() ; OnBoundary.mark(boundaries,1)
OnInterface = OnInterface() ; OnInterface.mark(boundaries,2)

boundary = mp.MeshRestriction(mesh, OnBoundary)
interface = mp.MeshRestriction(mesh, OnInterface)

# Function space
V = FunctionSpace(mesh, "Lagrange", 2)
# Block function space
W = mp.BlockFunctionSpace([V, V], restrict=[None, interface])

# TRIAL/TEST Functions
_cb = mp.BlockTrialFunction(W)
(_c, _b) = mp.block_split(_cb)
cb = mp.BlockFunction(W)
(c, b) = mp.block_split(cb)
vcvb = mp.BlockTestFunction(W)
(v_c, v_b) = mp.block_split(vcvb)

# Measures
dV = Measure("dx")(subdomain_data=domains)
dL = Measure("ds")(subdomain_data=boundaries)
dL = dL(2) # Interface

# Initial Conditions
c_n = interpolate(Constant(1), W.sub(0))
b_n = interpolate(Constant(0), W.sub(1))

# Constants
t = 0 ; dt = 0.01
Δt = Constant(dt)

def g(c,b): return c*(1-b)-b

# Variational Form
F = [(c-c_n)/Δt*v_c*dV + dot(grad(c),grad(v_c))*dV + g(c,b)*v_c*dL ,
    (b-b_n)/Δt*v_b*dL - g(c,b)*v_b*dL]

J = mp.block_derivative(F, cb, _cb)

problem = mp.BlockNonlinearProblem(F, cb, None, J)

solver_parameters = {"nonlinear_solver": "newton",
                          "solver_params": {"linear_solver": "mumps",
                                          "maximum_iterations": 20,
                                          "report": True,
                                          "error_on_nonconvergence": False}}

solver = mp.BlockPETScSNESSolver(problem)
solver.parameters.update(solver_parameters["solver_params"])

solver.solve()

(c, b) = cb.block_split()
plt.figure()
plot(c)
plt.figure()
plot(b)
plt.show()

```

---

<div class="post-metadata">

**Author:** ![cdaversin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cdaversin/32/44_2.png) [@cdaversin](https://fenicsproject.discourse.group/u/cdaversin)\
**Post date:** [May 21, 2021, 11:35am UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/7 "2021-05-21T11:35:37Z")

</div>

It might indeed be due to the solver you are using, but unfortunately SNES is not supported in the mixed-dimensional branch. I can try to investigate your variational form and boundary conditions a bit deeper.

---

<div class="post-metadata">

**Author:** ![cm324](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cm324/32/2541_2.png) [@cm324](https://fenicsproject.discourse.group/u/cm324)\
**Post date:** [May 21, 2021, 12:09pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/8 "2021-05-21T12:09:41Z")

</div>

That would be phenomenal! Essentially I am trying to replicate the workings of this article [https://doi.org/10.1038/nbt1388](https://doi.org/10.1038/nbt1388). Briefly, I am solving the advection-diffusion equations coupled with a Langmuir reaction that happens on the surface. For the sake of completeness I will state the actual mathematical problem here (for the MWE I simply set all parameters to unity). @cdaversin if it helps, I can privately share my project with you.

# Problem

> - Langmuir Isotherm (First order reaction kinectics):
> 
> \frac{\partial b}{\partial t} = k\_{\text{on}}c\_s(b\_{max}-b)-k\_{\text{off}}b
> 
> > where c\_s is the (ambient) target concentration at the sensor surface.
> 
> - Differential Equation:
> 
> \frac{\partial c}{\partial t} = D\nabla^2 c - \mathbf{u}\cdot \nabla c
> 
> - Boundary Conditions:
> 
> \begin{aligned} \Gamma\_{sensor} &:& \mathbf{n}\cdot D\nabla c = \frac{\partial b}{\partial t}\\ \Gamma\_{walls} &:& \mathbf{n}\cdot D\nabla c = 0\\ \end{aligned}
> 
> - Initial Conditions:
> 
> c(\mathbf{x},t)|\_{t=0}=c\_0
> 
> - Prescribed Velocity Field
> 
> \vec{u}(y)=\frac{6Q}{W\_c H^3} y(H-y) ~\hat{x}

## Scaling

> Non-dimensionalize the problem by introducing the scaling parameters:
> 
> x = \frac{x'}{L} \quad , \quad y = \frac{y'}{H} \quad , \quad c=\frac{c'}{c\_0} \quad , \quad \lambda = \frac{H}{L} \quad , \quad t = \frac{D t'}{H^2}
> 
> For the reaction equation we scale the surface concentration of bound analyte by the binding site density (b= b'/b\_{\text{max}}). The following parameters arise naturally from non-dimensionalizing the Langmuir binding kinetics equation:
> 
> \text{Da} = \frac{k\_{\text{on}} b\_{\text{max}}H}{D} \quad , \quad \epsilon = \frac{c\_0 H}{b\_{\text{max}}} \quad , \quad \tilde{c} = \frac{k\_{\text{on}} c\_0}{k\_{\text{off}}}

## The Scaled Problem

> - Differential Equation:
> 
> \frac{\partial c}{\partial t} = \lambda^2 \frac{\partial^2 c}{\partial x^2} + \frac{\partial^2 c}{\partial y^2} - \vec{u}(y) \frac{\partial c}{\partial x}
> 
> \frac{\partial b}{\partial t} = \epsilon \text{Da} \left( c\_s (1-b) - \tilde{c}^{-1} b \right)
> 
> - Boundary Conditions:
> 
> \begin{aligned} \Gamma\_{sensor} &:& \frac{\partial b}{\partial t} = -\epsilon \frac{\partial c}{\partial y} \\ \Gamma\_{walls} &:& \mathbf{n}\cdot \nabla c = 0\\ \end{aligned}
> 
> - Initial Conditions:
> 
> c(\mathbf{x},t)|\_{t=0}=1
> 
> - Prescribed Velocity Field
> 
> \vec{u}(y)= 6 \lambda \text{Pe}\_H y(1-y) ~\hat{x}
> 
> > where \text{Pe}\_H = \frac{Q}{W\_c D} is the dimensionless Péclet number.

## Variational Formulation

> F\_1 = \int\_{\Omega} \left[\frac{1}{\Delta t}\left(c-c^n \right) v\_c + \lambda^2 \frac{\partial c}{\partial x}\frac{\partial v\_c}{\partial x} + \frac{\partial c}{\partial y}\frac{\partial v\_c}{\partial y} + u\_x \frac{\partial c}{\partial x} v\_c \right] d\Omega - \int\_{\Gamma\_s} \frac{\partial c}{\partial y} v\_c ~d\Gamma = 0
> 
> F\_2 = \int\_{\Gamma\_s} \frac{1}{\Delta t}(b-b^n) v\_b - \epsilon \text{Da} \left[c(1-b) - \frac{1}{\tilde c}b \right]v\_b ~d\Gamma = 0
> 
> the flux condition on the boundary of the sensor is weakly enforced through the variational formulation:
> 
> \frac{\partial c}{\partial y} = -D\_a \left[c(1-b) - \frac{1}{\tilde c}b \right]

---

<div class="post-metadata">

**Author:** ![MikeS](https://avatars.discourse-cdn.com/v4/letter/m/919ad9/32.png) [@MikeS](https://fenicsproject.discourse.group/u/MikeS)\
**Post date:** [November 25, 2022, 10:36pm UTC](https://fenicsproject.discourse.group/t/diffusion-coupled-with-boundary-equation-mixed-dimensional/5580/9 "2022-11-25T22:36:33Z")

</div>

I have a problem similar to yours, cm324, with a PDE on a region and with a lower-dimensoinal PDE on its boundary. I have not been able to express this using UFL with a mixed-element (elements defined on the interior and elements defined on the boundary) formulation. Could you please give me an idea of how you did it?