Merge pull request #5874 from gassmoeller/cell_wise_particle_update

Cell wise grain size particle update

include/aspect/particle/property/grain_size.h

@@ -61,13 +61,13 @@ namespace aspect
                                             std::vector<double> &particle_properties) const override;
 
           /**
-           * @copydoc aspect::Particle::Property::Interface::update_particle_property()
+           * @copydoc aspect::Particle::Property::Interface::update_particle_properties()
            */
           void
-          update_particle_property (const unsigned int data_position,
-                                    const Vector<double> &solution,
-                                    const std::vector<Tensor<1,dim>> &gradients,
-                                    typename ParticleHandler<dim>::particle_iterator &particle) const override;
+          update_particle_properties (const unsigned int data_position,
+                                      const std::vector<Vector<double>> &solution,
+                                      const std::vector<std::vector<Tensor<1,dim>>> &gradients,
+                                      typename ParticleHandler<dim>::particle_iterator_range &particles) const override;
 
           /**
            * Returns an enum, which determines how this particle property is

include/aspect/particle/property/interface.h

@@ -324,6 +324,39 @@ namespace aspect
           initialize_one_particle_property (const Point<dim> &position,
                                             std::vector<double> &particle_properties) const;
 
+          /**
+           * Update function. This function is called every time an update is
+           * requested by need_update() for every cell for every property.
+           * It is expected to update the properties of all particles in the
+           * given range @p particles, which are all in one cell.
+           * It is obvious that
+           * this function is called a lot, so its code should be efficient.
+           * The interface provides a default implementation that does nothing,
+           * therefore derived plugins that do not require an update do not
+           * need to implement this function.
+           *
+           * @param [in] data_position An unsigned integer that denotes which
+           * component of each particle property vector is associated with the
+           * current property. For properties that own several components it
+           * denotes the first component of this property, all other components
+           * fill consecutive entries in the properties vector.
+           *
+           * @param [in] solution A vector of values of the solution variables
+           * at the given particle positions.
+           *
+           * @param [in] gradients A vector of gradients of the solution
+           * variables at the given particle positions.
+           *
+           * @param [in,out] particles The particles that are to be updated
+           * within this function.
+           */
+          virtual
+          void
+          update_particle_properties (const unsigned int data_position,
+                                      const std::vector<Vector<double>> &solution,
+                                      const std::vector<std::vector<Tensor<1,dim>>> &gradients,
+                                      typename ParticleHandler<dim>::particle_iterator_range &particles) const;
+
           /**
            * Update function. This function is called every time an update is
            * request by need_update() for every particle for every property.
@@ -349,7 +382,10 @@ namespace aspect
            * the call of this function. The particle location can be accessed
            * using particle->get_location() and its properties using
            * particle->get_properties().
+           *
+           * @deprecated Use update_particle_properties() instead.
            */
+          DEAL_II_DEPRECATED
           virtual
           void
           update_particle_property (const unsigned int data_position,
@@ -526,12 +562,19 @@ namespace aspect
 
           /**
            * Update function for particle properties. This function is
-           * called once every time step for every particle.
+           * called once every time step for every cell.
+           *
+           * @param particles The particles that are to be updated within
+           * this function.
+           * @param solution The values of the solution variables at the
+           * given particle positions.
+           * @param gradients The gradients of the solution variables at
+           * the given particle positions.
            */
           void
-          update_one_particle (typename ParticleHandler<dim>::particle_iterator &particle,
-                               const Vector<double> &solution,
-                               const std::vector<Tensor<1,dim>> &gradients) const;
+          update_particles (typename ParticleHandler<dim>::particle_iterator_range &particles,
+                            const std::vector<Vector<double>> &solution,
+                            const std::vector<std::vector<Tensor<1,dim>>> &gradients) const;
 
           /**
            * Returns an enum, which denotes at what time this class needs to

include/aspect/particle/world.h

@@ -484,8 +484,6 @@ namespace aspect
          */
         void
         local_update_particles(const typename DoFHandler<dim>::active_cell_iterator &cell,
-                               const typename ParticleHandler<dim>::particle_iterator &begin_particle,
-                               const typename ParticleHandler<dim>::particle_iterator &end_particle,
                                internal::SolutionEvaluators<dim> &evaluators);
 
         /**

source/particle/property/grain_size.cc

@@ -66,39 +66,53 @@ namespace aspect
 
       template <int dim>
       void
-      GrainSize<dim>::update_particle_property(const unsigned int data_position,
-                                               const Vector<double> &solution,
-                                               const std::vector<Tensor<1,dim>> &gradients,
-                                               typename ParticleHandler<dim>::particle_iterator &particle) const
+      GrainSize<dim>::update_particle_properties(const unsigned int data_position,
+                                                 const std::vector<Vector<double>> &solution,
+                                                 const std::vector<std::vector<Tensor<1,dim>>> &gradients,
+                                                 typename ParticleHandler<dim>::particle_iterator_range &particles) const
       {
-        material_inputs.position[0] = particle->get_location();
-
-        material_inputs.current_cell = typename DoFHandler<dim>::active_cell_iterator(*particle->get_surrounding_cell(),
+        material_inputs  = MaterialModel::MaterialModelInputs<dim>(solution.size(), this->n_compositional_fields());
+        material_outputs = MaterialModel::MaterialModelOutputs<dim>(solution.size(), this->n_compositional_fields());
+        material_inputs.requested_properties = MaterialModel::MaterialProperties::reaction_terms;
+        material_inputs.current_cell = typename DoFHandler<dim>::active_cell_iterator(*particles.begin()->get_surrounding_cell(),
                                                                                       &(this->get_dof_handler()));
 
-        material_inputs.temperature[0] = solution[this->introspection().component_indices.temperature];
+        unsigned int i = 0;
+        for (auto particle: particles)
+          {
+            // Make sure all particles are in the same cell
+            Assert(particle.get_surrounding_cell() == particles.begin()->get_surrounding_cell(),
+                   ExcMessage("All particles must be in the same cell."));
 
-        material_inputs.pressure[0] = solution[this->introspection().component_indices.pressure];
+            material_inputs.position[i] = particle.get_location();
+            material_inputs.temperature[i] = solution[i][this->introspection().component_indices.temperature];
+            material_inputs.pressure[i] = solution[i][this->introspection().component_indices.pressure];
 
-        for (unsigned int d = 0; d < dim; ++d)
-          material_inputs.velocity[0][d] = solution[this->introspection().component_indices.velocities[d]];
+            for (unsigned int d = 0; d < dim; ++d)
+              material_inputs.velocity[i][d] = solution[i][this->introspection().component_indices.velocities[d]];
 
-        for (unsigned int n = 0; n < this->n_compositional_fields(); ++n)
-          material_inputs.composition[0][n] = solution[this->introspection().component_indices.compositional_fields[n]];
+            for (unsigned int n = 0; n < this->n_compositional_fields(); ++n)
+              material_inputs.composition[i][n] = solution[i][this->introspection().component_indices.compositional_fields[n]];
 
-        material_inputs.composition[0][grain_size_index] = particle->get_properties()[data_position];
+            material_inputs.composition[i][grain_size_index] = particle.get_properties()[data_position];
 
-        Tensor<2,dim> grad_u;
-        for (unsigned int d=0; d<dim; ++d)
-          grad_u[d] = gradients[d];
-        material_inputs.strain_rate[0] = symmetrize (grad_u);
+            Tensor<2,dim> grad_u;
+            for (unsigned int d=0; d<dim; ++d)
+              grad_u[d] = gradients[i][d];
+            material_inputs.strain_rate[i] = symmetrize (grad_u);
 
-        material_inputs.requested_properties = MaterialModel::MaterialProperties::reaction_terms;
+            ++i;
+          }
 
         this->get_material_model().evaluate(material_inputs,
                                             material_outputs);
 
-        particle->get_properties()[data_position] += material_outputs.reaction_terms[0][grain_size_index];
+        i = 0;
+        for (auto particle: particles)
+          {
+            particle.get_properties()[data_position] += material_outputs.reaction_terms[i][grain_size_index];
+            ++i;
+          }
       }
 
 

source/particle/property/interface.cc

@@ -206,6 +206,39 @@ namespace aspect
 
 
 
+      DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+      template <int dim>
+      void
+      Interface<dim>::update_particle_properties (const unsigned int data_position,
+                                                  const std::vector<Vector<double>> &solution,
+                                                  const std::vector<std::vector<Tensor<1,dim>>> &gradients,
+                                                  typename ParticleHandler<dim>::particle_iterator_range &particles) const
+      {
+        unsigned int i = 0;
+        for (typename ParticleHandler<dim>::particle_iterator particle = particles.begin();
+             particle != particles.end(); ++particle, ++i)
+          {
+            // call the deprecated version of this function
+            update_particle_property(data_position,
+                                     solution[i],
+                                     gradients[i],
+                                     particle);
+          }
+      }
+      DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
+
+
+
+      template <int dim>
+      void
+      Interface<dim>::update_particle_property (const unsigned int /*data_position*/,
+                                                const Vector<double> &/*solution*/,
+                                                const std::vector<Tensor<1,dim>> &/*gradients*/,
+                                                typename ParticleHandler<dim>::particle_iterator &/*particle*/) const
+      {}
+
+
+
       template <int dim>
       UpdateTimeFlags
       Interface<dim>::need_update () const
@@ -243,16 +276,6 @@ namespace aspect
 
 
 
-      template <int dim>
-      void
-      Interface<dim>::update_particle_property (const unsigned int /*data_position*/,
-                                                const Vector<double> &/*solution*/,
-                                                const std::vector<Tensor<1,dim>> &/*gradients*/,
-                                                typename ParticleHandler<dim>::particle_iterator &/*particle*/) const
-      {}
-
-
-
       template <int dim>
       std::vector<std::pair<std::string, unsigned int>>
       IntegratorProperties<dim>::get_property_information() const
@@ -518,18 +541,18 @@ namespace aspect
 
       template <int dim>
       void
-      Manager<dim>::update_one_particle (typename ParticleHandler<dim>::particle_iterator &particle,
-                                         const Vector<double> &solution,
-                                         const std::vector<Tensor<1,dim>> &gradients) const
+      Manager<dim>::update_particles (typename ParticleHandler<dim>::particle_iterator_range &particles,
+                                      const std::vector<Vector<double>> &solution,
+                                      const std::vector<std::vector<Tensor<1,dim>>> &gradients) const
       {
         unsigned int plugin_index = 0;
         for (typename std::list<std::unique_ptr<Interface<dim>>>::const_iterator
              p = property_list.begin(); p!=property_list.end(); ++p,++plugin_index)
           {
-            (*p)->update_particle_property(property_information.get_position_by_plugin_index(plugin_index),
-                                           solution,
-                                           gradients,
-                                           particle);
+            (*p)->update_particle_properties(property_information.get_position_by_plugin_index(plugin_index),
+                                             solution,
+                                             gradients,
+                                             particles);
           }
       }
 

source/particle/world.cc

@@ -444,17 +444,16 @@ namespace aspect
     template <int dim>
     void
     World<dim>::local_update_particles(const typename DoFHandler<dim>::active_cell_iterator &cell,
-                                       const typename ParticleHandler<dim>::particle_iterator &begin_particle,
-                                       const typename ParticleHandler<dim>::particle_iterator &end_particle,
                                        internal::SolutionEvaluators<dim> &evaluators)
     {
       const unsigned int n_particles_in_cell = particle_handler->n_particles_in_cell(cell);
+      typename ParticleHandler<dim>::particle_iterator_range particles = particle_handler->particles_in_cell(cell);
 
       std::vector<Point<dim>> positions;
       positions.reserve(n_particles_in_cell);
 
-      for (auto particle = begin_particle; particle!=end_particle; ++particle)
-        positions.push_back(particle->get_reference_location());
+      for (const auto &particle : particles)
+        positions.push_back(particle.get_reference_location());
 
       const UpdateFlags update_flags = property_manager->get_needed_update_flags();
 
@@ -467,29 +466,28 @@ namespace aspect
       if (update_flags & (update_values | update_gradients))
         evaluators.reinit(cell, positions, {solution_values.data(), solution_values.size()}, update_flags);
 
-      Vector<double> solution;
+      std::vector<Vector<double>> solution;
       if (update_flags & update_values)
-        solution.reinit(this->introspection().n_components);
+        solution.resize(n_particles_in_cell,Vector<double>(this->introspection().n_components));
 
-      std::vector<Tensor<1,dim>> gradients;
+      std::vector<std::vector<Tensor<1,dim>>> gradients;
       if (update_flags & update_gradients)
-        gradients.resize(this->introspection().n_components);
+        gradients.resize(n_particles_in_cell,std::vector<Tensor<1,dim>>(this->introspection().n_components));
 
-      auto particle = begin_particle;
-      for (unsigned int i = 0; particle!=end_particle; ++particle,++i)
+      for (unsigned int i = 0; i<n_particles_in_cell; ++i)
         {
           // Evaluate the solution, but only if it is requested in the update_flags
           if (update_flags & update_values)
-            evaluators.get_solution(i, solution);
+            evaluators.get_solution(i, solution[i]);
 
           // Evaluate the gradients, but only if they are requested in the update_flags
           if (update_flags & update_gradients)
-            evaluators.get_gradients(i, gradients);
-
-          property_manager->update_one_particle(particle,
-                                                solution,
-                                                gradients);
+            evaluators.get_gradients(i, gradients[i]);
         }
+
+      property_manager->update_particles(particles,
+                                         solution,
+                                         gradients);
     }
 
 
@@ -1067,15 +1065,10 @@ namespace aspect
           for (const auto &cell : this->get_dof_handler().active_cell_iterators())
             if (cell->is_locally_owned())
               {
-                typename ParticleHandler<dim>::particle_iterator_range
-                particles_in_cell = particle_handler->particles_in_cell(cell);
-
-                // Only update particles, if there are any in this cell
-                if (particles_in_cell.begin() != particles_in_cell.end())
+                // Only update particles if there are any in this cell
+                if (particle_handler->n_particles_in_cell(cell) > 0)
                   {
                     local_update_particles(cell,
-                                           particles_in_cell.begin(),
-                                           particles_in_cell.end(),
                                            *evaluators);
                   }
 

tests/grain_size_phase_function/screen-output

@@ -16,7 +16,7 @@ Number of degrees of freedom: 4,645 (2,178+289+1,089+1,089)
 *** Timestep 1:  t=10000 years, dt=10000 years
    Solving temperature system... 0 iterations.
    Advecting particles...  done.
-   Solving Stokes system... 141+0 iterations.
+   Solving Stokes system... 136+0 iterations.
 
    Postprocessing:
      Writing graphical output:                          output-grain_size_phase_function/solution/solution-00001

tests/grain_size_phase_function/statistics

@@ -20,4 +20,4 @@
 # 20: Total mass (kg)
 # 21: Average iterations for ODE solver
 0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 1089 0 138 140 140 output-grain_size_phase_function/solution/solution-00000 1.00000000e-03 1.00000000e-03 1.00000000e+07 10000 output-grain_size_phase_function/particles/particles-00000 1.00000000e+03 3.19873350e+19 1.00000000e+13 nan 
-1 1.000000000000e+04 1.000000000000e+04 256 2467 1089 1089 0 140 142 142 output-grain_size_phase_function/solution/solution-00001 1.00000000e-03 1.60801892e-03 1.30379556e+07 10000 output-grain_size_phase_function/particles/particles-00001 1.00000000e+03 6.97386958e+19 1.00000000e+13 nan 
+1 1.000000000000e+04 1.000000000000e+04 256 2467 1089 1089 0 135 137 137 output-grain_size_phase_function/solution/solution-00001 1.00000000e-03 1.60801892e-03 1.30380370e+07 10000 output-grain_size_phase_function/particles/particles-00001 1.00000000e+03 6.97407271e+19 1.00000000e+13 nan 

tests/grain_size_ridge/screen-output

@@ -69,18 +69,18 @@ Number of degrees of freedom: 57,556 (21,186+2,737+10,593+23,040)
    Solving temperature system... 7 iterations.
    Advecting particles...  done.
    Solving Stokes system... 1+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.07544e-05, 0, 0.000124181
-      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.000124181
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.07544e-05, 0, 0.00012418
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.00012418
 
    Solving temperature system... 5 iterations.
    Advecting particles...  done.
    Solving Stokes system... 1+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.44492e-07, 0, 7.78241e-05
-      Relative nonlinear residual (total system) after nonlinear iteration 2: 7.78241e-05
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.44489e-07, 0, 7.78237e-05
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 7.78237e-05
 
 
    Postprocessing:
-     Compositions min/max/mass:          0.0006565/0.005/8.175e+09
+     Compositions min/max/mass:          0.0006567/0.005/8.175e+09
      RMS, max velocity:                  0.0355 m/year, 0.118 m/year
      Mass fluxes through boundary parts: 0 kg/yr, 4.133e+07 kg/yr, -1.688e+07 kg/yr, 0 kg/yr
      Number of advected particles:       499998