From: Luca Heltai Date: Tue, 12 May 2020 22:31:42 +0000 (+0200) Subject: Some minor changes. Add support for output directory. X-Git-Tag: v9.2.0-rc2~3^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a7991f53f6d54e40d0b952bfbce414fd2021dc6d;p=dealii.git Some minor changes. Add support for output directory. --- diff --git a/examples/step-70/doc/intro.dox b/examples/step-70/doc/intro.dox index 31e51e68e7..15b8a1e9e2 100644 --- a/examples/step-70/doc/intro.dox +++ b/examples/step-70/doc/intro.dox @@ -1,6 +1,7 @@
-This program was contributed by Luca Heltai (SISSA, Trieste) and Bruno Blais (Polytechnique Montréal) +This program was contributed by Luca Heltai (SISSA, Trieste), Bruno Blais (Polytechnique Montréal), +and Rene Gassmoeller (UC Davis) Change this!!!! diff --git a/examples/step-70/doc/results.dox b/examples/step-70/doc/results.dox index 43a12c9c76..7ef60bbfab 100644 --- a/examples/step-70/doc/results.dox +++ b/examples/step-70/doc/results.dox @@ -1,14 +1,15 @@

Results

-The directory in which this program is run contains a number of sample parameter -files that you can use to reproduce the results presented in this section. If -you do not specify a parameter file as an argument on command line, the program -will try to read the file "parameters.prm" by default, and will execute the two -dimensional version of the code. As explained before, if your file name contains -the string "23", then the program will run a three dimensional problem, with -immersed solid of co-dimension one. If it contains the string "3", it will run -a three dimensional problem, with immersed solid of co-dimension zero, otherwise -it will run a two dimensional problem with immersed solid of co-dimension zero. +The directory in which this program is run contains a number of sample +parameter files that you can use to reproduce the results presented in this +section. If you do not specify a parameter file as an argument on the command +line, the program will try to read the file "parameters.prm" by default, and +will execute the two dimensional version of the code. As explained before, if +your file name contains the string "23", then the program will run a three +dimensional problem, with immersed solid of co-dimension one. If it contains +the string "3", it will run a three dimensional problem, with immersed solid of +co-dimension zero, otherwise it will run a two dimensional problem with +immersed solid of co-dimension zero. Regardless of the specific parameter file name, if the specified file does not exist, when you execute the program you will get an exception that no such file @@ -32,11 +33,11 @@ Aborting! @endcode However, as the error message already states, the code that triggers the -exception will also generate the specified file ("parameters.prm" in this -case) that simply contains the default values for all parameters this program -cares about (for the correct dimension and co-dimension, according to the -wether a string "23" or "3" is contained in the file name). -By inspection of the default parameter file, we see the following: +exception will also generate the specified file ("parameters.prm" in this case) +that simply contains the default values for all parameters this program cares +about (for the correct dimension and co-dimension, according to the wether a +string "23" or "3" is contained in the file name). By inspection of the default +parameter file, we see the following: @code # Listing of Parameters @@ -55,12 +56,13 @@ subsection Stokes Immersed Problem # Initial mesh refinement used for the solid domain Gamma set Initial solid refinement = 5 - set Nitsche penalty term = 10 - set Number of time steps = 11 + set Nitsche penalty term = 100 + set Number of time steps = 501 + set Output directory = . set Output frequency = 1 # Refinement of the volumetric mesh used to insert the particles - set Particle insertion refinement = 1 + set Particle insertion refinement = 3 set Velocity degree = 2 set Viscosity = 1 @@ -70,12 +72,12 @@ subsection Stokes Immersed Problem # that describes the function, rather than having to use its numeric value # everywhere the constant appears. These values can be defined using this # parameter, in the form `var1=value1, var2=value2, ...'. - # + # # A typical example would be to set this runtime parameter to # `pi=3.1415926536' and then use `pi' in the expression of the actual # formula. (That said, for convenience this class actually defines both # `pi' and `Pi' by default, but you get the idea.) - set Function constants = + set Function constants = # The formula that denotes the function you want to evaluate for # particular values of the independent variables. This expression may @@ -86,7 +88,7 @@ subsection Stokes Immersed Problem # true, and to the third argument otherwise. For a full overview of # possible expressions accepted see the documentation of the muparser # library at http://muparser.beltoforion.de/. - # + # # If the function you are describing represents a vector-valued function # with multiple components, then separate the expressions for individual # components by a semicolon. @@ -130,12 +132,12 @@ subsection Stokes Immersed Problem # that describes the function, rather than having to use its numeric value # everywhere the constant appears. These values can be defined using this # parameter, in the form `var1=value1, var2=value2, ...'. - # + # # A typical example would be to set this runtime parameter to # `pi=3.1415926536' and then use `pi' in the expression of the actual # formula. (That said, for convenience this class actually defines both # `pi' and `Pi' by default, but you get the idea.) - set Function constants = + set Function constants = # The formula that denotes the function you want to evaluate for # particular values of the independent variables. This expression may @@ -146,7 +148,7 @@ subsection Stokes Immersed Problem # true, and to the third argument otherwise. For a full overview of # possible expressions accepted see the documentation of the muparser # library at http://muparser.beltoforion.de/. - # + # # If the function you are describing represents a vector-valued function # with multiple components, then separate the expressions for individual # components by a semicolon. @@ -169,9 +171,12 @@ subsection Stokes Immersed Problem end @endcode -If you now run the program, you will get a file called `parameters_used.prm`, -containing a shorter version of the above parameters (without comments and -documentation), documenting all parameters that were used to run your program: +If you now run the program, you will get a file called `parameters_22.prm` in +the directory specified by the parameter `Output directory` (which defaults to +the current directory) containing a shorter version of the above parameters +(without comments and documentation), documenting all parameters that were used +to run your program: + @code subsection Stokes Immersed Problem set Final time = 1 @@ -181,13 +186,14 @@ subsection Stokes Immersed Problem set Initial solid refinement = 5 set Nitsche penalty term = 100 set Number of time steps = 501 + set Output directory = . set Output frequency = 1 set Particle insertion refinement = 3 set Velocity degree = 2 set Viscosity = 1 subsection Angular velocity - set Function constants = - set Function expression = t < .500001 ? 6.283185 : -6.283185 # default: 0 + set Function constants = + set Function expression = t < .500001 ? 6.283185 : -6.283185 set Variable names = x,y,t end subsection Grid generation @@ -208,7 +214,7 @@ subsection Stokes Immersed Problem set Refinement strategy = fixed_fraction end subsection Right hand side - set Function constants = + set Function constants = set Function expression = 0; 0; 0 set Variable names = x,y,t end @@ -216,9 +222,11 @@ end @endcode The rationale behind creating first `parameters.prm` file (the first time the -program is run) and then a `parameters_used.prm` (every other times you run the -program), is because you may want to leave most parameters to their default -values, and only modify a handful of them. +program is run) and then a `output/parameters_22.prm` (every other times you +run the program), is because you may want to leave most parameters to their +default values, and only modify a handful of them, while still beeing able to +reproduce the results and inspect what parameters where used for a scpeficic +simulation. For example, you could use the following (perfectly valid) parameter file with this tutorial program: @@ -231,18 +239,19 @@ subsection Stokes Immersed Problem end @endcode -and you would run the program with Q3/Q2 Taylor-Hood finite elements, -for 101 steps, using a Nistche penalty of `10`, and leaving all the other -parameters to their default value. +and you would run the program with Q3/Q2 Taylor-Hood finite elements, for 101 +steps, using a Nistche penalty of `10`, and leaving all the other parameters to +their default value. You could then inspect all the other parameters in the +produced file `parameters_22.prm`.

Two dimensional test case

The default problem generates a co-dimension zero impeller, consisting of a rotating rectangular grid, where the rotation is for half a second in one direction, and half a second in the opposite direction, with constant angular -velocity equal to $\approx 2\Pi rad/s$. Consequently, the impeller does half a rotation -and returns to it's original position. The following animation displays -the velocity magnitude, the motion of the solid impeller and of the +velocity equal to $\approx 2\Pi rad/s$. Consequently, the impeller does half a +rotation and returns to it's original position. The following animation +displays the velocity magnitude, the motion of the solid impeller and of the tracer particles. @@ -330,8 +339,8 @@ simulation domain:

We see that, generally, the tracer particles have somewhat returned to their -original position, although they have been distorted by the flow field. -The following image compares the initial and the final position of the particles +original position, although they have been distorted by the flow field. The +following image compares the initial and the final position of the particles after 1s of flow.

@@ -342,13 +351,14 @@ after 1s of flow.

-In this case, we see that the tracer particles that were outside of the swept volume of the -impeller have returned very close to their initial position, whereas those in the swept -volume were slightly more deformed. This deformation is non-physical. It is caused by -the numerical error induced by the explicit Euler scheme used to advect the particles, -by the loss of accuracy due to the fictious domain and, finally, by the discretization -error on the Stokes equations. The first two errors are the leading cause of this deformation -and they could be alleviated by the use of a finer mesh and a lower time step. +In this case, we see that the tracer particles that were outside of the swept +volume of the impeller have returned very close to their initial position, +whereas those in the swept volume were slightly more deformed. This deformation +is non-physical. It is caused by the numerical error induced by the explicit +Euler scheme used to advect the particles, by the loss of accuracy due to the +fictious domain and, finally, by the discretization error on the Stokes +equations. The first two errors are the leading cause of this deformation and +they could be alleviated by the use of a finer mesh and a lower time step.

Three dimensional test case

@@ -442,16 +452,16 @@ points. The structure of the code already allows one to implement a two-way coupling, by exploiting the possibility to read values of the fluid velocity on the -quadrature points of the solid grid. For this to be more efficient in terms -of MPI communication patterns, one should maintain ownership of the quadrature -points on the solid processor that owns them. In the current code, it is -sufficient to define the IndexSet of the vectors used to exchange information -of the quadrature points by using the solid partition instead of the initial -fluid partition. - -This allows the combination of the technique used in this tutorial program -with those presented in the tutorial step-60 to solve a fluid structure -interaction problem with distributed Lagrange multipliers, on +quadrature points of the solid grid. For this to be more efficient in terms of +MPI communication patterns, one should maintain ownership of the quadrature +points on the solid processor that owns the cells where they have been created. +In the current code, it is sufficient to define the IndexSet of the vectors +used to exchange information of the quadrature points by using the solid +partition instead of the initial fluid partition. + +This allows the combination of the technique used in this tutorial program with +those presented in the tutorial step-60 to solve a fluid structure interaction +problem with distributed Lagrange multipliers, on parallel::distributed::Triangulation objects. The timings above show that the current preconditioning strategy does not work diff --git a/examples/step-70/parameters.prm b/examples/step-70/parameters.prm index 3234de62c4..04b42e9ccb 100644 --- a/examples/step-70/parameters.prm +++ b/examples/step-70/parameters.prm @@ -8,6 +8,7 @@ subsection Stokes Immersed Problem set Number of time steps = 501 set Velocity degree = 2 set Viscosity = 1 + set Output directory = results subsection Angular velocity set Function constants = set Function expression = t < .500001 ? 6.283185 : -6.283185 # default: 0 diff --git a/examples/step-70/step-70.cc b/examples/step-70/step-70.cc index 6f4f76e5f1..dfb2189f27 100644 --- a/examples/step-70/step-70.cc +++ b/examples/step-70/step-70.cc @@ -14,7 +14,7 @@ * --------------------------------------------------------------------- * - * Authors: Luca Heltai, Bruno Blais, 2020 + * Authors: Luca Heltai, Bruno Blais, Rene Gassmoeller, 2020 */ // @sect3{Include files} @@ -184,6 +184,9 @@ namespace Step70 angular_velocity.set_time(time); } + // this is where we write all the output files + std::string output_directory = "."; + // 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 @@ -350,10 +353,10 @@ namespace Step70 }; // Similarly, we assume that the solid position can be computed explicitly at - // each time step, exploiting the knoweledge of the agnular velocity. We - // perform 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`. + // each time step, exploiting the knoweledge of the angular velocity. We + // compute the exact position of the solid particle assuming that the solid is + // rotated by an amount equal to the time step multiplied by the angular + // velocity computed at the point `p`: template class SolidPosition : public Function { @@ -365,20 +368,20 @@ namespace Step70 , time_step(time_step) { static_assert(spacedim > 1, - "Cannot instantiate SolidDisplacement for spacedim == 1"); + "Cannot instantiate SolidPosition for spacedim == 1"); } virtual double value(const Point &p, unsigned int component = 0) const override { - Tensor<1, spacedim> displacement = p; + Point new_position = p; double dtheta = angular_velocity.value(p) * time_step; - displacement[0] = std::cos(dtheta) * p[0] - std::sin(dtheta) * p[1]; - displacement[1] = std::sin(dtheta) * p[0] + std::cos(dtheta) * p[1]; + new_position[0] = std::cos(dtheta) * p[0] - std::sin(dtheta) * p[1]; + new_position[1] = std::sin(dtheta) * p[0] + std::cos(dtheta) * p[1]; - return displacement[component]; + return new_position[component]; } void set_time_step(const double new_time_step) @@ -722,6 +725,7 @@ namespace Step70 } #else (void)ids_and_cad_file_names; + AssertThrow(false, ExcNotImplemented("Generation of the grid failed.")); #endif } @@ -900,7 +904,7 @@ namespace Step70 // The Particles::ParticleHandler class has a way to transfer information // from a cell to its children or to its parent upon refinement, without the // need to reconstruct the entire data structure. This is done by - // "registering" two callback functions to the triangulation. These + // registering two callback functions to the triangulation. These // functions will receive a signal when refinement is about to happen, and // when it has just happened, and will take care of transferring all // information to the newly refined grid with minimal computational cost. @@ -979,7 +983,21 @@ namespace Step70 } // We proceed in the same way we did with the tracer particles, reusing the - // computed bounding boxes. + // computed bounding boxes. However, we first check that the + // global_fluid_bounding_boxes object has been actually filled. This should + // certainly be the case here, since this method is called after the one + // that initializes the tracer particles. However, we want to make sure that + // if in the future someone decides (for whatever reason) to initialize + // first the solid particle handler, or to copy just this part of the + // tutorial, a meaningful exception is thrown when things don't work as + // expected + AssertThrow(!global_fluid_bounding_boxes.empty(), + ExcInternalError( + "I was expecting the " + "global_fluid_bounding_boxes to be filled at this stage. " + "Make sure you fill this vector before trying to use it " + "here. Bailing out.")); + // Since we have already stored the position of the quadrature point, // we can use these positions to insert the particles directly using // the solid_particle_handler instead of having to go through a @@ -1004,10 +1022,10 @@ namespace Step70 // We set up the finite element space and the quadrature formula to be // used throughout the step. For the fluid, we use Taylor-Hood elements (e.g. - // Q2-Q1). Since we do not solve any equation on the solid domain, an empty - // finite element space is generated. A natural extension of this program - // would be to solve a fluid structure interaction problem, which would - // require that the solid_fe use a non-empty FiniteElement. + // Q(P)-Q(P-1)). Since we do not solve any equation on the solid domain, an + // empty finite element space is generated. A natural extension of this + // program would be to solve a fluid structure interaction problem, which + // would require that the solid_fe use a non-empty FiniteElement. template void StokesImmersedProblem::initial_setup() { @@ -1032,6 +1050,13 @@ namespace Step70 std::make_unique>(par.velocity_degree + 1); solid_quadrature_formula = std::make_unique>(par.velocity_degree + 1); + + // Save the current parameter file in the output directory, for + // reproducibility + par.prm.print_parameters(par.output_directory + "/" + "parameters_" + + std::to_string(dim) + std::to_string(spacedim) + + ".prm", + ParameterHandler::Short); } @@ -1268,7 +1293,7 @@ namespace Step70 // of the cell in which the particle lies and then loop over all particles // within that cell. This enables us to skip the cells which do not contain // particles, yet to assemble the local matrix and rhs of each cell to apply - // the Nitsche restriction + // the Nitsche restriction. auto particle = solid_particle_handler.begin(); while (particle != solid_particle_handler.end()) { @@ -1335,8 +1360,8 @@ namespace Step70 } - // This function solves the linear system with MINRES with a block diagonal - // preconditioner and AMG for the two diagonal blocks as used in step-55. The + // This function solves the linear system with FGMRES with a block diagonal + // preconditioner and AMG for the two diagonal blocks. The // preconditioner applies a v cycle to the 0,0 block and a CG with the mass // matrix for the 1,1 block (the Schur complement). template @@ -1491,25 +1516,6 @@ namespace Step70 DataOut::type_dof_data, data_component_interpretation); - LA::MPI::BlockVector interpolated; - interpolated.reinit(fluid_owned_dofs, MPI_COMM_WORLD); - VectorTools::interpolate(fluid_dh, - ConstantFunction(1.0, spacedim + 1), - interpolated); - - LA::MPI::BlockVector interpolated_relevant(fluid_owned_dofs, - fluid_relevant_dofs, - MPI_COMM_WORLD); - interpolated_relevant = interpolated; - { - std::vector solution_names(spacedim, "ref_u"); - solution_names.emplace_back("ref_p"); - data_out.add_data_vector(interpolated_relevant, - solution_names, - DataOut::type_dof_data, - data_component_interpretation); - } - Vector subdomain(fluid_tria.n_active_cells()); for (unsigned int i = 0; i < subdomain.size(); ++i) @@ -1520,11 +1526,12 @@ namespace Step70 const std::string filename = "solution-" + Utilities::int_to_string(cycle) + ".vtu"; - data_out.write_vtu_in_parallel(filename, mpi_communicator); + data_out.write_vtu_in_parallel(par.output_directory + "/" + filename, + mpi_communicator); static std::vector> times_and_names; times_and_names.push_back(std::make_pair(time, filename)); - std::ofstream ofile("solution.pvd"); + std::ofstream ofile(par.output_directory + "/" + "solution.pvd"); DataOutBase::write_pvd_record(ofile, times_and_names); } @@ -1543,7 +1550,8 @@ namespace Step70 particles_out.build_patches(particles); const std::string filename = (fprefix + "-" + Utilities::int_to_string(iter) + ".vtu"); - particles_out.write_vtu_in_parallel(filename, mpi_communicator); + particles_out.write_vtu_in_parallel(par.output_directory + "/" + filename, + mpi_communicator); static std::map>> @@ -1552,13 +1560,13 @@ namespace Step70 times_and_names[fprefix].push_back(std::make_pair(time, filename)); else times_and_names[fprefix] = {std::make_pair(time, filename)}; - std::ofstream ofile(fprefix + ".pvd"); + std::ofstream ofile(par.output_directory + "/" + fprefix + ".pvd"); DataOutBase::write_pvd_record(ofile, times_and_names[fprefix]); } // This function orchestrates the entire simulation. It is very similar - // to the other transient steps. + // to the other time dependent tutorial programs. template void StokesImmersedProblem::run() { @@ -1575,8 +1583,10 @@ namespace Step70 ComponentMask velocity_mask(spacedim + 1, true); velocity_mask.set(spacedim, false); - const double time_step = par.final_time / (par.number_of_time_steps - 1); - double time = 0; + const double time_step = par.final_time / (par.number_of_time_steps - 1); + double time = 0; + unsigned int output_cycle = 0; + for (unsigned int cycle = 0; cycle < par.number_of_time_steps; ++cycle, time += time_step) { @@ -1594,14 +1604,20 @@ namespace Step70 setup_solid_particles(); tracer_particle_velocities.reinit(owned_tracer_particles, mpi_communicator); - output_results(0, time); + output_results(output_cycle, time); { TimerOutput::Scope t(computing_timer, "Output tracer particles"); - output_particles(tracer_particle_handler, "tracer", 0, time); + output_particles(tracer_particle_handler, + "tracer", + output_cycle, + time); } { TimerOutput::Scope t(computing_timer, "Output solid particles"); - output_particles(solid_particle_handler, "solid", 0, time); + output_particles(solid_particle_handler, + "solid", + output_cycle, + time); } } // On the other cycle, we displace the solid body to take into account @@ -1652,7 +1668,6 @@ namespace Step70 // particles, the tracer particles and the fluid domain. if (cycle % par.output_frequency == 0) { - static unsigned int output_cycle = 0; output_results(output_cycle, time); { TimerOutput::Scope t(computing_timer, "Output tracer particles"); @@ -1694,6 +1709,8 @@ namespace Step70 add_parameter("Number of time steps", number_of_time_steps); add_parameter("Output frequency", output_frequency); + add_parameter("Output directory", output_directory); + add_parameter("Final time", final_time); add_parameter("Viscosity", viscosity); @@ -1724,8 +1741,8 @@ namespace Step70 "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: `Fluid grid` - // is used to define the fluid domain, `Solid grid` defines the + // various grids. We will need three different triangulations: `Fluid + // grid` is used to define the fluid domain, `Solid grid` 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. @@ -1780,16 +1797,15 @@ namespace Step70 } // namespace Step70 -// The remainder of the code, the main function, is standard, with the exception -// of the handling of input parameter files. -// We allow the user to specify an optional parameter file as an argument to the -// program. If nothing is specified, we use the default file "parameters.prm", -// which is created if non existent. -// The file name is scanned for the the string "23" first, and "3" afterwards. -// If the filename contains the string "23", a the problem classes are -// instantiated with template arguments 2 and 3 respectively. If only the string -// "3" is found, then both template arguments are set to 3, otherwise both are -// set to 2. +// The remainder of the code, the main function, is standard, with the +// exception of the handling of input parameter files. We allow the user to +// specify an optional parameter file as an argument to the program. If +// nothing is specified, we use the default file "parameters.prm", which is +// created if non existent. The file name is scanned for the the string "23" +// first, and "3" afterwards. If the filename contains the string "23", a the +// problem classes are instantiated with template arguments 2 and 3 +// respectively. If only the string "3" is found, then both template arguments +// are set to 3, otherwise both are set to 2. int main(int argc, char *argv[]) { using namespace Step70; @@ -1802,23 +1818,16 @@ int main(int argc, char *argv[]) // Interpret the std::string prm_file; if (argc > 1) - { - prm_file = argv[1]; - } + prm_file = argv[1]; else - { - prm_file = "parameters.prm"; - } - - std::string used_prm_file = prm_file; - used_prm_file.insert(used_prm_file.find_last_of("."), "_used"); + prm_file = "parameters.prm"; // deduce the dimension of the problem from the name of the // parameter file specified at the command line if (prm_file.find("23") != std::string::npos) { StokesImmersedProblemParameters<2, 3> par; - ParameterAcceptor::initialize(prm_file, used_prm_file); + ParameterAcceptor::initialize(prm_file); StokesImmersedProblem<2, 3> problem(par); problem.run(); @@ -1826,7 +1835,7 @@ int main(int argc, char *argv[]) else if (prm_file.find("3") != std::string::npos) { StokesImmersedProblemParameters<3> par; - ParameterAcceptor::initialize(prm_file, used_prm_file); + ParameterAcceptor::initialize(prm_file); StokesImmersedProblem<3> problem(par); problem.run(); @@ -1834,7 +1843,7 @@ int main(int argc, char *argv[]) else { StokesImmersedProblemParameters<2> par; - ParameterAcceptor::initialize(prm_file, used_prm_file); + ParameterAcceptor::initialize(prm_file); StokesImmersedProblem<2> problem(par); problem.run();