From: blaisb Date: Fri, 29 May 2020 14:33:53 +0000 (-0400) Subject: Temporary working version while I fix properties X-Git-Tag: v9.3.0-rc1~1101^2~15 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5dd44fa4855132b99065b47f0f5780719c628706;p=dealii.git Temporary working version while I fix properties --- diff --git a/examples/step-68/step-68.cc b/examples/step-68/step-68.cc index 64071cd753..1753e0974e 100644 --- a/examples/step-68/step-68.cc +++ b/examples/step-68/step-68.cc @@ -190,13 +190,15 @@ namespace Step68 Vortex() : Function(dim) {} - virtual void vector_value(const Point &point, - Vector & values) const override; + virtual void + vector_value(const Point &point, + Vector & values) const override; }; template - void Vortex::vector_value(const Point &point, - Vector & values) const + void + Vortex::vector_value(const Point &point, + Vector & 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::cell_iterator &cell, const typename parallel::distributed::Triangulation::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 - unsigned int ParticleTracking::cell_weight( + unsigned int + ParticleTracking::cell_weight( const typename parallel::distributed::Triangulation::cell_iterator & cell, - const typename parallel::distributed::Triangulation::CellStatus status) const + const typename parallel::distributed::Triangulation::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 - void ParticleTracking::particles_generation() + void + ParticleTracking::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 - void ParticleTracking::setup_background_dofs() + void + ParticleTracking::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 - void ParticleTracking::interpolate_function_to_field() + void + ParticleTracking::interpolate_function_to_field() { const MappingQ 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 - void ParticleTracking::euler_analytical(double dt) + void + ParticleTracking::euler_analytical(double dt) { Vector 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 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 - void ParticleTracking::euler_interpolated(double dt) + void + ParticleTracking::euler_interpolated(double dt) { std::vector dof_indices(fluid_fe.dofs_per_cell); Vector 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 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 - void ParticleTracking::output_particles(unsigned int it) + void + ParticleTracking::output_particles(unsigned int it) { Particles::DataOut particle_output; particle_output.build_patches(particle_handler); @@ -615,7 +651,8 @@ namespace Step68 } template - void ParticleTracking::output_background(unsigned int it) + void + ParticleTracking::output_background(unsigned int it) { std::vector solution_names(dim, "velocity"); std::vector @@ -657,7 +694,8 @@ namespace Step68 // straightforward. template - void ParticleTracking::run() + void + ParticleTracking::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;