Dimensonless transport model showing 0 values in output data

[griffda]

Hi everyone,

As well as advice from the documentation written here I have seen some examples of transport simulations using scaling operations to help increase the number of solves.

Below is my attempt to do this for a tritium transport problem in a lithium ceramic. The issue is that the output data for every output variable returns a 0 value for every time-step. I have attempted to brute force non-zero values through changes to input parameters in order to figure out where the issue is without success.

My aim is to model transport in a 2D block of lithium ceramic material. I have generated the mesh externally using gmsh in a python script (shown at the bottom of this post), where bulk grains are modelled as squares 2 microns in width and amorphous grain boundaries between the bulk are 0.1microns in width. See screenshot of mesh below with gbs shown in red.

below is the input file. Any help would be hugely appreciated, thank you in advance.

# TMAP8 Tritium Transport in Li2TiO3 with Grain Boundaries
# Dimensionless formulation for numerical stability

# Physics:             tritium trapping and detrapping, diffusion
# Species variable:    mobile, trapped
# BCs:                 Vacuum boundary conditions at all boundaries

# ============================================================================
# CHARACTERISTIC SCALES
# ============================================================================

# Characteristic scales for non-dimensionalization
scale_length = 1e-6              # 1 μm (grain size scale)
scale_time = 1                   # 1 s
scale_quantity = 1e18            # atoms (choose based on expected inventory)

# Convert: dimensional → dimensionless
# [L] → L/scale_length
# [t] → t/scale_time  
# [C] → C/scale_quantity * scale_length^3

# ============================================================================
# DIMENSIONAL PARAMETERS (for reference)
# ============================================================================

# Physical parameters
temperature = 773                    # K
k_B = 8.617e-5                      # eV/K

# Bulk material properties (Li2TiO3 crystalline)
D_0_bulk_dim = 9.5e-6               # m^2/s
E_D_bulk = 0.47                     # eV
D_bulk_dim = ${fparse D_0_bulk_dim * exp(-E_D_bulk / (k_B * temperature))}
# = 9.5e-6 * exp(-0.47/0.0666) ≈ 1.35e-9 m²/s

# Grain boundary material properties  
D_0_gb_dim = 9.3e-7                 # m^2/s
E_D_gb = 0.41                       # eV
D_gb_dim = ${fparse D_0_gb_dim * exp(-E_D_gb / (k_B * temperature))}
# = 9.3e-7 * exp(-0.41/0.0666) ≈ 1.71e-10 m²/s

# Trapping parameters
k_t_prefactor = 1e13                # 1/s
trapping_energy_eV = 0.55           # eV
k_t_dim = ${fparse k_t_prefactor * exp(-trapping_energy_eV / (k_B * temperature))}
# = 1e13 * exp(-0.55/0.0666) ≈ 4.89e3 s⁻¹

k_r_prefactor = 1e13                # 1/s  
detrapping_energy_eV = 1.0          # eV (reduced from 1.35)
k_r_dim = ${fparse k_r_prefactor * exp(-detrapping_energy_eV / (k_B * temperature))}
# = 1e13 * exp(-1.0/0.0666) ≈ 2.8e-7 s⁻¹

# Trap site density 
#dpa_total = 0.2
#defects_per_dpa = 1e25              # defects/m³/dpa
#annealing_fraction = 0.5
#N_defects_dim = ${fparse dpa_total * defects_per_dpa * (1 - annealing_fraction)}
# = 0.2 * 1e25 * 0.5 = 1e24 m⁻³

# Trap site density
trapping_site_fraction = 1e-3       # Fraction of lattice sites that are traps
titanate_density = 5.9e28           # atoms/m³ (Li2TiO3 lattice density)
N_defects_dim = ${fparse trapping_site_fraction * titanate_density}
# = 1e-3 * 5.9e28 = 5.9e25 m⁻³

# Source terms
source_bulk_dim = 1e13              # atoms/m³/s
source_gb_dim = 1e13                # atoms/m³/s

# ===========================================================================
# DIMENSIONLESS PARAMETERS
# ============================================================================

# Diffusivities: D* = D × scale_time / scale_length²
D_bulk = ${fparse D_bulk_dim * scale_time / (scale_length * scale_length)}
# = 1.35e-9 * 1 / 1e-12 = 1350 (dimensionless)

D_gb = ${fparse D_gb_dim * scale_time / (scale_length * scale_length)}
# = 1.71e-10 * 1 / 1e-12 = 171 (dimensionless)

# Rate constants: k* = k × scale_time
k_t = ${fparse k_t_dim * scale_time}
# = 4.89e3 * 1 = 4890 (dimensionless)

k_r = ${fparse k_r_dim * scale_time}
# = 2.8e-7 * 1 = 2.8e-7 (dimensionless)

# Trap density: N* = N × scale_length³ / scale_quantity
N_traps = ${fparse N_defects_dim * scale_length*scale_length*scale_length / scale_quantity}
# = 1e24 * 1e-18 / 1e18 = 1.0 (dimensionless, O(1))

# Trapping rate: alpha_t = k_t × N_traps
alpha_t = ${fparse k_t * N_traps}
# = 4890 * 1.0 = 4890

# Detrapping rate: alpha_r = k_r
alpha_r = ${k_r}
# = 2.8e-7

# Source terms: S* = S × scale_time × scale_length³ / scale_quantity
source_bulk = ${fparse source_bulk_dim * scale_time * scale_length*scale_length*scale_length / scale_quantity}
# = 1e13 * 1 * 1e-18 / 1e18 = 1e-5 (dimensionless)

source_gb = ${fparse source_gb_dim * scale_time * scale_length*scale_length*scale_length / scale_quantity}
# = 1e13 * 1 * 1e-18 / 1e18 = 1e-5 (dimensionless)

# Simulation time (dimensionless)
end_time = 100                     # 100 s (dimensionless = dimensional here)

# ============================================================================
# MESH
# ============================================================================

[Mesh]
  [file]
    type = FileMeshGenerator
    file = 'grain_grid_10x10.msh'
  []
  [rename_blocks]
    type = RenameBlockGenerator
    input = file
    old_block = '1 2'  # Use Gmsh names
    new_block = 'bulk gb'
  []
  [rename_boundaries]
    type = RenameBoundaryGenerator
    input = rename_blocks
    old_boundary = '11 12 13 14'
    new_boundary = 'left right bottom top'
  []
[]

# ============================================================================
# VARIABLES (Dimensionless)
# ============================================================================

[Variables]
  [mobile]
    #initial_condition = 0        # Start from zero
    scaling = 1                   # Expect O(1) at steady state
  []
  [trapped]
    #initial_condition = 0
    scaling = 1                   # Expect O(1) when traps fill
    #block = 'bulk'
  []
[]

# ============================================================================
# AUXVARIABLES
# ============================================================================

[AuxVariables]
  [temp]
    initial_condition = ${temperature}
  []
  
  # For bounds (prevent negative concentrations)
  [bounds_dummy_mobile]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_trapped]
    order = FIRST
    family = LAGRANGE
    block = 'bulk'
  []
[]

# ============================================================================
# BOUNDS (Prevent negative concentrations)
# ============================================================================

[Bounds]
  [mobile_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_mobile
    bounded_variable = mobile
    bound_type = lower
    bound_value = 0.0
  []
  [trapped_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_trapped
    bounded_variable = trapped
    bound_type = lower
    bound_value = 0.0
  []
[]

# ============================================================================
# KERNELS
# ============================================================================

[Kernels]
  # Mobile species
  [mobile_time]
    type = TimeDerivative
    variable = mobile
  []
  [mobile_diffusion]
    type = MatDiffusion
    variable = mobile
    diffusivity = diffusivity
  []
  [mobile_source_bulk]
    type = BodyForce
    variable = mobile
    value = ${source_bulk}  # Same source everywhere
    block = 'bulk'
  []
  [mobile_source_gb]
    type = BodyForce
    variable = mobile
    value = ${source_gb}  # Same source everywhere
    block = 'gb'
  []

  # Coupling: mobile loses to trapping
  [mobile_trapped_coupling]
    type = CoupledTimeDerivative
    variable = mobile
    v = trapped
    block = 'bulk'
  []
  
  # Trapped species (only in bulk)
  [trapped_time]
    type = TimeDerivative
    variable = trapped
    block = 'bulk'
  []
  [trapping_reaction]
    type = ADMatCoupledDefectAnnihilation
    variable = trapped
    v = mobile
    eq_concentration = ${N_traps}
    annihilation_rate = ${alpha_t}
    block = 'bulk'
  []
  [detrapping_reaction]
    type = ADMatReaction
    variable = trapped
    reaction_rate = ${alpha_r}
    block = 'bulk'
  []
[]

# ============================================================================
# BOUNDARY CONDITIONS
# ============================================================================

[BCs]
  [left_bc]
    type = DirichletBC
    variable = mobile
    boundary = 'left'
    value = 0
  []
  [right_bc]
    type = DirichletBC
    variable = mobile
    boundary = 'right'
    value = 0
  []
  [top_bc]
    type = DirichletBC
    variable = mobile
    boundary = 'top'
    value = 0
  []
  [bottom_bc]
    type = DirichletBC
    variable = mobile
    boundary = 'bottom'
    value = 0
  []
[]

# ============================================================================
# MATERIALS
# ============================================================================

[Materials]
  [diffusivity_bulk]
    type = GenericConstantMaterial
    prop_names = 'diffusivity'
    prop_values = '${D_bulk}'
    block = 'bulk'
  []
  [diffusivity_gb]
    type = GenericConstantMaterial
    prop_names = 'diffusivity'
    prop_values = '${D_gb}'
    block = 'gb'
  []
[]

# ============================================================================
# POSTPROCESSORS
# ============================================================================

[Postprocessors]
  [total_mobile]
    type = ElementIntegralVariablePostprocessor
    variable = mobile
  []
  [total_trapped]
    type = ElementIntegralVariablePostprocessor
    variable = trapped
    block = 'bulk'
  []
  [avg_mobile]
    type = ElementAverageValue
    variable = mobile
  []
  [avg_trapped]
    type = ElementAverageValue
    variable = trapped
    block = 'bulk'
  []
  [max_trapped]
    type = ElementExtremeValue
    variable = trapped
    block = 'bulk'
  []
  [trap_saturation]
    type = ParsedPostprocessor
    pp_names = 'avg_trapped'
    expression = 'avg_trapped / ${N_traps}'
  []
  [flux_total]
    type = ParsedPostprocessor
    pp_names = 'flux_left flux_right flux_top flux_bottom'
    expression = 'abs(flux_left) + abs(flux_right) + abs(flux_top) + abs(flux_bottom)'
  []
  [flux_right]
    type = SideDiffusiveFluxIntegral
    variable = mobile
    boundary = 'right'
    diffusivity = diffusivity
  []
  [flux_left]
    type = SideDiffusiveFluxIntegral
    variable = mobile
    boundary = 'left'
    diffusivity = diffusivity
  []
  [flux_top]
    type = SideDiffusiveFluxIntegral
    variable = mobile
    boundary = 'top'
    diffusivity = diffusivity
  []
  [flux_bottom]
    type = SideDiffusiveFluxIntegral
    variable = mobile
    boundary = 'bottom'
    diffusivity = diffusivity
  []
[]

# ============================================================================
# PRECONDITIONING
# ============================================================================

[Preconditioning]
  [Newtonlu]
    type = SMP
    full = true
    solve_type = 'NEWTON'
    petsc_options_iname = '-pc_type -sub_pc_type -snes_type'
    petsc_options_value = 'asm      lu           vinewtonrsls' # This petsc option helps prevent negative concentrations'
  []
[]

# ============================================================================
# EXECUTIONER
# ============================================================================

[Executioner]
  type = Transient
  scheme = bdf2

  
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
  nl_max_its = 14
  
  l_tol = 1e-8
  l_max_its = 100
  
  start_time = 0
  end_time = ${end_time}
  
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-4                       # Start small but not tiny
    dtmax = 10                      # Allow reasonable growth
    optimal_iterations = 10
    iteration_window = 2
    growth_factor = 1.2
    cutback_factor = 0.75
    cutback_factor_at_failure = 0.5
  []
[]


# ============================================================================
# OUTPUTS
# ============================================================================

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  
  [console]
    type = Console
    fit_mode = 80
  []
[]

from mpi4py import MPI  
import os
import gmsh
from dolfinx.io import gmsh as gmsh_io 


# Geometry parameters
grain_size = 2e-6       # 2 micrometers
gb_width = 0.11e-6       # 0.1 micrometers grain boundary width

grid_size = 10            # 10×10 grid for ~20μm total width

def create_large_grain_mesh(save_mesh=True, mesh_output_name=None):
    """
    Create NxN grain grid mesh
    """
    
    gmsh.initialize()
    gmsh.model.add(f"grain_grid_{grid_size}x{grid_size}")
    
    # Calculate dimensions for NxN grid (same formula as 3x3)
    total_width = grid_size * grain_size + (grid_size - 1) * gb_width
    total_height = grid_size * grain_size + (grid_size - 1) * gb_width
    
    print(f"Creating {grid_size}×{grid_size} grain grid:")
    print(f"  - Grain size: {grain_size*1e6:.1f} μm")
    print(f"  - GB width: {gb_width*1e6:.3f} μm")
    print(f"  - Total dimensions: {total_width*1e6:.1f} × {total_height*1e6:.1f} μm")

    # Use built-in factory
    factory = gmsh.model.occ

    # Create grains individually
    grain_tags = []
    grain_id = 1
    
    for row in range(grid_size):
        for col in range(grid_size):
            x_pos = col * (grain_size + gb_width)
            y_pos = row * (grain_size + gb_width)
            grain_tag = factory.addRectangle(x_pos, y_pos, 0, grain_size, grain_size)
            grain_tags.append(grain_tag)
            grain_id += 1

    # Create vertical grain boundaries
    vgb_tags = []
    for col in range(grid_size - 1):  # (grid_size - 1) vertical GB strips
        x_pos = (col + 1) * grain_size + col * gb_width
        for row in range(grid_size):
            y_pos = row * (grain_size + gb_width)
            vgb_tag = factory.addRectangle(x_pos, y_pos, 0, gb_width, grain_size)
            vgb_tags.append(vgb_tag)

    # Create horizontal grain boundaries
    hgb_tags = []
    for row in range(grid_size - 1):  # (grid_size - 1) horizontal GB strips
        y_pos = (row + 1) * grain_size + row * gb_width
        for col in range(grid_size):
            x_pos = col * (grain_size + gb_width)
            hgb_tag = factory.addRectangle(x_pos, y_pos, 0, grain_size, gb_width)
            hgb_tags.append(hgb_tag)

    # Create junction grain boundaries 
    junction_tags = []
    for row in range(grid_size - 1):
        for col in range(grid_size - 1):
            x_pos = (col + 1) * grain_size + col * gb_width
            y_pos = (row + 1) * grain_size + row * gb_width
            junction_tag = factory.addRectangle(x_pos, y_pos, 0, gb_width, gb_width)
            junction_tags.append(junction_tag)
    
    # Synchronize before adding physical groups
    factory.synchronize()
    
    # Set mesh size
    mesh_size = grain_size / 3 
    print(f"Using mesh size: {mesh_size*1e6:.3f} μm")
    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)
    
    # Physical groups all grains together
    gmsh.model.addPhysicalGroup(2, grain_tags, 1, "bulk_grains")
    
    # All grain boundaries as one group 
    all_gb = vgb_tags + hgb_tags + junction_tags
    gmsh.model.addPhysicalGroup(2, all_gb, 2, "grain_boundaries")
    
    print(f"Created geometry:")
    print(f"  - {len(grain_tags)} individual grains")
    print(f"  - {len(vgb_tags)} vertical grain boundaries") 
    print(f"  - {len(hgb_tags)} horizontal grain boundaries")
    print(f"  - {len(junction_tags)} junction grain boundaries")
    print(f"  - {len(all_gb)} total grain boundary elements")
    
    # Get boundary curves (same logic as 3x3)
    all_curves = gmsh.model.getEntities(1)
    
    left_curves = []
    right_curves = []
    bottom_curves = []
    top_curves = []
    
    # Boundary detection
    for curve_dim, curve_tag in all_curves:
        try:
            bbox = gmsh.model.getBoundingBox(curve_dim, curve_tag)
            x_min, y_min, z_min, x_max, y_max, z_max = bbox
            
            tol = grain_size / 5  
            
            if x_min < tol and x_max < tol:
                left_curves.append(curve_tag)
            elif (total_width - x_min) < tol and (total_width - x_max) < tol:
                right_curves.append(curve_tag)
            elif y_min < tol and y_max < tol:
                bottom_curves.append(curve_tag)
            elif (total_height - y_min) < tol and (total_height - y_max) < tol:
                top_curves.append(curve_tag)
        except:
            continue
    
    # Add boundary physical groups 
    print(f"Boundary detection results:")
    print(f"  - Total curves: {len(all_curves)}")
    
    if left_curves:
        gmsh.model.addPhysicalGroup(1, left_curves, 11, "left_boundary")
        print(f"  - Left boundary: {len(left_curves)} curves")
    else:
        print("  - WARNING: No left boundary curves found!")
        
    if right_curves:
        gmsh.model.addPhysicalGroup(1, right_curves, 12, "right_boundary")
        print(f"  - Right boundary: {len(right_curves)} curves")
    else:
        print("  - WARNING: No right boundary curves found!")
        
    if bottom_curves:
        gmsh.model.addPhysicalGroup(1, bottom_curves, 13, "bottom_boundary")
        print(f"  - Bottom boundary: {len(bottom_curves)} curves")
    else:
        print("  - WARNING: No bottom boundary curves found!")
        
    if top_curves:
        gmsh.model.addPhysicalGroup(1, top_curves, 14, "top_boundary")
        print(f"  - Top boundary: {len(top_curves)} curves")
    else:
        print("  - WARNING: No top boundary curves found!")
    
    # Mesh generation
    gmsh.option.setNumber("Mesh.Algorithm", 6)  
    print("Generating mesh (this may take longer for large grids)...")
    
    try:
        gmsh.model.mesh.generate(2)
        gmsh.model.mesh.removeDuplicateNodes()
        
        # Verify mesh creation
        node_tags, coord, parametric_coord = gmsh.model.mesh.getNodes()
        print(f"✓ Mesh generated successfully with {len(node_tags)} nodes")
        
    except Exception as e:
        print(f"Mesh generation failed: {e}")
        # Try with slightly coarser mesh if it fails
        print("Trying with slightly coarser mesh...")
        coarser_mesh_size = grain_size / 2  # Less fine than /30
        gmsh.model.mesh.setSize(gmsh.model.getEntities(0), coarser_mesh_size)
        gmsh.model.mesh.generate(2)
        print(f"Generated with coarser mesh size: {coarser_mesh_size*1e6:.3f} μm")
    
    # Save mesh file
    mesh_filename = f"temp_{grid_size}x{grid_size}_grid.msh"
    gmsh.write(mesh_filename)
    
    # Save permanent copy BEFORE finalizing
    if save_mesh:
        if mesh_output_name is None:
            mesh_output_name = f"grain_grid_{grid_size}x{grid_size}.msh"
        
        import shutil
        shutil.copy(mesh_filename, mesh_output_name)
        print(f"✓ Mesh saved to: {mesh_output_name}")

            # Save in VTK format for ParaView
    vtk_filename = f"grain_grid_{grid_size}x{grid_size}.vtk"
    gmsh.write(vtk_filename)
    print(f"✓ VTK file saved to: {vtk_filename}")
    
    gmsh.finalize()
    
    # Load mesh with DOLFINx - returns 2 values
    mesh_data = gmsh_io.read_from_msh(
        mesh_filename, MPI.COMM_WORLD, 0, gdim=2
    )

    # Unpack the returned tuple
    mesh_2d = mesh_data[0]
    cell_tags = mesh_data[1]
    facet_tags = mesh_data[2]
    
    os.remove(mesh_filename)  # Remove temp file

    print(f"✓ {grid_size}×{grid_size} grain grid mesh created!")
    print(f"  - Total grains: {grid_size**2}")
    print(f"  - Actual dimensions: {total_width*1e6:.1f} × {total_height*1e6:.1f} μm")
    print(f"  - Mesh cells: {mesh_2d.topology.index_map(2).size_local}")
    print(f"  - Mesh vertices: {mesh_2d.topology.index_map(0).size_local}")
    
    return mesh_2d, total_width, total_height, cell_tags, facet_tags

# Create and save the mesh
mesh_2d, total_width, total_height, cell_tags, facet_tags = create_large_grain_mesh(
    save_mesh=True, 
    mesh_output_name=f"grain_grid_{grid_size}x{grid_size}.msh"
)

[simopier]

Hi @griffda, sorry for the delay on this discussion. I have a couple comments:

• We advise using nodal kernels for trapped species in general since trapped species do not diffuse. See item 5 in General TMAP Questions #7. Ver-1d is a good example of how nodal kernels can be used. That being said, it shouldn't be the cause of your issue here.
• Does the same input file work on a mesh with only bulk?

[simopier]

Playing with the input file a bit more, I have found that:

• The final concentration is equal to the value the user sets on the Dirichlet BCs, which makes me think that the Diffusivity might be too high (and the generation term too low), and the simulation immediately finds the steady state where the concentration field is flat and, in practice, imposed by the value set at the boundaries. By reducing the diffusivity by several orders of magnitude, I managed to get a non-uniform concentration field.
• The model currently imposes a generation term everywhere, including at the boundaries, and also imposes a set value at the BCs using a DirichletBC. This can be contradictory (how can a value increase and remain 0) and can lead to convergence challenges in some cases. This might not be the main issue here, but I suggest looking into modeling the reactions at the boundaries rather than imposing a value that can conflict with the generation term. This validation case can serve as an example: https://mooseframework.inl.gov/TMAP8/verification_and_validation/val-2f.html.
• you should move dtmax from the TimeSteppers block and put it in the overall Executioner block around it. It is currently unused in TimeSteppers

I hope that helps.

[griffda]

Hi @simopier, thank you for coming back to me.

• I have used the trapping nodal kernel in some other test simulations but attempted to apply this trapping logic to my own problem. I will retry with the example you have linked.
• In regard to the use of Dirichlet BCs, I agree that these are making the problem stiff. I am trying to resolve this by implementing surface reaction bcs instead. This is obviously also a more physical representation of the problem. I have seen help pages for this here, with ratedep conditions being most applicable to my problem: i.e., tritium leaving (recombination from the surface) of the ceramic by: 2T(surface) → T₂(gas) and tritium gas dissociating and re-entering the ceramic: T₂(gas) → 2T(surface). The forward reaction rate is: K_r = k_r0 * exp(-E_kr / (k_B * T)) and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T)). According to the docs this is a good example of its implementation. However, I see you have instead suggsted another validation case to capture this. Could you please explain why your suggestion should be used instead of the one I have linked?

[simopier]

Hi @griffda,

Did you manage to figure out if it was indeed the diffusivity that was too high in the model? That seemed to be the case from my limited testing with the input file, as described above. But let me know if that's not the case and we can look more into it.

In regard to the use of Dirichlet BCs, I agree that these are making the problem stiff. I am trying to resolve this by implementing surface reaction bcs instead. This is obviously also a more physical representation of the problem. I have seen help pages for this here, with ratedep conditions being most applicable to my problem: i.e., tritium leaving (recombination from the surface) of the ceramic by: 2T(surface) → T₂(gas) and tritium gas dissociating and re-entering the ceramic: T₂(gas) → 2T(surface). The forward reaction rate is: K_r = k_r0 * exp(-E_kr / (k_B * T)) and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T)). According to the docs this is a good example of its implementation. However, I see you have instead suggsted another validation case to capture this. Could you please explain why your suggestion should be used instead of the one I have linked?

Great question! The verification case that was listed in the theory manual is a good example of the implementation, but the reactions happen throughout the whole 2D domain. Physically, the simulation domain represents the surface.
However, the validation case I listed above (val-2f) relates more to the case you are modeling here. The surface equations are the same (althought the dissociation rate is assumed to be negligeable), but val-2f models a sample with the surface reactions happening at the boundaries (as boundary conditions, as opposed to everywhere in the domain like in ver-1ia).
If you want to see another example with similar surface reaction, but with the dissociation rate included too, you can check out val-2e, and val-2ee in particular.

[griffda]

Hi @simopier

Thank you for your repsonse.

I did play around with the diffusivity, however reducing it by several orders of magnitude as you suggest means that we are way off the value we calculate in our atomistic simulations for T in the ceramic.

I will work on implementing the surface reaction in the val-2f case.

Regarding the use of scaling operations, could you advise on whether you deem this necessary here?

Thanks!

[simopier]

Hi @griffda,

I wanted to reach out again to see if you were still having some issues about this simulation.

The diffusivity indeed has to be provided a reasonable value. I was changing it a lot solely for the purpose of debugging, not to imply that you should use a different diffusivity. My point was that with Dirichlet BCs, a small domain size, and such diffusivity values, then I was not surprised that you were getting a constant field equal to the value provided as a BC.

So the main solutions are to either increase the domain of the simulation (which can lead to increase in computational costs), or use a surface reaction BC, which you were hopefully able to deploy.

I am sorry I did not respond to your question about the scaling operation earlier, I somehow missed this question.
We indeed recommend the use of scaling operation if your concentration of trapped tritium tends to be orders of magnitude different than the concentration of mobile tritium. However, at the debugging stage, I would recommend starting without and only use it if you get numerical convergence issues, or to optimize convergence later on.

I hope that helps.