Add support for not executing integrated BCs when the variable is no longer defined near that boundary

[GiudGiud]

Reason

See Discussion
One gets a segfault when the variable is no longer defined near the boundary set to the integratedBC
@jmeier

Design

Add a boolean to skip the contribution when executing the BC res/jacobian loops and the variable does not have local DOFs anymore.
This will have to be optional. It is just as likely that the user intended for the boundary to "follow" the variable's domain of definition, and just did not set it up properly

Impact

Simpler setups when using mesh modification. Can allow to naturally "stop" fluxes when the variable is no longer defined.

Discussed in #29131

^{Originally posted by jmeier November 25, 2024}
Dear Moose Community,

I'm trying to realize a hydro-mechanical coupled simulation using the PorousFlowSink to define boundary conditions:

• The boundary condition is defined on one of my model boundaries
• The first calculation steps all is fine and Moose works as expected.
• The moment the boundary (or parts of it) my PorousFlowSink is active on is re-assinged to an 'inactive' subdomain via a TimedSubdomainModifier, I get the following message in DBG-mode and the model will not converge:

We caught a libMesh error in ThreadedElementLoopBase:Assertion `i < _val.size()' failed.
i = 0
_val.size() = 0

[4] /home/mjg/mambaforge3/envs/moose/libmesh/include/libmesh/dense_vector.h, line 422, compiled Nov 16 2024 at 18:32:16

To recover, the solution will fail and then be re-attempted with a reduced time step.

  Nonlinear solve did not converge due to DIVERGED_FUNCTION_DOMAIN iterations 0
 Solve Did NOT Converge!
Aborting as solve did not converge

OPT-mode will work as expected and the model will converge.

I can circumvent this message (and non-convergence) in DBG-mode by deactivating the PorousFlowSink using a TimesEnableControl. But this feels a little unnecessary and adds some complexity to the model I'd like to avoid.

Furthermore, other boundary conditions (e.g. DirichletBC) will not produce any problems even if the boundary (or parts of it) they live on is re-assigned to an 'inactive' subdomain.

I'd tried to understand what exactly is the problem for Moose but did not succeed. It is even unclear to me if this is a problem in Moose or in libmesh.

Could you please comment on this behaviour:

• Is this a problem in Moose or in libmesh?
• Should Moose avoid this message/error? I am in favour of ‘yes’. But if no, then Moose should also fail to converge in OPT-mode.
• If Moose should avoid this message/error: What has to be changed to do so?

Thanks in advance!

Please find a simple test-model below (run-time should be fast)

input-file

This model will fail in DBG-mode. Uncomment the TimesEnableControl to make it converge.

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  use_displaced_mesh = false
  PorousFlowDictator = dictator
[]

Box1_inactive_name = 'Box1_inactive'
inactive_domain_block_names = '${Box1_inactive_name}'

[Problem]
  solve = true
  kernel_coverage_check = SKIP_LIST
  kernel_coverage_block_list = '${inactive_domain_block_names}'
  material_coverage_check = SKIP_LIST
  material_coverage_block_list = '${inactive_domain_block_names}'
[]

[Mesh]
  [BaseMesh]
    type = GeneratedMeshGenerator
    subdomain_name = 'BaseMesh'
    elem_type = 'HEX8'
    dim = 3
    nx = 6
    ny = 6
    nz = 2
    xmin = -3
    xmax = +3
    ymin = -3
    ymax = +3
    zmin = -2
    zmax = +2
  []

  [Box1]
    type = SubdomainBoundingBoxGenerator
    block_name = 'Box1'
    input = "BaseMesh"
    block_id = 1
    location = "INSIDE"
    bottom_left = "-1.0 -1.0 -0.0"
    top_right = "+1.0 +1.0 +2.0"
  []

  add_subdomain_names = '${inactive_domain_block_names}'
  active_block_names = 'BaseMesh Box1'
[]

[Variables]
  [disp_x]
    block = '${Mesh/active_block_names}'
  []
  [disp_y]
    block = '${Mesh/active_block_names}'
  []
  [disp_z]
    block = '${Mesh/active_block_names}'
  []
  [porepressure]
    scaling = 1e-5
    block = '${Mesh/active_block_names}'
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        incremental = true
        add_variables = false
        generate_output = 'stress_xx stress_yy stress_zz'
        block = 'BaseMesh Box1'
      []
    []
  []
[]

# ===== Kernels: PorousFlow =====
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    block = '${Mesh/active_block_names}'
    fluid_component = 0
    variable = porepressure
  []
  [flux]
    type = PorousFlowFullySaturatedDarcyFlow
    block = '${Mesh/active_block_names}'
    variable = porepressure
    gravity = '0 0 0'
    fluid_component = 0
  []
  [poro_vol_exp]
    type = PorousFlowMassVolumetricExpansion
    block = '${Mesh/active_block_names}'
    variable = porepressure
    fluid_component = 0
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_x disp_y disp_z'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2E3
    density0 = 1000
    thermal_expansion = 0
    viscosity = 9.0E-10 #MPas -2.1e-12 * exp(1808/T) --> T = 298
  []
[]

[AuxVariables]
  [dummy]
    type = MooseVariableFVReal
  []
[]

# fix the lower model boundary in y and z direction
[BCs]
  [back_fix_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'back'
    value = 0.0
  []
  [back_fix_z]
    type = DirichletBC
    variable = disp_z
    boundary = 'back'
    value = 0.0
  []
[]

# fix the left model boundary in x direction
[BCs]
  [left_fix_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
[]

# put some pressure on the right model boundary
[BCs]
  [right_Dirichlet]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 'right'
    function = right_pressure_function
  []
[]
[Functions]
  [right_pressure_function]
    type = ParsedFunction
    expression = '-0.001 * t'
  []
[]

# PorousFlowSink at 'front'
[BCs]
  [front_pfs]
    type = PorousFlowSink
    boundary = 'front'
    variable = 'porepressure'
    flux_function = 0.0
  []
[]
# [Times]
#   [times_on]
#     type = InputTimes
#     times = '0.2'
#   []
# []
# [Controls]
#   [times]
#     type = TimesEnableControl
#     times = times_on
#     disable_objects = 'BCs::front_pfs'
#     execute_on = 'INITIAL TIMESTEP_BEGIN'
#   []
# []

# Material: Volume Elements
[Materials]

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 'BaseMesh Box1'
    youngs_modulus = 1e6
    poissons_ratio = 0.25
  []

  [stress]
    type = ComputeFiniteStrainElasticStress
    block = 'BaseMesh Box1'
  []


  [temperature]
    type = PorousFlowTemperature
    block = '${Mesh/active_block_names}'
  []
  [eff_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
    block = '${Mesh/active_block_names}'
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
    block = '${Mesh/active_block_names}'
  []
  [ppss]
    type = PorousFlow1PhaseFullySaturated
    porepressure = 'porepressure'
    block = '${Mesh/active_block_names}'
  []

  [massfrac]
    type = PorousFlowMassFraction
    block = '${Mesh/active_block_names}'
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
    block = '${Mesh/active_block_names}'
  []

  [porosity_bulk]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    ensure_positive = true
    porosity_zero = 0.11
    solid_bulk = 1.3333E10
    block = '${Mesh/active_block_names}'
  []

  [undrained_density_0]
    type = GenericConstantMaterial
    block = '${Mesh/active_block_names}'
    prop_names = density
    prop_values = 2500
  []

  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.1 0 0 0 1.1 0 0 0 1.1'
  []

[]

# move elements between subdomains back and forth
[MeshModifiers]
  [GlobalSubdomainModifier]
    type = TimedSubdomainModifier
    times = '0.2'
    blocks_from = 'Box1'
    blocks_to = 'Box1_inactive'
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]

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

[Executioner]
  type = Transient

  end_time = 1.0
  dtmin = 0.001
  [TimeSteppers]
    [BlockEventTimeStepper]
      type = TimeSequenceStepper
      time_sequence = '0.05 0.1 0.2 0.4 1.0'
    []
  []

  solve_type = 'PJFNK'
  petsc_options = '-snes_converged_reason'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       mumps'

  l_tol = 1E-5
  l_max_its = 20

  nl_abs_tol = 1E-3
  nl_rel_tol = 1e-7
  nl_max_its = 20
[]

[Outputs]
  exodus = true
[]

The stack trace from bdg looks like:

Time Step 3, time = 0.2, dt = 0.1

Thread 1 "shale-dbg" hit Catchpoint 2 (exception thrown), __cxxabiv1::__cxa_throw (obj=0x555556adc420, tinfo=0x7ffff7fb73c8 <typeinfo for libMesh::LogicError>, dest=0x7ffff45470c2 <libMesh::LogicError::~LogicError()>) at ../../../../libstdc++-v3/libsupc++/eh_throw.cc:80
80      in ../../../../libstdc++-v3/libsupc++/eh_throw.cc
(gdb) bt
#0  __cxxabiv1::__cxa_throw (obj=0x555556adc420, tinfo=0x7ffff7fb73c8 <typeinfo for libMesh::LogicError>, dest=0x7ffff45470c2 <libMesh::LogicError::~LogicError()>)
    at ../../../../libstdc++-v3/libsupc++/eh_throw.cc:80
#1  0x00007ffff450116a in libMesh::DenseVector<double>::operator() (this=0x5555565b3828, i=0) at /home/mjg/mambaforge3/envs/moose/libmesh/include/libmesh/dense_vector.h:422
#2  0x00007ffff5842cd8 in IntegratedBC::computeResidual (this=0x5555565b3440) at /home/mjg/projects/moose/framework/src/bcs/IntegratedBC.C:94
#3  0x00007ffff5c45ee1 in ComputeResidualThread::compute (this=0x7fffffff87a0, ro=...) at /home/mjg/projects/moose/framework/src/loops/ComputeResidualThread.C:36
#4  0x00007ffff5c4d448 in NonlinearThread::compute (this=0x7fffffff87a0, bc=...) at /home/mjg/projects/moose/framework/src/loops/NonlinearThread.C:281
#5  0x00007ffff5c4cb46 in NonlinearThread::computeOnBoundary (this=0x7fffffff87a0, bnd_id=5) at /home/mjg/projects/moose/framework/src/loops/NonlinearThread.C:180
#6  0x00007ffff5c4c9d6 in NonlinearThread::onBoundary (this=0x7fffffff87a0, elem=0x555555f98990, side=5, bnd_id=5, lower_d_elem=0x0)
    at /home/mjg/projects/moose/framework/src/loops/NonlinearThread.C:164
#7  0x00007ffff5c68116 in ThreadedElementLoopBase<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> >::operator() (this=0x7fffffff87a0, range=..., 
    bypass_threading=false) at /home/mjg/projects/moose/framework/build/header_symlinks/ThreadedElementLoopBase.h:277
#8  0x00007ffff5c4bf14 in NonlinearThread::operator() (this=0x7fffffff87a0, range=..., bypass_threading=false) at /home/mjg/projects/moose/framework/src/loops/NonlinearThread.C:60
#9  0x00007ffff596e80f in libMesh::Threads::parallel_reduce<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>, ComputeResidualThread> (range=..., body=...)
    at /home/mjg/mambaforge3/envs/moose/libmesh/include/libmesh/threads_pthread.h:367
#10 0x00007ffff5931661 in NonlinearSystemBase::computeResidualInternal (this=0x5555560053e0, tags=...) at /home/mjg/projects/moose/framework/src/systems/NonlinearSystemBase.C:1735
#11 0x00007ffff5928b22 in NonlinearSystemBase::computeResidualTags (this=0x5555560053e0, tags=...) at /home/mjg/projects/moose/framework/src/systems/NonlinearSystemBase.C:826
#12 0x00007ffff6392c82 in FEProblemBase::computeResidualTags (this=0x555555fd79e0, tags=...) at /home/mjg/projects/moose/framework/src/problems/FEProblemBase.C:7002
#13 0x00007ffff638fef7 in FEProblemBase::computeResidualInternal (this=0x555555fd79e0, soln=..., residual=..., tags=...) at /home/mjg/projects/moose/framework/src/problems/FEProblemBase.C:6842
#14 0x00007ffff638e27a in FEProblemBase::computeResidual (this=0x555555fd79e0, soln=..., residual=..., nl_sys_num=0) at /home/mjg/projects/moose/framework/src/problems/FEProblemBase.C:6646
#15 0x00007ffff638ddad in FEProblemBase::computeResidualSys (this=0x555555fd79e0, sys=..., soln=..., residual=...) at /home/mjg/projects/moose/framework/src/problems/FEProblemBase.C:6617
#16 0x00007ffff59121bb in ComputeResidualFunctor::residual (this=0x555556006cd8, soln=..., residual=..., sys=...) at /home/mjg/projects/moose/framework/src/systems/ComputeResidualFunctor.C:27
#17 0x00007fffe9339337 in libmesh_petsc_snes_residual () from /home/mjg/mambaforge3/envs/moose/libmesh/lib/libmesh_dbg.so.0
#18 0x00007fffe1f3e3a1 in SNESComputeFunction () from /home/mjg/mambaforge3/envs/moose/petsc/lib/libpetsc.so.3.21
#19 0x00007fffe1f18481 in SNESSolve_NEWTONLS () from /home/mjg/mambaforge3/envs/moose/petsc/lib/libpetsc.so.3.21
#20 0x00007fffe1f46065 in SNESSolve () from /home/mjg/mambaforge3/envs/moose/petsc/lib/libpetsc.so.3.21
#21 0x00007fffe9340d9e in libMesh::PetscNonlinearSolver<double>::solve(libMesh::SparseMatrix<double>&, libMesh::NumericVector<double>&, libMesh::NumericVector<double>&, double, unsigned int) ()
   from /home/mjg/mambaforge3/envs/moose/libmesh/lib/libmesh_dbg.so.0
#22 0x00007fffe93de6a8 in libMesh::NonlinearImplicitSystem::solve() () from /home/mjg/mambaforge3/envs/moose/libmesh/lib/libmesh_dbg.so.0
#23 0x00007ffff591e7d0 in NonlinearSystem::solve (this=0x5555560053e0) at /home/mjg/projects/moose/framework/src/systems/NonlinearSystem.C:197
#24 0x00007ffff63840a5 in FEProblemBase::solve (this=0x555555fd79e0, nl_sys_num=0) at /home/mjg/projects/moose/framework/src/problems/FEProblemBase.C:6208
#25 0x00007ffff5b3bbfa in FEProblemSolve::solve (this=0x5555560606c8) at /home/mjg/projects/moose/framework/src/executioners/FEProblemSolve.C:294
#26 0x00007ffff5b41126 in FixedPointSolve::solveStep (this=0x555556060d10, begin_norm=@0x5555567da990: 0, end_norm=@0x5555567e8a60: 0, transformed_dofs=...)
    at /home/mjg/projects/moose/framework/src/executioners/FixedPointSolve.C:439
#27 0x00007ffff5b3ffb2 in FixedPointSolve::solve (this=0x555556060d10) at /home/mjg/projects/moose/framework/src/executioners/FixedPointSolve.C:291
#28 0x00007ffff4edd150 in TimeStepper::step (this=0x555556122b40) at /home/mjg/projects/moose/framework/src/timesteppers/TimeStepper.C:168
#29 0x00007ffff4edb22b in TimeSequenceStepperBase::step (this=0x555556122b40) at /home/mjg/projects/moose/framework/src/timesteppers/TimeSequenceStepperBase.C:124
#30 0x00007ffff5b503f6 in Transient::takeStep (this=0x555556060310, input_dt=-1) at /home/mjg/projects/moose/framework/src/executioners/Transient.C:418
#31 0x00007ffff5b4ee76 in Transient::execute (this=0x555556060310) at /home/mjg/projects/moose/framework/src/executioners/Transient.C:293
#32 0x00007ffff65a9794 in MooseApp::executeExecutioner (this=0x555555946ae0) at /home/mjg/projects/moose/framework/src/base/MooseApp.C:1180
#33 0x00007ffff65b0b64 in MooseApp::run (this=0x555555946ae0) at /home/mjg/projects/moose/framework/src/base/MooseApp.C:1562
#34 0x0000555555558e94 in Moose::main<shaleTestApp> (argc=3, argv=0x7fffffffbf28) at /home/mjg/projects/moose/framework/build/header_symlinks/MooseMain.h:47
#35 0x00005555555586fd in main (argc=3, argv=0x7fffffffbf28) at /home/mjg/projects/shale/src/main.C:17

PS: here a similar error is described #9659