# Interface condition issue in multi-material diffusion with mixed solubility laws

**URL:** <https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158>\
**Category:** User Support\
**Tags:** festim2\
**Created:** [February 6, 2026, 3:35pm UTC](https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158 "2026-02-06T15:35:01Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![ChrisSo](https://avatars.discourse-cdn.com/v4/letter/c/91b2a8/32.png) [@ChrisSo](https://festim.discourse.group/u/ChrisSo)\
**Post date:** [February 6, 2026, 3:35pm UTC](https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158/1 "2026-02-06T15:35:01Z")

</div>

Hi all,  
I recently switched to **FESTIM2 (v2.0b1)** and ran into an issue when running multi-material simulations using the discontinuous formulations. While benchmarking against **FESTIM1** , I found that in certain cases the interface condition does not appear to be enforced correctly when using `HydrogenTransportProblemDiscontinuous()`.  
To investigate this, I benchmarked the following three implementations against each other:

- **FESTIM1**
- `HydrogenTransportProblemDiscontinuous()`
- `HydrogenTransportProblemDiscontinuousChangeVar()`

The physical problem was kept simple. I am looking at a 2D multi-material three-body setup:  
**Wall | Fluid | Wall**  
The physics is diffusion only (no advection). I apply a fixed concentration of 1e20 on the left boundary and zero concentration on the right boundary. The wall material is _Eurofer_, and the fluid material is either _Helium_ or _LiPb_. In addition, I also tested dummy material properties (_Dummy1_ for the fluid and _Dummy2_ for the walls) to check the robustness of the formulations.  
The solubility-law combinations I tested are:

- Sieverts | Henry | Sieverts (_Eurofer_ with _Helium_)
- Sieverts | Sieverts | Sieverts (_Eurofer_ with _LiPb_ and _Dummy1_ and _Dummy2_).  
What I observe is that **FESTIM1** and `HydrogenTransportProblemDiscontinuousChangeVar()` give very comparable results, while `HydrogenTransportProblemDiscontinuous()` produces qualitatively different ones.  
The table below shows three test cases I ran. The shown figures describe the concetration profile between the two boundaries.

| Material | FESTIM1 | `DiscontinuousChangeVar()` | `Discontinuous()` |
| --- | --- | --- | --- |
| Eurofer w. LiPb | ![image800x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/dc34e0235e3386536bac0685ebc73e0f03231816.png) | ![F2 changevar LiPb600x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/16c4cf73b01196a4411c73ad8e69faece76962a5.png) | ![F2 nitsche LiPb613x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/760f53116c2260264cf1cd958fc07c4d45f78eeb.png) |
| Eurofer w. Helium | ![F1 Helium613x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/3a01a130dafa69254949155ab22657df86b75776.png) | ![F2 changevar Helium613x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/09aca5c2f8e94af3dbd22a48b764c34c81835aaf.png) | ![F2 nitsche Helium600x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/1b94cd2bf99ac456fb6b04895f9772bbf9d217c8.png) |
| Dummy material w. Sieverts only | ![F1 dummy S613x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/7600cf521715ce0447a377ea45eeeedd7e2b19f8.png) | ![F2 changevar dummy S613x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/3c1202441d55be77069e63184c951756ef5035c7.png) | ![F2 nitsche dummy S613x391](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/d33300cee207b8e48613dc51cbae4c171ad157a0.png) |

For the cases for `HydrogenTransportProblemDiscontinuous()`, the change of concentration across the interface is inverse to that of the other two cases when realistic material properties are used. All 3 approaches coinside for the dummy materials values.  
**FESTIM1** and the `ChangeVar()` are similar especially for the first and third case which makes sense as they are solving, as I understand, the same transport equations. In the second case they differ as `HydrogenTransportProblemDiscontinuousChangeVar()` is not considering the solubility law difference.  
I explicitly checked the concentration jump at the interface against the theoretical jump expected from the solubility ratios. Using the simulated concentrations on both sides of the interface, the results from **FESTIM1** and `HydrogenTransportProblemDiscontinuousChangeVar()` match expectations, whereas the results from `HydrogenTransportProblemDiscontinuous()` do not.  
I tried global mesh refinement, but this did not resolve the issue (I am not refining specifically at the interface which might be something to try next).  
To isolate the problem further, I tested the same setups using dummy material properties (case 3). In that case, `HydrogenTransportProblemDiscontinuous()` matches the results from the other two implementations for purely Sieverts-type materials.

What I find particularly confusing is that once realistic material properties are used, `HydrogenTransportProblemDiscontinuous()` seems unable to impose the correct interface condition, even though it behaves as expected with dummy properties.  
I also encountered issues when using the penalty interface method. For cases with different solubility laws, the solver does not converge, regardless of whether realistic or dummy material properties are used. Mesh refinement did not help in this case either.  
At this point, I am unsure what the cause of this is and how to resolve it.  
If anyone has insight into this behavior, knows of constraints on material parameters, or can suggest a recommended strategy, I would appreciate the help.

A minimal working script with all material properties is provided below for FESTIM2.

```auto
import festim as F
from dolfinx.mesh import create_rectangle, create_unit_square
import dolfinx
from dolfinx import plot
import numpy as np
from mpi4py import MPI
import pyvista as pv
from pathlib import Path
import h_transport_materials as htm # Material properties

inlet_id=1
outlet_id=2
wall_id=3
axis_id=4
fluid_id=5
left_wall_id = 6
right_wall_id = 7
interface_l_id = 8
interface_r_id = 9
pipe_id=10
wall_l_thickness = 0.005
wall_r_thickness = 0.005
r_max=0.02
z_max=0.1
nx=100
ny=250
r_min=0.005
z_min=0

#set simulation typ
my_sim = F.HydrogenTransportProblemDiscontinuous()
#my_sim = F.HydrogenTransportProblemDiscontinuousChangeVar()

#generate Mesh
mesh_fenics = create_rectangle(MPI.COMM_WORLD,[[r_min-wall_l_thickness,z_min],[r_max+wall_r_thickness,z_max]],[nx,ny])

#define Material properties
Dummy1 = F.Material(D_0=1,E_D=0,K_S_0=2,E_K_S=0,solubility_law="SIEVERT") #HENRY SIEVERT
Dummy2 = F.Material(D_0=2,E_D=0,K_S_0=3,E_K_S=0,solubility_law="SIEVERT")

S_0=9.0537131e19
D_0=8.83757e-06
He = F.Material(D_0=D_0, E_D=0, K_S_0=S_0, E_K_S=0,solubility_law='HENRY')

D_pbli = (htm.diffusivities.filter(material="lipb").filter(isotope="h").filter(author="reiter"))[0]
S_pbli = (htm.solubilities.filter(material="lipb").filter(isotope="h").filter(author="aiello"))[0]

PbLi = F.Material(
    D_0=D_pbli.pre_exp.magnitude,
    E_D=D_pbli.act_energy.magnitude,
    K_S_0=S_pbli.pre_exp.magnitude,
    E_K_S=S_pbli.act_energy.magnitude,
    solubility_law="SIEVERT"
)

D_eurofer = htm.diffusivities.filter(material="eurofer_97").filter(author="chen")[0]
S_eurofer = htm.solubilities.filter(material="eurofer_97").filter(author="chen")[0]

Eurofer = F.Material(
    D_0=D_eurofer.pre_exp.magnitude,
    E_D=D_eurofer.act_energy.magnitude,
    K_S_0=S_eurofer.pre_exp.magnitude,
    E_K_S=S_eurofer.act_energy.magnitude,
    solubility_law="SIEVERT"
)

#assign surfaces
inlet_boundary = F.SurfaceSubdomain(id=inlet_id, locator=lambda x: (np.isclose(x[1], z_min) & (x[0] >= r_min) & (x[0] <= r_max)))
outlet_boundary = F.SurfaceSubdomain(id=outlet_id,locator=lambda x: (np.isclose(x[1], z_max) & (x[0] >= r_min) & (x[0] <= r_max)))
inner_wall_boundary = F.SurfaceSubdomain(id=axis_id, locator=lambda x: np.isclose(x[0], r_min-wall_l_thickness))
outer_wall_boundary = F.SurfaceSubdomain(id=wall_id, locator=lambda x: np.isclose(x[0], r_max+wall_r_thickness))

#assign Volumes
fluid_subdomain = F.VolumeSubdomain(id=fluid_id,material=He,locator=lambda x: (x[0]>=r_min) & (x[0]<=r_max))
#fluid_subdomain = F.VolumeSubdomain(id=fluid_id,material=PbLi,locator=lambda x: (x[0]>=r_min) & (x[0]<=r_max))
#fluid_subdomain = F.VolumeSubdomain(id=fluid_id,material=Dummy1,locator=lambda x: (x[0]>=r_min) & (x[0]<=r_max))

left_wall = F.VolumeSubdomain(
    id=left_wall_id, material=Eurofer, #Dummy2 Eurofer
    locator=lambda x: x[0] <= r_min
    )
right_wall = F.VolumeSubdomain(id=right_wall_id, material=Eurofer #Dummy2 Eurofer
                               ,locator=lambda x: x[0] >= r_max)

#setup Mesh and interfaces
my_sim.mesh=F.Mesh(mesh_fenics)
my_sim.subdomains = [
    outer_wall_boundary,
    inner_wall_boundary,
    inlet_boundary,
    outlet_boundary,
    fluid_subdomain,
    right_wall,
    left_wall,
    ]

my_sim.interfaces = [
    F.Interface(interface_l_id, (fluid_subdomain, left_wall), penalty_term=1000),
    F.Interface(interface_r_id, (fluid_subdomain, right_wall), penalty_term=1000),
]
my_sim.method_interface = "nitsche" #penalty nitsche

my_sim.surface_to_volume = {
    inlet_boundary: fluid_subdomain,
    outlet_boundary: fluid_subdomain,
    inner_wall_boundary: left_wall,
    outer_wall_boundary: right_wall,
}

#define Species
H = F.Species("H")
my_sim.species = [H]

for species in my_sim.species:
    species.subdomains=[fluid_subdomain, left_wall, 
                        right_wall
                        ]

my_sim.define_meshtags_and_measures()

# define Conditions and settings
my_sim.temperature=800

my_sim.boundary_conditions= [
    F.FixedConcentrationBC(species=H,subdomain=inner_wall_boundary,value=1e20),
    F.FixedConcentrationBC(species=H,subdomain=outer_wall_boundary,value=0),
]

my_sim.settings = F.Settings(atol=1e10, rtol=1e-10, transient=False,max_iterations=500)

# run Simulation
my_sim.initialise()
my_sim.run()

## Post-process ##
def make_ugrid(solution):
    topology, cell_types, geometry = plot.vtk_mesh(solution.function_space)
    u_grid = pv.UnstructuredGrid(topology, cell_types, geometry)
    u_grid.point_data["c"] = solution.x.array.real
    u_grid.set_active_scalars("c")
    return u_grid

if isinstance(my_sim, F.HydrogenTransportProblemDiscontinuous):
    # plotting for HydrogenTransportProblemDiscontinuous
    pv.start_xvfb()
    #pv.set_jupyter_backend("html")

    u_plotter = pv.Plotter()
    u_grid_fluid = make_ugrid(H.subdomain_to_post_processing_solution[fluid_subdomain])
    u_grid_wall_l = make_ugrid(H.subdomain_to_post_processing_solution[left_wall])
    u_grid_wall_r = make_ugrid(H.subdomain_to_post_processing_solution[right_wall])
    u_plotter.add_mesh(u_grid_fluid, show_edges=False)
    u_plotter.add_mesh(u_grid_wall_l,show_edges=False)
    u_plotter.add_mesh(u_grid_wall_r,show_edges=False)
    u_plotter.view_xy()
    u_plotter.add_text("Hydrogen concentration in multi-material problem for Discontinuous", font_size=12)

    if not pv.OFF_SCREEN:
        u_plotter.show()
    else:
        figure = u_plotter.screenshot("concentration.png")

elif isinstance(my_sim, F.HydrogenTransportProblemDiscontinuousChangeVar):
    # plotting for HydrogenTransportProblemDiscontinuousChangeVar
    u_grid = make_ugrid(H.post_processing_solution)
    u_plotter = pv.Plotter()
    u_plotter.add_mesh(u_grid, show_edges=False)
    u_plotter.view_xy()
    u_plotter.add_text("Hydrogen concentration in multi-material problem for ChangeVar", font_size=12)
    u_plotter.show()
else:
    raise TypeError(
        f"Unsupported FESTIM problem type: {type(my_sim)}"
    )

# ChangeVar line profiles: concentration vs x[0] at multiple x[1]
import numpy as np
import matplotlib.pyplot as plt
from dolfinx import geometry

def eval_on_mesh(u, points):
    mesh = u.function_space.mesh
    gdim = mesh.geometry.dim

    pts = np.asarray(points, dtype=np.float64)
    if pts.ndim != 2 or pts.shape[1] != gdim:
        raise ValueError(f"points should be (N, {gdim})")

    # pad to 3D for dolfinx geometry/eval
    if gdim == 2:
        pts3 = np.zeros((pts.shape[0], 3), dtype=np.float64)
        pts3[:, :2] = pts
    else:
        pts3 = pts

    pts3 = np.ascontiguousarray(pts3)

    bb = geometry.bb_tree(mesh, mesh.topology.dim)
    candidates = geometry.compute_collisions_points(bb, pts3)
    colliding = geometry.compute_colliding_cells(mesh, candidates, pts3)

    npts = pts3.shape[0]
    cells = np.full(npts, -1, dtype=np.int32)

    if hasattr(colliding, "links"):
        for i in range(npts):
            links = colliding.links(i)
            if len(links):
                cells[i] = links[0]
    elif hasattr(colliding, "array") and hasattr(colliding, "offsets"):
        arr = colliding.array
        off = colliding.offsets
        for i in range(npts):
            start, end = off[i], off[i + 1]
            if end > start:
                cells[i] = arr[start]
    else:
        raise TypeError("Unsupported AdjacencyList type from dolfinx")

    values = np.full(len(pts), np.nan, dtype=float)
    mask = cells >= 0
    if np.any(mask):
        vals = u.eval(pts3[mask], cells[mask])
        values[mask] = vals.reshape(-1)
    return values

# ---- plotting ----
z_slices = [z_min + 0.02, z_min + 0.05, z_min + 0.08]  
x0 = np.linspace(r_min - wall_l_thickness, r_max + wall_r_thickness, 600)

if isinstance(my_sim, F.HydrogenTransportProblemDiscontinuous):
    # grab discontinuous concentration functions for each subdomain
    u_left = H.subdomain_to_post_processing_solution[left_wall]
    u_fluid = H.subdomain_to_post_processing_solution[fluid_subdomain]
    u_right = H.subdomain_to_post_processing_solution[right_wall]

    fig, ax = plt.subplots(figsize=(7,4))
    for z in z_slices:
        y = np.full_like(x0, z, dtype=float)
        c = np.full_like(x0, np.nan, dtype=float)

        mask = x0 <= r_min
        c[mask] = eval_on_mesh(u_left, np.column_stack([x0[mask], y[mask]]))

        mask = (x0 >= r_min) & (x0 <= r_max)
        c[mask] = eval_on_mesh(u_fluid, np.column_stack([x0[mask], y[mask]]))

        mask = x0 >= r_max
        c[mask] = eval_on_mesh(u_right, np.column_stack([x0[mask], y[mask]]))

        ax.plot(x0, c, label=f"x[1]={z:.3g}")

    ax.axvline(r_min, color="k", lw=0.8, ls="--")
    ax.axvline(r_max, color="k", lw=0.8, ls="--")
    ax.set_xlabel("x[0]")
    #ax.set_yscale("log")
    ax.set_ylabel("Concentration")
    ax.set_title("Concentration profile across walls/fluid")
    ax.legend()
    plt.show()

elif isinstance(my_sim, F.HydrogenTransportProblemDiscontinuousChangeVar):
    # this is already concentration for ChangeVar
    u_c = H.post_processing_solution

    fig, ax = plt.subplots(figsize=(7,4))
    for z in z_slices:
        y = np.full_like(x0, z, dtype=float)
        pts = np.column_stack([x0, y])
        c = eval_on_mesh(u_c, pts)
        ax.plot(x0, c, label=f"x[1]={z:.3g}")

    ax.axvline(r_min, color="k", lw=0.8, ls="--")
    ax.axvline(r_max, color="k", lw=0.8, ls="--")
    ax.set_xlabel("x[0]")
    ax.set_ylabel("Concentration")
    #ax.set_yscale("log")
    ax.set_title("ChangeVar concentration profile")
    ax.legend()
    plt.show()
else:   
    raise TypeError(
        f"Unsupported FESTIM problem type: {type(my_sim)}"
    )

```

Best regards,  
Chris

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [February 6, 2026, 7:59pm UTC](https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158/2 "2026-02-06T19:59:30Z")

</div>

Hi Chris, thanks for reporting!

It seems to me that the penalty term isn’t large enough compared to the magnitude of the residual.

Consider using a penalty of say 1e20. You could also consider scaling the units of the problem to moles instead of atoms

---

<div class="post-metadata">

**Author:** ![ChrisSo](https://avatars.discourse-cdn.com/v4/letter/c/91b2a8/32.png) [@ChrisSo](https://festim.discourse.group/u/ChrisSo)\
**Post date:** [February 6, 2026, 9:20pm UTC](https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158/3 "2026-02-06T21:20:28Z")

</div>

Hey,  
thanks for the answer!  
Changing the penalty term made the sim converge and match the FESTIM1 benchmark, great!

My follow up question would be if there is any way to get an idea on how big the penalty value should be? I changed it before posting but I think the maximum value I tried was like 1e6.

You also mentioned that the penalty was small compared to the residual. Is there a way in FESTIM2 to actually look at or inspect the residual (similar to what was possible in FESTIM1 with log\_level) and then take a decision on the penalty term?

Scalling to moles sounds interesting I will definitly check that out. The reason for doing so would be to decrease the magnitude of the values involved and thereby make the numerics easier. Is that correct?

---

<div class="post-metadata">

**Author:** ![ChrisSo](https://avatars.discourse-cdn.com/v4/letter/c/91b2a8/32.png) [@ChrisSo](https://festim.discourse.group/u/ChrisSo)\
**Post date:** [February 11, 2026, 2:24pm UTC](https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158/4 "2026-02-11T14:24:07Z")

</div>

Hey Remi,

after your initial suggestion worked, I increased the domain to a 5 body problem (Wall | Fluid | Wall | Fluid | Wall), which is closer to my target setup.

Changing the penalty (tested from 1 to 1e35) did not resolve the issue. With identical solubility laws across the domain the simulation converges, but the result looks physically incorrect (Eurofer and LiPb). With mixed solubility laws it does not converge at all (Eurofer and He\_Purge).

Do you have any idea what could cause this?

I also experimented with scaling to moles but that didn’t help.

Below is a minimal working script. I reduced the problem to 1D.

```auto
import festim as F
from dolfinx.mesh import create_rectangle, create_unit_square
import dolfinx
from dolfinx import plot
import numpy as np
from mpi4py import MPI
import pyvista as pv
from pathlib import Path
import h_transport_materials as htm # Material properties

wall_id=3
axis_id=4
fluid_id=5
left_wall_id = 6
right_wall_id = 7
interface_l_id = 8
interface_r_id = 9
purge_id = 11
right_p_wall_id=12
wall_l_thickness = 0.005
wall_p_r_thickness =0.0025
wall_r_thickness = 0.005
r_max=0.02
z_max=0.1
nx=120
ny=200
r_min=0.005
z_min=0
height_pc =0.0025

my_sim = F.HydrogenTransportProblemDiscontinuous()
#my_sim = F.HydrogenTransportProblemDiscontinuousChangeVar()

n = 20
vertices = np.concatenate(
    [
        np.linspace(0,r_min,num=n*2),
        np.linspace(r_min,r_max, num=n*6),
        np.linspace(r_max,r_max+wall_p_r_thickness,num=n),
        np.linspace(r_max+wall_p_r_thickness,r_max+wall_p_r_thickness+height_pc,num=n),
        np.linspace(r_max+wall_p_r_thickness+height_pc,r_max+wall_p_r_thickness+height_pc+wall_r_thickness,num=n*2)
    ]
)

mesh = F.Mesh1D(vertices=vertices)

inner_wall_boundary = F.SurfaceSubdomain1D(id=axis_id,x=0)
outer_wall_boundary = F.SurfaceSubdomain1D(id=wall_id, x = r_max + wall_p_r_thickness + height_pc + wall_r_thickness)

my_sim.mesh=mesh

S_0=9.0537131e19
D_0=8.83757e-06
He_Purge = F.Material(D_0=D_0, E_D=0, K_S_0=S_0, E_K_S=0,
    solubility_law='HENRY'
    )

D_eurofer = htm.diffusivities.filter(material="eurofer_97").filter(author="chen")[0]
S_eurofer = htm.solubilities.filter(material="eurofer_97").filter(author="chen")[0]

Eurofer = F.Material(
    D_0=D_eurofer.pre_exp.magnitude,
    E_D=D_eurofer.act_energy.magnitude,
    K_S_0=S_eurofer.pre_exp.magnitude,
    E_K_S=S_eurofer.act_energy.magnitude,
    solubility_law="SIEVERT"
)

D_pbli = (htm.diffusivities.filter(material="lipb").filter(isotope="h").filter(author="reiter"))[0]
S_pbli = (htm.solubilities.filter(material="lipb").filter(isotope="h").filter(author="aiello"))[0]

PbLi = F.Material(
    D_0=D_pbli.pre_exp.magnitude,
    E_D=D_pbli.act_energy.magnitude,
    K_S_0=S_pbli.pre_exp.magnitude,
    E_K_S=S_pbli.act_energy.magnitude,
    solubility_law="SIEVERT"
)

Dummy1 = F.Material(D_0=1,E_D=0,K_S_0=2,E_K_S=0,solubility_law="SIEVERT") #HENRY SIEVERT
Dummy2 = F.Material(D_0=2,E_D=0,K_S_0=3,E_K_S=0,solubility_law="SIEVERT")

# Explicit x-locations of material interfaces (left -> right)
x_fluid_left = r_min
x_fluid_right = r_max
x_purge_wall_right = r_max + wall_p_r_thickness
x_purge_right = x_purge_wall_right + height_pc

left_wall = F.VolumeSubdomain1D(
    id=left_wall_id,
    material=Eurofer,#Dummy2 Eurofer
    borders=[0,x_fluid_left],
)

fluid_subdomain = F.VolumeSubdomain1D(
    id=fluid_id,
    material=He_Purge, #Dummy1 He_Purge PbLi
    borders= [x_fluid_left, x_fluid_right],
)

right_p_wall = F.VolumeSubdomain1D(
    id=right_p_wall_id,
    material=Eurofer, #Dummy2 Eurofer
    borders=[x_fluid_right, x_purge_wall_right]
)

purge_subdomain = F.VolumeSubdomain1D(
    id=purge_id,
    material=He_Purge, #Eurofer He_Purge Dummy1 PbLi
    borders=[x_purge_wall_right, x_purge_right]
)

right_wall = F.VolumeSubdomain1D(
    id=right_wall_id,
    material=Eurofer, #Dummy2 Eurofer
    borders=[x_purge_right, x_purge_right+wall_r_thickness],
)

my_sim.subdomains = [
    outer_wall_boundary,
    inner_wall_boundary,
    left_wall,
    fluid_subdomain,
    right_p_wall,
    purge_subdomain,
    right_wall,
]

pt = 1e20

my_sim.interfaces = [
    F.Interface(interface_l_id, (fluid_subdomain, left_wall), penalty_term=pt),
    F.Interface(interface_r_id, (fluid_subdomain, right_p_wall), penalty_term=pt),
    F.Interface(99, (right_p_wall,purge_subdomain), penalty_term=pt),
    F.Interface(98,(purge_subdomain,right_wall),penalty_term=pt)
]

my_sim.surface_to_volume = {
    inner_wall_boundary: left_wall,
    outer_wall_boundary: right_wall,
}

H = F.Species("H")
my_sim.species = [H]
H.subdomains = [fluid_subdomain, left_wall, right_p_wall, purge_subdomain,right_wall]

my_sim.temperature = 800

my_sim.boundary_conditions = [
    F.FixedConcentrationBC(species=H, subdomain=inner_wall_boundary, value=1e20),
    F.FixedConcentrationBC(species=H, subdomain=outer_wall_boundary, value=0),
]

my_sim.settings = F.Settings(
    transient=False,
    atol=1e10,
    rtol=1e-8,
    max_iterations=600,
)

my_sim.initialise()
my_sim.run()

import numpy as np
import matplotlib.pyplot as plt

plot_domains = [
    (left_wall, "left_wall"),
    (fluid_subdomain, "fluid"),
    (right_p_wall, "right_p_wall"),
    (purge_subdomain, "purge"),
    (right_wall,"right_wall")
]

plt.figure(figsize=(9, 4))
for volume, label in plot_domains:
    u = H.subdomain_to_post_processing_solution[volume]
    x = u.function_space.tabulate_dof_coordinates()[:, 0]
    c = u.x.array.real.copy()
    order = np.argsort(x)
    plt.plot(x[order], c[order], label=label)

for x_interface in [x_fluid_left, x_fluid_right, x_purge_wall_right]:
    plt.axvline(x_interface, color="k", linestyle="--", linewidth=0.8)

plt.xlabel("x [m]")
plt.ylabel("H concentration")
plt.title("1D discontinuous solution")
plt.legend()
plt.tight_layout()
plt.show()

```

Best,  
Chris

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [February 11, 2026, 4:19pm UTC](https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158/5 "2026-02-11T16:19:36Z")

</div>

I’m not sure simualting the He purge like this is the way to go. Plus, if you’re mixing solubility laws, you should use the `"penalty"` interface method and not the default “nitsche”.

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [April 13, 2026, 12:46pm UTC](https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158/6 "2026-04-13T12:46:31Z")

</div>

@ChrisSo we moved the discussion to this issue

> <https://github.com/festim-dev/FESTIM/issues/1092>
>
> Hey,
> 
> I would like to be able to choose a rate-limited interface condition in ad…dition to the current equilibrium interface conditions based on Henry's or Sieverts' law.
> 
> At the moment, recombination and dissociation can already be modelled through boundary conditions, but in my case the gas domain is solved explicitly. Therefore, I would need an analogous formulation as an internal interface condition between two solved domains rather than as a boundary condition.
> 
> This would be especially useful for gas/solid transport problems where the transfer across the interface is governed by dissociation and recombination kinetics.
> 
> It would also be important that such an interface condition can account for multi-species transport, for example a gas containing H2, HT, and T2 adjacent to a metallic wall containing only H and T.
> 
> A possible form of the interface flux for tritium would be:
> 
> \`\`\`math
> \- D\_{s,T} \\nabla c\_{s,T} \\cdot n\_s
> =
> 2 \\left( k\_B T K\_{d,T} c\_{g,T\_2} - K\_{r,T} c\_{s,T}^2 \\right)
> +
> \\left( k\_B T K\_{d,HT} c\_{g,HT} - K\_{r,HT} c\_{s,T} c\_{s,H} \\right)
> \`\`\`
> 
> For a single-species case, \`c\_{s,H}\` and \`c\_{g,HT}\` would be zero.
> 
> A reference for this type of rate-limited gas/solid interface model is \[Pasler et al.\](https://doi.org/10.3390/app11083481), who implemented dissociation/recombination-controlled interface conditions for hydrogen isotope transport in OpenFOAM.
> 
> This would be particularly useful for applications such as HCPB tritium transport, where surface-limited behaviour can have a significant influence on permeation and inventory.