I am trying to implement a test program using the matrix-free method,
combined with a geometric multigrid preconditioner. In my case my
vmult/apply_add-function needs access to the solution from the last step.
Thus I added a vector
LinearAlgebra::distributed::Vector<number> old_solution;
which is filled using
template <int dim, int fe_degree, typename number>
void LaplaceOperator<dim, fe_degree, number>::update_old_solution(const
LinearAlgebra::distributed::Vector<number> &solution){
this->old_solution =
LinearAlgebra::distributed::Vector<number>(solution);
}
For each level of grids I need to fill this vector with the current
solution, interpolated onto the correct mesh. Thus, I added the following
code:
const unsigned int nlevels = triangulation.n_global_levels();
mg_matrices.resize(0, nlevels - 1);
std::set<types::boundary_id> dirichlet_boundary;
dirichlet_boundary.insert(0);
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
dirichlet_boundary);
MGTransferMatrixFree<dim, float> mg_transfer(mg_constrained_dofs);
const unsigned int max_level =
dof_handler.get_triangulation().n_global_levels() - 1;
const unsigned int min_level = 0;
MGLevelObject<LinearAlgebra::distributed::Vector<float>>
level_projection(min_level, max_level);
pcout << "Interpolate to mg\n";
mg_transfer.interpolate_to_mg(dof_handler, level_projection, solution);
pcout << "Interpolated to mg\n";
for (unsigned int level = 0; level < nlevels; ++level)
{
pcout << "Level: " << level << '\n';
IndexSet relevant_dofs;
DoFTools::extract_locally_relevant_level_dofs(dof_handler,
level,
relevant_dofs);
AffineConstraints<double> level_constraints;
level_constraints.reinit(relevant_dofs);
level_constraints.add_lines(
mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
typename MatrixFree<dim, float>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree<dim, float>::AdditionalData::partition_partition;
additional_data.mapping_update_flags =
(update_gradients | update_JxW_values |
update_quadrature_points);
additional_data.level_mg_handler = level;
std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(
new MatrixFree<dim, float>());
mg_mf_storage_level->reinit(dof_handler,
level_constraints,
QGauss<1>(fe.degree + 1),
additional_data);
mg_matrices[level].initialize(mg_mf_storage_level,
mg_constrained_dofs,
level);
mg_matrices[level].set_time_step(time_step);
mg_matrices[level].update_old_solution(level_projection[level]);
}
My program works without crash (even though I get wrong results, but that
is another problem) when compiling it in debug mode, but crashes
immediately after pcout << "Interpolate to mg\n" with a segfault; when
running in release mode. Why is that? Or are there other methods I could
use for interpolation?
Thanks!
--
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 dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/c8e39843-b9d1-4694-ac8c-bef018b42782%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2009 - 2019 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.
*
* ---------------------------------------------------------------------
*
* Authors: Katharina Kormann, Martin Kronbichler, Uppsala University,
* 2009-2012, updated to MPI version with parallel vectors in 2016
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/operators.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <iostream>
#include <fstream>
namespace Step37
{
using namespace dealii;
const unsigned int degree_finite_element = 3;
const unsigned int dimension = 2;
template <int dim>
class Coefficient : public Function<dim>
{
public:
Coefficient()
: Function<dim>()
{}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
template <typename number>
number value(const Point<dim, number> &p,
const unsigned int component = 0) const;
};
template <int dim>
template <typename number>
number Coefficient<dim>::value(const Point<dim, number> &p,
const unsigned int /*component*/) const
{
return 2 * M_PI * M_PI * sin(M_PI * p[0]) * sin(M_PI * p[1]);
}
template <int dim>
double Coefficient<dim>::value(const Point<dim> & p,
const unsigned int component) const
{
return value<double>(p, component);
}
template <int dim>
class Solution : public Function<dim>
{
public:
Solution(const double cur_time) : Function<dim>(1), cur_time(cur_time)
{
}
virtual double value(const Point<dim> &p, const unsigned int component) const override;
//virtual Tensor<1, dim> gradient(const Point<dim> &p, const unsigned int component) const override;
private:
const double cur_time;
};
template <int dim>
double Solution<dim>::value(const Point<dim> &p, const unsigned int /*component*/) const
{
const double x = p[0];
const double y = p[1];
return /*exp(-cur_time) * */sin(M_PI * x) * sin(M_PI * y);
}
//template<int dim>
//Tensor<1, dim> Solution<dim>::gradient(const Point<dim> &p, const unsigned int /*component*/) const
//{
// Tensor<1, dim> return_value;
// AssertThrow(dim == 2, ExcNotImplemented());
// const double x = p[0];
// const double y = p[1];
// return_value[0] = exp(-cur_time) * M_PI * cos(M_PI * x) * sin(M_PI * y);
// return_value[1] = exp(-cur_time) * M_PI * sin(M_PI * y) * cos(M_PI * x);
// return return_value;
//}
using vector_t = LinearAlgebra::distributed::Vector<double>;
template <int dim, typename number>
void adjust_this_ghost_range_if_necessary(
const MatrixFree<dim, number> & data,
const LinearAlgebra::distributed::Vector<number> &vec)
{
if (vec.get_partitioner().get() ==
data.get_dof_info(0).vector_partitioner.get())
return;
std::cout << "vec: " << vec.local_size() << '\n';
LinearAlgebra::distributed::Vector<number> copy_vec(vec);
std::cout << "copy_vec: " << copy_vec.local_size() << '\n';
const_cast<LinearAlgebra::distributed::Vector<number> &>(vec).reinit(data.get_dof_info(0).vector_partitioner);
std::cout << "vec: " << vec.local_size() << '\n';
const_cast<LinearAlgebra::distributed::Vector<number> &>(vec).copy_locally_owned_data_from(copy_vec);
}
template <int dim, int fe_degree, typename number>
class LaplaceOperator
: public MatrixFreeOperators::
Base<dim, LinearAlgebra::distributed::Vector<number>>
{
public:
using value_type = number;
LaplaceOperator();
void clear() override;
virtual void compute_diagonal() override;
void set_time_step(const double time_step){
this->time_step = time_step;
};
double time_step = 0;
void update_old_solution(const LinearAlgebra::distributed::Vector<number> &solution);
private:
virtual void apply_add(
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src) const override;
void
local_apply(const MatrixFree<dim, number> & data,
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;
void local_compute_diagonal(
const MatrixFree<dim, number> & data,
LinearAlgebra::distributed::Vector<number> & dst,
const unsigned int & dummy,
const std::pair<unsigned int, unsigned int> &cell_range) const;
Table<2, VectorizedArray<number>> coefficient;
LinearAlgebra::distributed::Vector<number> old_solution;
};
template <int dim, int fe_degree, typename number>
LaplaceOperator<dim, fe_degree, number>::LaplaceOperator()
: MatrixFreeOperators::Base<dim,
LinearAlgebra::distributed::Vector<number>>()
{}
template <int dim, int fe_degree, typename number>
void LaplaceOperator<dim, fe_degree, number>::clear()
{
coefficient.reinit(0, 0);
MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>::
clear();
}
template <int dim, int fe_degree, typename number>
void LaplaceOperator<dim, fe_degree, number>::local_apply(
const MatrixFree<dim, number> & data,
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> & cell_range) const
{
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(data), phi_old(data);
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
{
phi.reinit(cell);
phi_old.reinit(cell);
LinearAlgebra::distributed::Vector<number> local_solution(old_solution), local_src(src);
local_src *= 1e-6;
// std::cout << "Local_src: " << local_src.local_size() << '\n';
// std::cout << "Solution: " << local_solution.local_size() << '\n';
local_solution += local_src;
phi.read_dof_values(local_solution);
phi_old.read_dof_values(old_solution);
phi.evaluate(true, true);
phi_old.evaluate(true, true);
for (unsigned int q = 0; q < phi.n_q_points; ++q){
phi.submit_gradient(1e-6 * (phi.get_gradient(q) - phi_old.get_gradient(q)), q);
// //phi.submit_value(phi.get_value(q) + make_vectorized_array<number>(0.0) * phi_old.get_value(q), q);
// phi.submit_gradient(phi.get_gradient(q),
// q);
}
phi.integrate(false, true);
phi.distribute_local_to_global(dst);
}
}
template <int dim, int fe_degree, typename number>
void LaplaceOperator<dim, fe_degree, number>::apply_add(
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src) const
{
this->data->cell_loop(&LaplaceOperator::local_apply, this, dst, src);
}
template <int dim, int fe_degree, typename number>
void LaplaceOperator<dim, fe_degree, number>::compute_diagonal()
{
this->inverse_diagonal_entries.reset(
new DiagonalMatrix<LinearAlgebra::distributed::Vector<number>>());
LinearAlgebra::distributed::Vector<number> &inverse_diagonal =
this->inverse_diagonal_entries->get_vector();
this->data->initialize_dof_vector(inverse_diagonal);
unsigned int dummy = 0;
this->data->cell_loop(&LaplaceOperator::local_compute_diagonal,
this,
inverse_diagonal,
dummy);
this->set_constrained_entries_to_one(inverse_diagonal);
for (unsigned int i = 0; i < inverse_diagonal.local_size(); ++i)
{
Assert(inverse_diagonal.local_element(i) > 0.,
ExcMessage("No diagonal entry in a positive definite operator "
"should be zero"));
inverse_diagonal.local_element(i) =
1. / inverse_diagonal.local_element(i);
}
}
template <int dim, int fe_degree, typename number>
void LaplaceOperator<dim, fe_degree, number>::local_compute_diagonal(
const MatrixFree<dim, number> & data,
LinearAlgebra::distributed::Vector<number> &dst,
const unsigned int &,
const std::pair<unsigned int, unsigned int> &cell_range) const
{
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(data), phi_old(data);
AlignedVector<VectorizedArray<number>> diagonal(phi.dofs_per_cell);
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
{
phi.reinit(cell);
phi_old.reinit(cell);
for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < phi.dofs_per_cell; ++j){
phi.submit_dof_value(VectorizedArray<number>(), j);
phi_old.submit_dof_value(VectorizedArray<number>(), j);
}
phi.submit_dof_value(make_vectorized_array<number>(1. + 1e-6), i);
phi_old.submit_dof_value(make_vectorized_array<number>(1.), i);
phi.evaluate(false, true);
phi_old.evaluate(false, true);
for (unsigned int q = 0; q < phi.n_q_points; ++q){
phi.submit_gradient((phi.get_gradient(q) - phi_old.get_gradient(q)) * 1e-6,
q);
}
phi.integrate(false, true);
diagonal[i] = phi.get_dof_value(i);
}
for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
phi.submit_dof_value(diagonal[i], i);
phi.distribute_local_to_global(dst);
}
}
template <int dim, int fe_degree, typename number>
void LaplaceOperator<dim, fe_degree, number>::update_old_solution(const LinearAlgebra::distributed::Vector<number> &solution){
this->old_solution = LinearAlgebra::distributed::Vector<number>(solution);
std::cout << "Solution-size is " << this->old_solution.local_size() << '\n';
}
template <int dim>
class LaplaceProblem
{
public:
LaplaceProblem();
void run();
private:
void setup_system();
void assemble_rhs();
void solve();
void output_results(const unsigned int cycle);
void refine_mesh();
#ifdef DEAL_II_WITH_P4EST
parallel::distributed::Triangulation<dim> triangulation;
#else
Triangulation<dim> triangulation;
#endif
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> boundary_constraints,
newton_constraints;
using SystemMatrixType =
LaplaceOperator<dim, degree_finite_element, double>;
SystemMatrixType system_matrix;
MGConstrainedDoFs mg_constrained_dofs;
using LevelMatrixType = LaplaceOperator<dim, degree_finite_element, float>;
MGLevelObject<LevelMatrixType> mg_matrices;
LinearAlgebra::distributed::Vector<double> solution;
LinearAlgebra::distributed::Vector<double> system_rhs;
LinearAlgebra::distributed::Vector<double> newton_update;
double setup_time;
ConditionalOStream pcout;
ConditionalOStream time_details;
const double time_step = 0.001;
};
template <int dim>
LaplaceProblem<dim>::LaplaceProblem()
:
#ifdef DEAL_II_WITH_P4EST
triangulation(
MPI_COMM_WORLD,
Triangulation<dim>::limit_level_difference_at_vertices,
parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy)
,
#else
triangulation(Triangulation<dim>::limit_level_difference_at_vertices)
,
#endif
fe(degree_finite_element)
, dof_handler(triangulation)
, setup_time(0.)
, pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
,
time_details(std::cout,
false && Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
{}
template <int dim>
void LaplaceProblem<dim>::setup_system()
{
Timer time;
setup_time = 0;
system_matrix.clear();
mg_matrices.clear_elements();
dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
boundary_constraints.clear();
boundary_constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, boundary_constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
Solution<dim>(0.),
boundary_constraints);
boundary_constraints.close();
newton_constraints.clear();
newton_constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, newton_constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(),
newton_constraints);
newton_constraints.close();
setup_time += time.wall_time();
time_details << "Distribute DoFs & B.C. (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s" << std::endl;
time.restart();
{
typename MatrixFree<dim, double>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree<dim, double>::AdditionalData::partition_partition;
additional_data.mapping_update_flags =
(update_gradients | update_JxW_values | update_quadrature_points);
std::shared_ptr<MatrixFree<dim, double>> system_mf_storage(
new MatrixFree<dim, double>());
system_mf_storage->reinit(dof_handler,
boundary_constraints,
QGauss<1>(fe.degree + 1),
additional_data);
system_matrix.initialize(system_mf_storage);
}
system_matrix.set_time_step(time_step);
system_matrix.initialize_dof_vector(solution);
vector_t local_solution;
system_matrix.initialize_dof_vector(local_solution);
local_solution = solution;
system_matrix.update_old_solution(local_solution);
system_matrix.initialize_dof_vector(newton_update);
system_matrix.initialize_dof_vector(system_rhs);
setup_time += time.wall_time();
time_details << "Setup matrix-free system (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s" << std::endl;
time.restart();
const unsigned int nlevels = triangulation.n_global_levels();
mg_matrices.resize(0, nlevels - 1);
std::set<types::boundary_id> dirichlet_boundary;
dirichlet_boundary.insert(0);
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
dirichlet_boundary);
MGTransferMatrixFree<dim, float> mg_transfer(mg_constrained_dofs);
const unsigned int max_level =
dof_handler.get_triangulation().n_global_levels() - 1;
const unsigned int min_level = 0;
MGLevelObject<LinearAlgebra::distributed::Vector<float>>
level_projection(min_level, max_level);
pcout << "Interpolate to mg\n";
mg_transfer.interpolate_to_mg(dof_handler, level_projection, solution);
pcout << "Interpolated to mg\n";
for (unsigned int level = 0; level < nlevels; ++level)
{
pcout << "Level: " << level << '\n';
IndexSet relevant_dofs;
DoFTools::extract_locally_relevant_level_dofs(dof_handler,
level,
relevant_dofs);
AffineConstraints<double> level_constraints;
level_constraints.reinit(relevant_dofs);
level_constraints.add_lines(
mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
typename MatrixFree<dim, float>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree<dim, float>::AdditionalData::partition_partition;
additional_data.mapping_update_flags =
(update_gradients | update_JxW_values | update_quadrature_points);
additional_data.level_mg_handler = level;
std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(
new MatrixFree<dim, float>());
mg_mf_storage_level->reinit(dof_handler,
level_constraints,
QGauss<1>(fe.degree + 1),
additional_data);
mg_matrices[level].initialize(mg_mf_storage_level,
mg_constrained_dofs,
level);
mg_matrices[level].set_time_step(time_step);
// LinearAlgebra::distributed::Vector<float> local_solution;
// mg_matrices[level].initialize_dof_vector(local_solution);
//std::cout << "In level " << level << " the solution size is " << level_projection[level].local_size() << '\n';
//local_solution = solution;
mg_matrices[level].update_old_solution(level_projection[level]);
}
setup_time += time.wall_time();
time_details << "Setup matrix-free levels (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s" << std::endl;
}
template <int dim>
void LaplaceProblem<dim>::assemble_rhs()
{
Timer time;
system_rhs = 0;
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
LinearAlgebra::distributed::Vector<double> copy_vec(solution);
system_matrix.initialize_dof_vector(solution);
solution.copy_locally_owned_data_from(copy_vec);
boundary_constraints.distribute(solution);
solution.update_ghost_values();
FEEvaluation<dim, degree_finite_element> phi(
*system_matrix.get_matrix_free());
for (unsigned int cell = 0;
cell < system_matrix.get_matrix_free()->n_macro_cells();
++cell)
{
phi.reinit(cell);
phi.read_dof_values(solution);
phi.evaluate(true, true);
for (unsigned int q = 0; q < phi.n_q_points; ++q){
phi.submit_gradient(make_vectorized_array<double>(1.) * phi.get_gradient(q), q);
phi.submit_value(2. * M_PI * M_PI * sin(M_PI * phi.quadrature_point(q)[0]) * sin(M_PI * phi.quadrature_point(q)[1]), q);
}
phi.integrate(true, true);
phi.distribute_local_to_global(system_rhs);
}
system_rhs.compress(VectorOperation::add);
setup_time += time.wall_time();
time_details << "Assemble right hand side (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s" << std::endl;
}
template <int dim>
void LaplaceProblem<dim>::refine_mesh()
{
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
LinearAlgebra::distributed::Vector<double> copy_vec(solution);
solution.reinit(dof_handler.locally_owned_dofs(),
locally_relevant_dofs,
triangulation.get_communicator());
solution.copy_locally_owned_data_from(copy_vec);
boundary_constraints.distribute(solution);
solution.update_ghost_values();
// LinearAlgebra::distributed::Vector<double> intermediate_solution(dof_handler.locally_owned_dofs(),
// locally_relevant_dofs,
// MPI_COMM_WORLD);
// intermediate_solution = solution;
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler,
QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>> soltrans(dof_handler);
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
triangulation,
estimated_error_per_cell,
1, 0);
triangulation.prepare_coarsening_and_refinement();
soltrans.prepare_for_coarsening_and_refinement(solution);
triangulation.execute_coarsening_and_refinement();
setup_system();
// LinearAlgebra::distributed::Vector<double> interpolated_solution(dof_handler.locally_owned_dofs(),
// locally_relevant_dofs,
// MPI_COMM_WORLD);
soltrans.interpolate(solution);
//solution = interpolated_solution;
}
template <int dim>
void LaplaceProblem<dim>::solve()
{
Timer time;
MGTransferMatrixFree<dim, float> mg_transfer(mg_constrained_dofs);
mg_transfer.build(dof_handler);
setup_time += time.wall_time();
time_details << "MG build transfer time (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s\n";
time.restart();
using SmootherType =
PreconditionChebyshev<LevelMatrixType,
LinearAlgebra::distributed::Vector<float>>;
mg::SmootherRelaxation<SmootherType,
LinearAlgebra::distributed::Vector<float>>
mg_smoother;
MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
smoother_data.resize(0, triangulation.n_global_levels() - 1);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level)
{
if (level > 0)
{
smoother_data[level].smoothing_range = 15.;
smoother_data[level].degree = 5;
smoother_data[level].eig_cg_n_iterations = 10;
}
else
{
smoother_data[0].smoothing_range = 1e-3;
smoother_data[0].degree = numbers::invalid_unsigned_int;
smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
}
mg_matrices[level].compute_diagonal();
smoother_data[level].preconditioner =
mg_matrices[level].get_matrix_diagonal_inverse();
}
mg_smoother.initialize(mg_matrices, smoother_data);
MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<float>>
mg_coarse;
mg_coarse.initialize(mg_smoother);
mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(
mg_matrices);
MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<LevelMatrixType>>
mg_interface_matrices;
mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level)
mg_interface_matrices[level].initialize(mg_matrices[level]);
mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_interface(
mg_interface_matrices);
Multigrid<LinearAlgebra::distributed::Vector<float>> mg(
mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
mg.set_edge_matrices(mg_interface, mg_interface);
PreconditionMG<dim,
LinearAlgebra::distributed::Vector<float>,
MGTransferMatrixFree<dim, float>>
preconditioner(dof_handler, mg, mg_transfer);
SolverControl solver_control(100, 1e-12 * system_rhs.l2_norm());
SolverGMRES<LinearAlgebra::distributed::Vector<double>> cg(solver_control);
setup_time += time.wall_time();
time_details << "MG build smoother time (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s\n";
pcout << "Total setup time (wall) " << setup_time << "s\n";
time.reset();
time.start();
newton_update = 0.;
cg.solve(system_matrix, newton_update, system_rhs, preconditioner);
newton_constraints.distribute(newton_update);
// pcout << "Updating solution\n";
solution += newton_update;
// pcout << "Solution updated\n";
pcout << "Time solve (" << solver_control.last_step() << " iterations)"
<< (solver_control.last_step() < 10 ? " " : " ") << "(CPU/wall) "
<< time.cpu_time() << "s/" << time.wall_time() << "s\n";
}
template <int dim>
void LaplaceProblem<dim>::output_results(const unsigned int cycle)
{
Timer time;
if (triangulation.n_global_active_cells() > 1000000)
return;
DataOut<dim> data_out;
solution.update_ghost_values();
newton_update.update_ghost_values();
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.add_data_vector(newton_update, "newton_update");
data_out.build_patches();
std::ofstream output(
"solution-" + std::to_string(cycle) + "." +
std::to_string(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)) +
".vtu");
DataOutBase::VtkFlags flags;
flags.compression_level = DataOutBase::VtkFlags::best_speed;
data_out.set_flags(flags);
data_out.write_vtu(output);
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
{
std::vector<std::string> filenames;
for (unsigned int i = 0;
i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
++i)
filenames.emplace_back("solution-" + std::to_string(cycle) + "." +
std::to_string(i) + ".vtu");
std::string master_name =
"solution-" + Utilities::to_string(cycle) + ".pvtu";
std::ofstream master_output(master_name);
data_out.write_pvtu_record(master_output, filenames);
}
time_details << "Time write output (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s\n";
pcout << "Solution expected to be " << exp(-(this->time_step) * cycle) << " at time " << (this->time_step) * cycle << '\n';
}
template <int dim>
void LaplaceProblem<dim>::run()
{
{
const unsigned int n_vect_doubles =
VectorizedArray<double>::n_array_elements;
const unsigned int n_vect_bits = 8 * sizeof(double) * n_vect_doubles;
pcout << "Vectorization over " << n_vect_doubles
<< " doubles = " << n_vect_bits << " bits ("
<< Utilities::System::get_current_vectorization_level()
<< "), VECTORIZATION_LEVEL=" << DEAL_II_COMPILER_VECTORIZATION_LEVEL
<< std::endl;
}
for (unsigned int cycle = 0; cycle < 20; ++cycle)
{
pcout << "Cycle " << cycle << std::endl;
if (cycle == 0)
{
GridGenerator::hyper_cube(triangulation, 0., 1.);
triangulation.refine_global(4);
setup_system();
VectorTools::project(dof_handler,
boundary_constraints,
QGauss<dim>(fe.degree + 2),
//Solution<dim>(0.),
ZeroFunction<dim>(),
solution);
}
else{
if(cycle < 3)
refine_mesh();
}
if(cycle == 0)
output_results(0);
assemble_rhs();
solve();
output_results(cycle + 1);
pcout << std::endl;
};
}
} // namespace Step37
int main(int argc, char *argv[])
{
try
{
using namespace Step37;
Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
LaplaceProblem<dimension> laplace_problem;
laplace_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;
}
