Hello All,
I am solving a 1D two-domain transport problem in legacy FEniCS. The physical problem is a Poisson-Nernst-Planck type system with one domain representing a gel and the other representing a supernatant. The two domains are represented on the same interval mesh, with the interface at x=1. The supernatant coordinate is reflected/scaled, so the supernatant top is at x=0 and the interface is also at x=1.
At the interface, I want to impose continuity of concentration and continuity of normal flux. For one species, schematically,
c_g = c_s
and
J_g · n_g + J_s · n_s = 0.
The flux is not prescribed external data; it is the constitutive Nernst–Planck flux, for example diffusion plus electromigration. In a simpler diffusion-only diagnostic,
J_g = -D_g grad(c_g) and J_s = -D_s grad(c_s).
So far, I have been using Dirichlet–Neumann type interface treatment. I impose concentration continuity strongly on the gel-side variable, and I impose the corresponding flux as a weak boundary integral in the supernatant equation. A simplified, diffusion-only version of the code looks like this:
# gel concentration: u1
# supernatant concentration: u2
# k_g/k_s: diffusion constants in gel/supernatant
# interface is ds(1); supernatant top is ds(0)
g1 = -k_g*dot(grad(u1), n)
F1 = k_g*inner(grad(u1), grad(v1))*dx
F2 = k_s*inner(grad(u2), grad(v2))*dx - ccr*g1*v2*ds(1)
bc_u1 = DirichletBC(V.sub(0), u(1.0)[1], bx1) # u1 = u2 at interface
bcs = [bc_u1]
solver.solve()
This formulation solves, but when I check the total species amount over both domains, I observe a small but systematic drift over time. The diagnostic is something like:
total_species = assemble(gel_free*dx) + assemble(supernatant_free*dx)
I then tried replacing the strong interface Dirichlet condition and one-sided flux term by a Nitsche-type internal interface formulation. In this version, I remove the strong interface DirichletBC, remove the one-sided interface flux term, and add Nitsche terms for the interface continuity. A simplified version for one diffusing species is:
# gel concentration: u1
# supernatant concentration: u2
# interface is ds(1); supernatant top is ds(0)
k_g = Constant(1.0) # or gel diffusivity
k_s = Constant(1.0) # or supernatant diffusivity
h = CellDiameter(mesh)
gamma_nitsche = Constant(20.0)
jump_u = u1 - u2
jump_v = v1 - v2
# n is the normal of the reference interval.
# The supernatant coordinate has top at x=0 and interface at x=1.
# Therefore, in a common gel-to-supernatant orientation, I use a minus sign
# for the supernatant-side normal flux.
q_g = dot(k_g*grad(u1), n)
q_s = -dot(k_s*grad(u2), n)
qv_g = dot(k_g*grad(v1), n)
qv_s = -dot(k_s*grad(v2), n)
avg_q = 0.5*(q_g + q_s)
avg_qv = 0.5*(qv_g + qv_s)
penalty = gamma_nitsche*((k_g/h) + (k_s/h))
F1 = k_g*inner(grad(u1), grad(v1))*dx
F2 = k_s*inner(grad(u2), grad(v2))*dx
F_nitsche = -avg_q*jump_v*ds(1) \
-avg_qv*jump_u*ds(1) \
+penalty*jump_u*jump_v*ds(1)
F = F1 + F2 + F_nitsche
# No strong DirichletBC at the interface.
# bcs should contain only external/physical Dirichlet conditions.
bcs = []
solver.solve()
With this Nitsche interface form, the total species/charge conservation appears much better in my tests. I removed the electric potential from the code above to keep it simple, but in my actual system, I also solve for the electric field using Poisson’s equation in both domains.
My questions are:
-
Is the first Dirichlet-Neumann interface treatment expected to be non-conservative at the discrete level, even though the continuous interface conditions look correct?
-
Is the Nitsche form above a legitimate way to impose internal interface concentration continuity while obtaining a conservative numerical interface flux?
-
Am I correct that this is not “using Nitsche to impose a prescribed Neumann condition,” but rather using Nitsche for an internal interface problem where both traces are unknown, and the flux continuity is part of the consistent interface formulation?
Any comments or experience with conservative interface coupling in FEniCS would be greatly appreciated.