Vortex()
: Function<dim>(dim)
{}
- virtual void vector_value(const Point<dim> &point,
- Vector<double> & values) const override;
+ virtual void
+ vector_value(const Point<dim> &point,
+ Vector<double> & values) const override;
};
template <int dim>
- void Vortex<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
public:
ParticleTracking(const ParticleTrackingParameters &par,
const bool interpolated_velocity);
- void run();
+ void
+ run();
private:
// The particles_generation function is responsible for the initial
// generation of the particles on top of the background grid
- void particles_generation();
+ void
+ particles_generation();
// When the velocity profile is interpolated to the position of the
// particles, it must first be stored using degrees of freedom.
// Consequently, as is the case for other parallel case (e.g. step-40) we
// initialize the degrees of freedom on the background grid
- void setup_background_dofs();
+ void
+ setup_background_dofs();
- void interpolate_function_to_field();
+ void
+ interpolate_function_to_field();
// The next two functions are responsible for carrying out explicit Euler
// time integration for the cases where the velocity field is interpolated
// at the positions of the particles or calculated analytically,
// respectively
- void euler_interpolated(double dt);
- void euler_analytical(double dt);
+ void
+ euler_interpolated(double dt);
+ void
+ euler_analytical(double dt);
// The cell_weight() function indicates to the triangulation how much
// computational work is expected to happen on this cell, and consequently
// roughly equal amount of work (potentially not an equal number of cells).
// While the function is called from the outside, it is connected to the
// corresponding signal from inside this class, therefore it can be private.
- unsigned int cell_weight(
+ unsigned int
+ cell_weight(
const typename parallel::distributed::Triangulation<dim>::cell_iterator
&cell,
const typename parallel::distributed::Triangulation<dim>::CellStatus
// The following two functions are responsible for outputting the simulation
// results for the particles and for the velocity profile on the background
// mesh, respectively.
- void output_particles(unsigned int it);
- void output_background(unsigned int it);
+ void
+ output_particles(unsigned int it);
+ void
+ output_background(unsigned int it);
// The private members of this class are similar to other parallel deal.II
// examples. The parameters are stored as a const member. It is important
// between ranks (the connection is created inside the
// particles_generation() function of this class).
template <int dim>
- unsigned int ParticleTracking<dim>::cell_weight(
+ unsigned int
+ ParticleTracking<dim>::cell_weight(
const typename parallel::distributed::Triangulation<dim>::cell_iterator
& cell,
- const typename parallel::distributed::Triangulation<dim>::CellStatus status) const
+ const typename parallel::distributed::Triangulation<dim>::CellStatus status)
+ const
{
// Assign no weight to cells we do not own.
if (!cell->is_locally_owned())
return 0;
// This determines how important particle work is compared to cell
- // work (by default every cell has a weight of 1000).
+ // work (by default every cell has a weight of 1000).
// We set the weight per particle much higher to indicate that
// the particle load is the only one that is important to distribute
// in this example. The optimal value of this number depends on the
// This function generates the tracer particles and the background
// triangulation on which these particles evolve.
template <int dim>
- void ParticleTracking<dim>::particles_generation()
+ void
+ ParticleTracking<dim>::particles_generation()
{
// We create an hyper_cube triangulation which we globally define. This
// triangulation englobes the full trajectory of the particles.
false));
// Establish the background triangulation where the particles are living
- particle_handler.initialize(background_triangulation, mapping);
+ // and the number of properties of the particles
+ particle_handler.initialize(background_triangulation, mapping, 1 + dim);
+ // pcout << "Number of properties "
+ // << particle_handler.n_properties_per_particle() << std::endl;
// We create a particle triangulation which is solely used to generate
// the points which will be used to insert the particles. This
// 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()
+ void
+ ParticleTracking<dim>::setup_background_dofs()
{
fluid_dh.distribute_dofs(fluid_fe);
IndexSet locally_owned_dofs = fluid_dh.locally_owned_dofs();
// Interpolates the Vortex velocity field to the field vector
template <int dim>
- void ParticleTracking<dim>::interpolate_function_to_field()
+ void
+ ParticleTracking<dim>::interpolate_function_to_field()
{
const MappingQ<dim> mapping(fluid_fe.degree);
// 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)
+ void
+ ParticleTracking<dim>::euler_analytical(double dt)
{
Vector<double> particle_velocity(dim);
particle_location[d] += particle_velocity[d] * dt;
particle->set_location(particle_location);
+
+ // Store the processor id and the particle velocity in the particle
+ // properties
+ ArrayView<double> properties = particle->get_properties();
+ properties[0] = Utilities::MPI::this_mpi_process(mpi_communicator);
+ for (int d = 0; d < dim; ++d)
+ properties[1 + d] += particle_velocity[d];
}
}
// 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)
+ 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);
for (int d = 0; d < dim; ++d)
particle_location[d] += particle_velocity[d] * dt;
particle->set_location(particle_location);
+
+ // Store the particle velocity and the processor id in the particle
+ // properties
+ ArrayView<double> properties = particle->get_properties();
+ properties[0] = Utilities::MPI::this_mpi_process(mpi_communicator);
+ for (int d = 0; d < dim; ++d)
+ properties[1 + d] += particle_velocity[d];
}
}
}
// and the background mesh to vtu with a pvtu record
template <int dim>
- void ParticleTracking<dim>::output_particles(unsigned int it)
+ void
+ ParticleTracking<dim>::output_particles(unsigned int it)
{
Particles::DataOut<dim, dim> particle_output;
particle_output.build_patches(particle_handler);
}
template <int dim>
- void ParticleTracking<dim>::output_background(unsigned int it)
+ void
+ ParticleTracking<dim>::output_background(unsigned int it)
{
std::vector<std::string> solution_names(dim, "velocity");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
// straightforward.
template <int dim>
- void ParticleTracking<dim>::run()
+ void
+ ParticleTracking<dim>::run()
{
DiscreteTime discrete_time(0, par.final_time, par.time_step);
pcout << "Repartitioning triangulation after particle generation"
<< std::endl;
- background_triangulation.repartition();
+ // background_triangulation.repartition();
setup_background_dofs();
interpolate_function_to_field();
// The remainder of the code, the `main()` function, is standard.
// We note that we run the particle tracking with the analytical velocity
// and the interpolated velocity and produce both results
-int main(int argc, char *argv[])
+int
+main(int argc, char *argv[])
{
using namespace Step68;
using namespace dealii;