diff --git a/doc/content/source/kernels/ReleasingKernel.md b/doc/content/source/kernels/ReleasingKernel.md new file mode 100644 index 000000000..0b858e188 --- /dev/null +++ b/doc/content/source/kernels/ReleasingKernel.md @@ -0,0 +1,12 @@ +# ReleasingKernel + +!syntax description /Kernels/ReleasingKernel + +This kernel is the volumetric equivalent of the [ReleasingNodalKernel.md]. See this nodal +kernel documentation for information on the releasing mechanism. + +!syntax parameters /Kernels/ReleasingKernel + +!syntax inputs /Kernels/ReleasingKernel + +!syntax children /Kernels/ReleasingKernel diff --git a/doc/content/source/kernels/TrappingKernel.md b/doc/content/source/kernels/TrappingKernel.md new file mode 100644 index 000000000..be38529b9 --- /dev/null +++ b/doc/content/source/kernels/TrappingKernel.md @@ -0,0 +1,12 @@ +# TrappingKernel + +!syntax description /Kernels/TrappingKernel + +This kernel is the volumetric equivalent of the [TrappingNodalKernel.md]. See this nodal +kernel documentation for information on the trapping mechanism. + +!syntax parameters /Kernels/TrappingKernel + +!syntax inputs /Kernels/TrappingKernel + +!syntax children /Kernels/TrappingKernel diff --git a/include/kernels/ReleasingKernel.h b/include/kernels/ReleasingKernel.h new file mode 100644 index 000000000..1c38eae77 --- /dev/null +++ b/include/kernels/ReleasingKernel.h @@ -0,0 +1,41 @@ +/************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* TMAP8: Tritium Migration Analysis Program, Version 8 */ +/* */ +/* Copyright 2021 - 2024 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/************************************************************/ + +#pragma once + +#include "Kernel.h" + +/** + * Implements the contribution to the residual and Jacobian of a species release using a volumetric + * integration + */ +class ReleasingKernel : public Kernel +{ +public: + ReleasingKernel(const InputParameters & parameters); + + static InputParameters validParams(); + +protected: + Real computeQpResidual() override; + Real computeQpJacobian() override; + Real computeQpOffDiagJacobian(unsigned int jvar) override; + + /// Release coefficient + const Real _alpha_r; + /// Energy from detrapping + const Real _detrapping_energy; + /// Local temperature + const VariableValue & _temperature; + /// Species concentration + const VariableValue & _v; + /// Index of the species variable + const unsigned int _v_index; + /// Whether the v variable is the kernel's variable parameter variable + const bool _v_is_u; +}; diff --git a/include/kernels/TrappingKernel.h b/include/kernels/TrappingKernel.h new file mode 100644 index 000000000..94c67b6e9 --- /dev/null +++ b/include/kernels/TrappingKernel.h @@ -0,0 +1,44 @@ +/************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* TMAP8: Tritium Migration Analysis Program, Version 8 */ +/* */ +/* Copyright 2021 - 2024 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/************************************************************/ + +#pragma once + +#include "ADKernel.h" + +class Function; + +/** + * Implements the contribution to the residual and Jacobian using a volumetric integration and + * automatic differentiation + */ +class TrappingKernel : public ADKernel +{ +public: + TrappingKernel(const InputParameters & parameters); + + static InputParameters validParams(); + +protected: + ADReal computeQpResidual() override; + + /// Trapping coefficient + const Real _alpha_t; + /// Trapping energy + const Real _trapping_energy; + const Real _N; + const Function & _Ct0; + /// Concentration of the mobile species + const ADVariableValue & _mobile_concentration; + /// Number of other species competing for the same traps + unsigned int _n_other_concs; + /// Concentration of the other species + std::vector _trapped_concentrations; + const Real _trap_per_free; + /// Local temperature + const ADVariableValue & _temperature; +}; diff --git a/include/nodal_kernels/ReleasingNodalKernel.h b/include/nodal_kernels/ReleasingNodalKernel.h index 56eb13b56..53347c6cc 100644 --- a/include/nodal_kernels/ReleasingNodalKernel.h +++ b/include/nodal_kernels/ReleasingNodalKernel.h @@ -20,8 +20,24 @@ class ReleasingNodalKernel : public NodalKernel protected: Real computeQpResidual() override; Real computeQpJacobian() override; + Real computeQpOffDiagJacobian(unsigned int jvar) override; + /// Release coefficient const Real _alpha_r; + /// Energy from detrapping const Real _detrapping_energy; + /// Local temperature const VariableValue & _temperature; + /// Species concentration + const VariableValue & _v; + /// Index of the species variable + const unsigned int _v_index; + /// Whether the v variable is the kernel's variable parameter variable + const bool _v_is_u; + /// Whether the kernels are mass lumped to make it compatible + const bool _mass_lumped; + /// Local node mass + const VariableValue & _nodal_mass; + /// An array with ones for convenience + const VariableValue _one; }; diff --git a/include/nodal_kernels/TrappingNodalKernel.h b/include/nodal_kernels/TrappingNodalKernel.h index 9e1daec59..10eb0b0a7 100644 --- a/include/nodal_kernels/TrappingNodalKernel.h +++ b/include/nodal_kernels/TrappingNodalKernel.h @@ -44,6 +44,12 @@ class TrappingNodalKernel : public NodalKernel const Real _trap_per_free; const VariableValue & _temperature; LocalDN _jacobian; + /// Whether the kernels are mass lumped to make it compatible + const bool _mass_lumped; + /// Local node mass + const VariableValue & _nodal_mass; + /// An array with ones for convenience + const VariableValue _one; private: void ADHelper(); diff --git a/include/physics/SpeciesRadioactiveDecay.h b/include/physics/SpeciesRadioactiveDecay.h new file mode 100644 index 000000000..4409f11fc --- /dev/null +++ b/include/physics/SpeciesRadioactiveDecay.h @@ -0,0 +1,51 @@ +/********************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* TMAP8: Tritium Migration Analysis Program, Version 8 */ +/* */ +/* Copyright 2021 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/********************************************************/ + +#pragma once + +#include "SpeciesPhysicsBase.h" + +class ActionComponent; + +/** + * Creates all the objects needed to solve for the radioactive decay of the species + */ +class SpeciesRadioactiveDecay : public SpeciesPhysicsBase +{ +public: + static InputParameters validParams(); + + SpeciesRadioactiveDecay(const InputParameters & parameters); + + void addComponent(const ActionComponent & component) override; + +protected: + /// Return the name of the species variable + /// @param c_i index of the component + /// @param s_j index of the species + VariableName getSpeciesVariableName(unsigned int c_i, unsigned int s_j) const; + + /// Decay products on each component, for each species, for each decay reaction + std::vector>>> _decay_products; + /// Decay constants on each component, for each species, for each decay reaction + std::vector>> _decay_constants; + /// Branching rations on each component, for each species, for each decay reaction, for each product species + std::vector>> _branching_ratios; + + /// Whether to define a single variable for each species for all components, or a different one for each component + const bool _single_variable_set; + /// Whether to define initial conditions for the decaying species + const bool _add_initial_conditions; + /// Whether to use a nodal or a volumetric formulation for the definition of kernels, indexed per component + std::vector _use_nodal_formulations; + +private: + virtual void addSolverVariables() override; + virtual void addInitialConditions() override; + virtual void addFEKernels() override; +}; diff --git a/include/physics/SpeciesTrappingPhysics.h b/include/physics/SpeciesTrappingPhysics.h index 7e6096712..df99d57ae 100644 --- a/include/physics/SpeciesTrappingPhysics.h +++ b/include/physics/SpeciesTrappingPhysics.h @@ -31,6 +31,8 @@ class SpeciesTrappingPhysics : public SpeciesPhysicsBase /// @param s_j index of the species VariableName getSpeciesVariableName(unsigned int c_i, unsigned int s_j) const; + /// The discretization method for the trapping equations + const MooseEnum _discretization; /// The mobile species of interest std::vector> _mobile_species_names; @@ -54,6 +56,7 @@ class SpeciesTrappingPhysics : public SpeciesPhysicsBase const bool _single_variable_set; private: + bool isNodal() const { return _discretization == "nodal"; } virtual void addSolverVariables() override; virtual void addInitialConditions() override; virtual void addFEKernels() override; diff --git a/moose b/moose index 201ada841..0f5a59842 160000 --- a/moose +++ b/moose @@ -1 +1 @@ -Subproject commit 201ada841f9ac9d9ed52bc1a21536739a77d6abb +Subproject commit 0f5a59842032ebabbe3bd418857801820caed136 diff --git a/src/base/TMAP8App.C b/src/base/TMAP8App.C index a9352f027..66c4db3c7 100644 --- a/src/base/TMAP8App.C +++ b/src/base/TMAP8App.C @@ -42,6 +42,8 @@ TMAP8App::registerAll(Factory & f, ActionFactory & af, Syntax & syntax) // TMAP8 specific Physics registerSyntax("SorptionExchangePhysics", "Physics/SorptionExchange/*"); registerSyntax("SpeciesTrappingPhysics", "Physics/SpeciesTrapping/*"); + registerSyntax("SpeciesDiffusionReactionCG", "Physics/SpeciesDiffusionReaction/*"); + registerSyntax("SpeciesRadioactiveDecay", "Physics/SpeciesRadioactiveDecay/*"); // Shorter syntax for MOOSE Physics used by TMAP8 registerSyntax("DiffusionCG", "Physics/Diffusion/*"); diff --git a/src/kernels/ReleasingKernel.C b/src/kernels/ReleasingKernel.C new file mode 100644 index 000000000..1c1f671c1 --- /dev/null +++ b/src/kernels/ReleasingKernel.C @@ -0,0 +1,69 @@ +/************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* TMAP8: Tritium Migration Analysis Program, Version 8 */ +/* */ +/* Copyright 2021 - 2024 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/************************************************************/ + +#include "ReleasingKernel.h" + +registerMooseObject("TMAP8App", ReleasingKernel); + +InputParameters +ReleasingKernel::validParams() +{ + InputParameters params = Kernel::validParams(); + params.addClassDescription( + "Implements a residual describing the release of trapped species in a material."); + params.addRequiredParam("alpha_r", "The release rate coefficient (1/s)"); + params.addCoupledVar("trapped_concentration", + "If specified, variable to compute the release rate with. " + "Else, uses the 'variable' argument"); + params.addParam("detrapping_energy", 0, "The detrapping energy (K)"); + params.addRequiredCoupledVar("temperature", "The temperature (K)"); + return params; +} + +ReleasingKernel::ReleasingKernel(const InputParameters & parameters) + : Kernel(parameters), + _alpha_r(getParam("alpha_r")), + _detrapping_energy(getParam("detrapping_energy")), + _temperature(coupledValue("temperature")), + _v(isParamValid("trapped_concentration") ? coupledValue("trapped_concentration") : _u), + _v_index(coupled("trapped_concentration")), + _v_is_u(!isCoupled("trapped_concentration") || + (coupled("trapped_concentration") == variable().number())) +{ +} + +Real +ReleasingKernel::computeQpResidual() +{ + // std::cout << "Volumetric " + // << _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _v[_qp] << + // std::endl; + return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _v[_qp] * _test[_i][_qp]; +} + +Real +ReleasingKernel::computeQpJacobian() +{ + if (_v_is_u) + return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _phi[_j][_qp] * + _test[_i][_qp]; + else + return 0; +} + +Real +ReleasingKernel::computeQpOffDiagJacobian(unsigned int jvar) +{ + if (!_v_is_u && jvar == _v_index) + return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _phi[_j][_qp] * + _test[_i][_qp]; + else + return 0; + + // TODO: add temperature off-diagonal term +} diff --git a/src/kernels/TrappingKernel.C b/src/kernels/TrappingKernel.C new file mode 100644 index 000000000..d4409dd7c --- /dev/null +++ b/src/kernels/TrappingKernel.C @@ -0,0 +1,77 @@ +/************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* TMAP8: Tritium Migration Analysis Program, Version 8 */ +/* */ +/* Copyright 2021 - 2024 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/************************************************************/ + +#include "TrappingKernel.h" +#include "Function.h" + +registerMooseObject("TMAP8App", TrappingKernel); + +InputParameters +TrappingKernel::validParams() +{ + InputParameters params = ADKernel::validParams(); + params.addClassDescription( + "Implements a residual describing the trapping of a species in a material."); + params.addRequiredParam("alpha_t", + "The trapping rate coefficient. This has units of 1/s (e.g. no " + "number densities are involved)"); + params.addParam("trapping_energy", 0, "The trapping energy (K)"); + params.addRequiredParam("N", "The atomic number density of the host material (1/m^3)"); + params.addRequiredParam( + "Ct0", "The fraction of host sites that can contribute to trapping as a function (-)"); + params.addParam( + "trap_per_free", + 1., + "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 (-)"); + params.addRequiredCoupledVar( + "mobile_concentration", + "The variable representing the mobile concentration of solute particles (1/m^3)"); + params.addCoupledVar("trapped_concentration", + "The variable representing the trapped concentration of solute particles " + "(1/m^3). If unspecified, defaults to the 'variable' parameter"); + params.addCoupledVar("other_trapped_concentration_variables", + "Other variables representing trapped particle concentrations."); + params.addRequiredCoupledVar("temperature", "The temperature (K)"); + return params; +} + +TrappingKernel::TrappingKernel(const InputParameters & parameters) + : ADKernel(parameters), + _alpha_t(getParam("alpha_t")), + _trapping_energy(getParam("trapping_energy")), + _N(getParam("N")), + _Ct0(getFunction("Ct0")), + _mobile_concentration(adCoupledValue("mobile_concentration")), + _trap_per_free(getParam("trap_per_free")), + _temperature(adCoupledValue("temperature")) +{ + _n_other_concs = coupledComponents("other_trapped_concentration_variables"); + + // Resize to n_other_concs plus the concentration corresponding to this Kernel's variable + _trapped_concentrations.resize(1 + _n_other_concs); + + // Keep pointers to references of the trapped species concentrations + for (MooseIndex(_n_other_concs) i = 0; i < _n_other_concs; ++i) + _trapped_concentrations[i] = + &adCoupledValue("other_trapped_concentration_variables", /*comp*/ 0); + _trapped_concentrations[_n_other_concs] = + isParamValid("trapped_concentration") ? &adCoupledValue("trapped_concentration") : &_u; +} + +ADReal +TrappingKernel::computeQpResidual() +{ + // Remove filled trapping sites from total number of trapping sites + ADReal empty_trapping_sites = _Ct0.value(_t, _q_point[_qp]) * _N; + for (const auto & trap_conc : _trapped_concentrations) + empty_trapping_sites -= (*trap_conc)[_qp] * _trap_per_free; + + return -_alpha_t * std::exp(-_trapping_energy / _temperature[_qp]) * empty_trapping_sites * + _mobile_concentration[_qp] / (_N * _trap_per_free) * _test[_i][_qp]; +} diff --git a/src/nodal_kernels/ReleasingNodalKernel.C b/src/nodal_kernels/ReleasingNodalKernel.C index faf123253..bba579345 100644 --- a/src/nodal_kernels/ReleasingNodalKernel.C +++ b/src/nodal_kernels/ReleasingNodalKernel.C @@ -17,8 +17,20 @@ ReleasingNodalKernel::validParams() params.addClassDescription( "Implements a residual describing the release of trapped species in a material."); params.addRequiredParam("alpha_r", "The release rate coefficient (1/s)"); + params.addCoupledVar("v", + "If specified, variable to compute the release rate with. " + "Else, uses the 'variable' argument"); + params.deprecateCoupledVar("v", "trapped_concentration", "1/1/2026"); params.addParam("detrapping_energy", 0, "The detrapping energy (K)"); params.addRequiredCoupledVar("temperature", "The temperature (K)"); + + // Optional mass lumping + params.addParam( + "use_mass_lumping", + false, + "Whether to use mass lumping to make this kernel compatible with volumetric kernels"); + params.addCoupledVar("nodal_mass", "The local nodal mass"); + params.addParamNamesToGroup("use_mass_lumping nodal_mass", "Mass lumping"); return params; } @@ -26,18 +38,42 @@ ReleasingNodalKernel::ReleasingNodalKernel(const InputParameters & parameters) : NodalKernel(parameters), _alpha_r(getParam("alpha_r")), _detrapping_energy(getParam("detrapping_energy")), - _temperature(coupledValue("temperature")) + _temperature(coupledValue("temperature")), + _v(isParamValid("trapped_concentration") ? coupledValue("trapped_concentration") : _u), + _v_index(coupled("trapped_concentration")), + _v_is_u(!isCoupled("trapped_concentration") || (_v_index == variable().number())), + _mass_lumped(getParam("use_mass_lumping")), + _nodal_mass(_mass_lumped ? coupledValue("nodal_mass") : _one), + _one(1) { } Real ReleasingNodalKernel::computeQpResidual() { - return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _u[_qp]; + const auto mass = _mass_lumped ? _nodal_mass[_qp] : 1.; + // std::cout << "Nodal " << _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _v[_qp] + // << std::endl; + return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _v[_qp] * mass; } Real ReleasingNodalKernel::computeQpJacobian() { - return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]); + const auto mass = _mass_lumped ? _nodal_mass[_qp] : 1.; + if (_v_is_u) + return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * mass; + else + return 0; +} + +Real +ReleasingNodalKernel::computeQpOffDiagJacobian(unsigned int jvar) +{ + const auto mass = _mass_lumped ? _nodal_mass[_qp] : 1.; + if (!_v_is_u && jvar == _v_index) + return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * mass; + else + return 0; + // TODO: add temperature off-diagonal term } diff --git a/src/nodal_kernels/TrappingNodalKernel.C b/src/nodal_kernels/TrappingNodalKernel.C index 0823fec2c..4f29a3915 100644 --- a/src/nodal_kernels/TrappingNodalKernel.C +++ b/src/nodal_kernels/TrappingNodalKernel.C @@ -32,9 +32,21 @@ TrappingNodalKernel::validParams() params.addRequiredCoupledVar( "mobile_concentration", "The variable representing the mobile concentration of solute particles (1/m^3)"); + params.addCoupledVar("trapped_concentration", + "The variable representing the trapped concentration of solute particles " + "(1/m^3). If unspecified, defaults to the 'variable' parameter"); params.addCoupledVar("other_trapped_concentration_variables", "Other variables representing trapped particle concentrations."); params.addRequiredCoupledVar("temperature", "The temperature (K)"); + + // Mass lumping + params.addParam( + "use_mass_lumping", + false, + "Whether to use mass lumping to make this kernel compatible with volumetric kernels"); + params.addCoupledVar("nodal_mass", "The local nodal mass"); + params.addParamNamesToGroup("use_mass_lumping nodal_mass", "Mass lumping"); + return params; } @@ -47,7 +59,10 @@ TrappingNodalKernel::TrappingNodalKernel(const InputParameters & parameters) _mobile_concentration(coupledValue("mobile_concentration")), _last_node(nullptr), _trap_per_free(getParam("trap_per_free")), - _temperature(coupledValue("temperature")) + _temperature(coupledValue("temperature")), + _mass_lumped(getParam("use_mass_lumping")), + _nodal_mass(_mass_lumped ? coupledValue("nodal_mass") : _one), + _one(1) { _n_other_concs = coupledComponents("other_trapped_concentration_variables"); @@ -62,8 +77,10 @@ TrappingNodalKernel::TrappingNodalKernel(const InputParameters & parameters) _trapped_concentrations[i] = &coupledValue("other_trapped_concentration_variables", /*comp=*/i); _var_numbers[i] = coupled("other_trapped_concentration_variables", i); } - _trapped_concentrations[_n_other_concs] = &_u; - _var_numbers[_n_other_concs] = _var.number(); + _trapped_concentrations[_n_other_concs] = + isCoupled("trapped_concentration") ? &coupledValue("trapped_concentration") : &_u; + _var_numbers[_n_other_concs] = + isParamValid("trapped_concentration") ? coupled("trapped_concentration") : _var.number(); _var_numbers[_n_other_concs + 1] = coupled("mobile_concentration"); } @@ -73,9 +90,10 @@ TrappingNodalKernel::computeQpResidual() auto empty_trapping_sites = _Ct0.value(_t, (*_current_node)) * _N; for (const auto & trap_conc : _trapped_concentrations) empty_trapping_sites -= (*trap_conc)[_qp] * _trap_per_free; + const auto mass = _mass_lumped ? _nodal_mass[_qp] : 1; return -_alpha_t * std::exp(-_trapping_energy / _temperature[_qp]) * empty_trapping_sites * - _mobile_concentration[_qp] / (_N * _trap_per_free); + _mobile_concentration[_qp] / (_N * _trap_per_free) * mass; } void @@ -98,9 +116,10 @@ TrappingNodalKernel::ADHelper() LocalDN mobile_concentration = _mobile_concentration[_qp]; mobile_concentration.derivatives().insert(_var_numbers.back()) = 1.; + const auto mass = _mass_lumped ? _nodal_mass[_qp] : 1; _jacobian = -_alpha_t * std::exp(-_trapping_energy / _temperature[_qp]) * empty_trapping_sites * - mobile_concentration / (_N * _trap_per_free); + mobile_concentration / (_N * _trap_per_free) * mass; } Real diff --git a/src/physics/SpeciesRadioactiveDecay.C b/src/physics/SpeciesRadioactiveDecay.C new file mode 100644 index 000000000..c6347eafe --- /dev/null +++ b/src/physics/SpeciesRadioactiveDecay.C @@ -0,0 +1,363 @@ +/********************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* TMAP8: Tritium Migration Analysis Program, Version 8 */ +/* */ +/* Copyright 2021 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/********************************************************/ + +#include "SpeciesRadioactiveDecay.h" +#include "ActionComponent.h" +#include "MooseUtils.h" +#include "FEProblemBase.h" + +// Component interaction +#include "Enclosure0D.h" + +// Register the actions for the objects actually used +registerMooseAction("TMAP8App", SpeciesRadioactiveDecay, "init_physics"); +registerMooseAction("TMAP8App", SpeciesRadioactiveDecay, "init_component_physics"); +registerMooseAction("TMAP8App", SpeciesRadioactiveDecay, "copy_vars_physics"); +registerMooseAction("TMAP8App", SpeciesRadioactiveDecay, "check_integrity"); +registerMooseAction("TMAP8App", SpeciesRadioactiveDecay, "check_integrity_early_physics"); +registerMooseAction("TMAP8App", SpeciesRadioactiveDecay, "add_variable"); +registerMooseAction("TMAP8App", SpeciesRadioactiveDecay, "add_ic"); +registerMooseAction("TMAP8App", SpeciesRadioactiveDecay, "add_kernel"); + +InputParameters +SpeciesRadioactiveDecay::validParams() +{ + InputParameters params = SpeciesPhysicsBase::validParams(); + params.addClassDescription("Add Physics for the radioactive decay of species."); + + params.addParam("separate_variables_per_component", + false, + "Whether to create new variables for each trapped species on every " + "component, or whether to only create variables."); + + // Manual input of decay products, constants and branching ratios + // We could consider just loading a file from ENDF! + params.addParam>>>( + "decay_products", + {}, + "Decay products. Use '|' to separate inputs for each decaying species. Use ';' to separate " + "each decay reaction for a given species. Use 'NA' " + "if tracking the product is not desired. If specified in the Physics, applies to every " + "component."); + params.addParam>>( + "decay_constants", + {}, + "The decay constants for each decay reaction. Use ';' to separate input for each decaying " + "species. If specified in the Physics, applies to every component."); + params.addParam>>( + "branching_ratios", + {}, + "Branching ratios. Use ';' to separate inputs for each decay reaction within each decaying " + "species. Defaults to 1 if not specified for a given species " + "or group. If specified in the Physics, applies to every component."); + + params.addRequiredParam( + "add_decaying_species_initial_conditions", + "Whether to set initial conditions from this Physics for the decaying species. This should " + "be set to false if another Physics is already defining the initial conditions"); + + return params; +} + +SpeciesRadioactiveDecay::SpeciesRadioactiveDecay(const InputParameters & parameters) + : SpeciesPhysicsBase(parameters), + // If specified in the Physics block, all parameters are retrieved here + _decay_products( + {getParam>>>("decay_products")}), + _decay_constants({getParam>>("decay_constants")}), + _branching_ratios({getParam>>("branching_ratios")}), + _single_variable_set(!getParam("separate_variables_per_component")), + _add_initial_conditions(getParam("add_decaying_species_initial_conditions")) +{ + // All the other parameters can vary on each component + if (_single_variable_set) + checkVectorParamNotEmpty("species"); + if (!_add_initial_conditions) + errorDependentParameter("add_initial_condition", "true", {"species_initial_concentrations"}); + + // Only set the other parameters if setting the species + checkSecondParamSetOnlyIfFirstOneSet("species", "decay_products"); + checkSecondParamSetOnlyIfFirstOneSet("species", "decay_constants"); + checkSecondParamSetOnlyIfFirstOneSet("species", "branching_ratios"); + + // Check sizes + checkVectorParamsSameLengthIfSet( + "species", "species_initial_concentrations", /*ignore_empty_second*/ true); + checkVectorParamsSameLengthIfSet>( + "species", "decay_constants", true); + // TODO: add a 3D vector check to the input parameter check util + if (isParamSetByUser("branching_ratios")) + checkTwoDVectorParamsSameLength, Real>("decay_products", + "branching_ratios"); + checkTwoDVectorParamsSameLength, Real>("decay_products", + "decay_constants"); +} + +void +SpeciesRadioactiveDecay::addComponent(const ActionComponent & component) +{ + for (const auto & block : component.blocks()) + _blocks.push_back(block); + _components.push_back(component.name()); + + // Index of the component in all the component-indexed vectors + const auto comp_index = _components.size() - 1; + + // Process each of the component's parameters, adding defaults to avoid breaking the double-vector + // indexing when acceptable + // These parameters are known to be defined for a Structure1D, so we retrieve them from the + // component's parameters. If they are not defined on the Physics or the component, we error + processComponentParameters>( + "species", component.name(), comp_index, _species, "species", false, {}); + processComponentParameters>("species_scaling_factors", + component.name(), + comp_index, + _scaling_factors, + "species_scaling_factors", + true, + std::vector(_species[comp_index].size(), 1)); + processComponentParameters>("species_initial_concentrations", + component.name(), + comp_index, + _initial_conditions, + "species_initial_concentrations", + false, + {}); + + // We dont support Enclosures yet + if (dynamic_cast(&component)) + mooseError("This Physics has not been implemented for 0D enclosures"); + mooseWarning( + "Processing parameters from Components has not been implemented yet for this Physics"); + + // These parameters should be defined as material properties by the user on the Component + // or on the Physics. + // We only support Real numbers for now as the consuming kernels only support Real + // TODO: decay constants, branching ratios and decay products + // processComponentMatprop>( + // "decay_constants", component, comp_index, _species.back(), _decay_constants); +} + +VariableName +SpeciesRadioactiveDecay::getSpeciesVariableName(unsigned int c_i, unsigned int s_j) const +{ + mooseAssert(c_i < _species.size(), "component index higher than number of components"); + mooseAssert(s_j < _species[c_i].size(), "species index higher than number of species"); + if (_single_variable_set) + return _species[0][s_j]; + else + // Add the component name if defining variables on a per-component basis + return _species[c_i][s_j] + "_" + _components[c_i]; +} + +void +SpeciesRadioactiveDecay::addSolverVariables() +{ + // TODO: this should be a scalar variable for a 0D component + + const std::string variable_type = "MooseVariable"; + InputParameters params = getFactory().getValidParams(variable_type); + params.set("family") = "LAGRANGE"; + params.set("order") = FIRST; + + // Allow using blocks even with the loops on components + if (_components.empty()) + { + if (_single_variable_set) + { + _components.push_back(""); + if (_species[0].empty()) + paramError("species", "Should not be empty if not using Components"); + if (_scaling_factors[0].empty()) + _scaling_factors[0] = std::vector(_species.size(), 1); + } + else + paramError("separate_variables_per_component", + "Physics is not defined on any Component, this parameter should be set to false"); + } + else + { + // Check component-indexed parameters + checkSizeComponentSpeciesIndexedVectorOfVector( + _scaling_factors, "species_scaling_factors", true); + checkSizeComponentSpeciesIndexedVectorOfVector(_decay_products, "decay_products", true); + } + + for (const auto c_i : index_range(_components)) + { + // Use the whole phyiscs block restriction if using the same species variable everywhere + if (_single_variable_set) + assignBlocks(params, _blocks); + else + assignBlocks(params, getActionComponent(_components[c_i]).blocks()); + + for (const auto s_j : index_range(_species[c_i])) + { + const auto species_name = getSpeciesVariableName(c_i, s_j); + if (isParamSetByUser("species_scaling_factor") || !_single_variable_set) + params.set>("scaling") = { + (_scaling_factors.size() > 1) + ? _scaling_factors[c_i][s_j] + : ((_scaling_factors.size() == 1) ? _scaling_factors[0][s_j] : 1)}; + params.set("solver_sys") = getSolverSystem(species_name); + + // TODO: add a check in PhysicsBase for avoiding defining variables with different scaling + // factors from multiple Physics Species is already likely added by another Physics + if (!getProblem().hasVariable(species_name)) + getProblem().addVariable(variable_type, species_name, params); + + // Keep track of variables + saveSolverVariableName(species_name); + } + if (_single_variable_set) + break; + } +} + +void +SpeciesRadioactiveDecay::addInitialConditions() +{ + if (!_add_initial_conditions) + return; + + const std::string ic_type = "FunctorIC"; + InputParameters params = getFactory().getValidParams(ic_type); + + // Check component-indexed parameters + if (_components.size()) + checkSizeComponentSpeciesIndexedVectorOfVector( + _initial_conditions, "species_initial_concentrations", true); + + for (const auto c_i : index_range(_components)) + { + // Use the whole phyiscs block restriction if using the same species variable everywhere + if (_single_variable_set) + if (isParamSetByUser("species_initial_concentrations")) + assignBlocks(params, _blocks); + else + break; + else + assignBlocks(params, getActionComponent(_components[c_i]).blocks()); + + // Decaying species + for (const auto s_j : index_range(_species[c_i])) + { + const auto species_name = getSpeciesVariableName(c_i, s_j); + params.set("variable") = species_name; + params.set("functor") = + ((_initial_conditions.size() > 1) + ? _initial_conditions[c_i][s_j] + : ((_initial_conditions.size() == 1) ? _initial_conditions[0][s_j] : "0")); + + // TODO: add a check in PhysicsBase for avoiding defining variables with different initial + // conditions from multiple Physics. The species IC is already likely added by another Physics + if (_add_initial_conditions) + getProblem().addInitialCondition( + ic_type, "IC_" + species_name + "_" + Moose::stringify(_blocks), params); + } + // TODO: create product species if they dont exist + + if (_single_variable_set) + break; + } +} + +void +SpeciesRadioactiveDecay::addFEKernels() +{ + // Prefill branching ratios with 1s if it is not specified + if (_single_variable_set) + if (!isParamSetByUser("branching_ratios")) + { + std::cout << "Resizing " << std::endl; + _branching_ratios.resize(_decay_products.size()); + for (const auto i_c : index_range(_decay_products)) + { + std::cout << "Resizing for species " << _decay_products[i_c].size() << std::endl; + _branching_ratios[i_c].resize(_decay_products[i_c].size()); + for (const auto i_s : index_range(_decay_products[i_c])) + { + _branching_ratios[i_c][i_s].resize(_decay_products[i_c][i_s].size(), 1.); + } + } + } + + // Check component-indexed parameters + checkSizeComponentSpeciesIndexedVectorOfVector(_decay_products, "decay_products", true); + checkSizeComponentSpeciesIndexedVectorOfVector(_decay_constants, "decay_constants", true); + checkSizeComponentSpeciesIndexedVectorOfVector(_branching_ratios, "branching_ratios", true); + + for (const auto c_i : index_range(_components)) + { + // Use the whole phyiscs block restriction if using the same species variable everywhere + const auto blocks = + _single_variable_set ? _blocks : getActionComponent(_components[c_i]).blocks(); + const auto & comp_name = _components[c_i]; + + // Create the kernel for each species + for (const auto s_j : index_range(_species[c_i])) + { + const auto species_name = getSpeciesVariableName(c_i, s_j); + + // Time derivative + if (isTransient()) + { + // TODO: add volumetric discretization option + const std::string kernel_type = "TimeDerivativeNodalKernel"; + InputParameters params = getFactory().getValidParams(kernel_type); + params.set("variable") = species_name; + assignBlocks(params, blocks); + // TODO: add a check in PhysicsBase for avoiding defining time derivatives for variable from + // multiple Physics on the same blocks, as this term is already likely added by another + // Physics + getProblem().addNodalKernel( + kernel_type, prefix() + comp_name + "_" + species_name + "_time", params); + } + + // Decay reactions for each species + for (const auto & r_k : index_range(_decay_products[c_i][s_j])) + { + /* Decay of species */ + // TODO: add volumetric discretization option + const std::string kernel_type = "ReactionNodalKernel"; + auto params = _factory.getValidParams(kernel_type); + assignBlocks(params, blocks); + params.set("variable") = species_name; + params.set("coeff") = _decay_constants[c_i][s_j][r_k]; + + // Name should be unique + getProblem().addNodalKernel(kernel_type, + prefix() + comp_name + "_" + species_name + "_decay_to_" + + Moose::stringify(_decay_products[c_i][s_j][r_k]), + params); + + /* Production of species */ + for (const auto & p_l : index_range(_decay_products[c_i][s_j][r_k])) + { + // TODO: add volumetric discretization option + const std::string kernel_type = "CoupledForceNodalKernel"; + auto params = _factory.getValidParams(kernel_type); + params.set("variable") = _decay_products[c_i][s_j][r_k][p_l]; + params.set>("v") = {species_name}; + params.set("coef") = + _branching_ratios[c_i][s_j][r_k] * _decay_constants[c_i][s_j][r_k]; + + // Add indices to be sure the name is unique + getProblem().addNodalKernel(kernel_type, + prefix() + comp_name + "_production_of_" + + _decay_products[c_i][s_j][r_k][p_l] + "_from_decay_of_" + + species_name + std::to_string(s_j) + "_" + + std::to_string(r_k) + "_" + std::to_string(p_l), + params); + } + } + } + if (_single_variable_set) + break; + } +} diff --git a/src/physics/SpeciesTrappingPhysics.C b/src/physics/SpeciesTrappingPhysics.C index 7d2cf99ac..349815412 100644 --- a/src/physics/SpeciesTrappingPhysics.C +++ b/src/physics/SpeciesTrappingPhysics.C @@ -18,8 +18,10 @@ registerMooseAction("TMAP8App", SpeciesTrappingPhysics, "copy_vars_physics"); registerMooseAction("TMAP8App", SpeciesTrappingPhysics, "check_integrity"); registerMooseAction("TMAP8App", SpeciesTrappingPhysics, "check_integrity_early_physics"); registerMooseAction("TMAP8App", SpeciesTrappingPhysics, "add_variable"); +registerMooseAction("TMAP8App", SpeciesTrappingPhysics, "add_aux_variable"); registerMooseAction("TMAP8App", SpeciesTrappingPhysics, "add_ic"); registerMooseAction("TMAP8App", SpeciesTrappingPhysics, "add_kernel"); +registerMooseAction("TMAP8App", SpeciesTrappingPhysics, "add_aux_kernel"); InputParameters SpeciesTrappingPhysics::validParams() @@ -27,6 +29,11 @@ SpeciesTrappingPhysics::validParams() InputParameters params = SpeciesPhysicsBase::validParams(); params.addClassDescription( "Add Physics for the trapping of species on multi-dimensional components."); + MooseEnum discr("nodal CG", "nodal"); + params.addParam("discretization", + discr, + "The discretization to use for the equations. We recommend 'CG' when " + "coupling to other volumetric equations"); params.addParam>( "mobile", {}, @@ -82,6 +89,7 @@ SpeciesTrappingPhysics::validParams() SpeciesTrappingPhysics::SpeciesTrappingPhysics(const InputParameters & parameters) : SpeciesPhysicsBase(parameters), + _discretization(getParam("discretization")), // If specified in the Physics block, all parameters are retrieved here _mobile_species_names({getParam>("mobile")}), _alpha_ts({getParam>("alpha_t")}), @@ -342,94 +350,166 @@ SpeciesTrappingPhysics::addFEKernels() // Time derivative if (isTransient()) { - const std::string kernel_type = "TimeDerivativeNodalKernel"; + const std::string kernel_type = isNodal() ? "TimeDerivativeNodalKernel" : "TimeDerivative"; InputParameters params = getFactory().getValidParams(kernel_type); params.set("variable") = species_name; assignBlocks(params, blocks); - getProblem().addNodalKernel(kernel_type, prefix() + species_name + "_time", params); + + const auto kernel_name = prefix() + species_name + "_time"; + if (isNodal()) + getProblem().addNodalKernel(kernel_type, kernel_name, params); + else + getProblem().addKernel(kernel_type, kernel_name, params); } // Trapping term { - const std::string kernel_type = "TrappingNodalKernel"; - auto params = _factory.getValidParams(kernel_type); - assignBlocks(params, blocks); - params.set("variable") = species_name; - params.set>("mobile_concentration") = {mobile_species_name}; - mooseAssert(c_i < _component_temperatures.size(), "Should not happen"); - // The default coupled value will not have been created by the Builder since we created - // the parameter as a MooseFunctorName in the Physics - if (MooseUtils::parsesToReal(_component_temperatures[c_i])) + auto get_trapping_params = [this, &c_i, &s_j, &blocks, &species_name, &mobile_species_name]( + const bool for_mobile_equation, + std::string & kernel_type) -> InputParameters { - std::istringstream ss(_component_temperatures[c_i]); - Real value; - ss >> value; - params.defaultCoupledValue("temperature", value, 0); - params.set>("temperature") = {}; - } - else if (getProblem().hasVariable(_component_temperatures[c_i])) - params.set>("temperature") = {_component_temperatures[c_i]}; + // We can't use a TrappingNodalKernel for the mobile equation. Even with nodal mass + // lumping, we are not getting the same results at this time. + kernel_type = for_mobile_equation + ? "TrappingKernel" + : (isNodal() ? "TrappingNodalKernel" : "TrappingKernel"); + auto params = _factory.getValidParams(kernel_type); + + assignBlocks(params, blocks); + if (!for_mobile_equation) + params.set("variable") = species_name; + else + { + params.set("variable") = mobile_species_name; + params.set>("trapped_concentration") = {species_name}; + } + + params.set>("mobile_concentration") = {mobile_species_name}; + mooseAssert(c_i < _component_temperatures.size(), "Should not happen"); + // The default coupled value will not have been created by the Builder since we created + // the parameter as a MooseFunctorName in the Physics + if (MooseUtils::parsesToReal(_component_temperatures[c_i])) + { + std::istringstream ss(_component_temperatures[c_i]); + Real value; + ss >> value; + params.defaultCoupledValue("temperature", value, 0); + params.set>("temperature") = {}; + } + else if (getProblem().hasVariable(_component_temperatures[c_i])) + params.set>("temperature") = {_component_temperatures[c_i]}; + else + paramError("temperature", "Should be a constant or the name of a variable"); + + if (!for_mobile_equation) + params.set("alpha_t") = _alpha_ts[c_i][s_j]; + else + params.set("alpha_t") = -_alpha_ts[c_i][s_j]; + params.set("trapping_energy") = _trapping_energies[c_i][s_j]; + params.set("N") = _Ns[c_i]; + params.set("Ct0") = _Ct0s[c_i][s_j]; + params.set("trap_per_free") = _trap_per_frees[c_i]; + + // Add the other species as occupying traps + std::vector copy_species; + for (const auto & sp_name : _species[c_i]) + if (sp_name != species_name) + copy_species.push_back(sp_name); + if (copy_species.size() && !getParam("different_traps_for_each_species")) + params.set>("other_trapped_concentration_variables") = + copy_species; + + return params; + }; + + std::string kernel_type; + auto params = get_trapping_params(/*for mobile equation*/ false, kernel_type); + + const auto kernel_name = + prefix() + "trapping_of_" + mobile_species_name + "_to_" + species_name; + if (isNodal()) + getProblem().addNodalKernel(kernel_type, kernel_name, params); else - paramError("temperature", "Should be a constant or the name of a variable"); - - params.set("alpha_t") = _alpha_ts[c_i][s_j]; - params.set("trapping_energy") = _trapping_energies[c_i][s_j]; - params.set("N") = _Ns[c_i]; - params.set("Ct0") = _Ct0s[c_i][s_j]; - params.set("trap_per_free") = _trap_per_frees[c_i]; - - // Add the other species as occupying traps - std::vector copy_species; - for (const auto & sp_name : _species[c_i]) - if (sp_name != species_name) - copy_species.push_back(sp_name); - if (copy_species.size() && !getParam("different_traps_for_each_species")) - params.set>("other_trapped_concentration_variables") = - copy_species; - - getProblem().addNodalKernel(kernel_type, prefix() + species_name + "_enc_trapping", params); + getProblem().addKernel(kernel_type, kernel_name, params); + + // Release term in the mobile species conservation equation + // We avoid just using the time derivative of the trapped species because another term (such + // as decay) may be acting on the trapped species concentration + { + auto params = get_trapping_params(/*for mobile equation*/ true, kernel_type); + + const auto kernel_name = + prefix() + "trapping_loss_of_" + mobile_species_name + "_to_" + species_name; + getProblem().addKernel(kernel_type, kernel_name, params); + } } // Release term { - const std::string kernel_type = "ReleasingNodalKernel"; - auto params = _factory.getValidParams(kernel_type); - assignBlocks(params, blocks); - params.set("variable") = species_name; - params.set("alpha_r") = _alpha_rs[c_i][s_j]; - params.set("detrapping_energy") = _detrapping_energies[c_i][s_j]; - // The default coupled value will not have been created by the Builder since we created - // the parameter as a MooseFunctorName in the Physics - if (MooseUtils::parsesToReal(_component_temperatures[c_i])) + auto get_releasing_params = + [this, &c_i, &s_j, &blocks, &species_name, &mobile_species_name]( + const bool for_mobile_equation, std::string & kernel_type) -> InputParameters { - std::istringstream ss(_component_temperatures[c_i]); - Real value; - ss >> value; - params.defaultCoupledValue("temperature", value, 0); - params.set>("temperature") = {}; - } - else if (getProblem().hasVariable(_component_temperatures[c_i])) - params.set>("temperature") = {_component_temperatures[c_i]}; + kernel_type = for_mobile_equation + ? "ReleasingKernel" + : (isNodal() ? "ReleasingNodalKernel" : "ReleasingKernel"); + auto params = _factory.getValidParams(kernel_type); + assignBlocks(params, blocks); + if (!for_mobile_equation) + params.set("variable") = species_name; + else + { + params.set("variable") = mobile_species_name; + params.set>("trapped_concentration") = {species_name}; + } + if (!for_mobile_equation) + params.set("alpha_r") = _alpha_rs[c_i][s_j]; + else + params.set("alpha_r") = -_alpha_rs[c_i][s_j]; + + params.set("detrapping_energy") = _detrapping_energies[c_i][s_j]; + // The default coupled value will not have been created by the Builder since we created + // the parameter as a MooseFunctorName in the Physics + if (MooseUtils::parsesToReal(_component_temperatures[c_i])) + { + std::istringstream ss(_component_temperatures[c_i]); + Real value; + ss >> value; + params.defaultCoupledValue("temperature", value, 0); + params.set>("temperature") = {}; + } + else if (getProblem().hasVariable(_component_temperatures[c_i])) + params.set>("temperature") = {_component_temperatures[c_i]}; + else + paramError("temperature", "Should be a constant or the name of a variable"); + + return params; + }; + + std::string kernel_type; + auto params = get_releasing_params(/* for mobile equation*/ false, kernel_type); + const auto kernel_name = + prefix() + "release_loss_of_" + species_name + "_to_" + mobile_species_name; + if (isNodal()) + getProblem().addNodalKernel(kernel_type, kernel_name, params); else - // there is an error right above in the trapping term - mooseAssert(false, "Should be a constant or the name of a variable"); + getProblem().addKernel(kernel_type, kernel_name, params); - getProblem().addNodalKernel(kernel_type, prefix() + species_name + "_enc_release", params); - } - - // Release term in the mobile species conservation equation - { - const std::string kernel_type = "ScaledCoupledTimeDerivative"; - auto params = _factory.getValidParams(kernel_type); - assignBlocks(params, blocks); - params.set("variable") = mobile_species_name; - params.set>("v") = {species_name}; - params.set("factor") = _trap_per_frees[c_i]; + // Release term in the mobile species conservation equation + // We avoid just using the time derivative of the trapped species because another term (such + // as decay) may be acting on the trapped species concentration + { + auto params = get_releasing_params(/*for mobile equation*/ true, kernel_type); - getProblem().addKernel( - kernel_type, prefix() + mobile_species_name + "_from_" + species_name, params); + const auto kernel_name = + prefix() + "release_of_" + mobile_species_name + "_from_" + species_name; + getProblem().addKernel(kernel_type, kernel_name, params); + } } } + + // Only need kernels once with the combined block restriction if using a single variable for all + // the components that we will be solving this Physics on if (_single_variable_set) break; } diff --git a/test/tests/actioncomponents/two_components/structure1D.i b/test/tests/actioncomponents/two_components/structure1D.i index 154e4a813..bac18d0f1 100644 --- a/test/tests/actioncomponents/two_components/structure1D.i +++ b/test/tests/actioncomponents/two_components/structure1D.i @@ -64,6 +64,7 @@ N = ${fparse 3.1622e22/cl} [trapped] species = 'trapped' mobile = 'mobile' + discretization = 'nodal' verbose = true [] [] diff --git a/test/tests/physics/more_species_on_component/species_trapping.i b/test/tests/physics/more_species_on_component/species_trapping.i index 36bcba2ae..c477cbd19 100644 --- a/test/tests/physics/more_species_on_component/species_trapping.i +++ b/test/tests/physics/more_species_on_component/species_trapping.i @@ -45,6 +45,8 @@ N = ${fparse 3.1622e22/cl} species = 'trapped' mobile = 'mobile' verbose = true + + # discretization = 'nodal' [] [] [] @@ -59,7 +61,7 @@ N = ${fparse 3.1622e22/cl} [Executioner] type = Transient - num_steps = 2 + num_steps = 2 end_time = 3 dt = .01 dtmin = .01 diff --git a/test/tests/ver-1d/ver-1d-diffusion-component.i b/test/tests/ver-1d/ver-1d-diffusion-component.i index 5ca25c2fe..7e12e1165 100644 --- a/test/tests/ver-1d/ver-1d-diffusion-component.i +++ b/test/tests/ver-1d/ver-1d-diffusion-component.i @@ -1,5 +1,5 @@ cl = 3.1622e18 -N = ${fparse 3.1622e22/cl} +N = '${fparse 3.1622e22/cl}' [ActionComponents] [structure] @@ -41,13 +41,71 @@ N = ${fparse 3.1622e22/cl} mobile = 'mobile' verbose = true + discretization = 'nodal' + dont_create_kernels = true + # Trapping parameters are specified using functors on the structure component # Releasing parameters are specified using functors on the structure component - temperature = 'temp' # we can use component ICs to set the temperature on the component + temperature = 'temp' # we can use component ICs to set the temperature on the component [] [] [] +[Kernels] + [trapped_trapping_loss_of_mobile_to_trapped] + type = TrappingKernel + Ct0 = 0.1 + N = 10000 + alpha_t = -1e+15 + block = structure + mobile_concentration = mobile + temperature = temp + trap_per_free = 1 + trapped_concentration = trapped + trapping_energy = 0 + variable = mobile + [] + [trapped_release_of_mobile_from_trapped] + type = ReleasingKernel + alpha_r = -1e+13 + block = structure + detrapping_energy = 100 + temperature = temp + trapped_concentration = trapped + variable = mobile + [] +[] + +[NodalKernels] + [trapped_trapped_time] + type = TimeDerivativeNodalKernel + block = structure + matrix_tags = 'system time' + variable = trapped + vector_tags = time + [] + [trapped_trapping_of_mobile_to_trapped] + type = TrappingNodalKernel + Ct0 = 0.1 + N = 10000 + alpha_t = 1e+15 + block = structure + mobile_concentration = mobile + temperature = temp + trap_per_free = 1 + trapping_energy = 0 + variable = trapped + [] + [trapped_release_loss_of_trapped_to_mobile] + type = ReleasingNodalKernel + alpha_r = 1e+13 + block = structure + detrapping_energy = 100 + temperature = temp + variable = trapped + [] +[] + [AuxVariables] [temp] initial_condition = 1000 @@ -78,7 +136,7 @@ N = ${fparse 3.1622e22/cl} [Executioner] type = Transient - num_steps = 2 + num_steps = 2 end_time = 3 dt = .01 dtmin = .01 @@ -86,6 +144,8 @@ N = ${fparse 3.1622e22/cl} petsc_options_iname = '-pc_type' petsc_options_value = 'lu' automatic_scaling = true + + line_search = 'none' [] [Outputs] diff --git a/test/tests/ver-1d/ver-1d-diffusion.i b/test/tests/ver-1d/ver-1d-diffusion.i index 20ebae80a..89914ba94 100644 --- a/test/tests/ver-1d/ver-1d-diffusion.i +++ b/test/tests/ver-1d/ver-1d-diffusion.i @@ -1,4 +1,4 @@ -cl=3.1622e18 +cl = 3.1622e18 temperature = 1000 [Mesh] @@ -32,15 +32,22 @@ temperature = 1000 variable = mobile extra_vector_tags = ref [] - [coupled_time] - type = CoupledTimeDerivative + # [trapping_minus_release] + # type = CoupledTimeDerivative + # variable = mobile + # v = trapped + # extra_vector_tags = ref + # [] +[] + +[NodalKernels] + [time2] + type = TimeDerivativeNodalKernel variable = mobile v = trapped - extra_vector_tags = ref + nodal_mass = mass [] -[] -[NodalKernels] [time] type = TimeDerivativeNodalKernel variable = trapped @@ -79,6 +86,15 @@ temperature = 1000 [] [] +[AuxVariables] + [mass] + [AuxKernel] + type = VolumeAux + execute_on = 'INITIAL' + [] + [] +[] + [Postprocessors] [outflux] type = SideDiffusiveFluxAverage @@ -106,11 +122,11 @@ temperature = 1000 dt = .01 dtmin = .01 solve_type = NEWTON - petsc_options_iname = '-pc_type' - petsc_options_value = 'lu' + petsc_options_iname = '-pc_type -pc_factor_shift_type' + petsc_options_value = 'lu NONZERO' automatic_scaling = true verbose = true - compute_scaling_once = false + nl_abs_tol = 1e-9 [] [Outputs] diff --git a/test/tests/ver-1dc/comparison_ver-1dc.py b/test/tests/ver-1dc/comparison_ver-1dc.py index ff4f9e2fa..74e2734bd 100644 --- a/test/tests/ver-1dc/comparison_ver-1dc.py +++ b/test/tests/ver-1dc/comparison_ver-1dc.py @@ -91,7 +91,7 @@ def summation_term(num_terms, time): label=r"Analytical", c='k', linestyle='--') ax.plot(tmap_time, tmap_prediction, label=r"TMAP8", c='tab:gray') # numerical solution ax.plot([tau_be, tau_be], [0, 3.15e18], label=r"Analytical breakthrough time", c='tab:brown', linestyle='--') -ax.plot([tmap_intercept, 17.5], [0, (17.5 - tmap_intercept) * np.max(tmap_slope)], +ax.plot([tmap_intercept, tmap_intercept], [0, 3.15e18], label=r"Numerical breakthrough time", c='tab:brown') ax.set_xlabel(u'Time (s)') ax.set_ylabel(u"Flux (atom/m$^2$s)") diff --git a/test/tests/ver-1dc/ver-1dc-components.i b/test/tests/ver-1dc/ver-1dc-components.i index addcddd76..8e6a8b6db 100644 --- a/test/tests/ver-1dc/ver-1dc-components.i +++ b/test/tests/ver-1dc/ver-1dc-components.i @@ -1,26 +1,25 @@ # This input shows the Component-Physics syntax for the ver-1dc case. # The ver-1dc case is formed without this syntax by combining ver-1dc_base.i and ver-1dc.i -epsilon_1 = ${units 100 K} -epsilon_2 = ${units 500 K} -epsilon_3 = ${units 800 K} -temperature = ${units 1000 K} +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 = 1 # m^2/s -cl = ${units 3.1622e18 atom/m^3} -N = ${units 3.1622e22 atom/m^3} +cl = '${units 3.1622e18 atom/m^3}' +N = '${units 3.1622e22 atom/m^3}' nx_num = 1000 # (-) -trapping_rate_coefficient = ${units 1e15 1/s} -release_rate_coefficient = ${units 1e13 1/s} -simulation_time = ${units 60 s} -time_interval_max = ${units 0.3 s} -time_step = ${units 1e-6 s} +trapping_rate_coefficient = '${units 1e15 1/s}' +release_rate_coefficient = '${units 1e13 1/s}' +simulation_time = '${units 60 s}' +time_interval_max = '${units 0.3 s}' +time_step = '${units 1e-6 s}' scheme = BDF2 - [ActionComponents] [structure] type = Structure1D @@ -47,30 +46,41 @@ scheme = BDF2 # Does not work for the species trapping preconditioning = 'none' + # extra_tag_vectors = 'ref' [] [] [SpeciesTrapping] [trapping] - species = 'trapped_1 trapped_2 trapped_3' - mobile = 'mobile mobile mobile' - species_initial_concentrations = '0 0 0' #'${units 1.0e-15 m^-3} ${units 1.0e-15 m^-3} ${units 1.0e-15 m^-3}' - separate_variables_per_component = false - - temperature = '${temperature}' - - alpha_t = '${trapping_rate_coefficient} ${trapping_rate_coefficient} ${trapping_rate_coefficient}' - trapping_energy = '0 0 0' - N = '${fparse N / cl}' - Ct0 = '${trapping_site_fraction_1} ${trapping_site_fraction_2} ${trapping_site_fraction_3}' - trap_per_free = 1.0e0 - different_traps_for_each_species = true - - alpha_r = '${release_rate_coefficient} ${release_rate_coefficient} ${release_rate_coefficient}' - detrapping_energy = '${epsilon_1} ${epsilon_2} ${epsilon_3}' + species = 'trapped_1 trapped_2 trapped_3' + mobile = 'mobile mobile mobile' + species_initial_concentrations = '0 0 0' #'${units 1.0e-15 m^-3} ${units 1.0e-15 m^-3} ${units 1.0e-15 m^-3}' + separate_variables_per_component = false + discretization = 'nodal' + + temperature = '${temperature}' + + alpha_t = '${trapping_rate_coefficient} ${trapping_rate_coefficient} ${trapping_rate_coefficient}' + trapping_energy = '0 0 0' + N = '${fparse N / cl}' + Ct0 = '${trapping_site_fraction_1} ${trapping_site_fraction_2} ${trapping_site_fraction_3}' + trap_per_free = 1.0e0 + # in the original ver-1dc specification + different_traps_for_each_species = true + + alpha_r = '${release_rate_coefficient} ${release_rate_coefficient} ${release_rate_coefficient}' + detrapping_energy = '${epsilon_1} ${epsilon_2} ${epsilon_3}' + + # extra_tag_vectors = 'ref' [] [] [] +# [Problem] +# type = ReferenceResidualProblem +# extra_tag_vectors = 'ref' +# reference_vector = 'ref' +# [] + [Preconditioning] [smp] type = SMP @@ -89,7 +99,9 @@ scheme = BDF2 line_search = 'none' automatic_scaling = true - nl_abs_tol = 5e-8 + nl_abs_tol = 1e-8 + # If not using reference residual + nl_forced_its = 1 [] [Postprocessors] diff --git a/test/tests/ver-1gc/ver-1gc-physics.i b/test/tests/ver-1gc/ver-1gc-physics.i new file mode 100644 index 000000000..e144e2264 --- /dev/null +++ b/test/tests/ver-1gc/ver-1gc-physics.i @@ -0,0 +1,80 @@ +concentration_A_0 = '${units 2.415e14 at/m^3 -> at/mum^3}' # atoms/microns^3 initial concentration of species A +k_1 = 0.0125 # 1/s reaction rate for first reaction +k_2 = 0.0025 # 1/s reaction rate for second reaction +end_time = 1500 # s + +[Mesh] + type = GeneratedMesh + dim = 1 +[] + +[Physics] + [SpeciesRadioactiveDecay] + [all] + species = 'c_A c_B c_C' + species_initial_concentrations = '${concentration_A_0} 0 0' + add_decaying_species_initial_conditions = true + + # | between species, ; between decay reactions + decay_products = 'c_B | c_C |' + # ; between species, ' ' between decay reactions + decay_constants = '${k_1}; ${k_2};' + [] + [] +[] + +# mK1, mK2 are used to flip the signs for the coefficient of reaction (production vs destruction) +[Materials] + [K1] + type = ADGenericConstantMaterial + prop_names = 'K1 mK1' + prop_values = '${k_1} -${k_1}' + [] + [K2] + type = ADGenericConstantMaterial + prop_names = 'K2 mK2' + prop_values = '${k_2} -${k_2}' + [] +[] + +[Postprocessors] + [concentration_A] + type = ElementAverageValue + variable = c_A + execute_on = 'INITIAL TIMESTEP_END' + [] + [concentration_B] + type = ElementAverageValue + variable = c_B + execute_on = 'INITIAL TIMESTEP_END' + [] + [concentration_C] + type = ElementAverageValue + variable = c_C + execute_on = 'INITIAL TIMESTEP_END' + [] +[] + +[Executioner] + type = Transient + scheme = bdf2 + nl_rel_tol = 1e-11 + nl_abs_tol = 1e-50 + l_tol = 1e-10 + solve_type = 'NEWTON' + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + end_time = ${end_time} + dtmax = 50 + # Ensures the time steps taken are the same + [TimeStepper] + type = ExodusTimeSequenceStepper + mesh = 'gold/ver-1gc_out.e' + [] +[] + +[Outputs] + exodus = true + csv = true +[] diff --git a/test/tests/ver-1jb/gold/ver-1jb-physics_time_dependent_out.csv b/test/tests/ver-1jb/gold/ver-1jb-physics_time_dependent_out.csv new file mode 120000 index 000000000..b0f0f16dc --- /dev/null +++ b/test/tests/ver-1jb/gold/ver-1jb-physics_time_dependent_out.csv @@ -0,0 +1 @@ +ver-1jb_time_dependent_out.csv \ No newline at end of file diff --git a/test/tests/ver-1jb/tests b/test/tests/ver-1jb/tests index 13a4c9f57..c4b3c6f4a 100644 --- a/test/tests/ver-1jb/tests +++ b/test/tests/ver-1jb/tests @@ -55,4 +55,15 @@ when modeling decay of tritium and associated growth of He in a diffusion segment with distributed traps.' required_python_packages = 'csv matplotlib numpy pandas os' [] + + [ver-1jb-physics_csvdiff] + type = CSVDiff + input = ver-1jb-physics.i + csvdiff = ver-1jb-physics_time_dependent_out.csv + requirement = 'The system shall be able to model decay of tritium and associated growth of He in a diffusion segment with + distributed traps using the shorthand Physics syntax.' + # scaled mobile is around 1e-25 so we cant trust it + # you can tighten the tolerances on ver-1jb.i to see this quantity diff as well + ignore_columns = 'tritium_mobile_inventory' + [] [] diff --git a/test/tests/ver-1jb/ver-1jb-physics.i b/test/tests/ver-1jb/ver-1jb-physics.i new file mode 100644 index 000000000..6396b5f0f --- /dev/null +++ b/test/tests/ver-1jb/ver-1jb-physics.i @@ -0,0 +1,317 @@ +# 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} +[] + +[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' + [] +[] + +[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 + # Scaling here + normal_factor = ${fparse density_scalar * trap_per_free} + source_variable = tritium_trapped_concentration_scaled + [] + [helium_concentration] + type = NormalizationAux + variable = helium_concentration + normal_factor = ${density_scalar} + source_variable = helium_concentration_scaled + [] +[] + +[Physics] + [SpeciesDiffusionReaction] + [decays] + species = 'tritium_mobile_concentration_scaled helium_concentration_scaled' + initial_conditions_species = '${fparse tritium_mobile_concentration_initial / density_scalar} 0' + diffusivity_matprops = 'diffusivity 0' + + # we dont set the decay for the trapped helium concentration because it's using nodal + # kernels in the original input which we are trying to reproduce + reacting_species = 'tritium_mobile_concentration_scaled' + product_species = 'helium_concentration_scaled' + reaction_coefficients = '${decay_rate_constant}' + [] + [] + [SpeciesTrapping] + [trapped] + species = 'tritium_trapped_concentration_scaled' + mobile = 'tritium_mobile_concentration_scaled' + species_initial_concentrations = 'tritium_trapped_concentration_initial_function_scaled' + + # Since we are not using components, we specify all the species trapping parameters here + temperature = 'temperature' + trap_per_free = ${trap_per_free} + 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) + alpha_r = ${fparse tritium_release_prefactor} # (1/s) + detrapping_energy = ${fparse tritium_release_energy / boltzmann_constant} # (K) + [] + [] +[] + + +# Alternatively, we could create a nodal decay Physics OR add decay reactions to SpeciesTrapping +# Note that this approach, used by the reference input, is not conservative because of the differences +# in integration of the flux from one species to the other between the nodal and volumetric kernels. +[Kernels] + [decay_trapped_to_helium] + type = ADMatReactionFlexible + reaction_rate_name = ${fparse decay_rate_constant} + variable = helium_concentration_scaled + vs = tritium_trapped_concentration_scaled + [] +[] +[NodalKernels] + [decay_trapped_to_helium] + 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 = ADDerivativeParsedMaterial + 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 = ADParsedMaterial + 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 = ADParsedMaterial + property_name = 'alpha_r_tot' + coupled_variables = 'temperature' + expression = '${tritium_release_prefactor} * exp(- ${tritium_release_energy} / ${boltzmann_constant}/temperature)' + outputs = 'all' + [] +[] + +[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-physics_out + [time_dependent_out] + type = CSV + execute_vector_postprocessors_on = NONE + file_base = ver-1jb-physics_time_dependent_out + [] + [profile_out] + type = CSV + sync_only = true + sync_times = '${units 45 year -> s}' + execute_postprocessors_on = NONE + file_base=ver-1jb-physics_profile_out + [] +[] + +[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' + [] +[] diff --git a/test/tests/ver-1jb/ver-1jb.i b/test/tests/ver-1jb/ver-1jb.i index 05cb32ae1..e00c2102f 100644 --- a/test/tests/ver-1jb/ver-1jb.i +++ b/test/tests/ver-1jb/ver-1jb.i @@ -174,19 +174,6 @@ dt_max = ${fparse end_time/100} # s v = 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 @@ -222,13 +209,21 @@ dt_max = ${fparse end_time/100} # s temperature = temperature # (K) trap_per_free = ${trap_per_free} [] - [release] + [release_from_trapped] 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) [] + [release_into_mobile] + type = ReleasingNodalKernel + variable = tritium_mobile_concentration_scaled # (atoms/m^3) / density_scalar = (-) + trapped_concentration = tritium_trapped_concentration_scaled + 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 = (-)