Pore-Scale Simulations of Tritium Transport Using TMAP8

This example demonstrates TMAP8's capability to (1) generate pore structures from input images, and (2) perform pore-scale simulations of tritium transport on these pore structures based on the model described in Simon et al. (2022). Simon et al. (2022) provides additional context and details about this example, including a description of the calibration efforts. Note that this approach can be extended to three-dimensional microstructures.

This example, as stated in Simon et al. (2022), showcases how TMAP8 can simulate complex geometries and help investigate the effects of pore microstructure on tritium absorption.

General Description of the Simulation Case and Corresponding Input Files

Introduction

The sections below describe the simulation details and explain how they translate into the example TMAP8 input files. Not every part of the input file is explained here. If the interested reader has more questions about this example case and how to modify this input file to adapt it to a different case, feel free to reach out to the TMAP8 development team on the TMAP8 GitHub discussion page.

In this example, TMAP8 is used to simulate tritium transport at the pore scale in ceramic breeder materials. The model tracks tritium in three different forms: freely moving in solid solution in the ceramic material ( $var element = document.getElementById("moose-equation-2f1a0f5e-792a-4f91-aa9d-b66c35122679");katex.render("T_s", element, {displayMode:false,throwOnError:false});$ ), trapped in the ceramic material ( $var element = document.getElementById("moose-equation-44ad1d2f-828a-4efe-b21e-c847f3d76851");katex.render("T_t", element, {displayMode:false,throwOnError:false});$ ), or in gaseous form in pores and outside the ceramic ( $var element = document.getElementById("moose-equation-8ef96403-2559-456e-80f3-3772d3ef5ad0");katex.render("T_2", element, {displayMode:false,throwOnError:false});$ ). This example is based on the simulations presented in Section V of Simon et al. (2022).

This example is divided into the two steps required to reproduce the results presented in Section V of Simon et al. (2022). First, it describes the simulation used to input an image and automatically generate a usable pore microstructure. Second, it shows how this microstructure can be imported to perform tritium transport calculations on pore microstructure. This process is illustrated in Figure 1.

Figure 1: Illustration of the process for pore scale simulation of tritium transport in TMAP8 based on a microstructure image. Step 1 takes an image and creates a corresponding microstructure, and step 2 performs tritium transport calculations on this microstructure (figures are reproduced from Simon et al. (2022)).

note:TMAP8 can also automatically generate microstructures, or import other simulation results

TMAP8 can also automatically generate microstructures using the initial condition system, as MOOSE is equipped with the user-friendly, built-in MeshGenerators system for creating meshes based on simple geometries. It can also import EBSD data or import microstructures generated via other means.

Generating the pore microstructure based on an image

In this first step, we use the following input file to import the image and create the smooth microstructure:

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Info:
# - This simulation predicts phase evolution based on a grain growth phase field model.
# - The pore and the ceramics represent two different phases, and the initial conditions
#   are provided with a black and white picture, which is read by ImageFunction.
# - The simulation smoothens the interface between the two phases.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_microstructure_reader_smoothing_base.i pore_structure_open.params

# mesh information
num_nodes_x = ${fparse 2*120} # (-)
num_nodes_y = ${fparse 2*120} # (-)
domain_start_x = ${units 0 mum}
domain_start_y = ${units 0 mum}
domain_end_x = ${units 5425 mum}
domain_end_y = ${units 5425 mum}

# grain growth parameters (since the point of the simulation is only to get a smooth interface, the model parameters can be selected arbitrarily)
op_num = 2 # Number of grains (-)
mobility_prefactor = ${units 2.5e-6 mum^4/J/s}
GB_energy = ${units 0.7 J/mum^2}
activation_energy = ${units 0.23 eV}
temperature = ${units 700 K}
width_diffuse_GB = ${units 50 mum} # Width of the diffuse GB

# image function option
threshold_image_function = 255.5 # (-)

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  [gen]
    type = DistributedRectilinearMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/DistributedRectilinearMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_nodes_x}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_nodes_y}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${domain_start_x}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${domain_end_x}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = ${domain_start_y}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${domain_end_y}
  []
[]

[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
  op_num = ${op_num}
  var_name_base = gr # Base name of grains
[]

[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
  [grain_tracker]
    type = GrainTracker<<<{"description": "Grain Tracker object for running reduced order parameter simulations without grain coalescence.", "href": "../../source/postprocessors/GrainTracker.html"}>>>
    flood_entity_type<<<{"description": "Determines whether the flood algorithm runs on nodes or elements"}>>> = ELEMENTAL
    compute_halo_maps<<<{"description": "Instruct the Postprocessor to communicate proper halo information to all ranks"}>>> = true
  []
[]

[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
  [image_func0]
    type = ImageFunction<<<{"description": "Function with values sampled from an image or image stack.", "href": "../../source/functions/ImageFunction.html"}>>>
    file<<<{"description": "Name of single image file to extract mesh parameters from.  If provided, a 2D mesh is created."}>>> = ${input_name}
    threshold<<<{"description": "The threshold value"}>>> = ${threshold_image_function}
    upper_value<<<{"description": "The value to set for data greater than the threshold value"}>>> = 1
    lower_value<<<{"description": "The value to set for data less than the threshold value"}>>> = 0
  []
  [image_func1]
    type = ImageFunction<<<{"description": "Function with values sampled from an image or image stack.", "href": "../../source/functions/ImageFunction.html"}>>>
    file<<<{"description": "Name of single image file to extract mesh parameters from.  If provided, a 2D mesh is created."}>>> = ${input_name}
    threshold<<<{"description": "The threshold value"}>>> = ${threshold_image_function}
    upper_value<<<{"description": "The value to set for data greater than the threshold value"}>>> = 0
    lower_value<<<{"description": "The value to set for data less than the threshold value"}>>> = 1
  []
[]

[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [gr0_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    function<<<{"description": "The initial condition function."}>>> = image_func0
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = gr0
  []
  [gr1_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    function<<<{"description": "The initial condition function."}>>> = image_func1
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = gr1
  []
[]

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [PolycrystalVariables<<<{"href": "../../syntax/Variables/PolycrystalVariables/index.html"}>>>]
    # Custom action that creates all of the phase variables and sets their initial condition
  []
[]

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [bnds]
    # Variable used to visualize the interfaces in the simulation
  []
  [unique_grains]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
[]

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [PolycrystalKernel<<<{"href": "../../syntax/Kernels/PolycrystalKernel/index.html"}>>>]
    # Custom action creating all necessary kernels for phase evolution.  All input parameters are up in GlobalParams
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [bnds_aux]
    type = BndsCalcAux<<<{"description": "Calculate location of grain boundaries in a polycrystalline sample", "href": "../../source/auxkernels/BndsCalcAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = bnds
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'initial timestep_end'
  []
  [unique_grains]
    type = FeatureFloodCountAux<<<{"description": "Feature detection by connectivity analysis", "href": "../../source/auxkernels/FeatureFloodCountAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = unique_grains
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    flood_counter<<<{"description": "The FeatureFloodCount UserObject to get values from."}>>> = grain_tracker
    field_display<<<{"description": "Determines how the auxilary field should be colored. (UNIQUE_REGION and VARIABLE_COLORING are nodal, CENTROID is elemental, default: UNIQUE_REGION)"}>>> = UNIQUE_REGION
  []
[]

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [InterfaceEvolution]
    type = GBEvolution<<<{"description": "Computes necessary material properties for the isotropic grain growth model", "href": "../../source/materials/GBEvolution.html"}>>>
    GBmob0<<<{"description": "Grain boundary mobility prefactor in m^4/(J*s)"}>>> = ${mobility_prefactor}
    GBenergy<<<{"description": "Grain boundary energy in J/m^2"}>>> = ${GB_energy}
    Q<<<{"description": "Grain boundary migration activation energy in eV"}>>> = ${activation_energy}
    T<<<{"description": "Temperature in Kelvin"}>>> = ${temperature}
    wGB<<<{"description": "Diffuse GB width in the length scale of the model"}>>> = ${width_diffuse_GB}
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  scheme = bdf2
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  l_max_its = 20
  nl_max_its = 9
  nl_rel_tol = 1e-10
  [TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    cutback_factor = 0.9
    dt = 3e-2
    growth_factor = 1.1
    optimal_iterations = 7
  []
  dtmax = 25

  start_time = 0.0
  end_time = 100
[]

[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = ${output_name}
  [exodus]
    type = Exodus<<<{"description": "Object for output data in the Exodus format", "href": "../../source/outputs/Exodus.html"}>>>
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'FINAL'
  []
[]

This section describes this input file, what it does, and how to run it to replicate the results from Simon et al. (2022). Note that the input file also contains extensive comments meant to help users.

Import the pore microstructure image

Figure 2 and Figure 3 show the initial images used to generate the microstructures from Fig. 5 in Simon et al. (2022). Although both samples have a density of 80%, they exhibit very different pore interconnectivities, which was demonstrated to affect tritium absorption (Simon et al., 2022).

Figure 2: Initial image used to generate the microstructure with closed pores. This corresponds to the original image to generate Fig. 5.a in Simon et al. (2022).

Figure 3: Initial image used to generate the microstructure with open pores. This corresponds to the original image to generate Fig. 5.b in Simon et al. (2022).

TMAP8 can import these images using the ImageFunction object and, using a threshold, attribute different contrast values to given phases. Here, the bright areas are identified as ceramic material, while dark areas are identified as pores.

The size of the microstructure is not directly provided by the images in Figure 2 and Figure 3. Users can scale, crop, and adjust the image to the desired size. In this example, the sample being 4.5 mm in radius, the domain size is set to 5.425 mm $var element = document.getElementById("moose-equation-50726c59-8822-47b4-a106-73bf11c1eadb");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 5.425 mm.

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  [gen]
    type = DistributedRectilinearMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/DistributedRectilinearMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_nodes_x}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_nodes_y}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${domain_start_x}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${domain_end_x}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = ${domain_start_y}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${domain_end_y}
  []
[]

[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
  [image_func0]
    type = ImageFunction<<<{"description": "Function with values sampled from an image or image stack.", "href": "../../source/functions/ImageFunction.html"}>>>
    file<<<{"description": "Name of single image file to extract mesh parameters from.  If provided, a 2D mesh is created."}>>> = ${input_name}
    threshold<<<{"description": "The threshold value"}>>> = ${threshold_image_function}
    upper_value<<<{"description": "The value to set for data greater than the threshold value"}>>> = 1
    lower_value<<<{"description": "The value to set for data less than the threshold value"}>>> = 0
  []
  [image_func1]
    type = ImageFunction<<<{"description": "Function with values sampled from an image or image stack.", "href": "../../source/functions/ImageFunction.html"}>>>
    file<<<{"description": "Name of single image file to extract mesh parameters from.  If provided, a 2D mesh is created."}>>> = ${input_name}
    threshold<<<{"description": "The threshold value"}>>> = ${threshold_image_function}
    upper_value<<<{"description": "The value to set for data greater than the threshold value"}>>> = 0
    lower_value<<<{"description": "The value to set for data less than the threshold value"}>>> = 1
  []
[]

Perform a phase field simulation to smooth the interfaces

The rest of the input file is set up as a phase field simulation to smoothen the interface between the pore and the ceramic material. As discussed in Simon et al. (2022), using a phase field approach with a continuous interface enables TMAP8 to model arbitrarily complex geometries that cannot be easily meshed. Since TMAP8 inherits the phase field module from MOOSE, it can perform phase field simulations.

In this case, we use a two-phase model with arbitrary parameter values to quickly obtain, from the pore structures illustrated in Figure 2 and Figure 3, corresponding microstructures with the desired interface width, i.e., $var element = document.getElementById("moose-equation-eb52ef2a-784e-43cb-a081-bb1a8fed3b10");katex.render("l", element, {displayMode:false,throwOnError:false});$ =50 $var element = document.getElementById("moose-equation-5f93b5d5-fb09-4a08-aa3d-874c1208ac40");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m in this case. The simulation time is selected to be large enough to obtain a smooth interface, while staying small enough to prevent significant changes in the microstructure.

Run simulations

To run the input file to reproduce the microstructures shown in Fig. 5 in Simon et al. (2022), users need to first install TMAP8. Then, to generate the microstructure with closed pores using 4 CPU cores in parallel with mpirun, run:


cd ~/projects/TMAP8/test/tests/pore_scale_transport/
mpirun -np 4 ~/projects/TMAP8/tmap8-opt -i 2D_microstructure_reader_smoothing_base.i pore_structure_closed.params

and to generate the microstructure with open pores using 4 CPU cores in parallel with mpirun, run:


cd ~/projects/TMAP8/test/tests/pore_scale_transport/
mpirun -np 4 ~/projects/TMAP8/tmap8-opt -i 2D_microstructure_reader_smoothing_base.i pore_structure_open.params

Note that the pore_structure_closed.params and pore_structure_open.params files complement the base input file 2D_microstructure_reader_smoothing_base.i by providing input and output names.

Resulting smooth pore microstructures

Figure 4 shows the smooth pore microstructures obtained with the simulation described above. These pore structures correspond to the ones shown in Fig. 5 in Simon et al. (2022) and can be found in ~/projects/TMAP8/test/tests/pore_scale_transport/ as pore_structure_close/pore_structure_close.e and pore_structure_open/pore_structure_open.e once the simulations are completed. These can be viewed in a visualization software such as ParaView.

Figure 4: Smooth pore microstructures obtained with the simulation described above. This corresponds Fig. 5 in Simon et al. (2022).

Model tritium transport on pore-scale microstructures

In this second step, we use the following input file to import the smooth microstructures from Figure 4 and model pore scale tritium transport during absorption:

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Phase Field Model:   Isotropic bulk diffusion, pore surface reactions, isotropic pore diffusion
# Type:                Transient
# Phase structure:     Contains two phases: the bulk of a material, and the pores. An interface separates the two
# Physics:             tritium trapping and detrapping, surface reactions, diffusion
# Species variable:    tritium_s (Ts), tritium_trap (trapped T), tritium_2g (T2(g))
# Species AuxVariable: --
# Chemistry:           at the pore surface: (s) means in the solid, (g) means in gas form
#    Diatomic molecule dissolving at surface
#                                 T2(g) -> 2 T(s) rate K [T2] (1-theta)^2
#    Diatomic molecule combining at surface
#                                 2 T(s) -> T2(g) rate K [T]^2
#
# BCs:                 No flux for tritium and tritium_trap, zero at pore center for gaseous tritium in every form
#
#
# Info:
# - This input file uses the microstructure generated by 2D_microstructure_reader_smoothing_base.i
#   and performs tritium transport simulation on this microstructure.
# - Basic 2D non-dimensional input file for the transport of a solute in a microstructure with pores.
# - The goal of this input base and associated parameter files is to have the same simulation with different pore configurations
#   to show the effect of microstructure on tritium release.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_absorption_base.i pore_structure_open_absorption.params

# Scaling factors
scale_quantity = ${units 1e-18 moles}
scale_time = ${units 1 s}
scale_length = ${units 1 mum}
tritium_trap_scaling = 1e3 # (-)
tritium_s_scaling = 1e2 # (-)

# Conditions
tritium_2g_concentration = ${units 0.1121 mol/mum^3}

# Material properties
tritium_trap_0 = ${units 3.472805925e-18 mol/mum^3}
tritium_trap_0_scaled = ${fparse tritium_trap_0*scale_length^3/scale_quantity} # non-dimensional
available_sites_mult = 4.313540691 # (-) should be >1, otherwise there are fewer 'available sites' than 'trapping sites'
Diffusion_min = ${units 0 mum^2/s} # used as diffusion coefficients where I don't want the species to move (more stable than 0)
Diffusion_coefficient_Db_D0 = ${units 1649969.728 mum^2/s}
Diffusion_coefficient_D_pore_D0 = ${units 224767738.6 mum^2/s}
detrapping_rate_0 = ${units -4.065104516 1/s} # should be negative
trapping_rate_0 = ${units 8.333146645e+16 1/s/mol} # should be positive
reactionRateSurface_sXYg_ref_0 = ${units 1549103878000000.0 1/s/mol}
reactionRateSurface_gXYs_ref_0 = ${units 34.69205021 1/s}

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  file = ${input_name}
[]

[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
  [initial_microstructure]
    type = SolutionUserObject<<<{"description": "Reads a variable from a mesh in one simulation to another", "href": "../../source/userobjects/SolutionUserObject.html"}>>>
    mesh<<<{"description": "The name of the mesh file (must be xda/xdr or exodusII file)."}>>> = ${input_name}
    timestep<<<{"description": "Index of the single timestep used or \"LATEST\" for the last timestep (exodusII only).  If not supplied, time interpolation will occur."}>>> = LATEST
  []
[]

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [tritium_s]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = ${tritium_s_scaling}
  []
  [tritium_trap]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = ${tritium_trap_scaling}
  []
  [tritium_2g]
  []
[]

[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [tritium_2g_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = tritium_2g
    function<<<{"description": "The initial condition function."}>>> = 'gaseous_concentration_1_function_ic'
  []
[]

[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
  [pore_position_func]
    type = SolutionFunction<<<{"description": "Function for reading a solution from file.", "href": "../../source/functions/SolutionFunction.html"}>>>
    solution<<<{"description": "The SolutionUserObject to extract data from."}>>> = initial_microstructure
    from_variable<<<{"description": "The name of the variable in the file that is to be extracted"}>>> = gr0
  []
  [gaseous_concentration_1_function_ic]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'len_pore  position_pore_int'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = '25   4521'
    expression<<<{"description": "The user defined function."}>>> = '${tritium_2g_concentration}*0.5*(1.0+tanh((sqrt(x*x+y*y)-position_pore_int)/(len_pore/2.2)))'
  []
[]

[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [tritium_2g_sides_d] # Fix concentration on the sides
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'right top'
    value<<<{"description": "Value of the BC"}>>> = ${tritium_2g_concentration}
  []
[]

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [pore]
    initial_from_file_var<<<{"description": "Gives the name of a variable for which to read an initial condition from a mesh file"}>>> = gr0
    initial_from_file_timestep<<<{"description": "Gives the timestep (or \"LATEST\") for which to read a solution from a file for a given variable. (Default: LATEST)"}>>> = LATEST
  []
  # Used to prevent negative concentrations
  [bounds_dummy_tritium_s]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = LAGRANGE
  []
  [bounds_dummy_tritium_2g]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = LAGRANGE
  []
[]

[Bounds<<<{"href": "../../syntax/Bounds/index.html"}>>>]
  # To prevent negative concentrations
  [tritium_s_lower_bound]
    type = ConstantBounds<<<{"description": "Provides constant bound of a variable for the PETSc's variational inequalities solver", "href": "../../source/bounds/ConstantBounds.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = bounds_dummy_tritium_s
    bounded_variable<<<{"description": "The variable to be bounded"}>>> = tritium_s
    bound_type<<<{"description": "Type of bound. 'upper' refers to the upper bound. 'lower' refers to the lower value."}>>> = lower
    bound_value<<<{"description": "The value of bound for the variable"}>>> = 0.
  []
  [tritium_g_lower_bound]
    type = ConstantBounds<<<{"description": "Provides constant bound of a variable for the PETSc's variational inequalities solver", "href": "../../source/bounds/ConstantBounds.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = bounds_dummy_tritium_2g
    bounded_variable<<<{"description": "The variable to be bounded"}>>> = tritium_2g
    bound_type<<<{"description": "Type of bound. 'upper' refers to the upper bound. 'lower' refers to the lower value."}>>> = lower
    bound_value<<<{"description": "The value of bound for the variable"}>>> = 0.
  []
[]

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  # ========================================== Time derivative for each variable
  [dtritium_sdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_s'
  []
  [dtritium_2gdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_2g'
  []
  [dtritium_trapdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
  []
  # ============================================================= Bulk diffusion
  [Diff_tritium]
    type = ADMatDiffusion<<<{"description": "Diffusion equation kernel that takes an isotropic diffusivity from a material property", "href": "../../source/kernels/ADMatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = D_tritium
  []
  # ============================================================= Pore diffusion
  [Diff_tritium2g]
    type = ADMatDiffusion<<<{"description": "Diffusion equation kernel that takes an isotropic diffusivity from a material property", "href": "../../source/kernels/ADMatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = D_pore
  []
  # ================================================ Tritium trapping/detrapping
  [Trapping]
    type = CoupledTimeDerivative<<<{"description": "Time derivative Kernel that acts on a coupled variable. Weak form: $(\\psi_i, \\frac{\\partial v_h}{\\partial t})$.", "href": "../../source/kernels/CoupledTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    v<<<{"description": "Coupled variable"}>>> = 'tritium_trap'
  []
  [Detrapping_trap]
    type = ADMatReaction<<<{"description": "Kernel representing the contribution of the PDE term $-L*v$, where $L$ is a reaction rate material property, $v$ is a scalar variable (nonlinear or coupled), and whose Jacobian contribution is calculated using automatic differentiation.", "href": "../../source/kernels/ADMatReaction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
    reaction_rate<<<{"description": "The reaction_rate used with the kernel."}>>> = ${fparse scale_time * detrapping_rate_0} # from 1/s to non-dimensional
  []
  [Trapping_trap]
    type = ADMatCoupledDefectAnnihilation<<<{"description": "Kernel to add K*v*(u0-u), where K=annihilation rate, u=variable, u0=equilibrium, v=coupled variable", "href": "../../source/kernels/ADMatCoupledDefectAnnihilation.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
    eq_concentration<<<{"description": "The equilibrium concentration of the variable used with the kernel"}>>> = ${tritium_trap_0_scaled}
    v<<<{"description": "Set this to make v a coupled variable, otherwise it will use the kernel's nonlinear variable for v"}>>> = tritium_s
    annihilation_rate<<<{"description": "The annihilation rate used with the kernel"}>>> = ${fparse scale_time*scale_quantity/scale_length^3*trapping_rate_0} # from um^3/s/mol to non-dimensional
  []
  # ========================= Surface Reactions ================================
  # ======================================= Diatomic molecule dissolve at surface
  # T2(g) -> 2 T(s) rate K [T2] (1-theta) ======================================
  [Surface_Reaction_dissolution_tritium_2g] # T2(g) -> 2 T(s) rate -K*tritium_2g
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_2g'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = -1
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_gXYs
  []
  [Surface_Reaction_dissolution_tritium_2g_tritium_s] # T2(g) -> 2 T(s) rate 2*K*tritium_2g
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_2g'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = 2
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_gXYs
  []
  # ======================================= Diatomic molecule combine at surface
  # 2 T(s) -> T2(g) rate K [T]^2 ===============================================
  [Surface_Reaction_release_tritium_s_tritium_2g] # 2 T(s) -> T2(g) rate -2*K*tritium_s^2
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_s tritium_s'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = -2
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_sXYg
  []
  [Surface_Reaction_formation_tritium_2g] # 2 T(s) -> T2(g) rate K*tritium_s^2
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_s tritium_s'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = 1
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_sXYg
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [pore_aux]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = pore
    function<<<{"description": "The function to use as the value"}>>> = 'pore_position_func'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL'
  []
[]

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  #===================================================== Interpolation functions
  [hsurfaceAD] # equal to 1 at pore surface, 0 elsewhere. The values of hsurface should be smaller than the ones for hmat or hpore to allow the diffusion of the species out of the surface
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '16*(1-pore)^2*pore^2'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hsurfaceAD'
  []
  [hmatAD] # equal to 1 in material, and 0 in pore.
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '(1-pore)*(1-pore)'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hmatAD'
  []
  [hporeAD] # equal to 1 in pore, and 1 in material. # notice that hpore overlaps a bit with hmat to prevent some agglomeration of chemicals on the edges of the surface.
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'pore*pore'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hporeAD'
  []
  #================================================= Local species concentration
  [Tritium_ceramic]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'tritium_ceramic'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'tritium_s tritium_trap'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'tritium_s+tritium_trap'
  []
  [Tritium_pore]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'tritium_pore'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'tritium_2g'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '2*tritium_2g'
  []
  [Tritium_total]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'tritium_tot'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'tritium_s tritium_trap tritium_2g'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'tritium_s+tritium_trap+2*tritium_2g'
  []
  #=============================================== Occupied sites on the surface
  [theta]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'theta'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'tritium_s tritium_trap'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '(tritium_s + tritium_trap)/(${available_sites_mult}*${tritium_trap_0_scaled})'
  []
  #====================================================== Diffusion coefficients
  [Diffusion_coefficient_D] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'D_tritium'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hmatAD(pore)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}/${scale_length}^2*(hmatAD*${Diffusion_coefficient_Db_D0} + (1-hmatAD)*${Diffusion_min})'
    derivative_order<<<{"description": "Maximum order of derivatives taken"}>>> = 2
  []
  [Diffusion_coefficient_D_pore] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'D_pore'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hporeAD(pore)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}/${scale_length}^2*(hporeAD*${Diffusion_coefficient_D_pore_D0} + (1-hporeAD)*${Diffusion_min})'
  []
  #============================================ Surface reaction rate of tritium
  [ReactionRateSurface_sXYg_ref] # from surface to gaseous for diatomic molecules # from um^3/s/mol to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = ReactionRateSurface_sXYg
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hsurfaceAD(pore)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}*${scale_quantity}/${scale_length}^3*${reactionRateSurface_sXYg_ref_0}*hsurfaceAD'
  []
  [ReactionRateSurface_gXYs_ref] # from gaseous to surface for diatomic molecules # from 1/s to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = ReactionRateSurface_gXYs
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore tritium_s tritium_trap'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hsurfaceAD(pore) theta(tritium_s,tritium_trap)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}*${reactionRateSurface_gXYs_ref_0}*hsurfaceAD*(1-theta)^2'
  []
[]

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [Tritium_amount_free]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_s
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_trapped]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_trap
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_ceramic]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_ceramic
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium2g_amount]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_2g
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_pore]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_pore
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_total]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_tot
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Surface_tot]
    type = ElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = 1
  []
  [dt]
    type = TimestepSize<<<{"description": "Reports the timestep size", "href": "../../source/postprocessors/TimestepSize.html"}>>>
  []
[]

[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
  [Newtonlu]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    solve_type<<<{"description": "PJFNK: Preconditioned Jacobian-Free Newton Krylov JFNK: Jacobian-Free Newton Krylov NEWTON: Full Newton Solve FD: Use finite differences to compute Jacobian LINEAR: Solving a linear problem"}>>> = 'NEWTON'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -sub_pc_type -snes_type'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'asm      lu           vinewtonrsls' # This petsc option helps prevent negative concentrations'
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  scheme = bdf2

  nl_rel_tol = 1e-10
  end_time = 2000
  dtmax = 50
  nl_max_its = 14
  [TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    optimal_iterations = 12
    iteration_window = 1
    growth_factor = 1.2
    dt = 2.37e-7
    cutback_factor = 0.75
    cutback_factor_at_failure = 0.75
  []
[]

[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  [exodus]
    type = Exodus<<<{"description": "Object for output data in the Exodus format", "href": "../../source/outputs/Exodus.html"}>>>
    time_step_interval<<<{"description": "The interval (number of time steps) at which output occurs. Unless explicitly set, the default value of this parameter is set to infinity if the wall_time_interval is explicitly set."}>>> = 5
  []
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = ${output_name}
[]

This section describes this input file, what it does, and how to run it. Note that the input file also contains extensive comments meant to help users.

Background and Governing Equations

In this step, we model tritium transport in the pores, surface reactions at the ceramics interface, and diffusion, trapping, and detrapping within the ceramics.

The evolution of tritium concentration in solid solution within the ceramic material, $var element = document.getElementById("moose-equation-3b46fc3f-0f91-4b44-ba78-30c3ff896992");katex.render("c_s", element, {displayMode:false,throwOnError:false});$ (mol⋅μm $var element = document.getElementById("moose-equation-cc77b69b-febb-43e5-a172-d2b9ab4f3124");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$ ), is governed by the following equation:

$(1)var element = document.getElementById("moose-equation-3a12a573-78b7-43ee-b66f-18085375341e");katex.render(" \\frac{\\partial c_s}{\\partial t} = \\nabla \\cdot (h_c D_c \\nabla c_s) - \\frac{\\partial c_t}{\\partial t} + 2(K_{gs} - K_{sg}) h_s", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-33ef57e7-8c73-4499-b572-e057b64d932d");katex.render("t", element, {displayMode:false,throwOnError:false});$ is time (s), $var element = document.getElementById("moose-equation-10de42f9-0f94-44dc-857c-b9b162033e52");katex.render("D_c", element, {displayMode:false,throwOnError:false});$ is the diffusion coefficient of tritium in the ceramic material ( $var element = document.getElementById("moose-equation-5b89dff8-1ffe-4d96-b3b4-e86b9c8866a6");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-693b1268-810c-4a23-9e88-9087c031082e");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s), $var element = document.getElementById("moose-equation-7d7a961a-bf00-4be0-9bc3-ef8100b451c5");katex.render("c_t", element, {displayMode:false,throwOnError:false});$ is the concentration of trapped tritium (mol/ $var element = document.getElementById("moose-equation-04b65746-8075-4608-9e5f-9166e53b007a");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-94453c76-4398-4f19-9d0c-f32d9b7d414f");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ ), $var element = document.getElementById("moose-equation-17f94b7e-7b64-41c6-ab3d-15d1a9cf5eb8");katex.render("K_{gs}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-2d929c02-6d7d-4819-a47f-be30f7cb890c");katex.render("K_{sg}", element, {displayMode:false,throwOnError:false});$ are the reaction rates at the surface (mol/ $var element = document.getElementById("moose-equation-0edca083-9298-49a4-9d55-cb300c636fd5");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-f399ece5-491b-43b1-8758-0efcec67f789");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ /s), and $var element = document.getElementById("moose-equation-deba5a5a-13c3-4adb-bfb6-677f61823535");katex.render("h_c", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-2bc7fd2c-4af3-44cb-95ac-c57f5e24b758");katex.render("h_s", element, {displayMode:false,throwOnError:false});$ are dimensionless interpolation functions that delimit the ceramic material and pore surface, respectively (see the next section below).

The concentration of trapped tritium is governed by:

$(2)var element = document.getElementById("moose-equation-d489633f-3518-42b6-aba4-1e5cdc1ff1f4");katex.render(" \\frac{\\partial c_t}{\\partial t} = K_T c_s (c_t^0 - c_t) - K_D c_t", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-b0ebb629-0280-43d9-9496-faca4c605ca0");katex.render("K_T", element, {displayMode:false,throwOnError:false});$ is the rate constant for trapping ( $var element = document.getElementById("moose-equation-df0da395-bc1a-4dbe-9f59-e9a7de3fca9c");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-a97fdc2a-4113-47e0-bf87-c03bda56c569");katex.render("^3", element, {displayMode:false,throwOnError:false});$ /s/mol), $var element = document.getElementById("moose-equation-e263540b-28bb-4a03-824c-6aa239ddc818");katex.render("c_t^0", element, {displayMode:false,throwOnError:false});$ is the density of available trapping sites (mol/ $var element = document.getElementById("moose-equation-83e15bd1-5c02-48d1-8d25-2a33a9d625af");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-1efc5e6d-02b9-4872-8bcd-9d9963837f94");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ ), and $var element = document.getElementById("moose-equation-34265ffc-0f93-44d9-bd5b-f7b5285d8dbd");katex.render("K_D", element, {displayMode:false,throwOnError:false});$ is the rate constant for detrapping (1/s).

The evolution of tritium concentration in pores and outside the ceramic ( $var element = document.getElementById("moose-equation-de012bf2-4478-48fa-b925-da005e877e3a");katex.render("c_2", element, {displayMode:false,throwOnError:false});$ ) is defined by:

$(3)var element = document.getElementById("moose-equation-fac5ca85-cfe4-4033-bfc3-2d96f91cea4d");katex.render(" \\frac{\\partial c_2}{\\partial t} = \\nabla \\cdot (h_p D_p \\nabla c_2) + (K_{sg} - K_{gs}) h_s", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-e6eba011-7d4e-4d86-9129-f780d09902d2");katex.render("D_p", element, {displayMode:false,throwOnError:false});$ is the diffusion coefficient of $var element = document.getElementById("moose-equation-ee2d6212-580d-4011-a930-aa120554abb8");katex.render("T_2", element, {displayMode:false,throwOnError:false});$ in the pores ( $var element = document.getElementById("moose-equation-aa4422a3-4dc1-4023-8812-c171d4efef08");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-44ef96e4-99f8-4498-814d-348080d6e910");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s) and $var element = document.getElementById("moose-equation-8d46eac2-b139-40c6-9a2d-fe50c75116f0");katex.render("h_p", element, {displayMode:false,throwOnError:false});$ is a dimensionless interpolation function that delimits the pore.

These concentration variables are defined in the input file as:

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [tritium_s]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = ${tritium_s_scaling}
  []
  [tritium_trap]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = ${tritium_trap_scaling}
  []
  [tritium_2g]
  []
[]

The equations are encoded into the input file as kernels, which each implement a specific term from the equations

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  # ========================================== Time derivative for each variable
  [dtritium_sdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_s'
  []
  [dtritium_2gdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_2g'
  []
  [dtritium_trapdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
  []
  # ============================================================= Bulk diffusion
  [Diff_tritium]
    type = ADMatDiffusion<<<{"description": "Diffusion equation kernel that takes an isotropic diffusivity from a material property", "href": "../../source/kernels/ADMatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = D_tritium
  []
  # ============================================================= Pore diffusion
  [Diff_tritium2g]
    type = ADMatDiffusion<<<{"description": "Diffusion equation kernel that takes an isotropic diffusivity from a material property", "href": "../../source/kernels/ADMatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = D_pore
  []
  # ================================================ Tritium trapping/detrapping
  [Trapping]
    type = CoupledTimeDerivative<<<{"description": "Time derivative Kernel that acts on a coupled variable. Weak form: $(\\psi_i, \\frac{\\partial v_h}{\\partial t})$.", "href": "../../source/kernels/CoupledTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    v<<<{"description": "Coupled variable"}>>> = 'tritium_trap'
  []
  [Detrapping_trap]
    type = ADMatReaction<<<{"description": "Kernel representing the contribution of the PDE term $-L*v$, where $L$ is a reaction rate material property, $v$ is a scalar variable (nonlinear or coupled), and whose Jacobian contribution is calculated using automatic differentiation.", "href": "../../source/kernels/ADMatReaction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
    reaction_rate<<<{"description": "The reaction_rate used with the kernel."}>>> = '${fparse scale_time * detrapping_rate_0}' # from 1/s to non-dimensional
  []
  [Trapping_trap]
    type = ADMatCoupledDefectAnnihilation<<<{"description": "Kernel to add K*v*(u0-u), where K=annihilation rate, u=variable, u0=equilibrium, v=coupled variable", "href": "../../source/kernels/ADMatCoupledDefectAnnihilation.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
    eq_concentration<<<{"description": "The equilibrium concentration of the variable used with the kernel"}>>> = ${tritium_trap_0_scaled}
    v<<<{"description": "Set this to make v a coupled variable, otherwise it will use the kernel's nonlinear variable for v"}>>> = tritium_s
    annihilation_rate<<<{"description": "The annihilation rate used with the kernel"}>>> = '${fparse scale_time*scale_quantity/scale_length^3*trapping_rate_0}' # from um^3/s/mol to non-dimensional
  []
  # ========================= Surface Reactions ================================
  # ======================================= Diatomic molecule dissolve at surface
  # T2(g) -> 2 T(s) rate K [T2] (1-theta) ======================================
  [Surface_Reaction_dissolution_tritium_2g]
    # T2(g) -> 2 T(s) rate -K*tritium_2g
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_2g'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = -1
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_gXYs
  []
  [Surface_Reaction_dissolution_tritium_2g_tritium_s]
    # T2(g) -> 2 T(s) rate 2*K*tritium_2g
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_2g'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = 2
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_gXYs
  []
  # ======================================= Diatomic molecule combine at surface
  # 2 T(s) -> T2(g) rate K [T]^2 ===============================================
  [Surface_Reaction_release_tritium_s_tritium_2g]
    # 2 T(s) -> T2(g) rate -2*K*tritium_s^2
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_s tritium_s'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = -2
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_sXYg
  []
  [Surface_Reaction_formation_tritium_2g]
    # 2 T(s) -> T2(g) rate K*tritium_s^2
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_s tritium_s'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = 1
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_sXYg
  []
[]

Reaction Rates

The reaction rates for tritium dissociation and formation on the ceramic surface are:

$(4)var element = document.getElementById("moose-equation-edf40e6e-39e9-45fb-9c0b-1c506b852fc1");katex.render(" K_{gs} = k_{gs} c_2 (1 - \\theta)^2,", element, {displayMode:true,throwOnError:false});$ and $(5)var element = document.getElementById("moose-equation-3a315659-278c-4879-ae9e-0f280c710387");katex.render(" K_{sg} = k_{sg} c_s^2,", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-0f52cabd-25f3-44d5-be95-4b869b970fca");katex.render("k_{gs}", element, {displayMode:false,throwOnError:false});$ (1/s) and $var element = document.getElementById("moose-equation-6b2b55ba-5230-46b5-8905-7e4b1e9bac94");katex.render("k_{sg}", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-c7f65590-d6d0-4d07-99ff-b0b923416a74");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-6e9b1a71-11a2-4af5-8512-cfde115c21f1");katex.render("^3", element, {displayMode:false,throwOnError:false});$ /s/mol) are reaction rate constants, and $var element = document.getElementById("moose-equation-60bdb88d-7fd1-4203-9a9e-cf12e6bf156f");katex.render("\\theta", element, {displayMode:false,throwOnError:false});$ is the fraction of available sites for $var element = document.getElementById("moose-equation-e2b28420-0138-4b49-8efd-4e556ce5bb56");katex.render("T_s", element, {displayMode:false,throwOnError:false});$ , defined as:

$(6)var element = document.getElementById("moose-equation-ecee3940-9148-4654-ad00-a107c217ba14");katex.render(" \\theta = \\frac{c_s + c_t}{\\theta_0 c_t^0},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-e7b84d43-7343-4003-b0d9-de9c530f30bf");katex.render("\\theta_0 \\geq 1", element, {displayMode:false,throwOnError:false});$ is a dimensionless multiplier that ensures the concentration of available sites for tritium exceeds the concentration of trapping sites.

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [ReactionRateSurface_sXYg_ref]
    # from surface to gaseous for diatomic molecules # from um^3/s/mol to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = ReactionRateSurface_sXYg
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hsurfaceAD(pore)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}*${scale_quantity}/${scale_length}^3*${reactionRateSurface_sXYg_ref_0}*hsurfaceAD'
  []

  [ReactionRateSurface_gXYs_ref]
    # from gaseous to surface for diatomic molecules # from 1/s to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = ReactionRateSurface_gXYs
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore tritium_s tritium_trap'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hsurfaceAD(pore) theta(tritium_s,tritium_trap)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}*${reactionRateSurface_gXYs_ref_0}*hsurfaceAD*(1-theta)^2'
  []
[]

Interpolation Functions

The microstructure is defined by a dimensionless order parameter $var element = document.getElementById("moose-equation-21b24f83-ad81-4d43-8ff3-08009536c11a");katex.render("\\eta", element, {displayMode:false,throwOnError:false});$ , which equals 1 in the pores, 0 in the ceramic, and varies continuously at the pore surface. This is commonly used in phase field modeling to describe continuous microstructures. The interpolation functions, which help ensure that diffusion coefficients (or other material properties) are restricted to the appropriate domains, are defined as:

$(7)var element = document.getElementById("moose-equation-7b9cf6c3-2dce-49f5-a3b3-317b6a90bb6f");katex.render(" h_c = (1 - \\eta)^2", element, {displayMode:true,throwOnError:false});$

for the ceramic material,

$(8)var element = document.getElementById("moose-equation-5a3d5eea-a589-440e-ac10-91c2b4afe6a1");katex.render(" h_p = \\eta^2", element, {displayMode:true,throwOnError:false});$

for the pore area, and

$(9)var element = document.getElementById("moose-equation-4bddbd56-ebf8-48d6-b404-fcfca2f24646");katex.render(" h_s = 16 \\eta^2 (1 - \\eta)^2", element, {displayMode:true,throwOnError:false});$

to describe the interface between the ceramic material and the pore. Note that the factor 16 ensures that at the middle of the interface, $var element = document.getElementById("moose-equation-678e892a-9a1b-47fc-9d37-da80c28ad80c");katex.render("h_s(\\eta = 0.5) = 1", element, {displayMode:false,throwOnError:false});$ , as desired.

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [hmatAD]
    # equal to 1 in material, and 0 in pore.
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '(1-pore)*(1-pore)'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hmatAD'
  []

  [hporeAD]
    # equal to 1 in pore, and 1 in material. # notice that hpore overlaps a bit with hmat to prevent some agglomeration of chemicals on the edges of the surface.
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'pore*pore'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hporeAD'
  []

  [hsurfaceAD]
    # equal to 1 at pore surface, 0 elsewhere. The values of hsurface should be smaller than the ones for hmat or hpore to allow the diffusion of the species out of the surface
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '16*(1-pore)^2*pore^2'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hsurfaceAD'
  []
[]

Dimensionless Modeling

To improve convergence, the model is rendered dimensionless using characteristic time $var element = document.getElementById("moose-equation-e5a77532-db1f-4c0f-b1d8-d1bae3408814");katex.render("t_c=1", element, {displayMode:false,throwOnError:false});$ s, length $var element = document.getElementById("moose-equation-266f4b6c-cd7c-4122-8b1b-e4fb01974d3b");katex.render("l_c=1", element, {displayMode:false,throwOnError:false});$ μm, and quantity $var element = document.getElementById("moose-equation-511a1895-c84c-4f23-a5b3-0142f53cdfb1");katex.render("n_c=10^{-18}", element, {displayMode:false,throwOnError:false});$ mol.

Calibration Against Experimental Data

In Simon et al. (2022), the one-dimensional version of this case is used to calibrate the model against experimental data. In this example, like in Section V of Simon et al. (2022), we use this calibrated model on the two-dimensional microstructure.

The resulting model parameter values are summarized in Table 1.

Table 1: Calibrated model’s parameter values from Simon et al. (2022).

Model Parameter	Calibrated Value	Units
$var element = document.getElementById("moose-equation-3e1b2fd8-9db2-4ffc-8a4d-977f3354a84e");katex.render("D_c", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-358c0d30-21f5-4b67-b5a2-c03ef809e336");katex.render("1.65 \\times 10^6", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-a016bb53-3160-497f-817b-4fd64de7d4a5");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-0e58ec11-e3e8-4c49-95cd-172e8add1d69");katex.render("^2\\cdot", element, {displayMode:false,throwOnError:false});$ s $var element = document.getElementById("moose-equation-6c3050a2-2fcc-4553-849f-9650d037146f");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-11fab172-ab8f-4d12-934a-a13c31d23ae5");katex.render("D_p", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-58dc85ef-1b68-40ba-82c6-961250d0a382");katex.render("2.25 \\times 10^8", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-4b0e993b-b813-447b-a5ec-20fea5211dea");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-3d23fbf1-330c-41c6-8e52-2868bcff0c0d");katex.render("^2\\cdot", element, {displayMode:false,throwOnError:false});$ s $var element = document.getElementById("moose-equation-4bcd5499-56cc-4714-8a8b-3a1702c3a200");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-b494902d-4d58-4859-af05-6ea9cefe646b");katex.render("K_D", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-554a7836-6a2e-4f0d-8eca-228b975bcc47");katex.render("4.06", element, {displayMode:false,throwOnError:false});$	s $var element = document.getElementById("moose-equation-9928033a-45c8-47f1-b37e-ba589db48881");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-db8dd8ae-8834-40b5-a718-4223646dd54c");katex.render("K_T", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-97fe7ab0-1309-4268-805e-9e86c0949897");katex.render("8.33 \\times 10^{16}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-695cf6ff-c9e5-4f07-ae0a-771d0c11036f");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-b7f33dc1-5092-4a65-88dd-5edd1c5f5342");katex.render("^{3} \\cdot", element, {displayMode:false,throwOnError:false});$ s $var element = document.getElementById("moose-equation-253dcd94-69c5-421d-97ab-6932fb290600");katex.render("^{-1}\\cdot", element, {displayMode:false,throwOnError:false});$ mol $var element = document.getElementById("moose-equation-ead38b28-4c05-4b3d-ba33-cba59a983165");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-c948ee85-6638-412a-97bb-95f6871629e5");katex.render("c_t^0", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-9d0caea5-f2f3-4285-8902-b7d1627f7956");katex.render("3.47 \\times 10^{-18}", element, {displayMode:false,throwOnError:false});$	mol $var element = document.getElementById("moose-equation-a2a85eef-0d45-452a-b38c-04bd53d56494");katex.render("\\cdot \\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-f7eacda1-bfc5-4eea-b98d-8f7b41a7a33e");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-38945643-e722-4f04-a723-9e84136aaaae");katex.render("\\theta^0", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-d1308e20-1108-4fa9-9d17-6766670cfdf1");katex.render("4.31", element, {displayMode:false,throwOnError:false});$	(-)
$var element = document.getElementById("moose-equation-7c62e4b0-8054-4074-85e3-66d10cddaf8a");katex.render("k_{sg}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-8056ac0b-a737-4dfd-a7bc-33b6bb82c8ef");katex.render("1.55 \\times 10^{15}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-130b73f0-e572-4fce-ab5e-fe506a0cce35");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-dcbc86f0-28c9-466c-a288-d7b94f72f431");katex.render("^3\\cdot", element, {displayMode:false,throwOnError:false});$ s $var element = document.getElementById("moose-equation-ad2b0676-e1d8-44f4-bae8-d2a0129f4d0e");katex.render("^{-1}\\cdot", element, {displayMode:false,throwOnError:false});$ mol $var element = document.getElementById("moose-equation-7e613662-191b-40a2-89ee-08b5de9872bf");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-e4d9ae88-9389-4081-b594-754671c2a165");katex.render("k_{gs}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-05a7f4bf-ec31-4c13-82c6-aa41dc2e546c");katex.render("34.69", element, {displayMode:false,throwOnError:false});$	s $var element = document.getElementById("moose-equation-00209fff-9a4e-4e0e-8ae4-dfb53fd71c30");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$

Numerical Method, boundary conditions, and initial condition

As shown in Figure 4, the simulated domain represents a quarter of the sample, whose radius is 4.5 mm.

The same numerical methods as described in Simon et al. (2022) are used here. The numerical methods, including time step, time adaptivity, time integration method, preconditioning, and others are defined in the Preconditioning and Executioner blocks, as:

[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
  [Newtonlu]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    solve_type<<<{"description": "PJFNK: Preconditioned Jacobian-Free Newton Krylov JFNK: Jacobian-Free Newton Krylov NEWTON: Full Newton Solve FD: Use finite differences to compute Jacobian LINEAR: Solving a linear problem"}>>> = 'NEWTON'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -sub_pc_type -snes_type'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'asm      lu           vinewtonrsls' # This petsc option helps prevent negative concentrations'
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  scheme = bdf2

  nl_rel_tol = 1e-10
  end_time = 2000
  dtmax = 50
  nl_max_its = 14
  [TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    optimal_iterations = 12
    iteration_window = 1
    growth_factor = 1.2
    dt = 2.37e-7
    cutback_factor = 0.75
    cutback_factor_at_failure = 0.75
  []
[]

Concerning boundary conditions, we apply no-flux boundary conditions except for $var element = document.getElementById("moose-equation-ac28a96c-0bd4-4b1e-b7f5-159a6ec78115");katex.render("c_2", element, {displayMode:false,throwOnError:false});$ on the right-hand and top sides, which is fixed as a Dirichlet condition according to:

$(10)var element = document.getElementById("moose-equation-23a9861f-dd09-4375-80b6-41ce6500cb7a");katex.render(" c_{s,BC,pore} = \\frac{P}{RT}", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-05a4527b-83e6-4016-b261-fb3c3d2ee172");katex.render("P", element, {displayMode:false,throwOnError:false});$ is the tritium pressure (Pa), $var element = document.getElementById("moose-equation-36745618-c390-450f-a2eb-5adf8fec3015");katex.render("R", element, {displayMode:false,throwOnError:false});$ is the gas constant from PhysicalConstants, and $var element = document.getElementById("moose-equation-327c4bb3-f7ca-472a-bb0d-44385a406ce8");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature (K).

[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [tritium_2g_sides_d]
    # Fix concentration on the sides
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'right top'
    value<<<{"description": "Value of the BC"}>>> = ${tritium_2g_concentration}
  []
[]

As initial conditions, the ceramic initially contains no tritium, and the pore is filled with $var element = document.getElementById("moose-equation-bb9fe5bb-6814-4140-b4af-c5782296108f");katex.render("T_2", element, {displayMode:false,throwOnError:false});$ with the value from Eq. (10).

[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [tritium_2g_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = tritium_2g
    function<<<{"description": "The initial condition function."}>>> = 'gaseous_concentration_1_function_ic'
  []
[]

[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
  [pore_position_func]
    type = SolutionFunction<<<{"description": "Function for reading a solution from file.", "href": "../../source/functions/SolutionFunction.html"}>>>
    solution<<<{"description": "The SolutionUserObject to extract data from."}>>> = initial_microstructure
    from_variable<<<{"description": "The name of the variable in the file that is to be extracted"}>>> = gr0
  []
  [gaseous_concentration_1_function_ic]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'len_pore  position_pore_int'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = '25   4521'
    expression<<<{"description": "The user defined function."}>>> = '${tritium_2g_concentration}*0.5*(1.0+tanh((sqrt(x*x+y*y)-position_pore_int)/(len_pore/2.2)))'
  []
[]

Run simulations

To run the input file to reproduce the absorption results shown in Fig. 6 in Simon et al. (2022), users need to first run the previous simulations from the first step. As a result, TMAP8 will have generated pore_structure_open/pore_structure_open.e and pore_structure_close/pore_structure_close.e. Then, to model tritium transport in the microstructure with closed pores, run:


~/projects/TMAP8/test/tests/pore_scale_transport/
mpirun -np 4 ~/projects/TMAP8/tmap8-opt -i 2D_absorption_base.i pore_structure_closed_absorption.params

and to model tritium transport in the microstructure with open pores, run:


~/projects/TMAP8/test/tests/pore_scale_transport/
mpirun -np 4 ~/projects/TMAP8/tmap8-opt -i 2D_absorption_base.i pore_structure_open_absorption.params

Note that pore_structure_closed_absorption.params and pore_structure_open_absorption.params complement the base input file 2D_absorption_base.i by providing input and output names.

Results and Discussion

Figure 6 in Simon et al. (2022) shows the tritium absorption profiles for samples with different pore interconnectivities. The results indicate that higher pore interconnectivity leads to faster tritium transport, which is the type of insight that one-dimensional simulations cannot provide.

These outputs are obtained via the postprocessor blocks:

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [Tritium_amount_free]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_s
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_trapped]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_trap
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_ceramic]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_ceramic
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium2g_amount]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_2g
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_pore]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_pore
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_total]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_tot
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Surface_tot]
    type = ElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = 1
  []
  [dt]
    type = TimestepSize<<<{"description": "Reports the timestep size", "href": "../../source/postprocessors/TimestepSize.html"}>>>
  []
[]

Note that these simulations are not run to generate Fig. 6 from Simon et al. (2022) on-the-fly in the documentation like it is custom to do in other examples and V&V cases due to the significant time these simulations can take when using a fine mesh. However, using the instructions and details provided above, users can replicate the results presented in Fig. 6 from Simon et al. (2022) and reproduced in Figure 1. If any challenge arise or if you have a question while working through this example case, feel free to reach out to the TMAP8 development team on the TMAP8 GitHub discussion page.

Complete input file

Below are the complete input files, which can be run reliably with approximately 4 processor cores. Note that this input file has not been optimized to reduce computational costs.

Input file used to import the image and create the smooth microstructure:

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Info:
# - This simulation predicts phase evolution based on a grain growth phase field model.
# - The pore and the ceramics represent two different phases, and the initial conditions
#   are provided with a black and white picture, which is read by ImageFunction.
# - The simulation smoothens the interface between the two phases.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_microstructure_reader_smoothing_base.i pore_structure_open.params

# mesh information
num_nodes_x = ${fparse 2*120} # (-)
num_nodes_y = ${fparse 2*120} # (-)
domain_start_x = ${units 0 mum}
domain_start_y = ${units 0 mum}
domain_end_x = ${units 5425 mum}
domain_end_y = ${units 5425 mum}

# grain growth parameters (since the point of the simulation is only to get a smooth interface, the model parameters can be selected arbitrarily)
op_num = 2 # Number of grains (-)
mobility_prefactor = ${units 2.5e-6 mum^4/J/s}
GB_energy = ${units 0.7 J/mum^2}
activation_energy = ${units 0.23 eV}
temperature = ${units 700 K}
width_diffuse_GB = ${units 50 mum} # Width of the diffuse GB

# image function option
threshold_image_function = 255.5 # (-)

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  [gen]
    type = DistributedRectilinearMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/DistributedRectilinearMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_nodes_x}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_nodes_y}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${domain_start_x}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${domain_end_x}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = ${domain_start_y}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${domain_end_y}
  []
[]

[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
  op_num = ${op_num}
  var_name_base = gr # Base name of grains
[]

[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
  [grain_tracker]
    type = GrainTracker<<<{"description": "Grain Tracker object for running reduced order parameter simulations without grain coalescence.", "href": "../../source/postprocessors/GrainTracker.html"}>>>
    flood_entity_type<<<{"description": "Determines whether the flood algorithm runs on nodes or elements"}>>> = ELEMENTAL
    compute_halo_maps<<<{"description": "Instruct the Postprocessor to communicate proper halo information to all ranks"}>>> = true
  []
[]

[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
  [image_func0]
    type = ImageFunction<<<{"description": "Function with values sampled from an image or image stack.", "href": "../../source/functions/ImageFunction.html"}>>>
    file<<<{"description": "Name of single image file to extract mesh parameters from.  If provided, a 2D mesh is created."}>>> = ${input_name}
    threshold<<<{"description": "The threshold value"}>>> = ${threshold_image_function}
    upper_value<<<{"description": "The value to set for data greater than the threshold value"}>>> = 1
    lower_value<<<{"description": "The value to set for data less than the threshold value"}>>> = 0
  []
  [image_func1]
    type = ImageFunction<<<{"description": "Function with values sampled from an image or image stack.", "href": "../../source/functions/ImageFunction.html"}>>>
    file<<<{"description": "Name of single image file to extract mesh parameters from.  If provided, a 2D mesh is created."}>>> = ${input_name}
    threshold<<<{"description": "The threshold value"}>>> = ${threshold_image_function}
    upper_value<<<{"description": "The value to set for data greater than the threshold value"}>>> = 0
    lower_value<<<{"description": "The value to set for data less than the threshold value"}>>> = 1
  []
[]

[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [gr0_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    function<<<{"description": "The initial condition function."}>>> = image_func0
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = gr0
  []
  [gr1_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    function<<<{"description": "The initial condition function."}>>> = image_func1
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = gr1
  []
[]

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [PolycrystalVariables<<<{"href": "../../syntax/Variables/PolycrystalVariables/index.html"}>>>]
    # Custom action that creates all of the phase variables and sets their initial condition
  []
[]

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [bnds]
    # Variable used to visualize the interfaces in the simulation
  []
  [unique_grains]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
[]

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [PolycrystalKernel<<<{"href": "../../syntax/Kernels/PolycrystalKernel/index.html"}>>>]
    # Custom action creating all necessary kernels for phase evolution.  All input parameters are up in GlobalParams
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [bnds_aux]
    type = BndsCalcAux<<<{"description": "Calculate location of grain boundaries in a polycrystalline sample", "href": "../../source/auxkernels/BndsCalcAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = bnds
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'initial timestep_end'
  []
  [unique_grains]
    type = FeatureFloodCountAux<<<{"description": "Feature detection by connectivity analysis", "href": "../../source/auxkernels/FeatureFloodCountAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = unique_grains
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    flood_counter<<<{"description": "The FeatureFloodCount UserObject to get values from."}>>> = grain_tracker
    field_display<<<{"description": "Determines how the auxilary field should be colored. (UNIQUE_REGION and VARIABLE_COLORING are nodal, CENTROID is elemental, default: UNIQUE_REGION)"}>>> = UNIQUE_REGION
  []
[]

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [InterfaceEvolution]
    type = GBEvolution<<<{"description": "Computes necessary material properties for the isotropic grain growth model", "href": "../../source/materials/GBEvolution.html"}>>>
    GBmob0<<<{"description": "Grain boundary mobility prefactor in m^4/(J*s)"}>>> = ${mobility_prefactor}
    GBenergy<<<{"description": "Grain boundary energy in J/m^2"}>>> = ${GB_energy}
    Q<<<{"description": "Grain boundary migration activation energy in eV"}>>> = ${activation_energy}
    T<<<{"description": "Temperature in Kelvin"}>>> = ${temperature}
    wGB<<<{"description": "Diffuse GB width in the length scale of the model"}>>> = ${width_diffuse_GB}
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  scheme = bdf2
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  l_max_its = 20
  nl_max_its = 9
  nl_rel_tol = 1e-10
  [TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    cutback_factor = 0.9
    dt = 3e-2
    growth_factor = 1.1
    optimal_iterations = 7
  []
  dtmax = 25

  start_time = 0.0
  end_time = 100
[]

[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = ${output_name}
  [exodus]
    type = Exodus<<<{"description": "Object for output data in the Exodus format", "href": "../../source/outputs/Exodus.html"}>>>
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'FINAL'
  []
[]

Input file utilizing the smooth microstructure and modeling pore scale tritium transport during absorption:

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Phase Field Model:   Isotropic bulk diffusion, pore surface reactions, isotropic pore diffusion
# Type:                Transient
# Phase structure:     Contains two phases: the bulk of a material, and the pores. An interface separates the two
# Physics:             tritium trapping and detrapping, surface reactions, diffusion
# Species variable:    tritium_s (Ts), tritium_trap (trapped T), tritium_2g (T2(g))
# Species AuxVariable: --
# Chemistry:           at the pore surface: (s) means in the solid, (g) means in gas form
#    Diatomic molecule dissolving at surface
#                                 T2(g) -> 2 T(s) rate K [T2] (1-theta)^2
#    Diatomic molecule combining at surface
#                                 2 T(s) -> T2(g) rate K [T]^2
#
# BCs:                 No flux for tritium and tritium_trap, zero at pore center for gaseous tritium in every form
#
#
# Info:
# - This input file uses the microstructure generated by 2D_microstructure_reader_smoothing_base.i
#   and performs tritium transport simulation on this microstructure.
# - Basic 2D non-dimensional input file for the transport of a solute in a microstructure with pores.
# - The goal of this input base and associated parameter files is to have the same simulation with different pore configurations
#   to show the effect of microstructure on tritium release.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_absorption_base.i pore_structure_open_absorption.params

# Scaling factors
scale_quantity = ${units 1e-18 moles}
scale_time = ${units 1 s}
scale_length = ${units 1 mum}
tritium_trap_scaling = 1e3 # (-)
tritium_s_scaling = 1e2 # (-)

# Conditions
tritium_2g_concentration = ${units 0.1121 mol/mum^3}

# Material properties
tritium_trap_0 = ${units 3.472805925e-18 mol/mum^3}
tritium_trap_0_scaled = ${fparse tritium_trap_0*scale_length^3/scale_quantity} # non-dimensional
available_sites_mult = 4.313540691 # (-) should be >1, otherwise there are fewer 'available sites' than 'trapping sites'
Diffusion_min = ${units 0 mum^2/s} # used as diffusion coefficients where I don't want the species to move (more stable than 0)
Diffusion_coefficient_Db_D0 = ${units 1649969.728 mum^2/s}
Diffusion_coefficient_D_pore_D0 = ${units 224767738.6 mum^2/s}
detrapping_rate_0 = ${units -4.065104516 1/s} # should be negative
trapping_rate_0 = ${units 8.333146645e+16 1/s/mol} # should be positive
reactionRateSurface_sXYg_ref_0 = ${units 1549103878000000.0 1/s/mol}
reactionRateSurface_gXYs_ref_0 = ${units 34.69205021 1/s}

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  file = ${input_name}
[]

[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
  [initial_microstructure]
    type = SolutionUserObject<<<{"description": "Reads a variable from a mesh in one simulation to another", "href": "../../source/userobjects/SolutionUserObject.html"}>>>
    mesh<<<{"description": "The name of the mesh file (must be xda/xdr or exodusII file)."}>>> = ${input_name}
    timestep<<<{"description": "Index of the single timestep used or \"LATEST\" for the last timestep (exodusII only).  If not supplied, time interpolation will occur."}>>> = LATEST
  []
[]

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [tritium_s]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = ${tritium_s_scaling}
  []
  [tritium_trap]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = ${tritium_trap_scaling}
  []
  [tritium_2g]
  []
[]

[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [tritium_2g_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = tritium_2g
    function<<<{"description": "The initial condition function."}>>> = 'gaseous_concentration_1_function_ic'
  []
[]

[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
  [pore_position_func]
    type = SolutionFunction<<<{"description": "Function for reading a solution from file.", "href": "../../source/functions/SolutionFunction.html"}>>>
    solution<<<{"description": "The SolutionUserObject to extract data from."}>>> = initial_microstructure
    from_variable<<<{"description": "The name of the variable in the file that is to be extracted"}>>> = gr0
  []
  [gaseous_concentration_1_function_ic]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'len_pore  position_pore_int'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = '25   4521'
    expression<<<{"description": "The user defined function."}>>> = '${tritium_2g_concentration}*0.5*(1.0+tanh((sqrt(x*x+y*y)-position_pore_int)/(len_pore/2.2)))'
  []
[]

[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [tritium_2g_sides_d] # Fix concentration on the sides
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'right top'
    value<<<{"description": "Value of the BC"}>>> = ${tritium_2g_concentration}
  []
[]

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [pore]
    initial_from_file_var<<<{"description": "Gives the name of a variable for which to read an initial condition from a mesh file"}>>> = gr0
    initial_from_file_timestep<<<{"description": "Gives the timestep (or \"LATEST\") for which to read a solution from a file for a given variable. (Default: LATEST)"}>>> = LATEST
  []
  # Used to prevent negative concentrations
  [bounds_dummy_tritium_s]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = LAGRANGE
  []
  [bounds_dummy_tritium_2g]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = LAGRANGE
  []
[]

[Bounds<<<{"href": "../../syntax/Bounds/index.html"}>>>]
  # To prevent negative concentrations
  [tritium_s_lower_bound]
    type = ConstantBounds<<<{"description": "Provides constant bound of a variable for the PETSc's variational inequalities solver", "href": "../../source/bounds/ConstantBounds.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = bounds_dummy_tritium_s
    bounded_variable<<<{"description": "The variable to be bounded"}>>> = tritium_s
    bound_type<<<{"description": "Type of bound. 'upper' refers to the upper bound. 'lower' refers to the lower value."}>>> = lower
    bound_value<<<{"description": "The value of bound for the variable"}>>> = 0.
  []
  [tritium_g_lower_bound]
    type = ConstantBounds<<<{"description": "Provides constant bound of a variable for the PETSc's variational inequalities solver", "href": "../../source/bounds/ConstantBounds.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = bounds_dummy_tritium_2g
    bounded_variable<<<{"description": "The variable to be bounded"}>>> = tritium_2g
    bound_type<<<{"description": "Type of bound. 'upper' refers to the upper bound. 'lower' refers to the lower value."}>>> = lower
    bound_value<<<{"description": "The value of bound for the variable"}>>> = 0.
  []
[]

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  # ========================================== Time derivative for each variable
  [dtritium_sdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_s'
  []
  [dtritium_2gdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_2g'
  []
  [dtritium_trapdt]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
  []
  # ============================================================= Bulk diffusion
  [Diff_tritium]
    type = ADMatDiffusion<<<{"description": "Diffusion equation kernel that takes an isotropic diffusivity from a material property", "href": "../../source/kernels/ADMatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = D_tritium
  []
  # ============================================================= Pore diffusion
  [Diff_tritium2g]
    type = ADMatDiffusion<<<{"description": "Diffusion equation kernel that takes an isotropic diffusivity from a material property", "href": "../../source/kernels/ADMatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = D_pore
  []
  # ================================================ Tritium trapping/detrapping
  [Trapping]
    type = CoupledTimeDerivative<<<{"description": "Time derivative Kernel that acts on a coupled variable. Weak form: $(\\psi_i, \\frac{\\partial v_h}{\\partial t})$.", "href": "../../source/kernels/CoupledTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    v<<<{"description": "Coupled variable"}>>> = 'tritium_trap'
  []
  [Detrapping_trap]
    type = ADMatReaction<<<{"description": "Kernel representing the contribution of the PDE term $-L*v$, where $L$ is a reaction rate material property, $v$ is a scalar variable (nonlinear or coupled), and whose Jacobian contribution is calculated using automatic differentiation.", "href": "../../source/kernels/ADMatReaction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
    reaction_rate<<<{"description": "The reaction_rate used with the kernel."}>>> = ${fparse scale_time * detrapping_rate_0} # from 1/s to non-dimensional
  []
  [Trapping_trap]
    type = ADMatCoupledDefectAnnihilation<<<{"description": "Kernel to add K*v*(u0-u), where K=annihilation rate, u=variable, u0=equilibrium, v=coupled variable", "href": "../../source/kernels/ADMatCoupledDefectAnnihilation.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 'tritium_trap'
    eq_concentration<<<{"description": "The equilibrium concentration of the variable used with the kernel"}>>> = ${tritium_trap_0_scaled}
    v<<<{"description": "Set this to make v a coupled variable, otherwise it will use the kernel's nonlinear variable for v"}>>> = tritium_s
    annihilation_rate<<<{"description": "The annihilation rate used with the kernel"}>>> = ${fparse scale_time*scale_quantity/scale_length^3*trapping_rate_0} # from um^3/s/mol to non-dimensional
  []
  # ========================= Surface Reactions ================================
  # ======================================= Diatomic molecule dissolve at surface
  # T2(g) -> 2 T(s) rate K [T2] (1-theta) ======================================
  [Surface_Reaction_dissolution_tritium_2g] # T2(g) -> 2 T(s) rate -K*tritium_2g
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_2g'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = -1
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_gXYs
  []
  [Surface_Reaction_dissolution_tritium_2g_tritium_s] # T2(g) -> 2 T(s) rate 2*K*tritium_2g
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_2g'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = 2
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_gXYs
  []
  # ======================================= Diatomic molecule combine at surface
  # 2 T(s) -> T2(g) rate K [T]^2 ===============================================
  [Surface_Reaction_release_tritium_s_tritium_2g] # 2 T(s) -> T2(g) rate -2*K*tritium_s^2
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_s
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_s tritium_s'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = -2
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_sXYg
  []
  [Surface_Reaction_formation_tritium_2g] # 2 T(s) -> T2(g) rate K*tritium_s^2
    type = ADMatReactionFlexible<<<{"description": "Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product", "href": "../../source/kernels/ADMatReactionFlexible.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tritium_2g
    vs<<<{"description": "Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v"}>>> = 'tritium_s tritium_s'
    coeff<<<{"description": "A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species."}>>> = 1
    reaction_rate_name<<<{"description": "The reaction rate used with the kernel"}>>> = ReactionRateSurface_sXYg
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [pore_aux]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = pore
    function<<<{"description": "The function to use as the value"}>>> = 'pore_position_func'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL'
  []
[]

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  #===================================================== Interpolation functions
  [hsurfaceAD] # equal to 1 at pore surface, 0 elsewhere. The values of hsurface should be smaller than the ones for hmat or hpore to allow the diffusion of the species out of the surface
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '16*(1-pore)^2*pore^2'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hsurfaceAD'
  []
  [hmatAD] # equal to 1 in material, and 0 in pore.
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '(1-pore)*(1-pore)'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hmatAD'
  []
  [hporeAD] # equal to 1 in pore, and 1 in material. # notice that hpore overlaps a bit with hmat to prevent some agglomeration of chemicals on the edges of the surface.
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'pore*pore'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'hporeAD'
  []
  #================================================= Local species concentration
  [Tritium_ceramic]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'tritium_ceramic'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'tritium_s tritium_trap'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'tritium_s+tritium_trap'
  []
  [Tritium_pore]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'tritium_pore'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'tritium_2g'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '2*tritium_2g'
  []
  [Tritium_total]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'tritium_tot'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'tritium_s tritium_trap tritium_2g'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'tritium_s+tritium_trap+2*tritium_2g'
  []
  #=============================================== Occupied sites on the surface
  [theta]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'theta'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'tritium_s tritium_trap'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '(tritium_s + tritium_trap)/(${available_sites_mult}*${tritium_trap_0_scaled})'
  []
  #====================================================== Diffusion coefficients
  [Diffusion_coefficient_D] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'D_tritium'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hmatAD(pore)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}/${scale_length}^2*(hmatAD*${Diffusion_coefficient_Db_D0} + (1-hmatAD)*${Diffusion_min})'
    derivative_order<<<{"description": "Maximum order of derivatives taken"}>>> = 2
  []
  [Diffusion_coefficient_D_pore] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'D_pore'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hporeAD(pore)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}/${scale_length}^2*(hporeAD*${Diffusion_coefficient_D_pore_D0} + (1-hporeAD)*${Diffusion_min})'
  []
  #============================================ Surface reaction rate of tritium
  [ReactionRateSurface_sXYg_ref] # from surface to gaseous for diatomic molecules # from um^3/s/mol to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = ReactionRateSurface_sXYg
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hsurfaceAD(pore)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}*${scale_quantity}/${scale_length}^3*${reactionRateSurface_sXYg_ref_0}*hsurfaceAD'
  []
  [ReactionRateSurface_gXYs_ref] # from gaseous to surface for diatomic molecules # from 1/s to non-dimensional
    type = ADDerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../source/materials/DerivativeParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = ReactionRateSurface_gXYs
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'pore tritium_s tritium_trap'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'hsurfaceAD(pore) theta(tritium_s,tritium_trap)'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${scale_time}*${reactionRateSurface_gXYs_ref_0}*hsurfaceAD*(1-theta)^2'
  []
[]

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [Tritium_amount_free]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_s
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_trapped]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_trap
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_ceramic]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_ceramic
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium2g_amount]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = tritium_2g
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_pore]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_pore
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_total]
    type = ADElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = tritium_tot
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'INITIAL TIMESTEP_END'
  []
  [Surface_tot]
    type = ElementIntegralMaterialProperty<<<{"description": "Compute the integral of the material property over the domain", "href": "../../source/postprocessors/ElementIntegralMaterialProperty.html"}>>>
    mat_prop<<<{"description": "The name of the material property"}>>> = 1
  []
  [dt]
    type = TimestepSize<<<{"description": "Reports the timestep size", "href": "../../source/postprocessors/TimestepSize.html"}>>>
  []
[]

[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
  [Newtonlu]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    solve_type<<<{"description": "PJFNK: Preconditioned Jacobian-Free Newton Krylov JFNK: Jacobian-Free Newton Krylov NEWTON: Full Newton Solve FD: Use finite differences to compute Jacobian LINEAR: Solving a linear problem"}>>> = 'NEWTON'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -sub_pc_type -snes_type'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'asm      lu           vinewtonrsls' # This petsc option helps prevent negative concentrations'
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  scheme = bdf2

  nl_rel_tol = 1e-10
  end_time = 2000
  dtmax = 50
  nl_max_its = 14
  [TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    optimal_iterations = 12
    iteration_window = 1
    growth_factor = 1.2
    dt = 2.37e-7
    cutback_factor = 0.75
    cutback_factor_at_failure = 0.75
  []
[]

[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  [exodus]
    type = Exodus<<<{"description": "Object for output data in the Exodus format", "href": "../../source/outputs/Exodus.html"}>>>
    time_step_interval<<<{"description": "The interval (number of time steps) at which output occurs. Unless explicitly set, the default value of this parameter is set to infinity if the wall_time_interval is explicitly set."}>>> = 5
  []
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = ${output_name}
[]

References

Pierre-Clément A. Simon, Paul W. Humrickhouse, and Alexander D. Lindsay. Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8. IEEE Transactions on Plasma Science, pages 1–7, 2022. URL: https://ieeexplore.ieee.org/document/9831077/, doi:10.1109/TPS.2022.3183525.

@article{Simon2022,
    author = "Simon, Pierre-Clément A. and Humrickhouse, Paul W. and Lindsay, Alexander D.",
    doi = "10.1109/TPS.2022.3183525",
    issn = "0093-3813",
    journal = "IEEE Transactions on Plasma Science",
    pages = "1-7",
    title = "{Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8}",
    url = "https://ieeexplore.ieee.org/document/9831077/",
    year = "2022"
}

(test/tests/pore_scale_transport/2D_microstructure_reader_smoothing_base.i)

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Info:
# - This simulation predicts phase evolution based on a grain growth phase field model.
# - The pore and the ceramics represent two different phases, and the initial conditions
#   are provided with a black and white picture, which is read by ImageFunction.
# - The simulation smoothens the interface between the two phases.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_microstructure_reader_smoothing_base.i pore_structure_open.params

# mesh information
num_nodes_x = ${fparse 2*120} # (-)
num_nodes_y = ${fparse 2*120} # (-)
domain_start_x = ${units 0 mum}
domain_start_y = ${units 0 mum}
domain_end_x = ${units 5425 mum}
domain_end_y = ${units 5425 mum}

# grain growth parameters (since the point of the simulation is only to get a smooth interface, the model parameters can be selected arbitrarily)
op_num = 2 # Number of grains (-)
mobility_prefactor = ${units 2.5e-6 mum^4/J/s}
GB_energy = ${units 0.7 J/mum^2}
activation_energy = ${units 0.23 eV}
temperature = ${units 700 K}
width_diffuse_GB = ${units 50 mum} # Width of the diffuse GB

# image function option
threshold_image_function = 255.5 # (-)

[Mesh]
  [gen]
    type = DistributedRectilinearMeshGenerator
    dim = 2
    nx = ${num_nodes_x}
    ny = ${num_nodes_y}
    xmin = ${domain_start_x}
    xmax = ${domain_end_x}
    ymin = ${domain_start_y}
    ymax = ${domain_end_y}
  []
[]

[GlobalParams]
  op_num = ${op_num}
  var_name_base = gr # Base name of grains
[]

[UserObjects]
  [grain_tracker]
    type = GrainTracker
    flood_entity_type = ELEMENTAL
    compute_halo_maps = true
  []
[]

[Functions]
  [image_func0]
    type = ImageFunction
    file = ${input_name}
    threshold = ${threshold_image_function}
    upper_value = 1
    lower_value = 0
  []
  [image_func1]
    type = ImageFunction
    file = ${input_name}
    threshold = ${threshold_image_function}
    upper_value = 0
    lower_value = 1
  []
[]

[ICs]
  [gr0_ic]
    type = FunctionIC
    function = image_func0
    variable = gr0
  []
  [gr1_ic]
    type = FunctionIC
    function = image_func1
    variable = gr1
  []
[]

[Variables]
  [PolycrystalVariables]
    # Custom action that creates all of the phase variables and sets their initial condition
  []
[]

[AuxVariables]
  [bnds]
    # Variable used to visualize the interfaces in the simulation
  []
  [unique_grains]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [PolycrystalKernel]
    # Custom action creating all necessary kernels for phase evolution.  All input parameters are up in GlobalParams
  []
[]

[AuxKernels]
  [bnds_aux]
    type = BndsCalcAux
    variable = bnds
    execute_on = 'initial timestep_end'
  []
  [unique_grains]
    type = FeatureFloodCountAux
    variable = unique_grains
    execute_on = timestep_end
    flood_counter = grain_tracker
    field_display = UNIQUE_REGION
  []
[]

[Materials]
  [InterfaceEvolution]
    type = GBEvolution
    GBmob0 = ${mobility_prefactor}
    GBenergy = ${GB_energy}
    Q = ${activation_energy}
    T = ${temperature}
    wGB = ${width_diffuse_GB}
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  l_max_its = 20
  nl_max_its = 9
  nl_rel_tol = 1e-10
  [TimeStepper]
    type = IterationAdaptiveDT
    cutback_factor = 0.9
    dt = 3e-2
    growth_factor = 1.1
    optimal_iterations = 7
  []
  dtmax = 25

  start_time = 0.0
  end_time = 100
[]

[Outputs]
  file_base = ${output_name}
  [exodus]
    type = Exodus
    execute_on = 'FINAL'
  []
[]

(test/tests/pore_scale_transport/2D_microstructure_reader_smoothing_base.i)

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Info:
# - This simulation predicts phase evolution based on a grain growth phase field model.
# - The pore and the ceramics represent two different phases, and the initial conditions
#   are provided with a black and white picture, which is read by ImageFunction.
# - The simulation smoothens the interface between the two phases.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_microstructure_reader_smoothing_base.i pore_structure_open.params

# mesh information
num_nodes_x = ${fparse 2*120} # (-)
num_nodes_y = ${fparse 2*120} # (-)
domain_start_x = ${units 0 mum}
domain_start_y = ${units 0 mum}
domain_end_x = ${units 5425 mum}
domain_end_y = ${units 5425 mum}

# grain growth parameters (since the point of the simulation is only to get a smooth interface, the model parameters can be selected arbitrarily)
op_num = 2 # Number of grains (-)
mobility_prefactor = ${units 2.5e-6 mum^4/J/s}
GB_energy = ${units 0.7 J/mum^2}
activation_energy = ${units 0.23 eV}
temperature = ${units 700 K}
width_diffuse_GB = ${units 50 mum} # Width of the diffuse GB

# image function option
threshold_image_function = 255.5 # (-)

[Mesh]
  [gen]
    type = DistributedRectilinearMeshGenerator
    dim = 2
    nx = ${num_nodes_x}
    ny = ${num_nodes_y}
    xmin = ${domain_start_x}
    xmax = ${domain_end_x}
    ymin = ${domain_start_y}
    ymax = ${domain_end_y}
  []
[]

[GlobalParams]
  op_num = ${op_num}
  var_name_base = gr # Base name of grains
[]

[UserObjects]
  [grain_tracker]
    type = GrainTracker
    flood_entity_type = ELEMENTAL
    compute_halo_maps = true
  []
[]

[Functions]
  [image_func0]
    type = ImageFunction
    file = ${input_name}
    threshold = ${threshold_image_function}
    upper_value = 1
    lower_value = 0
  []
  [image_func1]
    type = ImageFunction
    file = ${input_name}
    threshold = ${threshold_image_function}
    upper_value = 0
    lower_value = 1
  []
[]

[ICs]
  [gr0_ic]
    type = FunctionIC
    function = image_func0
    variable = gr0
  []
  [gr1_ic]
    type = FunctionIC
    function = image_func1
    variable = gr1
  []
[]

[Variables]
  [PolycrystalVariables]
    # Custom action that creates all of the phase variables and sets their initial condition
  []
[]

[AuxVariables]
  [bnds]
    # Variable used to visualize the interfaces in the simulation
  []
  [unique_grains]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [PolycrystalKernel]
    # Custom action creating all necessary kernels for phase evolution.  All input parameters are up in GlobalParams
  []
[]

[AuxKernels]
  [bnds_aux]
    type = BndsCalcAux
    variable = bnds
    execute_on = 'initial timestep_end'
  []
  [unique_grains]
    type = FeatureFloodCountAux
    variable = unique_grains
    execute_on = timestep_end
    flood_counter = grain_tracker
    field_display = UNIQUE_REGION
  []
[]

[Materials]
  [InterfaceEvolution]
    type = GBEvolution
    GBmob0 = ${mobility_prefactor}
    GBenergy = ${GB_energy}
    Q = ${activation_energy}
    T = ${temperature}
    wGB = ${width_diffuse_GB}
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  l_max_its = 20
  nl_max_its = 9
  nl_rel_tol = 1e-10
  [TimeStepper]
    type = IterationAdaptiveDT
    cutback_factor = 0.9
    dt = 3e-2
    growth_factor = 1.1
    optimal_iterations = 7
  []
  dtmax = 25

  start_time = 0.0
  end_time = 100
[]

[Outputs]
  file_base = ${output_name}
  [exodus]
    type = Exodus
    execute_on = 'FINAL'
  []
[]

(test/tests/pore_scale_transport/2D_absorption_base.i)

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Phase Field Model:   Isotropic bulk diffusion, pore surface reactions, isotropic pore diffusion
# Type:                Transient
# Phase structure:     Contains two phases: the bulk of a material, and the pores. An interface separates the two
# Physics:             tritium trapping and detrapping, surface reactions, diffusion
# Species variable:    tritium_s (Ts), tritium_trap (trapped T), tritium_2g (T2(g))
# Species AuxVariable: --
# Chemistry:           at the pore surface: (s) means in the solid, (g) means in gas form
#    Diatomic molecule dissolving at surface
#                                 T2(g) -> 2 T(s) rate K [T2] (1-theta)^2
#    Diatomic molecule combining at surface
#                                 2 T(s) -> T2(g) rate K [T]^2
#
# BCs:                 No flux for tritium and tritium_trap, zero at pore center for gaseous tritium in every form
#
#
# Info:
# - This input file uses the microstructure generated by 2D_microstructure_reader_smoothing_base.i
#   and performs tritium transport simulation on this microstructure.
# - Basic 2D non-dimensional input file for the transport of a solute in a microstructure with pores.
# - The goal of this input base and associated parameter files is to have the same simulation with different pore configurations
#   to show the effect of microstructure on tritium release.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_absorption_base.i pore_structure_open_absorption.params

# Scaling factors
scale_quantity = ${units 1e-18 moles}
scale_time = ${units 1 s}
scale_length = ${units 1 mum}
tritium_trap_scaling = 1e3 # (-)
tritium_s_scaling = 1e2 # (-)

# Conditions
tritium_2g_concentration = ${units 0.1121 mol/mum^3}

# Material properties
tritium_trap_0 = ${units 3.472805925e-18 mol/mum^3}
tritium_trap_0_scaled = ${fparse tritium_trap_0*scale_length^3/scale_quantity} # non-dimensional
available_sites_mult = 4.313540691 # (-) should be >1, otherwise there are fewer 'available sites' than 'trapping sites'
Diffusion_min = ${units 0 mum^2/s} # used as diffusion coefficients where I don't want the species to move (more stable than 0)
Diffusion_coefficient_Db_D0 = ${units 1649969.728 mum^2/s}
Diffusion_coefficient_D_pore_D0 = ${units 224767738.6 mum^2/s}
detrapping_rate_0 = ${units -4.065104516 1/s} # should be negative
trapping_rate_0 = ${units 8.333146645e+16 1/s/mol} # should be positive
reactionRateSurface_sXYg_ref_0 = ${units 1549103878000000.0 1/s/mol}
reactionRateSurface_gXYs_ref_0 = ${units 34.69205021 1/s}

[Mesh]
  file = ${input_name}
[]

[UserObjects]
  [initial_microstructure]
    type = SolutionUserObject
    mesh = ${input_name}
    timestep = LATEST
  []
[]

[Variables]
  [tritium_s]
    scaling = ${tritium_s_scaling}
  []
  [tritium_trap]
    scaling = ${tritium_trap_scaling}
  []
  [tritium_2g]
  []
[]

[ICs]
  [tritium_2g_ic]
    type = FunctionIC
    variable = tritium_2g
    function = 'gaseous_concentration_1_function_ic'
  []
[]

[Functions]
  [pore_position_func]
    type = SolutionFunction
    solution = initial_microstructure
    from_variable = gr0
  []
  [gaseous_concentration_1_function_ic]
    type = ParsedFunction
    symbol_names = 'len_pore  position_pore_int'
    symbol_values = '25   4521'
    expression = '${tritium_2g_concentration}*0.5*(1.0+tanh((sqrt(x*x+y*y)-position_pore_int)/(len_pore/2.2)))'
  []
[]

[BCs]
  [tritium_2g_sides_d] # Fix concentration on the sides
    type = DirichletBC
    variable = tritium_2g
    boundary = 'right top'
    value = ${tritium_2g_concentration}
  []
[]

[AuxVariables]
  [pore]
    initial_from_file_var = gr0
    initial_from_file_timestep = LATEST
  []
  # Used to prevent negative concentrations
  [bounds_dummy_tritium_s]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_tritium_2g]
    order = FIRST
    family = LAGRANGE
  []
[]

[Bounds]
  # To prevent negative concentrations
  [tritium_s_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_tritium_s
    bounded_variable = tritium_s
    bound_type = lower
    bound_value = 0.
  []
  [tritium_g_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_tritium_2g
    bounded_variable = tritium_2g
    bound_type = lower
    bound_value = 0.
  []
[]

[Kernels]
  # ========================================== Time derivative for each variable
  [dtritium_sdt]
    type = TimeDerivative
    variable = 'tritium_s'
  []
  [dtritium_2gdt]
    type = TimeDerivative
    variable = 'tritium_2g'
  []
  [dtritium_trapdt]
    type = TimeDerivative
    variable = 'tritium_trap'
  []
  # ============================================================= Bulk diffusion
  [Diff_tritium]
    type = ADMatDiffusion
    variable = tritium_s
    diffusivity = D_tritium
  []
  # ============================================================= Pore diffusion
  [Diff_tritium2g]
    type = ADMatDiffusion
    variable = tritium_2g
    diffusivity = D_pore
  []
  # ================================================ Tritium trapping/detrapping
  [Trapping]
    type = CoupledTimeDerivative
    variable = tritium_s
    v = 'tritium_trap'
  []
  [Detrapping_trap]
    type = ADMatReaction
    variable = 'tritium_trap'
    reaction_rate = ${fparse scale_time * detrapping_rate_0} # from 1/s to non-dimensional
  []
  [Trapping_trap]
    type = ADMatCoupledDefectAnnihilation
    variable = 'tritium_trap'
    eq_concentration = ${tritium_trap_0_scaled}
    v = tritium_s
    annihilation_rate = ${fparse scale_time*scale_quantity/scale_length^3*trapping_rate_0} # from um^3/s/mol to non-dimensional
  []
  # ========================= Surface Reactions ================================
  # ======================================= Diatomic molecule dissolve at surface
  # T2(g) -> 2 T(s) rate K [T2] (1-theta) ======================================
  [Surface_Reaction_dissolution_tritium_2g] # T2(g) -> 2 T(s) rate -K*tritium_2g
    type = ADMatReactionFlexible
    variable = tritium_2g
    vs = 'tritium_2g'
    coeff = -1
    reaction_rate_name = ReactionRateSurface_gXYs
  []
  [Surface_Reaction_dissolution_tritium_2g_tritium_s] # T2(g) -> 2 T(s) rate 2*K*tritium_2g
    type = ADMatReactionFlexible
    variable = tritium_s
    vs = 'tritium_2g'
    coeff = 2
    reaction_rate_name = ReactionRateSurface_gXYs
  []
  # ======================================= Diatomic molecule combine at surface
  # 2 T(s) -> T2(g) rate K [T]^2 ===============================================
  [Surface_Reaction_release_tritium_s_tritium_2g] # 2 T(s) -> T2(g) rate -2*K*tritium_s^2
    type = ADMatReactionFlexible
    variable = tritium_s
    vs = 'tritium_s tritium_s'
    coeff = -2
    reaction_rate_name = ReactionRateSurface_sXYg
  []
  [Surface_Reaction_formation_tritium_2g] # 2 T(s) -> T2(g) rate K*tritium_s^2
    type = ADMatReactionFlexible
    variable = tritium_2g
    vs = 'tritium_s tritium_s'
    coeff = 1
    reaction_rate_name = ReactionRateSurface_sXYg
  []
[]

[AuxKernels]
  [pore_aux]
    type = FunctionAux
    variable = pore
    function = 'pore_position_func'
    execute_on = 'INITIAL'
  []
[]

[Materials]
  #===================================================== Interpolation functions
  [hsurfaceAD] # equal to 1 at pore surface, 0 elsewhere. The values of hsurface should be smaller than the ones for hmat or hpore to allow the diffusion of the species out of the surface
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = '16*(1-pore)^2*pore^2'
    property_name = 'hsurfaceAD'
  []
  [hmatAD] # equal to 1 in material, and 0 in pore.
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = '(1-pore)*(1-pore)'
    property_name = 'hmatAD'
  []
  [hporeAD] # equal to 1 in pore, and 1 in material. # notice that hpore overlaps a bit with hmat to prevent some agglomeration of chemicals on the edges of the surface.
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = 'pore*pore'
    property_name = 'hporeAD'
  []
  #================================================= Local species concentration
  [Tritium_ceramic]
    type = ADParsedMaterial
    property_name = 'tritium_ceramic'
    coupled_variables = 'tritium_s tritium_trap'
    expression = 'tritium_s+tritium_trap'
  []
  [Tritium_pore]
    type = ADParsedMaterial
    property_name = 'tritium_pore'
    coupled_variables = 'tritium_2g'
    expression = '2*tritium_2g'
  []
  [Tritium_total]
    type = ADParsedMaterial
    property_name = 'tritium_tot'
    coupled_variables = 'tritium_s tritium_trap tritium_2g'
    expression = 'tritium_s+tritium_trap+2*tritium_2g'
  []
  #=============================================== Occupied sites on the surface
  [theta]
    type = ADParsedMaterial
    property_name = 'theta'
    coupled_variables = 'tritium_s tritium_trap'
    expression = '(tritium_s + tritium_trap)/(${available_sites_mult}*${tritium_trap_0_scaled})'
  []
  #====================================================== Diffusion coefficients
  [Diffusion_coefficient_D] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = 'D_tritium'
    coupled_variables = 'pore'
    material_property_names = 'hmatAD(pore)'
    expression = '${scale_time}/${scale_length}^2*(hmatAD*${Diffusion_coefficient_Db_D0} + (1-hmatAD)*${Diffusion_min})'
    derivative_order = 2
  []
  [Diffusion_coefficient_D_pore] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = 'D_pore'
    material_property_names = 'hporeAD(pore)'
    expression = '${scale_time}/${scale_length}^2*(hporeAD*${Diffusion_coefficient_D_pore_D0} + (1-hporeAD)*${Diffusion_min})'
  []
  #============================================ Surface reaction rate of tritium
  [ReactionRateSurface_sXYg_ref] # from surface to gaseous for diatomic molecules # from um^3/s/mol to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = ReactionRateSurface_sXYg
    material_property_names = 'hsurfaceAD(pore)'
    expression = '${scale_time}*${scale_quantity}/${scale_length}^3*${reactionRateSurface_sXYg_ref_0}*hsurfaceAD'
  []
  [ReactionRateSurface_gXYs_ref] # from gaseous to surface for diatomic molecules # from 1/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = ReactionRateSurface_gXYs
    coupled_variables = 'pore tritium_s tritium_trap'
    material_property_names = 'hsurfaceAD(pore) theta(tritium_s,tritium_trap)'
    expression = '${scale_time}*${reactionRateSurface_gXYs_ref_0}*hsurfaceAD*(1-theta)^2'
  []
[]

[Postprocessors]
  [Tritium_amount_free]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_s
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_trapped]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_trap
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_ceramic]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_ceramic
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium2g_amount]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_2g
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_pore]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_pore
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_total]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_tot
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Surface_tot]
    type = ElementIntegralMaterialProperty
    mat_prop = 1
  []
  [dt]
    type = TimestepSize
  []
[]

[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]
  type = Transient
  scheme = bdf2

  nl_rel_tol = 1e-10
  end_time = 2000
  dtmax = 50
  nl_max_its = 14
  [TimeStepper]
    type = IterationAdaptiveDT
    optimal_iterations = 12
    iteration_window = 1
    growth_factor = 1.2
    dt = 2.37e-7
    cutback_factor = 0.75
    cutback_factor_at_failure = 0.75
  []
[]

[Outputs]
  [exodus]
    type = Exodus
    time_step_interval = 5
  []
  csv = true
  file_base = ${output_name}
[]

(test/tests/pore_scale_transport/pore_structure_closed_absorption.params)

input_name=pore_structure_closed/pore_structure_closed.e
output_name=pore_structure_closed/pore_structure_closed_absorption

(test/tests/pore_scale_transport/2D_absorption_base.i)

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Phase Field Model:   Isotropic bulk diffusion, pore surface reactions, isotropic pore diffusion
# Type:                Transient
# Phase structure:     Contains two phases: the bulk of a material, and the pores. An interface separates the two
# Physics:             tritium trapping and detrapping, surface reactions, diffusion
# Species variable:    tritium_s (Ts), tritium_trap (trapped T), tritium_2g (T2(g))
# Species AuxVariable: --
# Chemistry:           at the pore surface: (s) means in the solid, (g) means in gas form
#    Diatomic molecule dissolving at surface
#                                 T2(g) -> 2 T(s) rate K [T2] (1-theta)^2
#    Diatomic molecule combining at surface
#                                 2 T(s) -> T2(g) rate K [T]^2
#
# BCs:                 No flux for tritium and tritium_trap, zero at pore center for gaseous tritium in every form
#
#
# Info:
# - This input file uses the microstructure generated by 2D_microstructure_reader_smoothing_base.i
#   and performs tritium transport simulation on this microstructure.
# - Basic 2D non-dimensional input file for the transport of a solute in a microstructure with pores.
# - The goal of this input base and associated parameter files is to have the same simulation with different pore configurations
#   to show the effect of microstructure on tritium release.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_absorption_base.i pore_structure_open_absorption.params

# Scaling factors
scale_quantity = ${units 1e-18 moles}
scale_time = ${units 1 s}
scale_length = ${units 1 mum}
tritium_trap_scaling = 1e3 # (-)
tritium_s_scaling = 1e2 # (-)

# Conditions
tritium_2g_concentration = ${units 0.1121 mol/mum^3}

# Material properties
tritium_trap_0 = ${units 3.472805925e-18 mol/mum^3}
tritium_trap_0_scaled = ${fparse tritium_trap_0*scale_length^3/scale_quantity} # non-dimensional
available_sites_mult = 4.313540691 # (-) should be >1, otherwise there are fewer 'available sites' than 'trapping sites'
Diffusion_min = ${units 0 mum^2/s} # used as diffusion coefficients where I don't want the species to move (more stable than 0)
Diffusion_coefficient_Db_D0 = ${units 1649969.728 mum^2/s}
Diffusion_coefficient_D_pore_D0 = ${units 224767738.6 mum^2/s}
detrapping_rate_0 = ${units -4.065104516 1/s} # should be negative
trapping_rate_0 = ${units 8.333146645e+16 1/s/mol} # should be positive
reactionRateSurface_sXYg_ref_0 = ${units 1549103878000000.0 1/s/mol}
reactionRateSurface_gXYs_ref_0 = ${units 34.69205021 1/s}

[Mesh]
  file = ${input_name}
[]

[UserObjects]
  [initial_microstructure]
    type = SolutionUserObject
    mesh = ${input_name}
    timestep = LATEST
  []
[]

[Variables]
  [tritium_s]
    scaling = ${tritium_s_scaling}
  []
  [tritium_trap]
    scaling = ${tritium_trap_scaling}
  []
  [tritium_2g]
  []
[]

[ICs]
  [tritium_2g_ic]
    type = FunctionIC
    variable = tritium_2g
    function = 'gaseous_concentration_1_function_ic'
  []
[]

[Functions]
  [pore_position_func]
    type = SolutionFunction
    solution = initial_microstructure
    from_variable = gr0
  []
  [gaseous_concentration_1_function_ic]
    type = ParsedFunction
    symbol_names = 'len_pore  position_pore_int'
    symbol_values = '25   4521'
    expression = '${tritium_2g_concentration}*0.5*(1.0+tanh((sqrt(x*x+y*y)-position_pore_int)/(len_pore/2.2)))'
  []
[]

[BCs]
  [tritium_2g_sides_d] # Fix concentration on the sides
    type = DirichletBC
    variable = tritium_2g
    boundary = 'right top'
    value = ${tritium_2g_concentration}
  []
[]

[AuxVariables]
  [pore]
    initial_from_file_var = gr0
    initial_from_file_timestep = LATEST
  []
  # Used to prevent negative concentrations
  [bounds_dummy_tritium_s]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_tritium_2g]
    order = FIRST
    family = LAGRANGE
  []
[]

[Bounds]
  # To prevent negative concentrations
  [tritium_s_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_tritium_s
    bounded_variable = tritium_s
    bound_type = lower
    bound_value = 0.
  []
  [tritium_g_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_tritium_2g
    bounded_variable = tritium_2g
    bound_type = lower
    bound_value = 0.
  []
[]

[Kernels]
  # ========================================== Time derivative for each variable
  [dtritium_sdt]
    type = TimeDerivative
    variable = 'tritium_s'
  []
  [dtritium_2gdt]
    type = TimeDerivative
    variable = 'tritium_2g'
  []
  [dtritium_trapdt]
    type = TimeDerivative
    variable = 'tritium_trap'
  []
  # ============================================================= Bulk diffusion
  [Diff_tritium]
    type = ADMatDiffusion
    variable = tritium_s
    diffusivity = D_tritium
  []
  # ============================================================= Pore diffusion
  [Diff_tritium2g]
    type = ADMatDiffusion
    variable = tritium_2g
    diffusivity = D_pore
  []
  # ================================================ Tritium trapping/detrapping
  [Trapping]
    type = CoupledTimeDerivative
    variable = tritium_s
    v = 'tritium_trap'
  []
  [Detrapping_trap]
    type = ADMatReaction
    variable = 'tritium_trap'
    reaction_rate = ${fparse scale_time * detrapping_rate_0} # from 1/s to non-dimensional
  []
  [Trapping_trap]
    type = ADMatCoupledDefectAnnihilation
    variable = 'tritium_trap'
    eq_concentration = ${tritium_trap_0_scaled}
    v = tritium_s
    annihilation_rate = ${fparse scale_time*scale_quantity/scale_length^3*trapping_rate_0} # from um^3/s/mol to non-dimensional
  []
  # ========================= Surface Reactions ================================
  # ======================================= Diatomic molecule dissolve at surface
  # T2(g) -> 2 T(s) rate K [T2] (1-theta) ======================================
  [Surface_Reaction_dissolution_tritium_2g] # T2(g) -> 2 T(s) rate -K*tritium_2g
    type = ADMatReactionFlexible
    variable = tritium_2g
    vs = 'tritium_2g'
    coeff = -1
    reaction_rate_name = ReactionRateSurface_gXYs
  []
  [Surface_Reaction_dissolution_tritium_2g_tritium_s] # T2(g) -> 2 T(s) rate 2*K*tritium_2g
    type = ADMatReactionFlexible
    variable = tritium_s
    vs = 'tritium_2g'
    coeff = 2
    reaction_rate_name = ReactionRateSurface_gXYs
  []
  # ======================================= Diatomic molecule combine at surface
  # 2 T(s) -> T2(g) rate K [T]^2 ===============================================
  [Surface_Reaction_release_tritium_s_tritium_2g] # 2 T(s) -> T2(g) rate -2*K*tritium_s^2
    type = ADMatReactionFlexible
    variable = tritium_s
    vs = 'tritium_s tritium_s'
    coeff = -2
    reaction_rate_name = ReactionRateSurface_sXYg
  []
  [Surface_Reaction_formation_tritium_2g] # 2 T(s) -> T2(g) rate K*tritium_s^2
    type = ADMatReactionFlexible
    variable = tritium_2g
    vs = 'tritium_s tritium_s'
    coeff = 1
    reaction_rate_name = ReactionRateSurface_sXYg
  []
[]

[AuxKernels]
  [pore_aux]
    type = FunctionAux
    variable = pore
    function = 'pore_position_func'
    execute_on = 'INITIAL'
  []
[]

[Materials]
  #===================================================== Interpolation functions
  [hsurfaceAD] # equal to 1 at pore surface, 0 elsewhere. The values of hsurface should be smaller than the ones for hmat or hpore to allow the diffusion of the species out of the surface
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = '16*(1-pore)^2*pore^2'
    property_name = 'hsurfaceAD'
  []
  [hmatAD] # equal to 1 in material, and 0 in pore.
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = '(1-pore)*(1-pore)'
    property_name = 'hmatAD'
  []
  [hporeAD] # equal to 1 in pore, and 1 in material. # notice that hpore overlaps a bit with hmat to prevent some agglomeration of chemicals on the edges of the surface.
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = 'pore*pore'
    property_name = 'hporeAD'
  []
  #================================================= Local species concentration
  [Tritium_ceramic]
    type = ADParsedMaterial
    property_name = 'tritium_ceramic'
    coupled_variables = 'tritium_s tritium_trap'
    expression = 'tritium_s+tritium_trap'
  []
  [Tritium_pore]
    type = ADParsedMaterial
    property_name = 'tritium_pore'
    coupled_variables = 'tritium_2g'
    expression = '2*tritium_2g'
  []
  [Tritium_total]
    type = ADParsedMaterial
    property_name = 'tritium_tot'
    coupled_variables = 'tritium_s tritium_trap tritium_2g'
    expression = 'tritium_s+tritium_trap+2*tritium_2g'
  []
  #=============================================== Occupied sites on the surface
  [theta]
    type = ADParsedMaterial
    property_name = 'theta'
    coupled_variables = 'tritium_s tritium_trap'
    expression = '(tritium_s + tritium_trap)/(${available_sites_mult}*${tritium_trap_0_scaled})'
  []
  #====================================================== Diffusion coefficients
  [Diffusion_coefficient_D] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = 'D_tritium'
    coupled_variables = 'pore'
    material_property_names = 'hmatAD(pore)'
    expression = '${scale_time}/${scale_length}^2*(hmatAD*${Diffusion_coefficient_Db_D0} + (1-hmatAD)*${Diffusion_min})'
    derivative_order = 2
  []
  [Diffusion_coefficient_D_pore] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = 'D_pore'
    material_property_names = 'hporeAD(pore)'
    expression = '${scale_time}/${scale_length}^2*(hporeAD*${Diffusion_coefficient_D_pore_D0} + (1-hporeAD)*${Diffusion_min})'
  []
  #============================================ Surface reaction rate of tritium
  [ReactionRateSurface_sXYg_ref] # from surface to gaseous for diatomic molecules # from um^3/s/mol to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = ReactionRateSurface_sXYg
    material_property_names = 'hsurfaceAD(pore)'
    expression = '${scale_time}*${scale_quantity}/${scale_length}^3*${reactionRateSurface_sXYg_ref_0}*hsurfaceAD'
  []
  [ReactionRateSurface_gXYs_ref] # from gaseous to surface for diatomic molecules # from 1/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = ReactionRateSurface_gXYs
    coupled_variables = 'pore tritium_s tritium_trap'
    material_property_names = 'hsurfaceAD(pore) theta(tritium_s,tritium_trap)'
    expression = '${scale_time}*${reactionRateSurface_gXYs_ref_0}*hsurfaceAD*(1-theta)^2'
  []
[]

[Postprocessors]
  [Tritium_amount_free]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_s
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_trapped]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_trap
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_ceramic]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_ceramic
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium2g_amount]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_2g
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_pore]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_pore
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_total]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_tot
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Surface_tot]
    type = ElementIntegralMaterialProperty
    mat_prop = 1
  []
  [dt]
    type = TimestepSize
  []
[]

[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]
  type = Transient
  scheme = bdf2

  nl_rel_tol = 1e-10
  end_time = 2000
  dtmax = 50
  nl_max_its = 14
  [TimeStepper]
    type = IterationAdaptiveDT
    optimal_iterations = 12
    iteration_window = 1
    growth_factor = 1.2
    dt = 2.37e-7
    cutback_factor = 0.75
    cutback_factor_at_failure = 0.75
  []
[]

[Outputs]
  [exodus]
    type = Exodus
    time_step_interval = 5
  []
  csv = true
  file_base = ${output_name}
[]

(test/tests/pore_scale_transport/2D_microstructure_reader_smoothing_base.i)

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Info:
# - This simulation predicts phase evolution based on a grain growth phase field model.
# - The pore and the ceramics represent two different phases, and the initial conditions
#   are provided with a black and white picture, which is read by ImageFunction.
# - The simulation smoothens the interface between the two phases.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_microstructure_reader_smoothing_base.i pore_structure_open.params

# mesh information
num_nodes_x = ${fparse 2*120} # (-)
num_nodes_y = ${fparse 2*120} # (-)
domain_start_x = ${units 0 mum}
domain_start_y = ${units 0 mum}
domain_end_x = ${units 5425 mum}
domain_end_y = ${units 5425 mum}

# grain growth parameters (since the point of the simulation is only to get a smooth interface, the model parameters can be selected arbitrarily)
op_num = 2 # Number of grains (-)
mobility_prefactor = ${units 2.5e-6 mum^4/J/s}
GB_energy = ${units 0.7 J/mum^2}
activation_energy = ${units 0.23 eV}
temperature = ${units 700 K}
width_diffuse_GB = ${units 50 mum} # Width of the diffuse GB

# image function option
threshold_image_function = 255.5 # (-)

[Mesh]
  [gen]
    type = DistributedRectilinearMeshGenerator
    dim = 2
    nx = ${num_nodes_x}
    ny = ${num_nodes_y}
    xmin = ${domain_start_x}
    xmax = ${domain_end_x}
    ymin = ${domain_start_y}
    ymax = ${domain_end_y}
  []
[]

[GlobalParams]
  op_num = ${op_num}
  var_name_base = gr # Base name of grains
[]

[UserObjects]
  [grain_tracker]
    type = GrainTracker
    flood_entity_type = ELEMENTAL
    compute_halo_maps = true
  []
[]

[Functions]
  [image_func0]
    type = ImageFunction
    file = ${input_name}
    threshold = ${threshold_image_function}
    upper_value = 1
    lower_value = 0
  []
  [image_func1]
    type = ImageFunction
    file = ${input_name}
    threshold = ${threshold_image_function}
    upper_value = 0
    lower_value = 1
  []
[]

[ICs]
  [gr0_ic]
    type = FunctionIC
    function = image_func0
    variable = gr0
  []
  [gr1_ic]
    type = FunctionIC
    function = image_func1
    variable = gr1
  []
[]

[Variables]
  [PolycrystalVariables]
    # Custom action that creates all of the phase variables and sets their initial condition
  []
[]

[AuxVariables]
  [bnds]
    # Variable used to visualize the interfaces in the simulation
  []
  [unique_grains]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [PolycrystalKernel]
    # Custom action creating all necessary kernels for phase evolution.  All input parameters are up in GlobalParams
  []
[]

[AuxKernels]
  [bnds_aux]
    type = BndsCalcAux
    variable = bnds
    execute_on = 'initial timestep_end'
  []
  [unique_grains]
    type = FeatureFloodCountAux
    variable = unique_grains
    execute_on = timestep_end
    flood_counter = grain_tracker
    field_display = UNIQUE_REGION
  []
[]

[Materials]
  [InterfaceEvolution]
    type = GBEvolution
    GBmob0 = ${mobility_prefactor}
    GBenergy = ${GB_energy}
    Q = ${activation_energy}
    T = ${temperature}
    wGB = ${width_diffuse_GB}
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  l_max_its = 20
  nl_max_its = 9
  nl_rel_tol = 1e-10
  [TimeStepper]
    type = IterationAdaptiveDT
    cutback_factor = 0.9
    dt = 3e-2
    growth_factor = 1.1
    optimal_iterations = 7
  []
  dtmax = 25

  start_time = 0.0
  end_time = 100
[]

[Outputs]
  file_base = ${output_name}
  [exodus]
    type = Exodus
    execute_on = 'FINAL'
  []
[]

(test/tests/pore_scale_transport/2D_absorption_base.i)

# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Phase Field Model:   Isotropic bulk diffusion, pore surface reactions, isotropic pore diffusion
# Type:                Transient
# Phase structure:     Contains two phases: the bulk of a material, and the pores. An interface separates the two
# Physics:             tritium trapping and detrapping, surface reactions, diffusion
# Species variable:    tritium_s (Ts), tritium_trap (trapped T), tritium_2g (T2(g))
# Species AuxVariable: --
# Chemistry:           at the pore surface: (s) means in the solid, (g) means in gas form
#    Diatomic molecule dissolving at surface
#                                 T2(g) -> 2 T(s) rate K [T2] (1-theta)^2
#    Diatomic molecule combining at surface
#                                 2 T(s) -> T2(g) rate K [T]^2
#
# BCs:                 No flux for tritium and tritium_trap, zero at pore center for gaseous tritium in every form
#
#
# Info:
# - This input file uses the microstructure generated by 2D_microstructure_reader_smoothing_base.i
#   and performs tritium transport simulation on this microstructure.
# - Basic 2D non-dimensional input file for the transport of a solute in a microstructure with pores.
# - The goal of this input base and associated parameter files is to have the same simulation with different pore configurations
#   to show the effect of microstructure on tritium release.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_absorption_base.i pore_structure_open_absorption.params

# Scaling factors
scale_quantity = ${units 1e-18 moles}
scale_time = ${units 1 s}
scale_length = ${units 1 mum}
tritium_trap_scaling = 1e3 # (-)
tritium_s_scaling = 1e2 # (-)

# Conditions
tritium_2g_concentration = ${units 0.1121 mol/mum^3}

# Material properties
tritium_trap_0 = ${units 3.472805925e-18 mol/mum^3}
tritium_trap_0_scaled = ${fparse tritium_trap_0*scale_length^3/scale_quantity} # non-dimensional
available_sites_mult = 4.313540691 # (-) should be >1, otherwise there are fewer 'available sites' than 'trapping sites'
Diffusion_min = ${units 0 mum^2/s} # used as diffusion coefficients where I don't want the species to move (more stable than 0)
Diffusion_coefficient_Db_D0 = ${units 1649969.728 mum^2/s}
Diffusion_coefficient_D_pore_D0 = ${units 224767738.6 mum^2/s}
detrapping_rate_0 = ${units -4.065104516 1/s} # should be negative
trapping_rate_0 = ${units 8.333146645e+16 1/s/mol} # should be positive
reactionRateSurface_sXYg_ref_0 = ${units 1549103878000000.0 1/s/mol}
reactionRateSurface_gXYs_ref_0 = ${units 34.69205021 1/s}

[Mesh]
  file = ${input_name}
[]

[UserObjects]
  [initial_microstructure]
    type = SolutionUserObject
    mesh = ${input_name}
    timestep = LATEST
  []
[]

[Variables]
  [tritium_s]
    scaling = ${tritium_s_scaling}
  []
  [tritium_trap]
    scaling = ${tritium_trap_scaling}
  []
  [tritium_2g]
  []
[]

[ICs]
  [tritium_2g_ic]
    type = FunctionIC
    variable = tritium_2g
    function = 'gaseous_concentration_1_function_ic'
  []
[]

[Functions]
  [pore_position_func]
    type = SolutionFunction
    solution = initial_microstructure
    from_variable = gr0
  []
  [gaseous_concentration_1_function_ic]
    type = ParsedFunction
    symbol_names = 'len_pore  position_pore_int'
    symbol_values = '25   4521'
    expression = '${tritium_2g_concentration}*0.5*(1.0+tanh((sqrt(x*x+y*y)-position_pore_int)/(len_pore/2.2)))'
  []
[]

[BCs]
  [tritium_2g_sides_d] # Fix concentration on the sides
    type = DirichletBC
    variable = tritium_2g
    boundary = 'right top'
    value = ${tritium_2g_concentration}
  []
[]

[AuxVariables]
  [pore]
    initial_from_file_var = gr0
    initial_from_file_timestep = LATEST
  []
  # Used to prevent negative concentrations
  [bounds_dummy_tritium_s]
    order = FIRST
    family = LAGRANGE
  []
  [bounds_dummy_tritium_2g]
    order = FIRST
    family = LAGRANGE
  []
[]

[Bounds]
  # To prevent negative concentrations
  [tritium_s_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_tritium_s
    bounded_variable = tritium_s
    bound_type = lower
    bound_value = 0.
  []
  [tritium_g_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy_tritium_2g
    bounded_variable = tritium_2g
    bound_type = lower
    bound_value = 0.
  []
[]

[Kernels]
  # ========================================== Time derivative for each variable
  [dtritium_sdt]
    type = TimeDerivative
    variable = 'tritium_s'
  []
  [dtritium_2gdt]
    type = TimeDerivative
    variable = 'tritium_2g'
  []
  [dtritium_trapdt]
    type = TimeDerivative
    variable = 'tritium_trap'
  []
  # ============================================================= Bulk diffusion
  [Diff_tritium]
    type = ADMatDiffusion
    variable = tritium_s
    diffusivity = D_tritium
  []
  # ============================================================= Pore diffusion
  [Diff_tritium2g]
    type = ADMatDiffusion
    variable = tritium_2g
    diffusivity = D_pore
  []
  # ================================================ Tritium trapping/detrapping
  [Trapping]
    type = CoupledTimeDerivative
    variable = tritium_s
    v = 'tritium_trap'
  []
  [Detrapping_trap]
    type = ADMatReaction
    variable = 'tritium_trap'
    reaction_rate = ${fparse scale_time * detrapping_rate_0} # from 1/s to non-dimensional
  []
  [Trapping_trap]
    type = ADMatCoupledDefectAnnihilation
    variable = 'tritium_trap'
    eq_concentration = ${tritium_trap_0_scaled}
    v = tritium_s
    annihilation_rate = ${fparse scale_time*scale_quantity/scale_length^3*trapping_rate_0} # from um^3/s/mol to non-dimensional
  []
  # ========================= Surface Reactions ================================
  # ======================================= Diatomic molecule dissolve at surface
  # T2(g) -> 2 T(s) rate K [T2] (1-theta) ======================================
  [Surface_Reaction_dissolution_tritium_2g] # T2(g) -> 2 T(s) rate -K*tritium_2g
    type = ADMatReactionFlexible
    variable = tritium_2g
    vs = 'tritium_2g'
    coeff = -1
    reaction_rate_name = ReactionRateSurface_gXYs
  []
  [Surface_Reaction_dissolution_tritium_2g_tritium_s] # T2(g) -> 2 T(s) rate 2*K*tritium_2g
    type = ADMatReactionFlexible
    variable = tritium_s
    vs = 'tritium_2g'
    coeff = 2
    reaction_rate_name = ReactionRateSurface_gXYs
  []
  # ======================================= Diatomic molecule combine at surface
  # 2 T(s) -> T2(g) rate K [T]^2 ===============================================
  [Surface_Reaction_release_tritium_s_tritium_2g] # 2 T(s) -> T2(g) rate -2*K*tritium_s^2
    type = ADMatReactionFlexible
    variable = tritium_s
    vs = 'tritium_s tritium_s'
    coeff = -2
    reaction_rate_name = ReactionRateSurface_sXYg
  []
  [Surface_Reaction_formation_tritium_2g] # 2 T(s) -> T2(g) rate K*tritium_s^2
    type = ADMatReactionFlexible
    variable = tritium_2g
    vs = 'tritium_s tritium_s'
    coeff = 1
    reaction_rate_name = ReactionRateSurface_sXYg
  []
[]

[AuxKernels]
  [pore_aux]
    type = FunctionAux
    variable = pore
    function = 'pore_position_func'
    execute_on = 'INITIAL'
  []
[]

[Materials]
  #===================================================== Interpolation functions
  [hsurfaceAD] # equal to 1 at pore surface, 0 elsewhere. The values of hsurface should be smaller than the ones for hmat or hpore to allow the diffusion of the species out of the surface
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = '16*(1-pore)^2*pore^2'
    property_name = 'hsurfaceAD'
  []
  [hmatAD] # equal to 1 in material, and 0 in pore.
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = '(1-pore)*(1-pore)'
    property_name = 'hmatAD'
  []
  [hporeAD] # equal to 1 in pore, and 1 in material. # notice that hpore overlaps a bit with hmat to prevent some agglomeration of chemicals on the edges of the surface.
    type = ADDerivativeParsedMaterial
    coupled_variables = 'pore'
    expression = 'pore*pore'
    property_name = 'hporeAD'
  []
  #================================================= Local species concentration
  [Tritium_ceramic]
    type = ADParsedMaterial
    property_name = 'tritium_ceramic'
    coupled_variables = 'tritium_s tritium_trap'
    expression = 'tritium_s+tritium_trap'
  []
  [Tritium_pore]
    type = ADParsedMaterial
    property_name = 'tritium_pore'
    coupled_variables = 'tritium_2g'
    expression = '2*tritium_2g'
  []
  [Tritium_total]
    type = ADParsedMaterial
    property_name = 'tritium_tot'
    coupled_variables = 'tritium_s tritium_trap tritium_2g'
    expression = 'tritium_s+tritium_trap+2*tritium_2g'
  []
  #=============================================== Occupied sites on the surface
  [theta]
    type = ADParsedMaterial
    property_name = 'theta'
    coupled_variables = 'tritium_s tritium_trap'
    expression = '(tritium_s + tritium_trap)/(${available_sites_mult}*${tritium_trap_0_scaled})'
  []
  #====================================================== Diffusion coefficients
  [Diffusion_coefficient_D] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = 'D_tritium'
    coupled_variables = 'pore'
    material_property_names = 'hmatAD(pore)'
    expression = '${scale_time}/${scale_length}^2*(hmatAD*${Diffusion_coefficient_Db_D0} + (1-hmatAD)*${Diffusion_min})'
    derivative_order = 2
  []
  [Diffusion_coefficient_D_pore] # from um^2/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = 'D_pore'
    material_property_names = 'hporeAD(pore)'
    expression = '${scale_time}/${scale_length}^2*(hporeAD*${Diffusion_coefficient_D_pore_D0} + (1-hporeAD)*${Diffusion_min})'
  []
  #============================================ Surface reaction rate of tritium
  [ReactionRateSurface_sXYg_ref] # from surface to gaseous for diatomic molecules # from um^3/s/mol to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = ReactionRateSurface_sXYg
    material_property_names = 'hsurfaceAD(pore)'
    expression = '${scale_time}*${scale_quantity}/${scale_length}^3*${reactionRateSurface_sXYg_ref_0}*hsurfaceAD'
  []
  [ReactionRateSurface_gXYs_ref] # from gaseous to surface for diatomic molecules # from 1/s to non-dimensional
    type = ADDerivativeParsedMaterial
    property_name = ReactionRateSurface_gXYs
    coupled_variables = 'pore tritium_s tritium_trap'
    material_property_names = 'hsurfaceAD(pore) theta(tritium_s,tritium_trap)'
    expression = '${scale_time}*${reactionRateSurface_gXYs_ref_0}*hsurfaceAD*(1-theta)^2'
  []
[]

[Postprocessors]
  [Tritium_amount_free]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_s
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_trapped]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_trap
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_ceramic]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_ceramic
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium2g_amount]
    type = ElementIntegralVariablePostprocessor
    variable = tritium_2g
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_pore]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_pore
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tritium_amount_total]
    type = ADElementIntegralMaterialProperty
    mat_prop = tritium_tot
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Surface_tot]
    type = ElementIntegralMaterialProperty
    mat_prop = 1
  []
  [dt]
    type = TimestepSize
  []
[]

[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]
  type = Transient
  scheme = bdf2

  nl_rel_tol = 1e-10
  end_time = 2000
  dtmax = 50
  nl_max_its = 14
  [TimeStepper]
    type = IterationAdaptiveDT
    optimal_iterations = 12
    iteration_window = 1
    growth_factor = 1.2
    dt = 2.37e-7
    cutback_factor = 0.75
    cutback_factor_at_failure = 0.75
  []
[]

[Outputs]
  [exodus]
    type = Exodus
    time_step_interval = 5
  []
  csv = true
  file_base = ${output_name}
[]

General Description of the Simulation Case and Corresponding Input Files
Complete input file
References

Pore-Scale Simulations of Tritium Transport Using TMAP8

General Description of the Simulation Case and Corresponding Input Files

Introduction

Generating the pore microstructure based on an image

Import the pore microstructure image

Perform a phase field simulation to smooth the interfaces

Run simulations

Resulting smooth pore microstructures

Model tritium transport on pore-scale microstructures

Background and Governing Equations

Reaction Rates

Interpolation Functions

Dimensionless Modeling

Calibration Against Experimental Data

Numerical Method, boundary conditions, and initial condition

Run simulations

Results and Discussion

Complete input file

References

(test/tests/pore_scale_transport/2D_microstructure_reader_smoothing_base.i)

(test/tests/pore_scale_transport/pore_structure_closed.params)

(test/tests/pore_scale_transport/pore_structure_open.params)

(test/tests/pore_scale_transport/2D_microstructure_reader_smoothing_base.i)

(test/tests/pore_scale_transport/2D_absorption_base.i)

(test/tests/pore_scale_transport/pore_structure_closed_absorption.params)

(test/tests/pore_scale_transport/pore_structure_open_absorption.params)

(test/tests/pore_scale_transport/2D_absorption_base.i)

(test/tests/pore_scale_transport/2D_microstructure_reader_smoothing_base.i)

(test/tests/pore_scale_transport/2D_absorption_base.i)