Bug Description

Hello!

I am trying to reproduce permeation data from an experiment (reference: https://doi.org/10.1016/j.fusengdes.2008.05.016) for validation purposes. The problem can be reconducted to a two-layers composite solid structure, where on one side a tritium partial pressure is imposed and on the other side the concentration is null. The two layers one models a nickel plate, where Sievert's law applied, and the other a layer of stationary molten FLiBe, where Henry's law applies.

My issue is with the interface BC between the two regions. I want to use the recently updated InterfaceSoprtion BC but I get a segmentation fault on the first timestep of the simulation. It probably is a machinie-specific bug since the same input runs fine on another pc. It shouldn't be an environment problem since the problem persists even after conda and TMAP( have been reinstalled.

This is the trace I get on MPI_Abort, running in debug mode:

 Thread 1 "tmap8-dbg" hit Breakpoint 1, 0x00007fffdf53e4f0 in PMPI_Abort () from /home/federico/miniforge/envs/moose/lib/libmpi.so.12
(gdb) bt
#0  0x00007fffdf53e4f0 in PMPI_Abort () from /home/federico/miniforge/envs/moose/lib/libmpi.so.12
#1  0x00007ffff5fec3bc in moose::internal::mooseErrorRaw (msg=..., prefix=...)
    at /home/federico/projects/TMAP8/moose/framework/src/base/MooseError.C:88
#2  0x00007ffff5fead92 in MooseBase::callMooseError (this=0x555555f79e90, msg=..., with_prefix=true)
    at /home/federico/projects/TMAP8/moose/framework/src/base/MooseBase.C:39
#3  0x00007ffff4a3f5a7 in MooseBaseErrorInterface::mooseError<char const (&) [61]> (this=0x555555f79f10)
    at /home/federico/projects/TMAP8/moose/framework/build/header_symlinks/MooseBaseErrorInterface.h:33
#4  0x00007ffff5f3edb0 in IterationAdaptiveDT::computeFailedDT (this=0x555555f79e90)
    at /home/federico/projects/TMAP8/moose/framework/src/timesteppers/IterationAdaptiveDT.C:328
#5  0x00007ffff5f45680 in TimeStepper::computeStep (this=0x555555f79e90)
    at /home/federico/projects/TMAP8/moose/framework/src/timesteppers/TimeStepper.C:87
#6  0x00007ffff5be9040 in Transient::computeDT (this=0x555555ed27a0)
    at /home/federico/projects/TMAP8/moose/framework/src/executioners/Transient.C:337
#7  0x00007ffff5be81dd in Transient::execute (this=0x555555ed27a0)
    at /home/federico/projects/TMAP8/moose/framework/src/executioners/Transient.C:290
#8  0x00007ffff5fd712e in MooseApp::executeExecutioner (this=0x555555907590)
    at /home/federico/projects/TMAP8/moose/framework/src/base/MooseApp.C:1172
#9  0x00007ffff5fde4f7 in MooseApp::run (this=0x555555907590)
    at /home/federico/projects/TMAP8/moose/framework/src/base/MooseApp.C:1554
#10 0x0000555555558e7e in Moose::main<TMAP8TestApp> (argc=3, argv=0x7fffffffcc88)
    at /home/federico/projects/TMAP8/moose/framework/build/header_symlinks/MooseMain.h:47
#11 0x00005555555586e6 in main (argc=3, argv=0x7fffffffcc88) at /home/federico/projects/TMAP8/src/main.C:15
(gdb)

I tested the same input swapping the BC with the ADPenaltyInterfaceDiffusion BC and the case runs. There might be a better way to model the BC, any idea is appreciated

Thanks!

Steps to Reproduce

What I did was simply a fresh TMAP8 installation and I run the case n°1 below. It doesn't run on my machine. Then I swapped the BC in the case n°2 and it runs.

Impact

This issue is preventing me from proposing two validation cases for TMAP8 on the V&V list.

Case n°1 - InterfaceSorption BC

endtime = 140000

R = 8.31446261815324 # Gas constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)

T = '${units 550 degC -> K}' # temperature

p_bnd = 1210 # pressure

L_Ni = ${units 2 mm -> m} # nickel thickness
L_salt = ${units 8.1 mm -> m} # salt thickness

num_nodes = 200 # (-)

[Mesh]
  [whole_domain]
    type = GeneratedMeshGenerator
    xmin = 0
    xmax = ${fparse L_Ni + L_salt}
    dim = 1
    nx = ${num_nodes}
  []
  [block_1]
    type = ParsedSubdomainMeshGenerator
    input = whole_domain
    combinatorial_geometry = 'x <= ${L_Ni}'
    block_id = 0
  []
  [block_2]
    type = ParsedSubdomainMeshGenerator
    input = block_1
    combinatorial_geometry = 'x > ${L_Ni}'
    block_id = 1
  [] 
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = block_2
    primary_block = '0' # Ni
    paired_block = '1' # salt
    new_boundary = 'interface'
  []
  [interface_other_side]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = '1' # salt
    paired_block = '0' # Ni
    new_boundary = 'interface_other'
  []
[]

[Variables]
  [conc_Ni]
    initial_condition = 0.9e-0
    block = 0
  []
  [conc_salt]
    initial_condition = 0.3e-0
    block = 1
  [] 
[]

[AuxVariables]
  [enclosure_pressure]
    family = SCALAR
    initial_condition = ${p_bnd}
  []
  [flux_x]
    order = FIRST
    family = MONOMIAL
  []
  [p_div_RT_salt]
    block = 1
    initial_condition = 0.001
  []
  [p_div_RT_Ni]
    block = 0
    initial_condition = 0.001
  []
[]

[Kernels]
  [diff_Ni]
    type = ADMatDiffusion
    variable = conc_Ni
    diffusivity = diffusivity_Ni
    block = 0
  []
  [diff_salt]
    type = ADMatDiffusion
    variable = conc_salt
    diffusivity = diffusivity_salt
    block = 1
  []
  [time_diff_Ni]
    type = TimeDerivative
    variable = conc_Ni
    block = 0
  []
  [time_diff_salt]
    type = TimeDerivative
    variable = conc_salt
    block = 1
  []
[]

[InterfaceKernels]
  [tied]
    type = InterfaceSorption
    K0 = 0.564
    Ea = 15800.0
    n_sorption = 0.5
    diffusivity = diffusivity_Ni_nonAD
    unit_scale = 1
    unit_scale_neighbor = 1
    temperature = ${T}
    variable = conc_Ni
    neighbor_var = p_div_RT_salt
    sorption_penalty = 0.1
    boundary = 'interface'
  []
[]

[AuxKernels]
  [flux_x_Ni]
    type = DiffusionFluxAux
    diffusivity = diffusivity_Ni
    variable = flux_x
    diffusion_variable = conc_Ni
    component = x
    block = 0
  []
  [flux_x_salt]
    type = DiffusionFluxAux
    diffusivity = diffusivity_salt
    variable = flux_x
    diffusion_variable = conc_salt
    component = x
    block = 1
  []
  [p_div_RT_Ni_kernel]
    variable = p_div_RT_Ni
    type = ParsedAux
    expression = '(conc_Ni / (${R} * ${T} * 0.564 * exp(-15800/(${R}*${T}))))^2'
    coupled_variables = 'conc_Ni'
  []
  [p_div_RT_salt_kernel]
    variable = p_div_RT_salt
    type = ParsedAux
    expression = '(conc_salt / (${R} * ${T} * 0.079 * exp(-35000/(${R}*${T}))))'
    coupled_variables = 'conc_salt'
  []
[]

[BCs]
  [left_flux]
    type = EquilibriumBC
    Ko = 0.564
    activation_energy = 15800.0
    boundary = left
    enclosure_var = enclosure_pressure
    temperature = ${T}
    variable = conc_Ni
    p = 0.5
  []
  [right_flux]
    type = ADDirichletBC
    boundary = right
    variable = conc_salt
    value = 0.0
  []
[]

[Functions]

  [diffusivity_Ni_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = '${T}'
    expression = '0.0000007*exp(-39500/(${R}*T))'
  []

  [diffusivity_salt_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = '${T}'
    expression = '0.00000093*exp(-42000/(${R}*T))'
  []

  [solubility_Ni_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = '${T}'
    expression = '0.564 * exp(-15800/(${R}*T))'
  []

  [solubility_salt_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = '${T}'
    expression = '0.079 * exp(-35000/(${R}*T))'
  []
[]

[Materials]
  [diff_solu]
    type = ADGenericFunctionMaterial
    prop_names = 'diffusivity_Ni diffusivity_salt solubility_Ni solubility_salt'
    prop_values = 'diffusivity_Ni_func diffusivity_salt_func solubility_Ni_func solubility_salt_func'
    outputs = all
  []
  [converter_to_regular]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_Ni diffusivity_salt'
    reg_props_out = 'diffusivity_Ni_nonAD diffusivity_salt_nonAD'
  []
[]

[Postprocessors]
  [avg_flux_right]
    type = SideDiffusiveFluxAverage
    variable = conc_salt
    boundary = right
    diffusivity = diffusivity_salt_nonAD
  []
[]

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

[Executioner]
  type = Transient
  steady_state_detection = true
  steady_state_start_time = 40000
  steady_state_tolerance = 1e-9

  scheme = bdf2 # bdf2 # crank-nicolson # explicit-euler
  solve_type = NEWTON # LINEAR # JFNK # NEWTON
  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  l_max_its = 10
  nl_max_its = 13
  nl_rel_tol = 1e-8 # nonlinear relative tolerance
  nl_abs_tol = 1e-20 #1e-30 # nonlinear absolute tolerance
  l_tol = 1e-6 # 1e-3 - 1e-5 # linear tolerance
  end_time = ${endtime}
  automatic_scaling = true
  line_search = none
  dtmax = 10.0
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-10
    optimal_iterations = 18 # 6-10 or 18
    growth_factor = 1.1
    cutback_factor = 0.5
  []
[]

[Outputs]
  execute_on = timestep_end
  exodus = true
  # checkpoint = true
  # 
  # hide = 'dt h0 h1 Fo0 Fo1 variable_ratio'
  [csv]
    type = CSV
    file_base = 'val-flibe_1210_550'
    time_step_interval = 500
  []
[]

[Dampers]
  [limit_salt]
    type = BoundingValueElementDamper
    variable = conc_salt
    max_value = 1e42
    min_value = -0.01
    min_damping = 0.001
  []
  [limit_Ni]
    type = BoundingValueElementDamper
    variable = conc_Ni
    max_value = 1e42
    min_value = -0.01
    min_damping = 0.001
  []
[]

Case n°2 - ADPenaltyInterfaceDiffusion BC

endtime = 140000

R = 8.31446261815324 # Gas constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)

T = '${units 550 degC -> K}' # temperature

p_bnd = 1210 # pressure

L_Ni = ${units 2 mm -> m} # nickel thickness
L_salt = ${units 8.1 mm -> m} # salt thickness

num_nodes = 200 # (-)

[Mesh]
  [whole_domain]
    type = GeneratedMeshGenerator
    xmin = 0
    xmax = ${fparse L_Ni + L_salt}
    dim = 1
    nx = ${num_nodes}
  []
  [block_1]
    type = ParsedSubdomainMeshGenerator
    input = whole_domain
    combinatorial_geometry = 'x <= ${L_Ni}'
    block_id = 0
  []
  [block_2]
    type = ParsedSubdomainMeshGenerator
    input = block_1
    combinatorial_geometry = 'x > ${L_Ni}'
    block_id = 1
  [] 
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = block_2
    primary_block = '0' # Ni
    paired_block = '1' # salt
    new_boundary = 'interface'
  []
  [interface_other_side]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = '1' # salt
    paired_block = '0' # Ni
    new_boundary = 'interface_other'
  []
[]

[Variables]
  [conc_Ni]
    initial_condition = 0.9e-0
    block = 0
  []
  [conc_salt]
    initial_condition = 0.3e-0
    block = 1
  [] 
[]

[AuxVariables]
  [enclosure_pressure]
    family = SCALAR
    initial_condition = ${p_bnd}
  []
  [flux_x]
    order = FIRST
    family = MONOMIAL
  []
  [p_div_RT_salt]
    block = 1
    initial_condition = 0.001
  []
  [p_div_RT_Ni]
    block = 0
    initial_condition = 0.001
  []
[]

[Kernels]
  [diff_Ni]
    type = ADMatDiffusion
    variable = conc_Ni
    diffusivity = diffusivity_Ni
    block = 0
  []
  [diff_salt]
    type = ADMatDiffusion
    variable = conc_salt
    diffusivity = diffusivity_salt
    block = 1
  []
  [time_diff_Ni]
    type = TimeDerivative
    variable = conc_Ni
    block = 0
  []
  [time_diff_salt]
    type = TimeDerivative
    variable = conc_salt
    block = 1
  []
[]

[InterfaceKernels]
  [tied]
    type = ADPenaltyInterfaceDiffusion
    variable = conc_Ni
    neighbor_var = conc_salt
    penalty = 0.06 # to avoid cutting time-step on civet testing
    jump_prop_name = solubility_ratio
    boundary = 'interface'
  []
[]

[AuxKernels]
  [flux_x_Ni]
    type = DiffusionFluxAux
    diffusivity = diffusivity_Ni
    variable = flux_x
    diffusion_variable = conc_Ni
    component = x
    block = 0
  []
  [flux_x_salt]
    type = DiffusionFluxAux
    diffusivity = diffusivity_salt
    variable = flux_x
    diffusion_variable = conc_salt
    component = x
    block = 1
  []
  [p_div_RT_Ni_kernel]
    variable = p_div_RT_Ni
    type = ParsedAux
    expression = '(conc_Ni / (${R} * ${T} * 0.564 * exp(-15800/(${R}*${T}))))^2'
    coupled_variables = 'conc_Ni'
  []
  [p_div_RT_salt_kernel]
    variable = p_div_RT_salt
    type = ParsedAux
    expression = '(conc_salt / (${R} * ${T} * 0.079 * exp(-35000/(${R}*${T}))))'
    coupled_variables = 'conc_salt'
  []
[]

[BCs]
  [left_flux]
    type = EquilibriumBC
    Ko = 0.564
    activation_energy = 15800.0
    boundary = left
    enclosure_var = enclosure_pressure
    temperature = ${T}
    variable = conc_Ni
    p = 0.5
  []
  [right_flux]
    type = ADDirichletBC
    boundary = right
    variable = conc_salt
    value = 0.0
  []
[]

[Functions]

  [diffusivity_Ni_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = '${T}'
    expression = '0.0000007*exp(-39500/(${R}*T))'
  []

  [diffusivity_salt_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = '${T}'
    expression = '0.00000093*exp(-42000/(${R}*T))'
  []

  [solubility_Ni_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = '${T}'
    expression = '0.564 * exp(-15800/(${R}*T))'
  []

  [solubility_salt_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = '${T}'
    expression = '0.079 * exp(-35000/(${R}*T))'
  []
[]

[Materials]
  [diff_solu]
    type = ADGenericFunctionMaterial
    prop_names = 'diffusivity_Ni diffusivity_salt solubility_Ni solubility_salt'
    prop_values = 'diffusivity_Ni_func diffusivity_salt_func solubility_Ni_func solubility_salt_func'
    outputs = all
  []
  [converter_to_regular]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_Ni diffusivity_salt'
    reg_props_out = 'diffusivity_Ni_nonAD diffusivity_salt_nonAD'
  []

  [interface_jump]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_Ni
    solubility_secondary = solubility_salt
    boundary = interface
    concentration_primary = conc_Ni
    concentration_secondary = conc_salt
  []
[]

[Postprocessors]
  [avg_flux_right]
    type = SideDiffusiveFluxAverage
    variable = conc_salt
    boundary = right
    diffusivity = diffusivity_salt_nonAD
  []
[]

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

[Executioner]
  type = Transient
  steady_state_detection = true
  steady_state_start_time = 40000
  steady_state_tolerance = 1e-9

  scheme = bdf2 # bdf2 # crank-nicolson # explicit-euler
  solve_type = NEWTON # LINEAR # JFNK # NEWTON
  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  l_max_its = 10
  nl_max_its = 13
  nl_rel_tol = 1e-8 # nonlinear relative tolerance
  nl_abs_tol = 1e-20 #1e-30 # nonlinear absolute tolerance
  l_tol = 1e-6 # 1e-3 - 1e-5 # linear tolerance
  end_time = ${endtime}
  automatic_scaling = true
  line_search = none
  dtmax = 10.0
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-10
    optimal_iterations = 18 # 6-10 or 18
    growth_factor = 1.1
    cutback_factor = 0.5
  []
[]

[Outputs]
  execute_on = timestep_end
  exodus = true
  # checkpoint = true
  # 
  # hide = 'dt h0 h1 Fo0 Fo1 variable_ratio'
  [csv]
    type = CSV
    file_base = 'val-flibe_1210_550'
    time_step_interval = 500
  []
[]

[Dampers]
  [limit_salt]
    type = BoundingValueElementDamper
    variable = conc_salt
    max_value = 1e42
    min_value = -0.01
    min_damping = 0.001
  []
  [limit_Ni]
    type = BoundingValueElementDamper
    variable = conc_Ni
    max_value = 1e42
    min_value = -0.01
    min_damping = 0.001
  []
[]

idaholab / TMAP8

Segmentation fault on first timestep when using InterfaceSorption InterfaceKernel #186

Bug Description

Steps to Reproduce

Impact

Case n°1 - InterfaceSorption BC

Case n°2 - ADPenaltyInterfaceDiffusion BC