From 5779d67c0fb940591458396971f0946ae2e4cf98 Mon Sep 17 00:00:00 2001 From: blaisb Date: Tue, 2 Jun 2020 13:39:08 -0400 Subject: [PATCH] - Removed dependency on trilinos or petsc to use deal.II distributed vectors - Finalized step-68 --- examples/step-68/CMakeLists.txt | 13 +--- examples/step-68/step-68.cc | 113 +++++++++++--------------------- 2 files changed, 39 insertions(+), 87 deletions(-) diff --git a/examples/step-68/CMakeLists.txt b/examples/step-68/CMakeLists.txt index fc66b05510..d704f7d531 100644 --- a/examples/step-68/CMakeLists.txt +++ b/examples/step-68/CMakeLists.txt @@ -37,22 +37,13 @@ ENDIF() # # Are all dependencies fulfilled? # -IF(NOT ((DEAL_II_WITH_PETSC AND NOT DEAL_II_PETSC_WITH_COMPLEX) OR DEAL_II_WITH_TRILINOS) OR NOT DEAL_II_WITH_P4EST) # keep in one line +IF(NOT DEAL_II_WITH_P4EST) # keep in one line MESSAGE(FATAL_ERROR " Error! This tutorial requires a deal.II library that was configured with the following options: - DEAL_II_WITH_PETSC = ON - DEAL_II_PETSC_WITH_COMPLEX = OFF - DEAL_II_WITH_P4EST = ON -or - DEAL_II_WITH_TRILINOS = ON DEAL_II_WITH_P4EST = ON However, the deal.II library found at ${DEAL_II_PATH} was configured with these options - DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC} - DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX} DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST} - DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS} -which conflict with the requirements. -One or both of the aforementioned combinations of prerequisites are not met by your installation, but at least one is required for this tutorial step." +which conflict with the requirements." ) ENDIF() diff --git a/examples/step-68/step-68.cc b/examples/step-68/step-68.cc index 73aa1a0af4..921d955914 100644 --- a/examples/step-68/step-68.cc +++ b/examples/step-68/step-68.cc @@ -44,9 +44,7 @@ #include #include -#include -#include -#include +#include #include #include @@ -71,24 +69,6 @@ // to commonly used parallel vtu format #include -// This step uses parallel vector to interpolate the velocity field -// at the position of the particles. This step supports the use of both -// Trilinos and PETSC distributed vectors -#define FORCE_USE_OF_TRILINOS - -namespace LA -{ -#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \ - !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS)) - using namespace dealii::LinearAlgebraPETSc; -# define USE_PETSC_LA -#elif defined(DEAL_II_WITH_TRILINOS) - using namespace dealii::LinearAlgebraTrilinos; -#else -# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required -#endif -} // namespace LA - #include #include @@ -185,15 +165,13 @@ 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 @@ -221,33 +199,27 @@ 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 @@ -255,8 +227,7 @@ 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 @@ -265,10 +236,8 @@ 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 @@ -281,11 +250,11 @@ namespace Step68 parallel::distributed::Triangulation background_triangulation; Particles::ParticleHandler particle_handler; - DoFHandler fluid_dh; - FESystem fluid_fe; - MappingQ mapping; - LA::MPI::Vector field_owned; - LA::MPI::Vector field_relevant; + DoFHandler fluid_dh; + FESystem fluid_fe; + MappingQ mapping; + LinearAlgebra::distributed::Vector field_owned; + LinearAlgebra::distributed::Vector field_relevant; Vortex velocity; @@ -299,7 +268,7 @@ namespace Step68 // @sect4{Constructor} // Constructors and destructors are rather trivial. They are very similar - // to what is done in step-40. we set the set of processors we want to work on + // to what is done in step-40. We set the processors we want to work on // to all machines available (MPI_COMM_WORLD) and // initialize the pcout variable to only allow processor zero // to output anything to the standard output. @@ -333,8 +302,7 @@ 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) @@ -387,8 +355,7 @@ 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. @@ -461,6 +428,10 @@ namespace Step68 const auto global_bounding_boxes = Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box); + // The quadrature points particle generator generates particles only + // on locally owned active cells. We therefore count how many of those + // are present in the triangulation, as this will be required to + // initialize the properties unsigned int n_locally_owned_cells = 0; for (const auto &cell : particle_triangulation.active_cell_iterators()) if (cell->is_locally_owned()) @@ -492,8 +463,7 @@ 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(); @@ -508,8 +478,7 @@ 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); @@ -523,8 +492,7 @@ 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); @@ -556,8 +524,7 @@ 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); @@ -635,8 +602,7 @@ 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; @@ -651,11 +617,9 @@ namespace Step68 data_component_interpretation.push_back( DataComponentInterpretation::component_is_scalar); - particle_output.build_patches(particle_handler); - - // particle_output.build_patches(particle_handler, - // solution_names, - // data_component_interpretation); + particle_output.build_patches(particle_handler, + solution_names, + data_component_interpretation); std::string output_folder(par.output_directory); std::string file_name(interpolated_velocity ? "interpolated-particles" : "analytical-particles"); @@ -668,8 +632,7 @@ 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 @@ -711,8 +674,7 @@ namespace Step68 // straightforward. template - void - ParticleTracking::run() + void ParticleTracking::run() { DiscreteTime discrete_time(0, par.final_time, par.time_step); @@ -774,8 +736,7 @@ 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; -- 2.39.5