Merge pull request #6068 from gassmoeller/test_load_balancing

Add test for load balancing and fix postprocessors

source/postprocess/load_balance_statistics.cc

@@ -45,8 +45,11 @@ namespace aspect
 
       if (this->n_particle_worlds() > 0)
         {
-          const unsigned int locally_owned_particles = this->get_particle_world(0).
-                                                       get_particle_handler().n_locally_owned_particles();
+          unsigned int locally_owned_particles = 0;
+          for (unsigned int particle_handler_index = 0; particle_handler_index < this->n_particle_worlds(); ++particle_handler_index)
+            locally_owned_particles += this->get_particle_world(particle_handler_index).
+                                       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

@@ -34,18 +34,21 @@ namespace aspect
     std::pair<std::string,std::string>
     ParticleCountStatistics<dim>::execute (TableHandler &statistics)
     {
-      const Particle::ParticleHandler<dim> &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(0).n_global_particles();
+
+      Particle::types::particle_index global_particles = 0;
+      for (unsigned int particle_handler_index = 0; particle_handler_index < this->n_particle_worlds(); ++particle_handler_index)
+        global_particles += this->get_particle_world(particle_handler_index).n_global_particles();
 
       // compute local min/max
       for (const auto &cell : this->get_dof_handler().active_cell_iterators())
         if (cell->is_locally_owned())
           {
-            const unsigned int particles_in_cell = particle_handler.n_particles_in_cell(cell);
+            unsigned int particles_in_cell = 0;
+            for (unsigned int particle_handler_index = 0; particle_handler_index < this->n_particle_worlds(); ++particle_handler_index)
+              particles_in_cell += this->get_particle_world(particle_handler_index).get_particle_handler().n_particles_in_cell(cell);
+
             local_min_particles = std::min(local_min_particles,particles_in_cell);
             local_max_particles = std::max(local_max_particles,particles_in_cell);
           }

source/postprocess/visualization/particle_count.cc

@@ -42,9 +42,6 @@ namespace aspect
       std::pair<std::string, std::unique_ptr<Vector<float>>>
       ParticleCount<dim>::execute() const
       {
-        const Particle::ParticleHandler<dim> &particle_handler =
-          this->get_particle_world(0).get_particle_handler();
-
         std::pair<std::string, std::unique_ptr<Vector<float>>>
         return_value ("particles_per_cell",
                       std::make_unique<Vector<float>>(this->get_triangulation().n_active_cells()));
@@ -53,7 +50,13 @@ namespace aspect
         for (const auto &cell : this->get_dof_handler().active_cell_iterators())
           if (cell->is_locally_owned())
             {
-              (*return_value.second)(cell->active_cell_index()) = static_cast<float> (particle_handler.n_particles_in_cell(cell));
+              unsigned int n_particles_in_cell = 0;
+              for (unsigned int particle_handler_index = 0;
+                   particle_handler_index < this->n_particle_worlds();
+                   ++particle_handler_index)
+                n_particles_in_cell += this->get_particle_world(particle_handler_index).get_particle_handler().n_particles_in_cell(cell);
+
+              (*return_value.second)(cell->active_cell_index()) = static_cast<float>(n_particles_in_cell);
             }
 
         return return_value;

tests/particle_load_balancing_repartition_multiple_worlds.prm

@@ -0,0 +1,147 @@
+# A test for the particle load balancing strategy 'repartition'
+# when multiple particle worlds are active.
+# The test is balanced to show that the cell weight used for
+# the repartition algorithm is exactly 1000 and only added once
+# from the first particle world. The test domain is split into
+# two particle worlds, one containing particles in the upper
+# right quarter of the domain, the other covering the rest of
+# the domain. The particle weights are chosen so that
+# for a cell weight of 1000 the integrated weight from the
+# upper right quarter of the domain equals the integrated
+# weight from the other three quarters, which leads to a
+# cell distribution between the two processes of 16:48,
+# which can be seen in the output.
+
+# MPI: 2
+
+set Dimension                              = 2
+set End time                               = 0
+set Use years in output instead of seconds = false
+set Number of particle worlds = 2
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent  = 1.0000
+    set Y extent  = 1.0000
+  end
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right
+  set Zero velocity boundary indicators       = bottom, top
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density             = 1010
+    set Viscosity                     = 1e2
+    set Thermal expansion coefficient = 0
+    set Density differential for compositional field 1 = -10
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10
+  end
+end
+
+############### Parameters describing the temperature field
+# Note: The temperature plays no role in this model
+
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0
+  end
+end
+
+############### Parameters describing the compositional field
+# Note: The compositional field is what drives the flow
+# in this example
+
+subsection Compositional fields
+  set Number of fields = 1
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,z
+    set Function constants  = pi=3.1415926
+    set Function expression = 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02))
+  end
+end
+
+############### Parameters describing the discretization
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 1
+  set Strategy                           = maximum refinement function,  minimum refinement function
+  set Initial global refinement          = 2
+  set Time steps between mesh refinement = 1
+  set Coarsening fraction                = 0.05
+  set Refinement fraction                = 0.3
+  subsection Maximum refinement function
+    set Function expression = 3
+  end
+  subsection Minimum refinement function
+    set Function expression = 3
+  end
+end
+
+############### Parameters describing what to do with the solution
+
+subsection Postprocess
+  set List of postprocessors = particles, particle count statistics, load balance statistics, visualization
+
+  subsection Visualization
+    set Time between graphical output = 0
+    set Output format = gnuplot
+    set List of output variables =  partition,  particle count
+  end
+
+  subsection Particles
+    set Time between data output = 100
+    set Data output format = gnuplot
+  end
+end
+
+subsection Particles
+  set Load balancing strategy = repartition
+  set Integration scheme = euler
+  set Particle weight = 5000
+  set Particle generator name = probability density function
+
+  subsection Generator
+    subsection Probability density function
+      set Number of particles = 16
+      set Function expression = x > 0.5 ? (y > 0.5 ? 1.0 : 1e-3) : 1e-3
+      set Random cell selection = false
+    end
+  end
+end
+
+subsection Particles 2
+  set Load balancing strategy = repartition
+  set Integration scheme = euler
+  set Particle weight = 1000
+  set Particle generator name = probability density function
+
+  subsection Generator
+    subsection Probability density function
+      set Number of particles = 48
+      set Function expression = x > 0.5 ? (y > 0.5 ? 1e-3 : 1.0) : 1.0
+      set Random cell selection = false
+    end
+  end
+end

tests/particle_load_balancing_repartition_multiple_worlds/screen-output

@@ -0,0 +1,32 @@
+
+Number of active cells: 16 (on 3 levels)
+Number of degrees of freedom: 349 (162+25+81+81)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Skipping temperature solve because RHS is zero.
+   Solving C_1 system ... 0 iterations.
+   Solving Stokes system (GMG)... 11+0 iterations.
+
+Number of active cells: 64 (on 4 levels)
+Number of degrees of freedom: 1,237 (578+81+289+289)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Skipping temperature solve because RHS is zero.
+   Advecting particles...  done.
+   Advecting particles...  done.
+   Solving C_1 system ... 0 iterations.
+   Solving Stokes system (GMG)... 12+0 iterations.
+
+   Postprocessing:
+     Writing particle output:             output-particle_load_balancing_repartition_multiple_worlds/particles/particles-00000
+     Particle count per cell min/avg/max: 0, 1, 3
+     Cells per process min/max/avg:       16/48/32
+     Writing graphical output:            output-particle_load_balancing_repartition_multiple_worlds/solution/solution-00000
+
+Number of active cells: 64 (on 4 levels)
+Number of degrees of freedom: 1,237 (578+81+289+289)
+
+Termination requested by criterion: end time
+
+
+