]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Temporary working version while I fix properties
authorblaisb <blais.bruno@gmail.com>
Fri, 29 May 2020 14:33:53 +0000 (10:33 -0400)
committerBruno Blais <blais.bruno@gmail.com>
Tue, 15 Sep 2020 19:47:34 +0000 (15:47 -0400)
examples/step-68/step-68.cc

index 64071cd753874531c10d90d9e8617f525e094441..1753e0974e944e20085cfb93f65fa86752efb9b4 100644 (file)
@@ -190,13 +190,15 @@ namespace Step68
     Vortex()
       : Function<dim>(dim)
     {}
-    virtual void vector_value(const Point<dim> &point,
-                              Vector<double> &  values) const override;
+    virtual void
+    vector_value(const Point<dim> &point,
+                 Vector<double> &  values) const override;
   };
 
   template <int dim>
-  void Vortex<dim>::vector_value(const Point<dim> &point,
-                                 Vector<double> &  values) const
+  void
+  Vortex<dim>::vector_value(const Point<dim> &point,
+                            Vector<double> &  values) const
   {
     const double T = 4;
     // Since the velocity profile is time dependant, the present time in the
@@ -224,27 +226,33 @@ namespace Step68
   public:
     ParticleTracking(const ParticleTrackingParameters &par,
                      const bool                        interpolated_velocity);
-    void run();
+    void
+    run();
 
   private:
     // The particles_generation function is responsible for the initial
     // generation of the particles on top of the background grid
-    void particles_generation();
+    void
+    particles_generation();
 
     // When the velocity profile is interpolated to the position of the
     // particles, it must first be stored using degrees of freedom.
     // Consequently, as is the case for other parallel case (e.g. step-40) we
     // initialize the degrees of freedom on the background grid
-    void setup_background_dofs();
+    void
+    setup_background_dofs();
 
-    void interpolate_function_to_field();
+    void
+    interpolate_function_to_field();
 
     // The next two functions are responsible for carrying out explicit Euler
     // time integration for the cases where the velocity field is interpolated
     // at the positions of the particles or calculated analytically,
     // respectively
-    void euler_interpolated(double dt);
-    void euler_analytical(double dt);
+    void
+    euler_interpolated(double dt);
+    void
+    euler_analytical(double dt);
 
     // The cell_weight() function indicates to the triangulation how much
     // computational work is expected to happen on this cell, and consequently
@@ -252,7 +260,8 @@ namespace Step68
     // roughly equal amount of work (potentially not an equal number of cells).
     // While the function is called from the outside, it is connected to the
     // corresponding signal from inside this class, therefore it can be private.
-    unsigned int cell_weight(
+    unsigned int
+    cell_weight(
       const typename parallel::distributed::Triangulation<dim>::cell_iterator
         &cell,
       const typename parallel::distributed::Triangulation<dim>::CellStatus
@@ -261,8 +270,10 @@ namespace Step68
     // The following two functions are responsible for outputting the simulation
     // results for the particles and for the velocity profile on the background
     // mesh, respectively.
-    void output_particles(unsigned int it);
-    void output_background(unsigned int it);
+    void
+    output_particles(unsigned int it);
+    void
+    output_background(unsigned int it);
 
     // The private members of this class are similar to other parallel deal.II
     // examples. The parameters are stored as a const member. It is important
@@ -327,17 +338,19 @@ namespace Step68
   // between ranks (the connection is created inside the
   // particles_generation() function of this class).
   template <int dim>
-  unsigned int ParticleTracking<dim>::cell_weight(
+  unsigned int
+  ParticleTracking<dim>::cell_weight(
     const typename parallel::distributed::Triangulation<dim>::cell_iterator
       &                                                                  cell,
-    const typename parallel::distributed::Triangulation<dim>::CellStatus status) const
+    const typename parallel::distributed::Triangulation<dim>::CellStatus status)
+    const
   {
     // Assign no weight to cells we do not own.
     if (!cell->is_locally_owned())
       return 0;
 
     // This determines how important particle work is compared to cell
-    // work (by default every cell has a weight of 1000). 
+    // work (by default every cell has a weight of 1000).
     // We set the weight per particle much higher to indicate that
     // the particle load is the only one that is important to distribute
     // in this example. The optimal value of this number depends on the
@@ -379,7 +392,8 @@ namespace Step68
   // This function generates the tracer particles and the background
   // triangulation on which these particles evolve.
   template <int dim>
-  void ParticleTracking<dim>::particles_generation()
+  void
+  ParticleTracking<dim>::particles_generation()
   {
     // We create an hyper_cube triangulation which we globally define. This
     // triangulation englobes the full trajectory of the particles.
@@ -419,7 +433,10 @@ namespace Step68
         false));
 
     // Establish the background triangulation where the particles are living
-    particle_handler.initialize(background_triangulation, mapping);
+    // and the number of properties of the particles
+    particle_handler.initialize(background_triangulation, mapping, 1 + dim);
+    // pcout << "Number of properties "
+    //      << particle_handler.n_properties_per_particle() << std::endl;
 
     // We create a particle triangulation which is solely used to generate
     // the points which will be used to insert the particles. This
@@ -474,7 +491,8 @@ namespace Step68
   // interpolation And allocate the field vector where the entire
   // solution of the velocity field is stored
   template <int dim>
-  void ParticleTracking<dim>::setup_background_dofs()
+  void
+  ParticleTracking<dim>::setup_background_dofs()
   {
     fluid_dh.distribute_dofs(fluid_fe);
     IndexSet locally_owned_dofs = fluid_dh.locally_owned_dofs();
@@ -489,7 +507,8 @@ namespace Step68
 
   // Interpolates the Vortex velocity field to the field vector
   template <int dim>
-  void ParticleTracking<dim>::interpolate_function_to_field()
+  void
+  ParticleTracking<dim>::interpolate_function_to_field()
   {
     const MappingQ<dim> mapping(fluid_fe.degree);
 
@@ -503,7 +522,8 @@ namespace Step68
   // using an analytically defined velocity field. This is a relatively trivial
   // usage of the particles.
   template <int dim>
-  void ParticleTracking<dim>::euler_analytical(double dt)
+  void
+  ParticleTracking<dim>::euler_analytical(double dt)
   {
     Vector<double> particle_velocity(dim);
 
@@ -522,6 +542,13 @@ namespace Step68
           particle_location[d] += particle_velocity[d] * dt;
 
         particle->set_location(particle_location);
+
+        // Store the processor id and the particle velocity in the particle
+        // properties
+        ArrayView<double> properties = particle->get_properties();
+        properties[0] = Utilities::MPI::this_mpi_process(mpi_communicator);
+        for (int d = 0; d < dim; ++d)
+          properties[1 + d] += particle_velocity[d];
       }
   }
 
@@ -529,7 +556,8 @@ namespace Step68
   // We integrate the particle trajectories by interpolating the value of the
   // velocity field at the degrees of freedom to the position of the particles.
   template <int dim>
-  void ParticleTracking<dim>::euler_interpolated(double dt)
+  void
+  ParticleTracking<dim>::euler_interpolated(double dt)
   {
     std::vector<types::global_dof_index> dof_indices(fluid_fe.dofs_per_cell);
     Vector<double> dof_data_per_cell(fluid_fe.dofs_per_cell);
@@ -588,6 +616,13 @@ namespace Step68
             for (int d = 0; d < dim; ++d)
               particle_location[d] += particle_velocity[d] * dt;
             particle->set_location(particle_location);
+
+            // Store the particle velocity and the processor id in the particle
+            // properties
+            ArrayView<double> properties = particle->get_properties();
+            properties[0] = Utilities::MPI::this_mpi_process(mpi_communicator);
+            for (int d = 0; d < dim; ++d)
+              properties[1 + d] += particle_velocity[d];
           }
       }
   }
@@ -599,7 +634,8 @@ namespace Step68
   // and the background mesh to vtu with a pvtu record
 
   template <int dim>
-  void ParticleTracking<dim>::output_particles(unsigned int it)
+  void
+  ParticleTracking<dim>::output_particles(unsigned int it)
   {
     Particles::DataOut<dim, dim> particle_output;
     particle_output.build_patches(particle_handler);
@@ -615,7 +651,8 @@ namespace Step68
   }
 
   template <int dim>
-  void ParticleTracking<dim>::output_background(unsigned int it)
+  void
+  ParticleTracking<dim>::output_background(unsigned int it)
   {
     std::vector<std::string> solution_names(dim, "velocity");
     std::vector<DataComponentInterpretation::DataComponentInterpretation>
@@ -657,7 +694,8 @@ namespace Step68
   // straightforward.
 
   template <int dim>
-  void ParticleTracking<dim>::run()
+  void
+  ParticleTracking<dim>::run()
   {
     DiscreteTime discrete_time(0, par.final_time, par.time_step);
 
@@ -665,7 +703,7 @@ namespace Step68
 
     pcout << "Repartitioning triangulation after particle generation"
           << std::endl;
-    background_triangulation.repartition();
+    //    background_triangulation.repartition();
 
     setup_background_dofs();
     interpolate_function_to_field();
@@ -709,7 +747,8 @@ namespace Step68
 // The remainder of the code, the `main()` function, is standard.
 // We note that we run the particle tracking with the analytical velocity
 // and the interpolated velocity and produce both results
-int main(int argc, char *argv[])
+int
+main(int argc, char *argv[])
 {
   using namespace Step68;
   using namespace dealii;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.