SolubilityRatioMaterial

Calculates the jump in concentration across an interface.

Overview

The SolubilityRatioMaterial object is used to calculate the jump in species concentration across a material interface, given different solubility values for those materials and species concentration on either side of the interface. Solubilities for the "primary" and "secondary" sides can be provided via the Materials System using automatic differentiation (a.k.a. an ADMaterialProperty).

The solubility ratio jump is calculated using the following relationship:

var element = document.getElementById("moose-equation-73099b7e-db3e-4e6b-98fd-f93c3ac634f0");katex.render("J = \\frac{c_p}{S_p} - \\frac{c_s}{S_s}", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-b65bfd6b-30d5-418e-b629-3c2ebe3810a5");katex.render("J", element, {displayMode:false,throwOnError:false});$ is the calculated solubility ratio jump (available as an InterfaceMaterial property, named the solubility_ratio), $var element = document.getElementById("moose-equation-8eae25d9-cf6a-499d-8d9e-d65b70462756");katex.render("c_i", element, {displayMode:false,throwOnError:false});$ is the species concentration on the $var element = document.getElementById("moose-equation-3ae2f161-7f9b-4c67-9d25-2db09de0b1eb");katex.render("i", element, {displayMode:false,throwOnError:false});$ th side of the interface, and $var element = document.getElementById("moose-equation-af313752-db94-40c7-bb09-551e0d32e993");katex.render("S_i", element, {displayMode:false,throwOnError:false});$ is the solubility on the $var element = document.getElementById("moose-equation-a32ee4d9-b014-4d59-a850-2da9481c60f5");katex.render("i", element, {displayMode:false,throwOnError:false});$ th side of the interface. The $var element = document.getElementById("moose-equation-25cd8c83-1b68-4afb-ac7f-dfcbf2c2b9d5");katex.render("p", element, {displayMode:false,throwOnError:false});$ subscript corresponds to the primary side and $var element = document.getElementById("moose-equation-a1de3cfc-a46a-4801-8145-e9a0e981bbd6");katex.render("s", element, {displayMode:false,throwOnError:false});$ corresponds to the secondary side. The ratio jump material property can then be used in InterfaceKernels, such as ADPenaltyInterfaceDiffusion, to solve for the appropriate concentrations on either side of the interface.

Example Input File Syntax

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [interface_jump]
    type = SolubilityRatioMaterial<<<{"description": "Calculates the jump in concentration across an interface.", "href": "SolubilityRatioMaterial.html"}>>>
    solubility_primary<<<{"description": "The material property on the primary side of the interface"}>>> = solubility_BeO
    solubility_secondary<<<{"description": "The material property on the secondary side of the interface"}>>> = solubility_Be
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = interface
    concentration_primary<<<{"description": "Primary side non-linear variable for jump computation"}>>> = deuterium_concentration_BeO
    concentration_secondary<<<{"description": "Secondary side non-linear variable for jump computation"}>>> = deuterium_concentration_Be
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'exodus'
  []
[]

Input Parameters

concentration_primaryPrimary side non-linear variable for jump computation
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Primary side non-linear variable for jump computation
concentration_secondarySecondary side non-linear variable for jump computation
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Secondary side non-linear variable for jump computation
solubility_primaryThe material property on the primary side of the interface
C++ Type:std::string
Controllable:No
Description:The material property on the primary side of the interface
solubility_secondaryThe material property on the secondary side of the interface
C++ Type:std::string
Controllable:No
Description:The material property on the secondary side of the interface

Required Parameters

blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
search_methodnearest_node_connected_sidesChoice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
Default:nearest_node_connected_sides
C++ Type:MooseEnum
Options:nearest_node_connected_sides, all_proximate_sides
Controllable:No
Description:Choice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object

Outputs Parameters

prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Material Property Retrieval Parameters

Input Files

(test/tests/divertor_monoblock/divertor_monoblock_physics.i)
(test/tests/divertor_monoblock/divertor_monoblock_physics-single-variable.i)
(test/tests/divertor_monoblock/divertor_monoblock_explicit_base.i)
(test/tests/val-2b/val-2b.i)

(test/tests/val-2b/val-2b.i)

# Validation Problem #2b from TMAP4/TMAP7 V&V document
# Thermal desorption of deuterium from BeO/Be with time-varying temperature and enclosure pressure
# No trapping terms

# Physical constants
R = ${units 8.31446261815324 J/mol/K} # ideal gas constant based on number used in include/utils/PhysicalConstants.h

# Pressure conditions
pressure_enclosure_init = ${units 13300 Pa}
pressure_enclosure_cooldown = ${units 1e-6 Pa} # vaccum
pressure_enclosure_desorption = ${units 1e-3 Pa} # vaccum, assumed

# Temperature conditions
temperature_initial = ${units 773 K}
temperature_desorption_min = ${units 300 K}
temperature_desorption_max = ${units 1073 K}
temperature_cooldown_min = ${temperature_desorption_min}
desorption_heating_rate = ${units ${fparse 3/60} K/s}

# Important times
charge_time = ${units 50 h -> s}
cooldown_time_constant = ${units ${fparse 45*60} s}
# TMAP4 and TMAP7 used 40 minutes for the cooldown duration,
# We use a 5 hour cooldown period to let the temperature decrease to around 300 K for the start of the desorption.
# R.G. Macaulay-Newcombe et al. (1991) is not very clear on how long samples cooled down.
cooldown_duration = ${units 5 h -> s}
desorption_duration = ${fparse (temperature_desorption_max-temperature_desorption_min)/desorption_heating_rate}
endtime = ${fparse charge_time + cooldown_duration + desorption_duration}

# Materials properties
concentration_scaling = 1e10 # (-)
diffusion_Be_preexponential = ${units 8.0e-9 m^2/s -> mum^2/s}
diffusion_Be_energy = ${units 4220 K}
diffusion_BeO_preexponential_charging = ${units 1.40e-4 m^2/s -> mum^2/s}
diffusion_BeO_energy_charging = ${units 24408 K}
diffusion_BeO_preexponential_desorption = ${units 7e-5 m^2/s -> mum^2/s}
diffusion_BeO_energy_desorption = ${units 27000 K}
solubility_order = .5 # order of the solubility law (Here, we use Sievert's law)
solubility_constant_BeO = ${fparse 5.00e20 / 1e18 / concentration_scaling} # at/m^3/Pa^0.5 -> at/mum^3/Pa^0.5}
solubility_energy_BeO = ${units -9377.7 K}
solubility_constant_Be = ${fparse 7.156e27 / 1e18 / concentration_scaling} # at/m^3/Pa^0.5 -> at/mum^3/Pa^0.5}
solubility_energy_Be = ${units 11606 K}
jump_penalty = 1e0 # (-)

# Numerical parameters
dt_max_large = ${units 100 s}
dt_max_small = ${units 10 s}
dt_start_charging = ${units 1 s}
dt_start_cooldown = ${units 10 s}
dt_start_desorption = ${units 1 s}

# Geometry and mesh
length_BeO = ${units 18 nm -> mum}
num_nodes_BeO = 18
node_length_BeO = ${fparse length_BeO / num_nodes_BeO}
length_Be = ${units 0.4 mm -> mum}
length_Be_modeled = ${fparse length_Be/2}
num_nodes_Be = 40
node_length_Be = ${fparse length_Be_modeled / num_nodes_Be}

[Mesh]
  [cmg]
    type = CartesianMeshGenerator
    dim = 1

    #     0                  1                  2                  3                  4                  5                  6                  7                  8                  9                  10                 11                 12                 13                 14                 15                 16                 17
    dx = '${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO}
          ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be}
          ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be}'
    #     18                19                20                21                22                23                24                25                26                27                28                29                30                31                32                33                34                35                36                37

    #               0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
    subdomain_id = '0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0
                    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
                    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1'
    #               18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37'

  []
  # two surfaces for interface
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = cmg
    primary_block = '0' #BeO
    paired_block = '1' # Be
    new_boundary = 'interface'
  []
  [interface_other_side]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = '1' #BeO
    paired_block = '0' # Be
    new_boundary = 'interface_other'
  []
[]

[Variables]
  # Mobile deuterium concentrations in BeO and Be domains
  [deuterium_concentration_Be] # (atoms/microns^3) / concentration_scaling
    block = 1
  []
  [deuterium_concentration_BeO] # (atoms/microns^3) / concentration_scaling
    block = 0
  []
[]

[AuxVariables]
  # Transient enclosure pressure and temperature used in source/BC functions
  [enclosure_pressure]
    family = SCALAR
    initial_condition = ${pressure_enclosure_init}
  []
  [temperature]
    initial_condition = ${temperature_initial}
  []
  # Flux diagnostic across both materials
  [flux_x]
    order = FIRST
    family = MONOMIAL
  []
[]

[Kernels]
  # Transient diffusion of deuterium in BeO
  [time_BeO]
    type = TimeDerivative
    variable = deuterium_concentration_BeO
    block = 0
  []
  [diffusion_BeO]
    type = ADMatDiffusion
    variable = deuterium_concentration_BeO
    diffusivity = diffusivity_BeO
    block = 0
  []
  # Transient diffusion of deuterium in Be
  [time_Be]
    type = TimeDerivative
    variable = deuterium_concentration_Be
    block = 1
  []
  [diffusion_Be]
    type = ADMatDiffusion
    variable = deuterium_concentration_Be
    diffusivity = diffusivity_Be
    block = 1
  []
[]

[InterfaceKernels]
  # Penalized continuity with solubility jump at BeO/Be interface
  [tied]
    type = ADPenaltyInterfaceDiffusion
    variable = deuterium_concentration_BeO
    neighbor_var = deuterium_concentration_Be
    penalty = ${jump_penalty}
    jump_prop_name = solubility_ratio
    boundary = 'interface'
  []
[]

[AuxScalarKernels]
  # Time-dependent enclosure pressure driving equilibrium BC
  [enclosure_pressure_aux]
    type = FunctionScalarAux
    variable = enclosure_pressure
    function = enclosure_pressure_func
  []
[]

[AuxKernels]
  # Temperature schedule applied as an auxiliary variable and used in functions
  [temperature_aux]
    type = FunctionAux
    variable = temperature
    function = temperature_bc_func
    execute_on = 'INITIAL LINEAR'
  []
  # Flux reconstruction for Be and BeO regions
  [flux_x_Be]
    type = DiffusionFluxAux
    diffusivity = diffusivity_Be
    variable = flux_x
    diffusion_variable = deuterium_concentration_Be
    component = x
    block = 1
  []
  [flux_x_BeO]
    type = DiffusionFluxAux
    diffusivity = diffusivity_BeO
    variable = flux_x
    diffusion_variable = deuterium_concentration_BeO
    component = x
    block = 0
  []
[]

[BCs]
  # Sieverts-law equilibrium on the left surface; symmetry on the right
  [left_flux]
    type = EquilibriumBC
    Ko = ${solubility_constant_BeO}
    activation_energy = '${fparse solubility_energy_BeO * R}'
    boundary = left
    enclosure_var = enclosure_pressure
    temperature = temperature
    variable = deuterium_concentration_BeO
    p = ${solubility_order}
  []
  [right_flux]
    type = ADNeumannBC
    boundary = right
    variable = deuterium_concentration_Be
    value = 0
  []
[]

[Functions]
  # Temperature, pressure, and property functions driving modeling
  [temperature_bc_func]
    type = ParsedFunction
    expression = 'if(t<${charge_time}, ${temperature_initial}, if(t<${fparse charge_time + cooldown_duration}, ${temperature_initial}-((1-exp(-(t-${charge_time})/${cooldown_time_constant}))*${fparse temperature_initial - temperature_cooldown_min}), ${temperature_desorption_min}+${desorption_heating_rate}*(t-${fparse charge_time + cooldown_duration})))'
  []
  [diffusivity_BeO_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = 'if(t<${fparse charge_time + cooldown_duration}, ${diffusion_BeO_preexponential_charging}*exp(-${diffusion_BeO_energy_charging}/T), ${diffusion_BeO_preexponential_desorption}*exp(-${diffusion_BeO_energy_desorption}/T))'
  []
  [diffusivity_Be_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${diffusion_Be_preexponential}*exp(-${diffusion_Be_energy}/T)'
  []
  [enclosure_pressure_func]
    type = ParsedFunction
    expression = 'if(t<${charge_time}, ${pressure_enclosure_init}, if(t<${fparse charge_time + cooldown_duration}, ${pressure_enclosure_cooldown}, ${pressure_enclosure_desorption}))'
  []
  [solubility_BeO_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${solubility_constant_BeO} * exp(-${solubility_energy_BeO}/T)'
  []
  [solubility_Be_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${solubility_constant_Be} * exp(-${solubility_energy_Be}/T)'
  []
  [max_time_step_size_func]
    type = ParsedFunction
    expression = 'if(t < ${fparse charge_time-dt_max_large}, ${dt_max_large}, ${dt_max_small})'
  []
[]

[Materials]
  # Temperature-dependent diffusivity and solubility materials plus interface jump ratio
  [diffusion_solubility]
    type = ADGenericFunctionMaterial
    prop_names = 'diffusivity_BeO diffusivity_Be solubility_Be solubility_BeO'
    prop_values = 'diffusivity_BeO_func diffusivity_Be_func solubility_Be_func solubility_BeO_func'
    outputs = 'exodus'
  []
  [converter_to_nonAD]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_Be diffusivity_BeO'
    reg_props_out = 'diffusivity_Be_nonAD diffusivity_BeO_nonAD'
    outputs = 'exodus'
  []
  [interface_jump]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_BeO
    solubility_secondary = solubility_Be
    boundary = interface
    concentration_primary = deuterium_concentration_BeO
    concentration_secondary = deuterium_concentration_Be
    outputs = 'exodus'
  []
[]

[Postprocessors]
  # Flux
  [avg_flux_left]
    type = SideDiffusiveFluxAverage
    variable = deuterium_concentration_BeO
    boundary = left
    diffusivity = diffusivity_BeO_nonAD
  []
  [avg_flux_total] # total flux coming out of the sample in atoms/microns^2/s
    type = ScalePostprocessor
    value = avg_flux_left
    scaling_factor = '${fparse 2 * concentration_scaling}'
    # Factor of 2 because symmetry is assumed and only one-half of the specimen is modeled.
    # Thus, the total flux coming out of the specimen (per unit area)
    # is twice the flux calculated at the left side of the domain.
    # The 'concentration_scaling' parameter is used to get a consistent concentration unit
  []
  # material properties
  [temperature]
    type = ElementAverageValue
    block = 0
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [diffusion_Be]
    type = ElementAverageValue
    block = 1
    variable = diffusivity_Be
    outputs = 'exodus csv'
  []
  [diffusion_BeO]
    type = ElementAverageValue
    block = 0
    variable = diffusivity_BeO
    outputs = 'exodus csv'
  []
  [solubility_Be]
    type = ElementAverageValue
    block = 1
    variable = solubility_Be
    outputs = 'exodus csv'
  []
  [solubility_BeO]
    type = ElementAverageValue
    block = 0
    variable = solubility_BeO
    outputs = 'exodus csv'
  []
  [gold_solubility_ratio]
    type = ParsedPostprocessor
    pp_names = 'solubility_BeO solubility_Be'
    expression = 'solubility_BeO / solubility_Be'
  []
  [BeO_interface]
    type = SideAverageValue
    boundary = interface
    variable = deuterium_concentration_BeO
  []
  [Be_interface]
    type = SideAverageValue
    boundary = interface_other
    variable = deuterium_concentration_Be
  []
  [variable_ratio]
    type = ParsedPostprocessor
    pp_names = 'BeO_interface Be_interface'
    expression = 'BeO_interface / Be_interface'
  []
  [dt]
    type = TimestepSize
  []
  [max_time_step_size_pp]
    type = FunctionValuePostprocessor
    function = max_time_step_size_func
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = none
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  nl_rel_tol = 1e-6
  nl_abs_tol = 5e-12
  end_time = ${endtime}
  automatic_scaling = true
  compute_scaling_once = false
  nl_max_its = 7
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${dt_start_charging}
    optimal_iterations = 5
    growth_factor = 1.1
    cutback_factor_at_failure = .9
    timestep_limiting_postprocessor = max_time_step_size_pp
    time_t = ' 0                    ${charge_time}        ${fparse charge_time + cooldown_duration}'
    time_dt = '${dt_start_charging} ${dt_start_cooldown}  ${dt_start_desorption}'
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[Outputs]
  csv = true
  [exodus]
    type = Exodus
    output_material_properties = true
  []
[]

(test/tests/divertor_monoblock/divertor_monoblock_physics.i)

### This input file is the complete input file for the divertor monoblock case using the physics
### syntax with several variable.
### This input uses the `!include` feature to incorporate other input files

### Nomenclatures
### C_mobile_j      mobile H concentration in "j" material, where j = CuCrZr, Cu, W
### C_trapped_j     trapped H concentration in "j" material, where j = CuCrZr, Cu, W
### C_total_j       total H concentration in "j" material, where j = CuCrZr, Cu, W
###
### S_empty_j       empty site concentration in "j" material, where j = CuCrZr, Cu, W
### S_trapped_j     trapped site concentration in "j" material, where j = CuCrZr, Cu, W
### S_total_j       total site H concentration in "j" material, where j = CuCrZr, Cu, W
###
### F_permeation    permeation flux
### F_recombination recombination flux
###
### Sc_             Scaled
### Int_            Integrated
### ScInt_          Scaled and integrated

C_mobile_CuCrZr_DirichletBC_Coolant = 1.0e-18 # at.fraction
C_mobile_W_init = 1.0e-20 # at.fraction
C_mobile_Cu_init = 5.0e-17 # at.fraction
C_mobile_CuCrZr_init = 1.0e-15 # at.fraction

# include sections of the input file shared with other inputs
!include divertor_monoblock_common_base.i
!include divertor_monoblock_mesh_base.i
!include divertor_monoblock_output_base.i
!include divertor_monoblock_multi_variable_base.i

[GlobalParams]
  species_scaling_factors = '1'
  # The heat conduction physics preconditioning does not apply well across the other physics
  preconditioning = 'defer'
[]

[Physics]
  [HeatConduction]
    [all]
      temperature_name = 'temperature'
      initial_temperature = ${temperature_initial}
      # if using AD, increase the size of the factorization space in petsc options
      # using -mat_mumps_icntl_14 300 or use superlu_dist over mumps
      use_automatic_differentiation = false

      thermal_conductivity = 'thermal_conductivity'
      specific_heat = 'specific_heat'

      heat_flux_boundaries = 'top'
      boundary_heat_fluxes = 'temperature_flux_bc_function'
      fixed_temperature_boundaries = '2to1'
      boundary_temperatures = 'temperature_inner_function'
    []
  []
  [Diffusion]
    [W]
      variable_name = 'C_mobile_W'
      block = '4'
      diffusivity_matprop = diffusivity_W
      initial_condition = ${C_mobile_W_init}

      neumann_boundaries = 'top'
      boundary_fluxes = 'mobile_flux_bc_function'
    []
    [Cu]
      variable_name = 'C_mobile_Cu'
      block = '3'
      diffusivity_matprop = diffusivity_Cu
      initial_condition = ${C_mobile_Cu_init}
    []
    [CuCrZr]
      variable_name = 'C_mobile_CuCrZr'
      block = '2'
      diffusivity_matprop = diffusivity_CuCrZr
      initial_condition = ${C_mobile_CuCrZr_init}

      dirichlet_boundaries = '2to1'
      boundary_values = '${C_mobile_CuCrZr_DirichletBC_Coolant}'
    []
  []
  [SpeciesTrapping]
    [W]
      species = 'C_trapped_W'
      mobile = 'C_mobile_W'
      block = '4'
      species_initial_concentrations = ${C_trapping_init}
      separate_variables_per_component = false

      temperature = temperature
      alpha_t = ${alpha_t}
      N = ${N_W}
      Ct0 = ${Ct0_W}
      trap_per_free = ${trap_per_free_W}
      trapping_energy = ${trapping_energy}

      alpha_r = ${alpha_r}
      detrapping_energy = ${detrapping_energy_W}
    []
    [Cu]
      species = 'C_trapped_Cu'
      mobile = 'C_mobile_Cu'
      block = '3'
      species_initial_concentrations = ${C_trapping_init}
      separate_variables_per_component = false
      temperature = temperature

      alpha_t = ${alpha_t}
      N = ${N_Cu}
      Ct0 = ${Ct0_Cu}
      trap_per_free = ${trap_per_free_Cu}
      trapping_energy = ${trapping_energy}

      alpha_r = ${alpha_r}
      detrapping_energy = ${detrapping_energy_Cu}
    []
    [CuCrZr]
      species = 'C_trapped_CuCrZr'
      mobile = 'C_mobile_CuCrZr'
      block = '2'
      species_initial_concentrations = ${C_trapping_init}
      separate_variables_per_component = false
      temperature = temperature

      alpha_t = ${alpha_t}
      N = ${N_CuCrZr}
      Ct0 = ${Ct0_CuCrZr}
      trap_per_free = ${trap_per_free_CuCrZr}
      trapping_energy = ${trapping_energy}

      alpha_r = ${alpha_r}
      detrapping_energy = ${detrapping_energy_CuCrZr}
    []
  []
[]

[InterfaceKernels]
  [tied_4to3]
    type = ADPenaltyInterfaceDiffusion
    variable = C_mobile_W
    neighbor_var = C_mobile_Cu
    penalty = 0.05 #  it will not converge with > 0.1, but it creates negative C_mobile _Cu with << 0.1
    # jump_prop_name = solubility_ratio_4to3
    jump_prop_name = solubility_ratio
    boundary = '4to3'
  []
  [tied_3to2]
    type = ADPenaltyInterfaceDiffusion
    variable = C_mobile_Cu
    neighbor_var = C_mobile_CuCrZr
    penalty = 0.05 #  it will not converge with > 0.1, but it creates negative C_mobile _Cu with << 0.1
    # jump_prop_name = solubility_ratio_3to2
    jump_prop_name = solubility_ratio
    boundary = '3to2'
  []
[]

[Materials]
  ############################## Materials for W (block = 4)
  [diffusivity_W]
    type = ADParsedMaterial
    property_name = diffusivity_W
    coupled_variables = 'temperature'
    block = 4
    expression = '${diffusivity_W_D0}*exp(-${diffusivity_W_Ea}/temperature)'
    outputs = all
  []
  [solubility_W]
    type = ADParsedMaterial
    property_name = solubility_W
    coupled_variables = 'temperature'
    block = 4
    expression = '${solubility_W_1_D0}*exp(-${solubility_W_1_Ea}/temperature) + ${solubility_W_2_D0}*exp(-${solubility_W_2_Ea}/temperature)'
    outputs = all
  []
  [heat_transfer_W]
    type = GenericConstantMaterial
    prop_names = 'density'
    prop_values = '${density_W}'
    block = 4
  []
  [specific_heat_W]
    type = ParsedMaterial
    property_name = specific_heat
    coupled_variables = 'temperature'
    block = 4
    expression = '1.16e2 + 7.11e-2 * temperature - 6.58e-5 * temperature^2 + 3.24e-8 * temperature^3 -5.45e-12 * temperature^4' # [J/kg-K]
    outputs = all
  []
  [thermal_conductivity_W]
    type = ParsedMaterial
    property_name = thermal_conductivity
    coupled_variables = 'temperature'
    block = 4
    expression = '2.41e2 - 2.90e-1 * temperature + 2.54e-4 * temperature^2 - 1.03e-7 * temperature^3 + 1.52e-11 * temperature^4' # [W/m-K]
    outputs = all
  []
  ############################## Materials for Cu (block = 3)
  [diffusivity_Cu]
    type = ADParsedMaterial
    property_name = diffusivity_Cu
    coupled_variables = 'temperature'
    block = 3
    expression = '${diffusivity_Cu_D0}*exp(-${diffusivity_Cu_Ea}/temperature)'
    outputs = all
  []
  [solubility_Cu]
    type = ADParsedMaterial
    property_name = solubility_Cu
    coupled_variables = 'temperature'
    block = 3
    expression = '${solubility_Cu_D0}*exp(-${solubility_Cu_Ea}/temperature)'
    outputs = all
  []
  [heat_transfer_Cu]
    type = GenericConstantMaterial
    prop_names = 'density'
    prop_values = '${density_Cu}'
    block = 3
  []
  [specific_heat_Cu]
    type = ParsedMaterial
    property_name = specific_heat
    coupled_variables = 'temperature'
    block = 3
    expression = '3.16e2 + 3.18e-1 * temperature - 3.49e-4 * temperature^2 + 1.66e-7 * temperature^3' # [J/kg-K]
    outputs = all
  []
  [thermal_conductivity_Cu]
    type = ParsedMaterial
    property_name = thermal_conductivity
    coupled_variables = 'temperature'
    block = 3
    # expression = '-3.9e-8 * temperature^3 + 3.8e-5 * temperature^2 - 7.9e-2 * temperature + 4.0e2'    # ~ 401.0  [ W/m-K] from R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038,
    expression = '4.21e2 - 6.85e-2 * temperature' # [W/m-K]
    outputs = all
  []
  ############################## Materials for CuCrZr (block = 2)
  [diffusivity_CuCrZr]
    type = ADParsedMaterial
    property_name = diffusivity_CuCrZr
    coupled_variables = 'temperature'
    block = 2
    expression = '${diffusivity_CuCrZr_D0}*exp(-${diffusivity_CuCrZr_Ea}/temperature)'
    outputs = all
  []
  [solubility_CuCrZr]
    type = ADParsedMaterial
    property_name = solubility_CuCrZr
    coupled_variables = 'temperature'
    block = 2
    expression = '${solubility_CuCrZr_D0}*exp(-${solubility_CuCrZr_Ea}/temperature)'
    outputs = all
  []
  [heat_transfer_CuCrZr]
    type = GenericConstantMaterial
    prop_names = 'density specific_heat'
    prop_values = '${density_CuCrZr} ${specific_heat_CuCrZr}'
    block = 2
  []
  [thermal_conductivity_CuCrZr]
    type = ParsedMaterial
    property_name = thermal_conductivity
    coupled_variables = 'temperature'
    block = 2
    expression = '3.87e2 - 1.28e-1 * temperature' # [W/m-K]
    outputs = all
  []
  ############################## Materials for others
  [interface_jump_4to3]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_W
    solubility_secondary = solubility_Cu
    boundary = '4to3'
    concentration_primary = C_mobile_W
    concentration_secondary = C_mobile_Cu
  []
  [interface_jump_3to2]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_Cu
    solubility_secondary = solubility_CuCrZr
    boundary = '3to2'
    concentration_primary = C_mobile_Cu
    concentration_secondary = C_mobile_CuCrZr
  []
[]

(test/tests/divertor_monoblock/divertor_monoblock_physics-single-variable.i)

### This input file is the complete input file for the divertor monoblock case using the physics
### syntax with a single variable across materials.
### This input uses the `!include` feature to incorporate other input files

### Nomenclatures
### C_mobile        mobile H concentration
### C_trapped       trapped H concentration

C_mobile_init = 1.0e-20 # at.fraction

# Include sections of the input file shared with other inputs
!include divertor_monoblock_common_base.i
!include divertor_monoblock_mesh_base.i
!include divertor_monoblock_output_base.i
!include divertor_monoblock_single_variable_base.i

[GlobalParams]
  species_scaling_factors = '1'
  # The heat conduction physics preconditioning does not apply well across the other physics
  preconditioning = 'defer'
[]

# Define the variable outside of the Physics to prevent the Physics from defining it
# with a block restriction
[Variables]
  [C_trapped]
  []
[]

[Physics]
  [HeatConduction]
    [all]
      temperature_name = 'temperature'
      initial_temperature = ${temperature_initial}

      thermal_conductivity = 'thermal_conductivity'
      specific_heat = 'specific_heat'

      heat_flux_boundaries = 'top'
      boundary_heat_fluxes = 'temperature_flux_bc_function'
      fixed_temperature_boundaries = '2to1'
      boundary_temperatures = 'temperature_inner_function'
    []
  []
  [Diffusion]
    [all]
      variable_name = 'C_mobile'
      diffusivity_matprop = diffusivity
      initial_condition = ${C_mobile_init}

      neumann_boundaries = 'top'
      boundary_fluxes = 'mobile_flux_bc_function'
    []
  []
  [SpeciesTrapping]
    [all]
      species = 'C_trapped'
      mobile = 'C_mobile'
      block = '4'
      species_initial_concentrations = ${C_trapping_init}
      separate_variables_per_component = false

      temperature = temperature
      alpha_t = ${alpha_t}
      N = ${N_W}
      Ct0 = ${Ct0_W}
      trap_per_free = ${trap_per_free_W}
      trapping_energy = ${trapping_energy}

      alpha_r = ${alpha_r}
      detrapping_energy = ${detrapping_energy_W}
    []
    [Cu]
      species = 'C_trapped'
      mobile = 'C_mobile'
      block = '3'
      species_initial_concentrations = ${C_trapping_init}
      separate_variables_per_component = false
      temperature = temperature

      alpha_t = ${alpha_t}
      N = ${N_Cu}
      Ct0 = ${Ct0_Cu}
      trap_per_free = ${trap_per_free_Cu}
      trapping_energy = ${trapping_energy}

      alpha_r = ${alpha_r}
      detrapping_energy = ${detrapping_energy_Cu}
    []
    [CuCrZr]
      species = 'C_trapped'
      mobile = 'C_mobile'
      block = '2'
      species_initial_concentrations = ${C_trapping_init}
      separate_variables_per_component = false
      temperature = temperature

      alpha_t = ${alpha_t}
      N = ${N_CuCrZr}
      Ct0 = ${Ct0_CuCrZr}
      trap_per_free = ${trap_per_free_CuCrZr}
      trapping_energy = ${trapping_energy}

      alpha_r = ${alpha_r}
      detrapping_energy = ${detrapping_energy_CuCrZr}
    []
  []
[]

[Materials]
  ############################## Materials for W (block = 4)
  [diffusivity_W]
    type = ADParsedMaterial
    property_name = diffusivity
    coupled_variables = 'temperature'
    block = 4
    expression = '${diffusivity_W_D0}*exp(-${diffusivity_W_Ea}/temperature)'
    outputs = all
  []
  [solubility_W]
    type = ADParsedMaterial
    property_name = solubility_W
    coupled_variables = 'temperature'
    block = 4
    expression = '${solubility_W_1_D0}*exp(-${solubility_W_1_Ea}/temperature) + ${solubility_W_2_D0}*exp(-${solubility_W_2_Ea}/temperature)'
    outputs = all
  []
  [heat_transfer_W]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = '${density_W}'
    block = 4
  []
  [specific_heat_W]
    type = ADParsedMaterial
    property_name = specific_heat
    coupled_variables = 'temperature'
    block = 4
    expression = '1.16e2 + 7.11e-2 * temperature - 6.58e-5 * temperature^2 + 3.24e-8 * temperature^3 -5.45e-12 * temperature^4' # [J/kg-K]
    outputs = all
  []
  [thermal_conductivity_W]
    type = ADParsedMaterial
    property_name = thermal_conductivity
    coupled_variables = 'temperature'
    block = 4
    expression = '2.41e2 - 2.90e-1 * temperature + 2.54e-4 * temperature^2 - 1.03e-7 * temperature^3 + 1.52e-11 * temperature^4' # [W/m-K]
    outputs = all
  []
  ############################## Materials for Cu (block = 3)
  [diffusivity_Cu]
    type = ADParsedMaterial
    property_name = diffusivity
    coupled_variables = 'temperature'
    block = 3
    expression = '${diffusivity_Cu_D0}*exp(-${diffusivity_Cu_Ea}/temperature)'
    outputs = all
  []
  [solubility_Cu]
    type = ADParsedMaterial
    property_name = solubility_Cu
    coupled_variables = 'temperature'
    block = 3
    expression = '${solubility_Cu_D0}*exp(-${solubility_Cu_Ea}/temperature)'
    outputs = all
  []
  [heat_transfer_Cu]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = '${density_Cu}'
    block = 3
  []
  [specific_heat_Cu]
    type = ADParsedMaterial
    property_name = specific_heat
    coupled_variables = 'temperature'
    block = 3
    expression = '3.16e2 + 3.18e-1 * temperature - 3.49e-4 * temperature^2 + 1.66e-7 * temperature^3' # [J/kg-K]
    outputs = all
  []
  [thermal_conductivity_Cu]
    type = ADParsedMaterial
    property_name = thermal_conductivity
    coupled_variables = 'temperature'
    block = 3
    # expression = '-3.9e-8 * temperature^3 + 3.8e-5 * temperature^2 - 7.9e-2 * temperature + 4.0e2'    # ~ 401.0  [ W/m-K] from R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038,
    expression = '4.21e2 - 6.85e-2 * temperature' # [W/m-K]
    outputs = all
  []
  ############################## Materials for CuCrZr (block = 2)
  [diffusivity_CuCrZr]
    type = ADParsedMaterial
    property_name = diffusivity
    coupled_variables = 'temperature'
    block = 2
    expression = '${diffusivity_CuCrZr_D0}*exp(-${diffusivity_CuCrZr_Ea}/temperature)'
    outputs = all
  []
  [solubility_CuCrZr]
    type = ADParsedMaterial
    property_name = solubility_CuCrZr
    coupled_variables = 'temperature'
    block = 2
    expression = '${solubility_CuCrZr_D0}*exp(-${solubility_CuCrZr_Ea}/temperature)'
    outputs = all
  []
  [heat_transfer_CuCrZr]
    type = ADGenericConstantMaterial
    prop_names = 'density specific_heat'
    prop_values = '${density_CuCrZr} ${specific_heat_CuCrZr}'
    block = 2
  []
  [thermal_conductivity_CuCrZr]
    type = ADParsedMaterial
    property_name = thermal_conductivity
    coupled_variables = 'temperature'
    block = 2
    expression = '3.87e2 - 1.28e-1 * temperature' # [W/m-K]
    outputs = all
  []
  ############################## Materials for others
  [interface_jump_4to3]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_W
    solubility_secondary = solubility_Cu
    boundary = '4to3'
    concentration_primary = C_mobile
    concentration_secondary = C_mobile
  []
  [interface_jump_3to2]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_Cu
    solubility_secondary = solubility_CuCrZr
    boundary = '3to2'
    concentration_primary = C_mobile
    concentration_secondary = C_mobile
  []
[]

(test/tests/divertor_monoblock/divertor_monoblock_explicit_base.i)

### This input file contains everything necessary for the divertor monoblock
### case that is not used in the physics-based inputs
### M. Shimada, P.-C. A. Simon, C. T. Icenhour, and G. Singh, “Toward a high-fidelity
### tritium transport modeling for retention and permeation experiments,” Fusion
### Engineering and Design, Volume 203, 2024, 114438, ISSN 0920-3796,
### https://doi.org/10.1016/j.fusengdes.2024.114438.

### This input uses the `!include` feature to incorporate other input files
### Nomenclatures
### C_mobile_j      mobile H concentration in "j" material, where j = CuCrZr, Cu, W
### C_trapped_j     trapped H concentration in "j" material, where j = CuCrZr, Cu, W
### C_total_j       total H concentration in "j" material, where j = CuCrZr, Cu, W
###
### Sc_             Scaled
### Int_            Integrated
### ScInt_          Scaled and integrated

C_mobile_CuCrZr_DirichletBC_Coolant = 1.0e-18 # at.fraction
C_mobile_W_init = 1.0e-20 # at.fraction
C_mobile_Cu_init = 5.0e-17 # at.fraction
C_mobile_CuCrZr_init = 1.0e-15 # at.fraction

!include divertor_monoblock_mesh_base.i
!include divertor_monoblock_common_base.i
!include divertor_monoblock_multi_variable_base.i

[Functions]
  [temperature_steady_state]
    type = ParsedFunction
    expression = "-1.59786e4*x^2  -1.11629611e4*x + 4.84297313e2 + 1.9491599e6*y^2 + 1.55723201e4*y "
                 "- 7.312884e5*x*y"
  []
[]

[Postprocessors]
  [Tritium_SideFluxIntegral]
    type = SideDiffusiveFluxIntegral
    boundary = '2to1'
    diffusivity = diffusivity_CuCrZr_nonAD
    variable = Sc_C_total_CuCrZr
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [Scaled_Tritium_Flux]
    type = ParsedPostprocessor
    expression = '5.01e-24 * Tritium_SideFluxIntegral'
    pp_names = Tritium_SideFluxIntegral
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [coolant_heat_flux]
    # units of W/m2
    type = SideDiffusiveFluxAverage
    boundary = '2to1'
    diffusivity = thermal_conductivity_CuCrZr
    variable = temperature
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [max_temperature_W]
    type = ElementExtremeValue
    block = 4
    variable = 'temperature'
    value_type = max
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [max_temperature_Cu]
    type = ElementExtremeValue
    block = 3
    variable = 'temperature'
    value_type = max
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [max_temperature_CuCrZr]
    type = ElementExtremeValue
    block = 2
    variable = 'temperature'
    value_type = max
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [avg_temperature_W]
    type = ElementAverageValue
    variable = temperature
    block = 4
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [avg_temperature_Cu]
    type = ElementAverageValue
    variable = temperature
    block = 3
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [avg_temperature_CuCrZr]
    type = ElementAverageValue
    variable = temperature
    block = 2
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [max_concentration_W]
    type = ElementExtremeValue
    variable = 'C_total_W'
    value_type = max
    block = 4
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [max_concentration_Cu]
    type = ElementExtremeValue
    variable = 'C_total_Cu'
    value_type = max
    block = 3
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [max_concentration_CuCrZr]
    type = ElementExtremeValue
    variable = 'C_total_CuCrZr'
    value_type = max
    block = 2
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [area_W]
    type = VolumePostprocessor
    block = 4
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [area_Cu]
    type = VolumePostprocessor
    block = 3
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [area_CuCrZr]
    type = VolumePostprocessor
    block = 2
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []
  [total_retention]
    type = SumPostprocessor
    values = 'ScInt_C_total_W ScInt_C_total_Cu ScInt_C_total_CuCrZr'
    execute_on = 'MULTIAPP_FIXED_POINT_END FINAL'
  []

  # TMAP8 will give an error if the parser does not find a use for defined parameters
  # It also seems that the parser does not deal well with multiple layers of includes,
  # so the definitions of the mesh parameters need to have a use somewhere in this input file
  # This is a dummy postprocessor that is used to satisfy the parser for known terms.
  [unused_parameters]
    type = ParsedPostprocessor
    expression = '${num_sectors} + ${rings_H2O} + ${rings_CuCrZr} + ${rings_Cu} + ${rings_W}'
    enable = False
    outputs = ''
  []
[]

[Problem]
  type = ReferenceResidualProblem
  extra_tag_vectors = 'ref'
  reference_vector = 'ref'
[]

[Variables]
  [temperature]
    order = FIRST
    family = LAGRANGE
  []
  ######################### Variables for W (block = 4)
  [C_mobile_W]
    order = FIRST
    family = LAGRANGE
    initial_condition = ${C_mobile_W_init}
    block = 4
  []
  [C_trapped_W]
    order = FIRST
    family = LAGRANGE
    initial_condition = ${C_trapping_init}
    block = 4
  []
  ######################### Variables for Cu (block = 3)
  [C_mobile_Cu]
    order = FIRST
    family = LAGRANGE
    initial_condition = ${C_mobile_Cu_init}
    block = 3
  []
  [C_trapped_Cu]
    order = FIRST
    family = LAGRANGE
    initial_condition = ${C_trapping_init}
    block = 3
  []
  ######################### Variables for CuCrZr (block = 2)
  [C_mobile_CuCrZr]
    order = FIRST
    family = LAGRANGE
    initial_condition = ${C_mobile_CuCrZr_init}
    block = 2
  []
  [C_trapped_CuCrZr]
    order = FIRST
    family = LAGRANGE
    initial_condition = ${C_trapping_init}
    block = 2
  []
[]

[Kernels]
  ############################## Kernels for W (block = 4)
  [diff_W]
    type = ADMatDiffusion
    variable = C_mobile_W
    diffusivity = diffusivity_W
    block = 4
    extra_vector_tags = ref
  []
  [time_diff_W]
    type = ADTimeDerivative
    variable = C_mobile_W
    block = 4
    extra_vector_tags = ref
  []
  [coupled_time_W]
    type = ScaledCoupledTimeDerivative
    variable = C_mobile_W
    v = C_trapped_W
    block = 4
    extra_vector_tags = ref
  []
  [heat_conduction_W]
    type = HeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity_W
    block = 4
    extra_vector_tags = ref
  []
  [time_heat_conduction_W]
    type = SpecificHeatConductionTimeDerivative
    variable = temperature
    specific_heat = specific_heat_W
    density = density_W
    block = 4
    extra_vector_tags = ref
  []
  ############################## Kernels for Cu (block = 3)
  [diff_Cu]
    type = ADMatDiffusion
    variable = C_mobile_Cu
    diffusivity = diffusivity_Cu
    block = 3
    extra_vector_tags = ref
  []
  [time_diff_Cu]
    type = ADTimeDerivative
    variable = C_mobile_Cu
    block = 3
    extra_vector_tags = ref
  []
  [coupled_time_Cu]
    type = ScaledCoupledTimeDerivative
    variable = C_mobile_Cu
    v = C_trapped_Cu
    block = 3
    extra_vector_tags = ref
  []
  [heat_conduction_Cu]
    type = HeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity_Cu
    block = 3
    extra_vector_tags = ref
  []
  [time_heat_conduction_Cu]
    type = SpecificHeatConductionTimeDerivative
    variable = temperature
    specific_heat = specific_heat_Cu
    density = density_Cu
    block = 3
    extra_vector_tags = ref
  []
  ############################## Kernels for CuCrZr (block = 2)
  [diff_CuCrZr]
    type = ADMatDiffusion
    variable = C_mobile_CuCrZr
    diffusivity = diffusivity_CuCrZr
    block = 2
    extra_vector_tags = ref
  []
  [time_diff_CuCrZr]
    type = ADTimeDerivative
    variable = C_mobile_CuCrZr
    block = 2
    extra_vector_tags = ref
  []
  [coupled_time_CuCrZr]
    type = ScaledCoupledTimeDerivative
    variable = C_mobile_CuCrZr
    v = C_trapped_CuCrZr
    block = 2
    extra_vector_tags = ref
  []
  [heat_conduction_CuCrZr]
    type = HeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity_CuCrZr
    block = 2
    extra_vector_tags = ref
  []
  [time_heat_conduction_CuCrZr]
    type = SpecificHeatConductionTimeDerivative
    variable = temperature
    specific_heat = specific_heat_CuCrZr
    density = density_CuCrZr
    block = 2
    extra_vector_tags = ref
  []
[]

[InterfaceKernels]
  [tied_4to3]
    type = ADPenaltyInterfaceDiffusion
    variable = C_mobile_W
    neighbor_var = C_mobile_Cu
    penalty = 0.05 #  it will not converge with > 0.1, but it creates negative C_mobile _Cu with << 0.1
    # jump_prop_name = solubility_ratio_4to3
    jump_prop_name = solubility_ratio
    boundary = '4to3'
  []
  [tied_3to2]
    type = ADPenaltyInterfaceDiffusion
    variable = C_mobile_Cu
    neighbor_var = C_mobile_CuCrZr
    penalty = 0.05 #  it will not converge with > 0.1, but it creates negative C_mobile _Cu with << 0.1
    # jump_prop_name = solubility_ratio_3to2
    jump_prop_name = solubility_ratio
    boundary = '3to2'
  []
[]

[NodalKernels]
  ############################## NodalKernels for W (block = 4)
  [time_W]
    type = TimeDerivativeNodalKernel
    variable = C_trapped_W
  []
  [trapping_W]
    type = TrappingNodalKernel
    variable = C_trapped_W
    temperature = temperature
    alpha_t = ${alpha_t}
    N = ${N_W}
    Ct0 = ${Ct0_W}
    trapping_energy = ${trapping_energy}
    trap_per_free = ${trap_per_free_W}
    mobile_concentration = 'C_mobile_W'
    extra_vector_tags = ref
  []
  [release_W]
    type = ReleasingNodalKernel
    alpha_r = ${alpha_r}
    temperature = temperature
    detrapping_energy = ${detrapping_energy_W}
    variable = C_trapped_W
  []
  ############################## NodalKernels for Cu (block = 3)
  [time_Cu]
    type = TimeDerivativeNodalKernel
    variable = C_trapped_Cu
  []
  [trapping_Cu]
    type = TrappingNodalKernel
    variable = C_trapped_Cu
    temperature = temperature
    alpha_t = ${alpha_t}
    N = ${N_Cu}
    Ct0 = ${Ct0_Cu}
    trapping_energy = ${trapping_energy}
    trap_per_free = ${trap_per_free_Cu}
    mobile_concentration = 'C_mobile_Cu'
    extra_vector_tags = ref
  []
  [release_Cu]
    type = ReleasingNodalKernel
    alpha_r = ${alpha_r}
    temperature = temperature
    detrapping_energy = ${detrapping_energy_Cu}
    variable = C_trapped_Cu
  []
  ############################## NodalKernels for CuCrZr (block = 2)
  [time_CuCrZr]
    type = TimeDerivativeNodalKernel
    variable = C_trapped_CuCrZr
  []
  [trapping_CuCrZr]
    type = TrappingNodalKernel
    variable = C_trapped_CuCrZr
    temperature = temperature
    alpha_t = ${alpha_t}
    N = ${N_CuCrZr}
    Ct0 = ${Ct0_CuCrZr}
    trapping_energy = ${trapping_energy}
    trap_per_free = ${trap_per_free_CuCrZr}
    mobile_concentration = 'C_mobile_CuCrZr'
    extra_vector_tags = ref
  []
  [release_CuCrZr]
    type = ReleasingNodalKernel
    alpha_r = ${alpha_r}
    temperature = temperature
    detrapping_energy = ${detrapping_energy_CuCrZr}
    variable = C_trapped_CuCrZr
  []
[]

[BCs]
  [C_mob_W_top_flux]
    type = FunctionNeumannBC
    variable = C_mobile_W
    boundary = 'top'
    function = mobile_flux_bc_function
  []
  [mobile_tube]
    type = DirichletBC
    variable = C_mobile_CuCrZr
    value = ${C_mobile_CuCrZr_DirichletBC_Coolant}
    boundary = '2to1'
  []
  [temperature_top]
    type = FunctionNeumannBC
    variable = temperature
    boundary = 'top'
    function = temperature_flux_bc_function
  []
  [temperature_tube]
    type = FunctionDirichletBC
    variable = temperature
    boundary = '2to1'
    function = temperature_inner_function
  []
[]

[Materials]
  ############################## Materials for W (block = 4)
  [diffusivity_W]
    type = ADParsedMaterial
    property_name = diffusivity_W
    coupled_variables = 'temperature'
    block = 4
    expression = '${diffusivity_W_D0}*exp(-${diffusivity_W_Ea}/temperature)'
    outputs = all
  []
  [solubility_W]
    type = ADParsedMaterial
    property_name = solubility_W
    coupled_variables = 'temperature'
    block = 4
    expression = '${solubility_W_1_D0}*exp(-${solubility_W_1_Ea}/temperature) + ${solubility_W_2_D0}*exp(-${solubility_W_2_Ea}/temperature)'
    outputs = all
  []
  [converter_to_regular_W]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_W'
    reg_props_out = 'diffusivity_W_nonAD'
    block = 4
  []
  [heat_transfer_W]
    type = GenericConstantMaterial
    prop_names = 'density_W'
    prop_values = '${density_W}'
    block = 4
  []
  [specific_heat_W]
    type = ParsedMaterial
    property_name = specific_heat_W
    coupled_variables = 'temperature'
    block = 4
    expression = '1.16e2 + 7.11e-2 * temperature - 6.58e-5 * temperature^2 + 3.24e-8 * temperature^3 -5.45e-12 * temperature^4' # [J/kg-K]
    outputs = all
  []
  [thermal_conductivity_W]
    type = ParsedMaterial
    property_name = thermal_conductivity_W
    coupled_variables = 'temperature'
    block = 4
    expression = '2.41e2 - 2.90e-1 * temperature + 2.54e-4 * temperature^2 - 1.03e-7 * temperature^3 + 1.52e-11 * temperature^4' # [W/m-K]
    outputs = all
  []
  ############################## Materials for Cu (block = 3)
  [diffusivity_Cu]
    type = ADParsedMaterial
    property_name = diffusivity_Cu
    coupled_variables = 'temperature'
    block = 3
    expression = '${diffusivity_Cu_D0}*exp(-${diffusivity_Cu_Ea}/temperature)'
    outputs = all
  []
  [solubility_Cu]
    type = ADParsedMaterial
    property_name = solubility_Cu
    coupled_variables = 'temperature'
    block = 3
    expression = '${solubility_Cu_D0}*exp(-${solubility_Cu_Ea}/temperature)'
    outputs = all
  []
  [converter_to_regular_Cu]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_Cu'
    reg_props_out = 'diffusivity_Cu_nonAD'
    block = 3
  []
  [heat_transfer_Cu]
    type = GenericConstantMaterial
    prop_names = 'density_Cu'
    prop_values = '${density_Cu}'
    block = 3
  []
  [specific_heat_Cu]
    type = ParsedMaterial
    property_name = specific_heat_Cu
    coupled_variables = 'temperature'
    block = 3
    expression = '3.16e2 + 3.18e-1 * temperature - 3.49e-4 * temperature^2 + 1.66e-7 * temperature^3' # [J/kg-K]
    outputs = all
  []
  [thermal_conductivity_Cu]
    type = ParsedMaterial
    property_name = thermal_conductivity_Cu
    coupled_variables = 'temperature'
    block = 3
    # expression = '-3.9e-8 * temperature^3 + 3.8e-5 * temperature^2 - 7.9e-2 * temperature + 4.0e2'    # ~ 401.0  [ W/m-K] from R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038,
    expression = '4.21e2 - 6.85e-2 * temperature' # [W/m-K]
    outputs = all
  []
  ############################## Materials for CuCrZr (block = 2)
  [diffusivity_CuCrZr]
    type = ADParsedMaterial
    property_name = diffusivity_CuCrZr
    coupled_variables = 'temperature'
    block = 2
    expression = '${diffusivity_CuCrZr_D0}*exp(-${diffusivity_CuCrZr_Ea}/temperature)'
    outputs = all
  []
  [solubility_CuCrZr]
    type = ADParsedMaterial
    property_name = solubility_CuCrZr
    coupled_variables = 'temperature'
    block = 2
    expression = '${solubility_CuCrZr_D0}*exp(-${solubility_CuCrZr_Ea}/temperature)'
    outputs = all
  []
  [converter_to_regular_CuCrZr]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_CuCrZr'
    reg_props_out = 'diffusivity_CuCrZr_nonAD'
    block = 2
  []
  [heat_transfer_CuCrZr]
    type = GenericConstantMaterial
    prop_names = 'density_CuCrZr specific_heat_CuCrZr'
    prop_values = '${density_CuCrZr} ${specific_heat_CuCrZr}'
    block = 2
  []
  [thermal_conductivity_CuCrZr]
    type = ParsedMaterial
    property_name = thermal_conductivity_CuCrZr
    coupled_variables = 'temperature'
    block = 2
    expression = '3.87e2 - 1.28e-1 * temperature' # [W/m-K]
    outputs = all
  []
  ############################## Materials for others
  [interface_jump_4to3]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_W
    solubility_secondary = solubility_Cu
    boundary = '4to3'
    concentration_primary = C_mobile_W
    concentration_secondary = C_mobile_Cu
  []
  [interface_jump_3to2]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_Cu
    solubility_secondary = solubility_CuCrZr
    boundary = '3to2'
    concentration_primary = C_mobile_Cu
    concentration_secondary = C_mobile_CuCrZr
  []
[]

(test/tests/val-2b/val-2b.i)

# Validation Problem #2b from TMAP4/TMAP7 V&V document
# Thermal desorption of deuterium from BeO/Be with time-varying temperature and enclosure pressure
# No trapping terms

# Physical constants
R = ${units 8.31446261815324 J/mol/K} # ideal gas constant based on number used in include/utils/PhysicalConstants.h

# Pressure conditions
pressure_enclosure_init = ${units 13300 Pa}
pressure_enclosure_cooldown = ${units 1e-6 Pa} # vaccum
pressure_enclosure_desorption = ${units 1e-3 Pa} # vaccum, assumed

# Temperature conditions
temperature_initial = ${units 773 K}
temperature_desorption_min = ${units 300 K}
temperature_desorption_max = ${units 1073 K}
temperature_cooldown_min = ${temperature_desorption_min}
desorption_heating_rate = ${units ${fparse 3/60} K/s}

# Important times
charge_time = ${units 50 h -> s}
cooldown_time_constant = ${units ${fparse 45*60} s}
# TMAP4 and TMAP7 used 40 minutes for the cooldown duration,
# We use a 5 hour cooldown period to let the temperature decrease to around 300 K for the start of the desorption.
# R.G. Macaulay-Newcombe et al. (1991) is not very clear on how long samples cooled down.
cooldown_duration = ${units 5 h -> s}
desorption_duration = ${fparse (temperature_desorption_max-temperature_desorption_min)/desorption_heating_rate}
endtime = ${fparse charge_time + cooldown_duration + desorption_duration}

# Materials properties
concentration_scaling = 1e10 # (-)
diffusion_Be_preexponential = ${units 8.0e-9 m^2/s -> mum^2/s}
diffusion_Be_energy = ${units 4220 K}
diffusion_BeO_preexponential_charging = ${units 1.40e-4 m^2/s -> mum^2/s}
diffusion_BeO_energy_charging = ${units 24408 K}
diffusion_BeO_preexponential_desorption = ${units 7e-5 m^2/s -> mum^2/s}
diffusion_BeO_energy_desorption = ${units 27000 K}
solubility_order = .5 # order of the solubility law (Here, we use Sievert's law)
solubility_constant_BeO = ${fparse 5.00e20 / 1e18 / concentration_scaling} # at/m^3/Pa^0.5 -> at/mum^3/Pa^0.5}
solubility_energy_BeO = ${units -9377.7 K}
solubility_constant_Be = ${fparse 7.156e27 / 1e18 / concentration_scaling} # at/m^3/Pa^0.5 -> at/mum^3/Pa^0.5}
solubility_energy_Be = ${units 11606 K}
jump_penalty = 1e0 # (-)

# Numerical parameters
dt_max_large = ${units 100 s}
dt_max_small = ${units 10 s}
dt_start_charging = ${units 1 s}
dt_start_cooldown = ${units 10 s}
dt_start_desorption = ${units 1 s}

# Geometry and mesh
length_BeO = ${units 18 nm -> mum}
num_nodes_BeO = 18
node_length_BeO = ${fparse length_BeO / num_nodes_BeO}
length_Be = ${units 0.4 mm -> mum}
length_Be_modeled = ${fparse length_Be/2}
num_nodes_Be = 40
node_length_Be = ${fparse length_Be_modeled / num_nodes_Be}

[Mesh]
  [cmg]
    type = CartesianMeshGenerator
    dim = 1

    #     0                  1                  2                  3                  4                  5                  6                  7                  8                  9                  10                 11                 12                 13                 14                 15                 16                 17
    dx = '${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO}
          ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be}
          ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be}'
    #     18                19                20                21                22                23                24                25                26                27                28                29                30                31                32                33                34                35                36                37

    #               0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
    subdomain_id = '0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0
                    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
                    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1'
    #               18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37'

  []
  # two surfaces for interface
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = cmg
    primary_block = '0' #BeO
    paired_block = '1' # Be
    new_boundary = 'interface'
  []
  [interface_other_side]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = '1' #BeO
    paired_block = '0' # Be
    new_boundary = 'interface_other'
  []
[]

[Variables]
  # Mobile deuterium concentrations in BeO and Be domains
  [deuterium_concentration_Be] # (atoms/microns^3) / concentration_scaling
    block = 1
  []
  [deuterium_concentration_BeO] # (atoms/microns^3) / concentration_scaling
    block = 0
  []
[]

[AuxVariables]
  # Transient enclosure pressure and temperature used in source/BC functions
  [enclosure_pressure]
    family = SCALAR
    initial_condition = ${pressure_enclosure_init}
  []
  [temperature]
    initial_condition = ${temperature_initial}
  []
  # Flux diagnostic across both materials
  [flux_x]
    order = FIRST
    family = MONOMIAL
  []
[]

[Kernels]
  # Transient diffusion of deuterium in BeO
  [time_BeO]
    type = TimeDerivative
    variable = deuterium_concentration_BeO
    block = 0
  []
  [diffusion_BeO]
    type = ADMatDiffusion
    variable = deuterium_concentration_BeO
    diffusivity = diffusivity_BeO
    block = 0
  []
  # Transient diffusion of deuterium in Be
  [time_Be]
    type = TimeDerivative
    variable = deuterium_concentration_Be
    block = 1
  []
  [diffusion_Be]
    type = ADMatDiffusion
    variable = deuterium_concentration_Be
    diffusivity = diffusivity_Be
    block = 1
  []
[]

[InterfaceKernels]
  # Penalized continuity with solubility jump at BeO/Be interface
  [tied]
    type = ADPenaltyInterfaceDiffusion
    variable = deuterium_concentration_BeO
    neighbor_var = deuterium_concentration_Be
    penalty = ${jump_penalty}
    jump_prop_name = solubility_ratio
    boundary = 'interface'
  []
[]

[AuxScalarKernels]
  # Time-dependent enclosure pressure driving equilibrium BC
  [enclosure_pressure_aux]
    type = FunctionScalarAux
    variable = enclosure_pressure
    function = enclosure_pressure_func
  []
[]

[AuxKernels]
  # Temperature schedule applied as an auxiliary variable and used in functions
  [temperature_aux]
    type = FunctionAux
    variable = temperature
    function = temperature_bc_func
    execute_on = 'INITIAL LINEAR'
  []
  # Flux reconstruction for Be and BeO regions
  [flux_x_Be]
    type = DiffusionFluxAux
    diffusivity = diffusivity_Be
    variable = flux_x
    diffusion_variable = deuterium_concentration_Be
    component = x
    block = 1
  []
  [flux_x_BeO]
    type = DiffusionFluxAux
    diffusivity = diffusivity_BeO
    variable = flux_x
    diffusion_variable = deuterium_concentration_BeO
    component = x
    block = 0
  []
[]

[BCs]
  # Sieverts-law equilibrium on the left surface; symmetry on the right
  [left_flux]
    type = EquilibriumBC
    Ko = ${solubility_constant_BeO}
    activation_energy = '${fparse solubility_energy_BeO * R}'
    boundary = left
    enclosure_var = enclosure_pressure
    temperature = temperature
    variable = deuterium_concentration_BeO
    p = ${solubility_order}
  []
  [right_flux]
    type = ADNeumannBC
    boundary = right
    variable = deuterium_concentration_Be
    value = 0
  []
[]

[Functions]
  # Temperature, pressure, and property functions driving modeling
  [temperature_bc_func]
    type = ParsedFunction
    expression = 'if(t<${charge_time}, ${temperature_initial}, if(t<${fparse charge_time + cooldown_duration}, ${temperature_initial}-((1-exp(-(t-${charge_time})/${cooldown_time_constant}))*${fparse temperature_initial - temperature_cooldown_min}), ${temperature_desorption_min}+${desorption_heating_rate}*(t-${fparse charge_time + cooldown_duration})))'
  []
  [diffusivity_BeO_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = 'if(t<${fparse charge_time + cooldown_duration}, ${diffusion_BeO_preexponential_charging}*exp(-${diffusion_BeO_energy_charging}/T), ${diffusion_BeO_preexponential_desorption}*exp(-${diffusion_BeO_energy_desorption}/T))'
  []
  [diffusivity_Be_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${diffusion_Be_preexponential}*exp(-${diffusion_Be_energy}/T)'
  []
  [enclosure_pressure_func]
    type = ParsedFunction
    expression = 'if(t<${charge_time}, ${pressure_enclosure_init}, if(t<${fparse charge_time + cooldown_duration}, ${pressure_enclosure_cooldown}, ${pressure_enclosure_desorption}))'
  []
  [solubility_BeO_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${solubility_constant_BeO} * exp(-${solubility_energy_BeO}/T)'
  []
  [solubility_Be_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${solubility_constant_Be} * exp(-${solubility_energy_Be}/T)'
  []
  [max_time_step_size_func]
    type = ParsedFunction
    expression = 'if(t < ${fparse charge_time-dt_max_large}, ${dt_max_large}, ${dt_max_small})'
  []
[]

[Materials]
  # Temperature-dependent diffusivity and solubility materials plus interface jump ratio
  [diffusion_solubility]
    type = ADGenericFunctionMaterial
    prop_names = 'diffusivity_BeO diffusivity_Be solubility_Be solubility_BeO'
    prop_values = 'diffusivity_BeO_func diffusivity_Be_func solubility_Be_func solubility_BeO_func'
    outputs = 'exodus'
  []
  [converter_to_nonAD]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_Be diffusivity_BeO'
    reg_props_out = 'diffusivity_Be_nonAD diffusivity_BeO_nonAD'
    outputs = 'exodus'
  []
  [interface_jump]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_BeO
    solubility_secondary = solubility_Be
    boundary = interface
    concentration_primary = deuterium_concentration_BeO
    concentration_secondary = deuterium_concentration_Be
    outputs = 'exodus'
  []
[]

[Postprocessors]
  # Flux
  [avg_flux_left]
    type = SideDiffusiveFluxAverage
    variable = deuterium_concentration_BeO
    boundary = left
    diffusivity = diffusivity_BeO_nonAD
  []
  [avg_flux_total] # total flux coming out of the sample in atoms/microns^2/s
    type = ScalePostprocessor
    value = avg_flux_left
    scaling_factor = '${fparse 2 * concentration_scaling}'
    # Factor of 2 because symmetry is assumed and only one-half of the specimen is modeled.
    # Thus, the total flux coming out of the specimen (per unit area)
    # is twice the flux calculated at the left side of the domain.
    # The 'concentration_scaling' parameter is used to get a consistent concentration unit
  []
  # material properties
  [temperature]
    type = ElementAverageValue
    block = 0
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [diffusion_Be]
    type = ElementAverageValue
    block = 1
    variable = diffusivity_Be
    outputs = 'exodus csv'
  []
  [diffusion_BeO]
    type = ElementAverageValue
    block = 0
    variable = diffusivity_BeO
    outputs = 'exodus csv'
  []
  [solubility_Be]
    type = ElementAverageValue
    block = 1
    variable = solubility_Be
    outputs = 'exodus csv'
  []
  [solubility_BeO]
    type = ElementAverageValue
    block = 0
    variable = solubility_BeO
    outputs = 'exodus csv'
  []
  [gold_solubility_ratio]
    type = ParsedPostprocessor
    pp_names = 'solubility_BeO solubility_Be'
    expression = 'solubility_BeO / solubility_Be'
  []
  [BeO_interface]
    type = SideAverageValue
    boundary = interface
    variable = deuterium_concentration_BeO
  []
  [Be_interface]
    type = SideAverageValue
    boundary = interface_other
    variable = deuterium_concentration_Be
  []
  [variable_ratio]
    type = ParsedPostprocessor
    pp_names = 'BeO_interface Be_interface'
    expression = 'BeO_interface / Be_interface'
  []
  [dt]
    type = TimestepSize
  []
  [max_time_step_size_pp]
    type = FunctionValuePostprocessor
    function = max_time_step_size_func
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = none
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  nl_rel_tol = 1e-6
  nl_abs_tol = 5e-12
  end_time = ${endtime}
  automatic_scaling = true
  compute_scaling_once = false
  nl_max_its = 7
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${dt_start_charging}
    optimal_iterations = 5
    growth_factor = 1.1
    cutback_factor_at_failure = .9
    timestep_limiting_postprocessor = max_time_step_size_pp
    time_t = ' 0                    ${charge_time}        ${fparse charge_time + cooldown_duration}'
    time_dt = '${dt_start_charging} ${dt_start_cooldown}  ${dt_start_desorption}'
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[Outputs]
  csv = true
  [exodus]
    type = Exodus
    output_material_properties = true
  []
[]

Overview
Example Input File Syntax
Input Parameters
Input Files