Material models do not check whether `Enable elasticity` is true even if they don't support it

[anne-glerum]

At the moment, when one sets Enable elasticity to true in Formulation, while using a Material model that does not support elasticity, ASPECT will crash with an unhelpful error while solving the Stokes solver (see below), instead of throwing a clear error before even starting the timestep. Would it be better to have each Material model check whether Enable elasticity is true and then throw, or to check whether a Material model that supports elasticity is used when Enable elasticity is read in as true? The former avoids having to add new Material model names to the Assert when they support elasticity in the future.

Happy to make a PR with the preferred fix.

This came up in #4370.

*** Timestep 0:  t=0 years, dt=0 years
   Solving temperature system... 0 iterations.
   Solving Stokes system... 
--------------------------------------------------------
An error occurred in line <1841> of file </var/folders/8z/hlb6vc015qjggytkxn84m6_c0000gn/T/heltai/spack-stage/spack-stage-dealii-9.3.0-zy7k3uwnakcqjvrajvacy5l4jrl7eaex/spack-src/include/deal.II/lac/la_parallel_vector.templates.h> in function
    typename Vector<Number, MemorySpaceType>::real_type dealii::LinearAlgebra::distributed::Vector<double, dealii::MemorySpace::Host>::norm_sqr_local() const
The violated condition was: 
    dealii::numbers::is_finite(sum)
Additional information: 
    In a significant number of places, deal.II checks that some
    intermediate value is a finite number (as opposed to plus or minus
    infinity, or NaN/Not a Number). In the current function, we
    encountered a number that is not finite (its value is (nan,0) and 
    therefore violates the current assertion).
    This may be due to the fact that some operation in this function
    created such a value, or because one of the arguments you passed to
    the function already had this value from some previous operation. In
    the latter case, this function only triggered the error but may not
    actually be responsible for the computation of the number that is not
    finite.

    There are two common cases where this situation happens. First, your
    code (or something in deal.II) divides by zero in a place where this
    should not happen. Or, you are trying to solve a linear system with an
    unsuitable solver (such as an indefinite or non-symmetric linear
    system using a Conjugate Gradient solver); such attempts oftentimes
    yield an operation somewhere that tries to divide by zero or take the
    square root of a negative value.

    In any case, when trying to find the source of the error, recall that
    the location where you are getting this error is simply the first
    place in the program where there is a check that a number (e.g., an
    element of a solution vector) is in fact finite, but that the actual
    error that computed the number may have happened far earlier. To find
    this location, you may want to add checks for finiteness in places of
    your program visited before the place where this error is produced.
    One way to check for finiteness is to use the 'AssertIsFinite' macro.

Stacktrace:
-----------
#0  2   libdeal_II.g.9.3.0.dylib            0x0000000127e8939c _ZNK6dealii13LinearAlgebra11distributed6VectorIdNS_11MemorySpace4HostEE14norm_sqr_localEv + 172: 2   libdeal_II.g.9.3.0.dylib            0x0000000127e8939c _ZNK6dealii13LinearAlgebra11distributed6VectorIdNS_11MemorySpace4HostEE14norm_sqr_localEv
#1  3   libdeal_II.g.9.3.0.dylib            0x0000000127ecd962 _ZNK6dealii13LinearAlgebra11distributed11BlockVectorIdE8norm_sqrEv + 82: 3   libdeal_II.g.9.3.0.dylib            0x0000000127ecd962 _ZNK6dealii13LinearAlgebra11distributed11BlockVectorIdE8norm_sqrEv
#2  4   libdeal_II.g.9.3.0.dylib            0x0000000127ecd909 _ZNK6dealii13LinearAlgebra11distributed11BlockVectorIdE7l2_normEv + 9: 4   libdeal_II.g.9.3.0.dylib            0x0000000127ecd909 _ZNK6dealii13LinearAlgebra11distributed11BlockVectorIdE7l2_normEv
#3  5   aspect                              0x000000010e953488 _ZN6aspect37StokesMatrixFreeHandlerImplementationILi2ELi2EE5solveEv + 3224: 5   aspect                              0x000000010e953488 _ZN6aspect37StokesMatrixFreeHandlerImplementationILi2ELi2EE5solveEv
#4  6   aspect                              0x000000010ebe35bd _ZN6aspect9SimulatorILi2EE12solve_stokesEv + 189: 6   aspect                              0x000000010ebe35bd _ZN6aspect9SimulatorILi2EE12solve_stokesEv
#5  7   aspect                              0x000000010ebe87f8 _ZN6aspect9SimulatorILi2EE25assemble_and_solve_stokesEbPd + 184: 7   aspect                              0x000000010ebe87f8 _ZN6aspect9SimulatorILi2EE25assemble_and_solve_stokesEbPd
#6  8   aspect                              0x000000010ebe8a81 _ZN6aspect9SimulatorILi2EE36solve_single_advection_single_stokesEv + 65: 8   aspect                              0x000000010ebe8a81 _ZN6aspect9SimulatorILi2EE36solve_single_advection_single_stokesEv
#7  9   aspect                              0x000000010e8286ab _ZN6aspect9SimulatorILi2EE14solve_timestepEv + 187: 9   aspect                              0x000000010e8286ab _ZN6aspect9SimulatorILi2EE14solve_timestepEv
#8  10  aspect                              0x000000010e8275d4 _ZN6aspect9SimulatorILi2EE3runEv + 1204: 10  aspect                              0x000000010e8275d4 _ZN6aspect9SimulatorILi2EE3runEv
#9  11  aspect                              0x000000010f3243cc _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb + 588: 11  aspect                              0x000000010f3243cc _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
#10  12  aspect                              0x000000010f323be9 main + 1097: 12  aspect                              0x000000010f323be9 main
#11  13  libdyld.dylib                       0x00007fff20957f3d start + 1: 13  libdyld.dylib                       0x00007fff20957f3d start
#12  14  ???                                 0x0000000000000002 0x0 + 2: 14  ???                                 0x0000000000000002 0x0

[gassmoeller]

I would like to avoid having material models check for all possible new features we introduce in the main code (it may be a growing list). It is equally unsustainable to keep a list of all models that support elasticity in the main code (it removes the idea of plugins that are independent of the main code). I think the right approach is to ask: What does a material model need to do to be able to support elasticity? I think the answer is it needs to create a certain additional material model outputs object (ElasticAdditionalOutputs) within create_additional_material_outputs (or create_additional_named_outputs I do not know the relation between the two at the moment). So if enable elasticity is set to true, we should check if the material model indeed creates such an object. Do you think that makes sense?

[anne-glerum]

Sure, that sounds good.

The ElasticAdditionalOutputs are created in create_additional_named_outputs.

I'm just wondering where to test if if enable_elasticity == true, whether ElasticAdditionalOutputs<dim> *elastic_out = out.template get_additional_output<ElasticAdditionalOutputs<dim>>() is not a null pointer. When parsing Enable elasticity, it's probably too early as the MaterialModelOutputs needs to be set up.

At the moment, rheology/elasticity.cc, viscoelastic.cc and latent_heat.cc check whether Enable elasticity is true or false. These checks could then be removed.

[bangerth]

I think after calling material_model->initialize() here would be a good place: https://github.com/geodynamics/aspect/blob/main/source/simulator/core.cc#L300-L303C35

[anne-glerum]

Addressed by #5218