<h1>Results</h1>
-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
@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
# 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
# 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
# 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.
# 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
# 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.
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
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
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
@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:
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`.
<h3> Two dimensional test case </h3>
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.
</p>
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.
<p align="center">
</div>
</p>
-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.
<h3> Three dimensional test case </h3>
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
* ---------------------------------------------------------------------
*
- * Authors: Luca Heltai, Bruno Blais, 2020
+ * Authors: Luca Heltai, Bruno Blais, Rene Gassmoeller, 2020
*/
// @sect3{Include files}
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
};
// 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 <int spacedim>
class SolidPosition : public Function<spacedim>
{
, 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<spacedim> &p,
unsigned int component = 0) const override
{
- Tensor<1, spacedim> displacement = p;
+ Point<spacedim> 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)
}
#else
(void)ids_and_cad_file_names;
+ AssertThrow(false, ExcNotImplemented("Generation of the grid failed."));
#endif
}
// 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.
}
// 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
// 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 <int dim, int spacedim>
void StokesImmersedProblem<dim, spacedim>::initial_setup()
{
std::make_unique<QGauss<spacedim>>(par.velocity_degree + 1);
solid_quadrature_formula =
std::make_unique<QGauss<dim>>(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);
}
// 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())
{
}
- // 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 <int dim, int spacedim>
DataOut<spacedim>::type_dof_data,
data_component_interpretation);
- LA::MPI::BlockVector interpolated;
- interpolated.reinit(fluid_owned_dofs, MPI_COMM_WORLD);
- VectorTools::interpolate(fluid_dh,
- ConstantFunction<spacedim>(1.0, spacedim + 1),
- interpolated);
-
- LA::MPI::BlockVector interpolated_relevant(fluid_owned_dofs,
- fluid_relevant_dofs,
- MPI_COMM_WORLD);
- interpolated_relevant = interpolated;
- {
- std::vector<std::string> solution_names(spacedim, "ref_u");
- solution_names.emplace_back("ref_p");
- data_out.add_data_vector(interpolated_relevant,
- solution_names,
- DataOut<spacedim>::type_dof_data,
- data_component_interpretation);
- }
-
Vector<float> subdomain(fluid_tria.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i)
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<std::pair<double, std::string>> 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);
}
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<std::string, std::vector<std::pair<double, std::string>>>
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 <int dim, int spacedim>
void StokesImmersedProblem<dim, spacedim>::run()
{
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)
{
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
// 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");
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);
"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.
} // 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;
// 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();
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();
else
{
StokesImmersedProblemParameters<2> par;
- ParameterAcceptor::initialize(prm_file, used_prm_file);
+ ParameterAcceptor::initialize(prm_file);
StokesImmersedProblem<2> problem(par);
problem.run();
projects / dealii.git / commitdiff