On Wednesday, August 10, 2022 at 1:06:58 AM UTC+2 Wolfgang Bangerth wrote:
>
> It's hard to say what the reason is without having a small test case that
> shows the issue. Since the issue happens during solution transfer, it is
> almost certainly not something related to whether or not you build
> solutions.
>
> As is generically the case, it is easiest to find the problem if you
> reduce
> the program that produces the error to something as minimal as possible --
> just throw out everything that can be thrown out. Do you think you can
> come up
> with such a program?
>
Thank you very much for your input! Your suggestion regarding the minimal
working example is much appreciated. I trimmed down a code snippet even
more and the current version is attached. This example is roughly 340 lines
long and the relevant part with refine_grid() start at line 204. I included
a dummy refinement criterion.
> Separately (and almost certainly unrelated): Did you mean to reinit() the
> 'solution' variable twice here?
> ```
> LinearAlgebra::distributed::Vector<Number> old_solution;
> old_solution.reinit(locally_owned_dofs,
> locally_relevant_dofs,
> MPI_COMM_WORLD);
> old_solution = solution;
> soltrans.prepare_for_coarsening_and_refinement(old_solution);
>
> solution.reinit(locally_owned_dofs, MPI_COMM_WORLD);
> soltrans.interpolate(solution);
> ```
>
For this section I just followed the documentation for the distributed
SolutionTransfer but I was also a little bit confused about this part, to
be honest. I am currently playing around with this part to see what the
different outcomes can be.
Best regards,
M. Bänsch
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/c92684d4-da96-4798-8825-640ee7ce3525n%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2020 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Martin Kronbichler, 2020
*/
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria.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/la_parallel_vector.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <deal.II/matrix_free/operators.h>
namespace Euler_DG
{
using namespace dealii;
constexpr unsigned int dimension = 2;
constexpr unsigned int n_global_refinements = 0;
constexpr unsigned int fe_degree = 1;
constexpr unsigned int n_q_points_1d = fe_degree + 2;
using Number = double;
template <int dim, int degree, int n_points_1d>
class EulerOperator
{
public:
static constexpr unsigned int n_quadrature_points_1d = n_points_1d;
EulerOperator(TimerOutput &timer_output);
void reinit(const Mapping<dim> & mapping,
const DoFHandler<dim> &dof_handler);
void
initialize_vector(LinearAlgebra::distributed::Vector<Number> &vector) const;
private:
MatrixFree<dim, Number> data;
TimerOutput &timer;
};
template <int dim, int degree, int n_points_1d>
EulerOperator<dim, degree, n_points_1d>::EulerOperator(TimerOutput &timer)
: timer(timer)
{}
template <int dim, int degree, int n_points_1d>
void EulerOperator<dim, degree, n_points_1d>::reinit(
const Mapping<dim> & mapping,
const DoFHandler<dim> &dof_handler)
{
const std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_handler};
const AffineConstraints<double> dummy;
const std::vector<const AffineConstraints<double> *> constraints = {&dummy};
const std::vector<Quadrature<1>> quadratures = {QGauss<1>(n_q_points_1d),
QGauss<1>(fe_degree + 1)};
typename MatrixFree<dim, Number>::AdditionalData additional_data;
additional_data.mapping_update_flags =
(update_gradients | update_JxW_values | update_quadrature_points |
update_values);
additional_data.mapping_update_flags_inner_faces =
(update_JxW_values | update_quadrature_points | update_normal_vectors |
update_values);
additional_data.mapping_update_flags_boundary_faces =
(update_JxW_values | update_quadrature_points | update_normal_vectors |
update_values);
additional_data.tasks_parallel_scheme =
MatrixFree<dim, Number>::AdditionalData::none;
data.reinit(
mapping, dof_handlers, constraints, quadratures, additional_data);
}
template <int dim, int degree, int n_points_1d>
void EulerOperator<dim, degree, n_points_1d>::initialize_vector(
LinearAlgebra::distributed::Vector<Number> &vector) const
{
data.initialize_dof_vector(vector);
}
template <int dim>
class EulerProblem
{
public:
EulerProblem();
void run();
private:
void make_grid_and_dofs();
void refine_grid();
void output_results(const unsigned int result_number);
LinearAlgebra::distributed::Vector<Number> solution;
ConditionalOStream pcout;
#ifdef DEAL_II_WITH_P4EST
parallel::distributed::Triangulation<dim> triangulation;
#else
Triangulation<dim> triangulation;
#endif
FESystem<dim> fe;
MappingQGeneric<dim> mapping;
DoFHandler<dim> dof_handler;
TimerOutput timer;
EulerOperator<dim, fe_degree, n_q_points_1d> euler_operator;
double time, time_step;
};
template <int dim>
EulerProblem<dim>::EulerProblem()
: pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
#ifdef DEAL_II_WITH_P4EST
, triangulation(MPI_COMM_WORLD)
#endif
, fe(FE_DGQ<dim>(fe_degree), dim + 2)
, mapping(fe_degree)
, dof_handler(triangulation)
, timer(pcout, TimerOutput::never, TimerOutput::wall_times)
, euler_operator(timer)
, time(0)
, time_step(0)
{}
template <int dim>
void EulerProblem<dim>::make_grid_and_dofs()
{
GridGenerator::channel_with_cylinder(
triangulation, 0.03, 1, 0, true);
triangulation.refine_global(n_global_refinements);
dof_handler.distribute_dofs(fe);
euler_operator.reinit(mapping, dof_handler);
euler_operator.initialize_vector(solution);
}
template <int dim>
void EulerProblem<dim>::refine_grid()
{
const int max_h_level = 5;
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell->clear_coarsen_flag();
cell->clear_refine_flag();
if ((cell->level() < max_h_level) &&
(cell->extent_in_direction(0) < 0.1))
cell->set_refine_flag();
}
// see distributed::SolutionTransfer documentation for details
parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<Number>> soltrans(dof_handler);
triangulation.prepare_coarsening_and_refinement();
soltrans.prepare_for_coarsening_and_refinement(solution);
triangulation.execute_coarsening_and_refinement();
dof_handler.distribute_dofs(fe);
IndexSet locally_owned_dofs, locally_relevant_dofs;
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
LinearAlgebra::distributed::Vector<Number> solution;
solution.reinit(locally_owned_dofs, MPI_COMM_WORLD);
LinearAlgebra::distributed::Vector<Number> old_solution;
old_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
MPI_COMM_WORLD);
old_solution = solution;
soltrans.prepare_for_coarsening_and_refinement(old_solution);
solution.reinit(locally_owned_dofs, MPI_COMM_WORLD);
soltrans.interpolate(solution);
}
template <int dim>
void EulerProblem<dim>::output_results(const unsigned int result_number)
{
{
TimerOutput::Scope t(timer, "output");
DataOut<dim> data_out;
DataOutBase::VtkFlags flags;
flags.write_higher_order_cells = true;
data_out.set_flags(flags);
data_out.attach_dof_handler(dof_handler);
Vector<double> mpi_owner(triangulation.n_active_cells());
mpi_owner = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
data_out.add_data_vector(mpi_owner, "owner");
data_out.build_patches(mapping,
fe.degree,
DataOut<dim>::curved_inner_cells);
const std::string filename =
"solution_" + Utilities::int_to_string(result_number, 3) + ".vtu";
data_out.write_vtu_in_parallel(filename, MPI_COMM_WORLD);
}
}
template <int dim>
void EulerProblem<dim>::run()
{
pcout << "Running with "
<< Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)
<< " MPI processes" << std::endl;
make_grid_and_dofs();
output_results(0);
refine_grid();
output_results(1);
}
} // namespace Euler_DG
int main(int argc, char **argv)
{
using namespace Euler_DG;
using namespace dealii;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
try
{
deallog.depth_console(0);
EulerProblem<dimension> euler_problem;
euler_problem.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
