From 0fe9ec48eaa61c84ec2d4964a6a99fe0aaa872bb Mon Sep 17 00:00:00 2001 From: Bruno Date: Fri, 18 Sep 2020 07:37:08 -0400 Subject: [PATCH] Finalized step-68 Fixed a bug related to deal.II distributed vector and to the particle loop in interpolated euler. --- examples/step-68/step-68.cc | 42 +++++++++++++--------------- include/deal.II/base/discrete_time.h | 6 ++-- 2 files changed, 22 insertions(+), 26 deletions(-) diff --git a/examples/step-68/step-68.cc b/examples/step-68/step-68.cc index 350349b821..9b951e8d96 100644 --- a/examples/step-68/step-68.cc +++ b/examples/step-68/step-68.cc @@ -280,8 +280,7 @@ namespace Step68 DoFHandler fluid_dh; FESystem fluid_fe; MappingQ mapping; - LinearAlgebra::distributed::Vector field_owned; - LinearAlgebra::distributed::Vector field_relevant; + LinearAlgebra::distributed::Vector velocity_field; Vortex velocity; @@ -416,15 +415,10 @@ namespace Step68 status) -> unsigned int { return this->cell_weight(cell, status); }); background_triangulation.signals.pre_distributed_repartition.connect( - std::bind( - &Particles::ParticleHandler::register_store_callback_function, - &particle_handler)); + [this]() { this->particle_handler.register_store_callback_function(); }); background_triangulation.signals.post_distributed_repartition.connect( - std::bind( - &Particles::ParticleHandler::register_load_callback_function, - &particle_handler, - false)); + [&]() { this->particle_handler.register_load_callback_function(false); }); // This initializes the background triangulation where the particles are // living and the number of properties of the particles. @@ -495,8 +489,7 @@ namespace Step68 IndexSet locally_relevant_dofs; DoFTools::extract_locally_relevant_dofs(fluid_dh, locally_relevant_dofs); - field_owned.reinit(locally_owned_dofs, mpi_communicator); - field_relevant.reinit(locally_owned_dofs, + velocity_field.reinit(locally_owned_dofs, locally_relevant_dofs, mpi_communicator); } @@ -511,8 +504,9 @@ namespace Step68 { const MappingQ mapping(fluid_fe.degree); - VectorTools::interpolate(mapping, fluid_dh, velocity, field_owned); - field_relevant = field_owned; + velocity_field.zero_out_ghosts(); + VectorTools::interpolate(mapping, fluid_dh, velocity, velocity_field); + velocity_field.update_ghost_values(); } @@ -569,7 +563,7 @@ namespace Step68 // Rather, 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 + // out of the `velocity_field` vector once and use them for all particles // that lie within the cell. auto particle = particle_handler.begin(); while (particle != particle_handler.end()) @@ -579,7 +573,7 @@ namespace Step68 const auto dh_cell = typename DoFHandler::cell_iterator(*cell, &fluid_dh); - dh_cell->get_dof_values(field_relevant, local_dof_values); + dh_cell->get_dof_values(velocity_field, local_dof_values); // Next, compute the velocity at the particle locations by evaluating // the finite element solution at the position of the particles. @@ -589,11 +583,12 @@ namespace Step68 // evaluation by hand, which is somewhat more efficient and only // matters for this tutorial, because the particle work is the // dominant cost of the whole program. - while (particle->get_surrounding_cell(background_triangulation) == cell) + const auto pic = particle_handler.particles_in_cell(cell); + Assert(pic.begin() == particle, ExcInternalError()); + for (auto &p : pic) { - const Point reference_location = - particle->get_reference_location(); - Tensor<1, dim> particle_velocity; + const Point reference_location = p.get_reference_location(); + Tensor<1, dim> particle_velocity; for (unsigned int j = 0; j < fluid_fe.dofs_per_cell; ++j) { const auto comp_j = fluid_fe.system_to_component_index(j); @@ -606,11 +601,11 @@ namespace Step68 Point particle_location = particle->get_location(); for (int d = 0; d < dim; ++d) particle_location[d] += particle_velocity[d] * dt; - particle->set_location(particle_location); + p.set_location(particle_location); // Again, we store the particle velocity and the processor id in the // particle properties for visualization purposes. - ArrayView properties = particle->get_properties(); + ArrayView properties = p.get_properties(); for (int d = 0; d < dim; ++d) properties[d] = particle_velocity[d]; @@ -675,7 +670,7 @@ namespace Step68 // Attach the solution data to data_out object data_out.attach_dof_handler(fluid_dh); - data_out.add_data_vector(field_relevant, + data_out.add_data_vector(velocity_field, solution_names, DataOut::type_dof_data, data_component_interpretation); @@ -732,7 +727,8 @@ namespace Step68 euler_step_analytical(0.); output_particles(discrete_time.get_step_number()); - output_background(discrete_time.get_step_number()); + if (interpolated_velocity) + output_background(discrete_time.get_step_number()); // The particles are advected by looping over time. while (!discrete_time.is_at_end()) diff --git a/include/deal.II/base/discrete_time.h b/include/deal.II/base/discrete_time.h index 84341d6f82..51c81f85f2 100644 --- a/include/deal.II/base/discrete_time.h +++ b/include/deal.II/base/discrete_time.h @@ -26,9 +26,9 @@ DEAL_II_NAMESPACE_OPEN * $T_{\text{start}}$ to an end time $T_{\text{end}}$. It also allows adjusting * the time step size during the simulation. This class provides the necessary * interface to be incorporated in any time-dependent simulation. As an - * example, the usage of this class is demonstrated in step-21 and step-68. - * This class attempts to replace the usage of TimestepControl with a better - * and more modern interface. + * example, the usage of this class is demonstrated in step-21. This class + * attempts to replace the usage of TimestepControl with a better and more + * modern interface. * * This class provides a number of invariants that are guaranteed to be * true at all times. -- 2.39.5