Merge pull request #5920 from MFraters/further_multipl_particle_world…

…_infrastructre

Add further multiple particle world instrastructure.

include/aspect/simulator.h

@@ -1980,7 +1980,7 @@ namespace aspect
       /**
        * The world holding the particles
        */
-      std::unique_ptr<Particle::World<dim>> particle_world;
+      std::vector<std::unique_ptr<Particle::World<dim>>> particle_worlds;
 
       /**
        * @}

include/aspect/simulator_access.h

@@ -937,20 +937,20 @@ namespace aspect
       n_particle_worlds() const;
 
       /**
-       * Returns a const reference to the particle world, in case anyone
+       * Returns a const reference to a single particle world, in case anyone
        * wants to query something about particles.
        */
       const Particle::World<dim> &
-      get_particle_world() const;
+      get_particle_world(const unsigned int particle_world_index) const;
 
       /**
-       * Returns a reference to the particle world, in case anyone wants to
+       * Returns a reference to a single particle world, in case anyone wants to
        * change something within the particle world. Use with care, usually
        * you want to only let the functions within the particle subsystem
        * change member variables of the particle world.
        */
       Particle::World<dim> &
-      get_particle_world();
+      get_particle_world(const unsigned int particle_world_index);
 
       /**
        * Return true if using the block GMG Stokes solver.

source/material_model/rheology/strain_dependent.cc

@@ -412,7 +412,7 @@ namespace aspect
         // Currently this functionality only works in field composition
         if (healing_mechanism != no_healing && this->n_particle_worlds() > 0)
           {
-            const Particle::Property::Manager<dim> &particle_property_manager = this->get_particle_world().get_property_manager();
+            const Particle::Property::Manager<dim> &particle_property_manager = this->get_particle_world(0).get_property_manager();
             AssertThrow(particle_property_manager.plugin_name_exists("viscoplastic strain invariants") == false, ExcMessage("This healing mechanism currently does not work if the strain is tracked on particles."));
           }
 

source/mesh_refinement/particle_density.cc

@@ -36,7 +36,7 @@ namespace aspect
                              "postprocessor plugin `particles' to be selected. Please activate the "
                              "particles or deactivate this mesh refinement plugin."));
 
-      const Particle::ParticleHandler<dim> &particle_handler = this->get_particle_world().get_particle_handler();
+      const Particle::ParticleHandler<dim> &particle_handler = this->get_particle_world(0).get_particle_handler();
 
       for (const auto &cell : this->get_dof_handler().active_cell_iterators())
         if (cell->is_locally_owned())

source/particle/integrator/rk_2.cc

@@ -41,7 +41,7 @@ namespace aspect
       void
       RK2<dim>::initialize ()
       {
-        const auto &property_information = this->get_particle_world().get_property_manager().get_data_info();
+        const auto &property_information = this->get_particle_world(this->get_particle_world_index()).get_property_manager().get_data_info();
         property_index_old_location = property_information.get_position_by_field_name("internal: integrator properties");
       }
 

source/particle/integrator/rk_4.cc

@@ -41,7 +41,7 @@ namespace aspect
       void
       RK4<dim>::initialize ()
       {
-        const auto &property_information = this->get_particle_world().get_property_manager().get_data_info();
+        const auto &property_information = this->get_particle_world(this->get_particle_world_index()).get_property_manager().get_data_info();
 
         property_indices[0] = property_information.get_position_by_field_name("internal: integrator properties");
         property_indices[1] = property_indices[0] + dim;

source/particle/interpolator/bilinear_least_squares.cc

@@ -281,7 +281,7 @@ namespace aspect
         {
           prm.enter_subsection("Bilinear least squares");
           {
-            const auto &particle_property_information = this->get_particle_world().get_property_manager().get_data_info();
+            const auto &particle_property_information = this->get_particle_world(this->get_particle_world_index()).get_property_manager().get_data_info();
             const unsigned int n_property_components = particle_property_information.n_components();
             const unsigned int n_internal_components = particle_property_information.get_components_by_field_name("internal: integrator properties");
 

source/particle/interpolator/quadratic_least_squares.cc

@@ -558,7 +558,7 @@ namespace aspect
         {
           prm.enter_subsection("Quadratic least squares");
           {
-            const auto &particle_property_information = this->get_particle_world().get_property_manager().get_data_info();
+            const auto &particle_property_information = this->get_particle_world(this->get_particle_world_index()).get_property_manager().get_data_info();
             const unsigned int n_property_components = particle_property_information.n_components();
             const unsigned int n_internal_components = particle_property_information.get_components_by_field_name("internal: integrator properties");
 

source/particle/property/cpo_bingham_average.cc

@@ -35,7 +35,7 @@ namespace aspect
       {
         const unsigned int my_rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
         this->random_number_generator.seed(random_number_seed+my_rank);
-        const auto &manager = this->get_particle_world().get_property_manager();
+        const auto &manager = this->get_particle_world(this->get_particle_world_index()).get_property_manager();
         AssertThrow(manager.plugin_name_exists("crystal preferred orientation"),
                     ExcMessage("No crystal preferred orientation property plugin found."));
 
@@ -252,7 +252,7 @@ namespace aspect
         {
           // Get a pointer to the CPO particle property.
           cpo_particle_property = std::make_unique<const Particle::Property::CrystalPreferredOrientation<dim>> (
-                                    this->get_particle_world().get_property_manager().template get_matching_property<Particle::Property::CrystalPreferredOrientation<dim>>());
+                                    this->get_particle_world(this->get_particle_world_index()).get_property_manager().template get_matching_property<Particle::Property::CrystalPreferredOrientation<dim>>());
 
           random_number_seed = prm.get_integer ("Random number seed");
           n_grains = cpo_particle_property->get_number_of_grains();

source/particle/property/cpo_elastic_tensor.cc

@@ -66,7 +66,7 @@ namespace aspect
       void
       CpoElasticTensor<dim>::initialize ()
       {
-        const auto &manager = this->get_particle_world().get_property_manager();
+        const auto &manager = this->get_particle_world(this->get_particle_world_index()).get_property_manager();
         AssertThrow(manager.plugin_name_exists("crystal preferred orientation"),
                     ExcMessage("No crystal preferred orientation property plugin found."));
 
@@ -147,7 +147,7 @@ namespace aspect
       {
         // Get a reference to the CPO particle property.
         const Particle::Property::CrystalPreferredOrientation<dim> &cpo_particle_property =
-          this->get_particle_world().get_property_manager().template get_matching_property<Particle::Property::CrystalPreferredOrientation<dim>>();
+          this->get_particle_world(this->get_particle_world_index()).get_property_manager().template get_matching_property<Particle::Property::CrystalPreferredOrientation<dim>>();
 
 
         const SymmetricTensor<2,6> C_average = voigt_average_elastic_tensor(cpo_particle_property,

source/particle/property/elastic_tensor_decomposition.cc

@@ -270,7 +270,7 @@ namespace aspect
       void
       ElasticTensorDecomposition<dim>::initialize ()
       {
-        const Particle::Property::Manager<dim> &manager = this->get_particle_world().get_property_manager();
+        const Particle::Property::Manager<dim> &manager = this->get_particle_world(this->get_particle_world_index()).get_property_manager();
         AssertThrow(manager.plugin_name_exists("crystal preferred orientation"),
                     ExcMessage("No cpo property plugin found."));
         AssertThrow(manager.plugin_name_exists("cpo elastic tensor"),

source/postprocess/crystal_preferred_orientation.cc

@@ -173,8 +173,8 @@ namespace aspect
     CrystalPreferredOrientation<dim>::execute (TableHandler &statistics)
     {
 
-      const Particle::Property::Manager<dim> &manager = this->get_particle_world().get_property_manager();
-      const Particle::Property::ParticleHandler<dim> &particle_handler = this->get_particle_world().get_particle_handler();
+      const Particle::Property::Manager<dim> &manager = this->get_particle_world(0).get_property_manager();
+      const Particle::Property::ParticleHandler<dim> &particle_handler = this->get_particle_world(0).get_particle_handler();
 
       const bool cpo_elastic_decomposition_plugin_exists = manager.plugin_name_exists("elastic tensor decomposition");
 
@@ -220,7 +220,7 @@ namespace aspect
                                                                         + "isotropic_norm_square") : "") << std::endl;
 
       // get particle data
-      const Particle::Property::ParticlePropertyInformation &property_information = this->get_particle_world().get_property_manager().get_data_info();
+      const Particle::Property::ParticlePropertyInformation &property_information = this->get_particle_world(0).get_property_manager().get_data_info();
 
       AssertThrow(property_information.fieldname_exists("cpo mineral 0 type") ,
                   ExcMessage("No CPO particle properties found. Make sure that the CPO particle property plugin is selected."));

source/postprocess/load_balance_statistics.cc

@@ -45,7 +45,7 @@ namespace aspect
 
       if (this->n_particle_worlds() > 0)
         {
-          const unsigned int locally_owned_particles = this->get_particle_world().
+          const unsigned int locally_owned_particles = this->get_particle_world(0).
                                                        get_particle_handler().n_locally_owned_particles();
           const dealii::Utilities::MPI::MinMaxAvg particles_per_process =
             dealii::Utilities::MPI::min_max_avg(locally_owned_particles,this->get_mpi_communicator());

source/postprocess/particle_count_statistics.cc

@@ -35,11 +35,11 @@ namespace aspect
     ParticleCountStatistics<dim>::execute (TableHandler &statistics)
     {
       const Particle::ParticleHandler<dim> &particle_handler =
-        this->get_particle_world().get_particle_handler();
+        this->get_particle_world(0).get_particle_handler();
 
       unsigned int local_min_particles = std::numeric_limits<unsigned int>::max();
       unsigned int local_max_particles = 0;
-      const Particle::types::particle_index global_particles = this->get_particle_world().n_global_particles();
+      const Particle::types::particle_index global_particles = this->get_particle_world(0).n_global_particles();
 
       // compute local min/max
       for (const auto &cell : this->get_dof_handler().active_cell_iterators())

source/postprocess/particles.cc

@@ -323,7 +323,7 @@ namespace aspect
       if (std::isnan(last_output_time))
         last_output_time = this->get_time() - output_interval;
 
-      const Particle::World<dim> &world = this->get_particle_world();
+      const Particle::World<dim> &world = this->get_particle_world(0);
 
       statistics.add_value("Number of advected particles",world.n_global_particles());
 
@@ -562,7 +562,7 @@ namespace aspect
       std::ostringstream os;
       aspect::oarchive oa (os);
 
-      this->get_particle_world().save(os);
+      this->get_particle_world(0).save(os);
       oa << (*this);
 
       status_strings["Particles"] = os.str();
@@ -580,7 +580,7 @@ namespace aspect
           aspect::iarchive ia (is);
 
           // Load the particle world
-          this->get_particle_world().load(is);
+          this->get_particle_world(0).load(is);
 
           ia >> (*this);
         }

source/postprocess/visualization/particle_count.cc

@@ -43,7 +43,7 @@ namespace aspect
       ParticleCount<dim>::execute() const
       {
         const Particle::ParticleHandler<dim> &particle_handler =
-          this->get_particle_world().get_particle_handler();
+          this->get_particle_world(0).get_particle_handler();
 
         std::pair<std::string, std::unique_ptr<Vector<float>>>
         return_value ("particles_per_cell",

source/simulator/core.cc

@@ -189,7 +189,7 @@ namespace aspect
                    nullptr),
 #endif
     boundary_heat_flux (BoundaryHeatFlux::create_boundary_heat_flux<dim>(prm)),
-    particle_world(nullptr),
+    particle_worlds(),
     time (numbers::signaling_nan<double>()),
     time_step (numbers::signaling_nan<double>()),
     old_time_step (numbers::signaling_nan<double>()),
@@ -448,12 +448,16 @@ namespace aspect
 
     if (postprocess_manager.template has_matching_postprocessor<Postprocess::Particles<dim>>())
       {
-        particle_world = std::make_unique<Particle::World<dim>>();
-        if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(particle_world.get()))
-          sim->initialize_simulator (*this);
+        particle_worlds.emplace_back(std::move(std::make_unique<Particle::World<dim>>()));
+        for (unsigned int particle_world_index = 0 ; particle_world_index < particle_worlds.size(); ++particle_world_index)
+          {
+            if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(particle_worlds[particle_world_index].get()))
+              sim->initialize_simulator (*this);
+
+            particle_worlds.back()->parse_parameters(prm,particle_world_index);
+            particle_worlds.back()->initialize();
 
-        particle_world->parse_parameters(prm,0);
-        particle_world->initialize();
+          }
       }
 
     mesh_refinement_manager.initialize_simulator (*this);
@@ -568,7 +572,9 @@ namespace aspect
     // is destroyed after the latter. But it stores a pointer to the
     // triangulation and uses it during destruction. This results in
     // trouble. So destroy it first.
-    particle_world.reset();
+
+    for (auto &particle_world : particle_worlds)
+      particle_world.reset();
 
     // wait if there is a thread that's still writing the statistics
     // object (set from the output_statistics() function)
@@ -608,9 +614,10 @@ namespace aspect
 
     // Copy particle handler to restore particle location and properties
     // before repeating a timestep
-    if (particle_world.get() != nullptr)
+    for (auto &particle_world : particle_worlds)
       particle_world->backup_particles();
 
+
     // then interpolate the current boundary velocities. copy constraints
     // into current_constraints and then add to current_constraints
     compute_current_constraints ();
@@ -642,7 +649,7 @@ namespace aspect
     if (prescribed_stokes_solution.get())
       prescribed_stokes_solution->update();
 
-    if (particle_world.get() != nullptr)
+    for (auto &particle_world : particle_worlds)
       particle_world->update();
 
     // do the same for the traction boundary conditions and other things
@@ -2133,7 +2140,7 @@ namespace aspect
             // Restore particles through stored copy of particle handler,
             // created in start_timestep(),
             // but only if this timestep is to be repeated.
-            if (particle_world.get() != nullptr)
+            for (auto &particle_world : particle_worlds)
               particle_world->restore_particles();
 
             continue; // repeat time step loop

source/simulator/initial_conditions.cc

@@ -284,8 +284,8 @@ 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_world->get_interpolator();
-    const Particle::Property::Manager<dim> &particle_property_manager = particle_world->get_property_manager();
+    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();
 
     std::vector<unsigned int> particle_property_indices;
     ComponentMask property_mask  (particle_property_manager.get_data_info().n_components(),false);
@@ -356,7 +356,7 @@ namespace aspect
           try
             {
               particle_properties =
-                particle_interpolator.properties_at_points(particle_world->get_particle_handler(),
+                particle_interpolator.properties_at_points(particle_worlds[0]->get_particle_handler(),
                                                            quadrature_points,
                                                            property_mask,
                                                            cell);

source/simulator/simulator_access.cc

@@ -821,30 +821,27 @@ namespace aspect
   unsigned int
   SimulatorAccess<dim>::n_particle_worlds() const
   {
-    // operator () returns whether an object is managed by a unique ptr
-    return (simulator->particle_world ? 1 : 0);
+    return simulator->particle_worlds.size();
   }
 
 
 
   template <int dim>
   const Particle::World<dim> &
-  SimulatorAccess<dim>::get_particle_world() const
+  SimulatorAccess<dim>::get_particle_world(unsigned int particle_world_index) const
   {
-    Assert (simulator->particle_world.get() != nullptr,
-            ExcMessage("You can not call this function if there is no particle world."));
-    return *simulator->particle_world.get();
+    AssertThrow (particle_world_index < simulator->particle_worlds.size(), ExcInternalError());
+    return *simulator->particle_worlds[particle_world_index].get();
   }
 
 
 
   template <int dim>
   Particle::World<dim> &
-  SimulatorAccess<dim>::get_particle_world()
+  SimulatorAccess<dim>::get_particle_world(unsigned int particle_world_index)
   {
-    Assert (simulator->particle_world.get() != nullptr,
-            ExcMessage("You can not call this function if there is no particle world."));
-    return *simulator->particle_world.get();
+    AssertThrow (particle_world_index < simulator->particle_worlds.size(), ExcInternalError());
+    return *simulator->particle_worlds[particle_world_index].get();
   }
 
 

source/simulator/solver_schemes.cc

@@ -226,7 +226,7 @@ namespace aspect
   {
     // Advect the particles before they are potentially used to
     // set up the compositional fields.
-    if (particle_world.get() != nullptr)
+    for (auto &particle_world : particle_worlds)
       {
         // Do not advect the particles in the initial refinement stage
         const bool in_initial_refinement = (timestep_number == 0)
@@ -958,8 +958,10 @@ namespace aspect
         // Restore particles through stored copy of particle handler,
         // but only if they have already been displaced in a nonlinear
         // iteration (in the assemble_and_solve_composition call).
-        if ((particle_world.get() != nullptr) && (nonlinear_iteration > 0))
-          particle_world->restore_particles();
+
+        if (nonlinear_iteration > 0)
+          for (auto &particle_world : particle_worlds)
+            particle_world->restore_particles();
 
         const double relative_temperature_residual =
           assemble_and_solve_temperature(initial_temperature_residual,
@@ -1115,8 +1117,10 @@ namespace aspect
         // Restore particles through stored copy of particle handler,
         // but only if they have already been displaced in a nonlinear
         // iteration (in the assemble_and_solve_composition call).
-        if ((particle_world.get() != nullptr) && (nonlinear_iteration > 0))
-          particle_world->restore_particles();
+
+        if (nonlinear_iteration > 0)
+          for (auto &particle_world : particle_worlds)
+            particle_world->restore_particles();
 
         const double relative_temperature_residual =
           assemble_and_solve_temperature(initial_temperature_residual,

tests/mesh_deformation_particles.cc

@@ -34,7 +34,7 @@ void post_mesh_deformation (const SimulatorAccess<dim> &simulator_access)
   // Get the reference location of the one particle,
   // and check that it has the correct value.
   const Point<dim> predicted_particle_reference_location =
-    simulator_access.get_particle_world().get_particle_handler().begin()->get_reference_location();
+    simulator_access.get_particle_world(0).get_particle_handler().begin()->get_reference_location();
   // At timestep 0, the mesh is not deformed
   // other than through the initial topography/
   // initial mesh deformation plugins, which

tests/particle_handler_copy.cc

@@ -33,7 +33,7 @@ void duplicate_particle_handler(const SimulatorAccess<dim> &simulator_access,
                                 const unsigned int,
                                 const SolverControl &)
 {
-  const auto &particle_world = simulator_access.get_particle_world();
+  const auto &particle_world = simulator_access.get_particle_world(0);
   dealii::Particles::ParticleHandler<dim> particle_handler;
   std::cout << "duplicating particle handler" << std::endl;