Grain-level 2D multi-material interface

[griff]

Hi All,

I am having an issue running a multi-material simulation in 2D. The model represents a square grain of ceramic (Li2TiO3) with a narrower grain adjacent to it that mimics an ‘amorphous’ layer. This layer has different material properties to the ceramic grain.

The tritium diffuses through the bulk, but when it reaches the interface with the amorphous layer it does not enter it. I wonder if it is an issue with the way I assign the F.Interface class?

Here is the code with a fig to illustrate:

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

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

hydrogen_problem = F.HydrogenTransportProblem()

width_bulk = 10e-6  # 5 micrometers
width_amorphous = 3e-6  # 3 micrometers

def create_two_region_mesh():
    """
    Create mesh using gmsh for better control over interfaces
    """
    
    grain_size = width_bulk
    gb_width = width_amorphous
    
    gmsh.initialize()
    gmsh.model.add("two_regions")
    
    # Create points
    p1 = gmsh.model.geo.addPoint(0, 0, 0)
    p2 = gmsh.model.geo.addPoint(grain_size, 0, 0)
    p3 = gmsh.model.geo.addPoint(grain_size, grain_size, 0)
    p4 = gmsh.model.geo.addPoint(0, grain_size, 0)
    
    p5 = gmsh.model.geo.addPoint(grain_size + gb_width, 0, 0)
    p6 = gmsh.model.geo.addPoint(grain_size + gb_width, grain_size, 0)
    
    # Create lines for grain region
    l1 = gmsh.model.geo.addLine(p1, p2)
    l2 = gmsh.model.geo.addLine(p2, p3)
    l3 = gmsh.model.geo.addLine(p3, p4)
    l4 = gmsh.model.geo.addLine(p4, p1)
    
    # Create lines for grain boundary region
    l5 = gmsh.model.geo.addLine(p2, p5)
    l6 = gmsh.model.geo.addLine(p5, p6)
    l7 = gmsh.model.geo.addLine(p6, p3)
    
    # Create curve loops
    grain_loop = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
    gb_loop = gmsh.model.geo.addCurveLoop([l5, l6, l7, -l2])
    
    # Create surfaces
    grain_surface = gmsh.model.geo.addPlaneSurface([grain_loop])
    gb_surface = gmsh.model.geo.addPlaneSurface([gb_loop])
    
    # Set mesh sizes
    gmsh.model.geo.mesh.setSize([(0, p1), (0, p2), (0, p3), (0, p4)], grain_size/ 200)  # Fine in grain
    gmsh.model.geo.mesh.setSize([(0, p5), (0, p6)], gb_width/ 200)  # Finer in GB
    
    # Physical groups
    gmsh.model.addPhysicalGroup(2, [grain_surface], 1)
    gmsh.model.setPhysicalName(2, 1, "grain")
    
    gmsh.model.addPhysicalGroup(2, [gb_surface], 2)
    gmsh.model.setPhysicalName(2, 2, "grain_boundary")
    
    # Boundaries
    gmsh.model.addPhysicalGroup(1, [l4], 3)  # left boundary
    gmsh.model.setPhysicalName(1, 3, "left")
    
    gmsh.model.addPhysicalGroup(1, [l6], 4)  # right boundary
    gmsh.model.setPhysicalName(1, 4, "right")
    
    gmsh.model.addPhysicalGroup(1, [l3, l7], 5)  # top boundary
    gmsh.model.setPhysicalName(1, 5, "top")
    
    gmsh.model.addPhysicalGroup(1, [l1, l5], 6)  # bottom boundary
    gmsh.model.setPhysicalName(1, 6, "bottom")
    
    gmsh.model.addPhysicalGroup(1, [l2], 7)  # interface
    gmsh.model.setPhysicalName(1, 7, "interface")
    
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(2)
    
    # Convert to dolfinx mesh

    mesh_2d, cell_tags, facet_tags = gmshio.model_to_mesh(
        gmsh.model, MPI.COMM_WORLD, 0, gdim=2
    )
    
    gmsh.finalize()
    
    total_width = grain_size + gb_width
    
    print(f"✓ GMSH two-region mesh created!")
    print(f"  - Total dimensions: {total_width*1e6:.1f} μm × {grain_size*1e6:.1f} μm")
    
    return mesh_2d, total_width, grain_size, gb_width, cell_tags, facet_tags

# Update region locators for the new geometry
def bulk_grain_locator(x):
    """Square grain region (left side)"""
    return x[0] < grain_size

def grain_boundary_locator(x):
    """Rectangle grain boundary region (right side)"""
    return x[0] >= grain_size

# Update boundary locators for new geometry
def left_boundary_locator(x):
    """Left edge of square grain"""
    return np.isclose(x[0], 0, atol=1e-12)

def right_boundary_locator(x):
    """Right edge of grain boundary"""
    return np.isclose(x[0], grain_size + gb_width, atol=1e-12)

def top_boundary_locator(x):
    """Top edge (both regions)"""
    return np.isclose(x[1], grain_size, atol=1e-12)

def bottom_boundary_locator(x):
    """Bottom edge (both regions)"""
    return np.isclose(x[1], 0, atol=1e-12)

def gb_interface_locator(x):
    """Interface between square grain and grain boundary"""
    return np.isclose(x[0], grain_size, atol=1e-12)


mesh_2d, total_width, grain_size, gb_width, cell_tags, facet_tags = create_two_region_mesh()

festim_mesh = F.Mesh(mesh_2d)
hydrogen_problem.mesh = festim_mesh

# Update materials (same as before)
bulk_material = F.Material(
    name="LiTiO3",
    D_0=9.5E-6,               
    E_D=0.80,                 
    thermal_conductivity=2.5,  
    density=2881,              
    heat_capacity=1600         
)

gb_material = F.Material(
    name="LiTiO3_amorphous", 
    D_0=9.3e-7,  
    E_D=0.41,    
    thermal_conductivity=2.5,  
    density=2881,
    heat_capacity=1600
)

# Volume subdomains for two regions
vol_bulk_grain = F.VolumeSubdomain(id=1, material=bulk_material, locator=bulk_grain_locator)
vol_grain_boundary = F.VolumeSubdomain(id=2, material=gb_material, locator=grain_boundary_locator)

# Surface subdomains (external boundaries only)
surf_left = F.SurfaceSubdomain(id=3, locator=left_boundary_locator)
surf_right = F.SurfaceSubdomain(id=4, locator=right_boundary_locator)
surf_top = F.SurfaceSubdomain(id=5, locator=top_boundary_locator)
surf_bottom = F.SurfaceSubdomain(id=6, locator=bottom_boundary_locator)

# Updated subdomains list
hydrogen_problem.subdomains = [
    vol_bulk_grain, 
    vol_grain_boundary, 
    surf_left, 
    surf_right, 
    surf_top, 
    surf_bottom
]

# Interface between bulk grain and grain boundary
interface_gb = F.Interface(
    id=7,
    subdomains=[vol_bulk_grain, vol_grain_boundary],
    penalty_term=1e25
)

hydrogen_problem.interfaces = [interface_gb]
hydrogen_problem.surface_to_volume = {
    surf_left: vol_bulk_grain,
    surf_right: vol_grain_boundary
}

# Species in both regions
mobile_T = F.Species("T", subdomains=[vol_bulk_grain, vol_grain_boundary])
hydrogen_problem.species = [mobile_T]

# Boundary conditions - source from left surface
boundary_left = F.FixedConcentrationBC(
    subdomain=surf_left, species=mobile_T, value=1e15
)

hydrogen_problem.boundary_conditions = [boundary_left]

final_time = 100.0 # seconds

# Initial conditions
temperature = 700  # K
hydrogen_problem.temperature = temperature

hydrogen_problem.settings = F.Settings(
    atol=1e10,
    rtol=1e-10,
    final_time=final_time,  
    transient=True,
    max_iterations=10000  
)
hydrogen_problem.settings.stepsize = F.Stepsize(1)  # 1 second time step

print("Creating tritium transport problem")

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

print("Individual problem exports configured")

# Run the grain problem
print("Initializing hydrogen transport simulation in 1D grains")
hydrogen_problem.initialise()
print("Running simulation...")
hydrogen_problem.run()
print("Simulation completed!")

def plot_one_grain_results():
    """
    Plot results for single grain + amorphous layer geometry
    """
    try:
        concentration_field = F.read_function_from_file(
            filename="results/2D/grain_level/tritium_transport/grain_tritium.bp",
            name="T",
            timestamp=final_time
        )
        
        V = concentration_field.function_space
        coords = V.tabulate_dof_coordinates()
        x_coords = coords[:, 0]
        y_coords = coords[:, 1]
        conc_values = concentration_field.x.array
        
        fig, ax = plt.subplots(1, 1, figsize=(10, 8))
        
        scatter = ax.scatter(x_coords*1e6, y_coords*1e6, c=conc_values,
                            cmap='viridis', s=15, alpha=0.8)
        ax.set_title('Tritium Concentration: Grain → Amorphous Layer', fontsize=14)
        ax.set_xlabel('X Position (μm)')
        ax.set_ylabel('Y Position (μm)')
        ax.axvline(x=grain_size*1e6, color='red', linestyle='--', alpha=0.7, 
                   linewidth=2, label='Grain-Amorphous Interface')
        ax.legend()
        ax.set_aspect('equal')
        ax.grid(True, alpha=0.3)
        plt.colorbar(scatter, ax=ax, label='Concentration (atoms/m³)')
        
        plt.tight_layout()
        plt.savefig('results/2D/grain_level/one_grain_results.png', dpi=300, bbox_inches='tight')
        plt.show()
        
        print(f"\n=== RESULTS ===")
        print(f"Mobile tritium concentration range: {np.min(conc_values):.2e} to {np.max(conc_values):.2e} atoms/m³")
        print(f"Amorphous Layer width: {gb_width*1e6:.1f} μm")
        print(f"Total width: {total_width*1e6:.1f} μm")
        
    except Exception as e:
        print(f"Error in plotting: {e}")

plot_one_grain_results()

fyi I tried it in 1D first and it works:

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

hydrogen_problem = F.HydrogenTransportProblem()

width_bulk = 10e-6  # 5 micrometers
width_amorphous = 1e-6  # 0.9 micrometers

vertices = np.concatenate(
    [
        np.linspace(0, width_bulk, num=600),
        np.linspace(width_bulk, (width_bulk + width_amorphous), num=600)
    ]
)

# Create the two-grain mesh
my_mesh = F.Mesh1D(vertices=vertices)
hydrogen_problem.mesh = my_mesh

bulk_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 
    density=2881,              # Density [kg/m³]
    heat_capacity=1600         # Heat capacity [J/kg/K]
)

gb_material = F.Material(
    name="LiTiO3_amorphous", 
    D_0=9.3e-7,  
    E_D=0.41,    
    thermal_conductivity=2.5,  
    density=2881,
    heat_capacity=1600
)

vol_left_grain = F.VolumeSubdomain1D(id=1, material=bulk_material, borders=[0, width_bulk])
vol_amorphous_layer = F.VolumeSubdomain1D(id=2, material=gb_material, borders=[width_bulk, (width_bulk + width_amorphous)])

left_surface = F.SurfaceSubdomain1D(id=3, x=0)
right_surface = F.SurfaceSubdomain1D(id=4, x=(width_bulk + width_amorphous))

hydrogen_problem.subdomains = [
    vol_left_grain,
    vol_amorphous_layer,
    left_surface,
    right_surface,
]

bulk_gb_interface = F.Interface(
    id=5, 
    subdomains=[vol_left_grain, vol_amorphous_layer],
    penalty_term=1e25  
    )

hydrogen_problem.interfaces = [bulk_gb_interface]
hydrogen_problem.surface_to_volume = {
    left_surface : vol_left_grain,
    right_surface : vol_amorphous_layer
}

mobile_T = F.Species("T", subdomains=[vol_left_grain, vol_amorphous_layer])
hydrogen_problem.species = [mobile_T]

hydrogen_problem.boundary_conditions = [
    F.FixedConcentrationBC(
        subdomain=left_surface, species=mobile_T, value=1e15
    ),
]


hydrogen_problem.temperature = 700  # K 

hydrogen_problem.settings = F.Settings(
    transient=True,
    atol=1e10,
    rtol=1e-10,
    final_time=200,
    stepsize=1,
    max_iterations=10000
)

print("Creating tritium transport problem")

# Exports for hydrogen problem
hydrogen_problem.exports = [
    F.VTXSpeciesExport(
        filename="results/1D/grain_level/tritium_transport/bulk_tritium.bp",
        field=[mobile_T],
    ),
    F.XDMFExport(
        filename="results/1D/grain_level/tritium_transport/grain_tritium.xdmf",
        field=[mobile_T],
    )
]

print("Individual problem exports configured")

# Run the grain problem
print("Initializing hydrogen transport simulation in 1D grains")
hydrogen_problem.initialise()
print("Running simulation...")
hydrogen_problem.run()
print("Simulation completed!")

I have also tried these other trouble shooting methods:

• a single rectangluar dolfinx mesh that is split into two regions.
• refining the mesh into the amorphous region
• using a tritium source term for the bulk instead of a fixed concentration bc on the left surface

Thank you in advance for your help!

[remidm]

Hi @griff,

This is a meshing/tagging issue.

If you plot the meshtags with:

import pyvista
from dolfinx import plot

mesh = hydrogen_problem.mesh.mesh
ft = hydrogen_problem.facet_meshtags
ct = hydrogen_problem.volume_meshtags

fdim = mesh.topology.dim - 1
tdim = mesh.topology.dim
mesh.topology.create_connectivity(fdim, tdim)

topology, cell_types, x = plot.vtk_mesh(mesh, tdim, ct.indices)

p = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
grid.cell_data["Cell Marker"] = ct.values
grid.set_active_scalars("Cell Marker")

p.add_mesh(grid, show_edges=False)
if pyvista.OFF_SCREEN:
    figure = p.screenshot("cell_marker.png")
p.show()

You will see that some cells are tagged with 0. This is because your locators aren’t defined properly.

How to fix this:

1. add some tolerances to your locators:

# Update region locators for the new geometry
def bulk_grain_locator(x):
    """Square grain region (left side)"""
    return x[0] < grain_size + 1e-12


def grain_boundary_locator(x):
    """Rectangle grain boundary region (right side)"""
    return x[0] >= grain_size - 1e-12

2. Since you use GMSH and generate tags, you can directly pass them to your problem and not define locators in VolumeSubdomain and SurfaceSubdomain:

# Volume subdomains for two regions
vol_bulk_grain = F.VolumeSubdomain(id=1, material=bulk_material)
vol_grain_boundary = F.VolumeSubdomain(id=2, material=gb_material)

# Surface subdomains (external boundaries only)
surf_left = F.SurfaceSubdomain(id=3)
surf_right = F.SurfaceSubdomain(id=4)
surf_top = F.SurfaceSubdomain(id=5)
surf_bottom = F.SurfaceSubdomain(id=6)

hydrogen_problem.facet_meshtags = facet_tags
hydrogen_problem.volume_meshtags = cell_tags

[remidm]

Note:

• If you want to impose a concentration discontinuity at the interface, then you should use F.HydrogenTransportProblemDiscontinuous and not F.HydrogenTransportProblem
• If you don’t want to use Paraview, I would recommend using pyvista to plot your fields in python.

import pyvista
from dolfinx import plot

topology, cell_types, x = plot.vtk_mesh(
    mobile_T.post_processing_solution.function_space
)

p = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
grid.point_data["C"] = mobile_T.post_processing_solution.x.array.real
grid.set_active_scalars("C")

p.add_mesh(grid, show_edges=False)
if pyvista.OFF_SCREEN:
    figure = p.screenshot("c.png")
p.show()

• In order to make your transient problem converge in more than zero iterations, I had to adjust the tolerances to atol=1e0, rtol=1e-20. This is a bug in dolfinx 0.9 that will be fixed in the next release. See Finite relative residuals for 0th Newton iteration in time-dependent problem - dolfinx - FEniCS Project, this PR fixed it, + the newton solver in dolfinx 0.10 will be completely replaced by SNES which should solve this issue altogether.

[griff]

Problem solved, thank you!