val-2f
Modelling self-damaged tungsten effects on deuterium transport
Case Description
The case being used for validation here involves the use of recrystallized polycrystalline tungsten (PCW) samples, which are subjected to ion irradiation and subsequent analysis using thermal desorption spectroscopy (TDS). The primary objective is to determine how damage influences deuterium trapping and release. This case is drawn and updated from Dark et al. (2024).
The TDS process is simulated using TMAP8 in a 1D tungsten sample with a thickness of 0.8 mm. The TDS simulation consisted of three phases: implantation, cooldown, and desorption.
The sample temperature histories are shown in Figure 1.
Figure 1: Temperature history.
Model Description
1- Implantation phase:
During the charging phase, deuterium is continuously implanted into the tungsten sample over a period of 72 hours. The temperature is maintained at 370 K throughout this phase. The surface is exposed to a constant flux of atoms/m/s, corresponding to a total fluence of atoms/m. The implantation profile follows a Gaussian distribution centered at the mean implantation depth of nm, with a standard deviation of nm :
(1)
where the surface flux function is given by:
The implantation distribution is illustrated in Figure 2: most of the implanted atoms are found within a few standard deviations () of the mean implantation depth (). In this context, it means that the mesh in the region where deuterium implantation occurs should be refined to a size comparable to . In this TMAP8 simulation, the first mesh region is set to length of , divided into 100 elements. This allows to capture the majority of the implantation profile and ensure that the mesh is sufficiently refined in this region.
Figure 2: Deuterium implantation.
2- Cooldown phase:
After the implantation phase, the system enters the cooldown phase, lasting 12 hours. During this period, the sample temperature is rapidly reduced from 370 K to 295 K. No additional deuterium is introduced during this phase, meaning the source term is set to zero.
3- Desorption phase:
The final stage of the simulation is the desorption phase, during which the sample is gradually heated from 300 K to 1000 K at a constant rate of K/s.
Governing equations
• Damaged induced traps
The evolution of the trapping site density, , over time under irradiation is described by (Dark et al., 2024):
(2)
This equation consists of two terms. The first term, , represents the creation of trapping sites. It increases the trap density up to a saturation density, , due to the damage imposed on the sample. Here, is the damage rate in (dpa/s) and is the trap creation factor in (traps/m/dpa). The second term, , accounts for the annealing effect which decreases the density of trapping sites. is the trap annealing rate in (1/s) and follows an Arrhenius law given by (3) where is a pre-exponential factor, is the activation energy, is the Boltzmann constant, and is the temperature. The parameters for the evolution of the irradiation-induced traps are listed in Table 2.
The analytical solution for the ODE is given by:
(4)
The computed trap densities are provided in Table 3.
• Deuterium transport equations
The general form of the transport equations for the deuterium in tungsten is given by:
(5)
and, for representing the trapping sites:
(6)
and
where is the concentration of mobile tritium, is the time, is the source term in sample due to the deuterium implantation (Eq. (1)), is the trapped species in trap , and are the trapping and release rate coefficients for trap , is a user-defined numerical factor scaling to be closer to for better numerical convergence, is the fraction of host sites that can contribute to trapping, is the concentration of empty trapping sites, and is the host density, and is the deuterium diffusivity in tungsten, which is defined as:
(7)
where is the diffusion activation energy, is the temperature, and is the maximum diffusivity coefficient.
and are defined as:
(8)
and
(9)
where and are the pre-exponential factors of trapping and release. The trapping energy is equal to the diffusion activation energy .
• Boundary conditions
There are several ways to model the surface behavior of deuterium, and this study tested two options. One can assume that every deuterium atom at the surface is immediately desorbed from the materials, effectively assuming an infinite recombination rate at the surface and imposing a null Dirichlet boundary condition for the mobile deuterium concentration ( atoms/m). This is the assumption used in Dark et al. (2024) and reproduced in Figure 5.
However, it is also possible to capture the rate of deuterium recombination at the surface as it recombines into gas. It can be described by the following surface flux:
where represents the recombination flux exiting the sample on both the left and right sides, is the surface area, and is the deuterium recombination coefficient. The coefficient of 2 accounts for the fact that 2 deuterium atoms combine to form one D molecule. Using this condition, the surface concentration will not be imposed, but be governed by the concentration of deuterium coming to the surface from the bulk and the rate of recombination. This boundary condition is used in Figure 7.
Case and Model Parameters
In this section provides the model parameters for this validation case. Table 1 provides the base model parameters for the validation case. Then, we presents two sets of parameter values for the trapping sites. First, in Table 2, we lists the parameters related to the damage-induced traps. These parameters capture the irradiation-dependent evolution of trap properties. Second, in Table 3, we reproduce the calibrated values used in Dark et al. (2024) to reproduce the desorption data. As illustrated in Figure 3, the trap densities in these two approaches, i.e., analytical solution in Eq. (4) and calibrated values in Figure 3, differ.
Table 1: Values of material properties.
| Parameter | Description | Value | Units | Reference |
|---|---|---|---|---|
| Boltzmann constant | 1.380649 | J/K | PhysicalConstants.h | |
| Initial temperature | 370 | K | Dark et al. (2024) | |
| Cooldown temperature | 295 | K | Dark et al. (2024) | |
| Desorption start temperature | 300 | K | Dark et al. (2024) | |
| Desorption end temperature | 1000 | K | Dark et al. (2024) | |
| Desorption heating rate | 0.05 | K/s | Dark et al. (2024) | |
| Charging time | 72 | h | Dark et al. (2024) | |
| Cooldown duration | 12 | h | Dark et al. (2024) | |
| Diffusion pre-factor | 1.6 | m/s | Dark et al. (2024) | |
| Activation energy for deuterium diffusion | 0.28 | eV | Dark et al. (2024) | |
| Mean implantation depth | 0.7 | nm | Dark et al. (2024) | |
| Standard deviation of implantation profile | 0.5 | nm | Dark et al. (2024) | |
| Incident fluence | 1.5 | atoms/m | Dark et al. (2024) | |
| Incident flux | 5.79 | atoms/m/s | Dark et al. (2024) | |
| Length of the tungsten sample | 0.8 | mm | Dark et al. (2024) | |
| Tungsten density | 6.3222 | at/m | J-Dark-PhD |
Table 2: Damaged-induced traps parameters from Dark et al. (2024) used in Eq. (2) to capture the evolution of trap properties with irradiation.
| Trap | (dpa/s) | (traps/m/dpa) | (atoms/m) | (1/s) | (eV) | (K) |
|---|---|---|---|---|---|---|
| Trap 1 | 8.9 | 9.0 | 6.9 | 6.18 | 0.24 | 800 |
| Trap 2 | 8.9 | 4.2 | 7.0 | 6.18 | 0.24 | 800 |
| Trap 3 | 8.9 | 2.5 | 6.0 | 6.18 | 0.30 | 800 |
| Trap 4 | 8.9 | 5.0 | 4.7 | 6.18 | 0.30 | 800 |
| Trap 5 | 8.9 | 1.0 | 2.0 | 0 | - | 800 |
Table 3: Calibrated values of traps parameters for 0.1 dpa from Dark et al. (2024) to match the experimental desorption data.
| Parameter | Description | Value | Units | Reference |
|---|---|---|---|---|
| Pre-factor of trapping rate coefficient | atoms/s | J-Dark-PhD | ||
| Pre-factor of release rate coefficient | atoms/s | J-Dark-PhD | ||
| Trapping energy for all traps | 0.28 | eV | Dark et al. (2024) | |
| Release energy for trap 1 | 1.15 | eV | Dark et al. (2024) | |
| Release energy for trap 2 | 1.35 | eV | Dark et al. (2024) | |
| Release energy for trap 3 | 1.65 | eV | Dark et al. (2024) | |
| Release energy for trap 4 | 1.85 | eV | Dark et al. (2024) | |
| Release energy for trap 5 | 2.05 | eV | Dark et al. (2024) | |
| Release energy for intrinsic trap | 1.04 | eV | Dark et al. (2024) | |
| Density for trap 1 | 4.63 | atoms/m | Eq. (4) | |
| Density for trap 2 | 2.87 | atoms/m | Eq. (4) | |
| Density for trap 3 | 1.96 | atoms/m | Eq. (4) | |
| Density for trap 4 | 2.97 | atoms/m | Eq. (4) | |
| Density for trap 5 | 7.87 | atoms/m | Eq. (4) | |
| Density for intrinsic trap | 2.4 | atoms/m | Dark et al. (2024) |
 at 0.1 dpa.](figures/val-2f_trap_induced_density.png)
Figure 3: Trap densities as a function of damage. The analytical solutions are represented by the solid curve, while the points indicate the fitted values used in Dark et al. (2024) at 0.1 dpa.
As described in the boundary conditions section, we use two different approaches for the description of the surface recombination rates. We first impose atoms/m to reproduce the results from Dark et al. (2024). However, we also capture the finite kinetics of surface reaction using the recombination rate from Zhao et al. (2020), which has been increased in this study to better match the TDS data. The values used are listed in Table 4:
Table 4: Values of recombination rates.
| Parameter | Description | Value | Units | Reference |
|---|---|---|---|---|
| Deuterium recombination coefficient | 3.8 | m/at/s | Zhao et al. (2020) | |
| Adapted deuterium recombination coefficient | 3.8 | m/at/s | - |
warningwarning:Typo in formula from Zhao et al. (2020)
There is a typo in the expression for the deuterium recombination coefficient for clean tungsten surfaces from Zhao et al. (2020) where the minus sign in the exponential is missing, even though the data shows it should be present. Consequently, we used the corrected value in our simulations, which includes the minus sign.
Results
The figures below show the comparison of the TMAP8 calculation and the experimental data during desorption. The experimental data are provided by T. Schwarz-Selinger and are available here.
We first reproduce the results from Dark et al. (2024) using a Dirichlet boundary condition at the surface and using the calibrated trapping properties from Table 3. Figure 4 displays the quantities of mobile, trapped, and desorbing deuterium atoms during the desorption process. During desorption, the temperature increases from 300 K to 1000 K. The amount of deuterium trapped will decrease as the temperature rises and the various trapping energies are reached, meaning that deuterium will leave the traps, become mobile, and diffuse out. During desorption, no further implantation occurs, resulting in a decrease in the number of mobile and trapped deuterium atoms and an increase in the number of desorbed deuterium atoms.
Mass conservation in Figure 4 is maintained during desorption, with a 0.01% root mean squared percentage error (RMSPE) between the initial number of mobile and trapped deuterium atoms and the total number of deuterium atoms (mobile, trapped, and desorbed).
Figure 5 shows a good alignment with the experimental data and the results presented in Dark et al. (2024). With an infinite recombination rate, the amount of mobile deuterium is minimal compared to the trapped deuterium, as demonstrated in Figure 4. Consequently, the desorbed deuterium flux is almost solely due to the deuterium desorbing from the traps.
Figure 4: Quantity of deuterium atoms during the desorption process for an infinite recombination rate and calibrated trapping properties.
Figure 5: Comparison of TMAP8 calculations with experimental data on deuterium flux (atoms/m/s) for a damage of 0.1 dpa, an infinite recombination rate and calibrated trapping properties.
In this validation case, the recombination rate has been set to a finite value to better reflect real conditions. Specifically, the recombination rate is set to a value 10 orders of magnitude higher than that reported by Zhao et al. (2020). A recombination rate that is too low leads to the accumulation of a significant amount of mobile deuterium at the onset of the desorption phase. Consequently, this would result in desorption being attributed to mobile deuterium rather than trapped deuterium. The updated, increased recombination rate limits deuterium accumulation in the sample, without imposing a surface concentration. As shown in Figure 6, the quantity of mobile deuterium is nearly equal to zero, similar to the case with an infinite recombination rate, and mass conservation is maintained. The TMAP8 results in Figure 7 match well with the experimental data. There still exists an offset between the TMAP8 results and the experimental data, which can be reduced by better fitting the recombination rate and trapping parameters described in Table 3 using the MOOSE's stochastic tools module. It is important to mention here that the mesh for the infinite recombination rate requires a finer mesh as well as smaller time steps compared to the finite recombination rate case.
Figure 6: Quantity of deuterium atoms during the desorption process for an adapted recombination rate and uncalibrated trapping properties.
Figure 7: Comparison of TMAP8 calculations with experimental data on deuterium flux (atoms/m/s) for a damage of 0.1 dpa, an adapted recombination rate and uncalibrated trapping properties.
Input files
The input file for this case can be found at (test/tests/val-2f/val-2f.i). To minimize the length of the input file and organize it, it is divided into several parts:
(test/tests/val-2f/parameters_val-2f.params) lists the key values and model parameters used in this simulation
(test/tests/val-2f/val-2f_trapping_intrinsic.i) provides the blocks necessary to introduce the intrinsic traps in the simulation
(test/tests/val-2f/val-2f_trapping_1.i), (test/tests/val-2f/val-2f_trapping_2.i), (test/tests/val-2f/val-2f_trapping_3.i), (test/tests/val-2f/val-2f_trapping_4.i), (test/tests/val-2f/val-2f_trapping_5.i) provide the blocks necessary to introduce the trapping sites 1, 2, 3, 4, and 5, respectively, in the simulation.
To combine them into one input file when running the simulation, (test/tests/val-2f/val-2f.i) uses the !include feature.
schooltip:Input file include syntax information
To learn more about the !include feature, refer to the Input File Syntax page.
Note that both surface conditions can be modeled using (test/tests/val-2f/val-2f.i). By running it as is, the recombination condition with the updated recombination rate is utilized. (test/tests/val-2f/tests) used cli_args to modify (test/tests/val-2f/val-2f.i) into using the effectively infinite recombination rate at the surface and impose a null concentration at the surfaces.
To limit the computational costs of the test case, the test runs a version of the file with a smaller and coarser mesh, and fewer time steps. In addition, this validation case includes tests using two time integration schemes. The implicit-euler scheme (first-order, single-step) provides higher stability and consistency of the solution across different operating systems and processor core counts, making it better suited for use in Bayesian optimization workflows where numerical stability is important. The bdf2 scheme (second-order, multi-step) provides higher accuracy and is used for final results calculation. Both schemes are tested for this validation case in (test/tests/val-2f/tests). More information about the changes can be found in the test specification file for this case, namely (test/tests/val-2f/tests).
References
- James Dark, Rémi Delaporte-Mathurin, Thomas Schwarz-Selinger, Etienne A Hodille, Jonathan Mougenot, Yann Charles, and Christian Grisolia.
Modelling neutron damage effects on tritium transport in tungsten.
Nuclear Fusion, 64(8):086026, 2024.[Export]
BibTeX
@article{dark2024modelling, author = "Dark, James and Delaporte-Mathurin, R{\'e}mi and Schwarz-Selinger, Thomas and Hodille, Etienne A and Mougenot, Jonathan and Charles, Yann and Grisolia, Christian", title = "Modelling neutron damage effects on tritium transport in tungsten", journal = "Nuclear Fusion", volume = "64", number = "8", pages = "086026", year = "2024", publisher = "IOP Publishing" }RIS
TY - JOUR AU - Dark, James AU - Delaporte-Mathurin, Rémi AU - Schwarz-Selinger, Thomas AU - Hodille, Etienne A AU - Mougenot, Jonathan AU - Charles, Yann AU - Grisolia, Christian TI - Modelling neutron damage effects on tritium transport in tungsten JO - Nuclear Fusion PY - 2024 VL - 64 IS - 8 PB - IOP Publishing SP - 086026 ER -Plain Text
James Dark, Rémi Delaporte-Mathurin, Thomas Schwarz-Selinger, Etienne A Hodille, Jonathan Mougenot, Yann Charles, and Christian Grisolia. Modelling neutron damage effects on tritium transport in tungsten. Nuclear Fusion, 64(8):086026, 2024. - Mingzhong Zhao, Shota Yamazaki, Takuro Wada, Ayaka Koike, Fei Sun, Naoko Ashikawa, Yoji Someya, Tetsu Mieno, and Yasuhisa Oya.
Deuterium recombination coefficient on tungsten surface determined by plasma driven permeation.
Fusion Engineering and Design, 160:111853, 2020.[Export]
BibTeX
@article{zhao2020deuterium, author = "Zhao, Mingzhong and Yamazaki, Shota and Wada, Takuro and Koike, Ayaka and Sun, Fei and Ashikawa, Naoko and Someya, Yoji and Mieno, Tetsu and Oya, Yasuhisa", title = "Deuterium recombination coefficient on tungsten surface determined by plasma driven permeation", journal = "Fusion Engineering and Design", volume = "160", pages = "111853", year = "2020", publisher = "Elsevier" }RIS
TY - JOUR AU - Zhao, Mingzhong AU - Yamazaki, Shota AU - Wada, Takuro AU - Koike, Ayaka AU - Sun, Fei AU - Ashikawa, Naoko AU - Someya, Yoji AU - Mieno, Tetsu AU - Oya, Yasuhisa TI - Deuterium recombination coefficient on tungsten surface determined by plasma driven permeation JO - Fusion Engineering and Design PY - 2020 VL - 160 PB - Elsevier SP - 111853 ER -Plain Text
Mingzhong Zhao, Shota Yamazaki, Takuro Wada, Ayaka Koike, Fei Sun, Naoko Ashikawa, Yoji Someya, Tetsu Mieno, and Yasuhisa Oya. Deuterium recombination coefficient on tungsten surface determined by plasma driven permeation. Fusion Engineering and Design, 160:111853, 2020.
(test/tests/val-2f/val-2f.i)
# Validation Problem #2f
# Self-damaged Tungsten Effects on Deuterium Transport extended from an original model from
# Dark, J., Delaporte-Mathurin, R., Schwarz-Selinger, T., Hodille, E. A., Mougenot, J.,
# Charles, Y., & Grisolia, C. (2024). Modelling neutron damage effects on tritium transport
# in tungsten. Nuclear Fusion, 64(8), 086026.
!include parameters_val-2f.params
!include val-2f_trapping_intrinsic.i
!include val-2f_trapping_5.i
!include val-2f_trapping_4.i
!include val-2f_trapping_3.i
!include val-2f_trapping_2.i
!include val-2f_trapping_1.i
[Mesh]
[cartesian_mesh]
type = CartesianMeshGenerator
dim = 1
dx = '${dx1_hat} ${dx2_hat} ${dx3_hat} ${dx4_hat} ${dx5_hat}'
ix = '${ix1} ${ix2} ${ix3} ${ix4} ${ix5}'
subdomain_id = '0 0 0 0 0'
[]
[]
[Variables]
[deuterium_concentration_W]
[]
[]
[AuxVariables]
active = 'bounds_dummy temperature'
[bounds_dummy]
order = FIRST
family = LAGRANGE
[]
[temperature]
initial_condition = ${temperature_initial}
[]
[]
[Bounds]
[deuterium_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = deuterium_concentration_W
bound_type = lower
bound_value = 0
[]
[]
[Kernels]
[time_W]
type = TimeDerivative
variable = deuterium_concentration_W
[]
[diffusion_W]
type = ADMatDiffusion
variable = deuterium_concentration_W
diffusivity = diffusivity_W
[]
[source_deuterium]
type = BodyForce
variable = deuterium_concentration_W
function = source_deuterium
[]
[]
[AuxKernels]
[temperature_aux]
type = FunctionAux
variable = temperature
function = temperature_bc_func
execute_on = 'INITIAL LINEAR'
[]
[]
[BCs]
active = 'left_recombination_flux right_recombination_flux'
# Kinetic boundary conditions
[left_recombination_flux]
type = ADMatNeumannBC
variable = deuterium_concentration_W
boundary = left
value = 1
boundary_material = flux_recombination_surface
[]
[right_recombination_flux]
type = ADMatNeumannBC
variable = deuterium_concentration_W
boundary = right
value = 1
boundary_material = flux_recombination_surface
[]
# Sieverts boundary conditions
[left_concentration_sieverts]
type = ADDirichletBC
value = '${sieverts_boundary_hat}'
boundary = left
variable = deuterium_concentration_W
[]
[right_concentration_sieverts]
type = ADDirichletBC
value = '${sieverts_boundary_hat}'
boundary = right
variable = deuterium_concentration_W
[]
[]
[Functions]
[temperature_bc_func]
type = ParsedFunction
expression = 'if(t<${charge_time_hat}, ${temperature_initial},
if(t<${fparse charge_time_hat + cooldown_duration_hat}, ${temperature_cooldown},
${temperature_desorption_min}+${desorption_heating_rate_hat}*(t-${fparse charge_time_hat + cooldown_duration_hat})))'
[]
[source_distribution]
type = ParsedFunction
expression = '1 / (${sigma_hat} * sqrt(2 * pi)) * exp(-0.5 * ((x - ${R_p_hat}) / ${sigma_hat}) ^ 2)'
[]
[surface_flux_func]
type = ParsedFunction
expression = 'if(t<${charge_time_hat}, ${surface_flux_hat}, 0)'
[]
[source_deuterium]
type = ParsedFunction
symbol_names = 'source_distribution surface_flux_func'
symbol_values = 'source_distribution surface_flux_func'
expression = 'source_distribution * surface_flux_func'
[]
[max_dt_size_function]
type = ParsedFunction
expression = 'if(t<${fparse 5 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 8 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 12 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 20 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 35 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 450 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 5000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse 11000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 13000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse (charge_time + cooldown_duration + 4500) / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 313000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 315000 / time_reference}, ${fparse 1e1 / time_reference}, ${fparse 1e3 / time_reference}))))))))))))'
[]
[max_dt_size_function_inf]
type = ParsedFunction
expression = 'if(t<${fparse 5 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 8 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 12 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 20 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 35 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 450 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 5000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse 11000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 13000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse (charge_time + cooldown_duration + 4500) / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 315000 / time_reference}, ${fparse 1e1 / time_reference}, ${fparse 1e3 / time_reference})))))))))))'
[]
[max_dt_size_function_coarse]
type = ParsedFunction
expression = 'if(t<${fparse 1e-1 / time_reference}, ${fparse 1e4 / time_reference}, ${fparse 1e5 / time_reference})'
[]
[]
[Materials]
active = 'diffusivity_W_func diffusivity_nonAD recombination_rate_surface flux_recombination_surface'
[diffusivity_W_func]
type = ADDerivativeParsedMaterial
property_name = 'diffusivity_W'
functor_names = 'temperature_bc_func'
functor_symbols = 'temperature'
expression = '${diffusion_W_preexponential_hat} * exp(- ${diffusion_W_energy} / ${kb_eV} / temperature)'
[]
[diffusivity_nonAD]
type = MaterialADConverter
ad_props_in = 'diffusivity_W'
reg_props_out = 'diffusivity_W_nonAD'
[]
[recombination_rate_surface]
type = ADDerivativeParsedMaterial
property_name = 'Kr'
functor_names = 'temperature_bc_func'
functor_symbols = 'temperature'
expression = '${recombination_coefficient_hat} * exp(- ${recombination_energy} / ${kb_eV} / temperature)'
[]
[flux_recombination_surface]
type = ADDerivativeParsedMaterial
coupled_variables = 'deuterium_concentration_W'
property_name = 'flux_recombination_surface'
material_property_names = 'Kr'
expression = '- 2 * Kr * deuterium_concentration_W ^ 2'
[]
[]
[Postprocessors]
active = 'integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration
scaled_mobile_deuterium flux_surface_left scaled_flux_surface_left
flux_surface_right scaled_flux_surface_right temperature_pps diffusion_W_hat diffusion_W
max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1
integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3
integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5
integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic
spatial_max_mobile_d2 spatial_max_trapped_1 spatial_max_trapped_2 spatial_max_trapped_3 spatial_max_trapped_4 spatial_max_trapped_5 spatial_max_trapped_intrinsic
max_mobile_d2 max_trapped_1 max_trapped_2 max_trapped_3 max_trapped_4 max_trapped_5 max_trapped_intrinsic max_scaled_flux_surface_left max_scaled_flux_surface_right
max_scaled_mobile_deuterium max_scaled_trapped_deuterium_intrinsic'
[integral_source_deuterium]
type = FunctionElementIntegral
function = source_deuterium
outputs = none
[]
[scaled_implanted_deuterium]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = integral_source_deuterium
[]
[integral_deuterium_concentration]
type = ElementIntegralVariablePostprocessor
variable = deuterium_concentration_W
outputs = none
[]
[scaled_mobile_deuterium]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_deuterium_concentration
[]
[max_scaled_mobile_deuterium]
type = TimeExtremeValue
postprocessor = scaled_mobile_deuterium
value_type = max
outputs = 'console'
[]
[flux_surface_left]
type = ADSideAverageMaterialProperty
boundary = 'left'
property = flux_recombination_surface
outputs = none
[]
[scaled_flux_surface_left]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_left
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_left]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_left
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[flux_surface_left_sieverts]
type = SideDiffusiveFluxAverage
variable = deuterium_concentration_W
boundary = 'left'
diffusivity = 'diffusivity_W_nonAD'
outputs = none
[]
[scaled_flux_surface_left_sieverts]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_left_sieverts
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_left_sieverts]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_left_sieverts
value_type = max
[]
[flux_surface_right]
type = ADSideAverageMaterialProperty
boundary = 'right'
property = flux_recombination_surface
outputs = none
[]
[scaled_flux_surface_right]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_right
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_right]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_right
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[flux_surface_right_sieverts]
type = SideDiffusiveFluxAverage
variable = deuterium_concentration_W
boundary = 'right'
diffusivity = 'diffusivity_W_nonAD'
outputs = none
[]
[scaled_flux_surface_right_sieverts]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_right_sieverts
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_right_sieverts]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_right_sieverts
value_type = max
[]
[temperature_pps]
type = ElementAverageValue
variable = temperature
execute_on = 'initial timestep_end'
[]
[diffusion_W_hat]
type = ElementAverageMaterialProperty
mat_prop = diffusivity_W_nonAD
outputs = none
[]
[diffusion_W]
type = ScalePostprocessor
scaling_factor = '${fparse length_reference ^ 2 / time_reference}'
value = diffusion_W_hat
outputs = none
[]
[max_time_step_size]
type = FunctionValuePostprocessor
function = max_dt_size_function
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[max_time_step_size_coarse]
type = FunctionValuePostprocessor
function = max_dt_size_function_coarse
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[spatial_max_mobile_d2]
type = NodalExtremeValue
value_type = 'max'
variable = deuterium_concentration_W
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_1]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_1
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_2]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_3]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_3
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_4]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_4
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_5]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_5
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_intrinsic]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_intrinsic
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_mobile_d2]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_mobile_d2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_1]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_1
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_2]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_3]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_3
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_4]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_4
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_5]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_5
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_intrinsic]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_intrinsic
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = 'Newton'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -snes_type'
petsc_options_value = 'lu mumps vinewtonrsls'
end_time = ${endtime_hat}
line_search = 'none'
nl_rel_tol = 1e-8
nl_abs_tol = 4e-5
nl_max_its = 34
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse dt_init / time_reference}'
growth_factor = 1.1
timestep_limiting_postprocessor = max_time_step_size
[]
[Predictor]
type = SimplePredictor
scale = 1.0
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Outputs]
[csv]
type = CSV
file_base = 'val-2f_out'
[]
[csv_temperature_history]
type = CSV
show = 'temperature_pps'
file_base = 'val-2f_temperature_out'
[]
[exodus]
type = Exodus
output_material_properties = true
time_step_interval = 200
file_base = 'val-2f_out'
[]
[]
(test/tests/val-2f/parameters_val-2f.params)
# This input files provides the model parameter values for validation case val-2f.
# It is included in val-2f.i
# Physical constants
kb = '${units 1.380649e-23 J/K}' # Boltzmann constant J/K - from PhysicalConstants.h
eV_to_J = '${units 1.602176634e-19 J/eV}' # Conversion coefficient from eV to Joules - from PhysicalConstants.h
kb_eV = '${units ${fparse kb / eV_to_J} eV/K}' # Boltzmann constant eV/K
tungsten_density = '${units 6.3222e28 at/m^3 -> at/mum^3}'
# Temperature conditions
temperature_initial = '${units 370 K}'
temperature_cooldown = '${units 295 K}'
temperature_desorption_min = '${units 300 K}'
temperature_desorption_max = '${units 1000 K}'
desorption_heating_rate = '${units ${fparse 3/60} K/s}'
# Important times
charge_time = '${units 72 h -> s}'
cooldown_duration = '${units 12 h -> s}'
desorption_duration = '${fparse (temperature_desorption_max-temperature_desorption_min)/desorption_heating_rate}'
dt_init = 1e-8
# Materials properties
diffusion_W_preexponential = '${units 1.6e-7 m^2/s -> mum^2/s}'
diffusion_W_energy = '${units 0.28 eV}'
recombination_coefficient = '${units ${fparse 3.8e-16} m^4/at/s -> mum^4/at/s}'
recombination_energy = '${units 0.34 eV}'
trapping_energy = '${fparse ${units ${diffusion_W_energy} eV -> J} / kb}'
trapping_prefactor = '${units ${fparse ${units ${diffusion_W_preexponential} mum^2/s -> m^2/s} / (1.1e-10^2 * 6) } 1/s}'
detrapping_prefactor = '${units 1e13 1/s}'
# Source term parameters
sigma = '${units 0.5e-9 m -> mum}'
R_p = '${units 0.7e-9 m -> mum}'
flux = '${units ${fparse 5.79e19} at/m^2/s -> at/mum^2/s}'
# Fermi-Dirac trap density distribution
depth_center = '${units 2.5e-6 m -> mum}'
depth_width = '${units 5e-7 m -> mum}'
sample_thickness = '${units 0.8e-3 m -> mum}'
dx1 = '${fparse 6*sigma}'
dx2 = '${units 3e-6 m -> mum}'
dx3 = '${units 7.997e-6 m -> mum}'
dx4 = '${units 7.2e-5 m -> mum}'
ix1 = 100
ix2 = 100
ix3 = 100
ix4 = 100
ix5 = 100
# Damaged-induced traps parameters
A0 = '${units 6.18e-3 1/s}'
phi = '${units 8.9e-5 dpa/s}'
temperature_damage = '${units 800 K}'
damage = '${units 0.1 dpa}'
damage_duration = '${units ${fparse damage/phi} s}'
# Trap 1 parameters
trapping_energy_1 = '${fparse trapping_energy}'
trapping_prefactor_1 = '${fparse trapping_prefactor}'
detrapping_energy_1 = '${fparse ${units 1.15 eV -> J} / ${kb}}'
detrapping_prefactor_1 = '${fparse detrapping_prefactor}'
K_1 = '${units 9.0e26 traps/m^3/dpa}'
nmax_1 = '${units 6.9e25 1/m^3}'
Ea_1 = '${units 0.24 eV}'
A_1 = '${fparse A0 * exp(-Ea_1/(kb_eV*temperature_damage))}'
S_1 = '${fparse phi*K_1}'
r_1 = '${fparse phi*K_1/nmax_1 + A_1}'
trap_density_01dpa_1 = '${units ${fparse -S_1/r_1 * exp(-r_1*damage_duration) + S_1/r_1} at/m^3 -> at/mum^3}'
trapping_site_fraction_1 = '${fparse trap_density_01dpa_1 / tungsten_density}'
trap_density_01dpa_1_inf = '${units 4.8e25 at/m^3 -> at/mum^3}'
trapping_site_fraction_1_inf = '${fparse trap_density_01dpa_1_inf / tungsten_density}'
trap_concentration_reference_1 = '${fparse tungsten_density * 1e-8}'
# Trap 2 parameters
trapping_energy_2 = '${fparse trapping_energy}'
trapping_prefactor_2 = '${fparse trapping_prefactor}'
detrapping_energy_2 = '${fparse ${units 1.35 eV -> J} / kb}'
detrapping_prefactor_2 = '${fparse detrapping_prefactor}'
K_2 = '${units 4.2e26 traps/m^3/dpa}'
nmax_2 = '${units 7.0e25 1/m^3}'
Ea_2 = '${units 0.24 eV}'
A_2 = '${fparse A0 * exp(-Ea_2/(kb_eV*temperature_damage))}'
S_2 = '${fparse phi*K_2}'
r_2 = '${fparse phi*K_2/nmax_2 + A_2}'
trap_density_01dpa_2 = '${units ${fparse -S_2/r_2 * exp(-r_2*damage_duration) + S_2/r_2} at/m^3 -> at/mum^3}'
trapping_site_fraction_2 = '${fparse trap_density_01dpa_2 / tungsten_density}'
trap_density_01dpa_2_inf = '${units 3.8e25 at/m^3 -> at/mum^3}'
trapping_site_fraction_2_inf = '${fparse trap_density_01dpa_2_inf / tungsten_density}'
trap_concentration_reference_2 = '${fparse tungsten_density * 1e-8}'
# Trap 3 parameters
trapping_energy_3 = '${fparse trapping_energy}'
trapping_prefactor_3 = '${fparse trapping_prefactor}'
detrapping_energy_3 = '${fparse ${units 1.65 eV -> J} / kb}'
detrapping_prefactor_3 = '${fparse detrapping_prefactor}'
K_3 = '${units 2.5e26 traps/m^3/dpa}'
nmax_3 = '${units 6.0e25 1/m^3}'
Ea_3 = '${units 0.30 eV}'
A_3 = '${fparse A0 * exp(-Ea_3/(kb_eV*temperature_damage))}'
S_3 = '${fparse phi*K_3}'
r_3 = '${fparse phi*K_3/nmax_3 + A_3}'
trap_density_01dpa_3 = '${units ${fparse -S_3/r_3 * exp(-r_3*damage_duration) + S_3/r_3} at/m^3 -> at/mum^3}'
trapping_site_fraction_3 = '${fparse trap_density_01dpa_3 / tungsten_density}'
trap_density_01dpa_3_inf = '${units 2.6e25 at/m^3 -> at/mum^3}'
trapping_site_fraction_3_inf = '${fparse trap_density_01dpa_3_inf / tungsten_density}'
trap_concentration_reference_3 = '${fparse tungsten_density * 1e-8}'
# Trap 4 parameters
trapping_energy_4 = '${fparse trapping_energy}'
trapping_prefactor_4 = '${fparse trapping_prefactor}'
detrapping_energy_4 = '${fparse ${units 1.85 eV -> J} / kb}'
detrapping_prefactor_4 = '${fparse detrapping_prefactor}'
K_4 = '${units 5.0e26 traps/m^3/dpa}'
nmax_4 = '${units 4.7e25 1/m^3}'
Ea_4 = '${units 0.30 eV}'
A_4 = '${fparse A0 * exp(-Ea_4/(kb_eV*temperature_damage))}'
S_4 = '${fparse phi*K_4}'
r_4 = '${fparse phi*K_4/nmax_4 + A_4}'
trap_density_01dpa_4 = '${units ${fparse -S_4/r_4 * exp(-r_4*damage_duration) + S_4/r_4} at/m^3 -> at/mum^3}'
trapping_site_fraction_4 = '${fparse trap_density_01dpa_4 / tungsten_density}'
trap_density_01dpa_4_inf = '${units 3.6e25 at/m^3 -> at/mum^3}'
trapping_site_fraction_4_inf = '${fparse trap_density_01dpa_4_inf / tungsten_density}'
trap_concentration_reference_4 = '${fparse tungsten_density * 1e-8}'
# Trap 5 parameters
trapping_energy_5 = '${fparse trapping_energy}'
trapping_prefactor_5 = '${fparse trapping_prefactor}'
detrapping_energy_5 = '${fparse ${units 2.05 eV -> J} / kb}'
detrapping_prefactor_5 = '${fparse detrapping_prefactor}'
K_5 = '${units 1.0e26 traps/m^3/dpa}'
nmax_5 = '${units 2.0e25 1/m^3}'
S_5 = '${fparse phi*K_5}'
r_5 = '${fparse phi*K_5/nmax_5}'
trap_density_01dpa_5 = '${units ${fparse -S_5/r_5 * exp(-r_5*damage_duration) + S_5/r_5} at/m^3 -> at/mum^3}'
trapping_site_fraction_5 = '${fparse trap_density_01dpa_5 / tungsten_density}'
trap_density_01dpa_5_inf = '${units 1.1e25 at/m^3 -> at/mum^3}'
trapping_site_fraction_5_inf = '${fparse trap_density_01dpa_5_inf / tungsten_density}'
trap_concentration_reference_5 = '${fparse tungsten_density * 1e-8}'
# Trap intrinsic parameters
trapping_energy_intrinsic = '${fparse trapping_energy}'
trapping_prefactor_intrinsic = '${fparse trapping_prefactor}'
detrapping_energy_intrinsic = '${fparse ${units 1.04 eV -> J} / kb}'
detrapping_prefactor_intrinsic = '${fparse detrapping_prefactor}'
trap_density_01dpa_intrinsic = '${units 2.4e22 at/m^3 -> at/mum^3}'
trapping_site_fraction_intrinsic = '${fparse trap_density_01dpa_intrinsic / tungsten_density}'
trap_concentration_reference_intrinsic = '${fparse tungsten_density * 1e-11}'
mobile_concentration_reference = '${fparse tungsten_density * 1e-12}'
length_reference = ${units 1 mum}
time_reference = ${units 1 s}
# Dimensionless trap rate groups: k_t_hat = alpha_t * C_m_ref / N * t_ref
# (from nondimensionalizing dC_t/dt = alpha_t/N * C_m * (N*Ct0 - C_t))
dimensionless_trapping_rate_coefficient_intrinsic = '${fparse trapping_prefactor_intrinsic * time_reference * mobile_concentration_reference / tungsten_density}'
dimensionless_trapping_rate_coefficient_1 = '${fparse trapping_prefactor_1 * time_reference * mobile_concentration_reference / tungsten_density}'
dimensionless_trapping_rate_coefficient_2 = '${fparse trapping_prefactor_2 * time_reference * mobile_concentration_reference / tungsten_density}'
dimensionless_trapping_rate_coefficient_3 = '${fparse trapping_prefactor_3 * time_reference * mobile_concentration_reference / tungsten_density}'
dimensionless_trapping_rate_coefficient_4 = '${fparse trapping_prefactor_4 * time_reference * mobile_concentration_reference / tungsten_density}'
dimensionless_trapping_rate_coefficient_5 = '${fparse trapping_prefactor_5 * time_reference * mobile_concentration_reference / tungsten_density}'
dimensionless_release_rate_coefficient_intrinsic = '${fparse detrapping_prefactor_intrinsic * time_reference}'
dimensionless_release_rate_coefficient_1 = '${fparse detrapping_prefactor_1 * time_reference}'
dimensionless_release_rate_coefficient_2 = '${fparse detrapping_prefactor_2 * time_reference}'
dimensionless_release_rate_coefficient_3 = '${fparse detrapping_prefactor_3 * time_reference}'
dimensionless_release_rate_coefficient_4 = '${fparse detrapping_prefactor_4 * time_reference}'
dimensionless_release_rate_coefficient_5 = '${fparse detrapping_prefactor_5 * time_reference}'
# Dimensionless geometry and time quantities used by val-2f-dimensionless.
sigma_hat = '${fparse sigma / length_reference}'
R_p_hat = '${fparse R_p / length_reference}'
depth_center_hat = '${fparse depth_center / length_reference}'
depth_width_hat = '${fparse depth_width / length_reference}'
sample_thickness_hat = '${fparse sample_thickness / length_reference}'
dx1_hat = '${fparse dx1 / length_reference}'
dx2_hat = '${fparse dx2 / length_reference}'
dx3_hat = '${fparse dx3 / length_reference}'
dx4_hat = '${fparse dx4 / length_reference}'
dx5_hat = '${fparse sample_thickness_hat - dx1_hat - dx2_hat - dx3_hat - dx4_hat}'
charge_time_hat = '${fparse charge_time / time_reference}'
cooldown_duration_hat = '${fparse cooldown_duration / time_reference}'
desorption_duration_hat = '${fparse desorption_duration / time_reference}'
endtime_hat = '${fparse charge_time_hat + cooldown_duration_hat + desorption_duration_hat}'
desorption_heating_rate_hat = '${fparse desorption_heating_rate * time_reference}'
# Dimensionless mobile transport quantities.
diffusion_W_preexponential_hat = '${fparse diffusion_W_preexponential * time_reference / (length_reference^2)}'
surface_flux_hat = '${fparse flux * time_reference / (mobile_concentration_reference * length_reference)}'
recombination_coefficient_hat = '${fparse recombination_coefficient * mobile_concentration_reference * time_reference / length_reference}'
sieverts_boundary_hat = '${fparse 1e-10 / mobile_concentration_reference}'
(test/tests/val-2f/val-2f_trapping_intrinsic.i)
# This input file adds the intrinsic trapping sites for validation case val-2f.
# It is included in val-2f.i
[Bounds]
[trapped_intrinsic_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_intrinsic
bound_type = lower
bound_value = 0
[]
[]
[Physics]
[SpeciesTrapping]
[trapping_intrinsic]
species = 'trapped_intrinsic'
species_scaling_factors = '1'
species_initial_concentrations = '0'
mobile = 'deuterium_concentration_W'
dimensionless_trapping_rate_coefficient = '${dimensionless_trapping_rate_coefficient_intrinsic}'
trapping_energy = '${trapping_energy_intrinsic}'
N = ${tungsten_density}
Ct0 = '${trapping_site_fraction_intrinsic}'
trap_concentration_reference = '${trap_concentration_reference_intrinsic}'
mobile_concentration_reference = ${mobile_concentration_reference}
dimensionless_release_rate_coefficient = '${dimensionless_release_rate_coefficient_intrinsic}'
detrapping_energy = '${detrapping_energy_intrinsic}'
temperature = 'temperature'
dimensionless_species = true
[]
[]
[]
[Postprocessors]
[integral_trapped_concentration_intrinsic]
type = ElementIntegralVariablePostprocessor
variable = trapped_intrinsic
outputs = none
[]
[scaled_trapped_deuterium_intrinsic]
type = ScalePostprocessor
scaling_factor = '${fparse trap_concentration_reference_intrinsic * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_intrinsic
[]
[max_scaled_trapped_deuterium_intrinsic]
type = TimeExtremeValue
postprocessor = scaled_trapped_deuterium_intrinsic
value_type = max
outputs = 'console'
[]
[]
(test/tests/val-2f/val-2f_trapping_1.i)
# This input file adds the trapping sites 1 for validation case val-2f.
# It is included in val-2f.i
[Bounds]
[trapped_1_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_1
bound_type = lower
bound_value = 0
[]
[]
[Functions]
[trap_distribution_function_1]
type = ParsedFunction
expression = '${trapping_site_fraction_1} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[trap_distribution_function_1_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_1_inf} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[]
[Physics]
[SpeciesTrapping]
[trapping_1]
species = 'trapped_1'
species_scaling_factors = '1'
species_initial_concentrations = '0'
mobile = 'deuterium_concentration_W'
dimensionless_trapping_rate_coefficient = '${dimensionless_trapping_rate_coefficient_1}'
trapping_energy = '${trapping_energy_1}'
N = ${tungsten_density}
Ct0 = 'trap_distribution_function_1'
trap_concentration_reference = '${trap_concentration_reference_1}'
mobile_concentration_reference = ${mobile_concentration_reference}
dimensionless_release_rate_coefficient = '${dimensionless_release_rate_coefficient_1}'
detrapping_energy = '${detrapping_energy_1}'
temperature = 'temperature'
dimensionless_species = true
[]
[]
[]
[Postprocessors]
[integral_trapped_concentration_1]
type = ElementIntegralVariablePostprocessor
variable = trapped_1
outputs = none
[]
[scaled_trapped_deuterium_1]
type = ScalePostprocessor
scaling_factor = '${fparse trap_concentration_reference_1 * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_1
[]
[]
(test/tests/val-2f/val-2f_trapping_2.i)
# This input file adds the trapping sites 2 for validation case val-2f.
# It is included in val-2f.i
[Bounds]
[trapped_2_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_2
bound_type = lower
bound_value = 0
[]
[]
[Functions]
[trap_distribution_function_2]
type = ParsedFunction
expression = '${trapping_site_fraction_2} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[trap_distribution_function_2_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_2_inf} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[]
[Physics]
[SpeciesTrapping]
[trapping_2]
species = 'trapped_2'
species_scaling_factors = '1'
species_initial_concentrations = '0'
mobile = 'deuterium_concentration_W'
dimensionless_trapping_rate_coefficient = '${dimensionless_trapping_rate_coefficient_2}'
trapping_energy = '${trapping_energy_2}'
N = ${tungsten_density}
Ct0 = 'trap_distribution_function_2'
trap_concentration_reference = '${trap_concentration_reference_2}'
mobile_concentration_reference = ${mobile_concentration_reference}
dimensionless_release_rate_coefficient = '${dimensionless_release_rate_coefficient_2}'
detrapping_energy = '${detrapping_energy_2}'
temperature = 'temperature'
dimensionless_species = true
[]
[]
[]
[Postprocessors]
[integral_trapped_concentration_2]
type = ElementIntegralVariablePostprocessor
variable = trapped_2
outputs = none
[]
[scaled_trapped_deuterium_2]
type = ScalePostprocessor
scaling_factor = '${fparse trap_concentration_reference_2 * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_2
[]
[]
(test/tests/val-2f/val-2f_trapping_3.i)
# This input file adds the trapping sites 3 for validation case val-2f.
# It is included in val-2f.i
[Bounds]
[trapped_3_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_3
bound_type = lower
bound_value = 0
[]
[]
[Functions]
[trap_distribution_function_3]
type = ParsedFunction
expression = '${trapping_site_fraction_3} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[trap_distribution_function_3_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_3_inf} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[]
[Physics]
[SpeciesTrapping]
[trapping_3]
species = 'trapped_3'
species_scaling_factors = '1'
species_initial_concentrations = '0'
mobile = 'deuterium_concentration_W'
dimensionless_trapping_rate_coefficient = '${dimensionless_trapping_rate_coefficient_3}'
trapping_energy = '${trapping_energy_3}'
N = ${tungsten_density}
Ct0 = 'trap_distribution_function_3'
trap_concentration_reference = '${trap_concentration_reference_3}'
mobile_concentration_reference = ${mobile_concentration_reference}
dimensionless_release_rate_coefficient = '${dimensionless_release_rate_coefficient_3}'
detrapping_energy = '${detrapping_energy_3}'
temperature = 'temperature'
dimensionless_species = true
[]
[]
[]
[Postprocessors]
[integral_trapped_concentration_3]
type = ElementIntegralVariablePostprocessor
variable = trapped_3
outputs = none
[]
[scaled_trapped_deuterium_3]
type = ScalePostprocessor
scaling_factor = '${fparse trap_concentration_reference_3 * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_3
[]
[]
(test/tests/val-2f/val-2f_trapping_4.i)
# This input file adds the trapping sites 4 for validation case val-2f.
# It is included in val-2f.i
[Bounds]
[trapped_4_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_4
bound_type = lower
bound_value = 0
[]
[]
[Functions]
[trap_distribution_function_4]
type = ParsedFunction
expression = '${trapping_site_fraction_4} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[trap_distribution_function_4_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_4_inf} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[]
[Physics]
[SpeciesTrapping]
[trapping_4]
species = 'trapped_4'
species_scaling_factors = '1'
species_initial_concentrations = '0'
mobile = 'deuterium_concentration_W'
dimensionless_trapping_rate_coefficient = '${dimensionless_trapping_rate_coefficient_4}'
trapping_energy = '${trapping_energy_4}'
N = ${tungsten_density}
Ct0 = 'trap_distribution_function_4'
trap_concentration_reference = '${trap_concentration_reference_4}'
mobile_concentration_reference = ${mobile_concentration_reference}
dimensionless_release_rate_coefficient = '${dimensionless_release_rate_coefficient_4}'
detrapping_energy = '${detrapping_energy_4}'
temperature = 'temperature'
dimensionless_species = true
[]
[]
[]
[Postprocessors]
[integral_trapped_concentration_4]
type = ElementIntegralVariablePostprocessor
variable = trapped_4
outputs = none
[]
[scaled_trapped_deuterium_4]
type = ScalePostprocessor
scaling_factor = '${fparse trap_concentration_reference_4 * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_4
[]
[]
(test/tests/val-2f/val-2f_trapping_5.i)
# This input file adds the trapping sites 5 for validation case val-2f.
# It is included in val-2f.i
[Bounds]
[trapped_5_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_5
bound_type = lower
bound_value = 0
[]
[]
[Functions]
[trap_distribution_function_5]
type = ParsedFunction
expression = '${trapping_site_fraction_5} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[trap_distribution_function_5_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_5_inf} / (1 + exp((x - ${depth_center_hat}) / ${depth_width_hat}))'
[]
[]
[Physics]
[SpeciesTrapping]
[trapping_5]
species = 'trapped_5'
species_scaling_factors = '1'
species_initial_concentrations = '0'
mobile = 'deuterium_concentration_W'
dimensionless_trapping_rate_coefficient = '${dimensionless_trapping_rate_coefficient_5}'
trapping_energy = '${trapping_energy_5}'
N = ${tungsten_density}
Ct0 = 'trap_distribution_function_5'
trap_concentration_reference = '${trap_concentration_reference_5}'
mobile_concentration_reference = ${mobile_concentration_reference}
dimensionless_release_rate_coefficient = '${dimensionless_release_rate_coefficient_5}'
detrapping_energy = '${detrapping_energy_5}'
temperature = 'temperature'
dimensionless_species = true
[]
[]
[]
[Postprocessors]
[integral_trapped_concentration_5]
type = ElementIntegralVariablePostprocessor
variable = trapped_5
outputs = none
[]
[scaled_trapped_deuterium_5]
type = ScalePostprocessor
scaling_factor = '${fparse trap_concentration_reference_5 * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_5
[]
[]
(test/tests/val-2f/val-2f.i)
# Validation Problem #2f
# Self-damaged Tungsten Effects on Deuterium Transport extended from an original model from
# Dark, J., Delaporte-Mathurin, R., Schwarz-Selinger, T., Hodille, E. A., Mougenot, J.,
# Charles, Y., & Grisolia, C. (2024). Modelling neutron damage effects on tritium transport
# in tungsten. Nuclear Fusion, 64(8), 086026.
!include parameters_val-2f.params
!include val-2f_trapping_intrinsic.i
!include val-2f_trapping_5.i
!include val-2f_trapping_4.i
!include val-2f_trapping_3.i
!include val-2f_trapping_2.i
!include val-2f_trapping_1.i
[Mesh]
[cartesian_mesh]
type = CartesianMeshGenerator
dim = 1
dx = '${dx1_hat} ${dx2_hat} ${dx3_hat} ${dx4_hat} ${dx5_hat}'
ix = '${ix1} ${ix2} ${ix3} ${ix4} ${ix5}'
subdomain_id = '0 0 0 0 0'
[]
[]
[Variables]
[deuterium_concentration_W]
[]
[]
[AuxVariables]
active = 'bounds_dummy temperature'
[bounds_dummy]
order = FIRST
family = LAGRANGE
[]
[temperature]
initial_condition = ${temperature_initial}
[]
[]
[Bounds]
[deuterium_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = deuterium_concentration_W
bound_type = lower
bound_value = 0
[]
[]
[Kernels]
[time_W]
type = TimeDerivative
variable = deuterium_concentration_W
[]
[diffusion_W]
type = ADMatDiffusion
variable = deuterium_concentration_W
diffusivity = diffusivity_W
[]
[source_deuterium]
type = BodyForce
variable = deuterium_concentration_W
function = source_deuterium
[]
[]
[AuxKernels]
[temperature_aux]
type = FunctionAux
variable = temperature
function = temperature_bc_func
execute_on = 'INITIAL LINEAR'
[]
[]
[BCs]
active = 'left_recombination_flux right_recombination_flux'
# Kinetic boundary conditions
[left_recombination_flux]
type = ADMatNeumannBC
variable = deuterium_concentration_W
boundary = left
value = 1
boundary_material = flux_recombination_surface
[]
[right_recombination_flux]
type = ADMatNeumannBC
variable = deuterium_concentration_W
boundary = right
value = 1
boundary_material = flux_recombination_surface
[]
# Sieverts boundary conditions
[left_concentration_sieverts]
type = ADDirichletBC
value = '${sieverts_boundary_hat}'
boundary = left
variable = deuterium_concentration_W
[]
[right_concentration_sieverts]
type = ADDirichletBC
value = '${sieverts_boundary_hat}'
boundary = right
variable = deuterium_concentration_W
[]
[]
[Functions]
[temperature_bc_func]
type = ParsedFunction
expression = 'if(t<${charge_time_hat}, ${temperature_initial},
if(t<${fparse charge_time_hat + cooldown_duration_hat}, ${temperature_cooldown},
${temperature_desorption_min}+${desorption_heating_rate_hat}*(t-${fparse charge_time_hat + cooldown_duration_hat})))'
[]
[source_distribution]
type = ParsedFunction
expression = '1 / (${sigma_hat} * sqrt(2 * pi)) * exp(-0.5 * ((x - ${R_p_hat}) / ${sigma_hat}) ^ 2)'
[]
[surface_flux_func]
type = ParsedFunction
expression = 'if(t<${charge_time_hat}, ${surface_flux_hat}, 0)'
[]
[source_deuterium]
type = ParsedFunction
symbol_names = 'source_distribution surface_flux_func'
symbol_values = 'source_distribution surface_flux_func'
expression = 'source_distribution * surface_flux_func'
[]
[max_dt_size_function]
type = ParsedFunction
expression = 'if(t<${fparse 5 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 8 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 12 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 20 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 35 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 450 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 5000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse 11000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 13000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse (charge_time + cooldown_duration + 4500) / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 313000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 315000 / time_reference}, ${fparse 1e1 / time_reference}, ${fparse 1e3 / time_reference}))))))))))))'
[]
[max_dt_size_function_inf]
type = ParsedFunction
expression = 'if(t<${fparse 5 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 8 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 12 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 20 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 35 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 450 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 5000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse 11000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 13000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse (charge_time + cooldown_duration + 4500) / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 315000 / time_reference}, ${fparse 1e1 / time_reference}, ${fparse 1e3 / time_reference})))))))))))'
[]
[max_dt_size_function_coarse]
type = ParsedFunction
expression = 'if(t<${fparse 1e-1 / time_reference}, ${fparse 1e4 / time_reference}, ${fparse 1e5 / time_reference})'
[]
[]
[Materials]
active = 'diffusivity_W_func diffusivity_nonAD recombination_rate_surface flux_recombination_surface'
[diffusivity_W_func]
type = ADDerivativeParsedMaterial
property_name = 'diffusivity_W'
functor_names = 'temperature_bc_func'
functor_symbols = 'temperature'
expression = '${diffusion_W_preexponential_hat} * exp(- ${diffusion_W_energy} / ${kb_eV} / temperature)'
[]
[diffusivity_nonAD]
type = MaterialADConverter
ad_props_in = 'diffusivity_W'
reg_props_out = 'diffusivity_W_nonAD'
[]
[recombination_rate_surface]
type = ADDerivativeParsedMaterial
property_name = 'Kr'
functor_names = 'temperature_bc_func'
functor_symbols = 'temperature'
expression = '${recombination_coefficient_hat} * exp(- ${recombination_energy} / ${kb_eV} / temperature)'
[]
[flux_recombination_surface]
type = ADDerivativeParsedMaterial
coupled_variables = 'deuterium_concentration_W'
property_name = 'flux_recombination_surface'
material_property_names = 'Kr'
expression = '- 2 * Kr * deuterium_concentration_W ^ 2'
[]
[]
[Postprocessors]
active = 'integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration
scaled_mobile_deuterium flux_surface_left scaled_flux_surface_left
flux_surface_right scaled_flux_surface_right temperature_pps diffusion_W_hat diffusion_W
max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1
integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3
integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5
integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic
spatial_max_mobile_d2 spatial_max_trapped_1 spatial_max_trapped_2 spatial_max_trapped_3 spatial_max_trapped_4 spatial_max_trapped_5 spatial_max_trapped_intrinsic
max_mobile_d2 max_trapped_1 max_trapped_2 max_trapped_3 max_trapped_4 max_trapped_5 max_trapped_intrinsic max_scaled_flux_surface_left max_scaled_flux_surface_right
max_scaled_mobile_deuterium max_scaled_trapped_deuterium_intrinsic'
[integral_source_deuterium]
type = FunctionElementIntegral
function = source_deuterium
outputs = none
[]
[scaled_implanted_deuterium]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = integral_source_deuterium
[]
[integral_deuterium_concentration]
type = ElementIntegralVariablePostprocessor
variable = deuterium_concentration_W
outputs = none
[]
[scaled_mobile_deuterium]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_deuterium_concentration
[]
[max_scaled_mobile_deuterium]
type = TimeExtremeValue
postprocessor = scaled_mobile_deuterium
value_type = max
outputs = 'console'
[]
[flux_surface_left]
type = ADSideAverageMaterialProperty
boundary = 'left'
property = flux_recombination_surface
outputs = none
[]
[scaled_flux_surface_left]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_left
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_left]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_left
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[flux_surface_left_sieverts]
type = SideDiffusiveFluxAverage
variable = deuterium_concentration_W
boundary = 'left'
diffusivity = 'diffusivity_W_nonAD'
outputs = none
[]
[scaled_flux_surface_left_sieverts]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_left_sieverts
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_left_sieverts]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_left_sieverts
value_type = max
[]
[flux_surface_right]
type = ADSideAverageMaterialProperty
boundary = 'right'
property = flux_recombination_surface
outputs = none
[]
[scaled_flux_surface_right]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_right
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_right]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_right
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[flux_surface_right_sieverts]
type = SideDiffusiveFluxAverage
variable = deuterium_concentration_W
boundary = 'right'
diffusivity = 'diffusivity_W_nonAD'
outputs = none
[]
[scaled_flux_surface_right_sieverts]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_right_sieverts
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_right_sieverts]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_right_sieverts
value_type = max
[]
[temperature_pps]
type = ElementAverageValue
variable = temperature
execute_on = 'initial timestep_end'
[]
[diffusion_W_hat]
type = ElementAverageMaterialProperty
mat_prop = diffusivity_W_nonAD
outputs = none
[]
[diffusion_W]
type = ScalePostprocessor
scaling_factor = '${fparse length_reference ^ 2 / time_reference}'
value = diffusion_W_hat
outputs = none
[]
[max_time_step_size]
type = FunctionValuePostprocessor
function = max_dt_size_function
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[max_time_step_size_coarse]
type = FunctionValuePostprocessor
function = max_dt_size_function_coarse
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[spatial_max_mobile_d2]
type = NodalExtremeValue
value_type = 'max'
variable = deuterium_concentration_W
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_1]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_1
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_2]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_3]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_3
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_4]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_4
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_5]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_5
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_intrinsic]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_intrinsic
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_mobile_d2]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_mobile_d2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_1]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_1
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_2]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_3]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_3
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_4]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_4
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_5]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_5
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_intrinsic]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_intrinsic
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = 'Newton'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -snes_type'
petsc_options_value = 'lu mumps vinewtonrsls'
end_time = ${endtime_hat}
line_search = 'none'
nl_rel_tol = 1e-8
nl_abs_tol = 4e-5
nl_max_its = 34
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse dt_init / time_reference}'
growth_factor = 1.1
timestep_limiting_postprocessor = max_time_step_size
[]
[Predictor]
type = SimplePredictor
scale = 1.0
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Outputs]
[csv]
type = CSV
file_base = 'val-2f_out'
[]
[csv_temperature_history]
type = CSV
show = 'temperature_pps'
file_base = 'val-2f_temperature_out'
[]
[exodus]
type = Exodus
output_material_properties = true
time_step_interval = 200
file_base = 'val-2f_out'
[]
[]
(test/tests/val-2f/val-2f.i)
# Validation Problem #2f
# Self-damaged Tungsten Effects on Deuterium Transport extended from an original model from
# Dark, J., Delaporte-Mathurin, R., Schwarz-Selinger, T., Hodille, E. A., Mougenot, J.,
# Charles, Y., & Grisolia, C. (2024). Modelling neutron damage effects on tritium transport
# in tungsten. Nuclear Fusion, 64(8), 086026.
!include parameters_val-2f.params
!include val-2f_trapping_intrinsic.i
!include val-2f_trapping_5.i
!include val-2f_trapping_4.i
!include val-2f_trapping_3.i
!include val-2f_trapping_2.i
!include val-2f_trapping_1.i
[Mesh]
[cartesian_mesh]
type = CartesianMeshGenerator
dim = 1
dx = '${dx1_hat} ${dx2_hat} ${dx3_hat} ${dx4_hat} ${dx5_hat}'
ix = '${ix1} ${ix2} ${ix3} ${ix4} ${ix5}'
subdomain_id = '0 0 0 0 0'
[]
[]
[Variables]
[deuterium_concentration_W]
[]
[]
[AuxVariables]
active = 'bounds_dummy temperature'
[bounds_dummy]
order = FIRST
family = LAGRANGE
[]
[temperature]
initial_condition = ${temperature_initial}
[]
[]
[Bounds]
[deuterium_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = deuterium_concentration_W
bound_type = lower
bound_value = 0
[]
[]
[Kernels]
[time_W]
type = TimeDerivative
variable = deuterium_concentration_W
[]
[diffusion_W]
type = ADMatDiffusion
variable = deuterium_concentration_W
diffusivity = diffusivity_W
[]
[source_deuterium]
type = BodyForce
variable = deuterium_concentration_W
function = source_deuterium
[]
[]
[AuxKernels]
[temperature_aux]
type = FunctionAux
variable = temperature
function = temperature_bc_func
execute_on = 'INITIAL LINEAR'
[]
[]
[BCs]
active = 'left_recombination_flux right_recombination_flux'
# Kinetic boundary conditions
[left_recombination_flux]
type = ADMatNeumannBC
variable = deuterium_concentration_W
boundary = left
value = 1
boundary_material = flux_recombination_surface
[]
[right_recombination_flux]
type = ADMatNeumannBC
variable = deuterium_concentration_W
boundary = right
value = 1
boundary_material = flux_recombination_surface
[]
# Sieverts boundary conditions
[left_concentration_sieverts]
type = ADDirichletBC
value = '${sieverts_boundary_hat}'
boundary = left
variable = deuterium_concentration_W
[]
[right_concentration_sieverts]
type = ADDirichletBC
value = '${sieverts_boundary_hat}'
boundary = right
variable = deuterium_concentration_W
[]
[]
[Functions]
[temperature_bc_func]
type = ParsedFunction
expression = 'if(t<${charge_time_hat}, ${temperature_initial},
if(t<${fparse charge_time_hat + cooldown_duration_hat}, ${temperature_cooldown},
${temperature_desorption_min}+${desorption_heating_rate_hat}*(t-${fparse charge_time_hat + cooldown_duration_hat})))'
[]
[source_distribution]
type = ParsedFunction
expression = '1 / (${sigma_hat} * sqrt(2 * pi)) * exp(-0.5 * ((x - ${R_p_hat}) / ${sigma_hat}) ^ 2)'
[]
[surface_flux_func]
type = ParsedFunction
expression = 'if(t<${charge_time_hat}, ${surface_flux_hat}, 0)'
[]
[source_deuterium]
type = ParsedFunction
symbol_names = 'source_distribution surface_flux_func'
symbol_values = 'source_distribution surface_flux_func'
expression = 'source_distribution * surface_flux_func'
[]
[max_dt_size_function]
type = ParsedFunction
expression = 'if(t<${fparse 5 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 8 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 12 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 20 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 35 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 450 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 5000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse 11000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 13000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse (charge_time + cooldown_duration + 4500) / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 313000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 315000 / time_reference}, ${fparse 1e1 / time_reference}, ${fparse 1e3 / time_reference}))))))))))))'
[]
[max_dt_size_function_inf]
type = ParsedFunction
expression = 'if(t<${fparse 5 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 8 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 12 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 20 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 35 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 450 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 5000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse 11000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 13000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse (charge_time + cooldown_duration + 4500) / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 315000 / time_reference}, ${fparse 1e1 / time_reference}, ${fparse 1e3 / time_reference})))))))))))'
[]
[max_dt_size_function_coarse]
type = ParsedFunction
expression = 'if(t<${fparse 1e-1 / time_reference}, ${fparse 1e4 / time_reference}, ${fparse 1e5 / time_reference})'
[]
[]
[Materials]
active = 'diffusivity_W_func diffusivity_nonAD recombination_rate_surface flux_recombination_surface'
[diffusivity_W_func]
type = ADDerivativeParsedMaterial
property_name = 'diffusivity_W'
functor_names = 'temperature_bc_func'
functor_symbols = 'temperature'
expression = '${diffusion_W_preexponential_hat} * exp(- ${diffusion_W_energy} / ${kb_eV} / temperature)'
[]
[diffusivity_nonAD]
type = MaterialADConverter
ad_props_in = 'diffusivity_W'
reg_props_out = 'diffusivity_W_nonAD'
[]
[recombination_rate_surface]
type = ADDerivativeParsedMaterial
property_name = 'Kr'
functor_names = 'temperature_bc_func'
functor_symbols = 'temperature'
expression = '${recombination_coefficient_hat} * exp(- ${recombination_energy} / ${kb_eV} / temperature)'
[]
[flux_recombination_surface]
type = ADDerivativeParsedMaterial
coupled_variables = 'deuterium_concentration_W'
property_name = 'flux_recombination_surface'
material_property_names = 'Kr'
expression = '- 2 * Kr * deuterium_concentration_W ^ 2'
[]
[]
[Postprocessors]
active = 'integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration
scaled_mobile_deuterium flux_surface_left scaled_flux_surface_left
flux_surface_right scaled_flux_surface_right temperature_pps diffusion_W_hat diffusion_W
max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1
integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3
integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5
integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic
spatial_max_mobile_d2 spatial_max_trapped_1 spatial_max_trapped_2 spatial_max_trapped_3 spatial_max_trapped_4 spatial_max_trapped_5 spatial_max_trapped_intrinsic
max_mobile_d2 max_trapped_1 max_trapped_2 max_trapped_3 max_trapped_4 max_trapped_5 max_trapped_intrinsic max_scaled_flux_surface_left max_scaled_flux_surface_right
max_scaled_mobile_deuterium max_scaled_trapped_deuterium_intrinsic'
[integral_source_deuterium]
type = FunctionElementIntegral
function = source_deuterium
outputs = none
[]
[scaled_implanted_deuterium]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = integral_source_deuterium
[]
[integral_deuterium_concentration]
type = ElementIntegralVariablePostprocessor
variable = deuterium_concentration_W
outputs = none
[]
[scaled_mobile_deuterium]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_deuterium_concentration
[]
[max_scaled_mobile_deuterium]
type = TimeExtremeValue
postprocessor = scaled_mobile_deuterium
value_type = max
outputs = 'console'
[]
[flux_surface_left]
type = ADSideAverageMaterialProperty
boundary = 'left'
property = flux_recombination_surface
outputs = none
[]
[scaled_flux_surface_left]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_left
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_left]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_left
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[flux_surface_left_sieverts]
type = SideDiffusiveFluxAverage
variable = deuterium_concentration_W
boundary = 'left'
diffusivity = 'diffusivity_W_nonAD'
outputs = none
[]
[scaled_flux_surface_left_sieverts]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_left_sieverts
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_left_sieverts]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_left_sieverts
value_type = max
[]
[flux_surface_right]
type = ADSideAverageMaterialProperty
boundary = 'right'
property = flux_recombination_surface
outputs = none
[]
[scaled_flux_surface_right]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_right
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_right]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_right
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[flux_surface_right_sieverts]
type = SideDiffusiveFluxAverage
variable = deuterium_concentration_W
boundary = 'right'
diffusivity = 'diffusivity_W_nonAD'
outputs = none
[]
[scaled_flux_surface_right_sieverts]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_right_sieverts
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_right_sieverts]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_right_sieverts
value_type = max
[]
[temperature_pps]
type = ElementAverageValue
variable = temperature
execute_on = 'initial timestep_end'
[]
[diffusion_W_hat]
type = ElementAverageMaterialProperty
mat_prop = diffusivity_W_nonAD
outputs = none
[]
[diffusion_W]
type = ScalePostprocessor
scaling_factor = '${fparse length_reference ^ 2 / time_reference}'
value = diffusion_W_hat
outputs = none
[]
[max_time_step_size]
type = FunctionValuePostprocessor
function = max_dt_size_function
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[max_time_step_size_coarse]
type = FunctionValuePostprocessor
function = max_dt_size_function_coarse
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[spatial_max_mobile_d2]
type = NodalExtremeValue
value_type = 'max'
variable = deuterium_concentration_W
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_1]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_1
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_2]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_3]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_3
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_4]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_4
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_5]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_5
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_intrinsic]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_intrinsic
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_mobile_d2]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_mobile_d2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_1]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_1
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_2]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_3]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_3
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_4]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_4
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_5]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_5
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_intrinsic]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_intrinsic
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = 'Newton'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -snes_type'
petsc_options_value = 'lu mumps vinewtonrsls'
end_time = ${endtime_hat}
line_search = 'none'
nl_rel_tol = 1e-8
nl_abs_tol = 4e-5
nl_max_its = 34
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse dt_init / time_reference}'
growth_factor = 1.1
timestep_limiting_postprocessor = max_time_step_size
[]
[Predictor]
type = SimplePredictor
scale = 1.0
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Outputs]
[csv]
type = CSV
file_base = 'val-2f_out'
[]
[csv_temperature_history]
type = CSV
show = 'temperature_pps'
file_base = 'val-2f_temperature_out'
[]
[exodus]
type = Exodus
output_material_properties = true
time_step_interval = 200
file_base = 'val-2f_out'
[]
[]
(test/tests/val-2f/tests)
[Tests]
design = 'MatNeumannBC.md CoefCoupledTimeDerivative.md TimeDerivativeNodalKernel.md TrappingNodalKernel.md ReleasingNodalKernel.md'
issues = '#12'
validation = 'val-2f.md'
[val-2f_light_csv]
type = CSVDiff
input = val-2f.i
cli_args = "sample_thickness='${units 0.9e-4 m -> mum}'
ix1=8 ix2=8 ix3=8 ix4=8 ix5=8
dt_init=2e-3
charge_time='${units 1e-1 s}' cooldown_duration='${units 1e-1 s}'
Executioner/TimeStepper/timestep_limiting_postprocessor=max_time_step_size_coarse
Outputs/csv/file_base=val-2f_light_out
Outputs/exodus/file_base=val-2f_light_out
Outputs/csv/start_time=0.2
Executioner/TimeStepper/growth_factor=1.05
Executioner/scheme=implicit-euler"
csvdiff = val-2f_light_out.csv
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport and generate CSV data output with a short runtime and coarse mesh testing.'
rel_err = 1e-3
override_columns = 'scaled_flux_surface_right scaled_flux_surface_left scaled_mobile_deuterium'
override_abs_zero = '1e4 1e6 1e7'
override_rel_err = '1e-3 1e-3 1e-3'
# Some relevant scales when considering diffs
#
# scaled_flux_surface_right O(1e12)
# scaled_flux_surface_left O(1e19)
# scaled_mobile_deuterium O(1e13)
# scaled_trapped_deuterium_intrinsic O(1e14)
[]
[val-2f_light_exodus]
type = Exodiff
input = val-2f.i
exodiff = val-2f_light_out.e
prereq = val-2f_light_csv
should_execute = false # this test relies on the output files from val-2f_light_csv, so it shouldn't be run twice
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport with a short runtime and coarse mesh testing.'
rel_err = 1e-3
abs_zero = 1e-5
[]
# Check test data consistency across operating systems and core counts for heavy tests using `Executioner/scheme=implicit-euler` due to its high stability
[val-2f_heavy_csv_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv/file_base=val-2f_heavy_desorption_out
Outputs/exodus/file_base=val-2f_heavy_desorption_out
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
csvdiff = val-2f_heavy_desorption_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output. Only desorption phase data is compared with every other timestep.'
max_time = 2000
[]
[val-2f_heavy_exodus_implicit-euler]
type = Exodiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv/file_base=val-2f_heavy_desorption_out
Outputs/exodus/file_base=val-2f_heavy_desorption_out
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
exodiff = val-2f_heavy_desorption_out.e
heavy = true
prereq = val-2f_heavy_csv_implicit-euler
should_execute = false # this test relies on the output files from val-2f_heavy_csv_implicit-euler, so it shouldn't be run twice
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability. Only desorption phase data is compared with every other timestep.'
[]
[val-2f_heavy_temperature_csv_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
prereq = val-2f_heavy_csv_implicit-euler
should_execute = false # this test relies on the output files from val-2f_heavy_csv_implicit-euler, so it shouldn't be run twice
csvdiff = val-2f_temperature_implicit_euler_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output about the temperature for documentation plots.'
max_time = 2000
[]
[val-2f_heavy_csv_inf_recombination_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "ix1='${fparse 500}' ix4='${fparse 500}'
AuxVariables/active='bounds_dummy temperature'
BCs/active='left_concentration_sieverts right_concentration_sieverts'
Physics/SpeciesTrapping/trapping_1/Ct0='trap_distribution_function_1_inf'
Physics/SpeciesTrapping/trapping_2/Ct0='trap_distribution_function_2_inf'
Physics/SpeciesTrapping/trapping_3/Ct0='trap_distribution_function_3_inf'
Physics/SpeciesTrapping/trapping_4/Ct0='trap_distribution_function_4_inf'
Physics/SpeciesTrapping/trapping_5/Ct0='trap_distribution_function_5_inf'
Materials/active='diffusivity_W_func diffusivity_nonAD'
Postprocessors/active='integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration scaled_mobile_deuterium flux_surface_left_sieverts scaled_flux_surface_left_sieverts flux_surface_right_sieverts scaled_flux_surface_right_sieverts temperature_pps max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1 integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3 integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5 integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic'
Postprocessors/max_time_step_size/function=max_dt_size_function_inf
Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/csv/file_base=val-2f_heavy_desorption_inf_recombination_out
Outputs/exodus/file_base=val-2f_heavy_desorption_inf_recombination_out
Executioner/abort_on_solve_fail=false
Executioner/nl_abs_tol=5e-4
Executioner/scheme=implicit-euler"
csvdiff = val-2f_heavy_desorption_inf_recombination_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output, for the infinite recombination case. Only desorption phase data is compared with every other timestep.'
max_time = 3600
[]
# Check the test data for the same heavy tests using 'Executioner/scheme=bdf2' for documentation. Since consistency has already been verified with `Executioner/scheme=implicit-euler` tests, the error tolerance can be relaxed.
[val-2f_heavy_csv_bdf2]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/sync_times='${fparse charge_time_hat + cooldown_duration_hat}'"
csvdiff = val-2f_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the bdf2 time integration scheme for accuracy and generate CSV data output.'
max_time = 2000
rel_err = 1e-4
override_columns = 'scaled_flux_surface_right scaled_flux_surface_left scaled_mobile_deuterium scaled_trapped_deuterium_1 scaled_trapped_deuterium_2'
override_abs_zero = '1e4 1e6 1e4 1e9 1e4'
override_rel_err = '5e-3 1e-3 9e-3 5e-3 5e-3'
[]
[val-2f_heavy_csv_inf_recombination_bdf2]
type = CSVDiff
input = val-2f.i
cli_args = "ix1='${fparse 500}' ix4='${fparse 500}'
AuxVariables/active='bounds_dummy temperature'
BCs/active='left_concentration_sieverts right_concentration_sieverts'
Physics/SpeciesTrapping/trapping_1/Ct0='trap_distribution_function_1_inf'
Physics/SpeciesTrapping/trapping_2/Ct0='trap_distribution_function_2_inf'
Physics/SpeciesTrapping/trapping_3/Ct0='trap_distribution_function_3_inf'
Physics/SpeciesTrapping/trapping_4/Ct0='trap_distribution_function_4_inf'
Physics/SpeciesTrapping/trapping_5/Ct0='trap_distribution_function_5_inf'
Materials/active='diffusivity_W_func diffusivity_nonAD'
Postprocessors/active='integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration scaled_mobile_deuterium flux_surface_left_sieverts scaled_flux_surface_left_sieverts flux_surface_right_sieverts scaled_flux_surface_right_sieverts temperature_pps diffusion_W_hat diffusion_W max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1 integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3 integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5 integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic spatial_max_mobile_d2 spatial_max_trapped_1 spatial_max_trapped_2 spatial_max_trapped_3 spatial_max_trapped_4 spatial_max_trapped_5 spatial_max_trapped_intrinsic max_mobile_d2 max_trapped_1 max_trapped_2 max_trapped_3 max_trapped_4 max_trapped_5 max_trapped_intrinsic max_scaled_flux_surface_left_sieverts max_scaled_flux_surface_right_sieverts max_scaled_mobile_deuterium max_scaled_trapped_deuterium_intrinsic'
Postprocessors/max_time_step_size/function=max_dt_size_function_inf
Outputs/csv/file_base='val-2f_out_inf_recombination'
Outputs/exodus/file_base='val-2f_out_inf_recombination'
Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/sync_times='${fparse charge_time_hat + cooldown_duration_hat}'
Executioner/abort_on_solve_fail=false
Executioner/nl_abs_tol=5e-4"
csvdiff = val-2f_out_inf_recombination.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the bdf2 time integration scheme for accuracy and generate CSV data output, for the infinite recombination case.'
max_time = 3600
rel_err = 1e-5
abs_zero = 1e4
override_columns = 'scaled_flux_surface_right_sieverts scaled_flux_surface_left_sieverts scaled_mobile_deuterium scaled_trapped_deuterium_1 scaled_trapped_deuterium_2 scaled_trapped_deuterium_3 scaled_trapped_deuterium_4 scaled_trapped_deuterium_5'
override_abs_zero = '1e6 1e9 1e6 1e9 1e9 1e9 1e9 1e10'
override_rel_err = '1e-3 3e-3 5e-3 1e-3 1e-3 3e-3 1e-3 5e-5'
[]
[val-2f_comparison]
type = RunCommand
command = 'python3 comparison_val-2f.py'
requirement = 'The system shall be able to generate comparison plots between simulated solutions and experimental data of validation case val-2f, modeling self-damaged tungsten effects on deuterium transport.'
required_python_packages = 'matplotlib numpy pandas scipy os'
[]
[]
(test/tests/val-2f/val-2f.i)
# Validation Problem #2f
# Self-damaged Tungsten Effects on Deuterium Transport extended from an original model from
# Dark, J., Delaporte-Mathurin, R., Schwarz-Selinger, T., Hodille, E. A., Mougenot, J.,
# Charles, Y., & Grisolia, C. (2024). Modelling neutron damage effects on tritium transport
# in tungsten. Nuclear Fusion, 64(8), 086026.
!include parameters_val-2f.params
!include val-2f_trapping_intrinsic.i
!include val-2f_trapping_5.i
!include val-2f_trapping_4.i
!include val-2f_trapping_3.i
!include val-2f_trapping_2.i
!include val-2f_trapping_1.i
[Mesh]
[cartesian_mesh]
type = CartesianMeshGenerator
dim = 1
dx = '${dx1_hat} ${dx2_hat} ${dx3_hat} ${dx4_hat} ${dx5_hat}'
ix = '${ix1} ${ix2} ${ix3} ${ix4} ${ix5}'
subdomain_id = '0 0 0 0 0'
[]
[]
[Variables]
[deuterium_concentration_W]
[]
[]
[AuxVariables]
active = 'bounds_dummy temperature'
[bounds_dummy]
order = FIRST
family = LAGRANGE
[]
[temperature]
initial_condition = ${temperature_initial}
[]
[]
[Bounds]
[deuterium_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = deuterium_concentration_W
bound_type = lower
bound_value = 0
[]
[]
[Kernels]
[time_W]
type = TimeDerivative
variable = deuterium_concentration_W
[]
[diffusion_W]
type = ADMatDiffusion
variable = deuterium_concentration_W
diffusivity = diffusivity_W
[]
[source_deuterium]
type = BodyForce
variable = deuterium_concentration_W
function = source_deuterium
[]
[]
[AuxKernels]
[temperature_aux]
type = FunctionAux
variable = temperature
function = temperature_bc_func
execute_on = 'INITIAL LINEAR'
[]
[]
[BCs]
active = 'left_recombination_flux right_recombination_flux'
# Kinetic boundary conditions
[left_recombination_flux]
type = ADMatNeumannBC
variable = deuterium_concentration_W
boundary = left
value = 1
boundary_material = flux_recombination_surface
[]
[right_recombination_flux]
type = ADMatNeumannBC
variable = deuterium_concentration_W
boundary = right
value = 1
boundary_material = flux_recombination_surface
[]
# Sieverts boundary conditions
[left_concentration_sieverts]
type = ADDirichletBC
value = '${sieverts_boundary_hat}'
boundary = left
variable = deuterium_concentration_W
[]
[right_concentration_sieverts]
type = ADDirichletBC
value = '${sieverts_boundary_hat}'
boundary = right
variable = deuterium_concentration_W
[]
[]
[Functions]
[temperature_bc_func]
type = ParsedFunction
expression = 'if(t<${charge_time_hat}, ${temperature_initial},
if(t<${fparse charge_time_hat + cooldown_duration_hat}, ${temperature_cooldown},
${temperature_desorption_min}+${desorption_heating_rate_hat}*(t-${fparse charge_time_hat + cooldown_duration_hat})))'
[]
[source_distribution]
type = ParsedFunction
expression = '1 / (${sigma_hat} * sqrt(2 * pi)) * exp(-0.5 * ((x - ${R_p_hat}) / ${sigma_hat}) ^ 2)'
[]
[surface_flux_func]
type = ParsedFunction
expression = 'if(t<${charge_time_hat}, ${surface_flux_hat}, 0)'
[]
[source_deuterium]
type = ParsedFunction
symbol_names = 'source_distribution surface_flux_func'
symbol_values = 'source_distribution surface_flux_func'
expression = 'source_distribution * surface_flux_func'
[]
[max_dt_size_function]
type = ParsedFunction
expression = 'if(t<${fparse 5 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 8 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 12 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 20 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 35 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 450 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 5000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse 11000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 13000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse (charge_time + cooldown_duration + 4500) / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 313000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 315000 / time_reference}, ${fparse 1e1 / time_reference}, ${fparse 1e3 / time_reference}))))))))))))'
[]
[max_dt_size_function_inf]
type = ParsedFunction
expression = 'if(t<${fparse 5 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 8 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 12 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 20 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 35 / time_reference}, ${fparse 1e-2 / time_reference},
if(t<${fparse 450 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 5000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse 11000 / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 13000 / time_reference}, ${fparse 1e1 / time_reference},
if(t<${fparse (charge_time + cooldown_duration + 4500) / time_reference}, ${fparse 1e2 / time_reference},
if(t<${fparse 315000 / time_reference}, ${fparse 1e1 / time_reference}, ${fparse 1e3 / time_reference})))))))))))'
[]
[max_dt_size_function_coarse]
type = ParsedFunction
expression = 'if(t<${fparse 1e-1 / time_reference}, ${fparse 1e4 / time_reference}, ${fparse 1e5 / time_reference})'
[]
[]
[Materials]
active = 'diffusivity_W_func diffusivity_nonAD recombination_rate_surface flux_recombination_surface'
[diffusivity_W_func]
type = ADDerivativeParsedMaterial
property_name = 'diffusivity_W'
functor_names = 'temperature_bc_func'
functor_symbols = 'temperature'
expression = '${diffusion_W_preexponential_hat} * exp(- ${diffusion_W_energy} / ${kb_eV} / temperature)'
[]
[diffusivity_nonAD]
type = MaterialADConverter
ad_props_in = 'diffusivity_W'
reg_props_out = 'diffusivity_W_nonAD'
[]
[recombination_rate_surface]
type = ADDerivativeParsedMaterial
property_name = 'Kr'
functor_names = 'temperature_bc_func'
functor_symbols = 'temperature'
expression = '${recombination_coefficient_hat} * exp(- ${recombination_energy} / ${kb_eV} / temperature)'
[]
[flux_recombination_surface]
type = ADDerivativeParsedMaterial
coupled_variables = 'deuterium_concentration_W'
property_name = 'flux_recombination_surface'
material_property_names = 'Kr'
expression = '- 2 * Kr * deuterium_concentration_W ^ 2'
[]
[]
[Postprocessors]
active = 'integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration
scaled_mobile_deuterium flux_surface_left scaled_flux_surface_left
flux_surface_right scaled_flux_surface_right temperature_pps diffusion_W_hat diffusion_W
max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1
integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3
integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5
integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic
spatial_max_mobile_d2 spatial_max_trapped_1 spatial_max_trapped_2 spatial_max_trapped_3 spatial_max_trapped_4 spatial_max_trapped_5 spatial_max_trapped_intrinsic
max_mobile_d2 max_trapped_1 max_trapped_2 max_trapped_3 max_trapped_4 max_trapped_5 max_trapped_intrinsic max_scaled_flux_surface_left max_scaled_flux_surface_right
max_scaled_mobile_deuterium max_scaled_trapped_deuterium_intrinsic'
[integral_source_deuterium]
type = FunctionElementIntegral
function = source_deuterium
outputs = none
[]
[scaled_implanted_deuterium]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = integral_source_deuterium
[]
[integral_deuterium_concentration]
type = ElementIntegralVariablePostprocessor
variable = deuterium_concentration_W
outputs = none
[]
[scaled_mobile_deuterium]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2}}'
value = integral_deuterium_concentration
[]
[max_scaled_mobile_deuterium]
type = TimeExtremeValue
postprocessor = scaled_mobile_deuterium
value_type = max
outputs = 'console'
[]
[flux_surface_left]
type = ADSideAverageMaterialProperty
boundary = 'left'
property = flux_recombination_surface
outputs = none
[]
[scaled_flux_surface_left]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_left
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_left]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_left
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[flux_surface_left_sieverts]
type = SideDiffusiveFluxAverage
variable = deuterium_concentration_W
boundary = 'left'
diffusivity = 'diffusivity_W_nonAD'
outputs = none
[]
[scaled_flux_surface_left_sieverts]
type = ScalePostprocessor
scaling_factor = '${fparse mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_left_sieverts
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_left_sieverts]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_left_sieverts
value_type = max
[]
[flux_surface_right]
type = ADSideAverageMaterialProperty
boundary = 'right'
property = flux_recombination_surface
outputs = none
[]
[scaled_flux_surface_right]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_right
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_right]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_right
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[flux_surface_right_sieverts]
type = SideDiffusiveFluxAverage
variable = deuterium_concentration_W
boundary = 'right'
diffusivity = 'diffusivity_W_nonAD'
outputs = none
[]
[scaled_flux_surface_right_sieverts]
type = ScalePostprocessor
scaling_factor = '${fparse -1 * mobile_concentration_reference * length_reference * ${units 1 m^2 -> mum^2} / time_reference}'
value = flux_surface_right_sieverts
execute_on = 'initial nonlinear linear timestep_end'
[]
[max_scaled_flux_surface_right_sieverts]
type = TimeExtremeValue
postprocessor = scaled_flux_surface_right_sieverts
value_type = max
[]
[temperature_pps]
type = ElementAverageValue
variable = temperature
execute_on = 'initial timestep_end'
[]
[diffusion_W_hat]
type = ElementAverageMaterialProperty
mat_prop = diffusivity_W_nonAD
outputs = none
[]
[diffusion_W]
type = ScalePostprocessor
scaling_factor = '${fparse length_reference ^ 2 / time_reference}'
value = diffusion_W_hat
outputs = none
[]
[max_time_step_size]
type = FunctionValuePostprocessor
function = max_dt_size_function
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[max_time_step_size_coarse]
type = FunctionValuePostprocessor
function = max_dt_size_function_coarse
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[spatial_max_mobile_d2]
type = NodalExtremeValue
value_type = 'max'
variable = deuterium_concentration_W
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_1]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_1
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_2]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_3]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_3
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_4]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_4
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_5]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_5
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[spatial_max_trapped_intrinsic]
type = NodalExtremeValue
value_type = 'max'
variable = trapped_intrinsic
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_mobile_d2]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_mobile_d2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_1]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_1
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_2]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_2
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_3]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_3
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_4]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_4
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_5]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_5
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[max_trapped_intrinsic]
type = TimeExtremeValue
value_type = 'max'
postprocessor = spatial_max_trapped_intrinsic
execute_on = 'initial timestep_end'
outputs = 'console'
[]
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = 'Newton'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -snes_type'
petsc_options_value = 'lu mumps vinewtonrsls'
end_time = ${endtime_hat}
line_search = 'none'
nl_rel_tol = 1e-8
nl_abs_tol = 4e-5
nl_max_its = 34
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse dt_init / time_reference}'
growth_factor = 1.1
timestep_limiting_postprocessor = max_time_step_size
[]
[Predictor]
type = SimplePredictor
scale = 1.0
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Outputs]
[csv]
type = CSV
file_base = 'val-2f_out'
[]
[csv_temperature_history]
type = CSV
show = 'temperature_pps'
file_base = 'val-2f_temperature_out'
[]
[exodus]
type = Exodus
output_material_properties = true
time_step_interval = 200
file_base = 'val-2f_out'
[]
[]
(test/tests/val-2f/tests)
[Tests]
design = 'MatNeumannBC.md CoefCoupledTimeDerivative.md TimeDerivativeNodalKernel.md TrappingNodalKernel.md ReleasingNodalKernel.md'
issues = '#12'
validation = 'val-2f.md'
[val-2f_light_csv]
type = CSVDiff
input = val-2f.i
cli_args = "sample_thickness='${units 0.9e-4 m -> mum}'
ix1=8 ix2=8 ix3=8 ix4=8 ix5=8
dt_init=2e-3
charge_time='${units 1e-1 s}' cooldown_duration='${units 1e-1 s}'
Executioner/TimeStepper/timestep_limiting_postprocessor=max_time_step_size_coarse
Outputs/csv/file_base=val-2f_light_out
Outputs/exodus/file_base=val-2f_light_out
Outputs/csv/start_time=0.2
Executioner/TimeStepper/growth_factor=1.05
Executioner/scheme=implicit-euler"
csvdiff = val-2f_light_out.csv
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport and generate CSV data output with a short runtime and coarse mesh testing.'
rel_err = 1e-3
override_columns = 'scaled_flux_surface_right scaled_flux_surface_left scaled_mobile_deuterium'
override_abs_zero = '1e4 1e6 1e7'
override_rel_err = '1e-3 1e-3 1e-3'
# Some relevant scales when considering diffs
#
# scaled_flux_surface_right O(1e12)
# scaled_flux_surface_left O(1e19)
# scaled_mobile_deuterium O(1e13)
# scaled_trapped_deuterium_intrinsic O(1e14)
[]
[val-2f_light_exodus]
type = Exodiff
input = val-2f.i
exodiff = val-2f_light_out.e
prereq = val-2f_light_csv
should_execute = false # this test relies on the output files from val-2f_light_csv, so it shouldn't be run twice
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport with a short runtime and coarse mesh testing.'
rel_err = 1e-3
abs_zero = 1e-5
[]
# Check test data consistency across operating systems and core counts for heavy tests using `Executioner/scheme=implicit-euler` due to its high stability
[val-2f_heavy_csv_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv/file_base=val-2f_heavy_desorption_out
Outputs/exodus/file_base=val-2f_heavy_desorption_out
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
csvdiff = val-2f_heavy_desorption_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output. Only desorption phase data is compared with every other timestep.'
max_time = 2000
[]
[val-2f_heavy_exodus_implicit-euler]
type = Exodiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv/file_base=val-2f_heavy_desorption_out
Outputs/exodus/file_base=val-2f_heavy_desorption_out
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
exodiff = val-2f_heavy_desorption_out.e
heavy = true
prereq = val-2f_heavy_csv_implicit-euler
should_execute = false # this test relies on the output files from val-2f_heavy_csv_implicit-euler, so it shouldn't be run twice
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability. Only desorption phase data is compared with every other timestep.'
[]
[val-2f_heavy_temperature_csv_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
prereq = val-2f_heavy_csv_implicit-euler
should_execute = false # this test relies on the output files from val-2f_heavy_csv_implicit-euler, so it shouldn't be run twice
csvdiff = val-2f_temperature_implicit_euler_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output about the temperature for documentation plots.'
max_time = 2000
[]
[val-2f_heavy_csv_inf_recombination_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "ix1='${fparse 500}' ix4='${fparse 500}'
AuxVariables/active='bounds_dummy temperature'
BCs/active='left_concentration_sieverts right_concentration_sieverts'
Physics/SpeciesTrapping/trapping_1/Ct0='trap_distribution_function_1_inf'
Physics/SpeciesTrapping/trapping_2/Ct0='trap_distribution_function_2_inf'
Physics/SpeciesTrapping/trapping_3/Ct0='trap_distribution_function_3_inf'
Physics/SpeciesTrapping/trapping_4/Ct0='trap_distribution_function_4_inf'
Physics/SpeciesTrapping/trapping_5/Ct0='trap_distribution_function_5_inf'
Materials/active='diffusivity_W_func diffusivity_nonAD'
Postprocessors/active='integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration scaled_mobile_deuterium flux_surface_left_sieverts scaled_flux_surface_left_sieverts flux_surface_right_sieverts scaled_flux_surface_right_sieverts temperature_pps max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1 integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3 integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5 integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic'
Postprocessors/max_time_step_size/function=max_dt_size_function_inf
Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/csv/file_base=val-2f_heavy_desorption_inf_recombination_out
Outputs/exodus/file_base=val-2f_heavy_desorption_inf_recombination_out
Executioner/abort_on_solve_fail=false
Executioner/nl_abs_tol=5e-4
Executioner/scheme=implicit-euler"
csvdiff = val-2f_heavy_desorption_inf_recombination_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output, for the infinite recombination case. Only desorption phase data is compared with every other timestep.'
max_time = 3600
[]
# Check the test data for the same heavy tests using 'Executioner/scheme=bdf2' for documentation. Since consistency has already been verified with `Executioner/scheme=implicit-euler` tests, the error tolerance can be relaxed.
[val-2f_heavy_csv_bdf2]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/sync_times='${fparse charge_time_hat + cooldown_duration_hat}'"
csvdiff = val-2f_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the bdf2 time integration scheme for accuracy and generate CSV data output.'
max_time = 2000
rel_err = 1e-4
override_columns = 'scaled_flux_surface_right scaled_flux_surface_left scaled_mobile_deuterium scaled_trapped_deuterium_1 scaled_trapped_deuterium_2'
override_abs_zero = '1e4 1e6 1e4 1e9 1e4'
override_rel_err = '5e-3 1e-3 9e-3 5e-3 5e-3'
[]
[val-2f_heavy_csv_inf_recombination_bdf2]
type = CSVDiff
input = val-2f.i
cli_args = "ix1='${fparse 500}' ix4='${fparse 500}'
AuxVariables/active='bounds_dummy temperature'
BCs/active='left_concentration_sieverts right_concentration_sieverts'
Physics/SpeciesTrapping/trapping_1/Ct0='trap_distribution_function_1_inf'
Physics/SpeciesTrapping/trapping_2/Ct0='trap_distribution_function_2_inf'
Physics/SpeciesTrapping/trapping_3/Ct0='trap_distribution_function_3_inf'
Physics/SpeciesTrapping/trapping_4/Ct0='trap_distribution_function_4_inf'
Physics/SpeciesTrapping/trapping_5/Ct0='trap_distribution_function_5_inf'
Materials/active='diffusivity_W_func diffusivity_nonAD'
Postprocessors/active='integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration scaled_mobile_deuterium flux_surface_left_sieverts scaled_flux_surface_left_sieverts flux_surface_right_sieverts scaled_flux_surface_right_sieverts temperature_pps diffusion_W_hat diffusion_W max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1 integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3 integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5 integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic spatial_max_mobile_d2 spatial_max_trapped_1 spatial_max_trapped_2 spatial_max_trapped_3 spatial_max_trapped_4 spatial_max_trapped_5 spatial_max_trapped_intrinsic max_mobile_d2 max_trapped_1 max_trapped_2 max_trapped_3 max_trapped_4 max_trapped_5 max_trapped_intrinsic max_scaled_flux_surface_left_sieverts max_scaled_flux_surface_right_sieverts max_scaled_mobile_deuterium max_scaled_trapped_deuterium_intrinsic'
Postprocessors/max_time_step_size/function=max_dt_size_function_inf
Outputs/csv/file_base='val-2f_out_inf_recombination'
Outputs/exodus/file_base='val-2f_out_inf_recombination'
Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/sync_times='${fparse charge_time_hat + cooldown_duration_hat}'
Executioner/abort_on_solve_fail=false
Executioner/nl_abs_tol=5e-4"
csvdiff = val-2f_out_inf_recombination.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the bdf2 time integration scheme for accuracy and generate CSV data output, for the infinite recombination case.'
max_time = 3600
rel_err = 1e-5
abs_zero = 1e4
override_columns = 'scaled_flux_surface_right_sieverts scaled_flux_surface_left_sieverts scaled_mobile_deuterium scaled_trapped_deuterium_1 scaled_trapped_deuterium_2 scaled_trapped_deuterium_3 scaled_trapped_deuterium_4 scaled_trapped_deuterium_5'
override_abs_zero = '1e6 1e9 1e6 1e9 1e9 1e9 1e9 1e10'
override_rel_err = '1e-3 3e-3 5e-3 1e-3 1e-3 3e-3 1e-3 5e-5'
[]
[val-2f_comparison]
type = RunCommand
command = 'python3 comparison_val-2f.py'
requirement = 'The system shall be able to generate comparison plots between simulated solutions and experimental data of validation case val-2f, modeling self-damaged tungsten effects on deuterium transport.'
required_python_packages = 'matplotlib numpy pandas scipy os'
[]
[]
(test/tests/val-2f/tests)
[Tests]
design = 'MatNeumannBC.md CoefCoupledTimeDerivative.md TimeDerivativeNodalKernel.md TrappingNodalKernel.md ReleasingNodalKernel.md'
issues = '#12'
validation = 'val-2f.md'
[val-2f_light_csv]
type = CSVDiff
input = val-2f.i
cli_args = "sample_thickness='${units 0.9e-4 m -> mum}'
ix1=8 ix2=8 ix3=8 ix4=8 ix5=8
dt_init=2e-3
charge_time='${units 1e-1 s}' cooldown_duration='${units 1e-1 s}'
Executioner/TimeStepper/timestep_limiting_postprocessor=max_time_step_size_coarse
Outputs/csv/file_base=val-2f_light_out
Outputs/exodus/file_base=val-2f_light_out
Outputs/csv/start_time=0.2
Executioner/TimeStepper/growth_factor=1.05
Executioner/scheme=implicit-euler"
csvdiff = val-2f_light_out.csv
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport and generate CSV data output with a short runtime and coarse mesh testing.'
rel_err = 1e-3
override_columns = 'scaled_flux_surface_right scaled_flux_surface_left scaled_mobile_deuterium'
override_abs_zero = '1e4 1e6 1e7'
override_rel_err = '1e-3 1e-3 1e-3'
# Some relevant scales when considering diffs
#
# scaled_flux_surface_right O(1e12)
# scaled_flux_surface_left O(1e19)
# scaled_mobile_deuterium O(1e13)
# scaled_trapped_deuterium_intrinsic O(1e14)
[]
[val-2f_light_exodus]
type = Exodiff
input = val-2f.i
exodiff = val-2f_light_out.e
prereq = val-2f_light_csv
should_execute = false # this test relies on the output files from val-2f_light_csv, so it shouldn't be run twice
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport with a short runtime and coarse mesh testing.'
rel_err = 1e-3
abs_zero = 1e-5
[]
# Check test data consistency across operating systems and core counts for heavy tests using `Executioner/scheme=implicit-euler` due to its high stability
[val-2f_heavy_csv_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv/file_base=val-2f_heavy_desorption_out
Outputs/exodus/file_base=val-2f_heavy_desorption_out
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
csvdiff = val-2f_heavy_desorption_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output. Only desorption phase data is compared with every other timestep.'
max_time = 2000
[]
[val-2f_heavy_exodus_implicit-euler]
type = Exodiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv/file_base=val-2f_heavy_desorption_out
Outputs/exodus/file_base=val-2f_heavy_desorption_out
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
exodiff = val-2f_heavy_desorption_out.e
heavy = true
prereq = val-2f_heavy_csv_implicit-euler
should_execute = false # this test relies on the output files from val-2f_heavy_csv_implicit-euler, so it shouldn't be run twice
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability. Only desorption phase data is compared with every other timestep.'
[]
[val-2f_heavy_temperature_csv_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/exodus/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/exodus/time_step_interval=400
Outputs/csv_temperature_history/file_base=val-2f_temperature_implicit_euler_out
Executioner/scheme=implicit-euler"
prereq = val-2f_heavy_csv_implicit-euler
should_execute = false # this test relies on the output files from val-2f_heavy_csv_implicit-euler, so it shouldn't be run twice
csvdiff = val-2f_temperature_implicit_euler_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output about the temperature for documentation plots.'
max_time = 2000
[]
[val-2f_heavy_csv_inf_recombination_implicit-euler]
type = CSVDiff
input = val-2f.i
cli_args = "ix1='${fparse 500}' ix4='${fparse 500}'
AuxVariables/active='bounds_dummy temperature'
BCs/active='left_concentration_sieverts right_concentration_sieverts'
Physics/SpeciesTrapping/trapping_1/Ct0='trap_distribution_function_1_inf'
Physics/SpeciesTrapping/trapping_2/Ct0='trap_distribution_function_2_inf'
Physics/SpeciesTrapping/trapping_3/Ct0='trap_distribution_function_3_inf'
Physics/SpeciesTrapping/trapping_4/Ct0='trap_distribution_function_4_inf'
Physics/SpeciesTrapping/trapping_5/Ct0='trap_distribution_function_5_inf'
Materials/active='diffusivity_W_func diffusivity_nonAD'
Postprocessors/active='integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration scaled_mobile_deuterium flux_surface_left_sieverts scaled_flux_surface_left_sieverts flux_surface_right_sieverts scaled_flux_surface_right_sieverts temperature_pps max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1 integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3 integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5 integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic'
Postprocessors/max_time_step_size/function=max_dt_size_function_inf
Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/time_step_interval=2
Outputs/csv/file_base=val-2f_heavy_desorption_inf_recombination_out
Outputs/exodus/file_base=val-2f_heavy_desorption_inf_recombination_out
Executioner/abort_on_solve_fail=false
Executioner/nl_abs_tol=5e-4
Executioner/scheme=implicit-euler"
csvdiff = val-2f_heavy_desorption_inf_recombination_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the implicit-euler time integration scheme for numerical stability and generate CSV data output, for the infinite recombination case. Only desorption phase data is compared with every other timestep.'
max_time = 3600
[]
# Check the test data for the same heavy tests using 'Executioner/scheme=bdf2' for documentation. Since consistency has already been verified with `Executioner/scheme=implicit-euler` tests, the error tolerance can be relaxed.
[val-2f_heavy_csv_bdf2]
type = CSVDiff
input = val-2f.i
cli_args = "Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/sync_times='${fparse charge_time_hat + cooldown_duration_hat}'"
csvdiff = val-2f_out.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the bdf2 time integration scheme for accuracy and generate CSV data output.'
max_time = 2000
rel_err = 1e-4
override_columns = 'scaled_flux_surface_right scaled_flux_surface_left scaled_mobile_deuterium scaled_trapped_deuterium_1 scaled_trapped_deuterium_2'
override_abs_zero = '1e4 1e6 1e4 1e9 1e4'
override_rel_err = '5e-3 1e-3 9e-3 5e-3 5e-3'
[]
[val-2f_heavy_csv_inf_recombination_bdf2]
type = CSVDiff
input = val-2f.i
cli_args = "ix1='${fparse 500}' ix4='${fparse 500}'
AuxVariables/active='bounds_dummy temperature'
BCs/active='left_concentration_sieverts right_concentration_sieverts'
Physics/SpeciesTrapping/trapping_1/Ct0='trap_distribution_function_1_inf'
Physics/SpeciesTrapping/trapping_2/Ct0='trap_distribution_function_2_inf'
Physics/SpeciesTrapping/trapping_3/Ct0='trap_distribution_function_3_inf'
Physics/SpeciesTrapping/trapping_4/Ct0='trap_distribution_function_4_inf'
Physics/SpeciesTrapping/trapping_5/Ct0='trap_distribution_function_5_inf'
Materials/active='diffusivity_W_func diffusivity_nonAD'
Postprocessors/active='integral_source_deuterium scaled_implanted_deuterium integral_deuterium_concentration scaled_mobile_deuterium flux_surface_left_sieverts scaled_flux_surface_left_sieverts flux_surface_right_sieverts scaled_flux_surface_right_sieverts temperature_pps diffusion_W_hat diffusion_W max_time_step_size max_time_step_size_coarse integral_trapped_concentration_1 scaled_trapped_deuterium_1 integral_trapped_concentration_2 scaled_trapped_deuterium_2 integral_trapped_concentration_3 scaled_trapped_deuterium_3 integral_trapped_concentration_4 scaled_trapped_deuterium_4 integral_trapped_concentration_5 scaled_trapped_deuterium_5 integral_trapped_concentration_intrinsic scaled_trapped_deuterium_intrinsic spatial_max_mobile_d2 spatial_max_trapped_1 spatial_max_trapped_2 spatial_max_trapped_3 spatial_max_trapped_4 spatial_max_trapped_5 spatial_max_trapped_intrinsic max_mobile_d2 max_trapped_1 max_trapped_2 max_trapped_3 max_trapped_4 max_trapped_5 max_trapped_intrinsic max_scaled_flux_surface_left_sieverts max_scaled_flux_surface_right_sieverts max_scaled_mobile_deuterium max_scaled_trapped_deuterium_intrinsic'
Postprocessors/max_time_step_size/function=max_dt_size_function_inf
Outputs/csv/file_base='val-2f_out_inf_recombination'
Outputs/exodus/file_base='val-2f_out_inf_recombination'
Outputs/csv/start_time='${fparse charge_time_hat + cooldown_duration_hat}'
Outputs/csv/sync_times='${fparse charge_time_hat + cooldown_duration_hat}'
Executioner/abort_on_solve_fail=false
Executioner/nl_abs_tol=5e-4"
csvdiff = val-2f_out_inf_recombination.csv
heavy = true
requirement = 'The system shall be able to model self-damaged tungsten effects on deuterium transport using the bdf2 time integration scheme for accuracy and generate CSV data output, for the infinite recombination case.'
max_time = 3600
rel_err = 1e-5
abs_zero = 1e4
override_columns = 'scaled_flux_surface_right_sieverts scaled_flux_surface_left_sieverts scaled_mobile_deuterium scaled_trapped_deuterium_1 scaled_trapped_deuterium_2 scaled_trapped_deuterium_3 scaled_trapped_deuterium_4 scaled_trapped_deuterium_5'
override_abs_zero = '1e6 1e9 1e6 1e9 1e9 1e9 1e9 1e10'
override_rel_err = '1e-3 3e-3 5e-3 1e-3 1e-3 3e-3 1e-3 5e-5'
[]
[val-2f_comparison]
type = RunCommand
command = 'python3 comparison_val-2f.py'
requirement = 'The system shall be able to generate comparison plots between simulated solutions and experimental data of validation case val-2f, modeling self-damaged tungsten effects on deuterium transport.'
required_python_packages = 'matplotlib numpy pandas scipy os'
[]
[]