--- /dev/null
+<br>
+
+<i>
+
+Bruno Blais (Polytechnique Montréal),
+Toni El Geitani Nehme (Polytechnique Montreal),
+Rene Gassmöller (University of California Davis),
+and Peter Munch
+</i>
+
+@dealiiTutorialDOI{10.5281/zenodo.3829064,https://zenodo.org/badge/DOI/10.5281/zenodo.3829064.svg}
+
+
+<h1>Introduction</h1>
+
+<h3>Simulation of the motion of massless tracer particles in a vortical flow</h3>
+
+Particles play an important part in numerical models for a large
+ number of applications. Particles are routinely used
+ as massless tracer to visualize the dynamic of a transient flow. They
+ can also play an intrinsic role as part of a more complex finite element
+ model, as is the case of the Particle-In-Cell (PIC) method (Gassmöller et al. 2018)
+ or they can even be used to simulate the motion of granular matter, as is
+ the case with the Discrete Element Method (DEM) (Blais et al. 2019). In the case
+ of DEM, the resulting model is not related to the finite element method anymore,
+ but just leads to a serie of ordinary differential equation which describes
+ the motion of the particles and the dynamic of their collisions. All of
+ these models can be built using deal.II particle handling capabilities.
+
+In the present step, we use particles as massless tracer to illustrate
+the dynamic of a vortical flows. Since the particles are massless tracers,
+the position of each particle $i$ is described by the
+following ordinary differential equation (ODE):
+@f[
+\frac{d \textbf{x}_i}{dt} =\textbf{u}(\textbf{x}_i)
+@f]
+
+where $\textbf{x}_i$ is the postion of particle $i$. In the present step,
+this ODE is solved using the explicit Euler method. The resulting scheme is:
+@f[
+\textbf{x}_{i}^{t+\Delta t} = \textbf{x}_{i}^{t} + \Delta t \; \textbf{u}(\textbf{x}_{i}^{t})
+@f]
+
+where $\textbf{x}_{i}^{t+\Delta t}$ and $\textbf{x}_{i}^{t}$ are the position
+of particle $i$ at time $t+\Delta t$ and $t$, respectively and where $\Delta t$
+is the time step. In the present step, the velocity at the location of particles
+is obtainedin two different fashions:
+- By evaluating the velocity function at the location of the particles
+- By evaluating the velocity function on a background triangulation and, using
+a finite element support, interpolating at the position of the particle.
+
+The first approach is not really practical, since in general, the velocity profile
+is not known analytically. The second approach, based on interpolating a solution
+at the position of the particles, mimics exactly what would be done in a
+realistic computational fluid dynamic simulation. In this step, we illustrate both strategies.
+
+
+We note that much greater accuracy could be obtained by using a fourth
+order Runge-Kutta method or another appropriate scheme for the time integration
+of the motion of the particles. However, this does not alter which capacities
+are displayed in the present step.
+
+<h3>Particles in deal.II</h3>
+
+In deal.II, Particles::Particle are very simple and flexible entities that can be used
+to build PIC, DEM or any type of particle-based models. Particles have a location
+in real space, a location in the reference space of the element in which they
+lie and a unique ID. In the majority of cases, simulations that include
+particles require a significant number of them. Thus, it becomes interesting
+to handle all particles through an entity which agglomerates all particles.
+In deal.II, this is achieved through the use of the Particles::ParticleHandler class.
+
+By default, particles do not have a diameter,
+a mass or any other physical properties which we would generally expect of physical particles. Howevever, through
+a ParticleHandler, particles have access to a Particles::PropertyPool. This PropertyPool is
+an array which can be used to store any arbitrary number of properties
+associated with the particles. Consequently, users can build their own
+particle solver and attribute the desired properties to the particles (e.g. mass,
+ diameter, temperature, etc.). In the present tutorial, this is used to
+ store the value of the fluid velocity and the processor id to which the particles
+ belong.
+
+
+<h3>The testcase</h3>
+
+In the present step, we use particles as massless tracer to illustrate
+the dynamic of a vortical flow : the Rayleigh-Kotte Vortex. This flow pattern
+is generally used as a complex test case for interface tracking methods (e.g.
+ volume-of-fluid and level set approches) since
+it leads to strong rotation and elongation of the fluid (Blais 2014).
+
+The
+stream function $\Psi$ of this Rayleigh-Kotte vortex is defined as:
+
+@f[
+\Psi = \frac{1}{\pi} sin^2 (\pi x) \sin^2 (\pi y) \cos \left( \pi \frac{t}{T} \right)
+@f]
+where $T$ is the period of the flow. The velocity profile in 2D ($\textbf{u}=[u,v]^T$) is :
+@f{eqnarray*}
+ u &=& - \frac{\partial\Psi}{\partial y} = \frac{1}{\pi} \sin^2 (\pi x) \sin (\pi y) \cos (\pi y) \cos \left( \pi \frac{t}{T} \right)\\
+ v &=& \frac{\partial\Psi}{\partial x} = \frac{1}{\pi} \cos(\pi x) \sin(\pi x) \sin^2 (\pi y) \cos \left( \pi \frac{t}{T} \right)
+@f}
+
+
+
+
+
+<h3>References</h3>
+
+<ul>
+<li> Blais, Bruno, et al. (2019) "Experimental Methods in Chemical Engineering: Discrete
+ Element Method—DEM." The Canadian Journal of Chemical Engineering 97.7 : 1964-1973.
+
+<li>Gassmöller, Rene, et al. (2018). "Flexible and Scalable Particle‐in‐Cell Methods With
+ Adaptive Mesh Refinement for Geodynamic Computations." Geochemistry, Geophysics,
+ Geosystems 19.9 : 3596-3604.
+
+ <li>Blais, Bruno, et al. (2013) "Dealing with more than two materials in the FVCF–ENIP method."
+ European Journal of Mechanics-B/Fluids 42 1-9.
+ </ul>
--- /dev/null
+<h1>Results</h1>
+
+The directory in which this program is run contains a example parameter file by defualt.
+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.
+
+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
+can be found:
+
+@code
+----------------------------------------------------
+Exception on processing:
+
+--------------------------------------------------------
+An error occurred in line <74> of file <../source/base/parameter_acceptor.cc> in function
+ static void dealii::ParameterAcceptor::initialize(const std::string &, const std::string &, const ParameterHandler::OutputStyle, dealii::ParameterHandler &)
+The violated condition was:
+ false
+Additional information:
+ You specified <parameters.prm> as input parameter file, but it does not exist. We created it for you.
+--------------------------------------------------------
+
+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).
+
+On any number of core, the simulation output will look like:
+
+@code
+bash$ mpirun -np 4 ./step-68 parameters.prm
+Number of particles inserted: 606
+Repartitioning triangulation after particle generation
+Writing particle output file: analytical-particles-0
+Writing background field file: background-0
+Writing particle output file: analytical-particles-10
+Writing background field file: background-10
+Writing particle output file: analytical-particles-20
+Writing background field file: background-20
+Writing particle output file: analytical-particles-30
+Writing background field file: background-30
+...
+Number of particles inserted: 606
+Repartitioning triangulation after particle generation
+Writing particle output file: analytical-particles-0
+Writing background field file: background-0
+Writing particle output file: analytical-particles-10
+Writing background field file: background-10
+Writing particle output file: analytical-particles-20
+Writing background field file: background-20
+Writing particle output file: analytical-particles-30
+Writing background field file: background-30
+...
+Writing particle output file: interpolated-particles-1980
+Writing background field file: background-1980
+Writing particle output file: interpolated-particles-1990
+Writing background field file: background-1990
+Writing particle output file: interpolated-particles-2000
+Writing background field file: background-2000
+@endcode
+
+We notice that, by default, the simulation runs the particle tracking with
+an analytical velocity for 2000 iterations, then runs the particle tracking with
+velocity interpolationn for the same duration. The results are written every
+10 iterations.
+
+<h3> Motion of the particles </h3>
+
+The following animation displays the trajectory of the particles as they
+are advected by the flow field. We see that after the complete duration of the
+flow, the particle go back to their initial configuration as is expected.
+
+@htmlonly
+<p align="center">
+ <iframe width="560" height="500" src="https://youtu.be/EbgS5Ch35Xs"
+ frameborder="0"
+ allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
+ allowfullscreen></iframe>
+ </p>
+@endhtmlonly
+
+<h3> Dynamic load balancing </h3>
+
+The following animation shows the impact of dynamic load balancing. We clearly
+see that the subdomains adapt themselves to balance the number of particles per
+subdomain. However, a perfect load balancing is not reached, in part due to
+the coarseness of the background mesh.
+
+@htmlonly
+<p align="center">
+ <iframe width="560" height="500" src="https://youtu.be/ubUcsR4ECj4"
+ frameborder="0"
+ allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
+ allowfullscreen></iframe>
+ </p>
+@endhtmlonly
+
+
+<h3>Possibilities for extensions</h3>
+
+This steps highlights some of the main capabilities of particle, notably their
+capacity to be use in distributed parallel simulations. However, this step could
+be exteded in numerous manners:
+- High-order time integration (for example using a Runge-Kutta 4 method) could be
+used to increase the accuracy and allow for an increased time-steps.
+- The full equation of motion (with inertia) could be solved for the particles. In
+this case the particles would need to have additional properties such as their mass,
+and their diameter.
+- Coupling to a flow solver. This step could be straightforwardly coupled to any parallel
+steps in which the Stokes or the Navier-Stokes equations are solved (e.g. step-57)
* ---------------------------------------------------------------------
*
- * Authors: Bruno Blais, Toni El Geitani Nehme, Rene Gassmoeller, Luca Heltai,
- Wolfgang Banghert 2020
+ * Authors: Bruno Blais, Toni El Geitani Nehme, Rene Gassmoeller, Peter Munch
*/
+
+// @sect3{Include files}
+
+// The majority of the include files are generic
+
#include <deal.II/base/bounding_box.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/discrete_time.h>
-#include <deal.II/base/index_set.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/parameter_acceptor.h>
-#include <deal.II/base/std_cxx14/memory.h>
#include <deal.II/base/timer.h>
#include <deal.II/distributed/cell_weights.h>
-#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_tools.h>
-#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/petsc_vector.h>
-#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/particles/data_out.h>
-#include <deal.II/particles/generators.h>
+// From the following include file we import the ParticleHandler class
+// that allows you to manage
+// a collection of particles (objects of type Particles::Particle), representing
+// a collection of points with some attached properties (e.g., an id) floating
+// on a parallel::distributed::Triangulation. The methods and classes in the
+// namespace Particles allows one to easily implement Particle-In-Cell methods
+// and particle tracing on distributed triangulations
#include <deal.II/particles/particle_handler.h>
-#include <deal.II/particles/utilities.h>
+// We import the particles generator
+// which allow us to insert the particles. In the present step, the particle
+// are globally inserted using a non-matching hyper-shell triangulation
+#include <deal.II/particles/generators.h>
+
+// Since the particles do not form a triangulation, they have their
+// own specific data out class which will enable us to write them
+// to commonly used parallel vtu format
+#include <deal.II/particles/data_out.h>
+
+
+// This step uses parallel vector to interpolate the velocity field
+// at the position of the particles. This step supports the use of both
+// Trilinos and PETSC distributed vectors
#define FORCE_USE_OF_TRILINOS
namespace LA
} // namespace LA
#include <cmath>
-#include <fstream>
#include <iostream>
-#include <memory>
-namespace Stepx
+namespace Step68
{
using namespace dealii;
+ // @sect3{Run-time parameter handling}
+
+ // Similarly to what is done in step-60, we set up a class that holds
+ // all the parameters of our problem and derive it from the ParameterAcceptor
+ // class to simplify the management and creation of parameter files.
+ //
+ // The ParameterAcceptor paradigm requires all parameters to be writable by
+ // the ParameterAcceptor methods. In order to avoid bugs that would be very
+ // difficult to track down (such as writing things like `time = 0` instead of
+ // `time == 0`), we declare all the parameters in an external class, which is
+ // initialized before the actual `ParticleTracking` class, and pass it to
+ // the main class as a `const` reference.
+ //
+ // The constructor of the class 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.
class ParticleTrackingParameters : public ParameterAcceptor
{
public:
// There remains the task of declaring what run-time parameters we can accept
- // in input files.
+ // in input files. Since we have a very limited number of parameters, all
+ // parameters are declared in the same category.
ParticleTrackingParameters::ParticleTrackingParameters()
: ParameterAcceptor("Particle Tracking Problem/")
{
"Refinement of the volumetric mesh used to insert the particles");
}
- // Creating the function for the velocity profile.
+
+ // @sect3{Velocity profile}
+
+ // The velocity profile is provided as a Function object. We provide the
+ // velocity profile. In the present step, this function is hard-coded within
+ // the example. However, it could have been easily made using a ParsedFunction
template <int dim>
- class SingleVortex : public Function<dim>
+ class Vortex : public Function<dim>
{
public:
- SingleVortex()
+ Vortex()
: Function<dim>(dim)
{}
virtual void vector_value(const Point<dim> &point,
};
template <int dim>
- void SingleVortex<dim>::vector_value(const Point<dim> &point,
- Vector<double> & values) const
+ void Vortex<dim>::vector_value(const Point<dim> &point,
+ Vector<double> & values) const
{
const double T = 4;
+ // Since the velocity profile is time dependant, the present time in the
+ // simulation must be gathered from the Function object.
const double t = this->get_time();
const double px = numbers::PI * point(0);
}
}
- // Solver
+ // @sect3{The <code>PatricleTracking</code> class declaration}
+
+ // We are now ready to introduce the main class of our tutorial program.
+ // Contrarily to some other steps, there is an additional function that is
+ // left public other than the constructor and the `run()` method, which is the
+ // cell_weight function. This function is connected to the triangulation and
+ // must be callable from outside of the scope of this class. Everything else
+ // is left `private`, and accessed through the run method itself.
template <int dim>
class ParticleTracking
{
public:
ParticleTracking(const ParticleTrackingParameters &par,
const bool interpolated_velocity);
- void run_analytical_velocity();
+ void run();
+
+ // Rene you would be more proefficient than me to write this
unsigned int cell_weight(
const typename parallel::distributed::Triangulation<dim>::cell_iterator
&cell,
status);
private:
+ // The particles_generation function is responsible for the initial
+ // generation of the particles on top of the background grid
void particles_generation();
+
+ // When the velocity profile is interpolated to the position of the
+ // particles, it must first be stored using degrees of freedom.
+ // Consequently, as is the case for other parallel case (e.g. step-40) we
+ // initialize the degrees of freedom on the background grid
void setup_background_dofs();
+
+
void interpolate_function_to_field();
+
+ // The next two functions are responsible for carrying out explicit Euler
+ // time integration for the cases where the velocity field is interpolated
+ // at the positions of the particles or calculated analytically,
+ // respectively
void euler_interpolated(double dt);
void euler_analytical(double dt);
- void field_euler(double t, double dt, double T);
+
+ // The following two functions are responsible for outputting the simulation
+ // results for the particles and for the velocity profile on the background
+ // mesh, respectively.
void output_particles(unsigned int it);
void output_background(unsigned int it);
+ // The private member of this class are similar to other parallel deal.II
+ // examples. The parameters are stored as a const member. It is important
+ // to note that we keep the Vortex class as a member since its time
+ // must be modified as the simulation proceeds.
+
const ParticleTrackingParameters ∥
LA::MPI::Vector field_owned;
LA::MPI::Vector field_relevant;
- SingleVortex<dim> velocity;
+ Vortex<dim> velocity;
ConditionalOStream pcout;
- bool interpolated_velocity = false;
+ bool interpolated_velocity;
};
+ // @sect3{The <code>PatricleTracking</code> class implementation}
+
+ // @sect4{Constructor}
+
+ // Constructors and destructors are rather trivial. They are very similar
+ // to what is done in step-40. we set the set of processors we want to work on
+ // to all machines available (MPI_COMM_WORLD) and
+ // initialize the <code>pcout</code> variable to only allow processor zero
+ // to output anything to the standard output.
+
template <int dim>
ParticleTracking<dim>::ParticleTracking(const ParticleTrackingParameters &par,
const bool interpolated_velocity)
: par(par)
, mpi_communicator(MPI_COMM_WORLD)
- , background_triangulation(MPI_COMM_WORLD)
+ , background_triangulation(mpi_communicator)
, fluid_dh(background_triangulation)
, fluid_fe(FE_Q<dim>(par.velocity_degree), dim)
, mapping(par.velocity_degree)
- , pcout({std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0})
+ , pcout(
+ {std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0})
, interpolated_velocity(interpolated_velocity)
{}
+ // @sect4{Cell weight}
+
+ // To be completed by Rene
+ // This function is the key component that allow us to do dynamic load
+ // balancing. It attributes a weight to every cell that depends on the number
+ // of particles that lie within that cell.
+
template <int dim>
unsigned int ParticleTracking<dim>::cell_weight(
const typename parallel::distributed::Triangulation<dim>::cell_iterator
return 0;
}
- // Generation of particles using the grid where particles are generated at the
- // locations of the degrees of freedom.
+ // @sect4{Particles generation}
+
+ // This function generates the tracer particles and the background
+ // triangulation on which these particles evolve.
template <int dim>
void ParticleTracking<dim>::particles_generation()
{
- // Create a square triangulation
+ // We create an hyper_cube triangulation which we globally define. This
+ // triangulation englobes the full trajectory of the particles.
GridGenerator::hyper_cube(background_triangulation, 0, 1);
background_triangulation.refine_global(par.fluid_refinement);
&particle_handler,
false));
- // Establish where the particles are living
+ // Establish the background triangulation where the particles are living
particle_handler.initialize(background_triangulation, mapping);
- // Generate the necessary bounding boxes for the generator of the particles
- const auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
- background_triangulation, IteratorFilters::LocallyOwnedCell());
- const auto global_bounding_boxes =
- Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box);
+ // We create a particle triangulation which is solely used to generate
+ // the points which will be used to insert the particles. This
+ // triangulation is an hyper_shell which is off-set from the
+ // center of the simulation domain.
Point<dim> center;
center[0] = 0.5;
}
const double outer_radius = 0.15;
- const double inner_radius = 0.001;
+ const double inner_radius = 0.01;
- // Generation and refinement of the grid where the particles will be
- // created.
parallel::distributed::Triangulation<dim> particle_triangulation(
MPI_COMM_WORLD);
particles_dof_handler.distribute_dofs(particles_fe);
- // Generation of the particles using the Particles::Generators
+ // We generate the necessary bounding boxes for the particles generator
+ // These bounding boxes are required to quickly identify in which
+ // processors and which cell the inserted particle lies.
+ const auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
+ background_triangulation, IteratorFilters::LocallyOwnedCell());
+ const auto global_bounding_boxes =
+ Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box);
+
+ // We generate the particles at the position of the degree of
+ // freedom of the dummy particle triangulation
Particles::Generators::dof_support_points(particles_dof_handler,
global_bounding_boxes,
particle_handler);
<< particle_handler.n_global_particles() << std::endl;
}
- // Sets up the background degree of freedom using their interpolation
- // And allocated a vector where you can store the entire solution
- // of the velocity field
+ // @sect4{Background DOFs and interpolation}
+
+
+ // Sets up the background degree of freedom used for the velocity
+ // interpolation And allocate the field vector where the entire
+ // solution of the velocity field is stored
template <int dim>
void ParticleTracking<dim>::setup_background_dofs()
{
field_relevant.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
-
- pcout << "Number of degrees of freedom in background grid: "
- << fluid_dh.n_dofs() << std::endl;
}
+ // Interpolates the Vortex velocity field to the field vector
template <int dim>
void ParticleTracking<dim>::interpolate_function_to_field()
{
field_relevant = field_owned;
}
+ // @sect4{Time integration of the trajectories}
+
+ // We integrate the particle trajectories
+ // using an analytically defined velocity field. This is a relatively trivial
+ // usage of the particles.
+ template <int dim>
+ void ParticleTracking<dim>::euler_analytical(double dt)
+ {
+ Vector<double> particle_velocity(dim);
+
+ // Looping over all particles in the domain using a particle iterator
+ for (auto particle = particle_handler.begin();
+ particle != particle_handler.end();
+ ++particle)
+ {
+ // Get the velocity using the current location of particle
+ Point<dim> particle_location = particle->get_location();
+ velocity.vector_value(particle_location, particle_velocity);
+
+ // Updating the position of the particles and Setting the old position
+ // equal to the new position of the particle
+ for (int d = 0; d < dim; ++d)
+ particle_location[d] += particle_velocity[d] * dt;
+
+ particle->set_location(particle_location);
+ }
+ }
+
+
+ // We integrate the particle trajectories by interpolating the value of the
+ // velocity field at the degrees of freedom to the position of the particles.
template <int dim>
void ParticleTracking<dim>::euler_interpolated(double dt)
{
std::vector<types::global_dof_index> dof_indices(fluid_fe.dofs_per_cell);
-
Vector<double> dof_data_per_cell(fluid_fe.dofs_per_cell);
-
Tensor<1, dim> particle_velocity;
+ // We loop over all the local particles. Although this could be achieved
+ // directly by looping over all the cells, this would force us
+ // to loop over numerous cells which do not contain particles.
+ // Consequently, we loop over all the particles, but, we get the reference
+ // of the cell in which the particle lies and then loop over all particles
+ // within that cell. This enables us to gather the values of the velocity
+ // out of the field_relevant vector once and use them for all particles
+ // that lie within the cell. Once we are done with all particles on one
+ // cell, we advance the `particle` iterator to the particle past the end of
+ // the ones on the current cell (this is the last line of the `while` loop's
+ // body).
auto particle = particle_handler.begin();
while (particle != particle_handler.end())
{
typename DoFHandler<dim>::cell_iterator(*cell, &fluid_dh);
dh_cell->get_dof_indices(dof_indices);
- // Gather the DOF information in a local vector to prevent dynamically
- // re-accessing everything when there are multiple particles in a cell
+ // Gather the velocity information in a local vector to prevent
+ // dynamically re-accessing everything when there are multiple particles
+ // in a cell
for (unsigned int j = 0; j < fluid_fe.dofs_per_cell; ++j)
{
dof_data_per_cell[j] = field_relevant(dof_indices[j]);
}
}
- template <int dim>
- void ParticleTracking<dim>::euler_analytical(double dt)
- {
- Vector<double> particle_velocity(dim);
- // Looping over all particles in the domain using a particle iterator
- for (auto particle = particle_handler.begin();
- particle != particle_handler.end();
- ++particle)
- {
- // Get the velocity using the current location of particle
- velocity.vector_value(particle->get_location(), particle_velocity);
+ // @sect4{Data output}
- Point<dim> particle_location = particle->get_location();
- // Updating the position of the particles and Setting the old position
- // equal to the new position of the particle
- for (int d = 0; d < dim; ++d)
- particle_location[d] += particle_velocity[d] * dt;
-
- particle->set_location(particle_location);
- }
- }
-
- // template <int dim>
- // void
- // ParticleTracking<dim>::parallel_weight()
- // {
- // parallel::CellWeights<dim> cell_weights(background_dh);
-
- // }
+ // These two functions take care of writing both the particles
+ // and the background mesh to vtu with a pvtu record
template <int dim>
void ParticleTracking<dim>::output_particles(unsigned int it)
output_folder, file_name, it, mpi_communicator, 6);
}
+ // @sect4{Running the simulation}
+ // This function orchestrates the entire simulation. It is very similar
+ // to the other time dependent tutorial programs -- take step-21 or step-26 as
+ // an example. Note that we use the DiscreteTime class to monitor the time,
+ // the time-step and the step-number. This function is relatively
+ // straightforward.
+
template <int dim>
- void ParticleTracking<dim>::run_analytical_velocity()
+ void ParticleTracking<dim>::run()
{
DiscreteTime discrete_time(0, par.final_time, par.time_step);
if ((discrete_time.get_step_number() % par.repartition_frequency) == 0)
{
- pcout << "Repartitioning triangulation after particle advection"
- << std::endl;
background_triangulation.repartition();
setup_background_dofs();
}
}
}
-} // namespace Stepx
+} // namespace Step68
// @sect3{The main() function}
// The remainder of the code, the `main()` function, is standard.
+// We note that we run the particle tracking with the analytical velocity
+// and the interpolated velocity and produce both results
int main(int argc, char *argv[])
{
- using namespace Stepx;
+ using namespace Step68;
using namespace dealii;
deallog.depth_console(1);
ParticleTrackingParameters par;
ParameterAcceptor::initialize(prm_file);
{
- Stepx::ParticleTracking<2> particle_tracking(par, false);
- particle_tracking.run_analytical_velocity();
+ Step68::ParticleTracking<2> particle_tracking(par, false);
+ particle_tracking.run();
}
{
- Stepx::ParticleTracking<2> particle_tracking(par, true);
- particle_tracking.run_analytical_velocity();
+ Step68::ParticleTracking<2> particle_tracking(par, true);
+ particle_tracking.run();
}
}
catch (std::exception &exc)