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 andHydrogenTransportProblemDiscontinuousChangeVar()give very comparable results, whileHydrogenTransportProblemDiscontinuous()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 | |||
| Eurofer w. Helium | |||
| Dummy material w. Sieverts only |
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.
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