val-2b results don't match documentation

[RemDelaporteMathurin]

Hi all!

Been having a look at the V&V recently and wanted to plot the results in the val 2b folder.

In the documentation, the following figure is displayed:

When trying to plot the data in the val 2b folder (stored in the csv files, either this one or that one), the line doesn't seem to match (see below)

To obtain this, I scaled the experimental data csv by 1e15 and the TMAP8 data by 2*1e20 as done in your plotting script.

which one is correct (ie. produced by the val-2b.i script)?

[singhgp4321]

Hi @RemDelaporteMathurin,

Could you uncomment the lines 97 and 199 in the input file (https://github.com/idaholab/TMAP8/blob/devel/test/tests/val-2b/val-2b.i) and comment out 98 and 198 and rerun? We had used different parameters for generating the figure used in documentation; those parameters allow getting a more accurate solution but are not ideal for testing purpose, so those parameters are commented out in the input file. Once you reactivate those parameters you should be able to reproduce the figure shown in documentation.

[RemDelaporteMathurin]

I haven't run TMAP8 yet I just plotted the scripts available in the repo. I understand the differences for testing purposes.

On a side note, I noticed that there was an extra 15s for the high pressure stage:

https://github.com/idaholab/TMAP8/blob/894c92aac2961569ad4b3cb6d52aaaa202d1407e/test/tests/val-2b/val-2b.i#L179

These extra 15s don't appear in the temperature expression: https://github.com/idaholab/TMAP8/blob/894c92aac2961569ad4b3cb6d52aaaa202d1407e/test/tests/val-2b/val-2b.i#L159

Running this case with FESTIM I realised that these extra 15s were needed to obtain the results shown on the Figure in the docs. This should be clearly stated in the description of the test case otherwise people may struggle to reproduce the results (unless they go and dig in the input files like I did 😄 )

[singhgp4321]

I will add that information in the documentation. Thanks for your suggestion.

[RemDelaporteMathurin]

@singhgp4321 as discussed I ran the val-2b.i script by uncommenting/commenting the lines you suggested.

This is the script I ran:

```i endtime = 197860 scale = 1e20 [Mesh] [cmg] type = CartesianMeshGenerator dim = 1 # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 dx = '0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-9 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5 0.5e-5' # 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 subdomain_id = '0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1' # 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37' [] [interface] type = SideSetsBetweenSubdomainsGenerator input = cmg primary_block = '0' #BeO paired_block = '1' # Be new_boundary = 'interface' [] [interface_other_side] type = SideSetsBetweenSubdomainsGenerator input = interface primary_block = '1' #BeO paired_block = '0' # Be new_boundary = 'interface_other' [] [] [Variables] [conc_Be] order = FIRST family = LAGRANGE initial_condition = 0 block = 1 [] [conc_BeO] order = FIRST family = LAGRANGE initial_condition = 0 block = 0 [] [] [AuxVariables] [enclosure_pressure] family = SCALAR initial_condition = 13300.0 [] [temp] initial_condition = 300 [] [flux_x] order = FIRST family = MONOMIAL [] [] [Kernels] [diff_Be] type = ADMatDiffusion variable = conc_Be diffusivity = diffusivity_Be block = 1 [] [diff_BeO] type = ADMatDiffusion variable = conc_BeO diffusivity = diffusivity_BeO block = 0 [] [time_diff_Be] type = TimeDerivative variable = conc_Be block = 1 [] [time_diff_BeO] type = TimeDerivative variable = conc_BeO block = 0 [] [] [InterfaceKernels] [tied] type = ADPenaltyInterfaceDiffusion variable = conc_BeO neighbor_var = conc_Be penalty = 0.09 # for figure generation # penalty = 0.06 # to avoid cutting time-step on civet testing jump_prop_name = solubility_ratio boundary = 'interface' [] [] [AuxScalarKernels] [enclosure_pressure_aux] type = FunctionScalarAux variable = enclosure_pressure function = enclosure_pressure_func [] [] [AuxKernels] [temp_aux] type = FunctionAux variable = temp function = temp_bc_func execute_on = 'INITIAL LINEAR' [] [flux_x_Be] type = DiffusionFluxAux diffusivity = diffusivity_Be variable = flux_x diffusion_variable = conc_Be component = x block = 1 [] [flux_x_BeO] type = DiffusionFluxAux diffusivity = diffusivity_BeO variable = flux_x diffusion_variable = conc_BeO component = x block = 0 [] [] [BCs] [left_flux] type = EquilibriumBC Ko = 5.0 activation_energy = -77966.2 boundary = left enclosure_scalar_var = enclosure_pressure temp = temp variable = conc_BeO p = 0.5 [] [right_flux] type = ADNeumannBC boundary = right variable = conc_Be value = 0 [] [] [Functions] [temp_bc_func] type = ParsedFunction value = 'if(t<180000.0, 773.0, if(t<182400.0, 773.0-((1-exp(-(t-180000)/2700))*475), 300+0.05*(t-182400)))' [] [diffusivity_BeO_func] type = ParsedFunction vars = 'T' vals = 'temp_bc_func' # value = 'if(t<182400, 1.40e-4*exp(-24408/T), 7e-5*exp(-24408/T))' value = 'if(t<182400, 1.40e-4*exp(-24408/T), 7e-5*exp(-27000/T))' # TMAP7 [] [diffusivity_Be_func] type = ParsedFunction vars = 'T' vals = 'temp_bc_func' value = '8.0e-9*exp(-4220/T)' [] [enclosure_pressure_func] type = ParsedFunction value = 'if(t<180015.0, 13300.0, if(t<182400.0, 1e-6, 0.001))' [] [solubility_BeO_func] type = ParsedFunction vars = 'T' vals = 'temp_bc_func' value = '5.00e20 * exp(9377.7/T)/${scale}' [] [solubility_Be_func] type = ParsedFunction vars = 'T' vals = 'temp_bc_func' value = '7.156e27 * exp(-11606/T)/${scale}' [] [max_time_step_size_func] type = ParsedFunction # expression = 'if(t < 170000, 10000, 100)' expression = 'if(t < 170000, 100, 10)' # for figure generation [] [] [Materials] [diff_solu] type = ADGenericFunctionMaterial prop_names = 'diffusivity_BeO diffusivity_Be solubility_Be solubility_BeO' prop_values = 'diffusivity_BeO_func diffusivity_Be_func solubility_Be_func solubility_BeO_func' outputs = all [] [converter_to_regular] type = MaterialADConverter ad_props_in = 'diffusivity_Be diffusivity_BeO' reg_props_out = 'diffusivity_Be_nonAD diffusivity_BeO_nonAD' [] [interface_jump] type = SolubilityRatioMaterial solubility_primary = solubility_BeO solubility_secondary = solubility_Be boundary = interface concentration_primary = conc_BeO concentration_secondary = conc_Be [] [] [Postprocessors] [avg_flux_left] type = SideDiffusiveFluxAverage variable = conc_BeO boundary = left diffusivity = diffusivity_BeO_nonAD [] [Temp] type = ElementAverageValue block = 1 variable = temp [] [diff_Be] type = ElementAverageValue block = 1 variable = diffusivity_Be [] [diff_BeO] type = ElementAverageValue block = 0 variable = diffusivity_BeO [] [sol_Be] type = ElementAverageValue block = 1 variable = solubility_Be [] [sol_BeO] type = ElementAverageValue block = 0 variable = solubility_BeO [] [gold_solubility_ratio] type = ParsedPostprocessor function = 'sol_BeO / sol_Be' pp_names = 'sol_BeO sol_Be' [] [BeO_interface] type = SideAverageValue boundary = interface variable = conc_BeO [] [Be_interface] type = SideAverageValue boundary = interface_other variable = conc_Be [] [variable_ratio] type = ParsedPostprocessor function = 'BeO_interface / Be_interface' pp_names = 'BeO_interface Be_interface' [] [dt] type = TimestepSize [] [h0] type = AverageElementSize block = 0 [] [h1] type = AverageElementSize block = 1 [] [Fo0] type = ParsedPostprocessor function = 'diff_BeO * dt / h0^2' pp_names = 'dt h0 diff_BeO' [] [Fo1] type = ParsedPostprocessor function = 'diff_Be * dt / h1^2' pp_names = 'dt h1 diff_Be' [] [max_time_step_size_pp] type = FunctionValuePostprocessor function = max_time_step_size_func execute_on = 'INITIAL TIMESTEP_END' outputs = none [] [] [Preconditioning] [SMP] type = SMP full = true [] [] [Executioner] type = Transient scheme = bdf2 solve_type = NEWTON petsc_options_iname = '-pc_type' petsc_options_value = 'lu' nl_rel_tol = 1e-8 nl_abs_tol = 1e-7 l_tol = 1e-4 end_time = ${endtime} automatic_scaling = true line_search = 'none' [TimeStepper] type = IterationAdaptiveDT dt = 1 optimal_iterations = 4 growth_factor = 1.1 cutback_factor = 0.5 timestep_limiting_postprocessor = max_time_step_size_pp [] [] # [Debug] # show_var_residual_norms = true # [] [Outputs] # execute_on = FINAL # exodus = true # checkpoint = true csv = true hide = 'BeO_interface Be_interface sol_Be sol_BeO diff_BeO diff_Be dt h0 h1 Fo0 Fo1 variable_ratio gold_solubility_ratio' # [exodus] # type = Exodus # output_material_properties = true # [] [] ```

Unfortunately it didn't go through as it diverged at time step 362 (time 9446.46s). It cut the stepsize several times after eventually being killed because the stepsize was too small.

This is the last nonlinear solve:

``` Solve failed, cutting timestep. Time Step 362, time = 9446.46, dt = 1e-12 0 Nonlinear |R| = 6.599805e+04 0 Linear |R| = 6.599805e+04 1 Linear |R| = 1.653900e-11 1 Nonlinear |R| = 2.977724e+03 0 Linear |R| = 2.977724e+03 1 Linear |R| = 5.318095e-13 2 Nonlinear |R| = 2.432095e+03 0 Linear |R| = 2.432095e+03 1 Linear |R| = 6.544774e-13 3 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 4 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 5 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 6 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 7 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 8 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 9 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 10 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 11 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 12 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 13 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 14 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 15 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 16 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 17 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 18 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 19 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 20 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 21 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 22 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 23 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 24 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 25 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 26 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 27 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 28 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 29 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 30 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 31 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 32 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 33 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 34 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 35 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 36 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 37 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 38 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 39 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 40 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 41 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 42 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 43 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 44 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 45 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 46 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 47 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 48 Nonlinear |R| = 1.927877e+03 0 Linear |R| = 1.927877e+03 1 Linear |R| = 4.797107e-13 49 Nonlinear |R| = 1.784836e+03 0 Linear |R| = 1.784836e+03 1 Linear |R| = 2.706537e-13 50 Nonlinear |R| = 1.927877e+03 Nonlinear solve did not converge due to DIVERGED_MAX_IT iterations 50 Solve Did NOT Converge! Aborting as solve did not converge *** ERROR *** The following error occurred in the object "IterationAdaptiveDT", of type "IterationAdaptiveDT". Solve failed and timestep already at dtmin, cannot continue! Stack frames: 13 0: libMesh::print_trace(std::ostream&) 1: moose::internal::mooseErrorRaw(std::__cxx11::basic_string, std::allocator >, std::__cxx11::basic_string, std::allocator >) 2: callMooseErrorRaw(std::__cxx11::basic_string, std::allocator >&, MooseApp*) 3: /home/remidm/festim-review/comparison/tmap8/TMAP8/moose/framework/libmoose-opt.so.0(+0x103f491) [0x7f0dad777491] 4: IterationAdaptiveDT::computeFailedDT() 5: TimeStepper::computeStep() 6: Transient::execute() 7: MooseApp::executeExecutioner() 8: MooseApp::run() 9: main 10: /lib/x86_64-linux-gnu/libc.so.6(+0x29d90) [0x7f0da541ed90] 11: __libc_start_main 12: ./tmap8-opt(+0x33ef) [0x5582f82ae3ef] application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0 [unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1 : system msg for write_line failure : Bad file descriptor ```

Any idea?

[singhgp4321]

Could you try decreasing the penalty value slightly? Instead of 0.09 may be try 0.085 or 0.08. Ideally you would want the penalty value to be high for accuracy but low enough for converging.

[RemDelaporteMathurin]

Using 0.08 seems to have done the trick

[singhgp4321]

Closing this issue as it seems you have successfully reproduced the solution described in the documentation.

[RemDelaporteMathurin]

Closing this issue as it seems you have successfully reproduced the solution described in the documentation.

I haven't checked quantitatively, however I would suggest updating the documentation to avoid any confusion as it is a bit misleading as it is.