Conditional markers and interface facet handling in DG two-phase model

[Giuliano3110]

Hi,

I am new to Firedrake and trying to solve a 2D two-phase diffusion problem in Firedrake using DG(0) discretization for each phase. The two phases (electrolyte and active material) are distinguished by a phase marker, and fluxes are exchanged across their interfaces using dS terms. The problem is similar to coupled mass transfer with an interfacial reaction term.

However, my solver diverges quickly — even with very small time steps. I suspect the issue is in the way the fluxes (dS terms) are set up. Is anyone familiar with defining PDEs with UFL conditionals.

Below is a minimal reproducible example (no external files needed) that shows the problem.

from firedrake import *
import numpy as np
from ufl import conditional, eq, And

# --- Parameters ---
nx, ny = 10, 10
Lx, Ly = 1.0, 1.0 
Dx, Dy = Lx/nx, Ly/ny # voxel size

# --- Mesh ---
mesh = RectangleMesh(nx, ny, Lx, Ly, quadrilateral=True)
x, y = SpatialCoordinate(mesh)

# --- Function spaces ---
V = FunctionSpace(mesh, "DG", 0)
W = V * V
v_ce, v_cs = TestFunction(W)
w = Function(W, name="w")
w_old = Function(W, name="w_old")
c_e, c_s = split(w)
c_e_old, c_s_old = split(w_old)

# --- Phase marker: left half = phase 0, right half = phase 1 ---
phase_marker = Function(V)
phase_marker.interpolate(conditional(x < 0.5*Lx, 0, 1))

# Indicators
Omega_e = conditional(eq(phase_marker, 0), 1, 0)   # electrolyte
Omega_am = conditional(eq(phase_marker, 1), 1, 0)  # active material

# --- Parameters ---
dt = Constant(0.01)
D_e = Constant(0.01)
D_s = Constant(1.0)
alpha = Constant(1.0)
n = FacetNormal(mesh)
h = CellSize(mesh)

# --- Initial conditions ---
w.sub(0).assign(1000)  # c_e
w.sub(1).assign(21112)  # c_s
w_old.assign(w)

# --- Weak form ---

# Accumulation
F_c_e = conditional(eq(Omega_e, 1), (c_e - c_e_old)/dt * v_ce, 0) * dx
F_c_s = conditional(eq(Omega_am, 1), (c_s - c_s_old)/dt * v_cs, 0) * dx

# Diffusion terms
F_c_e += conditional(And(eq(Omega_e('+'), 1), eq(Omega_e('-'), 1)),
                     D_e * alpha / h('+') * dot(jump(c_e, n), jump(v_ce, n)),
                     0) * dS
F_c_s += conditional(And(eq(Omega_am('+'), 1), eq(Omega_am('-'), 1)),
                     D_s * alpha / h('+') * dot(jump(c_s, n), jump(v_cs, n)),
                     0) * dS

# Interface reaction flux
flux_e_to_s = conditional(
    And(eq(Omega_e('+'), 1), eq(Omega_am('-'), 1)),
    100 * c_e('+')/c_s('-'),
    0.0
)
F_c_e += -flux_e_to_s * v_ce('+') * dS
F_c_s +=  flux_e_to_s * v_cs('-') * dS

flux_s_to_e = conditional(
    And(eq(Omega_am('+'), 1), eq(Omega_e('-'), 1)),
    100 * c_e('-')/c_s('+'),
    0.0
)
F_c_e += -flux_s_to_e * v_ce('-') * dS
F_c_s +=  flux_s_to_e * v_cs('+') * dS

# BC -> Noflow

F = F_c_e + F_c_s
# --- Solver ---
solver_parameters = {
    "snes_monitor": None,
    "ksp_type": "gmres",
    "pc_type": "ilu",
    "snes_rtol": 1e-8,
}

# --- Time loop ---
t = 0.0
for step in range(5):
    solve(F == 0, w, solver_parameters=solver_parameters)
    w_old.assign(w)
    t += float(dt)
    print(f"Step {step} done")

[connorjward]

Can I suggest also posting this on the Firedrake Slack? That has a wider audience so you'll be more likely to find someone who can answer this for you.

[Giuliano3110]

Thanks for the advice. I've just posted it.