From: Bruno Blais Date: Sun, 2 Jul 2023 16:03:51 +0000 (-0400) Subject: Add Step-68 performance test X-Git-Tag: relicensing~786^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=075b8bcb09cc0cdef94fd5bae684990671a6dac9;p=dealii.git Add Step-68 performance test --- diff --git a/tests/performance/timing_step_68.cc b/tests/performance/timing_step_68.cc new file mode 100644 index 0000000000..03f7b0a519 --- /dev/null +++ b/tests/performance/timing_step_68.cc @@ -0,0 +1,510 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2020 - 2023 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * This test is based on a modified version of step-68. Modifications were + * made to ensure that particles were more evenly distributed + * between the cells to ensure that the computational load could be + * balanced in some way. This is a more realistic usage of particles. + * The performance test measures + * the interpolation of a finite element field to particle location, + * the displacement of particles and their localizaton within subdomains + * and cells. + */ + + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#define ENABLE_MPI + +#include "performance_test_driver.h" + +constexpr bool debug = false; + +namespace Step68 +{ + using namespace dealii; + + + struct ParticleTrackingParameters + { + ParticleTrackingParameters() + { + if (get_testing_environment() == TestingEnvironment::medium) + { + fluid_refinement = 8; + particle_insertion_refinement = 10; + } + else if (get_testing_environment() == TestingEnvironment::heavy) + { + fluid_refinement = 9; + particle_insertion_refinement = 11; + } + } + + unsigned int velocity_degree = 1; + double time_step = 0.0001; + double final_time = 0.0050; + unsigned int output_frequency = 40; + unsigned int repartition_frequency = 5; + + unsigned int fluid_refinement = 7; + unsigned int particle_insertion_refinement = 9; + }; + + + template + class Vortex : public Function + { + public: + Vortex() + : Function(dim) + {} + + + virtual void + vector_value(const Point &point, + Vector & values) const override; + }; + + + template + void + Vortex::vector_value(const Point &point, + Vector & values) const + { + const double T = 4; + const double t = this->get_time(); + + const double px = numbers::PI * point(0); + const double py = numbers::PI * point(1); + const double pt = numbers::PI / T * t; + + values[0] = -2 * cos(pt) * pow(sin(px), 2) * sin(py) * cos(py); + values[1] = 2 * cos(pt) * pow(sin(py), 2) * sin(px) * cos(px); + if (dim == 3) + { + values[2] = 0; + } + } + + + + template + class ParticleTracking + { + public: + ParticleTracking(const ParticleTrackingParameters &par, + const bool interpolated_velocity); + + Measurement + run(); + + private: + void + generate_particles(); + + void + setup_background_dofs(); + + void + interpolate_function_to_field(); + + void + euler_step_interpolated(const double dt); + + unsigned int + cell_weight( + const typename parallel::distributed::Triangulation::cell_iterator + &cell, + const typename parallel::distributed::Triangulation::CellStatus + status) const; + + void + output_particles(); + void + output_background(); + + + const ParticleTrackingParameters ∥ + + MPI_Comm mpi_communicator; + parallel::distributed::Triangulation background_triangulation; + Particles::ParticleHandler particle_handler; + + DoFHandler fluid_dh; + FESystem fluid_fe; + MappingQ1 mapping; + LinearAlgebra::distributed::Vector velocity_field; + + Vortex velocity; + + ConditionalOStream pcout; + + bool interpolated_velocity; + }; + + + + template + ParticleTracking::ParticleTracking(const ParticleTrackingParameters &par, + const bool interpolated_velocity) + : par(par) + , mpi_communicator(MPI_COMM_WORLD) + , background_triangulation(mpi_communicator) + , fluid_dh(background_triangulation) + , fluid_fe(FE_Q(par.velocity_degree), dim) + , pcout(std::cout, + debug && Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + , interpolated_velocity(interpolated_velocity) + + {} + + + + template + unsigned int + ParticleTracking::cell_weight( + const typename parallel::distributed::Triangulation::cell_iterator + & cell, + const typename parallel::distributed::Triangulation::CellStatus status) + const + { + const unsigned int base_weight = 1; + + const unsigned int particle_weight = 1; + + unsigned int n_particles_in_cell = 0; + switch (status) + { + case parallel::distributed::Triangulation::CELL_PERSIST: + case parallel::distributed::Triangulation::CELL_REFINE: + n_particles_in_cell = particle_handler.n_particles_in_cell(cell); + break; + + case parallel::distributed::Triangulation::CELL_INVALID: + break; + + case parallel::distributed::Triangulation::CELL_COARSEN: + for (const auto &child : cell->child_iterators()) + n_particles_in_cell += particle_handler.n_particles_in_cell(child); + break; + + default: + Assert(false, ExcInternalError()); + break; + } + + return base_weight + particle_weight * n_particles_in_cell; + } + + + + template + void + ParticleTracking::generate_particles() + { + GridGenerator::hyper_cube(background_triangulation, 0, 1); + background_triangulation.refine_global(par.fluid_refinement); + + background_triangulation.signals.weight.connect( + [&]( + const typename parallel::distributed::Triangulation::cell_iterator + &cell, + const typename parallel::distributed::Triangulation::CellStatus + status) -> unsigned int { return this->cell_weight(cell, status); }); + + particle_handler.initialize(background_triangulation, mapping, 1 + dim); + + Point center; + center[0] = 0.0; + center[1] = 0.0; + if (dim == 3) + center[2] = 0.0; + + const double outer_radius = 0.50; + const double inner_radius = 0.01; + + parallel::distributed::Triangulation particle_triangulation( + MPI_COMM_WORLD); + + GridGenerator::hyper_cube(particle_triangulation, 0, 1); + particle_triangulation.refine_global(par.particle_insertion_refinement); + + const auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box( + background_triangulation, IteratorFilters::LocallyOwnedCell()); + const auto global_bounding_boxes = + Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box); + + std::vector> properties( + particle_triangulation.n_locally_owned_active_cells(), + std::vector(dim + 1, 0.)); + + Particles::Generators::quadrature_points(particle_triangulation, + QMidpoint(), + global_bounding_boxes, + particle_handler, + mapping, + properties); + + pcout << "Number of particles inserted: " + << particle_handler.n_global_particles() << std::endl; + } + + + + template + void + ParticleTracking::setup_background_dofs() + { + fluid_dh.distribute_dofs(fluid_fe); + const IndexSet locally_owned_dofs = fluid_dh.locally_owned_dofs(); + const IndexSet locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(fluid_dh); + + velocity_field.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + } + + + + template + void + ParticleTracking::interpolate_function_to_field() + { + velocity_field.zero_out_ghost_values(); + VectorTools::interpolate(mapping, fluid_dh, velocity, velocity_field); + velocity_field.update_ghost_values(); + } + + + template + void + ParticleTracking::euler_step_interpolated(const double dt) + { + Vector local_dof_values(fluid_fe.dofs_per_cell); + + auto particle = particle_handler.begin(); + while (particle != particle_handler.end()) + { + const auto cell = particle->get_surrounding_cell(); + const auto dh_cell = + typename DoFHandler::cell_iterator(*cell, &fluid_dh); + + dh_cell->get_dof_values(velocity_field, local_dof_values); + + const auto pic = particle_handler.particles_in_cell(cell); + Assert(pic.begin() == particle, ExcInternalError()); + for (auto &p : pic) + { + const Point reference_location = p.get_reference_location(); + Tensor<1, dim> particle_velocity; + for (unsigned int j = 0; j < fluid_fe.dofs_per_cell; ++j) + { + const auto comp_j = fluid_fe.system_to_component_index(j); + + particle_velocity[comp_j.first] += + fluid_fe.shape_value(j, reference_location) * + local_dof_values[j]; + } + + Point particle_location = particle->get_location(); + for (int d = 0; d < dim; ++d) + particle_location[d] += particle_velocity[d] * dt; + p.set_location(particle_location); + + ArrayView properties = p.get_properties(); + for (int d = 0; d < dim; ++d) + properties[d] = particle_velocity[d]; + + properties[dim] = + Utilities::MPI::this_mpi_process(mpi_communicator); + + ++particle; + } + } + } + + + + template + void + ParticleTracking::output_particles() + { + Particles::DataOut particle_output; + + std::vector solution_names(dim, "velocity"); + solution_names.emplace_back("process_id"); + + std::vector + data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + + particle_output.build_patches(particle_handler, + solution_names, + data_component_interpretation); + } + + + + template + void + ParticleTracking::output_background() + { + std::vector solution_names(dim, "velocity"); + std::vector + data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + + DataOut data_out; + + data_out.attach_dof_handler(fluid_dh); + data_out.add_data_vector(velocity_field, + solution_names, + DataOut::type_dof_data, + data_component_interpretation); + Vector subdomain(background_triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = background_triangulation.locally_owned_subdomain(); + data_out.add_data_vector(subdomain, "subdomain"); + + data_out.build_patches(mapping); + } + + + + template + Measurement + ParticleTracking::run() + { + std::map timer; + + DiscreteTime discrete_time(0, par.final_time, par.time_step); + + timer["generate_particles"].start(); + generate_particles(); + timer["generate_particles"].stop(); + + pcout << "Repartitioning triangulation after particle generation" + << std::endl; + + particle_handler.prepare_for_coarsening_and_refinement(); + background_triangulation.repartition(); + particle_handler.unpack_after_coarsening_and_refinement(); + + setup_background_dofs(); + interpolate_function_to_field(); + euler_step_interpolated(0.); + + + while (!discrete_time.is_at_end()) + { + discrete_time.advance_time(); + velocity.set_time(discrete_time.get_previous_time()); + + if ((discrete_time.get_step_number() % par.repartition_frequency) == 0) + { + timer["load_balance"].start(); + particle_handler.prepare_for_coarsening_and_refinement(); + background_triangulation.repartition(); + particle_handler.unpack_after_coarsening_and_refinement(); + + setup_background_dofs(); + timer["load_balance"].stop(); + } + + + timer["advect"].start(); + interpolate_function_to_field(); + euler_step_interpolated(discrete_time.get_previous_step_size()); + timer["advect"].stop(); + + + timer["sort"].start(); + particle_handler.sort_particles_into_subdomains_and_cells(); + timer["sort"].stop(); + + if ((discrete_time.get_step_number() % par.output_frequency) == 0) + { + timer["output"].start(); + output_particles(); + output_background(); + timer["output"].stop(); + } + } + + return {timer["generate_particles"].wall_time(), + timer["load_balance"].wall_time(), + timer["advect"].wall_time(), + timer["sort"].wall_time(), + timer["output"].wall_time()}; + } + +} // namespace Step68 + + +std::tuple> +describe_measurements() +{ + return {Metric::timing, + 4, + {"generate_particles", "load_balance", "advect", "sort", "output"}}; +} + + +Measurement +perform_single_measurement() +{ + using namespace Step68; + using namespace dealii; + + ParticleTrackingParameters par; + + Step68::ParticleTracking<2> particle_tracking(par, true); + return particle_tracking.run(); +} diff --git a/tests/performance/timing_step_68.threads=1.mpirun=max.exclusive.release.run_only b/tests/performance/timing_step_68.threads=1.mpirun=max.exclusive.release.run_only new file mode 100644 index 0000000000..e69de29bb2