From b25c4c9b149fe8e511e82c8982940f0b80d6057b Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Thu, 19 Sep 2024 18:33:05 +0200 Subject: [PATCH] Support interpolation + repartition for worlds --- source/particle/world.cc | 34 +++++------ source/simulator/checkpoint_restart.cc | 9 +++ source/simulator/initial_conditions.cc | 79 +++++++++++++++----------- 3 files changed, 74 insertions(+), 48 deletions(-) diff --git a/source/particle/world.cc b/source/particle/world.cc index 969b2a36c1b..655b203a4e6 100644 --- a/source/particle/world.cc +++ b/source/particle/world.cc @@ -52,20 +52,6 @@ namespace aspect World::initialize() { CitationInfo::add("particles"); - if (particle_load_balancing & ParticleLoadBalancing::repartition) - this->get_triangulation().signals.weight.connect( -#if DEAL_II_VERSION_GTE(9,6,0) - [&] (const typename parallel::distributed::Triangulation::cell_iterator &cell, - const CellStatus status) - -> unsigned int -#else - [&] (const typename parallel::distributed::Triangulation::cell_iterator &cell, - const typename parallel::distributed::Triangulation::CellStatus status) - -> unsigned int -#endif - { - return this->cell_weight(cell, status); - }); // Create a particle handler that stores the future particles. // If we restarted from a checkpoint we will fill this particle handler @@ -373,7 +359,6 @@ namespace aspect if (cell->is_active() && !cell->is_locally_owned()) return 0; - const unsigned int base_weight = 1000; unsigned int n_particles_in_cell = 0; switch (status) { @@ -408,7 +393,7 @@ namespace aspect Assert(false, ExcInternalError()); break; } - return base_weight + n_particles_in_cell * particle_weight; + return n_particles_in_cell * particle_weight; } @@ -896,6 +881,23 @@ namespace aspect "are listed in the corresponding manual subsection.")); } + if (particle_load_balancing & ParticleLoadBalancing::repartition) + this->get_triangulation().signals.weight.connect( +#if DEAL_II_VERSION_GTE(9,6,0) + [ &, world_index] (const typename parallel::distributed::Triangulation::cell_iterator &cell, + const CellStatus status) + -> unsigned int +#else + [ &, world_index] (const typename parallel::distributed::Triangulation::cell_iterator &cell, + const typename parallel::distributed::Triangulation::CellStatus status) + -> unsigned int +#endif + { + // Only add the base weight of cells in particle world 0, because all weights will be summed + // across all particle worlds. + return (world_index == 0) ? 1000 + this->cell_weight(cell, status) : this->cell_weight(cell, status); + }); + TimerOutput::Scope timer_section(this->get_computing_timer(), "Particles: Initialization"); diff --git a/source/simulator/checkpoint_restart.cc b/source/simulator/checkpoint_restart.cc index 11a81f8a98f..18ad0674bdf 100644 --- a/source/simulator/checkpoint_restart.cc +++ b/source/simulator/checkpoint_restart.cc @@ -97,6 +97,7 @@ namespace aspect oa << parameters.names_of_compositional_fields; oa << parameters.normalized_fields; oa << parameters.mesh_deformation_enabled; + oa << parameters.n_particle_worlds; } @@ -260,6 +261,14 @@ namespace aspect "These need to be the same during restarting " "from a checkpoint.")); + unsigned int n_particle_worlds; + ia >> n_particle_worlds; + AssertThrow (n_particle_worlds == parameters.n_particle_worlds, + ExcMessage ("The number of particle worlds that were stored " + "in the checkpoint file is not the same as the one " + "you currently set in your input file. " + "These need to be the same during restarting " + "from a checkpoint.")); } } 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