Coupled simulation - application of tolerances

[griff]

Hi all,

I am running simple a coupled simulation for heat transfer and hydrogen transport in a ceramic pebble using fairly modest heat flux and tritium source rate values. I am seeing some non-physical results which vary depending on the tolerances - see first screen shot.

The sphere is represented as a 2D circular cross-section in Cartesian coordinates.

Vacuum BCs are applied to each hemisphere of the pebble for hydrogen transport.

Fixed temperature BCs are applied to each hemisphere for heat transfer (~room temp).

When I ‘loosen’ the tolerances, tritium concentrations show non-zero values but peak in just before the boundary edge – almost like it builds up before leaving the pebble surface (strange). See second screen shot.

fyi I have also tried:

• Using a dolphinx rectangular mesh

• Applying HeatFluxBCs to the surfaces (for convective cooling) instead of FixedTemperature ones for heat transfer.

• Applying SurfaceReactionBCs to the surfaces for tritium transport

• Using a finer (and coarser) mesh.

• Changing the material properties to Tungsten, see third screen shot

Would really appreciate any help and advice - thanks!

import festim as F
import numpy as np
import matplotlib.pyplot as plt
import dolfinx
import dolfinx.io
from mpi4py import MPI  
import os
import gmsh

if not os.path.exists('results'):
    os.makedirs('results')


pebble_material = F.Material(
    name="LiTiO3",
    D_0=9.5E-6,               # Diffusion pre-exponential factor [m²/s] 9.5E-6
    E_D=0.80,                 # Activation energy [eV] 0.80
    thermal_conductivity=2.5,  # Thermal conductivity at 300K [W/m/K]
    density=2881,              # Density [kg/m³]
    heat_capacity=1600         # Heat capacity [J/kg/K]
)

def create_circular_cross_section():
    """
    Create 2D circular disk with properly centered left/right boundaries
    """
    gmsh.initialize()
    gmsh.model.add("circular_pebble_boundaries")

    radius = 0.5e-3 # 0.5 mm radius (1 mm diameter pebble)

    # Create circle with more control points for smoother boundary
    center = gmsh.model.geo.addPoint(0.0, 0.0, 0.0)

    # Create more points around the circle for smoother boundary
    n_points = 10  # Use even number for symmetric division
    points = []
    arcs = []

    for i in range(n_points):
        angle = 2 * np.pi * i / n_points
        x = radius * np.cos(angle)
        y = radius * np.sin(angle)
        point = gmsh.model.geo.addPoint(x, y, 0.0)
        points.append(point)

    # Create arcs between consecutive points
    for i in range(n_points):
        next_i = (i + 1) % n_points
        arc = gmsh.model.geo.addCircleArc(points[i], center, points[next_i])
        arcs.append(arc)

    # Left boundary: from top (index 4) to bottom (index 12) going through left side
    left_start = n_points // 4      # Index 4 (top)
    left_end = 3 * n_points // 4    # Index 12 (bottom)
    left_arcs = arcs[left_start:left_end]

    # Right boundary: from bottom (index 12) to top (index 4) going through right side
    right_arcs = arcs[left_end:] + arcs[:left_start]

    # Create surface
    curve_loop = gmsh.model.geo.addCurveLoop(arcs)
    surface = gmsh.model.geo.addPlaneSurface([curve_loop])

    gmsh.model.geo.synchronize()

    # Physical groups
    gmsh.model.addPhysicalGroup(1, left_arcs, 1, "left_boundary")
    gmsh.model.addPhysicalGroup(1, right_arcs, 2, "right_boundary")
    gmsh.model.addPhysicalGroup(2, [surface], 3, "volume")

    # IMPROVED mesh settings
    char_length = radius / 50 # Finer mesh
    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), char_length)

    # Set mesh algorithm for better quality
    gmsh.option.setNumber("Mesh.Algorithm", 6)  
    gmsh.option.setNumber("Mesh.RecombineAll", 1)  
    gmsh.option.setNumber("Mesh.Smoothing", 5)  

    gmsh.model.mesh.generate(2)

    mesh_2d, cell_tags, facet_tags = dolfinx.io.gmshio.model_to_mesh(
        gmsh.model, comm=MPI.COMM_WORLD, rank=0, gdim=2
    )

    gmsh.finalize()
    return mesh_2d, cell_tags, facet_tags, radius

# FESTIM setup with boundaries
mesh_2d, cell_tags, facet_tags, radius = create_circular_cross_section()

festim_mesh = F.Mesh(mesh_2d)

# Volume subdomain
vol_sphere = F.VolumeSubdomain(id=3, material=pebble_material)

# CORRECTED boundary locator functions
def left_boundary_locator(x):
    """Locate left boundary (neutron-facing side) - properly centered"""
    r_coord = x[0]  # x-coordinate
    z_coord = x[1]  # y-coordinate
    distance_from_center = np.sqrt(r_coord**2 + z_coord**2)
    # Left side of circle (x < 0) and on boundary
    return np.logical_and(
        np.isclose(distance_from_center, radius, atol=1e-10),
        r_coord < 0.0 
    )

def right_boundary_locator(x):
    """Locate right boundary (coolant side) - properly centered"""
    r_coord = x[0]  # x-coordinate
    z_coord = x[1]  # y-coordinate
    distance_from_center = np.sqrt(r_coord**2 + z_coord**2)
    # Right side of circle (x > 0) and on boundary
    return np.logical_and(
        np.isclose(distance_from_center, radius, atol=1e-10),
        r_coord > 0.0  
    )

# Surface subdomains
surf_left = F.SurfaceSubdomain(id=1, locator=left_boundary_locator)      
surf_right = F.SurfaceSubdomain(id=2, locator=right_boundary_locator)    

subdomains = [vol_sphere, surf_left, surf_right]


# 1. CREATE HYDROGEN TRANSPORT PROBLEM
hydrogen_problem = F.HydrogenTransportProblem()
hydrogen_problem.mesh = festim_mesh
hydrogen_problem.subdomains = subdomains


mobile_T = F.Species("T")
hydrogen_problem.species = [mobile_T]

# Tritium source
tritium_source = F.ParticleSource(
    value=1e14,  # Source rate [atoms/m³/s]
    volume=vol_sphere,
    species=mobile_T
)
hydrogen_problem.sources = [tritium_source]


# Left boundary: 
left_bc = F.DirichletBC(
    subdomain=surf_left,
    species=mobile_T,
    value=0  # Vacuum 
)

# Right boundary: 
right_bc = F.DirichletBC(
    subdomain=surf_right,
    value=0,      # Vacuum
    species=mobile_T
)

hydrogen_problem.boundary_conditions = [left_bc, right_bc]

# 2. CREATE HEAT TRANSFER PROBLEM
heat_problem = F.HeatTransferProblem()
heat_problem.mesh = festim_mesh
heat_problem.subdomains = subdomains
heat_value = 1e3 # Volumetric heat source [W/m³]

heat_source = F.HeatSource(
    value=heat_value,  # Volumetric heat source [W/m³]
    volume=vol_sphere
)

heat_problem.sources = [heat_source]
chamber_temperature = 350  # K 

# initital conditions
heat_problem.initial_conditions = [
    F.InitialTemperature(
        value=chamber_temperature,  # Initial temperature [K]
        volume=vol_sphere
    )
]

# Make sure the boundary condition values make sense
left_bc_T = F.FixedTemperatureBC(
    subdomain=surf_left,
    value=chamber_temperature  # Fixed temperature [K] at boundary
)

right_bc_T = F.FixedTemperatureBC(
    subdomain=surf_right,
    value=chamber_temperature  # Fixed temperature [K] at boundary
) 

heat_problem.boundary_conditions = [left_bc_T, right_bc_T]

# simulation settings
# Run time:
run_time = 50

hydrogen_problem.settings = F.Settings(
    atol=1e12,
    rtol=1e-6,
    final_time=run_time,  
    transient=True  
)
hydrogen_problem.settings.stepsize = F.Stepsize(1)  # 1 second time step

heat_problem.settings = F.Settings(
    atol=1e-6,      
    rtol=1e-10,      
    final_time=run_time,
    transient=True  
)
heat_problem.settings.stepsize = F.Stepsize(1)  # 1 second time step

print("Creating coupled problem")

# 3. CREATE COUPLED PROBLEM
model = F.CoupledTransientHeatTransferHydrogenTransport(
    heat_problem=heat_problem,
    hydrogen_problem=hydrogen_problem
)

# Exports for hydrogen problem
hydrogen_problem.exports = [
    F.VTXSpeciesExport(
        filename="results/2D/coupled/tritium_transport/pebble_tritium.bp",
        field=[mobile_T],
        checkpoint=True,
    ),
    F.XDMFExport(
        filename="results/2D/coupled/tritium_transport/pebble_tritium.xdmf",
        field=[mobile_T]
    )
]

# Export temperature data
heat_problem.exports = [
    F.VTXTemperatureExport(
        filename="results/2D/coupled/heat_transfer/pebble_temperature.bp",
    ),
]

print("Individual problem exports configured")

# Run the coupled problem
print("Initializing coupled heat transfer and hydrogen transport simulation...")
model.initialise()
print("Running coupled simulation...")
model.run()
print("Coupled simulation completed!")

# Post-process: Read final concentration field
final_concentration = F.read_function_from_file(
    filename="results/2D/coupled/tritium_transport/pebble_tritium.bp",
    name="T",
    timestamp=hydrogen_problem.settings.final_time,
)

def plot_coupled_results():
    """
    Plot both temperature and tritium concentration from coupled simulation
    """
    # Get temperature and concentration from coupled model
    temperature_field = model.heat_problem.u  # Temperature solution
    concentration_field = model.hydrogen_problem.u  # Concentration solution

    V = concentration_field.function_space
    coords = V.tabulate_dof_coordinates()
    x_coords = coords[:, 0]  # X coordinates
    y_coords = coords[:, 1]  # Y coordinates
    conc_values = concentration_field.x.array
    temp_values = temperature_field.x.array

    # Create 2x2 subplot layout
    fig, axes = plt.subplots(2, 2, figsize=(15, 12))

    # Plot 1: Temperature distribution
    scatter1 = axes[0,0].scatter(x_coords*1000, y_coords*1000, c=temp_values,
                                cmap='plasma', s=20, alpha=0.8)
    axes[0,0].set_xlabel('X Position (mm)')
    axes[0,0].set_ylabel('Y Position (mm)')
    axes[0,0].set_title('Temperature Distribution')
    axes[0,0].set_aspect('equal')
    axes[0,0].grid(True, alpha=0.3)
    plt.colorbar(scatter1, ax=axes[0,0], label='Temperature (K)')

    # Plot 2: Tritium concentration distribution
    scatter2 = axes[0,1].scatter(x_coords*1000, y_coords*1000, c=conc_values,
                                cmap='viridis', s=20, alpha=0.8)
    axes[0,1].set_xlabel('X Position (mm)')
    axes[0,1].set_ylabel('Y Position (mm)')
    axes[0,1].set_title('Tritium Concentration Distribution')
    axes[0,1].set_aspect('equal')
    axes[0,1].grid(True, alpha=0.3)
    plt.colorbar(scatter2, ax=axes[0,1], label='Concentration (atoms/m³)')

    # Plot 3: Temperature profile (horizontal slice through center)
    y_tolerance = radius * 0.05
    center_mask = np.abs(y_coords) < y_tolerance
    if np.any(center_mask):
        x_center = x_coords[center_mask]
        temp_center = temp_values[center_mask]
        sort_idx = np.argsort(x_center)

        axes[1,0].plot(x_center[sort_idx]*1000, temp_center[sort_idx], 'r-o', linewidth=2, markersize=4)
        axes[1,0].set_xlabel('X Position (mm)')
        axes[1,0].set_ylabel('Temperature (K)')
        axes[1,0].set_title('Temperature Profile: Left → Right')
        axes[1,0].grid(True, alpha=0.3)
        axes[1,0].legend()

    # Plot 4: Concentration profile (horizontal slice through center)
    if np.any(center_mask):
        x_center = x_coords[center_mask]
        conc_center = conc_values[center_mask]
        sort_idx = np.argsort(x_center)

        axes[1,1].plot(x_center[sort_idx]*1000, conc_center[sort_idx], 'b-o', linewidth=2, markersize=4)
        axes[1,1].set_xlabel('X Position (mm)')
        axes[1,1].set_ylabel('Tritium Concentration (atoms/m³)')
        axes[1,1].set_title('Concentration Profile: Left → Right')
        axes[1,1].grid(True, alpha=0.3)
        axes[1,1].legend()

    plt.tight_layout()

    # Create output directory
    os.makedirs('results/2D/coupled', exist_ok=True)
    plt.savefig('results/2D/coupled/coupled_results_analysis.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Print summary
    print(f"\n=== COUPLED SIMULATION RESULTS ===")
    print(f"Temperature range: {np.min(temp_values):.1f} - {np.max(temp_values):.1f} K")
    print(f"Temperature gradient: {np.max(temp_values) - np.min(temp_values):.1f} K")
    print(f"Concentration range: {np.min(conc_values):.2e} - {np.max(conc_values):.2e} atoms/m³")
    print(f"Average temperature: {np.mean(temp_values):.1f} K")
    print(f"Average concentration: {np.mean(conc_values):.2e} atoms/m³")
    print(f"Mesh points: {len(conc_values)}")
    print(f"Simulation time: {run_time} seconds")

plot_coupled_results()

[remidm]

Hi @griff thanks for posting on the discourse!

This is a typical finite element mesh issue called the Gibbs Phenomenon. In a nutshell, you have high gradients at the boundary of your domain, and to be able to approximate these gradients you need local refinement. That’s the reason why you don’t see this on the heat transfer problem: gradients are smoother and easier to approximate.

To convince you further, try increasing (a lot) the diffusivity of tritium in the pebble and these overshoots will disappear.

It’s better explained in the dolfinx tutorial.

Solutions:

• increase the diffusivity: not possible if this is a fixed material property obviously

• local refinement: can be costly in a 2D geometry (even more in 3D)

Related to the last bullet point, you could simplify your problem by doing a 1D spherical/cylindrical geometry, that way you can massively refine the mesh near the boundary. See Non-uniform 1D meshes

[remidm]

Something else to consider: do you need a transient heat transfer simulation? I think it could be a simple steady state given the size of the pebble

[griff]

Hi @remidm,

Thanks for coming back to me.

I tried increasing the diffusivity and as you say, the issue is resolved. Sadly, local refinement in 2D does not seem to help (see first set of figs below).

In 1D the situation is similar even when using a refined mesh (see below script and second set of figs).

import festim as F
import numpy as np
import matplotlib.pyplot as plt
import os
if not os.path.exists('results'):
    os.makedirs('results')

# Create 1D mesh for pebble with cartesian coordinates
pebble_length = 0.5e-3 
n_elements = 10000

def create_exponential_refined_mesh(length, n_elements, grading_factor=2.0):
    """
    Create mesh with exponential grading toward boundaries
    grading_factor: higher values give more aggressive refinement
    """
    # Create symmetric grading from center to boundaries
    half_n = n_elements // 2
    
    # Generate graded spacing from center (0.5) to boundary (0)
    eta = np.linspace(0, 1, half_n + 1)
    spacing_ratio = np.exp(-grading_factor * eta)
    spacing_ratio = spacing_ratio / np.sum(spacing_ratio[:-1])  # Normalize
    
    # Create positions for left half
    left_positions = np.zeros(half_n + 1)
    for i in range(1, half_n + 1):
        left_positions[i] = left_positions[i-1] + spacing_ratio[i-1] * length/2
    
    # Mirror for right half
    right_positions = length - left_positions[::-1]
    
    # Combine and remove center duplicate
    vertices = np.concatenate([left_positions[:-1], right_positions])
    return vertices


refined_vertices = create_exponential_refined_mesh(pebble_length, n_elements)
mesh = F.Mesh1D(refined_vertices)

run_time = 10

# Define material properties
pebble_material = F.Material(
    name="LiTiO3",
    D_0=9.5E-6,           # Diffusion pre-exponential factor [m²/s] 9.5E-6
    E_D=0.80,              # Activation energy [eV] 0.80
    thermal_conductivity=2.5,  # Thermal conductivity at 300K [W/m/K]
    density=2881,              # Density [kg/m³]
    heat_capacity=1600         # Heat capacity [J/kg/K]
)

# Define subdomains
vol = F.VolumeSubdomain1D(id=1, borders=[0, pebble_length], material=pebble_material)
surf_left = F.SurfaceSubdomain1D(id=1, x=0)
surf_right = F.SurfaceSubdomain1D(id=2, x=pebble_length)
subdomains = [vol, surf_left, surf_right]

# 1. CREATE HYDROGEN TRANSPORT PROBLEM
hydrogen_problem = F.HydrogenTransportProblem()
hydrogen_problem.mesh = mesh
hydrogen_problem.subdomains = subdomains


mobile_T = F.Species("T")
hydrogen_problem.species = [mobile_T]

# Tritium source
tritium_source = F.ParticleSource(
    value=5.22e13 / run_time,  # Source rate [atoms/m³/s]
    volume=vol,
    species=mobile_T
)
hydrogen_problem.sources = [tritium_source]

# surface reaction parameters for 1D geometry
k_r0 = 1e-2               # Forward rate pre-factor [m⁴/(atom·s)]
E_kr = 0.5                # Forward activation energy [eV] - for T2 formation
k_d0 = 1e-2               # Backward rate pre-factor [Pa·m⁴/(atom·s)]
E_kd = 0.8                # Backward activation energy [eV] - for T2 dissociation

def purge_gas_T_pressure(t=None):
    """
    Time-dependent purge gas tritium pressure for 2D geometry
    """
    if t is None:
        t = 0
    base_pressure = 1e-3  # 1 mPa base pressure
    buildup = 1e-4 * t    # Gradual tritium accumulation
    #return min(base_pressure + buildup, 1e-2)  # Cap at 10 mPa
    return 1e-2

# surface reactions
left_bc = F.SurfaceReactionBC(
    reactant=[mobile_T],
    gas_pressure=purge_gas_T_pressure,
    k_r0=k_r0,       
    E_kr=E_kr,
    k_d0=k_d0,
    E_kd=E_kd,
    subdomain=surf_left,
)

right_bc = F.SurfaceReactionBC(
    reactant=[mobile_T],
    gas_pressure=purge_gas_T_pressure,
    k_r0=k_r0,            
    E_kr=E_kr,      
    k_d0=k_d0,
    E_kd=E_kd,
    subdomain=surf_right,
)

hydrogen_problem.boundary_conditions = [left_bc, right_bc]

# 2. CREATE HEAT TRANSFER PROBLEM
heat_problem = F.HeatTransferProblem()
heat_problem.mesh = mesh
heat_problem.subdomains = subdomains
heat_value = 139 # Volumetric heat source [W/m³]

heat_source = F.HeatSource(
    value=heat_value,  
    volume=vol
)

heat_problem.sources = [heat_source]
chamber_temperature = 300  # K 

# initital conditions
heat_problem.initial_conditions = [
    F.InitialTemperature(
        value=chamber_temperature,  # Initial temperature [K]
        volume=vol
    )
]

# Make sure the boundary condition values make sense
left_bc_T = F.FixedTemperatureBC(
    subdomain=surf_left,
    value=chamber_temperature # Fixed temperature [K] at boundary
)

right_bc_T = F.FixedTemperatureBC(
    subdomain=surf_right,
    value=chamber_temperature  # Fixed temperature [K] at boundary
) 

heat_problem.boundary_conditions = [left_bc_T, right_bc_T]

# settings for both problems
hydrogen_problem.settings = F.Settings(
    atol=1e7,
    rtol=1e-8,
    final_time=run_time,  
    transient=True   
)
hydrogen_problem.settings.stepsize = F.Stepsize(1)  # 1 second time step

heat_problem.settings = F.Settings(
    atol=1,      
    rtol=1e-5,      
    final_time=run_time,
    transient=True
)
heat_problem.settings.stepsize = F.Stepsize(1)  # 1 second time step

print("Creating coupled problem")

# 3. CREATE COUPLED PROBLEM
model = F.CoupledTransientHeatTransferHydrogenTransport(
    heat_problem=heat_problem,
    hydrogen_problem=hydrogen_problem
)

print("Setting up exports on individual problems...")

# Export hydrogen/tritium data
hydrogen_problem.exports = [
    F.VTXSpeciesExport(
        filename="results/1D/heat_trans+h_transport/surf_reactions/pebble_tritium.bp",
        field=[mobile_T],
        checkpoint=True,
    ),
    F.XDMFExport(
        filename="results/1D/heat_trans+h_transport/surf_reactions/pebble_tritium.xdmf",
        field=[mobile_T]
    )
]

# Export temperature data
heat_problem.exports = [
    F.VTXTemperatureExport(
        filename="results/1D/heat_trans+h_transport/surf_reactions/pebble_temperature.bp",
    )
]
print("Individual problem exports configured")

# Run the coupled problem
print("Initializing coupled heat transfer and hydrogen transport simulation...")
model.initialise()
print("Running coupled simulation...")
model.run()
print("Coupled simulation completed!")

# Post-process: Read final concentration field
final_concentration = F.read_function_from_file(
    filename="results/1D/heat_trans+h_transport/surf_reactions/pebble_tritium.bp",
    name="T",
    timestamp=hydrogen_problem.settings.final_time,
)

def analyze_and_plot_results(model, pebble_material, heat_value, run_time, pebble_length):
    """
    Simplified analysis without format string issues
    """
    print("Extracting data directly from model...")
    
    # Get concentration data from hydrogen problem
    final_conc_values = model.hydrogen_problem.u.x.array
    # Get temperature data from heat problem
    final_temp_values = model.heat_problem.u.x.array

    print(f"Successfully extracted data:")
    print(f"  Concentration range: {np.min(final_conc_values):.2e} to {np.max(final_conc_values):.2e} atoms/m³")
    print(f"  Temperature range: {np.min(final_temp_values):.1f} to {np.max(final_temp_values):.1f} K")

    x_coords = np.linspace(0, pebble_length, len(final_conc_values))
    
    plt.figure(figsize=(12, 6))
    
    # Plot 1: Temperature Profile
    plt.subplot(1, 2, 1)
    plt.plot(x_coords*1000, final_temp_values, 'r-', linewidth=3, label='Temperature')
    plt.xlabel('Position (mm)')
    plt.ylabel('Temperature (K)')
    plt.title('Temperature Profile')
    plt.grid(True, alpha=0.3)
    plt.legend()

    # Plot 2: Concentration Profile
    plt.subplot(1, 2, 2)
    plt.plot(x_coords*1000, final_conc_values, 'b-', linewidth=3, label='Tritium Concentration')
    plt.xlabel('Position (mm)')
    plt.ylabel('Tritium Concentration (atoms/m³)')
    plt.title('Tritium Concentration Profile')
    plt.grid(True, alpha=0.3)
    plt.legend()

    plt.tight_layout()
    
    # Save plot
    save_dir = 'results/1D/heat_trans+h_transport/surf_reactions/'
    os.makedirs(save_dir, exist_ok=True)
    
    save_path = os.path.join(save_dir, 'coupled_analysis_refined_mesh.png')
    plt.savefig(save_path, dpi=300, bbox_inches='tight')
    print(f"Plot saved to: {save_path}")

    plt.show()

    return final_temp_values, final_conc_values, x_coords

final_temp_values, final_conc_values, x_coords = analyze_and_plot_results(
        model, pebble_material, heat_value, run_time, pebble_length
    )

print("Coupled analysis complete!")

On your point re the need for a transient simulation, my understanding is that festim-2 problems must set this way in the settings, I get an error if it transient = False

raise TypeError(
    "Both the heat and hydrogen problems must be set to transient"
)

[remidm]

Hi Tom, a few things to consider:

• your problem seems extremely stiff. Diffusivity evolves as the inverse exponential of temperature and it seems like on your 1D case the temperature is ~300 K everywhere (that’s room temperature).
• for simplicity I replaced the temperature in your model by a homogeneous temperature
• I also simplified the local refinement with

zone_to_refine = 0.001e-2

refined_vertices = np.concatenate(
    [
        np.linspace(0, zone_to_refine, num=3000),
        np.linspace(zone_to_refine, pebble_length - zone_to_refine, num=1000),
        np.linspace(pebble_length - zone_to_refine, pebble_length, num=3000),
    ]
)
mesh = F.Mesh1D(refined_vertices)

• with T = 310 K and this mesh, I got a smoother field. But again, this is an extremely stiff problem.

If the temperature increases, diffusivity will increase, and fields will get smoother

Regarding the coupling:

• set the heat transfer problem to steady state
• solve the heat transfer problem (.initialise() then .run())
• set up your H transport problem and set hydrogen_problem.temperature = heat_problem.u
• solve the H transport problem