From 9be492777b41d7e443db8bbaca76c8ea60321649 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 28 Apr 2020 23:36:19 +0200 Subject: [PATCH] Some documentation. --- examples/step-70/step-70.cc | 701 ++++++++++++++++++++++++------------ 1 file changed, 471 insertions(+), 230 deletions(-) diff --git a/examples/step-70/step-70.cc b/examples/step-70/step-70.cc index b1507ebef4..1fdf21aa8f 100644 --- a/examples/step-70/step-70.cc +++ b/examples/step-70/step-70.cc @@ -17,6 +17,10 @@ * Authors: Luca Heltai, Bruno Blais, 2019 */ +// @sect3{Include files} +// Most of these have been introduced elsewhere, we'll comment only on the new +// ones. + #include #include #include @@ -26,8 +30,6 @@ #include #include -#include - #define FORCE_USE_OF_TRILINOS namespace LA @@ -88,6 +90,39 @@ namespace LA #include #include +// These are the only new include files w.r.t. step-60. In this tutorial, +// the non-matching coupling between the solid and the fluid is computed using +// an intermediate data structure that keeps track of how the quadrature points +// of the solid evolve w.r.t. the fluid mesh. This data structure needs to keep +// track of the position of the quadrature points on each cell describing the +// solid domain, of the quadrature weights, and possibly of the normal vector +// to each point, if the solid domain is of co-dimension one. +// +// Deal.II offers these facilities on the Particles namespace, through the +// ParticleHandler class. ParticleHandler is a class that allows you to manage +// a collection of particles (objects of type Particles::Particle), representing +// a collection of points with some attached properties floating on a +// parallel::distributed::Triangulation. The methods and classes on the +// namespace Particles allows one to easily implement Particle In Cell methods +// and particle tracing on distributed triangulations. +// +// We "abuse" this data structure to store information about the location of +// solid quadrature points w.r.t. to the surrounding fluid grid, including +// integration weights, and possibly surface normals. The reason why we use this +// additional data structure is related to the fact that the solid and the fluid +// grids are non-overlapping, and distributed independently among processes. +// +// In order to couple the two problems, we rely on the ParticleHandler class, +// storing in each particle the position of a solid quadrature point (which is +// in general not aligned to any of the fluid quadrature points), its weight, +// and any other information that may be required to couple the two problems. +// +// Ownership of the solid quadrature points is inherited by the MPI partitioning +// on the solid mesh itslef. The Particles so generated are later distributed to +// the fluid mesh using the methods of the ParticleHandler class. This allows +// transparent exchange of information between MPI processes about the +// overlapping pattern between fluid cells and solid quadrature points. +#include #include #include @@ -100,6 +135,7 @@ namespace Step70 { using namespace dealii; + // REMOVE THIS FUNCTION ONCE #9891 is merged. template class StokesImmersedProblemParameters : public ParameterAcceptor { public: - StokesImmersedProblemParameters() - : ParameterAcceptor("Stokes Immersed Problem/") - , rhs("Right hand side", spacedim + 1) - , angular_velocity("Angular velocity", spacedim == 3 ? spacedim : 1) - { - add_parameter("Velocity degree", - velocity_degree, - "", - this->prm, - Patterns::Integer(1)); - - add_parameter("Number of time steps", number_of_time_steps); - add_parameter("Output frequency", mod_output); - - add_parameter("Final time", final_time); - - add_parameter("Viscosity", viscosity); - - add_parameter("Nitsche penalty term", penalty_term); - - add_parameter("Initial fluid refinement", - initial_fluid_refinement, - "Initial mesh refinement used for the fluid domain Omega"); - - add_parameter("Initial solid refinement", - initial_solid_refinement, - "Initial mesh refinement used for the solid domain Gamma"); - - add_parameter( - "Particle insertion refinement", - particle_insertion_refinement, - "Refinement of the volumetric mesh used to insert the particles"); - - add_parameter( - "Homogeneous Dirichlet boundary ids", - homogeneous_dirichlet_ids, - "Boundary Ids over which homogeneous Dirichlet boundary conditions are applied"); - - enter_my_subsection(this->prm); - this->prm.enter_subsection("Grid generation"); - this->prm.add_parameter("Grid one generator", name_of_grid1); - this->prm.add_parameter("Grid one generator arguments", - arguments_for_grid1); - - this->prm.add_parameter("Grid two generator", name_of_grid2); - this->prm.add_parameter("Grid two generator arguments", - arguments_for_grid2); - - this->prm.add_parameter("Particle grid generator", name_of_particle_grid); - this->prm.add_parameter("Particle grid generator arguments", - arguments_for_particle_grid); - this->prm.leave_subsection(); - - leave_my_subsection(this->prm); - - - - enter_my_subsection(this->prm); - this->prm.enter_subsection("Refinement and remeshing"); - this->prm.add_parameter("Refinement step frequency", mod_refinement); - this->prm.add_parameter("Refinement maximal level", max_level_refinement); - this->prm.add_parameter("Refinement strategy", - refinement_strategy, - "", - Patterns::Selection( - "fixed_fraction|fixed_number")); - this->prm.add_parameter("Refinement coarsening fraction", - coarsening_fraction); - this->prm.add_parameter("Refinement fraction", refinement_fraction); - this->prm.add_parameter("Maximum number of cells", max_cells); - - this->prm.leave_subsection(); - leave_my_subsection(this->prm); - - // correct the default dimension for the functions - rhs.declare_parameters_call_back.connect([&]() { - Functions::ParsedFunction::declare_parameters(this->prm, - spacedim + 1); - }); - angular_velocity.declare_parameters_call_back.connect([&]() { - Functions::ParsedFunction::declare_parameters( - this->prm, spacedim == 3 ? spacedim : 1); - }); - } - + // The constructor is responsible for the connection between the members of + // this class and the corresponding entries in the ParameterHandler. Thanks + // to the use of the ParameterHandler::add_parameter() method, this + // connection is trivial, but requires all members of this class to be + // writeable + StokesImmersedProblemParameters(); + + // however, since this class will be passed as a const reference to the + // StokesImmersedProblem class, we have to make sure we can still set the + // time correctly in the objects derived by the Function class defined + // here. In order to do so, we declare both the + // StokesImmersedProblemParameters::rhs and + // StokesImmersedProblemParameters::angular_velocity members to be mutable, + // and define this little helper method that sets their time to the correct + // value. void set_time(const double &time) const { rhs.set_time(time); angular_velocity.set_time(time); } - unsigned int velocity_degree = 2; - unsigned int number_of_time_steps = 1; - double viscosity = 1.0; - double final_time = 1.0; - unsigned int initial_fluid_refinement = 3; - unsigned int initial_solid_refinement = 3; - unsigned int particle_insertion_refinement = 1; - double penalty_term = 1e3; - std::list homogeneous_dirichlet_ids{0, 1, 2, 3}; - std::string name_of_grid1 = "hyper_cube"; - std::string arguments_for_grid1 = "-1: 1: false"; - std::string name_of_grid2 = "hyper_rectangle"; - std::string arguments_for_grid2 = + // We will use a Taylor-Hood function space of arbitrary order. This + // parameter is used to initialize the FiniteElement space with the corret + // FESystem object + unsigned int velocity_degree = 2; + + // Instead of defining a time step increment, in this tutorial we prefer to + // let the user choose a final simulation time, and the number of steps in + // which we want to reach the final time + unsigned int number_of_time_steps = 1; + double final_time = 1.0; + + // Instead of producing an output at every time step, we allow the user to + // select the frequency at which output is produced: + unsigned int output_frequency = 1; + + // We allow every grid to be refined independently. In this tutorial, no + // physics is resolved on the solid grid, and its velocity is given as a + // datum. However it relatively straight forward to incorporate some + // elasticity model in this tutorial, and transform it in a fully fledged + // FSI solver. + unsigned int initial_fluid_refinement = 3; + unsigned int initial_solid_refinement = 3; + unsigned int particle_insertion_refinement = 1; + + // The only two parameters used in the equations are the viscosity of the + // fluid, and the penalty term used in the Nitsche formulation: + double viscosity = 1.0; + double penalty_term = 1e3; + + // By default, we create a hyper_cube without colorisation, and we use + // homogenous Dirichlet boundary conditions. In this set we store the + // boundary ids to use when setting the boundary conditions: + std::list homogeneous_dirichlet_ids{0}; + + // We illustrate here another way to create a Triangulation from a parameter + // file, using the method GridGenerator::generate_from_name_and_arguments(), + // that takes the name of a function in the GridGenerator namespace, and its + // arguments as a single string representing the arguments as a tuple. + // + // The mechanism with which the arguments are parsed from and to a string is + // explained in detail in the Patterns::Tools::Convert class, which is + // used to translate from strings to most of the basic STL types (vectors, + // maps, tuples) and basic dealii types (Point, Tensor, BoundingBox, etc.). + // + // In general objects that can be represented by rank 1 uniform elements + // (i.e., std::vector, Point, std::set, etc.) are comma + // separated. Additional ranks take a semicolon, allowing you to parse + // strings into objects of type `std::vector>`, or, + // for example, `std::vector>`, as `0.0, 0.1; 0.1, 0.2`. This + // string could be interpreted as a vector of two Point objects, or a vector + // of vector of doubles. + // + // When the entries are not uniform, as in the tuple case, we use a colon + // to separate the various entries. For example, a string like `5: 0.1, 0.2` + // could be used to parse an object of type `std::pair>` or a + // `std::tuple>`. + // + // In our case most of the arguments are Point objects (representing + // centers, corners, subdivision elements, etc.), integer values (number of + // subdivisions), double values (radius, lengths, etc.), or boolean options + // (such as the `colorize` option that many GridGenerator functions take). + // + // In the example below, we set reasonable default values, but these can be + // changed at run time by selecting any other supported function of the + // GridGenerator namespace. + // + // We do this for each of the generated grids, to be as generic as possible: + std::string name_of_grid1 = "hyper_cube"; + std::string arguments_for_grid1 = "-1: 1: false"; + std::string name_of_grid2 = "hyper_rectangle"; + std::string arguments_for_grid2 = dim == 2 ? "-.5, -.1: .5, .1: false" : "-.5, -.1, -.1: .5, .1, .1: false"; std::string name_of_particle_grid = "hyper_ball"; std::string arguments_for_particle_grid = dim == 2 ? "0.3, 0.3: 0.1: false" : "0.3, 0.3, 0.3 : 0.1: false"; - // Refinement parameters + // Similarly, we allow for different local refinement strategies. In + // particular, we limit the maximum number of refinement levels, in order + // to control the minimum size of the fluid grid, and guarantee that it is + // compatible with the solid grid, and we perform local refinement based + // on standard error estimators on the fluid velocity field. + // + // We permit the user to choose between the + // two most common refinement strategies, namely `fixed_number` or + // `fixed_fraction`, that refer to the methods + // GridRefinement::refine_and_coarsen_fixed_fraction() and + // GridRefinement::refine_and_coarsen_fixed_number(). + // + // Refinement may be done every few time steps, instead of continuosly, and + // we control this value by the `refinement_frequency` parameter: int max_level_refinement = 5; std::string refinement_strategy = "fixed_fraction"; double coarsening_fraction = 0.3; double refinement_fraction = 0.3; unsigned int max_cells = 1000; - int mod_refinement = 5; - int mod_output = 1; - + int refinement_frequency = 5; + + // These two functions are used to control the source term of Stokes flow + // and the angular velocity at which we move solid. In a more realistic + // simulation, the solid velocity or its deformation would come from the + // solution of an auxiliary problem on the solid domain. In this example + // step we leave this part aside, and simply impose a fixed rotational + // velocity field on the immersed solid, governed by function that can be + // specified in the parameter file: mutable ParameterAcceptorProxy> rhs; mutable ParameterAcceptorProxy> angular_velocity; }; // namespace Step70 + // Once the angular velocity is provided as a Function object, we reconstruct + // the pointwise solid velocity thrugh the following class. template class SolidVelocity : public Function { public: SolidVelocity(const Functions::ParsedFunction &angular_velocity) : angular_velocity(angular_velocity) - {} + { + static_assert(spacedim > 1, + "Cannot instatiate SolidVelocity for spacedim == 1"); + } virtual double value(const Point &p, unsigned int component = 0) const @@ -318,8 +376,7 @@ namespace Step70 velocity = cross_product_3d(p, omega); } - - if (spacedim == 2) + else if (spacedim == 2) { double omega = angular_velocity.value(p, 0); @@ -334,7 +391,10 @@ namespace Step70 const Functions::ParsedFunction &angular_velocity; }; - + // Similarly, we assume that the incremental solid displacement can be + // computed simply by a one step time integration process (here using a + // trivial forward Euler method), so that at each time step, the solid simply + // displaces by `v*dt`. template class SolidDisplacement : public Function { @@ -345,7 +405,10 @@ namespace Step70 : Function(spacedim) , angular_velocity(angular_velocity) , time_step(time_step) - {} + { + static_assert(spacedim > 1, + "Cannot instatiate SolidDisplacement for spacedim == 1"); + } virtual double value(const Point &p, unsigned int component = 0) const @@ -370,6 +433,7 @@ namespace Step70 double time_step; }; + // We are now ready to introduce the main class of our tutorial program. template class StokesImmersedProblem { @@ -377,18 +441,44 @@ namespace Step70 StokesImmersedProblem( const StokesImmersedProblemParameters &par); + // As usual, we leave a single public entry point to the user: the run + // method. Everything else is left private, and accessed through the run + // method itself. void run(); private: void make_grid(); + + // These two methods are new w.r.t. previous examples, and initiliaze the + // ParticleHandler objects used in this class. We have two such objects: one + // is a passive tracer, used to plot the trajectories of fluid particles, + // while the the other is composed of the actual solid quadrature points, + // and represent material particles of the solid. void setup_tracer_particles(); void setup_solid_particles(); + + // The setup is split in two parts: create all objects that are needed once + // per simulation, void initial_setup(); + // followed by all objects that need to be reinitialized at every refinement + // step. void setup_dofs(); + + // The assembly rutine is identical to other Stokes assembly rutines, void assemble_stokes_system(); + // with the exception of the Nistche restriction part, which exploits one of + // the particle handlers to integrate on a non-matching part of the fluid + // domain, corresponding to the position of the solid. void assemble_nitche_restriction(); + void solve(); - void refine_grid(); + + // The refine_and_transfer() method is called only every + // `refinement_frequency` steps, and makes sure that all the fields + // that were computed on the time step before refinement are transfered + // correctly to the new grid. This includes vector fields, as well as + // particle information. + void refine_and_transfer(); void output_results(const unsigned int cycle, const double time) const; void @@ -397,29 +487,52 @@ namespace Step70 const unsigned int iter, const double time) const; + // As noted before, make sure we cannot modify this object from within this + // class, by making it a const reference. const StokesImmersedProblemParameters ∥ MPI_Comm mpi_communicator; - std::unique_ptr> fe1; - std::unique_ptr> fe2; - - parallel::distributed::Triangulation tria1; - parallel::distributed::Triangulation tria2; - - DoFHandler dh1; - DoFHandler dh2; - - std::unique_ptr> mapping2; - - std::vector owned1; - std::vector owned2; - - std::vector relevant1; - std::vector relevant2; - - IndexSet owned_tracer_particles; - IndexSet relevant_tracer_particles; + // For the current implemenation, only `fluid_fe` would be really necessary. + // For completeness, and to allow easy extension, we keep also the + // `solid_fe` around, which is however initialized to a FE_Nothing finite + // element space, i.e., one that has no degrees of freedom. + // + // We declare both finite element spaces as unique pointers, to allow their + // generation after StokesImmersedProblemParameters has been initialized. In + // particular, we assume that they are filled only after initial_setup() has + // been called. + std::unique_ptr> fluid_fe; + std::unique_ptr> solid_fe; + + // This is one of the main novelty w.r.t. the tutorial step-60. Here we + // assume that both the solid and the fluid are fully distributed + // triangulations. This allows the problem to scale to a very large number + // of degrees of freedom, at the cost of communicating all the overlapping + // regions between non matching triangulations. This is especially tricky, + // since we make no assumptions on the relative position or distribution of + // the various subdomains. In particular, we assume that ever process owns + // only a part of the solid_tria, and only a part of the fluid_tria, not + // necessarily in the same physical region, and not necessarily overlapping. + // + // In order to couple the overlapping regions, we exploit the facilities + // implemented in the ParticleHandler class. + parallel::distributed::Triangulation fluid_tria; + parallel::distributed::Triangulation solid_tria; + + DoFHandler fluid_dh; + DoFHandler solid_dh; + + std::unique_ptr> solid_mapping; + + // Similarly to how things are done in step-32, we use a block system to + // treat the Stokes part of the problem, and follow very closely what was + // done there. + std::vector fluid_owned_dofs; + std::vector solid_owned_dofs; + + std::vector fluid_relevant_dofs; + std::vector solid_relevant_dofs; AffineConstraints constraints; @@ -431,9 +544,40 @@ namespace Step70 LA::MPI::BlockVector locally_relevant_solution; LA::MPI::BlockVector system_rhs; + // For every tracer particle, we need to compute the velocity field in its + // current position, and update its position using a discrete time stepping + // scheme. We do this using distributed linear algebra objects, where the + // owner of a particle is set to be equal to the process that generated that + // particle at time t=0. This information is stored for every process in the + // `owned_tracer_particles` IndexSet, that indicates which particles the + // current process owns. + // + // Once the particles have been distributed around to match the process that + // owns the region where the particle lives, we will need read access from + // that process on the corresponding velocity field. We achieve this by + // filling a read only velocity vector field, that contains the relevant + // information in ghost entries. This is achieved using the + // `relevant_tracer_particles` IndexSet, that keeps track of how things + // change during the simulation, i.e., it keeps track of where particles + // that I own have ended up being, and who owns the particles that ended up + // in my subdomain. + // + // While this is not the most efficient strategy, we keep it this way to + // illustrate how things would work in a real FSI problem. If a particle + // is linked to a specific solid degree of freedom, we are not free to + // choose who owns it, and we have to communicate this information around. + // We illustrate this here, and show that the communication pattern is + // point-to-point, and negligible in terms of total cost of the algorithm. + IndexSet owned_tracer_particles; + IndexSet relevant_tracer_particles; + + // These vectors are used to store the particles velocities (read-only, with + // ghost entries) and their displacement (read/write, no ghost entries). LA::MPI::Vector tracer_particle_velocities; LA::MPI::Vector relevant_tracer_particle_displacements; + // We fix once the quadrature formula that is used to integrate the solid + // domain. std::unique_ptr> quadrature_formula; Particles::ParticleHandler tracer_particle_handler; @@ -450,16 +594,16 @@ namespace Step70 const StokesImmersedProblemParameters &par) : par(par) , mpi_communicator(MPI_COMM_WORLD) - , tria1(mpi_communicator, - typename Triangulation::MeshSmoothing( - Triangulation::smoothing_on_refinement | - Triangulation::smoothing_on_coarsening)) - , tria2(mpi_communicator, - typename Triangulation::MeshSmoothing( - Triangulation::smoothing_on_refinement | - Triangulation::smoothing_on_coarsening)) - , dh1(tria1) - , dh2(tria2) + , fluid_tria(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , solid_tria(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fluid_dh(fluid_tria) + , solid_dh(solid_tria) , pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) , computing_timer(mpi_communicator, @@ -472,15 +616,15 @@ namespace Step70 template void StokesImmersedProblem::make_grid() { - GridGenerator::generate_from_name_and_arguments(tria1, + GridGenerator::generate_from_name_and_arguments(fluid_tria, par.name_of_grid1, par.arguments_for_grid1); - tria1.refine_global(par.initial_fluid_refinement); + fluid_tria.refine_global(par.initial_fluid_refinement); - GridGenerator::generate_from_name_and_arguments(tria2, + GridGenerator::generate_from_name_and_arguments(solid_tria, par.name_of_grid2, par.arguments_for_grid2); - tria2.refine_global(par.initial_solid_refinement); + solid_tria.refine_global(par.initial_solid_refinement); } template @@ -503,7 +647,7 @@ namespace Step70 particles_dof_handler.distribute_dofs(particles_fe); // Create the particle handler associated with the fluid triangulation - tracer_particle_handler.initialize(tria1, + tracer_particle_handler.initialize(fluid_tria, StaticMappingQ1::mapping); @@ -511,7 +655,7 @@ namespace Step70 // The generation of the global bounding boxes requires an all-to-all // communication auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box( - tria1, IteratorFilters::LocallyOwnedCell()); + fluid_tria, IteratorFilters::LocallyOwnedCell()); auto global_bounding_boxes = Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box); @@ -529,11 +673,11 @@ namespace Step70 relevant_tracer_particles = owned_tracer_particles; // Now make sure that upon refinement, particles are correctly transferred - tria1.signals.pre_distributed_refinement.connect(std::bind( + fluid_tria.signals.pre_distributed_refinement.connect(std::bind( &Particles::ParticleHandler::register_store_callback_function, &tracer_particle_handler)); - tria1.signals.post_distributed_refinement.connect(std::bind( + fluid_tria.signals.post_distributed_refinement.connect(std::bind( &Particles::ParticleHandler::register_load_callback_function, &tracer_particle_handler, @@ -546,28 +690,28 @@ namespace Step70 template void StokesImmersedProblem::setup_solid_particles() { - QGauss quadrature(fe1->degree + 1); + QGauss quadrature(fluid_fe->degree + 1); // In codimension one case, we store also the normal, else only the // quadrature weight. const unsigned int n_properties = (dim == spacedim) ? 1 : spacedim + 1; - solid_particle_handler.initialize(tria1, + solid_particle_handler.initialize(fluid_tria, StaticMappingQ1::mapping, n_properties); std::vector> quadrature_points_vec( - quadrature.size() * tria2.n_locally_owned_active_cells()); + quadrature.size() * solid_tria.n_locally_owned_active_cells()); std::vector> properties( - quadrature.size() * tria2.n_locally_owned_active_cells(), + quadrature.size() * solid_tria.n_locally_owned_active_cells(), std::vector(n_properties)); UpdateFlags flags = update_JxW_values | update_quadrature_points; if (spacedim > dim) flags |= update_normal_vectors; - FEValues fe_v(*fe2, quadrature, flags); + FEValues fe_v(*solid_fe, quadrature, flags); unsigned int cell_index = 0; - for (const auto &cell : dh2.active_cell_iterators()) + for (const auto &cell : solid_dh.active_cell_iterators()) if (cell->is_locally_owned()) { fe_v.reinit(cell); @@ -591,7 +735,7 @@ namespace Step70 // Distribute the local points to the processor that owns // them on the triangulation auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box( - tria1, IteratorFilters::LocallyOwnedCell()); + fluid_tria, IteratorFilters::LocallyOwnedCell()); auto global_bounding_boxes = Utilities::MPI::all_gather(mpi_communicator, my_bounding_box); @@ -603,11 +747,11 @@ namespace Step70 // Now make sure that upon refinement, particles are correctly transferred - tria1.signals.pre_distributed_refinement.connect(std::bind( + fluid_tria.signals.pre_distributed_refinement.connect(std::bind( &Particles::ParticleHandler::register_store_callback_function, &solid_particle_handler)); - tria1.signals.post_distributed_refinement.connect(std::bind( + fluid_tria.signals.post_distributed_refinement.connect(std::bind( &Particles::ParticleHandler::register_load_callback_function, &solid_particle_handler, @@ -622,7 +766,7 @@ namespace Step70 { TimerOutput::Scope t(computing_timer, "initial setup"); - fe1 = + fluid_fe = std::make_unique>(FE_Q(par.velocity_degree), spacedim, FE_Q(par.velocity_degree - @@ -630,8 +774,8 @@ namespace Step70 1); - fe2 = std::make_unique>(); - dh2.distribute_dofs(*fe2); + solid_fe = std::make_unique>(); + solid_dh.distribute_dofs(*solid_fe); quadrature_formula = std::make_unique>(par.velocity_degree + 1); } @@ -642,43 +786,44 @@ namespace Step70 { TimerOutput::Scope t(computing_timer, "setup dofs"); - dh1.distribute_dofs(*fe1); + fluid_dh.distribute_dofs(*fluid_fe); std::vector stokes_sub_blocks(dim + 1, 0); stokes_sub_blocks[dim] = 1; - DoFRenumbering::component_wise(dh1, stokes_sub_blocks); + DoFRenumbering::component_wise(fluid_dh, stokes_sub_blocks); auto dofs_per_block = - DoFTools::count_dofs_per_fe_block(dh1, stokes_sub_blocks); + DoFTools::count_dofs_per_fe_block(fluid_dh, stokes_sub_blocks); const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1]; - pcout << " Number of degrees of freedom: " << dh1.n_dofs() << " (" << n_u - << '+' << n_p << " -- " << solid_particle_handler.n_global_particles() - << '+' << tracer_particle_handler.n_global_particles() << ')' - << std::endl; + pcout << " Number of degrees of freedom: " << fluid_dh.n_dofs() << " (" + << n_u << '+' << n_p << " -- " + << solid_particle_handler.n_global_particles() << '+' + << tracer_particle_handler.n_global_particles() << ')' << std::endl; - owned1.resize(2); - owned1[0] = dh1.locally_owned_dofs().get_view(0, n_u); - owned1[1] = dh1.locally_owned_dofs().get_view(n_u, n_u + n_p); + fluid_owned_dofs.resize(2); + fluid_owned_dofs[0] = fluid_dh.locally_owned_dofs().get_view(0, n_u); + fluid_owned_dofs[1] = + fluid_dh.locally_owned_dofs().get_view(n_u, n_u + n_p); IndexSet locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs(dh1, locally_relevant_dofs); - relevant1.resize(2); - relevant1[0] = locally_relevant_dofs.get_view(0, n_u); - relevant1[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p); + DoFTools::extract_locally_relevant_dofs(fluid_dh, locally_relevant_dofs); + fluid_relevant_dofs.resize(2); + fluid_relevant_dofs[0] = locally_relevant_dofs.get_view(0, n_u); + fluid_relevant_dofs[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p); { constraints.reinit(locally_relevant_dofs); FEValuesExtractors::Vector velocities(0); - DoFTools::make_hanging_node_constraints(dh1, constraints); - VectorTools::interpolate_boundary_values(dh1, - 0, - ZeroFunction(spacedim + - 1), - constraints, - fe1->component_mask(velocities)); + DoFTools::make_hanging_node_constraints(fluid_dh, constraints); + VectorTools::interpolate_boundary_values( + fluid_dh, + 0, + ZeroFunction(spacedim + 1), + constraints, + fluid_fe->component_mask(velocities)); constraints.close(); } @@ -697,15 +842,16 @@ namespace Step70 BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block); - DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false); + DoFTools::make_sparsity_pattern( + fluid_dh, coupling, dsp, constraints, false); SparsityTools::distribute_sparsity_pattern( dsp, - dh1.compute_locally_owned_dofs_per_processor(), + fluid_dh.compute_locally_owned_dofs_per_processor(), mpi_communicator, locally_relevant_dofs); - system_matrix.reinit(owned1, dsp, mpi_communicator); + system_matrix.reinit(fluid_owned_dofs, dsp, mpi_communicator); } { @@ -721,18 +867,21 @@ namespace Step70 BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block); - DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false); + DoFTools::make_sparsity_pattern( + fluid_dh, coupling, dsp, constraints, false); SparsityTools::distribute_sparsity_pattern( dsp, - dh1.compute_locally_owned_dofs_per_processor(), + fluid_dh.compute_locally_owned_dofs_per_processor(), mpi_communicator, locally_relevant_dofs); - preconditioner_matrix.reinit(owned1, dsp, mpi_communicator); + preconditioner_matrix.reinit(fluid_owned_dofs, dsp, mpi_communicator); } - locally_relevant_solution.reinit(owned1, relevant1, mpi_communicator); - system_rhs.reinit(owned1, mpi_communicator); - solution.reinit(owned1, mpi_communicator); + locally_relevant_solution.reinit(fluid_owned_dofs, + fluid_relevant_dofs, + mpi_communicator); + system_rhs.reinit(fluid_owned_dofs, mpi_communicator); + solution.reinit(fluid_owned_dofs, mpi_communicator); } @@ -747,13 +896,13 @@ namespace Step70 TimerOutput::Scope t(computing_timer, "Stokes_assembly"); - FEValues fe_values(*fe1, + FEValues fe_values(*fluid_fe, *quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); - const unsigned int dofs_per_cell = fe1->dofs_per_cell; + const unsigned int dofs_per_cell = fluid_fe->dofs_per_cell; const unsigned int n_q_points = quadrature_formula->size(); FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); @@ -771,7 +920,7 @@ namespace Step70 const FEValuesExtractors::Vector velocities(0); const FEValuesExtractors::Scalar pressure(spacedim); - for (const auto &cell : dh1.active_cell_iterators()) + for (const auto &cell : fluid_dh.active_cell_iterators()) if (cell->is_locally_owned()) { cell_matrix = 0; @@ -805,7 +954,7 @@ namespace Step70 } const unsigned int component_i = - fe1->system_to_component_index(i).first; + fluid_fe->system_to_component_index(i).first; cell_rhs(i) += fe_values.shape_value(i, q) * rhs_values[q](component_i) * fe_values.JxW(q); } @@ -838,19 +987,20 @@ namespace Step70 SolidVelocity solid_velocity(par.angular_velocity); - std::vector dof_indices1(fe1->dofs_per_cell); + std::vector dof_indices1(fluid_fe->dofs_per_cell); - FullMatrix local_matrix(fe1->dofs_per_cell, fe1->dofs_per_cell); - dealii::Vector local_rhs(fe1->dofs_per_cell); + FullMatrix local_matrix(fluid_fe->dofs_per_cell, + fluid_fe->dofs_per_cell); + dealii::Vector local_rhs(fluid_fe->dofs_per_cell); auto particle = solid_particle_handler.begin(); while (particle != solid_particle_handler.end()) { local_matrix = 0; local_rhs = 0; - const auto &cell = particle->get_surrounding_cell(tria1); + const auto &cell = particle->get_surrounding_cell(fluid_tria); const auto &dh_cell = - typename DoFHandler::cell_iterator(*cell, &dh1); + typename DoFHandler::cell_iterator(*cell, &fluid_dh); dh_cell->get_dof_indices(dof_indices1); const auto pic = solid_particle_handler.particles_in_cell(cell); @@ -861,23 +1011,24 @@ namespace Step70 const auto real_q = p.get_location(); const auto properties = p.get_properties(); const auto &JxW = properties[0]; - for (unsigned int i = 0; i < fe1->dofs_per_cell; ++i) + for (unsigned int i = 0; i < fluid_fe->dofs_per_cell; ++i) { - const auto comp_i = fe1->system_to_component_index(i).first; + const auto comp_i = + fluid_fe->system_to_component_index(i).first; if (comp_i < spacedim) { - for (unsigned int j = 0; j < fe1->dofs_per_cell; ++j) + for (unsigned int j = 0; j < fluid_fe->dofs_per_cell; ++j) { const auto comp_j = - fe1->system_to_component_index(j).first; + fluid_fe->system_to_component_index(j).first; if (comp_i == comp_j) local_matrix(i, j) += - par.penalty_term * fe1->shape_value(i, ref_q) * - fe1->shape_value(j, ref_q) * JxW; + par.penalty_term * fluid_fe->shape_value(i, ref_q) * + fluid_fe->shape_value(j, ref_q) * JxW; } local_rhs(i) += par.penalty_term * solid_velocity.value(real_q, comp_i) * - fe1->shape_value(i, ref_q) * JxW; + fluid_fe->shape_value(i, ref_q) * JxW; } } } @@ -948,7 +1099,7 @@ namespace Step70 locally_relevant_solution = solution; const double mean_pressure = - VectorTools::compute_mean_value(dh1, + VectorTools::compute_mean_value(fluid_dh, QGauss(par.velocity_degree + 2), locally_relevant_solution, spacedim); @@ -959,23 +1110,23 @@ namespace Step70 template - void StokesImmersedProblem::refine_grid() + void StokesImmersedProblem::refine_and_transfer() { TimerOutput::Scope t(computing_timer, "refine"); const FEValuesExtractors::Vector velocity(0); - Vector error_per_cell(tria1.n_active_cells()); - KellyErrorEstimator::estimate(dh1, + Vector error_per_cell(fluid_tria.n_active_cells()); + KellyErrorEstimator::estimate(fluid_dh, QGauss(par.velocity_degree + 1), {}, locally_relevant_solution, error_per_cell, - fe1->component_mask(velocity)); + fluid_fe->component_mask(velocity)); if (par.refinement_strategy == "fixed_fraction") { parallel::distributed::GridRefinement:: - refine_and_coarsen_fixed_fraction(tria1, + refine_and_coarsen_fixed_fraction(fluid_tria, error_per_cell, par.refinement_fraction, par.coarsening_fraction); @@ -983,22 +1134,22 @@ namespace Step70 else if (par.refinement_strategy == "fixed_number") { parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( - tria1, + fluid_tria, error_per_cell, par.refinement_fraction, par.coarsening_fraction, par.max_cells); } - for (const auto &cell : tria1.active_cell_iterators()) + for (const auto &cell : fluid_tria.active_cell_iterators()) if (cell->refine_flag_set() && cell->level() == par.max_level_refinement) cell->clear_refine_flag(); parallel::distributed::SolutionTransfer transfer( - dh1); - tria1.prepare_coarsening_and_refinement(); + fluid_dh); + fluid_tria.prepare_coarsening_and_refinement(); transfer.prepare_for_coarsening_and_refinement(locally_relevant_solution); - tria1.execute_coarsening_and_refinement(); + fluid_tria.execute_coarsening_and_refinement(); setup_dofs(); transfer.interpolate(solution); constraints.distribute(solution); @@ -1023,20 +1174,20 @@ namespace Step70 DataComponentInterpretation::component_is_scalar); DataOut data_out; - data_out.attach_dof_handler(dh1); + data_out.attach_dof_handler(fluid_dh); data_out.add_data_vector(locally_relevant_solution, solution_names, DataOut::type_dof_data, data_component_interpretation); LA::MPI::BlockVector interpolated; - interpolated.reinit(owned1, MPI_COMM_WORLD); - VectorTools::interpolate(dh1, + interpolated.reinit(fluid_owned_dofs, MPI_COMM_WORLD); + VectorTools::interpolate(fluid_dh, ConstantFunction(1.0, spacedim + 1), interpolated); - LA::MPI::BlockVector interpolated_relevant(owned1, - relevant1, + LA::MPI::BlockVector interpolated_relevant(fluid_owned_dofs, + fluid_relevant_dofs, MPI_COMM_WORLD); interpolated_relevant = interpolated; { @@ -1049,9 +1200,9 @@ namespace Step70 } - Vector subdomain(tria1.n_active_cells()); + Vector subdomain(fluid_tria.n_active_cells()); for (unsigned int i = 0; i < subdomain.size(); ++i) - subdomain(i) = tria1.locally_owned_subdomain(); + subdomain(i) = fluid_tria.locally_owned_subdomain(); data_out.add_data_vector(subdomain, "subdomain"); data_out.build_patches(); @@ -1134,7 +1285,7 @@ namespace Step70 } { TimerOutput::Scope t(computing_timer, "Set tracer particle motion"); - interpolate_field_on_particles(dh1, + interpolate_field_on_particles(fluid_dh, tracer_particle_handler, locally_relevant_solution, tracer_particle_velocities, @@ -1160,7 +1311,7 @@ namespace Step70 assemble_nitche_restriction(); solve(); - if (cycle % par.mod_output == 0) + if (cycle % par.output_frequency == 0) { static unsigned int output_cycle = 0; output_results(output_cycle, time); @@ -1180,11 +1331,104 @@ namespace Step70 } ++output_cycle; } - if (cycle % par.mod_refinement == 0 && + if (cycle % par.refinement_frequency == 0 && cycle != par.number_of_time_steps - 1) - refine_grid(); + refine_and_transfer(); } } + + template + StokesImmersedProblemParameters::StokesImmersedProblemParameters() + : ParameterAcceptor("Stokes Immersed Problem/") + , rhs("Right hand side", spacedim + 1) + , angular_velocity("Angular velocity", spacedim == 3 ? spacedim : 1) + { + // We split the parameters in various cathegories, by putting them in + // different sections of the ParameterHandler class. We begin by declaring + // all the global parameters used by StokesImmersedProblem in the global + // scope: + add_parameter( + "Velocity degree", velocity_degree, "", this->prm, Patterns::Integer(1)); + + add_parameter("Number of time steps", number_of_time_steps); + add_parameter("Output frequency", output_frequency); + + add_parameter("Final time", final_time); + + add_parameter("Viscosity", viscosity); + + add_parameter("Nitsche penalty term", penalty_term); + + add_parameter("Initial fluid refinement", + initial_fluid_refinement, + "Initial mesh refinement used for the fluid domain Omega"); + + add_parameter("Initial solid refinement", + initial_solid_refinement, + "Initial mesh refinement used for the solid domain Gamma"); + + add_parameter( + "Particle insertion refinement", + particle_insertion_refinement, + "Refinement of the volumetric mesh used to insert the particles"); + + add_parameter( + "Homogeneous Dirichlet boundary ids", + homogeneous_dirichlet_ids, + "Boundary Ids over which homogeneous Dirichlet boundary conditions are applied"); + + // Next section is dedicated to the parameters used to create the various + // grids. We will need three different triangulations: `Grid one` is used + // to define the fluid domain, `Grid two` defines the solid domain, and + // `Particle grid` is used to distribute some tracer particles, that are + // advected with the velocity and only used as passive tracers. + enter_my_subsection(this->prm); + this->prm.enter_subsection("Grid generation"); + this->prm.add_parameter("Grid one generator", name_of_grid1); + this->prm.add_parameter("Grid one generator arguments", + arguments_for_grid1); + + this->prm.add_parameter("Grid two generator", name_of_grid2); + this->prm.add_parameter("Grid two generator arguments", + arguments_for_grid2); + + this->prm.add_parameter("Particle grid generator", name_of_particle_grid); + this->prm.add_parameter("Particle grid generator arguments", + arguments_for_particle_grid); + this->prm.leave_subsection(); + + leave_my_subsection(this->prm); + + + + enter_my_subsection(this->prm); + this->prm.enter_subsection("Refinement and remeshing"); + this->prm.add_parameter("Refinement step frequency", refinement_frequency); + this->prm.add_parameter("Refinement maximal level", max_level_refinement); + this->prm.add_parameter("Refinement strategy", + refinement_strategy, + "", + Patterns::Selection("fixed_fraction|fixed_number")); + this->prm.add_parameter("Refinement coarsening fraction", + coarsening_fraction); + this->prm.add_parameter("Refinement fraction", refinement_fraction); + this->prm.add_parameter("Maximum number of cells", max_cells); + + this->prm.leave_subsection(); + leave_my_subsection(this->prm); + + // correct the default dimension for the functions + rhs.declare_parameters_call_back.connect([&]() { + Functions::ParsedFunction::declare_parameters(this->prm, + spacedim + 1); + }); + angular_velocity.declare_parameters_call_back.connect([&]() { + Functions::ParsedFunction::declare_parameters( + this->prm, spacedim == 3 ? spacedim : 1); + }); + } + } // namespace Step70 @@ -1199,9 +1443,6 @@ int main(int argc, char *argv[]) Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); StokesImmersedProblemParameters<2> par; - par.declare_all_parameters(); - std::ofstream out("default.prm"); - par.prm.print_parameters(out, ParameterHandler::ShortText); ParameterAcceptor::initialize("parameters.prm", "used_parameters.prm"); StokesImmersedProblem<2> problem(par); -- 2.39.5