Mass balancing in FESTIM

User Support
1

Hello everyone,

I’m new to FESTIM and currently preparing some Tritium transport modelling work. Before starting, I discussed good practices with colleagues experienced in hydrogen transport (mostly using OpenFOAM), and they recommended regularly checking mass conservation. I wanted to implement this in FESTIM to build confidence in my setup.

To check mass balance for steady-state cases, I compared the hydrogen fluxes across all boundaries (Inflow - Outflow). When I used a volumetric source term, I added the generated hydrogen to the inflow side. I obtained the boundary fluxes using HydrogenFlux via DerivedQuantities.

In several cases, my mass-balance check suggests that the total fluxes do not match. Since I’m still new to FESTIM, I’m unsure whether my way of calculating the balance is complete, whether I am missing a detail in how fluxes should be interpreted, or whether there are limitations in how fluxes are currently exported. I would really appreciate some guidance from more experienced users.

Here is how I computed the relevant quantities:

# Before running the simulation

Flux_wall_l = F.HydrogenFlux(surface=wall_l_id)
Flux_wall_r = F.HydrogenFlux(surface=wall_r_id)
Flux_outlet = F.HydrogenFlux(surface=outlet_id)
Flux_inlet = F.HydrogenFlux(surface=inlet_id)

derived_quantities = F.DerivedQuantities(
    [Flux_wall_l, Flux_wall_r, Flux_outlet, Flux_inlet],
    show_units=True
)

my_sim.exports = [derived_quantities]

# After running the simulation

Difference = (Flux_wall_r.data[0]
              + Flux_outlet.data[0]
              + Flux_inlet.data[0]
              + Flux_wall_l.data[0]
              + Source)

Relative = Difference / (Flux_inlet.data[0] + Source)

print(Difference)
print(f'Relative difference = {Relative:.2%}')

For simple 2D Cartesian cases with Dirichlet-zero walls, flux conditions at the inlet, and/or a source term, the mass balance looks consistent.

However, when I switch to a cylindrical (axisymmetric) mesh and remove the Dirichlet BC on the wall representing the symmetry axis, I see a significant discrepancy. Just to clarify my interpretation: in my mesh, wall_l corresponds to the axis, so there should be no net mass transfer across it. Anything crossing the axis should remain inside the domain. Even with that understanding, I still observe that the inflow is ~50% higher than the outflow in these axisymmetric cases.

When I add an advection term to the transport equation, the discrepancy becomes even larger and appears to scale with the advection velocity (higher velocity → larger mismatch). This makes me suspect that I may be overlooking some detail about how fluxes should be combined or interpreted in FESTIM, especially in combination with advection.

I would be very grateful if someone could help me understand whether I am performing the mass balance correctly, or if there is something important I might be missing in this kind of setup.

Below is a minimal script that reproduces the behaviour I described.

Thanks a lot in advance for your help,
Chris

import festim as F
import numpy as np
from fenics import UnitSquareMesh, CompiledSubDomain, MeshFunction, plot, RectangleMesh, Point
import matplotlib.pyplot as plt


def create_mesh2D_rec(inlet_id,outlet_id,wall_left_id,wall_right_id, r_max, z_max, r_min=0, z_min=0,nx=40, ny=200):
    """Create a mesh for an annular flow mass exhanger
    
    Args: 
        inlet_id
        outlet_id
        wall_id
        axis_id
        r_max
        z_max
        r_min
        z_min        
    """
    
    # Create a mesh with fenix
    mesh_fenics = RectangleMesh(Point(r_min, z_min), Point(r_max, z_max), nx, ny)
    # --- Ensure topology is initialized ---
    mesh_fenics.init(mesh_fenics.topology().dim() - 1)
    # marking the physical groups (volumes and surfaces)
    volume_markers = MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim())
    volume_markers.set_all(1)
    tol = 1e-14
    inlet_surface = CompiledSubDomain('on_boundary && near(x[1],z_min,tol)',tol=tol, z_min=z_min)
    outlet_surface = CompiledSubDomain('on_boundary && near(x[1], z_max , tol)', tol=tol, z_max = z_max)
    wall_right_surface = CompiledSubDomain('on_boundary && near(x[0], r_max, tol)', tol=tol, r_max = r_max)
    wall_left_surface = CompiledSubDomain('on_boundary && near(x[0], r_min, tol)', tol=tol, r_min = r_min)
    
    surface_markers = MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim() - 1)
    surface_markers.set_all(0)
    
    inlet_surface.mark(surface_markers, inlet_id)
    outlet_surface.mark(surface_markers, outlet_id)
    wall_right_surface.mark(surface_markers, wall_right_id)
    wall_left_surface.mark(surface_markers, wall_left_id)
    plot(mesh_fenics)
    return mesh_fenics, volume_markers, surface_markers


inner_radius = 0
outer_radius = 1
length = 1
inlet_id = 1
outlet_id = 2
wall_r_id = 3
wall_l_id = 4
my_sim = F.Simulation()
mesh, volume_markers, surface_markers = create_mesh2D_rec(inlet_id=inlet_id,outlet_id=outlet_id
                                                      ,wall_right_id=wall_r_id,wall_left_id = wall_l_id,
                                                      r_max=outer_radius, z_max=length,
                                                      nx=500,ny=500)


my_sim.mesh = F.Mesh(
    mesh = mesh,
    volume_markers=volume_markers,
    surface_markers=surface_markers,
    #type="cylindrical",
)
my_sim.materials = F.Material(id=1, D_0=1, E_D=0)
my_sim.T = 500
my_sim.boundary_conditions = [
    F.DirichletBC(field=0, surfaces=wall_r_id, value=0),
    F.DirichletBC(field=0, surfaces=wall_l_id, value=0),
    F.FluxBC(field=0, surfaces=inlet_id, value=1e17)
]
my_sim.settings = F.Settings(
    absolute_tolerance=1e-10,
    relative_tolerance=1e-10,
    transient=False
)
Flux_wall_l = F.HydrogenFlux(surface=wall_l_id)
Flux_wall_r = F.HydrogenFlux(surface=wall_r_id)
Flux_outlet = F.HydrogenFlux(surface=outlet_id)
Flux_inlet = F.HydrogenFlux(surface=inlet_id)
derived_quantities = F.DerivedQuantities(
                                        [Flux_wall_l,Flux_wall_r,Flux_outlet,Flux_inlet],
                                        show_units=True)
my_sim.exports = [derived_quantities]


my_sim.initialise()
my_sim.run()


hydrogen_concentration = my_sim.h_transport_problem.mobile.solution
plt.figure()
cs=plot(hydrogen_concentration)
plt.colorbar(cs)
plt.show()


Difference_sym= Flux_wall_r.data[0] + Flux_outlet.data[0]  + Flux_inlet.data[0] + Flux_wall_l.data[0]
Relative_sym = Difference_sym/Flux_inlet.data[0]
print(Difference_sym)
print(f'Relative difference = {Relative_sym:.2%}')
print(Flux_inlet.data[0])
print(Flux_wall_r.data[0])
print(Flux_wall_l.data[0])
print(Flux_outlet.data[0])


from fenics import interpolate, Expression, VectorFunctionSpace
functionspace = VectorFunctionSpace(my_sim.mesh.mesh, "CG", 1)
velocity = Expression(("0", "4*v/pow(R,2)*x[0]*(R-x[0])"),v=20,R=outer_radius,degree=2)
velocity = interpolate(velocity, functionspace)
CS = plot(velocity)
plt.colorbar(CS)
plt.show()


Flux_wall_l = F.HydrogenFlux(surface=wall_l_id)
Flux_wall_r = F.HydrogenFlux(surface=wall_r_id)
Flux_outlet = F.HydrogenFlux(surface=outlet_id)
Flux_inlet = F.HydrogenFlux(surface=inlet_id)
derived_quantities = F.DerivedQuantities(
                                        [Flux_wall_l,Flux_wall_r,Flux_outlet,Flux_inlet],
                                        show_units=True)
my_sim.exports = [derived_quantities]


from fenics import inner, dot, grad
my_sim.initialise()
hydrogen_concentration = my_sim.h_transport_problem.mobile.solution
test_function_solute = my_sim.h_transport_problem.mobile.test_function
advection_term = inner(dot(grad(hydrogen_concentration),velocity),test_function_solute)*my_sim.mesh.dx
my_sim.h_transport_problem.F += advection_term
my_sim.run()
ch = plot(hydrogen_concentration)
plt.colorbar(ch)
plt.show()


Difference_sym= Flux_wall_r.data[0] + Flux_outlet.data[0]  + Flux_inlet.data[0] + Flux_wall_l.data[0]
Relative_sym = Difference_sym/Flux_inlet.data[0]
print(Difference_sym)
print(f'Relative difference = {Relative_sym:.2%}')
print(Flux_inlet.data[0])
print(Flux_wall_r.data[0])
print(Flux_wall_l.data[0])
print(Flux_outlet.data[0])
2

Hi @ChrisSo and welcome to the FESTIM community!

The HydrogenFlux derived quantity class is only suitable for cartesian coordinates. A warning statement should be thrown when using it in cylindrical or spherical systems. So that’s a first source of error. Instead, consider using SurfaceFlux.

Regarding advection-diffusion problem, this is a known issue with Galerkin methods. See this comment on the FEniCS discourse. Depending on your exact problem, you may have to use some stabilisation methods.

Finally, you may want to look into FESTIM2 instead of FESTIM1!

3

Hi @remidm,

thank you very much for the quick response and the clarifications.

I have now repeated the calculations using SurfaceFlux as suggested. The results remain the same as with HydrogenFlux, but I do now see the warning message related to cylindrical meshes, which I had previously overlooked. That is helpful to know.

Thank you as well for pointing me to the advection discussion on the FEniCS forum. I understand the issue better now. I was wondering whether these advection-related limitations are handled differently in FESTIM2, is mass conserved in advection cases without additionally modifying the code?

As this project is for my master thesis on component-level tritium transport modelling, I currently assume that I should stick with the released version(FESTIM1). Since FESTIM2 is (as far as I can tell) not completely finished or validated, I am unsure whether I can justify its use in my thesis before the professor without doing a decent amount of validation myself. Is there an approximate timeline for when FESTIM2 might be released or considered validated enough for research work?

Thanks again for the support and for pointing me in the right direction.
Best,
Chris

4

Hi Chris,

My apologies, I meant SurfaceFluxCylindrical.

Regarding the mass conservation in advection-diffusion problems, it comes down to several things such as mesh Peclet number and stepsize.

If your mesh Peclet number is too high, then Galerkin methods will require some stabilisation. Which is something that can be easily implemented in FESTIM2 and we can share some code to do so (not sure if there’s already a tutorial on it @ckhurana768 ?)

If the stepsize is too big, simply reduce the stepsize! :smiley:

Regarding FESTIM2, it is currently released as a beta version and the V&V book has started to be updated for FESTIM2 (see this verification case for example).

If you want to look into advection-diffusion problems in multi-material geometries (solid-fluid) then I would recommend FESTIM2 because of new features such as mixed-formulation.

5

Hi Remi,

Thanks for your reply, and apologies for the late response. I took some time to better understand the problem and now have a much clearer picture.

Regarding the time step: I am currently running steady-state simulations, so the time discretisation should not play a role. Do you see any benefit in running the problem transient?

Regarding the Peclet number: my domain-wide Peclet number is around ~120 for the case I investigated, so advection-dominated. My understanding is that the problem I am facing is how the continuous Galerkin formulation behaves in such regimes. As in CG, continuity of the concentration field is enforced across element interfaces, but continuity of the advective flux is not enforced. For flux-driven transport this can lead to significant mass-balance errors, which seems to be what I am observing.

I discussed this with my supervisors, and given the scope and timeline of my work, they recommended that I stick with FESTIM1 for now, mainly because it is already validated and would not require additional benchmarking effort.

I did try adding stabilisation terms directly to h_transport_problem.F, but I was not able to improve the mass balance, yet. I am not sure whether FESTIM1 is still supported in this regard, but if you have any guidance or example on how stabilisation has been implemented successfully in FESTIM1, that would be extremely helpful.

In terms of application, the goal is to model a HCPB blanket. Ultimately this involves solid-fluid interfaces, but initially I am focusing on transport in the pebble bed itself and plan to add complexity step by step once that is working well.

Best regards,
Chris

6

@ChrisSo we are indeed focusing development and maintenance efforts on FESTIM2.

Implementing this stabilisation term in FESTIM1 would be very tricky because you would need to change the diffusion coefficient.

I don’t see any particular benefit in running it transient so sticking to steady state seems like a good idea.

Have you had any luck using SurfaceFluxCylindrical?

7

Hi Rémi,

Yes, using SurfaceFluxCylindrical resolved the mass balance issue for axisymmetric diffusion only cases.

I also implemented an axisymmetric advective surface flux by following the same pattern as SurfaceFluxCylindrical and adapting the Cartesian example from Task 16 of the tutorial. With this implementation, mass conservation is satisfied.
Below is the implementation in case it is useful for others facing the same issue:

class SurfaceAdvectionFluxCylindrical(F.SurfaceFluxCylindrical):
    """
    Object to compute the advective flux J_adv of a field c through a surface
    J_adv = ∫ (c * u) · n dS_cyl
    """
    def __init__(self, field, surface, velocity_field, azimuth_range=(0, 2 * np.pi)):
        super().__init__(field=field, surface=surface, azimuth_range=azimuth_range)
        self.velocity_field = velocity_field
        self.r = None

    @property
    def title(self):
        quantity_title = f"{self.field} advective flux surface {self.surface}"
        if self.show_units:
            return quantity_title + f" ({self.export_unit})"
        else:
            return quantity_title

    def compute(self):
        if self.r is None:
            mesh = self.function.function_space().mesh()
            rthetaz = f.SpatialCoordinate(mesh)
            self.r = rthetaz[0]

        flux = f.assemble(
            self.r
            * f.dot(self.velocity_field * self.function, self.n)
            * self.ds(self.surface)
        )

        flux *= self.azimuth_range[1] - self.azimuth_range[0]
        return flux

Using this approach, the advective contribution is now consistent with the diffusive surface flux in cylindrical coordinates and the overall mass balance checks out.

1 Like