- Ct0The fraction of host sites that can contribute to trapping as a function (-)
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:The fraction of host sites that can contribute to trapping as a function (-)
- NThe atomic number density of the host material (1/m^3)
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The atomic number density of the host material (1/m^3)
- alpha_tThe trapping rate coefficient. This has units of 1/s (e.g. no number densities are involved)
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The trapping rate coefficient. This has units of 1/s (e.g. no number densities are involved)
- mobile_concentrationThe variable representing the mobile concentration of solute particles.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The variable representing the mobile concentration of solute particles.
- temperatureThe temperature (K)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The temperature (K)
- variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on
TrappingNodalKernel
Implements a residual describing the trapping of a species in a material.
Overview
This object implements the residual
where is the trapping rate coefficient, which has dimensions of 1/s, is the atomic number density of the host material, is the concentration of empty trapping sites, and is the concentration of the mobile species. is defined as with being the pre-exponential factor in units of 1/s, being the trapping energy in units of Kelvin, and being the temperature in units of Kelvin.
As outlined in Scaling, it is important to scale different species numerical concentrations to the same order of magnitude in order to have robust (non)linear solves. Unfortunately, incorporation of scaling tends to obfuscate residual-computing objects a bit, including TrappingNodalKernel. Hopefully we can make things clear here, however. Let's imagine that the trapped concentration is 1000 times greater than the mobile concentration. We can bring the trapped concentration to the same numerical level of the mobile concentration by changing the trapped concentration units from #/volume to k#/volume where 'k' denotes kilo. Now let's figure out how this transformation fits into TrappingNodalKernel. TrappingNodalKernel is producing trapped species, so consequently its residual must have units (in our example) of k#/(volume*s). Trapped concentrations have units of k#/volume and the mobile concentration and host density () have units of #/volume. TrappingNodalKernel computes the empty trapping sites concentration, in the following way:
Real empty_trapping_sites = _Ct0.value(_t, (*_current_node)) * _N;
for (const auto i : index_range(_occupancy_concentrations))
empty_trapping_sites -= (*_occupancy_concentrations[i])[_qp] * _occupancy_weights[i];
return -_trapping_rate_coefficient * std::exp(-_trapping_energy / _temperature[_qp]) *
empty_trapping_sites * _mobile_concentration[_qp] / _residual_denominator;
The trapping concentration, in units of k#/volume, is converted to units of #/volume by multiplying by trap_per_free which, in this example of k#/volume trapping concentration and #/volume mobile concentration, has a value of 1000 #/(k#). For the dimensional TrappingNodalKernel, the residual denominator is N * trap_per_free, so the residual is computed with the code
return -_trapping_rate_coefficient * std::exp(-_trapping_energy / _temperature[_qp]) *
empty_trapping_sites * _mobile_concentration[_qp] / _residual_denominator;
Let's carry through the units: 1/s * #/volume * #/volume / (#/volume * 1000 #/(k#)) -> k#/(s*volume) which is exactly the units that we needed.
Input 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
- other_trapped_concentration_variablesOther variables representing trapped particle concentrations.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Other variables representing trapped particle concentrations.
- trap_per_free1An estimate for the ratio of the concentration magnitude of trapped species to free species. Setting a value for this can be helpful in producing a well-scaled matrix (-)
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:An estimate for the ratio of the concentration magnitude of trapped species to free species. Setting a value for this can be helpful in producing a well-scaled matrix (-)
- trapping_energy0The trapping energy (K)
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The trapping energy (K)
Optional Parameters
- absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
- extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
- extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
- matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
- matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
- vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill
Contribution To Tagged Field Data 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.
- diag_save_inThe name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- 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
- save_inThe name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- 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
Input Files
- (test/tests/divertor_monoblock/divertor_monoblock_explicit_base.i)
- (test/tests/val-2g/val-2g_trapping.i)
- (test/tests/ver-1d/ver-1d-diffusion.i)
- (test/tests/ver-1d/ver-1d-trapping_profile.i)
- (test/tests/ver-1dc/ver-1dc_base.i)
- (test/tests/val-2j/val-2j_base.i)
- (test/tests/ver-1jb/ver-1jb.i)
- (test/tests/val-2d/val-2d.i)
- (test/tests/ver-1d/ver-1d-trapping.i)
(src/nodal_kernels/TrappingNodalKernelBase.C)
/************************************************************/
/* DO NOT MODIFY THIS HEADER */
/* TMAP8: Tritium Migration Analysis Program, Version 8 */
/* */
/* Copyright 2021 - 2025 Battelle Energy Alliance, LLC */
/* ALL RIGHTS RESERVED */
/************************************************************/
#include "TrappingNodalKernelBase.h"
#include "Function.h"
InputParameters
TrappingNodalKernelBase::validParams()
{
InputParameters params = NodalKernel::validParams();
params.addParam<Real>("trapping_energy", 0, "The trapping energy (K)");
params.addRequiredParam<Real>("N", "The atomic number density of the host material (1/m^3)");
params.addRequiredParam<FunctionName>(
"Ct0", "The fraction of host sites that can contribute to trapping as a function (-)");
params.addRequiredCoupledVar("mobile_concentration",
"The variable representing the mobile concentration of solute "
"particles.");
params.addCoupledVar("other_trapped_concentration_variables",
"Other variables representing trapped particle concentrations.");
params.addRequiredCoupledVar("temperature", "The temperature (K)");
return params;
}
TrappingNodalKernelBase::TrappingNodalKernelBase(const InputParameters & parameters,
Real trapping_rate_coefficient,
Real residual_denominator)
: NodalKernel(parameters),
_trapping_rate_coefficient(trapping_rate_coefficient),
_trapping_energy(getParam<Real>("trapping_energy")),
_N(getParam<Real>("N")),
_Ct0(getFunction("Ct0")),
_mobile_concentration(coupledValue("mobile_concentration")),
_n_other_concs(coupledComponents("other_trapped_concentration_variables")),
_last_node(nullptr),
_temperature(coupledValue("temperature")),
_residual_denominator(residual_denominator)
{
}
void
TrappingNodalKernelBase::initializeOccupancyTracking(const std::vector<Real> & other_weights,
Real self_weight)
{
mooseAssert(other_weights.size() == _n_other_concs,
"other occupancy weight count must match coupled concentration count");
_occupancy_concentrations.resize(1 + _n_other_concs);
_occupancy_weights.resize(1 + _n_other_concs);
_var_numbers.resize(2 + _n_other_concs);
for (const auto i : make_range(_n_other_concs))
{
_occupancy_concentrations[i] =
&coupledValue("other_trapped_concentration_variables", /*comp=*/i);
_occupancy_weights[i] = other_weights[i];
_var_numbers[i] = coupled("other_trapped_concentration_variables", i);
}
_occupancy_concentrations[_n_other_concs] = &_u;
_occupancy_weights[_n_other_concs] = self_weight;
_var_numbers[_n_other_concs] = _var.number();
_var_numbers[_n_other_concs + 1] = coupled("mobile_concentration");
}
Real
TrappingNodalKernelBase::computeQpResidual()
{
Real empty_trapping_sites = _Ct0.value(_t, (*_current_node)) * _N;
for (const auto i : index_range(_occupancy_concentrations))
empty_trapping_sites -= (*_occupancy_concentrations[i])[_qp] * _occupancy_weights[i];
return -_trapping_rate_coefficient * std::exp(-_trapping_energy / _temperature[_qp]) *
empty_trapping_sites * _mobile_concentration[_qp] / _residual_denominator;
}
void
TrappingNodalKernelBase::ADHelper()
{
if (_current_node == _last_node)
return;
_last_node = _current_node;
LocalDN empty_trapping_sites = _Ct0.value(_t, (*_current_node)) * _N;
for (const auto i : index_range(_occupancy_concentrations))
{
LocalDN trap_conc_dn = (*_occupancy_concentrations[i])[_qp];
trap_conc_dn.derivatives().insert(_var_numbers[i]) = 1.;
empty_trapping_sites -= trap_conc_dn * _occupancy_weights[i];
}
LocalDN mobile_concentration = _mobile_concentration[_qp];
mobile_concentration.derivatives().insert(_var_numbers.back()) = 1.;
_jacobian = -_trapping_rate_coefficient * std::exp(-_trapping_energy / _temperature[_qp]) *
empty_trapping_sites * mobile_concentration / _residual_denominator;
}
Real
TrappingNodalKernelBase::computeQpJacobian()
{
ADHelper();
return _jacobian.derivatives()[_var.number()];
}
Real
TrappingNodalKernelBase::computeQpOffDiagJacobian(unsigned int jvar)
{
if (std::find(_var_numbers.begin(), _var_numbers.end(), jvar) == _var_numbers.end())
return 0;
ADHelper();
return _jacobian.derivatives()[jvar];
}
(src/nodal_kernels/TrappingNodalKernelBase.C)
/************************************************************/
/* DO NOT MODIFY THIS HEADER */
/* TMAP8: Tritium Migration Analysis Program, Version 8 */
/* */
/* Copyright 2021 - 2025 Battelle Energy Alliance, LLC */
/* ALL RIGHTS RESERVED */
/************************************************************/
#include "TrappingNodalKernelBase.h"
#include "Function.h"
InputParameters
TrappingNodalKernelBase::validParams()
{
InputParameters params = NodalKernel::validParams();
params.addParam<Real>("trapping_energy", 0, "The trapping energy (K)");
params.addRequiredParam<Real>("N", "The atomic number density of the host material (1/m^3)");
params.addRequiredParam<FunctionName>(
"Ct0", "The fraction of host sites that can contribute to trapping as a function (-)");
params.addRequiredCoupledVar("mobile_concentration",
"The variable representing the mobile concentration of solute "
"particles.");
params.addCoupledVar("other_trapped_concentration_variables",
"Other variables representing trapped particle concentrations.");
params.addRequiredCoupledVar("temperature", "The temperature (K)");
return params;
}
TrappingNodalKernelBase::TrappingNodalKernelBase(const InputParameters & parameters,
Real trapping_rate_coefficient,
Real residual_denominator)
: NodalKernel(parameters),
_trapping_rate_coefficient(trapping_rate_coefficient),
_trapping_energy(getParam<Real>("trapping_energy")),
_N(getParam<Real>("N")),
_Ct0(getFunction("Ct0")),
_mobile_concentration(coupledValue("mobile_concentration")),
_n_other_concs(coupledComponents("other_trapped_concentration_variables")),
_last_node(nullptr),
_temperature(coupledValue("temperature")),
_residual_denominator(residual_denominator)
{
}
void
TrappingNodalKernelBase::initializeOccupancyTracking(const std::vector<Real> & other_weights,
Real self_weight)
{
mooseAssert(other_weights.size() == _n_other_concs,
"other occupancy weight count must match coupled concentration count");
_occupancy_concentrations.resize(1 + _n_other_concs);
_occupancy_weights.resize(1 + _n_other_concs);
_var_numbers.resize(2 + _n_other_concs);
for (const auto i : make_range(_n_other_concs))
{
_occupancy_concentrations[i] =
&coupledValue("other_trapped_concentration_variables", /*comp=*/i);
_occupancy_weights[i] = other_weights[i];
_var_numbers[i] = coupled("other_trapped_concentration_variables", i);
}
_occupancy_concentrations[_n_other_concs] = &_u;
_occupancy_weights[_n_other_concs] = self_weight;
_var_numbers[_n_other_concs] = _var.number();
_var_numbers[_n_other_concs + 1] = coupled("mobile_concentration");
}
Real
TrappingNodalKernelBase::computeQpResidual()
{
Real empty_trapping_sites = _Ct0.value(_t, (*_current_node)) * _N;
for (const auto i : index_range(_occupancy_concentrations))
empty_trapping_sites -= (*_occupancy_concentrations[i])[_qp] * _occupancy_weights[i];
return -_trapping_rate_coefficient * std::exp(-_trapping_energy / _temperature[_qp]) *
empty_trapping_sites * _mobile_concentration[_qp] / _residual_denominator;
}
void
TrappingNodalKernelBase::ADHelper()
{
if (_current_node == _last_node)
return;
_last_node = _current_node;
LocalDN empty_trapping_sites = _Ct0.value(_t, (*_current_node)) * _N;
for (const auto i : index_range(_occupancy_concentrations))
{
LocalDN trap_conc_dn = (*_occupancy_concentrations[i])[_qp];
trap_conc_dn.derivatives().insert(_var_numbers[i]) = 1.;
empty_trapping_sites -= trap_conc_dn * _occupancy_weights[i];
}
LocalDN mobile_concentration = _mobile_concentration[_qp];
mobile_concentration.derivatives().insert(_var_numbers.back()) = 1.;
_jacobian = -_trapping_rate_coefficient * std::exp(-_trapping_energy / _temperature[_qp]) *
empty_trapping_sites * mobile_concentration / _residual_denominator;
}
Real
TrappingNodalKernelBase::computeQpJacobian()
{
ADHelper();
return _jacobian.derivatives()[_var.number()];
}
Real
TrappingNodalKernelBase::computeQpOffDiagJacobian(unsigned int jvar)
{
if (std::find(_var_numbers.begin(), _var_numbers.end(), jvar) == _var_numbers.end())
return 0;
ADHelper();
return _jacobian.derivatives()[jvar];
}
(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-2g/val-2g_trapping.i)
# Validation Problem #2g for TMAP8
# Deuterium Transport in Proton-Conducting Ceramics
# Diffusion, Trapping, Surface Reaction under Wet and Dry Considered, No Soret effect
# Physical constants
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant based on number used in include/utils/PhysicalConstants.h
eV_to_J = '${units 1.602176634e-19 eV/J}' # eV to J conversion factor based on number used in include/utils/PhysicalConstants.h
N_a = '${units 6.02214076e23 at/mol}' # Avogadro's number based on number used in include/utils/PhysicalConstants.h
k_B = '${units 8.61733e-5 eV/K}' # Boltzmann constant in eV
# thermal parameters
temperature_low = '${units 300 K}'
temperature_initial = '${units 873 K}'
temperature_high = '${units 1400 K}'
temperature_rate = '${units 0.5 K/s}'
# Model parameters
dissolve_duration = '${units 1 h -> s}'
cooldown_time_constant = '${units ${fparse 10*60} s}'
cooldown_duration = '${units 1 h -> s}'
desorption_duration = '${fparse (temperature_high - temperature_low) / temperature_rate}'
endtime = '${units ${fparse dissolve_duration + cooldown_duration + desorption_duration} s}'
dt_start_charging = '${units 1e-4 s}'
bound_value_min = '${units -1e-20 at/mum^3}'
# Geometry and mesh
length = '${units 0.5 mm -> mum}'
edge_number = 200
boundary_mesh = '${units 12 mum}'
num_nodes = 600
Area = '${units ${fparse 7.7e-3 * 2.2e-3} m^2 -> mum^2}'
# Material properties
density_BZY20 = '${units 5.98 g/cm^3 -> g/m^3}'
molar_mass_BZY20 = '${units 276.085 g/mol}'
N = '${units ${fparse density_BZY20 / molar_mass_BZY20 * N_a} at/m^3 -> at/mum^3}' # 1.3043649601e10
# Initial concentrations
OT_concentration_initial = 0
hydration_limit_S = 0.1
oxygen_vacancy_concentration_initial = '${units ${fparse hydration_limit_S / 2 * N} at/mum^3}'
oxygen_concentration_initial = '${units ${fparse 3 * N - oxygen_vacancy_concentration_initial - OT_concentration_initial} at/mum^3}'
electron_concentration_initial = '${units ${fparse 10 ^ electron_concentration_initial_expo * N} at/mum^3}' # 0.001
# Traps parameters
initial_concentration_trap_1 = 0 # (-)
detrapping_energy_1 = '${fparse detrapping_energy_1_ev / k_B}'
trapping_site_fraction_1 = ${fparse 3 * 10 ^ trapping_site_fraction_1_expo} # (-)
trapping_rate_prefactor = '${units ${fparse 4.8 * 10 ^ trapping_rate_prefactor_expo} 1/s}' # 9.1316e12
release_rate_profactor = '${units ${fparse 2.6 * 10 ^ release_rate_profactor_expo }1/s}' # 8.4e12
trapping_energy = '${fparse trapping_energy_ev / k_B}'
trap_per_free_1 = 1e0 # (-)
##### Dry Pressure conditions
pressure_T2_high = '${units 1.33e3 Pa}'
pressure_T2_low = '${units 1e-5 Pa}'
pressure_T2O_dry = '${units 0 Pa}' # We assume the pressure of T2O is 0
##### Wet Pressure conditions
pressure_T2O_high = '${units 2.8e3 Pa}'
pressure_T2O_low = '${units 1e-5 Pa}'
pressure_T2_wet = '${units 0 Pa}' # We assume the pressure of T2 is 0
# chemical_reaction
T2O_reaction_forward_value = '${units ${fparse 2 * 10 ^ T2O_reaction_forward_value_expo} m^4/at/s -> mum^4/at/s}'
T2_reaction_forward_value = '${units ${fparse 2 * 10 ^ T2_reaction_forward_value_expo} m^4/at/s -> mum^4/at/s}'
# Materials diffusivities (Deuterium: diffusivity and solubility data from Hossain 2020)
diffusivity_OT_prefactor = '${units ${fparse diffusivity_OT_prefactor_m2s * sqrt(3/2)} m^2/s -> mum^2/s}'
diffusivity_OT_energy = '${units ${fparse diffusivity_OT_energy_ev * eV_to_J * N_a} J/mol}'
diffusivity_V_O_prefactor = '${units ${diffusivity_V_O_prefactor_m2s} m^2/s -> mum^2/s}'
diffusivity_e_prefactor = '${units ${diffusivity_e_prefactor_m2s} m^2/s -> mum^2/s}'
[Mesh]
active = 'cmg_edge'
[cmg_edge]
type = CartesianMeshGenerator
dim = 1
dx = '${boundary_mesh} ${fparse length - 2 * boundary_mesh} ${boundary_mesh}'
ix = '${edge_number} ${fparse num_nodes - 2*edge_number} ${edge_number}'
subdomain_id = '0 0 0'
[]
[cmg]
type = CartesianMeshGenerator
dim = 1
dx = '${fparse length}'
ix = '${fparse num_nodes}'
subdomain_id = '0'
[]
[]
[Variables]
#### Dry variable
[OT_concentration_dry] # (atoms/microns^3)
initial_condition = ${OT_concentration_initial}
[]
[Oxygen_vacancy_concentration_dry]
initial_condition = ${oxygen_vacancy_concentration_initial}
[]
[electron_concentration_dry]
initial_condition = ${electron_concentration_initial}
[]
[trapped_1_dry]
order = FIRST
family = LAGRANGE
initial_condition = '${fparse initial_concentration_trap_1 * trapping_site_fraction_1 * N}'
[]
#### Wet variable
[OT_concentration_wet] # (atoms/microns^3)
initial_condition = ${OT_concentration_initial}
[]
[Oxygen_vacancy_concentration_wet]
initial_condition = ${oxygen_vacancy_concentration_initial}
[]
[electron_concentration_wet]
initial_condition = ${electron_concentration_initial}
[]
[trapped_1_wet]
order = FIRST
family = LAGRANGE
initial_condition = '${fparse initial_concentration_trap_1 * trapping_site_fraction_1 * N}'
[]
[]
[Bounds]
[concentration_dry_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = OT_concentration_dry
bound_type = lower
bound_value = ${bound_value_min}
[]
[trap1_dry_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_1_dry
bound_type = lower
bound_value = ${bound_value_min}
[]
[concentration_wet_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = OT_concentration_wet
bound_type = lower
bound_value = ${bound_value_min}
[]
[trap1_wet_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_1_wet
bound_type = lower
bound_value = ${bound_value_min}
[]
[]
[AuxVariables]
[bounds_dummy]
order = FIRST
family = LAGRANGE
[]
[temperature]
initial_condition = ${temperature_initial}
[]
#### Dry auxvariable
[pressure_T2_dry]
initial_condition = ${pressure_T2_high}
[]
[Oxygen_concentration_dry]
initial_condition = ${oxygen_concentration_initial}
[]
#### Wet auxvariable
[pressure_T2O_wet]
initial_condition = ${pressure_T2O_high}
[]
[Oxygen_concentration_wet]
initial_condition = ${oxygen_concentration_initial}
[]
[]
[AuxKernels]
[temperature_Aux]
type = FunctionAux
variable = temperature
function = Temperature_function
[]
#### Dry auxkernels
[pressure_T2_dry_Aux]
type = FunctionAux
variable = pressure_T2_dry
function = Pressure_T2_dry_function
[]
[Oxygen_concentration_dry_Aux] # at/mum^3
type = ParsedAux
variable = Oxygen_concentration_dry
coupled_variables = 'Oxygen_vacancy_concentration_dry OT_concentration_dry'
expression = '3 * ${N} - Oxygen_vacancy_concentration_dry - OT_concentration_dry'
[]
#### Wet auxkernels
[pressure_T2O_wet_Aux]
type = FunctionAux
variable = pressure_T2O_wet
function = Pressure_T2O_wet_function
[]
[Oxygen_concentration_wet_Aux] # at/mum^3
type = ParsedAux
variable = Oxygen_concentration_wet
coupled_variables = 'Oxygen_vacancy_concentration_wet OT_concentration_wet'
expression = '3 * ${N} - Oxygen_vacancy_concentration_wet - OT_concentration_wet'
[]
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Kernels]
#### Dry kernels
[time_OT_dry]
type = ADTimeDerivative
variable = OT_concentration_dry
extra_vector_tags = ref
[]
[diffusion_OT_dry]
type = ADMatDiffusion
variable = OT_concentration_dry
diffusivity = diffusivity_OT
extra_vector_tags = ref
[]
[time_V_O_dry]
type = ADTimeDerivative
variable = Oxygen_vacancy_concentration_dry
extra_vector_tags = ref
[]
[diffusion_V_O_dry]
type = ADMatDiffusion
variable = Oxygen_vacancy_concentration_dry
diffusivity = diffusivity_V_O
extra_vector_tags = ref
[]
[time_e_dry]
type = ADTimeDerivative
variable = electron_concentration_dry
extra_vector_tags = ref
[]
[diffusion_e_dry]
type = ADMatDiffusion
variable = electron_concentration_dry
diffusivity = diffusivity_e
extra_vector_tags = ref
[]
# trapping kernel
[coupled_time_trap_1_dry]
type = ADCoefCoupledTimeDerivative
variable = OT_concentration_dry
v = trapped_1_dry
coef = ${trap_per_free_1}
block = 0
extra_vector_tags = ref
[]
#### Wet kernels
[time_OT_wet]
type = ADTimeDerivative
variable = OT_concentration_wet
extra_vector_tags = ref
[]
[diffusion_OT_wet]
type = ADMatDiffusion
variable = OT_concentration_wet
diffusivity = diffusivity_OT
extra_vector_tags = ref
[]
[time_V_O_wet]
type = ADTimeDerivative
variable = Oxygen_vacancy_concentration_wet
extra_vector_tags = ref
[]
[diffusion_V_O_wet]
type = ADMatDiffusion
variable = Oxygen_vacancy_concentration_wet
diffusivity = diffusivity_V_O
extra_vector_tags = ref
[]
[time_e_wet]
type = ADTimeDerivative
variable = electron_concentration_wet
extra_vector_tags = ref
[]
[diffusion_e_wet]
type = ADMatDiffusion
variable = electron_concentration_wet
diffusivity = diffusivity_e
extra_vector_tags = ref
[]
# trapping kernel
[coupled_time_trap_1_wet]
type = ADCoefCoupledTimeDerivative
variable = OT_concentration_wet
v = trapped_1_wet
coef = ${trap_per_free_1}
block = 0
extra_vector_tags = ref
[]
[]
[NodalKernels]
#### First traps under dry
[time_1_dry]
type = TimeDerivativeNodalKernel
variable = trapped_1_dry
[]
[trapping_1_dry]
type = TrappingNodalKernel
variable = trapped_1_dry
mobile_concentration = OT_concentration_dry
alpha_t = '${trapping_rate_prefactor}'
trapping_energy = '${trapping_energy}'
N = '${N}'
Ct0 = '${trapping_site_fraction_1}'
temperature = 'temperature'
trap_per_free = '${trap_per_free_1}'
extra_vector_tags = ref
[]
[release_1_dry]
type = ReleasingNodalKernel
variable = trapped_1_dry
alpha_r = '${release_rate_profactor}'
detrapping_energy = '${detrapping_energy_1}'
temperature = 'temperature'
[]
#### First traps under wet
[time_1_wet]
type = TimeDerivativeNodalKernel
variable = trapped_1_wet
[]
[trapping_1_wet]
type = TrappingNodalKernel
variable = trapped_1_wet
mobile_concentration = OT_concentration_wet
alpha_t = '${trapping_rate_prefactor}'
trapping_energy = '${trapping_energy}'
N = '${N}'
Ct0 = '${trapping_site_fraction_1}'
temperature = 'temperature'
trap_per_free = '${trap_per_free_1}'
extra_vector_tags = ref
[]
[release_1_wet]
type = ReleasingNodalKernel
variable = trapped_1_wet
alpha_r = '${release_rate_profactor}'
detrapping_energy = '${detrapping_energy_1}'
temperature = 'temperature'
[]
[]
[BCs]
#### Dry BCs
[left_OT_dry]
type = ADMatNeumannBC
variable = OT_concentration_dry
boundary = left
value = 1
boundary_material = flux_on_OT_dry
[]
[right_OT_dry]
type = ADMatNeumannBC
variable = OT_concentration_dry
boundary = right
value = 1
boundary_material = flux_on_OT_dry
[]
[left_V_O_dry]
type = ADMatNeumannBC
variable = Oxygen_vacancy_concentration_dry
boundary = left
value = 1
boundary_material = flux_on_V_O_dry
[]
[right_V_O_dry]
type = ADMatNeumannBC
variable = Oxygen_vacancy_concentration_dry
boundary = right
value = 1
boundary_material = flux_on_V_O_dry
[]
[left_e_dry]
type = ADMatNeumannBC
variable = electron_concentration_dry
boundary = left
value = 1
boundary_material = flux_on_e_dry
[]
[right_e_dry]
type = ADMatNeumannBC
variable = electron_concentration_dry
boundary = right
value = 1
boundary_material = flux_on_e_dry
[]
#### Wet BCs
[left_OT_wet]
type = ADMatNeumannBC
variable = OT_concentration_wet
boundary = left
value = 1
boundary_material = flux_on_OT_wet
[]
[right_OT_wet]
type = ADMatNeumannBC
variable = OT_concentration_wet
boundary = right
value = 1
boundary_material = flux_on_OT_wet
[]
[left_V_O_wet]
type = ADMatNeumannBC
variable = Oxygen_vacancy_concentration_wet
boundary = left
value = 1
boundary_material = flux_on_V_O_wet
[]
[right_V_O_wet]
type = ADMatNeumannBC
variable = Oxygen_vacancy_concentration_wet
boundary = right
value = 1
boundary_material = flux_on_V_O_wet
[]
[left_e_wet]
type = ADMatNeumannBC
variable = electron_concentration_wet
boundary = left
value = 1
boundary_material = flux_on_e_wet
[]
[right_e_wet]
type = ADMatNeumannBC
variable = electron_concentration_wet
boundary = right
value = 1
boundary_material = flux_on_e_wet
[]
[]
[Functions]
[Temperature_function]
type = ParsedFunction
expression = 'if(t<${dissolve_duration},
${temperature_initial},
if(t<${dissolve_duration} + ${cooldown_duration},
${temperature_initial}-((1-exp(-(t-${dissolve_duration})/${cooldown_time_constant}))*${fparse temperature_initial - temperature_low}),
if(t<${dissolve_duration} + ${cooldown_duration} + ${desorption_duration},
${temperature_low} + ${temperature_rate} * (t - ${dissolve_duration} - ${cooldown_duration}),
${temperature_high})))'
[]
[Pressure_T2_dry_function]
type = ParsedFunction
expression = 'if(t<${dissolve_duration} + ${cooldown_duration} - 1000, ${pressure_T2_high},
if(t<${dissolve_duration} + ${cooldown_duration}, ${pressure_T2_high} - (1-exp(-(t - ${dissolve_duration} - ${cooldown_duration} + 1000)/10)) * ${fparse pressure_T2_high - pressure_T2_low},
${pressure_T2_low}))'
[]
[Pressure_T2O_wet_function]
type = ParsedFunction
expression = 'if(t<${dissolve_duration} + ${cooldown_duration} - 1000, ${pressure_T2O_high},
if(t<${dissolve_duration} + ${cooldown_duration}, ${pressure_T2O_high} - (1-exp(-(t - ${dissolve_duration} - ${cooldown_duration} + 1000)/10)) * ${fparse pressure_T2O_high - pressure_T2O_low},
${pressure_T2O_low}))'
[]
[max_dt_size_function]
type = ParsedFunction
expression = 'if(t<${dissolve_duration} + 200, 50,
if(t<${dissolve_duration} + ${cooldown_duration} + 100, 10, 10))'
[]
#### Optimization for dry
[T_flux_T2O_dry_function] # T2O * 2
type = ParsedFunction
symbol_names = 'recombination_flux_T2O_dry_left'
symbol_values = 'recombination_flux_T2O_dry_left'
expression = 'if(t<7200, 0,
if(t<9165, 2 * recombination_flux_T2O_dry_left, 0))'
[]
[T_flux_T2_dry_function] # T2 * 2
type = ParsedFunction
symbol_names = 'recombination_flux_T2_dry_left'
symbol_values = 'recombination_flux_T2_dry_left'
expression = 'if(t<7200, 0,
if(t<9165, 2 * recombination_flux_T2_dry_left, 0))'
[]
[experiment_data_T2O_dry_interpolation_function]
type = PiecewiseLinearFromVectorPostprocessor
argument_column = 'Time'
value_column = 'Flux'
vectorpostprocessor_name = experiment_data_D2O_dry
[]
[experiment_data_T2_dry_interpolation_function]
type = PiecewiseLinearFromVectorPostprocessor
argument_column = 'Time'
value_column = 'Flux'
vectorpostprocessor_name = experiment_data_D2_dry
[]
[experiment_data_T2O_dry_scale_function]
type = ParsedFunction
symbol_names = 'experiment_data_T2O_dry_interpolation_function'
symbol_values = 'experiment_data_T2O_dry_interpolation_function'
expression = 'if(t<7200, 0,
if(t<9165, experiment_data_T2O_dry_interpolation_function * ${N_a} / ${Area}, 0))'
[]
[experiment_data_T2_dry_scale_function]
type = ParsedFunction
symbol_names = 'experiment_data_T2_dry_interpolation_function'
symbol_values = 'experiment_data_T2_dry_interpolation_function'
expression = 'if(t<7200, 0,
if(t<9165, experiment_data_T2_dry_interpolation_function * ${N_a} / ${Area}, 0))'
[]
[difference_square_T2O_dry_function]
type = ParsedFunction
symbol_names = 'T_flux_T2O_dry_function experiment_data_T2O_dry_scale_function'
symbol_values = 'T_flux_T2O_dry_function experiment_data_T2O_dry_scale_function'
expression = '(T_flux_T2O_dry_function - experiment_data_T2O_dry_scale_function) ^ 2'
[]
[difference_square_T2_dry_function]
type = ParsedFunction
symbol_names = 'T_flux_T2_dry_function experiment_data_T2_dry_scale_function'
symbol_values = 'T_flux_T2_dry_function experiment_data_T2_dry_scale_function'
expression = '(T_flux_T2_dry_function - experiment_data_T2_dry_scale_function) ^ 2'
[]
#### Optimization for wet
[T_flux_T2O_wet_function] # T2O * 2
type = ParsedFunction
symbol_names = 'recombination_flux_T2O_wet_left'
symbol_values = 'recombination_flux_T2O_wet_left'
expression = 'if(t<7200, 0,
if(t<9165, 2 * recombination_flux_T2O_wet_left, 0))'
[]
[T_flux_T2_wet_function] # T2 * 2
type = ParsedFunction
symbol_names = 'recombination_flux_T2_wet_left'
symbol_values = 'recombination_flux_T2_wet_left'
expression = 'if(t<7200, 0,
if(t<9165, 2 * recombination_flux_T2_wet_left, 0))'
[]
[experiment_data_T2O_wet_interpolation_function]
type = PiecewiseLinearFromVectorPostprocessor
argument_column = 'Time'
value_column = 'Flux'
vectorpostprocessor_name = experiment_data_D2O_wet
[]
[experiment_data_T2_wet_interpolation_function]
type = PiecewiseLinearFromVectorPostprocessor
argument_column = 'Time'
value_column = 'Flux'
vectorpostprocessor_name = experiment_data_D2_wet
[]
[experiment_data_T2O_wet_scale_function]
type = ParsedFunction
symbol_names = 'experiment_data_T2O_wet_interpolation_function'
symbol_values = 'experiment_data_T2O_wet_interpolation_function'
expression = 'if(t<7200, 0,
if(t<9165, experiment_data_T2O_wet_interpolation_function * ${N_a} / ${Area}, 0))'
[]
[experiment_data_T2_wet_scale_function]
type = ParsedFunction
symbol_names = 'experiment_data_T2_wet_interpolation_function'
symbol_values = 'experiment_data_T2_wet_interpolation_function'
expression = 'if(t<7200, 0,
if(t<9165, experiment_data_T2_wet_interpolation_function * ${N_a} / ${Area}, 0))'
[]
[difference_square_T2O_wet_function]
type = ParsedFunction
symbol_names = 'T_flux_T2O_wet_function experiment_data_T2O_wet_scale_function'
symbol_values = 'T_flux_T2O_wet_function experiment_data_T2O_wet_scale_function'
expression = '(T_flux_T2O_wet_function - experiment_data_T2O_wet_scale_function) ^ 2'
[]
[difference_square_T2_wet_function]
type = ParsedFunction
symbol_names = 'T_flux_T2_wet_function experiment_data_T2_wet_scale_function'
symbol_values = 'T_flux_T2_wet_function experiment_data_T2_wet_scale_function'
expression = '(T_flux_T2_wet_function - experiment_data_T2_wet_scale_function) ^ 2'
[]
[]
[Materials]
[diffusivity_OT]
type = ADParsedMaterial
property_name = 'diffusivity_OT'
coupled_variables = 'temperature'
expression = '${diffusivity_OT_prefactor} * exp(-${diffusivity_OT_energy} / ${R} / temperature)'
[]
[diffusivity_V_O]
type = ADParsedMaterial
property_name = 'diffusivity_V_O'
coupled_variables = 'temperature'
expression = '${diffusivity_V_O_prefactor} * exp(-${diffusivity_V_O_energy} / ${R} / temperature)'
[]
[diffusivity_e]
type = ADParsedMaterial
property_name = 'diffusivity_e'
coupled_variables = 'temperature'
expression = '${diffusivity_e_prefactor} * exp(-${diffusivity_e_energy} / ${R} / temperature)'
[]
[converter_to_nonAD]
type = MaterialADConverter
ad_props_in = 'diffusivity_OT diffusivity_V_O diffusivity_e'
reg_props_out = 'diffusivity_OT_nonAD diffusivity_V_O_nonAD diffusivity_e_nonAD'
[]
[reaction_equilibrium_constant_T2O]
type = ADParsedMaterial
property_name = 'T2O_K_eq'
coupled_variables = 'temperature'
expression = 'exp( (temperature * ${delta_S_T2O} - ${delta_H_T2O}) / ${R} / temperature )'
[]
[reaction_forward_T2O]
type = ADParsedMaterial
property_name = 'T2O_K_forward'
expression = '${T2O_reaction_forward_value}'
[]
[reaction_reverse_T2O]
type = ADParsedMaterial
property_name = 'T2O_K_reverse'
material_property_names = 'T2O_K_forward T2O_K_eq'
expression = 'T2O_K_forward / T2O_K_eq'
[]
[reaction_equilibrium_constant_T2]
type = ADParsedMaterial
property_name = 'T2_K_eq'
coupled_variables = 'temperature'
expression = 'exp( (temperature * ${delta_S_T2} - ${delta_H_T2}) / ${R} / temperature )'
[]
[reaction_forward_T2]
type = ADParsedMaterial
property_name = 'T2_K_forward'
expression = '${T2_reaction_forward_value}'
[]
[reaction_reverse_T2]
type = ADParsedMaterial
property_name = 'T2_K_reverse'
material_property_names = 'T2_K_forward T2_K_eq'
expression = 'T2_K_forward / T2_K_eq'
[]
#### Reaction for dry
[flux_base_on_T2_dry] # T2 + 2 O -> 2 OT + 2 e
type = ADDerivativeParsedMaterial
coupled_variables = 'OT_concentration_dry pressure_T2_dry Oxygen_concentration_dry electron_concentration_dry'
property_name = 'flux_base_on_T2_dry'
material_property_names = 'T2_K_forward T2_K_reverse'
expression = '(T2_K_forward * pressure_T2_dry * Oxygen_concentration_dry^2 - T2_K_reverse * OT_concentration_dry^2 * electron_concentration_dry^2)'
[]
[flux_base_on_T2O_dry] # T2O + V_O + O -> 2 OT
type = ADDerivativeParsedMaterial
coupled_variables = 'OT_concentration_dry Oxygen_concentration_dry Oxygen_vacancy_concentration_dry'
property_name = 'flux_base_on_T2O_dry'
material_property_names = 'T2O_K_forward T2O_K_reverse'
expression = '(T2O_K_forward * ${pressure_T2O_dry} * Oxygen_concentration_dry * Oxygen_vacancy_concentration_dry - T2O_K_reverse * OT_concentration_dry^2)'
[]
#### Reaction for wet
[flux_base_on_T2O_wet] # T2O + V_O + O -> 2 OT
type = ADDerivativeParsedMaterial
coupled_variables = 'OT_concentration_wet pressure_T2O_wet Oxygen_concentration_wet Oxygen_vacancy_concentration_wet'
property_name = 'flux_base_on_T2O_wet'
material_property_names = 'T2O_K_forward T2O_K_reverse'
expression = '(T2O_K_forward * pressure_T2O_wet * Oxygen_concentration_wet * Oxygen_vacancy_concentration_wet - T2O_K_reverse * OT_concentration_wet^2)'
[]
[flux_base_on_T2_wet] # T2 + 2 O -> 2 OT + 2 e
type = ADDerivativeParsedMaterial
coupled_variables = 'OT_concentration_wet Oxygen_concentration_wet electron_concentration_wet'
property_name = 'flux_base_on_T2_wet'
material_property_names = 'T2_K_forward T2_K_reverse'
expression = '(T2_K_forward * ${pressure_T2_wet} * Oxygen_concentration_wet^2 - T2_K_reverse * OT_concentration_wet^2 * electron_concentration_wet^2)'
[]
#### Flux for dry
[flux_on_e_dry] # electron
type = ADDerivativeParsedMaterial
property_name = 'flux_on_e_dry'
material_property_names = 'flux_base_on_T2_dry'
expression = '2 * flux_base_on_T2_dry'
[]
[flux_on_OT_dry] # OT
type = ADDerivativeParsedMaterial
property_name = 'flux_on_OT_dry'
material_property_names = 'flux_base_on_T2_dry flux_base_on_T2O_dry'
expression = '2 * flux_base_on_T2_dry + 2 * flux_base_on_T2O_dry'
[]
[flux_on_T2_dry] # T2
type = ADDerivativeParsedMaterial
property_name = 'flux_on_T2_dry'
material_property_names = 'flux_base_on_T2_dry'
expression = '-1 * flux_base_on_T2_dry'
[]
[flux_on_V_O_dry] # V_O
type = ADDerivativeParsedMaterial
property_name = 'flux_on_V_O_dry'
material_property_names = 'flux_base_on_T2O_dry'
expression = '-1 * flux_base_on_T2O_dry'
[]
[flux_on_T2O_dry] # T2O
type = ADDerivativeParsedMaterial
property_name = 'flux_on_T2O_dry'
material_property_names = 'flux_base_on_T2O_dry'
expression = '-1 * flux_base_on_T2O_dry'
[]
#### Flux for wet
[flux_on_e_wet] # electron
type = ADDerivativeParsedMaterial
property_name = 'flux_on_e_wet'
material_property_names = 'flux_base_on_T2_wet'
expression = '2 * flux_base_on_T2_wet'
[]
[flux_on_OT_wet] # OT
type = ADDerivativeParsedMaterial
property_name = 'flux_on_OT_wet'
material_property_names = 'flux_base_on_T2O_wet'
expression = '2 * flux_base_on_T2O_wet'
[]
[flux_on_T2_wet] # T2
type = ADDerivativeParsedMaterial
property_name = 'flux_on_T2_wet'
material_property_names = 'flux_base_on_T2_wet'
expression = '-1 * flux_base_on_T2_wet'
[]
[flux_on_V_O_wet] # V_O
type = ADDerivativeParsedMaterial
property_name = 'flux_on_V_O_wet'
material_property_names = 'flux_base_on_T2O_wet'
expression = '-1 * flux_base_on_T2O_wet'
[]
[flux_on_T2O_wet] # T2O
type = ADDerivativeParsedMaterial
property_name = 'flux_on_T2O_wet'
material_property_names = 'flux_base_on_T2O_wet'
expression = '-1 * flux_base_on_T2O_wet'
[]
[]
[VectorPostprocessors]
#### Experiment data for dry
[experiment_data_D2O_dry]
type = CSVReaderVectorPostprocessor
csv_file = 'gold/BZY_873K_D2_exposed_D2O_flux.csv'
outputs = none
[]
[experiment_data_D2_dry]
type = CSVReaderVectorPostprocessor
csv_file = 'gold/BZY_873K_D2_exposed_D2_flux.csv'
outputs = none
[]
#### Experiment data for wet
[experiment_data_D2O_wet]
type = CSVReaderVectorPostprocessor
csv_file = 'gold/BZY_873K_D2O_exposed_D2O_flux.csv'
outputs = none
[]
[experiment_data_D2_wet]
type = CSVReaderVectorPostprocessor
csv_file = 'gold/BZY_873K_D2O_exposed_D2_flux.csv'
outputs = none
[]
[]
[Postprocessors]
#### Postprocessors for flux under dry and wet
[recombination_flux_T2O_dry_left]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_T2O_dry
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'console csv exodus'
[]
[recombination_flux_T2_dry_left]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_T2_dry
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'console csv exodus'
[]
[recombination_flux_T2O_wet_left]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_T2O_wet
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'console csv exodus'
[]
[recombination_flux_T2_wet_left]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_T2_wet
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'console csv exodus'
[]
# necessary parameters
[T2_K_eq_average]
type = ADElementAverageMaterialProperty
mat_prop = T2_K_eq
[]
[T2_K_forward_average]
type = ADElementAverageMaterialProperty
mat_prop = T2_K_forward
[]
[T2_K_reverse_average]
type = ADElementAverageMaterialProperty
mat_prop = T2_K_reverse
[]
[T2O_K_eq_average]
type = ADElementAverageMaterialProperty
mat_prop = T2O_K_eq
[]
[T2O_K_forward_average]
type = ADElementAverageMaterialProperty
mat_prop = T2O_K_forward
[]
[T2O_K_reverse_average]
type = ADElementAverageMaterialProperty
mat_prop = T2O_K_reverse
[]
[diffusivity_OT_average]
type = ADElementAverageMaterialProperty
mat_prop = diffusivity_OT
execute_on = 'INITIAL TIMESTEP_END'
[]
[diffusivity_V_O_average]
type = ADElementAverageMaterialProperty
mat_prop = diffusivity_V_O
execute_on = 'INITIAL TIMESTEP_END'
[]
[diffusivity_e_average]
type = ADElementAverageMaterialProperty
mat_prop = diffusivity_e
execute_on = 'INITIAL TIMESTEP_END'
[]
[temperature_average]
type = ElementAverageValue
variable = temperature
execute_on = 'INITIAL TIMESTEP_END'
[]
[pressure_T2_average]
type = ElementAverageValue
variable = pressure_T2_dry
execute_on = 'INITIAL TIMESTEP_END'
[]
[pressure_T2O_average]
type = ElementAverageValue
variable = pressure_T2O_wet
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_time_step_size]
type = FunctionValuePostprocessor
function = max_dt_size_function
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[timestep_number_pre]
type = ParsedPostprocessor
pp_names = pp_experiment_data_T2_dry_interpolation
expression = 'if(pp_experiment_data_T2_dry_interpolation > 0, 1, 1e-10)'
use_t = true
execute_on = 'INITIAL TIMESTEP_END'
[]
[timestep_number]
type = CumulativeValuePostprocessor
postprocessor = timestep_number_pre
execute_on = 'INITIAL TIMESTEP_END'
[]
#### Postprocessors optimization for dry
[pp_experiment_data_T2O_dry_interpolation]
type = FunctionValuePostprocessor
function = experiment_data_T2O_dry_scale_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[pp_simulation_data_T2O_dry]
type = FunctionValuePostprocessor
function = T_flux_T2O_dry_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[pp_experiment_data_T2_dry_interpolation]
type = FunctionValuePostprocessor
function = experiment_data_T2_dry_scale_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[pp_simulation_data_T2_dry]
type = FunctionValuePostprocessor
function = T_flux_T2_dry_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[differece_square_T2O_dry]
type = FunctionValuePostprocessor
function = difference_square_T2O_dry_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[differece_square_T2_dry]
type = FunctionValuePostprocessor
function = difference_square_T2_dry_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[sum_difference_square_T2O_dry]
type = CumulativeValuePostprocessor
postprocessor = differece_square_T2O_dry
execute_on = 'INITIAL TIMESTEP_END'
[]
[sum_difference_square_T2_dry]
type = CumulativeValuePostprocessor
postprocessor = differece_square_T2_dry
execute_on = 'INITIAL TIMESTEP_END'
[]
[sum_experiment_data_T2O_dry]
type = CumulativeValuePostprocessor
postprocessor = pp_experiment_data_T2O_dry_interpolation
execute_on = 'INITIAL TIMESTEP_END'
[]
[sum_experiment_data_T2_dry]
type = CumulativeValuePostprocessor
postprocessor = pp_experiment_data_T2_dry_interpolation
execute_on = 'INITIAL TIMESTEP_END'
[]
[RMSPE_T2O_dry]
type = ParsedPostprocessor
pp_names = 'timestep_number sum_difference_square_T2O_dry sum_experiment_data_T2O_dry'
expression = 'sqrt(sum_difference_square_T2O_dry / timestep_number) / (sum_experiment_data_T2O_dry / timestep_number + 1e-10)'
execute_on = 'TIMESTEP_END'
[]
[RMSPE_T2_dry]
type = ParsedPostprocessor
pp_names = 'timestep_number sum_difference_square_T2_dry sum_experiment_data_T2_dry'
expression = 'sqrt(sum_difference_square_T2_dry / timestep_number) / (sum_experiment_data_T2_dry / timestep_number + 1e-10)'
execute_on = 'TIMESTEP_END'
[]
#### Postprocessors optimization for wet
[pp_experiment_data_T2O_wet_interpolation]
type = FunctionValuePostprocessor
function = experiment_data_T2O_wet_scale_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[pp_simulation_data_T2O_wet]
type = FunctionValuePostprocessor
function = T_flux_T2O_wet_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[pp_experiment_data_T2_wet_interpolation]
type = FunctionValuePostprocessor
function = experiment_data_T2_wet_scale_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[pp_simulation_data_T2_wet]
type = FunctionValuePostprocessor
function = T_flux_T2_wet_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[differece_square_T2O_wet]
type = FunctionValuePostprocessor
function = difference_square_T2O_wet_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[differece_square_T2_wet]
type = FunctionValuePostprocessor
function = difference_square_T2_wet_function
execute_on = 'INITIAL TIMESTEP_END'
[]
[sum_difference_square_T2O_wet]
type = CumulativeValuePostprocessor
postprocessor = differece_square_T2O_wet
execute_on = 'INITIAL TIMESTEP_END'
[]
[sum_difference_square_T2_wet]
type = CumulativeValuePostprocessor
postprocessor = differece_square_T2_wet
execute_on = 'INITIAL TIMESTEP_END'
[]
[sum_experiment_data_T2O_wet]
type = CumulativeValuePostprocessor
postprocessor = pp_experiment_data_T2O_wet_interpolation
execute_on = 'INITIAL TIMESTEP_END'
[]
[sum_experiment_data_T2_wet]
type = CumulativeValuePostprocessor
postprocessor = pp_experiment_data_T2_wet_interpolation
execute_on = 'INITIAL TIMESTEP_END'
[]
[RMSPE_T2O_wet]
type = ParsedPostprocessor
pp_names = 'timestep_number sum_difference_square_T2O_wet sum_experiment_data_T2O_wet'
expression = 'sqrt(sum_difference_square_T2O_wet / timestep_number) / (sum_experiment_data_T2O_wet / timestep_number + 1e-10)'
execute_on = 'TIMESTEP_END'
[]
[RMSPE_T2_wet]
type = ParsedPostprocessor
pp_names = 'timestep_number sum_difference_square_T2_wet sum_experiment_data_T2_wet'
expression = 'sqrt(sum_difference_square_T2_wet / timestep_number) / (sum_experiment_data_T2_wet / timestep_number + 1e-10)'
execute_on = 'TIMESTEP_END'
[]
[log_inverse_error]
type = ParsedPostprocessor
pp_names = 'RMSPE_T2O_wet RMSPE_T2_wet RMSPE_T2O_dry RMSPE_T2_dry'
expression = 'if(RMSPE_T2O_wet>0,
if(RMSPE_T2_wet>0,
if(RMSPE_T2O_dry>0,
if(RMSPE_T2_dry>0, log(1 / (RMSPE_T2O_wet + RMSPE_T2_wet + RMSPE_T2O_dry + RMSPE_T2_dry)), -20), -20), -20), -20)'
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Controls]
[stochastic]
type = SamplerReceiver
[]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = NEWTON
petsc_options_iname = '-pc_type -snes_type'
petsc_options_value = 'lu vinewtonrsls'
nl_rel_tol = 1e-7
nl_abs_tol = 1e-10
end_time = ${endtime}
automatic_scaling = true
compute_scaling_once = true
line_search = none
error_on_dtmin = false
nl_max_its = 10
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt_start_charging}
optimal_iterations = 7
growth_factor = 1.1
cutback_factor = 0.9
cutback_factor_at_failure = 0.9
timestep_limiting_postprocessor = max_time_step_size
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Outputs]
[csv]
type = CSV
[]
[exodus]
type = Exodus
start_time = ${fparse dissolve_duration + cooldown_duration}
[]
[]
(test/tests/ver-1d/ver-1d-diffusion.i)
# Verification Problem #1d from TMAP4/TMAP7 V&V document
# Permeation Problem with Trapping in Diffusion-limited Case
# No Soret effect, or solubility included.
# Modeling parameters
node_num = 200
thickness = '${units 1 m}'
end_time = ${units 3 s}
temperature = '${units 1000 K}'
diffusivity = ${units 1 m^2/s}
# Trapping parameters
density = '${units 3.1622e22 at/m^3}'
cl = '${units 3.1622e18 at/m^3}'
trapping_prefactor = ${units 1e15 1/s}
release_prefactor = ${units 1e13 1/s}
release_energy = ${units 100 K}
trapping_fraction = 0.1 # -
[Mesh]
type = GeneratedMesh
dim = 1
nx = ${node_num}
xmax = ${thickness}
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
# mobile tritium variable
[mobile]
[]
# trapped tritium variable
[trapped]
[]
[]
[Kernels]
# kernel for mobile tritium
[diff]
type = MatDiffusion
variable = mobile
diffusivity = ${diffusivity}
extra_vector_tags = ref
[]
[time]
type = TimeDerivative
variable = mobile
extra_vector_tags = ref
[]
# kernel for trapped tritium
[coupled_time]
type = CoupledTimeDerivative
variable = mobile
v = trapped
extra_vector_tags = ref
[]
[]
[NodalKernels]
# kernel for trapped tritium in nodal sites
[time]
type = TimeDerivativeNodalKernel
variable = trapped
[]
[trapping]
type = TrappingNodalKernel
variable = trapped
alpha_t = ${trapping_prefactor}
N = '${fparse density / cl}'
Ct0 = ${trapping_fraction}
mobile_concentration = 'mobile'
temperature = ${temperature}
extra_vector_tags = ref
[]
[release]
type = ReleasingNodalKernel
alpha_r = ${release_prefactor}
temperature = ${temperature}
detrapping_energy = ${release_energy}
variable = trapped
[]
[]
[BCs]
[left]
type = DirichletBC
variable = mobile
value = '${fparse cl / cl}'
boundary = left
[]
[right]
type = DirichletBC
variable = mobile
value = 0
boundary = right
[]
[]
[Postprocessors]
[outflux]
type = SideDiffusiveFluxAverage
boundary = 'right'
diffusivity = ${diffusivity}
variable = mobile
[]
[scaled_outflux]
type = ScalePostprocessor
value = outflux
scaling_factor = ${cl}
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = ${end_time}
dt = .01
dtmin = .01
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
verbose = true
compute_scaling_once = false
[]
[Outputs]
csv = true
exodus = true
[dof]
type = DOFMap
execute_on = initial
[]
perf_graph = true
[]
(test/tests/ver-1d/ver-1d-trapping_profile.i)
# This input file is slightly adapted from ver-1d-trapping.i to include a spatially-dependent trap site density
# Modeling parameters
node_num = 20
thickness = '${units 1 m}'
end_time = '${units 1000 s}'
time_scaling = 1
temperature = '${units 1000 K}'
diffusivity = '${fparse ${units 1 m^2/s} / time_scaling}'
# Trapping parameters
cl = '${units 3.1622e18 at/m^3}'
N = '${units 3.1622e22 at/m^3}'
trapping_prefactor = '${fparse ${units 1e15 1/s} / time_scaling}'
release_prefactor = '${fparse ${units 1e13 1/s} / time_scaling}'
epsilon = '${units 10000 K}'
trap_per_free = 1e3
Ct0_surface = 0.1 # trapping site density close to the sample surface
Ct0_bulk = .01 # trapping site density in the sample bulk
trap_profile_depth = 0.2 # position of the transition from Ct0_surface to Ct0_bulk
[Mesh]
type = GeneratedMesh
dim = 1
nx = ${node_num}
xmax = ${thickness}
[]
[Functions]
# Assign different trapping site density in the sample
[trapping_sites_density_function]
type = ParsedFunction
expression = 'if(x<${trap_profile_depth}, ${Ct0_surface}, ${Ct0_bulk})'
execute_on = 'INITIAL LINEAR NONLINEAR TIMESTEP_END'
[]
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[mobile][]
[trapped][]
[]
[AuxVariables]
[empty_sites][]
[scaled_empty_sites][]
[trapped_sites][]
[total_sites][]
[]
[AuxKernels]
[empty_sites]
variable = empty_sites
type = EmptySitesAux
N = ${fparse N / cl}
Ct0 = trapping_sites_density_function
trap_per_free = ${trap_per_free}
trapped_concentration_variables = trapped
[]
[scaled_empty]
variable = scaled_empty_sites
type = NormalizationAux
normal_factor = ${cl}
source_variable = empty_sites
[]
[trapped_sites]
variable = trapped_sites
type = NormalizationAux
normal_factor = ${trap_per_free}
source_variable = trapped
[]
[total_sites]
variable = total_sites
type = ParsedAux
expression = 'trapped_sites + empty_sites'
coupled_variables = 'trapped_sites empty_sites'
[]
[]
[Kernels]
[diff]
type = MatDiffusion
variable = mobile
diffusivity = ${diffusivity}
extra_vector_tags = ref
[]
[time]
type = TimeDerivative
variable = mobile
extra_vector_tags = ref
[]
[coupled_time]
type = ScaledCoupledTimeDerivative
variable = mobile
v = trapped
factor = ${trap_per_free}
extra_vector_tags = ref
[]
[]
[NodalKernels]
[time]
type = TimeDerivativeNodalKernel
variable = trapped
[]
[trapping]
type = TrappingNodalKernel
variable = trapped
alpha_t = ${trapping_prefactor}
N = ${fparse N / cl}
Ct0 = trapping_sites_density_function
mobile_concentration = 'mobile'
temperature = ${temperature}
trap_per_free = ${trap_per_free}
extra_vector_tags = ref
[]
[release]
type = ReleasingNodalKernel
alpha_r = ${release_prefactor}
temperature = ${temperature}
detrapping_energy = ${epsilon}
variable = trapped
[]
[]
[BCs]
[left]
type = DirichletBC
variable = mobile
value = ${fparse cl / cl}
boundary = left
[]
[right]
type = DirichletBC
variable = mobile
value = 0
boundary = right
[]
[]
[Postprocessors]
[outflux]
type = SideDiffusiveFluxAverage
boundary = 'right'
diffusivity = 1
variable = mobile
[]
[scaled_outflux]
type = ScalePostprocessor
value = outflux
scaling_factor = ${cl}
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = ${end_time}
dt = 1
solve_type = NEWTON
line_search = 'none'
automatic_scaling = true
verbose = true
compute_scaling_once = false
l_tol = 1e-11
[]
[Outputs]
exodus = true
csv = true
[dof]
type = DOFMap
execute_on = initial
[]
perf_graph = true
[]
(test/tests/ver-1dc/ver-1dc_base.i)
# This is the base input file for the ver-1dc case.
# This input file is meant to be included within both the ver-1dc.i
# and ver-1dc_mms.i input files to be complete.
# Its purpose is to centralize the capability common to the two cases
# within one file to minimize redundancy and ease maintenance.
# It is not meant to be run on its own.
epsilon_1 = ${units 100 K}
epsilon_2 = ${units 500 K}
epsilon_3 = ${units 800 K}
temperature = ${units 1000 K}
trapping_site_fraction_1 = 0.10 # (-)
trapping_site_fraction_2 = 0.15 # (-)
trapping_site_fraction_3 = 0.20 # (-)
diffusivity = ${units 1 m^2/s}
[Mesh]
type = GeneratedMesh
dim = 1
nx = ${nx_num}
xmax = 1
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[mobile]
[]
[trapped_1]
[]
[trapped_2]
[]
[trapped_3]
[]
[]
[Kernels]
[diff]
type = MatDiffusion
variable = mobile
diffusivity = ${diffusivity}
extra_vector_tags = ref
[]
[time]
type = TimeDerivative
variable = mobile
extra_vector_tags = ref
[]
[coupled_time_1]
type = CoupledTimeDerivative
variable = mobile
v = trapped_1
extra_vector_tags = ref
[]
[coupled_time_2]
type = CoupledTimeDerivative
variable = mobile
v = trapped_2
extra_vector_tags = ref
[]
[coupled_time_3]
type = CoupledTimeDerivative
variable = mobile
v = trapped_3
extra_vector_tags = ref
[]
[]
[NodalKernels]
# For first traps
[time_1]
type = TimeDerivativeNodalKernel
variable = trapped_1
[]
[trapping_1]
type = TrappingNodalKernel
variable = trapped_1
alpha_t = '${trapping_rate_coefficient}'
N = '${fparse N / cl}'
Ct0 = '${trapping_site_fraction_1}'
mobile_concentration = 'mobile'
temperature = '${temperature}'
extra_vector_tags = ref
[]
[release_1]
type = ReleasingNodalKernel
alpha_r = '${release_rate_coefficient}'
temperature = '${temperature}'
detrapping_energy = '${epsilon_1}'
variable = trapped_1
[]
# For second traps
[time_2]
type = TimeDerivativeNodalKernel
variable = trapped_2
[]
[trapping_2]
type = TrappingNodalKernel
variable = trapped_2
alpha_t = '${trapping_rate_coefficient}'
N = '${fparse N / cl}'
Ct0 = '${trapping_site_fraction_2}'
mobile_concentration = 'mobile'
temperature = '${temperature}'
extra_vector_tags = ref
[]
[release_2]
type = ReleasingNodalKernel
alpha_r = '${release_rate_coefficient}'
temperature = '${temperature}'
detrapping_energy = '${epsilon_2}'
variable = trapped_2
[]
# For third traps
[time_3]
type = TimeDerivativeNodalKernel
variable = trapped_3
[]
[trapping_3]
type = TrappingNodalKernel
variable = trapped_3
alpha_t = '${trapping_rate_coefficient}'
N = '${fparse N / cl}'
Ct0 = '${trapping_site_fraction_3}'
mobile_concentration = 'mobile'
temperature = '${temperature}'
extra_vector_tags = ref
[]
[release_3]
type = ReleasingNodalKernel
alpha_r = '${release_rate_coefficient}'
temperature = '${temperature}'
detrapping_energy = '${epsilon_3}'
variable = trapped_3
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = ${simulation_time}
dtmax = ${time_interval_max}
solve_type = NEWTON
scheme = ${scheme}
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
line_search = 'none'
[]
(test/tests/val-2j/val-2j_base.i)
# Base input file for val-2j: Tritium TDS from Li2TiO3 solid breeder
# This file contains the shared simulation blocks for all val-2j cases.
# The including input file must define:
# D0, E_d, alpha_t, epsilon_t, alpha_r, epsilon_r (diffusion/trapping params)
# alpha_anneal, E_anneal (defect annihilation params)
# ============ Geometry ============
grain_radius = '${units 1.5e-6 m -> mum}' # average grain radius
num_elems = 100 # number of radial elements
# ============ Material Properties ============
# Li2TiO3 lattice density: rho=3.43 g/cm3, MW=109.75 g/mol
N_lattice = '${units 1.88e28 1/m^3 -> 1/mum^3}' # lattice molecule density
Ct0 = 0.018 # total defect fraction: Sample E (Table 1)
trap_per_free = 1 # concentration scaling factor
# ============ TDS Parameters ============
T_start = '${units 300 K}' # initial temperature
heating_rate = '${units ${fparse 5.0/60.0} K/s}' # 5 K/min
end_time = '${fparse (900.0 - 300.0) / ${heating_rate}}' # ramp duration (300 K to 900 K)
# ============ Initial Conditions ============
# Mobile and trapped concentrations start in local trapping/detrapping
# equilibrium at T_start. This avoids an initial transient from an
# imbalance between trapping and detrapping rates.
# Since TDS output is normalized, absolute concentrations don't matter.
C0_trapped = 1.0 # initial trapped concentration (arb. units, uniform)
C0_mobile = '${fparse alpha_r * exp(-epsilon_r / T_start) * C0_trapped / (alpha_t * exp(-epsilon_t / T_start) * Ct0)}'
[Mesh]
type = GeneratedMesh
dim = 1
nx = ${num_elems}
xmax = ${grain_radius}
coord_type = RSPHERICAL
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[mobile]
initial_condition = ${C0_mobile}
[]
[trapped]
initial_condition = ${C0_trapped}
[]
[defect_density]
initial_condition = ${Ct0}
[]
[]
[AuxVariables]
[temperature]
[]
[]
[AuxKernels]
[temperature_aux]
type = FunctionAux
variable = temperature
function = temperature_ramp
execute_on = 'initial timestep_end linear'
[]
[]
[Kernels]
[diffusion]
type = ADMatDiffusion
variable = mobile
diffusivity = diffusivity
extra_vector_tags = ref
[]
[time_mobile]
type = ADTimeDerivative
variable = mobile
extra_vector_tags = ref
[]
[coupled_time_trapped]
type = ScaledCoupledTimeDerivative
variable = mobile
v = trapped
factor = ${trap_per_free}
extra_vector_tags = ref
[]
[]
[NodalKernels]
# --- (O-)-center (hydroxyl) trap ---
[time_trapped]
type = TimeDerivativeNodalKernel
variable = trapped
[]
[trapping]
type = TrappingNodalKernel
variable = trapped
mobile_concentration = mobile
alpha_t = ${alpha_t}
trapping_energy = ${epsilon_t}
N = ${N_lattice}
Ct0 = 'Ct0_func'
temperature = temperature
trap_per_free = ${trap_per_free}
extra_vector_tags = ref
[]
[releasing]
type = ReleasingNodalKernel
variable = trapped
alpha_r = ${alpha_r}
detrapping_energy = ${epsilon_r}
temperature = temperature
extra_vector_tags = ref
[]
[releasing_defect_annealing]
type = ReleasingNodalKernel
variable = trapped
alpha_r = ${alpha_anneal}
detrapping_energy = ${E_anneal}
temperature = temperature
extra_vector_tags = ref
[]
# --- Defect annihilation ODE: dD_id/dt = -k_anneal * D_id ---
[time_defect_density]
type = TimeDerivativeNodalKernel
variable = defect_density
[]
[defect_annihilation]
type = ReleasingNodalKernel
variable = defect_density
alpha_r = ${alpha_anneal}
detrapping_energy = ${E_anneal}
temperature = temperature
extra_vector_tags = ref
[]
[]
[BCs]
# C = 0 at grain surface (fast surface release - Eq. 20)
# No BC at r=0; symmetry is automatic in RSPHERICAL
[surface]
type = ADDirichletBC
variable = mobile
boundary = right
value = 0
[]
[]
[Materials]
[diffusivity_mat]
type = ADDerivativeParsedMaterial
property_name = diffusivity
functor_names = 'temperature_ramp'
functor_symbols = 'T'
expression = '${D0} * exp(-${E_d} / T)'
[]
[converter_to_regular]
type = MaterialADConverter
ad_props_in = 'diffusivity'
reg_props_out = 'diffusivity_nonAD'
outputs = none
[]
[]
[Functions]
[temperature_ramp]
type = ParsedFunction
expression = '${T_start} + ${heating_rate} * t'
[]
[Ct0_func]
type = ParsedFunction
symbol_names = 'Ct0_value'
symbol_values = 'defect_density_pp'
expression = 'Ct0_value'
[]
[]
[Postprocessors]
# Tritium release rate = diffusive flux at grain surface (right boundary)
[release_rate]
type = SideDiffusiveFluxIntegral
variable = mobile
diffusivity = diffusivity_nonAD
boundary = right
[]
[temperature_pp]
type = FunctionValuePostprocessor
function = temperature_ramp
[]
[defect_density_pp]
type = AverageNodalVariableValue
variable = defect_density
execute_on = 'initial timestep_end'
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
line_search = none
automatic_scaling = true
compute_scaling_once = true
end_time = ${end_time}
dtmax = 30 # limit time step to ~2.5 K per step for TDS resolution
[TimeStepper]
type = IterationAdaptiveDT
dt = 1.0
optimal_iterations = 12
growth_factor = 1.2
cutback_factor = 0.8
[]
nl_rel_tol = 1e-10
nl_abs_tol = 1e-12
[]
(test/tests/ver-1jb/ver-1jb.i)
# Verification Problem #1jb from TMAP7 V&V document
# Radioactive Decay of Tritium in a Distributed Trap
# Physical Constants
ideal_gas_constant = ${units 8.31446261815324 J/K/mol} # from PhysicalConstants.h
boltzmann_constant = ${units 1.380649e-23 J/K -> eV/K } # from PhysicalConstants.h
# Case and model parameters (adapted from TMAP7)
slab_length = ${units 1.5 m }
tritium_mobile_concentration_initial = ${units 1 atoms/m3}
trapping_sites_atomic_fraction_max = 0.001 # (-)
trapping_sites_fraction_occupied_initial = 0.5 # (-)
normal_center_position = ${fparse slab_length/2} # m
normal_standard_deviation = ${fparse slab_length/4} # m
density_material = ${units 6.34e28 atoms/m^3} # for tungsten
density_scalar = ${units 1e25 atoms/m^3} # used to scale variables and use a dimensionless system
temperature = ${units 300 K} # assumed (TMAP7's input file lists 273 K)
tritium_diffusivity_prefactor = ${units 1.58e-4 m^2/s} # from TMAP7 V&V input file
tritium_diffusivity_energy = ${units 308000.0 J/K} # from TMAP7 V&V input file
tritium_release_prefactor = ${units 1.0e13 1/s} # from TMAP7 V&V input file
tritium_release_energy = ${units 4.2 eV}
tritium_trapping_prefactor = ${units 2.096e15 1/s} # from TMAP7 V&V input file
tritium_trapping_energy = ${tritium_diffusivity_energy} # J/K - from TMAP7 V&V input file
trap_per_free = 1.e-25 # (-)
half_life_s = ${units 12.3232 year -> s}
decay_rate_constant = ${fparse 0.693/half_life_s} # 1/s
# Simulation parameters
num_mesh_element = 50
end_time = ${units 100 year -> s} # s
dt_start = ${fparse end_time/250} # s
dt_max = ${fparse end_time/100} # s
[Mesh]
type = GeneratedMesh
dim = 1
nx = ${num_mesh_element}
xmax = ${slab_length}
[]
[Variables]
# tritium mobile concentration in atoms/m^3 / density_scalar = (-)
[tritium_mobile_concentration_scaled]
initial_condition = ${fparse tritium_mobile_concentration_initial / density_scalar}
[]
# tritium trapped concentration in (atoms/m^3) / density_scalar = (-)
[tritium_trapped_concentration_scaled]
[]
# helium concentration in (atoms/m^3) / density_scalar = (-)
[helium_concentration_scaled]
[]
[]
[Functions]
[trapping_sites_fraction_function] # (atomic fraction)
type = ParsedFunction
expression = '${trapping_sites_atomic_fraction_max} * exp(-1/2*((x-${normal_center_position})/${normal_standard_deviation})^2)'
[]
[density_material_function] # (atoms/m^3)
type = ConstantFunction
value = ${density_material}
[]
[trapping_sites_concentration_function] # (atoms/m^3)
type = CompositeFunction
functions = 'density_material_function trapping_sites_fraction_function'
[]
[initial_trapping_sites_occupied_function] # (-)
type = ConstantFunction
value = ${trapping_sites_fraction_occupied_initial}
[]
[tritium_trapped_concentration_initial_function] # (atoms/m^3)
type = CompositeFunction
functions = 'initial_trapping_sites_occupied_function trapping_sites_concentration_function'
[]
[density_scalar_inverse_function] # (atoms/m^3)^(-1)
type = ConstantFunction
value = ${fparse 1/density_scalar}
[]
[tritium_trapped_concentration_initial_function_scaled] # (-)
type = CompositeFunction
functions = 'density_scalar_inverse_function tritium_trapped_concentration_initial_function'
[]
[]
[ICs]
[tritium_trapped_concentration_IC]
type = FunctionIC
variable = tritium_trapped_concentration_scaled
function = 'tritium_trapped_concentration_initial_function_scaled'
[]
[]
[AuxVariables]
[empty_sites] # (-)
[]
[scaled_empty_sites] # atoms/m^3
[]
[trapped_sites] # (atoms/m^3)
[]
[total_sites] # (atoms/m^3)
[]
[temperature] # (K)
initial_condition = ${temperature}
[]
[tritium_mobile_concentration] # (atoms/m^3)
[]
[tritium_trapped_concentration] # (atoms/m^3)
[]
[helium_concentration] # (atoms/m^3)
[]
[]
[AuxKernels]
[empty_sites]
variable = empty_sites
type = EmptySitesAux
N = ${fparse density_material / density_scalar} # (-)
Ct0 = 'trapping_sites_fraction_function' # atomic fraction
trap_per_free = ${trap_per_free}
trapped_concentration_variables = tritium_trapped_concentration_scaled
[]
[scaled_empty_sites]
variable = scaled_empty_sites
type = NormalizationAux
normal_factor = ${density_scalar}
source_variable = empty_sites
[]
[trapped_sites]
variable = trapped_sites
type = NormalizationAux
normal_factor = ${density_scalar}
source_variable = tritium_trapped_concentration_scaled
[]
[total_sites]
variable = total_sites
type = ParsedAux
coupled_variables = 'trapped_sites scaled_empty_sites'
expression = 'trapped_sites + scaled_empty_sites'
[]
[tritium_mobile_concentration]
type = NormalizationAux
variable = tritium_mobile_concentration
normal_factor = ${density_scalar}
source_variable = tritium_mobile_concentration_scaled
[]
[tritium_trapped_concentration]
type = NormalizationAux
variable = tritium_trapped_concentration
normal_factor = ${density_scalar}
source_variable = tritium_trapped_concentration_scaled
[]
[helium_concentration]
type = NormalizationAux
variable = helium_concentration
normal_factor = ${density_scalar}
source_variable = helium_concentration_scaled
[]
[]
[Kernels]
# kernels for the tritium concentration equation
[time_tritium]
type = TimeDerivative
variable = tritium_mobile_concentration_scaled
[]
[diffusion]
type = MatDiffusion
variable = tritium_mobile_concentration_scaled
diffusivity = diffusivity
[]
[decay_tritium]
type = MatReaction
variable = tritium_mobile_concentration_scaled
reaction_rate = '${fparse -decay_rate_constant}'
[]
[coupled_time_tritium]
type = ScaledCoupledTimeDerivative
variable = tritium_mobile_concentration_scaled
v = tritium_trapped_concentration_scaled
factor = ${trap_per_free}
[]
# re-adding it to the equation of mobile tritium because it is accounted for in coupled_time_tritium, and needs to be removed
[decay_tritium_trapped]
type = MatReaction
variable = tritium_mobile_concentration_scaled
v = tritium_trapped_concentration_scaled
reaction_rate = '${fparse - decay_rate_constant * trap_per_free}'
[]
# kernels for the helium concentration equation
[time_helium]
type = TimeDerivative
variable = helium_concentration_scaled
[]
[decay_helium_mobile]
type = MatReaction
variable = helium_concentration_scaled
v = tritium_mobile_concentration_scaled
reaction_rate = '${fparse decay_rate_constant}'
[]
[decay_helium_trapped]
type = MatReaction
variable = helium_concentration_scaled
v = tritium_trapped_concentration_scaled
reaction_rate = '${fparse decay_rate_constant}'
[]
[]
[NodalKernels]
[time]
type = TimeDerivativeNodalKernel
variable = tritium_trapped_concentration_scaled # (atoms/m^3) / density_scalar = (-)
[]
[trapping]
type = TrappingNodalKernel
variable = tritium_trapped_concentration_scaled # (-)
alpha_t = ${tritium_trapping_prefactor} # (1/s)
trapping_energy = ${fparse tritium_trapping_energy/ideal_gas_constant} # (K)
N = ${fparse density_material / density_scalar} # (-)
Ct0 = 'trapping_sites_fraction_function' # (atomic fraction)
mobile_concentration = 'tritium_mobile_concentration_scaled' # (-)
temperature = temperature # (K)
trap_per_free = ${trap_per_free}
[]
[release]
type = ReleasingNodalKernel
variable = tritium_trapped_concentration_scaled # (atoms/m^3) / density_scalar = (-)
alpha_r = ${fparse tritium_release_prefactor} # (1/s)
detrapping_energy = ${fparse tritium_release_energy / boltzmann_constant} # (K)
temperature = temperature # (K)
[]
[decay]
type = ReleasingNodalKernel
variable = tritium_trapped_concentration_scaled # (atoms/m^3) / density_scalar = (-)
alpha_r = ${decay_rate_constant} # (1/s) - decay rate
detrapping_energy = 0 # (K) - The decay rate is independent of temperature
temperature = temperature # (K)
[]
[]
[Materials]
[diffusivity] # (m2/s) tritium diffusivity
type = DerivativeParsedMaterial
property_name = 'diffusivity'
coupled_variables = 'temperature'
expression = '${tritium_diffusivity_prefactor}*exp(-${tritium_diffusivity_energy}/${ideal_gas_constant}/temperature)'
[]
[alpha_t_tot] # (1/s) - trapping rate
type = ParsedMaterial
property_name = 'alpha_t_tot'
coupled_variables = 'temperature'
expression = '${tritium_trapping_prefactor} * exp(- ${tritium_trapping_energy}/${ideal_gas_constant}/temperature)'
outputs = 'all'
[]
[alpha_r_tot] # (1/s) - detrapping rate
type = ParsedMaterial
property_name = 'alpha_r_tot'
coupled_variables = 'temperature'
expression = '${tritium_release_prefactor} * exp(- ${tritium_release_energy} / ${boltzmann_constant}/temperature)'
outputs = 'all'
[]
[]
[Postprocessors]
# Amount of mobile tritium in the sample in (atoms/m^3) / density_scalar * m = (m)
[tritium_mobile_inventory_scaled]
type = ElementIntegralVariablePostprocessor
variable = tritium_mobile_concentration_scaled
execute_on = 'INITIAL TIMESTEP_END'
[]
[tritium_mobile_inventory] # (atoms/m^2)
type = ScalePostprocessor
value = tritium_mobile_inventory_scaled
scaling_factor = ${density_scalar}
execute_on = 'INITIAL TIMESTEP_END'
[]
# Amount of trapped tritium in the sample (atoms/m^3) / density_scalar * m = (m)
[tritium_trapped_inventory_scaled]
type = ElementIntegralVariablePostprocessor
variable = tritium_trapped_concentration_scaled
execute_on = 'INITIAL TIMESTEP_END'
[]
[tritium_trapped_inventory] # (atoms/m^2)
type = ScalePostprocessor
value = tritium_trapped_inventory_scaled
scaling_factor = ${density_scalar}
execute_on = 'INITIAL TIMESTEP_END'
[]
# Amount of helium in the sample in atoms/m^3 / density_scalar * m = (m)
[helium_inventory_scaled]
type = ElementIntegralVariablePostprocessor
variable = helium_concentration_scaled
execute_on = 'INITIAL TIMESTEP_END'
[]
[helium_inventory] # (atoms/m^2)
type = ScalePostprocessor
value = helium_inventory_scaled
scaling_factor = ${density_scalar}
execute_on = 'INITIAL TIMESTEP_END'
[]
# check mass conservation - this should remain constant - (atoms/m^2)
[total_inventory] # (atoms/m^2)
type = SumPostprocessor
values = 'tritium_mobile_inventory tritium_trapped_inventory helium_inventory'
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[VectorPostprocessors]
[line]
type = LineValueSampler
start_point = '0 0 0'
end_point = '${slab_length} 0 0'
num_points = ${num_mesh_element}
sort_by = 'x'
variable = 'tritium_mobile_concentration tritium_trapped_concentration helium_concentration'
execute_on = 'timestep_end'
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = ${end_time}
solve_type = NEWTON
scheme = 'bdf2'
dtmin = 1
dtmax = ${dt_max}
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt_start}
optimal_iterations = 9
growth_factor = 1.2
cutback_factor = 0.9
[]
[]
[Outputs]
perf_graph = true
file_base = ver-1jb_out
[time_dependent_out]
type = CSV
execute_vector_postprocessors_on = NONE
file_base = ver-1jb_time_dependent_out
[]
[profile_out]
type = CSV
sync_only = true
sync_times = '${units 45 year -> s}'
execute_postprocessors_on = NONE
file_base=ver-1jb_profile_out
[]
[]
(test/tests/val-2d/val-2d.i)
# Validation Problem #2d from TMAP7 V&V document
# Thermal Desorption Spectroscopy on Tungsten
# General parameters
kB = '${units 1.380649e-23 J/K}' # Boltzmann constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
# Model parameters
TDS_initial_time = '${units 5e3 s}'
TDS_critial_time_1 = '${units 5400 s}'
TDS_critial_time_2 = '${units 5404 s}'
simulation_time = '${units 6.8e3 s}'
outputs_initial_time = '${units 4000 s}'
step_interval_max = 50 # (-)
step_interval_mid = 15 # (-)
step_interval_min = 6 # (-)
bound_value_max = '${units 2e4 at/mum^3}'
bound_value_min = '${units -1e-10 at/mum^3}'
# Diffusion parameters
flux_high = '${units 1e19 at/m^2/s -> at/mum^2/s}'
flux_low = '${units 0 at/mum^2/s}'
diffusivity_coefficient = '${units 4.1e-7 m^2/s -> mum^2/s}'
E_D = '${units 0.39 eV -> J}'
initial_concentration = '${units 1e-10 at/m^3 -> at/mum^3}'
width_source = '${units 3e-9 m -> mum}'
depth_source = '${units 4.6e-9 m -> mum}'
# Traps parameters
N = '${units 6.25e28 at/m^3 -> at/mum^3}'
initial_concentration_trap_2 = 4.4e-10 # (-)
initial_concentration_trap_3 = 1.4e-10 # (-)
trapping_energy = '${fparse ${units 0.39 eV -> J} / ${kB}}'
detrapping_energy_1 = '${fparse ${units 1.2 eV -> J} / ${kB}}'
detrapping_energy_2 = '${fparse ${units 1.6 eV -> J} / ${kB}}'
detrapping_energy_3 = '${fparse ${units 3.1 eV -> J} / ${kB}}'
trapping_site_fraction_1 = 0.002156 # (-)
trapping_site_fraction_2 = 0.00175 # (-)
trapping_site_fraction_3 = 0.0020 # (-)
trapping_rate_prefactor = '${units 9.1316e12 1/s}'
release_rate_profactor = '${units 8.4e12 1/s}'
trap_per_free_1 = 1e6 # (-)
trap_per_free_2 = 1e4 # (-)
trap_per_free_3 = 1e4 # (-)
width_trap1 = '${units 10e-9 m -> mum}'
# thermal parameters
temperature_low = '${units 300 K}'
temperature_high = '${units 1273 K}'
temperature_rate = '${units ${fparse 50 / 60} K/s}'
[Mesh]
active = 'cartesian_mesh'
# finer mesh for heavy simulations
[cartesian_mesh]
nx_scale = 2
type = CartesianMeshGenerator
dim = 1
dx = '${fparse 10 * ${units 1.5e-9 m -> mum}}
${units 1e-9 m -> mum} ${units 1e-8 m -> mum} ${units 1e-7 m -> mum}
${units 4e-6 m -> mum} ${units 4.407e-6 m -> mum} ${fparse 11 * ${units 7.407e-6 m -> mum}}'
ix = '${fparse 10 * ${nx_scale}}
${fparse 1 * ${nx_scale}} ${fparse 1 * ${nx_scale}} ${fparse 1 * ${nx_scale}}
${fparse 50 * ${nx_scale}} ${fparse 2} ${fparse 1}'
subdomain_id = '0 1 1 1 1 1 1'
[]
# coarse mesh for tests
[cartesian_mesh_coarse]
nx_scale = 1
type = CartesianMeshGenerator
dim = 1
dx = '${fparse 10 * ${units 1.5e-9 m -> mum}}
${units 1e-9 m -> mum} ${units 1e-8 m -> mum} ${units 1e-7 m -> mum}
${units 4e-6 m -> mum} ${units 4.407e-6 m -> mum} ${fparse 11 * ${units 7.407e-6 m -> mum}}'
ix = '${fparse 10 * ${nx_scale}}
${fparse 1 * ${nx_scale}} ${fparse 1 * ${nx_scale}} ${fparse 1 * ${nx_scale}}
${fparse 6 * ${nx_scale}} ${fparse 1} ${fparse 1}'
subdomain_id = '0 1 1 1 1 1 1'
[]
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Bounds]
# Bounds to limit the maximum and minimum values for variables
[concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = concentration
bound_type = lower
bound_value = ${bound_value_min}
[]
[trapped_1_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_1
bound_type = lower
bound_value = ${bound_value_min}
[]
[trapped_2_upper_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_2
bound_type = upper
bound_value = ${bound_value_max}
[]
[trapped_2_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_2
bound_type = lower
bound_value = ${bound_value_min}
[]
[trapped_3_upper_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_3
bound_type = upper
bound_value = ${bound_value_max}
[]
[trapped_3_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_3
bound_type = lower
bound_value = ${bound_value_min}
[]
[]
[Variables]
# Mobile concentration variable
[concentration]
order = FIRST
family = LAGRANGE
initial_condition = ${initial_concentration}
[]
# Three trapped concentration variables
[trapped_1]
order = FIRST
family = LAGRANGE
block = 0
outputs = none
[]
[trapped_2]
order = FIRST
family = LAGRANGE
initial_condition = '${fparse initial_concentration_trap_2 * trapping_site_fraction_2 * N}'
outputs = none
[]
[trapped_3]
order = FIRST
family = LAGRANGE
initial_condition = '${fparse initial_concentration_trap_3 * trapping_site_fraction_3 * N}'
outputs = none
[]
[]
[AuxVariables]
[temperature]
[]
# necessary dummy variable for bounds
[bounds_dummy]
order = FIRST
family = LAGRANGE
[]
[]
[AuxKernels]
[temperature_Aux]
type = FunctionAux
variable = temperature
function = Temperature_function
execute_on = 'initial timestep_end linear'
[]
[]
[Kernels]
# mobile concentration kernel
[diffusion_implantation]
type = ADMatDiffusion
variable = concentration
diffusivity = Diffusivity
extra_vector_tags = ref
[]
[time_diffusion_implantation]
type = ADTimeDerivative
variable = concentration
extra_vector_tags = ref
[]
[source]
type = ADBodyForce
variable = concentration
function = concentration_source_norm_function
[]
# trapping kernel
[coupled_time_trap_1]
type = ADCoefCoupledTimeDerivative
variable = concentration
v = trapped_1
coef = ${trap_per_free_1}
block = 0
extra_vector_tags = ref
[]
[coupled_time_trap_2_implantation]
type = ADCoefCoupledTimeDerivative
variable = concentration
v = trapped_2
coef = ${trap_per_free_2}
extra_vector_tags = ref
[]
[coupled_time_trap_3_implantation]
type = ADCoefCoupledTimeDerivative
variable = concentration
v = trapped_3
coef = ${trap_per_free_3}
extra_vector_tags = ref
[]
[]
[NodalKernels]
# First traps
[time_1]
type = TimeDerivativeNodalKernel
variable = trapped_1
[]
[trapping_1]
type = TrappingNodalKernel
variable = trapped_1
mobile_concentration = concentration
alpha_t = '${trapping_rate_prefactor}'
trapping_energy = '${trapping_energy}'
N = '${N}'
Ct0 = 'trap_1_distribution_function'
temperature = 'temperature'
trap_per_free = ${trap_per_free_1}
extra_vector_tags = ref
[]
[release_1]
type = ReleasingNodalKernel
variable = trapped_1
alpha_r = '${release_rate_profactor}'
detrapping_energy = '${detrapping_energy_1}'
temperature = 'temperature'
[]
# Second traps
[time_2]
type = TimeDerivativeNodalKernel
variable = trapped_2
[]
[trapping_2_implantation]
type = TrappingNodalKernel
variable = trapped_2
mobile_concentration = concentration
alpha_t = '${trapping_rate_prefactor}'
trapping_energy = '${trapping_energy}'
N = '${N}'
Ct0 = '${trapping_site_fraction_2}'
temperature = 'temperature'
trap_per_free = ${trap_per_free_2}
extra_vector_tags = ref
[]
[release_2_implantation]
type = ReleasingNodalKernel
variable = trapped_2
alpha_r = '${release_rate_profactor}'
detrapping_energy = '${detrapping_energy_2}'
temperature = 'temperature'
[]
# Third traps
[time_3]
type = TimeDerivativeNodalKernel
variable = trapped_3
[]
[trapping_3_implantation]
type = TrappingNodalKernel
variable = trapped_3
mobile_concentration = concentration
alpha_t = '${trapping_rate_prefactor}'
trapping_energy = '${trapping_energy}'
N = '${N}'
Ct0 = '${trapping_site_fraction_3}'
temperature = 'temperature'
trap_per_free = ${trap_per_free_3}
extra_vector_tags = ref
[]
[release_3_implantation]
type = ReleasingNodalKernel
variable = trapped_3
alpha_r = '${release_rate_profactor}'
detrapping_energy = '${detrapping_energy_3}'
temperature = 'temperature'
[]
[]
[BCs]
# Constant concentration on both boundary
[left]
type = ADDirichletBC
variable = concentration
boundary = left
value = 0
[]
[right]
type = ADDirichletBC
variable = concentration
boundary = right
value = 0
[]
[]
[Materials]
[Diffusivity_implantation]
type = ADDerivativeParsedMaterial
property_name = 'Diffusivity'
functor_names = 'Temperature_function'
functor_symbols = 'Temperature_function'
expression = '${diffusivity_coefficient} * exp(- ${E_D} / ${kB} / Temperature_function)'
block = 0
output_properties = 'Diffusivity'
[]
[Diffusivity_Tungsten]
type = ADDerivativeParsedMaterial
property_name = 'Diffusivity'
functor_names = 'Temperature_function'
functor_symbols = 'Temperature_function'
expression = '${diffusivity_coefficient} * exp(- ${E_D} / ${kB} / Temperature_function) * 10'
block = 1
output_properties = 'Diffusivity'
[]
[converter_to_regular]
type = MaterialADConverter
ad_props_in = 'Diffusivity'
reg_props_out = 'Diffusivity_nonAD'
outputs = none
[]
[]
[Functions]
# Function for temperature and input flux
[Temperature_function]
type = ParsedFunction
expression = 'if(t<${TDS_initial_time}, ${temperature_low},
if(t<${TDS_initial_time} + (${temperature_high} - ${temperature_low}) /
${temperature_rate}, ${temperature_low} + ${temperature_rate} * (t - ${TDS_initial_time}),
${temperature_high}))'
[]
[surface_flux_function]
type = ParsedFunction
expression = 'if(t<${TDS_initial_time}, ${flux_high}, ${flux_low})'
[]
[source_distribution_function]
type = ParsedFunction
expression = '1 / ( ${width_source} * sqrt(2 * pi) ) * exp(-0.5 * ((x - ${depth_source}) / ${width_source} ) ^ 2)'
[]
[trap_1_distribution_function]
type = ParsedFunction
expression = ' ${trapping_site_fraction_1} / ( ${width_trap1} * sqrt(2 * pi) ) * exp(-0.5 * ((x - ${depth_source}) / ${width_trap1}) ^ 2)'
[]
[concentration_source_norm_function]
type = ParsedFunction
symbol_names = 'source_distribution_function surface_flux_function'
symbol_values = 'source_distribution_function surface_flux_function'
expression = 'source_distribution_function * surface_flux_function'
[]
# maximum time size during simulations
[max_dt_size_function]
type = ParsedFunction
expression = 'if(t<${TDS_initial_time} , ${step_interval_mid}, ${step_interval_min})'
[]
[max_dt_size_function_coarse]
type = ParsedFunction
expression = 'if(t<${TDS_initial_time} , ${step_interval_mid},
if(t<${TDS_critial_time_1}, ${step_interval_max},
if(t<${TDS_critial_time_2}, ${step_interval_min}, ${step_interval_max})))'
[]
[]
[Postprocessors]
# left flux
[flux_surface_left]
type = SideDiffusiveFluxIntegral
variable = concentration
diffusivity = 'Diffusivity_nonAD'
boundary = 'left'
outputs = none
[]
[scaled_flux_surface_left]
type = ScalePostprocessor
scaling_factor = '${units 1 m^2 -> mum^2}'
value = flux_surface_left
execute_on = 'initial nonlinear linear timestep_end'
outputs = 'console csv exodus'
[]
# right flux
[flux_surface_right]
type = SideDiffusiveFluxIntegral
variable = concentration
diffusivity = 'Diffusivity_nonAD'
boundary = 'right'
outputs = none
[]
[scaled_flux_surface_right]
type = ScalePostprocessor
scaling_factor = '${units 1 m^2 -> mum^2}'
value = flux_surface_right
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[max_time_step_size]
type = FunctionValuePostprocessor
function = max_dt_size_function
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = NEWTON
petsc_options_iname = '-pc_type -snes_type'
petsc_options_value = 'lu vinewtonrsls'
end_time = ${simulation_time}
line_search = 'none'
automatic_scaling = true
nl_rel_tol = 1e-10
nl_max_its = 34
[TimeStepper]
type = IterationAdaptiveDT
dt = 1.0
iteration_window = 5
optimal_iterations = 26
growth_factor = 1.1
cutback_factor = 0.9
cutback_factor_at_failure = 0.9
timestep_limiting_postprocessor = max_time_step_size
[]
[]
[Outputs]
file_base = 'val-2d_out'
[csv]
type = CSV
start_time = ${outputs_initial_time}
[]
[exodus]
type = Exodus
start_time = ${outputs_initial_time}
output_material_properties = true
time_step_interval = 20
[]
[]
(test/tests/ver-1d/ver-1d-trapping.i)
# Verification Problem #1d from TMAP4/TMAP7 V&V document
# Permeation Problem with Trapping in Trapping-limited Case
# No Soret effect, or solubility included.
# Modeling parameters
node_num = 1000
thickness = '${units 1 m}'
end_time = '${units 1000 s}'
time_scaling = 1
temperature = '${units 1000 K}'
diffusivity = '${fparse ${units 1 m^2/s} / time_scaling}'
# Trapping parameters
cl = '${units 3.1622e18 at/m^3}'
N = '${units 3.1622e22 at/m^3}'
trapping_prefactor = '${fparse ${units 1e15 1/s} / time_scaling}'
release_prefactor = '${fparse ${units 1e13 1/s} / time_scaling}'
epsilon=${units 10000 K}
trapping_fraction = 0.1 # -
trap_per_free = 1e3
[Mesh]
type = GeneratedMesh
dim = 1
nx = ${node_num}
xmax = ${thickness}
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[mobile]
[]
[trapped]
[]
[]
[AuxVariables]
[empty_sites]
[]
[scaled_empty_sites]
[]
[trapped_sites]
[]
[total_sites]
[]
[]
[AuxKernels]
[empty_sites]
variable = empty_sites
type = EmptySitesAux
N = '${fparse N / cl}'
Ct0 = ${trapping_fraction}
trap_per_free = ${trap_per_free}
trapped_concentration_variables = trapped
[]
[scaled_empty]
variable = scaled_empty_sites
type = NormalizationAux
normal_factor = ${cl}
source_variable = empty_sites
[]
[trapped_sites]
variable = trapped_sites
type = NormalizationAux
normal_factor = ${trap_per_free}
source_variable = trapped
[]
[total_sites]
variable = total_sites
type = ParsedAux
expression = 'trapped_sites + empty_sites'
coupled_variables = 'trapped_sites empty_sites'
[]
[]
[Kernels]
[diff]
type = MatDiffusion
variable = mobile
diffusivity = ${diffusivity}
extra_vector_tags = ref
[]
[time]
type = TimeDerivative
variable = mobile
extra_vector_tags = ref
[]
[coupled_time]
type = ScaledCoupledTimeDerivative
variable = mobile
v = trapped
factor = ${trap_per_free}
extra_vector_tags = ref
[]
[]
[NodalKernels]
[time]
type = TimeDerivativeNodalKernel
variable = trapped
[]
[trapping]
type = TrappingNodalKernel
variable = trapped
alpha_t = ${trapping_prefactor}
N = '${fparse N / cl}'
Ct0 = ${trapping_fraction}
mobile_concentration = 'mobile'
temperature = ${temperature}
trap_per_free = ${trap_per_free}
extra_vector_tags = ref
[]
[release]
type = ReleasingNodalKernel
alpha_r = ${release_prefactor}
temperature = ${temperature}
detrapping_energy = ${epsilon}
variable = trapped
[]
[]
[BCs]
[left]
type = FunctionDirichletBC
variable = mobile
function = 'BC_func'
boundary = left
[]
[right]
type = DirichletBC
variable = mobile
value = 0
boundary = right
[]
[]
[Functions]
[BC_func]
type = ParsedFunction
expression = '${fparse cl / cl}*tanh( 3 * t )'
[]
[]
[Postprocessors]
[outflux]
type = SideDiffusiveFluxAverage
boundary = 'right'
diffusivity = 1
variable = mobile
[]
[scaled_outflux]
type = ScalePostprocessor
value = outflux
scaling_factor = ${cl}
[]
[min_trapped]
type = NodalExtremeValue
value_type = MIN
variable = trapped
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = ${end_time}
dtmax = 5
solve_type = NEWTON
scheme = BDF2
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
line_search = 'none'
[TimeStepper]
type = IterationAdaptiveDT
dt = 1e-6
optimal_iterations = 9
growth_factor = 1.1
cutback_factor = 0.909
[]
[]
[Outputs]
exodus = true
csv = true
[dof]
type = DOFMap
execute_on = initial
[]
perf_graph = true
[]