PorousFlowSink: libMesh error in ThreadedElementLoopBase:Assertion `i < _val.size()' failed. #29131

jmeier · 2024-11-25T11:12:52Z

jmeier
Nov 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

jmeier · 2024-11-25T13:00:06Z

jmeier
Nov 25, 2024
Author

I've seen that in my test-model above the error seems to be triggered in line 94 of IntegratedBC.C:

moose/framework/src/bcs/IntegratedBC.C

Lines 90 to 95 in 37b999b

    
           for (_qp = 0; _qp < _qrule->n_points(); _qp++) 
        
           { 
        
             precalculateQpResidual(); 
        
             for (_i = 0; _i < _test.size(); _i++) 
        
               _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual(); 
        
           }

Here, _test.size() returns 8 while _local_re.size() is 0. Would it be an option to change line 93 from
for (_i = 0; _i < _test.size(); _i++)
to
for (_i = 0; _i < _local_re.size(); _i++)
?

Similar are ...

... line 117 + 118

moose/framework/src/bcs/IntegratedBC.C

Lines 116 to 120 in 37b999b

    
             precalculateQpJacobian(); 
        
             for (_i = 0; _i < _test.size(); _i++) 
        
               for (_j = 0; _j < _phi.size(); _j++) 
        
                 _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian(); 
        
           }

... line 161 + 162

moose/framework/src/bcs/IntegratedBC.C

Lines 158 to 164 in 37b999b

    
           for (_qp = 0; _qp < _qrule->n_points(); _qp++) 
        
           { 
        
             precalculateQpOffDiagJacobian(jvar); 
        
             for (_i = 0; _i < _test.size(); _i++) 
        
               for (_j = 0; _j < phi_size; _j++) 
        
                 _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar_num); 
        
           }

... line 176

moose/framework/src/bcs/IntegratedBC.C

Lines 175 to 178 in 37b999b

    
           for (_qp = 0; _qp < _qrule->n_points(); _qp++) 
        
             for (_i = 0; _i < _test.size(); _i++) 
        
               for (_j = 0; _j < jv.order(); _j++) 
        
                 _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobianScalar(jvar);

12 replies

jmeier Dec 4, 2024
Author

I'm not sure I understand the last question. I'll try to answer:

In the example attached in my initial post above, the BC is active from the beginning (t=0). The BC covers all of the model boundary at z=+2 (this boundary is named "front") as shown with the magenta dots in the image below.

At t=0.2 the GlobalSubdomainModifier re-assignes all elements of the subdomain "Box1" shown in red above to an 'inactive' subdomain "Box1_inactive". This is the moment, Moose will show the error reported.

GiudGiud Dec 4, 2024
Collaborator

ok this makes sense now!
We need the boundary condition to not be defined on any boundary next to this inactive block anymore.
The easiest way is to modify the boundary at the same time as the block.

jmeier Dec 5, 2024
Author

We need the boundary condition to not be defined
on any boundary next to this inactive block anymore.
The easiest way is to modify the boundary at the same time as the block.

The input file in the initial post uses a controller on the BC for this purpose (commented out to show the error).

The motivation of this discussion is to avoid those additional controller-blocks in the input file. 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. It would be great it we can change the behaviour of IntegratedBC etc. to allow this.

GiudGiud Dec 5, 2024
Collaborator

I dont think the fact that it works for the DirichletBC is intentional either.
I'll create an issue. Considering you already have a solution, I ll take a look end of December

lindsayad Dec 6, 2024
Maintainer

A failed assertion is definitely not the outcome we want. We would like a uniform behavior between the two BC types. An up-front error is probably best, or in this case I suppose an error at the time the mesh is changed. Probably "silently not executing a residual function" is a distant second place

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

PorousFlowSink: libMesh error in ThreadedElementLoopBase:Assertion `i < _val.size()' failed. #29131

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 12 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

PorousFlowSink: libMesh error in ThreadedElementLoopBase:Assertion `i < _val.size()' failed. #29131

jmeier Nov 25, 2024

Replies: 1 comment · 12 replies

jmeier Nov 25, 2024 Author

jmeier Dec 4, 2024 Author

GiudGiud Dec 4, 2024 Collaborator

jmeier Dec 5, 2024 Author

GiudGiud Dec 5, 2024 Collaborator

lindsayad Dec 6, 2024 Maintainer

jmeier
Nov 25, 2024

Replies: 1 comment 12 replies

jmeier
Nov 25, 2024
Author

jmeier Dec 4, 2024
Author

GiudGiud Dec 4, 2024
Collaborator

jmeier Dec 5, 2024
Author

GiudGiud Dec 5, 2024
Collaborator

lindsayad Dec 6, 2024
Maintainer