Out-of-bounds issue with certain types of fields

[LeeviT]

I recently added a feature to my model to read dominant phase from a lookup table file, using visco-plastic material model. That caused a crash and a non-explanatory error message specifying nothing but "out of memory". Started simplifying the model to narrow down the error source, and eventually by ditching a few chemical compositions, ASPECT produced much more informative error message about a certain composition key having wrong number of values.

Essentially the same error can be replicated by using the test https://github.com/geodynamics/aspect/blob/main/tests/visco_plastic_phases_discrete_comp.prm but modifying the compositional fields section from

subsection Compositional fields
  set Number of fields = 1
  set Names of fields = spharz
  set Compositional field methods = field
end

to

subsection Compositional fields
  set Number of fields = 2
  set Names of fields = plastic_strain, spharz
  set Types of fields = strain, chemical composition
  set Compositional field methods = field, field
end

The full error message by release build executable is

----------------------------------------------------
Exception 'ExcMessage("The key <" + field_name + "> in <"+ options.property_name + "> does not have " + "the expected number of values. It expects " + std::to_string(n_expected_values) + " or 1 values, but we found " + std::to_string(n_values) + " values.")' on rank 0 on processing: 

--------------------------------------------------------
An error occurred in line <335> of file <aspect/source/utilities.cc> in function
    void aspect::Utilities::MapParsing::{anonymous}::read_or_check_map_structure(std::multimap<std::__cxx11::basic_string<char>, double>&, aspect::Utilities::MapParsing::Options&)
The violated condition was: 
    (n_expected_values == n_values || n_values == 1)
Additional information: 
    The key <spharz> in <Prefactors for diffusion creep> does not have the
    expected number of values. It expects 226560816 or 1 values, but we
    found 2 values.

Stacktrace:
-----------
#0  aspect-dev/bin/aspect-release: aspect::Utilities::MapParsing::parse_map_to_double_array(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, aspect::Utilities::MapParsing::Options&)
#1  aspect-dev/bin/aspect-release: aspect::MaterialModel::Rheology::DiffusionCreep<2>::parse_parameters(dealii::ParameterHandler&, std::unique_ptr<std::vector<unsigned int, std::allocator<unsigned int> >, std::default_delete<std::vector<unsigned int, std::allocator<unsigned int> > > > const&)
#2  aspect-dev/bin/aspect-release: aspect::MaterialModel::Rheology::ViscoPlastic<2>::parse_parameters(dealii::ParameterHandler&, std::unique_ptr<std::vector<unsigned int, std::allocator<unsigned int> >, std::default_delete<std::vector<unsigned int, std::allocator<unsigned int> > > > const&)
#3  aspect-dev/bin/aspect-release: aspect::MaterialModel::ViscoPlastic<2>::parse_parameters(dealii::ParameterHandler&)
#4  aspect-dev/bin/aspect-release: aspect::Simulator<2>::Simulator(int, dealii::ParameterHandler&)
#5  aspect-dev/bin/aspect-release: ) [0x1bf06be]
#6  aspect-dev/bin/aspect-release: main
--------------------------------------------------------

while debug build produces only

/usr/include/c++/14/bits/stl_vector.h:1143: std::vector<_Tp, _Alloc>::reference std::vector<_Tp, _Alloc>::operator[](size_type) [with _Tp = unsigned int; _Alloc = std::allocator<unsigned int>; reference = unsigned int&; size_type = long unsigned int]: Assertion '__n < this->size()' failed.
Aborted (core dumped)

So the issue seems to be caused by addition of the strain field. If I interpret the relevant parts (at least material_model/utilities.cc) of the code right, the indices are picked from a vector that has all the fields but then later on applied on a vector, which contains only chemical compositions. Hence by adding other types of fields in the beginning of the field list causes index overflow, as chemical composition vector is at least one element shorter than the total field vector.

My original model had two strain fields and one generic in the beginning of the list, followed by nine chemical compositions, which was apparently enough to lead reaching outside of a memory block and sudden crash. I'll keep investigating this further on and trying to develop a patch.

[LeeviT]

OK, found a solution. At line 1413 in https://github.com/geodynamics/aspect/blob/main/source/material_model/utilities.cc the for-loop runs through the indices from the total field vector, i.e. in case of the test model index 1 which is field spharz, as the index 0 is reserved for plastic_strain field.

Inside the loop n_phases_per_composition vector is accessed. It has background field already set up at index 0 so everything is accessed at index i+1. Within the test model the loop however starts at i=1, so n_phases_per_composition is accessed only at i=1+1 which is already out of bounds. A simple fix would be modifying line 1413 from
for (auto i : this->introspection().chemical_composition_field_indices())
to
for (unsigned int i = 0; i < this->introspection().chemical_composition_field_indices().size(); ++i)

If this looks reasonable, I'll submit a PR.

[bangerth]

I don't know this part of ASPECT, but do you know where the value 226560816 comes from? It doesn't make sense to expect 226M entries for anything in an input file. Surely that can't be right. In the error message, that corresponds to the n_expected_values variable, which is declared as

                  const unsigned int n_expected_values = options.n_values_per_key[field_index];

If I understand you correctly, this may be an out-of-bounds access. If so (or even if you just suspect so), put the following in front of the line:

      AssertIndexRange(field_index, options.n_values_per_key.size());

This checks that the array subscript is valid. If it is not, the program will abort here with a decent error message rather than in some random downstream place with a message that doesn't make any sense.