From c510ad650e69f8ada80367e63f557859cb92330e Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Fri, 20 Sep 2024 17:04:39 +0200 Subject: [PATCH] prepare property lookup --- source/simulator/initial_conditions.cc | 79 +++++++++++++++----------- 1 file changed, 47 insertions(+), 32 deletions(-) diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc index e9d4abbe266..cfb0a254181 100644 --- a/source/simulator/initial_conditions.cc +++ b/source/simulator/initial_conditions.cc @@ -284,36 +284,51 @@ namespace aspect // need to write into it and we can not // write into vectors with ghost elements - const Particle::Interpolator::Interface &particle_interpolator = particle_worlds[0]->get_interpolator(); - const Particle::Property::Manager &particle_property_manager = particle_worlds[0]->get_property_manager(); + // For each particle world store the pairs of corresponding advection field index (first) and particle property index (second) + std::vector>> particle_property_indices; + // A property component mask indicating for each particle world which particle properties need to be interpolated + std::vector property_mask; - std::vector particle_property_indices; - ComponentMask property_mask (particle_property_manager.get_data_info().n_components(),false); - - for (const auto &advection_field: advection_fields) + for (unsigned int world_index = 0; world_index < particle_worlds.size(); ++world_index) { - if (parameters.mapped_particle_properties.size() != 0) - { - const std::pair particle_property_and_component = parameters.mapped_particle_properties.find(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; + const Particle::Property::Manager &particle_property_manager = particle_worlds[world_index]->get_property_manager(); - particle_property_indices.push_back(particle_property_index); + particle_property_indices.push_back(std::vector>()); + property_mask.push_back(ComponentMask(particle_property_manager.get_data_info().n_components(),false)); - property_mask.set(particle_property_index,true); - } - else + for (unsigned int advection_field=0; advection_field::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")); - - particle_property_indices.push_back(particle_property_index); - property_mask.set(particle_property_index,true); + if (parameters.mapped_particle_properties.size() != 0) + { + const std::pair 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 world. + // If not: assume we find it in another world. + if (particle_property_manager.get_data_info().fieldname_exists(particle_property_and_component.first)) + { + 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; + + particle_property_indices[world_index].push_back({advection_field, particle_property_index}); + property_mask[world_index].set(particle_property_index,true); + } + } + else + { + Assert(particle_worlds.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'.")); + + 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")); + + particle_property_indices[world_index].push_back({advection_field,particle_property_index}); + property_mask[world_index].set(particle_property_index,true); + } } } @@ -357,10 +372,10 @@ namespace aspect try { particle_properties = - particle_interpolator.properties_at_points(particle_worlds[world_index]->get_particle_handler(), - quadrature_points, - property_mask, - cell); + particle_worlds[world_index]->get_interpolator().properties_at_points(particle_worlds[world_index]->get_particle_handler(), + quadrature_points, + property_mask[world_index], + cell); } // interpolators that throw exceptions usually do not result in // anything good, because they result in an unwinding of the stack @@ -390,14 +405,14 @@ namespace aspect // to the particle field interpolated at these points cell->get_dof_indices (local_dof_indices); const unsigned int n_dofs_per_cell = finite_element.base_element(base_element_index).dofs_per_cell; - for (unsigned int j=0; j