diff --git a/source/particle/property/crust_and_lithosphere_formation.cc b/source/particle/property/crust_and_lithosphere_formation.cc index 143cd8519e8..e4c96e153a9 100644 --- a/source/particle/property/crust_and_lithosphere_formation.cc +++ b/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 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 requires the compositional field 'harzburgite' " + "to be mapped to the particle property with the same name.")); } diff --git a/source/particle/property/elastic_stress.cc b/source/particle/property/elastic_stress.cc index 62a228b1171..66f47c49749 100644 --- a/source/particle/property/elastic_stress.cc +++ b/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> 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 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 all_field_indices(this->n_compositional_fields()); std::iota (std::begin(all_field_indices), std::end(all_field_indices), 0); diff --git a/source/particle/property/grain_size.cc b/source/particle/property/grain_size.cc index ac90a06535d..d63d1509c1f 100644 --- a/source/particle/property/grain_size.cc +++ b/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 requires the compositional field 'grain_size' " + "to be mapped to the particle property with the same name.")); } diff --git a/source/particle/property/initial_composition.cc b/source/particle/property/initial_composition.cc index 8c452052041..76e79b74769 100644 --- a/source/particle/property/initial_composition.cc +++ b/source/particle/property/initial_composition.cc @@ -32,8 +32,14 @@ namespace aspect InitialComposition::initialize_one_particle_property(const Point &position, std::vector &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> InitialComposition::get_property_information() const { - AssertThrow(this->n_compositional_fields() > 0, - ExcMessage("You have requested the particle property , but the number of compositional fields is 0. " - "Please add compositional fields to your model, or remove " - "this particle property.")); - std::vector> 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 , 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; } } diff --git a/source/particle/property/viscoplastic_strain_invariants.cc b/source/particle/property/viscoplastic_strain_invariants.cc index 2c23c5a42f9..787c88563f9 100644 --- a/source/particle/property/viscoplastic_strain_invariants.cc +++ b/source/particle/property/viscoplastic_strain_invariants.cc @@ -52,17 +52,16 @@ namespace aspect material_inputs = MaterialModel::MaterialModelInputs(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 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 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, diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc index 25e9bf6f577..920a59bd8f4 100644 --- a/source/simulator/initial_conditions.cc +++ b/source/simulator/initial_conditions.cc @@ -294,40 +294,20 @@ namespace aspect for (unsigned int advection_field=0; advection_field 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 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::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); } } diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc index 4b948299999..5d55ae6921d 100644 --- a/source/simulator/parameters.cc +++ b/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)