Support interpolation + repartition for worlds

source/particle/world.cc

@@ -52,20 +52,6 @@ namespace aspect
  World<dim>::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<dim>::cell_iterator &cell,
- const CellStatus status)
- -> unsigned int
-#else
- [&] (const typename parallel::distributed::Triangulation<dim>::cell_iterator &cell,
- const typename parallel::distributed::Triangulation<dim>::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<dim>::cell_iterator &cell,
+ const CellStatus status)
+ -> unsigned int
+#else
+ [ &, world_index] (const typename parallel::distributed::Triangulation<dim>::cell_iterator &cell,
+ const typename parallel::distributed::Triangulation<dim>::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");
 

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."));
  }
  }
 

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<dim> &particle_interpolator = particle_worlds[0]->get_interpolator();
- const Particle::Property::Manager<dim> &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<std::vector<std::pair<unsigned int, unsigned int>>> particle_property_indices;
+ // A property component mask indicating for each particle world which particle properties need to be interpolated
+ std::vector<ComponentMask> property_mask;
 
- std::vector<unsigned int> 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<std::string,unsigned int> 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<dim> &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<std::pair<unsigned int, unsigned int>>());
+ 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<advection_fields.size(); ++advection_field)
  {
- const unsigned int particle_property_index = std::count(introspection.compositional_field_methods.begin(),
- introspection.compositional_field_methods.begin() + 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"));
-
- 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<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 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<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"));
+
+ 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<advection_fields.size(); ++j)
+ for (unsigned int j=0; j<particle_property_indices[world_index].size(); ++j)
  for (unsigned int i=0; i<n_dofs_per_cell; ++i)
  {
  const unsigned int system_local_dof
- = finite_element.component_to_system_index(advection_fields[j].component_index(introspection),
+ = finite_element.component_to_system_index(advection_fields[particle_property_indices[world_index][j].first].component_index(introspection),
  /*dof index within component=*/i);
 
- particle_solution(local_dof_indices[system_local_dof]) = particle_properties[i][particle_property_indices[j]];
+ particle_solution(local_dof_indices[system_local_dof]) = particle_properties[i][particle_property_indices[world_index][j].second];
  }
  }
  }