[WIP] Make sure that the compositional fields are mapped to the correct particle properties

source/particle/property/crust_and_lithosphere_formation.cc

@@ -52,6 +52,17 @@ namespace aspect
 
         basalt_index = this->introspection().compositional_index_for_name("basalt");
         harzburgite_index = this->introspection().compositional_index_for_name("harzburgite");
+
+        // Make sure that the compositional fields are mapped to the particle properties with the same names
+        const auto basalt_it = this->get_parameters().mapped_particle_properties.find(basalt_index);
+        AssertThrow(basalt_it != this->get_parameters().mapped_particle_properties.end() && basalt_it->second.first == "basalt",
+                    ExcMessage("Particle property plugin <crust and lithosphere formation> requires the compositional field 'basalt' "
+                               "to be mapped to the particle property with the same name."));
+
+        const auto harzburgite_it = this->get_parameters().mapped_particle_properties.find(harzburgite_index);
+        AssertThrow(harzburgite_it != this->get_parameters().mapped_particle_properties.end() && harzburgite_it->second.first == "harzburgite",
+                    ExcMessage("Particle property plugin <crust and lithosphere formation> requires the compositional field 'harzburgite' "
+                               "to be mapped to the particle property with the same name."));
       }
 
 

source/particle/property/elastic_stress.cc

@@ -68,6 +68,18 @@ namespace aspect
         AssertThrow((stress_field_indices.size() == 2*SymmetricTensor<2,dim>::n_independent_components),
                     ExcMessage("The number of stress tensor element fields in the 'elastic stress' plugin does not equal twice the number of independent components."));
 
+        // Make sure that the stress components are mapped to the corresponding particle properties
+        const std::vector<std::pair<std::string, unsigned int>> property_information = this->get_property_information();
+        for (unsigned int i = 0; i < 2 * SymmetricTensor<2, dim>::n_independent_components; ++i)
+          {
+            const auto it = this->get_parameters().mapped_particle_properties.find(stress_field_indices[i]);
+            AssertThrow(it != this->get_parameters().mapped_particle_properties.end() &&
+                        it->second.first == property_information[i].first && it->second.second == 0,
+                        ExcMessage("If the particle property <elastic stress> is requested, then there must be 3+3 (or 6+6) compositional fields "
+                                   "mapped to particle properties with the name ve_stress_* in 2D (or 3D). Please read the documentation of "
+                                   "this particle property for details."));
+          }
+
         // Get the indices of all compositions that do not correspond to stress tensor elements.
         std::vector<unsigned int> all_field_indices(this->n_compositional_fields());
         std::iota (std::begin(all_field_indices), std::end(all_field_indices), 0);

source/particle/property/grain_size.cc

@@ -51,6 +51,12 @@ namespace aspect
                                "there is a compositional field named 'grain_size'."));
 
         grain_size_index = this->introspection().compositional_index_for_name("grain_size");
+
+        // Make sure that the compositional field is mapped to the particle property with the same name
+        const auto it = this->get_parameters().mapped_particle_properties.find(grain_size_index);
+        AssertThrow(it != this->get_parameters().mapped_particle_properties.end() && it->second.first == "grain_size",
+                    ExcMessage("Particle property plugin <grain size> requires the compositional field 'grain_size' "
+                               "to be mapped to the particle property with the same name."));
       }
 
 

source/particle/property/initial_composition.cc

@@ -32,8 +32,14 @@ namespace aspect
       InitialComposition<dim>::initialize_one_particle_property(const Point<dim> &position,
                                                                 std::vector<double> &data) const
       {
-        for (unsigned int i = 0; i < this->n_compositional_fields(); ++i)
-          data.push_back(this->get_initial_composition_manager().initial_composition(position,i));
+        const std::string prefix = "initial ";
+        for (const auto &key_and_value : this->get_parameters().mapped_particle_properties)
+          {
+            const std::string &property_name = key_and_value.second.first;
+            if (property_name.size() > prefix.size() &&
+                property_name.compare(0, prefix.size(), prefix) == 0)
+              data.push_back(this->get_initial_composition_manager().initial_composition(position, key_and_value.first));
+          }
       }
 
 
@@ -63,21 +69,28 @@ namespace aspect
       std::vector<std::pair<std::string, unsigned int>>
       InitialComposition<dim>::get_property_information() const
       {
-        AssertThrow(this->n_compositional_fields() > 0,
-                    ExcMessage("You have requested the particle property <initial "
-                               "composition>, but the number of compositional fields is 0. "
-                               "Please add compositional fields to your model, or remove "
-                               "this particle property."));
-
         std::vector<std::pair<std::string,unsigned int>> property_information;
 
-        for (unsigned int i = 0; i < this->n_compositional_fields(); ++i)
+        // Find the compositional fields that are mapped to this particle property
+        const std::string prefix = "initial ";
+        for (const auto &key_and_value : this->get_parameters().mapped_particle_properties)
           {
-            std::ostringstream field_name;
-            field_name << "initial " << this->introspection().name_for_compositional_index(i);
-            property_information.emplace_back(field_name.str(),1);
+            const std::string &property_name = key_and_value.second.first;
+            if (property_name.size() > prefix.size() &&
+                property_name.compare(0, prefix.size(), prefix) == 0)
+              {
+                AssertThrow(key_and_value.second.second == 0, ExcNotImplemented());
+                property_information.emplace_back(key_and_value.second.first, 1);
+              }
           }
 
+        AssertThrow(property_information.size() > 0,
+                    ExcMessage("You have requested the particle property <initial "
+                               "composition>, but there is no compositional field "
+                               "mapped to this particle property. Please prefix the "
+                               "property names of the corresponding compositional fields "
+                               "by 'initial ', or remove this particle property."));
+
         return property_information;
       }
     }

source/particle/property/viscoplastic_strain_invariants.cc

@@ -52,17 +52,16 @@ namespace aspect
         material_inputs = MaterialModel::MaterialModelInputs<dim>(1,this->n_compositional_fields());
 
         // Find out which fields are used.
-        if (this->introspection().compositional_name_exists("plastic_strain"))
-          n_components += 1;
-
-        if (this->introspection().compositional_name_exists("viscous_strain"))
-          n_components += 1;
-
-        if (this->introspection().compositional_name_exists("total_strain"))
-          n_components += 1;
-
-        if (this->introspection().compositional_name_exists("noninitial_plastic_strain"))
-          n_components += 1;
+        const std::vector<std::string> possible_names = {"plastic_strain", "viscous_strain", "total_strain", "noninitial_plastic_strain"};
+        for (const auto &name : possible_names)
+          {
+            // Make sure that the compositional field is mapped to the particle property with the same name
+            const auto it = this->get_parameters().mapped_particle_properties.find(this->introspection().compositional_index_for_name(name));
+            AssertThrow(it != this->get_parameters().mapped_particle_properties.end() && it->second.first == name,
+                        ExcMessage("Particle property plugin <viscoplastic strain invariants> requires the compositional field '"
+                                   + name + "' to be mapped to the particle property with the same name."));
+            n_components += 1;
+          }
 
         if (n_components == 0)
           AssertThrow(false,

source/simulator/initial_conditions.cc

@@ -294,40 +294,20 @@ namespace aspect
 
         for (unsigned int advection_field=0; advection_field<advection_fields.size(); ++advection_field)
           {
-            if (parameters.mapped_particle_properties.size() != 0)
-              {
-                const std::pair<std::string,unsigned int> particle_property_and_component = parameters.mapped_particle_properties.find(advection_fields[advection_field].compositional_variable)->second;
-
-                // Check if the required particle property exists in the current particle manager.
-                // If not: assume we find it in another world.
-                if (particle_property_manager.get_data_info().fieldname_exists(particle_property_and_component.first))
-                  {
-                    Assert (advection_field_has_been_found[advection_field] == false,
-                            ExcMessage("The field " + advection_fields[advection_field].name(introspection) + " is mapped to particle properties in more than one particle manager. This is not supported."));
+            const std::pair<std::string,unsigned int> particle_property_and_component = parameters.mapped_particle_properties.find(advection_fields[advection_field].compositional_variable)->second;
 
-                    const unsigned int particle_property_index = particle_property_manager.get_data_info().get_position_by_field_name(particle_property_and_component.first)
-                                                                 + particle_property_and_component.second;
-
-                    advection_field_has_been_found[advection_field] = true;
-                    particle_property_indices[particle_manager].emplace_back(advection_field, particle_property_index);
-                    property_mask[particle_manager].set(particle_property_index,true);
-                  }
-              }
-            else
+            // Check if the required particle property exists in the current particle manager.
+            // If not: assume we find it in another world.
+            if (particle_property_manager.get_data_info().fieldname_exists(particle_property_and_component.first))
               {
-                Assert(particle_managers.size() == 1,
-                       ExcMessage("Automatically mapping particle properties to compositional fields is only supported if there is exactly one set of particles. "
-                                  "Please specify the particle properties manually in the parameter file using the parameter 'Compositional Fields/Mapped particle properties'."));
+                Assert (advection_field_has_been_found[advection_field] == false,
+                        ExcMessage("The field " + advection_fields[advection_field].name(introspection) + " is mapped to particle properties in more than one particle manager. This is not supported."));
 
-                const unsigned int particle_property_index = std::count(introspection.compositional_field_methods.begin(),
-                                                                        introspection.compositional_field_methods.begin() + advection_fields[advection_field].compositional_variable,
-                                                                        Parameters<dim>::AdvectionFieldMethod::particles);
-                AssertThrow(particle_property_index <= particle_property_manager.get_data_info().n_components(),
-                            ExcMessage("Can not automatically match particle properties to fields, because there are"
-                                       "more fields that are marked as particle advected than particle properties"));
+                const unsigned int particle_property_index = particle_property_manager.get_data_info().get_position_by_field_name(particle_property_and_component.first)
+                                                             + particle_property_and_component.second;
 
                 advection_field_has_been_found[advection_field] = true;
-                particle_property_indices[particle_manager].emplace_back(advection_field,particle_property_index);
+                particle_property_indices[particle_manager].emplace_back(advection_field, particle_property_index);
                 property_mask[particle_manager].set(particle_property_index,true);
               }
           }

source/simulator/parameters.cc

@@ -2246,9 +2246,8 @@ namespace aspect
       const unsigned int number_of_particle_fields =
         std::count(compositional_field_methods.begin(),compositional_field_methods.end(),AdvectionFieldMethod::particles);
 
-      AssertThrow ((x_mapped_particle_properties.size() == number_of_particle_fields)
-                   || (x_mapped_particle_properties.size() == 0),
-                   ExcMessage ("The list of names for the mapped particle property fields needs to either be empty or have a length equal to "
+      AssertThrow (x_mapped_particle_properties.size() == number_of_particle_fields,
+                   ExcMessage ("The list of names for the mapped particle property fields needs to have a length equal to "
                                "the number of compositional fields that are interpolated from particle properties."));
 
       if (std::find(compositional_field_methods.begin(), compositional_field_methods.end(), AdvectionFieldMethod::fem_darcy_field)