I am trying to port example 15 to a MPI-supported version, using the
geometric multigrid method as preconditioner (to get more familiar with the
implementation). I followed example 50 while doing that, but currently I
get the error
An error occurred in line <3866> of file
<~/Downloads/git-files/dealii/include/deal.II/dofs/dof_accessor.templates.h>
in function
void dealii::DoFCellAccessor<DoFHandlerType, lda>::get_dof_values(const
InputVector&, ForwardIterator, ForwardIterator) const [with InputVector =
dealii::TrilinosWrappers::MPI::Vector; ForwardIterator = double*;
DoFHandlerType = dealii::DoFHandler<2, 2>; bool level_dof_access = true]
The violated condition was:
this->is_artificial() == false
Additional information:
Can't ask for DoF indices on artificial cells.
as soon as I try to run my program using several MPI threads. I tried to
guard that part using if (cell->level_subdomain_id() ==
triangulation.locally_owned_subdomain()) as suggested in step-50, but that
did not help. Why am I still trying to request data from artificial cells,
even though I guarded that code part (at least I assume it guards it
accordingly).
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/d7eecf00-5fce-46c7-a191-a6d6261dac35%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2012 - 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.
*
* ---------------------------------------------------------------------
*
* Author: Sven Wetterauer, University of Heidelberg, 2012
*/
// @sect3{Include files}
// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/mg_transfer.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/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/meshworker/mesh_loop.h>
#include <fstream>
#include <iostream>
// We will use adaptive mesh refinement between Newton iterations. To do so,
// we need to be able to work with a solution on the new mesh, although it was
// computed on the old one. The SolutionTransfer class transfers the solution
// from the old to the new mesh:
#include <deal.II/numerics/solution_transfer.h>
// We then open a namespace for this program and import everything from the
// dealii namespace into it, as in previous programs:
namespace Step15
{
using namespace dealii;
namespace LA{
using namespace dealii::LinearAlgebraTrilinos;
}
template <int dim>
struct ScratchData
{
ScratchData(const Mapping<dim> & mapping,
const FiniteElement<dim> &fe,
const unsigned int quadrature_degree,
const UpdateFlags update_flags)
: fe_values(mapping, fe, QGauss<dim>(quadrature_degree), update_flags)
{}
ScratchData(const ScratchData<dim> &scratch_data)
: fe_values(scratch_data.fe_values.get_mapping(),
scratch_data.fe_values.get_fe(),
scratch_data.fe_values.get_quadrature(),
scratch_data.fe_values.get_update_flags())
{}
FEValues<dim> fe_values;
};
struct CopyData
{
unsigned int level;
FullMatrix<double> cell_matrix;
Vector<double> cell_rhs;
std::vector<types::global_dof_index> local_dof_indices;
template <class Iterator>
void reinit(const Iterator &cell, unsigned int dofs_per_cell)
{
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
cell_rhs.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
cell->get_active_or_mg_dof_indices(local_dof_indices);
level = cell->level();
}
};
// @sect3{The <code>MinimalSurfaceProblem</code> class template}
// The class template is basically the same as in step-6. Three additions
// are made:
// - There are two solution vectors, one for the Newton update
// $\delta u^n$, and one for the current iterate $u^n$.
// - The <code>setup_system</code> function takes an argument that denotes
// whether this is the first time it is called or not. The difference is
// that the first time around we need to distribute the degrees of freedom
// and set the solution vector for $u^n$ to the correct size. The following
// times, the function is called after we have already done these steps as
// part of refining the mesh in <code>refine_mesh</code>.
// - We then also need new functions: <code>set_boundary_values()</code>
// takes care of setting the boundary values on the solution vector
// correctly, as discussed at the end of the
// introduction. <code>compute_residual()</code> is a function that computes
// the norm of the nonlinear (discrete) residual. We use this function to
// monitor convergence of the Newton iteration. The function takes a step
// length $\alpha^n$ as argument to compute the residual of $u^n + \alpha^n
// \; \delta u^n$. This is something one typically needs for step length
// control, although we will not use this feature here. Finally,
// <code>determine_step_length()</code> computes the step length $\alpha^n$
// in each Newton iteration. As discussed in the introduction, we here use a
// fixed step length and leave implementing a better strategy as an
// exercise.
template <int dim>
class MinimalSurfaceProblem
{
public:
MinimalSurfaceProblem();
~MinimalSurfaceProblem();
void run();
private:
template <class Iterator>
void cell_worker(const Iterator & cell,
ScratchData<dim> &scratch_data,
CopyData & copy_data);
template <class Iterator>
void mg_cell_worker(const Iterator & cell,
ScratchData<dim> &scratch_data,
CopyData & copy_data);
void setup_system();
void assemble_system();
void assemble_multigrid();
void solve();
void refine_mesh();
void set_boundary_values();
double compute_residual(const double alpha) const;
double determine_step_length() const;
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FE_Q<dim> fe;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
AffineConstraints<double> hanging_node_constraints;
AffineConstraints<double> newton_constraints;
SparsityPattern sparsity_pattern;
using matrix_t = LA::MPI::SparseMatrix;
using vector_t = LA::MPI::Vector;
vector_t present_solution;
vector_t newton_update;
vector_t system_rhs;
matrix_t system_matrix;
MGLevelObject<SparsityPattern> mg_sparsity_patterns;
MGLevelObject<SparsityPattern> mg_interface_sparsity_patterns;
MGLevelObject<matrix_t> mg_matrices;
MGLevelObject<matrix_t> mg_interface_matrices;
MGConstrainedDoFs mg_constrained_dofs;
ConditionalOStream pcout;
TimerOutput computing_timer;
};
// @sect3{Boundary condition}
// The boundary condition is implemented just like in step-4. It is chosen
// as $g(x,y)=\sin(2 \pi (x+y))$:
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues()
: Function<dim>()
{}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
template <int dim>
double BoundaryValues<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
return std::sin(2 * numbers::PI * (p[0] + p[1]));
}
// @sect3{The <code>MinimalSurfaceProblem</code> class implementation}
// @sect4{MinimalSurfaceProblem::MinimalSurfaceProblem}
// The constructor and destructor of the class are the same as in the first
// few tutorials.
template <int dim>
MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
: mpi_communicator(MPI_COMM_WORLD),
triangulation(MPI_COMM_WORLD,
typename Triangulation<dim>::MeshSmoothing(
Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening |
Triangulation<dim>::limit_level_difference_at_vertices),
parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
dof_handler(triangulation),
fe(2),
pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
computing_timer(mpi_communicator,
pcout,
TimerOutput::summary,
TimerOutput::wall_times)
{}
template <int dim>
MinimalSurfaceProblem<dim>::~MinimalSurfaceProblem()
{
dof_handler.clear();
}
// @sect4{MinimalSurfaceProblem::setup_system}
// As always in the setup-system function, we setup the variables of the
// finite element method. There are same differences to step-6, because
// there we start solving the PDE from scratch in every refinement cycle
// whereas here we need to take the solution from the previous mesh onto the
// current mesh. Consequently, we can't just reset solution vectors. The
// argument passed to this function thus indicates whether we can
// distributed degrees of freedom (plus compute constraints) and set the
// solution vector to zero or whether this has happened elsewhere already
// (specifically, in <code>refine_mesh()</code>).
template <int dim>
void MinimalSurfaceProblem<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
present_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
hanging_node_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler,
hanging_node_constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
BoundaryValues<dim>(),
hanging_node_constraints);
hanging_node_constraints.close();
newton_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler,
newton_constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(),
newton_constraints);
newton_constraints.close();
//hanging_node_constraints.distribute(present_solution);
// The remaining parts of the function are the same as in step-6.
newton_update.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
system_rhs.reinit(locally_owned_dofs,
mpi_communicator);
DynamicSparsityPattern dsp(dof_handler.n_dofs(),
dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
hanging_node_constraints,
false);
SparsityTools::distribute_sparsity_pattern(
dsp,
dof_handler.compute_n_locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs);
sparsity_pattern.copy_from(dsp);
//hanging_node_constraints.condense(dsp);
// system_matrix.reinit(locally_owned_dofs,
// locally_owned_dofs,
// dsp,
// mpi_communicator);
system_matrix.reinit(sparsity_pattern);
std::set<types::boundary_id> dirichlet_boundary_ids = {0};
mg_constrained_dofs.clear();
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
dirichlet_boundary_ids);
const unsigned int n_levels = triangulation.n_levels();
mg_interface_matrices.resize(0, n_levels - 1);
mg_matrices.resize(0, n_levels - 1);
mg_sparsity_patterns.resize(0, n_levels - 1);
mg_interface_sparsity_patterns.resize(0, n_levels - 1);
for(size_t level = 0; level < n_levels; ++level){
DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
dof_handler.n_dofs(level));
MGTools::make_sparsity_pattern(dof_handler,
dsp,
level);
mg_sparsity_patterns[level].copy_from(dsp);
mg_matrices[level].reinit(dof_handler.locally_owned_mg_dofs(level),
dof_handler.locally_owned_mg_dofs(level),
dsp,
mpi_communicator,
true);
mg_interface_matrices[level].reinit(dof_handler.locally_owned_mg_dofs(level),
dof_handler.locally_owned_mg_dofs(level),
dsp,
mpi_communicator,
true);
}
}
template <int dim>
template <class Iterator>
void MinimalSurfaceProblem<dim>::cell_worker(const Iterator &cell,
ScratchData<dim> &scratch_data,
CopyData ©_data)
{
FEValues<dim> &fe_values = scratch_data.fe_values;
fe_values.reinit(cell);
const size_t dofs_per_cell = fe_values.get_fe().dofs_per_cell;
const size_t n_q_points = fe_values.get_quadrature().size();
copy_data.reinit(cell, dofs_per_cell);
const std::vector<double> &JxW = fe_values.get_JxW_values();
std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
vector_t local_evaluation_point(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
local_evaluation_point = present_solution;
fe_values.get_function_gradients(local_evaluation_point,
old_solution_gradients);
for (unsigned int q = 0; q < n_q_points; ++q)
{
const double coeff =
1.0 / std::sqrt(1 + old_solution_gradients[q] *
old_solution_gradients[q]);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j){
copy_data.cell_matrix(i, j) +=
(((fe_values.shape_grad(i, q) // ((\nabla \phi_i
* coeff // * a_n
* fe_values.shape_grad(j, q)) // * \nabla \phi_j)
- // -
(fe_values.shape_grad(i, q) // (\nabla \phi_i
* coeff * coeff * coeff // * a_n^3
* (fe_values.shape_grad(j, q) // * (\nabla \phi_j
* old_solution_gradients[q]) // * \nabla u_n)
* old_solution_gradients[q])) // * \nabla u_n)))
* JxW[q]); // * dx
}
copy_data.cell_rhs(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
* coeff // * a_n
* old_solution_gradients[q] // * u_n
* JxW[q]); // * dx
}
}
}
template <int dim>
template <class Iterator>
void MinimalSurfaceProblem<dim>::mg_cell_worker(const Iterator &cell,
ScratchData<dim> &scratch_data,
CopyData ©_data)
{
if (cell->level_subdomain_id() == triangulation.locally_owned_subdomain()){
FEValues<dim> &fe_values = scratch_data.fe_values;
fe_values.reinit(cell);
const size_t dofs_per_cell = fe_values.get_fe().dofs_per_cell;
const size_t n_q_points = fe_values.get_quadrature().size();
copy_data.reinit(cell, dofs_per_cell);
const std::vector<double> &JxW = fe_values.get_JxW_values();
std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
vector_t local_evaluation_point(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
local_evaluation_point = present_solution;
fe_values.get_function_gradients(present_solution,
old_solution_gradients);
for (unsigned int q = 0; q < n_q_points; ++q)
{
const double coeff =
1.0 / std::sqrt(1 + old_solution_gradients[q] *
old_solution_gradients[q]);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j){
copy_data.cell_matrix(i, j) +=
(((fe_values.shape_grad(i, q) // ((\nabla \phi_i
* coeff // * a_n
* fe_values.shape_grad(j, q)) // * \nabla \phi_j)
- // -
(fe_values.shape_grad(i, q) // (\nabla \phi_i
* coeff * coeff * coeff // * a_n^3
* (fe_values.shape_grad(j, q) // * (\nabla \phi_j
* old_solution_gradients[q]) // * \nabla u_n)
* old_solution_gradients[q])) // * \nabla u_n)))
* JxW[q]); // * dx
}
copy_data.cell_rhs(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
* coeff // * a_n
* old_solution_gradients[q] // * u_n
* JxW[q]); // * dx
}
}
}
}
// @sect4{MinimalSurfaceProblem::assemble_system}
// This function does the same as in the previous tutorials except that now,
// of course, the matrix and right hand side functions depend on the
// previous iteration's solution. As discussed in the introduction, we need
// to use zero boundary values for the Newton updates; we compute them at
// the end of this function.
//
// The top of the function contains the usual boilerplate code, setting up
// the objects that allow us to evaluate shape functions at quadrature
// points and temporary storage locations for the local matrices and
// vectors, as well as for the gradients of the previous solution at the
// quadrature points. We then start the loop over all cells:
template <int dim>
void MinimalSurfaceProblem<dim>::assemble_system()
{
TimerOutput::Scope t(computing_timer, "assembly");
MappingQ1<dim> mapping;
auto cell_worker = [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData<dim> &scratch_data,
CopyData ©_data){
this->cell_worker(cell, scratch_data, copy_data);
};
auto copier = [&](const CopyData &cd){
this->newton_constraints.distribute_local_to_global(cd.cell_matrix,
cd.cell_rhs,
cd.local_dof_indices,
system_matrix,
system_rhs);
};
const size_t n_gauss_points = fe.degree + 1;
ScratchData<dim> scratch_data(mapping,
fe,
n_gauss_points,
update_values
| update_gradients
| update_JxW_values
| update_quadrature_points);
MeshWorker::mesh_loop(dof_handler.begin_active(),
dof_handler.end(),
cell_worker,
copier,
scratch_data,
CopyData(),
MeshWorker::assemble_own_cells);
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
}
template <int dim>
void MinimalSurfaceProblem<dim>::assemble_multigrid(){
MappingQ1<dim> mapping;
const size_t n_levels = triangulation.n_levels();
std::vector<AffineConstraints<double>> boundary_constraints(n_levels);
for(size_t level = 0; level < n_levels; ++level){
IndexSet dofset;
DoFTools::extract_locally_relevant_level_dofs(dof_handler,
level,
dofset);
boundary_constraints[level].reinit(dofset);
boundary_constraints[level].add_lines(mg_constrained_dofs.get_refinement_edge_indices(level));
boundary_constraints[level].add_lines(mg_constrained_dofs.get_boundary_indices(level));
boundary_constraints[level].close();
}
auto mg_cell_worker = [&](const typename DoFHandler<dim>::level_cell_iterator &cell,
ScratchData<dim> &scratch_data,
CopyData ©_data){
this->mg_cell_worker(cell,
scratch_data,
copy_data);
};
auto copier = [&](const CopyData &cd){
boundary_constraints[cd.level].distribute_local_to_global(cd.cell_matrix,
cd.local_dof_indices,
mg_matrices[cd.level]);
const size_t dofs_per_cell = cd.local_dof_indices.size();
for(size_t i = 0; i < dofs_per_cell; ++i)
for(size_t j = 0; j < dofs_per_cell; ++j)
if(mg_constrained_dofs.is_interface_matrix_entry(cd.level,
cd.local_dof_indices[i],
cd.local_dof_indices[j])){
mg_interface_matrices[cd.level].add(cd.local_dof_indices[i],
cd.local_dof_indices[j],
cd.cell_matrix(i, j));
}
};
const size_t n_gauss_points = fe.degree + 1;
ScratchData<dim> scratch_data(mapping,
fe,
n_gauss_points,
update_values
| update_gradients
| update_JxW_values
| update_quadrature_points);
MeshWorker::mesh_loop(dof_handler.begin_mg(),
dof_handler.end_mg(),
mg_cell_worker,
copier,
scratch_data,
CopyData(),
MeshWorker::assemble_own_cells);
for (unsigned int i = 0; i < triangulation.n_global_levels(); ++i)
{
mg_matrices[i].compress(VectorOperation::add);
mg_interface_matrices[i].compress(VectorOperation::add);
}
}
// @sect4{MinimalSurfaceProblem::solve}
// The solve function is the same as always. At the end of the solution
// process we update the current solution by setting
// $u^{n+1}=u^n+\alpha^n\;\delta u^n$.
template <int dim>
void MinimalSurfaceProblem<dim>::solve()
{
MGTransferPrebuilt<vector_t> mg_transfer(mg_constrained_dofs);
mg_transfer.build_matrices(dof_handler);
matrix_t &coarse_matrix = mg_matrices[0];
SolverControl coarse_solver_control(1000, 1e-10, false, false);
SolverGMRES<vector_t> coarse_solver(coarse_solver_control);
PreconditionIdentity id;
MGCoarseGridIterativeSolver<vector_t,
SolverGMRES<vector_t>,
matrix_t,
PreconditionIdentity>
coarse_grid_solver(coarse_solver, coarse_matrix, id);
using Smoother = LA::MPI::PreconditionJacobi;
MGSmootherPrecondition<matrix_t, Smoother, vector_t> mg_smoother;
mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5));
mg_smoother.set_steps(2);
mg::Matrix<vector_t> mg_matrix(mg_matrices);
mg::Matrix<vector_t> mg_interface_up(mg_interface_matrices);
mg::Matrix<vector_t> mg_interface_down(mg_interface_matrices);
Multigrid<vector_t> mg(mg_matrix, coarse_grid_solver, mg_transfer, mg_smoother, mg_smoother);
mg.set_edge_matrices(mg_interface_down, mg_interface_up);
PreconditionMG<dim, vector_t, MGTransferPrebuilt<vector_t>> preconditioner(dof_handler,
mg,
mg_transfer);
vector_t completely_distributed_update(locally_owned_dofs,
mpi_communicator);
SolverControl solver_control(system_rhs.size(),
system_rhs.l2_norm() * 1e-6);
SolverGMRES<vector_t> solver(solver_control);
solver.solve(system_matrix, completely_distributed_update, system_rhs, preconditioner);
newton_constraints.distribute(completely_distributed_update);
newton_update = completely_distributed_update;
const double alpha = 0.1;//determine_step_length();
newton_update *= alpha;
present_solution += newton_update;
}
// @sect4{MinimalSurfaceProblem::refine_mesh}
// The first part of this function is the same as in step-6... However,
// after refining the mesh we have to transfer the old solution to the new
// one which we do with the help of the SolutionTransfer class. The process
// is slightly convoluted, so let us describe it in detail:
template <int dim>
void MinimalSurfaceProblem<dim>::refine_mesh()
{
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> *>(),
present_solution,
estimated_error_per_cell);
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(triangulation,
estimated_error_per_cell,
0.3,
0.03);
triangulation.prepare_coarsening_and_refinement();
parallel::distributed::SolutionTransfer<dim, vector_t> solution_transfer(dof_handler);
solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
triangulation.execute_coarsening_and_refinement();
// Vector<double> old_tmp(dof_handler.n_dofs());
// old_tmp = present_solution;
setup_system();
vector_t tmp(locally_owned_dofs,
mpi_communicator);
solution_transfer.interpolate(tmp);
hanging_node_constraints.distribute(tmp);
present_solution = tmp;
}
template <int dim>
double MinimalSurfaceProblem<dim>::compute_residual(const double alpha) const
{
vector_t residual(locally_owned_dofs,
mpi_communicator);
vector_t evaluation_point(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
vector_t local_newton_update(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
local_newton_update = newton_update;
local_newton_update *= alpha;
evaluation_point = present_solution;
evaluation_point += local_newton_update;
const QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_gradients | update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> cell_residual(dofs_per_cell);
std::vector<Tensor<1, dim>> gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
if(cell->is_locally_owned()){
cell_residual = 0;
fe_values.reinit(cell);
fe_values.get_function_gradients(evaluation_point, gradients);
for (unsigned int q = 0; q < n_q_points; ++q)
{
const double coeff = 1 / std::sqrt(1 + gradients[q] * gradients[q]);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_residual(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
* coeff // * a_n
* gradients[q] // * u_n
* fe_values.JxW(q)); // * dx
}
cell->get_dof_indices(local_dof_indices);
newton_constraints.distribute_local_to_global(cell_residual,
local_dof_indices,
residual);
}
}
residual.compress(VectorOperation::add);
return residual.l2_norm();
}
// @sect4{MinimalSurfaceProblem::determine_step_length}
// As discussed in the introduction, Newton's method frequently does not
// converge if we always take full steps, i.e., compute $u^{n+1}=u^n+\delta
// u^n$. Rather, one needs a damping parameter (step length) $\alpha^n$ and
// set $u^{n+1}=u^n+\alpha^n\delta u^n$. This function is the one called
// to compute $\alpha^n$.
//
// Here, we simply always return 0.1. This is of course a sub-optimal
// choice: ideally, what one wants is that the step size goes to one as we
// get closer to the solution, so that we get to enjoy the rapid quadratic
// convergence of Newton's method. We will discuss better strategies below
// in the results section.
template <int dim>
double MinimalSurfaceProblem<dim>::determine_step_length() const
{
auto residual_0 = compute_residual(0);
auto residual_25 = compute_residual(0.25);
auto residual_5 = compute_residual(0.5);
auto residual_75 = compute_residual(0.75);
auto residual_1 = compute_residual(1);
if(residual_1 < residual_0)
return 0.1;
else{
if(residual_1 == 0 || residual_1 == 0)
return 0.1;
else{
// const double eta_val = (residual_1 / residual_0) * (residual_1 / residual_0);
// return 0.1;//(sqrt(1 + 6 * eta_val) - 1) / (3 * eta_val);
if(residual_75 < residual_0)
return 0.1;
else if(residual_5 < residual_0)
return 0.5;
else if(residual_25 < residual_0)
return 0.25;
else
return 0.1;
}
}
}
// @sect4{MinimalSurfaceProblem::run}
// In the run function, we build the first grid and then have the top-level
// logic for the Newton iteration. The function has two variables, one that
// indicates whether this is the first time we solve for a Newton update and
// one that indicates the refinement level of the mesh:
template <int dim>
void MinimalSurfaceProblem<dim>::run()
{
unsigned int refinement = 0;
bool first_step = true;
// As described in the introduction, the domain is the unit disk around
// the origin, created in the same way as shown in step-6. The mesh is
// globally refined twice followed later on by several adaptive cycles:
GridGenerator::hyper_ball(triangulation);
triangulation.refine_global(2);
// The Newton iteration starts next. During the first step we do not have
// information about the residual prior to this step and so we continue
// the Newton iteration until we have reached at least one iteration and
// until residual is less than $10^{-3}$.
//
// At the beginning of the loop, we do a bit of setup work. In the first
// go around, we compute the solution on the twice globally refined mesh
// after setting up the basic data structures and ensuring that the first
// Newton iterate already has the correct boundary values. In all
// following mesh refinement loops, the mesh will be refined adaptively.
double previous_res = 0;
while (first_step || (previous_res > 1e-3))
{
if (first_step == true)
{
pcout << "******** Initial mesh "
<< " ********" << std::endl;
setup_system();
vector_t local_solution(locally_owned_dofs,
mpi_communicator);
hanging_node_constraints.distribute(local_solution);
present_solution = local_solution;
first_step = false;
}
else
{
++refinement;
pcout << "******** Refined mesh " << refinement << " ********"
<< std::endl;
refine_mesh();
}
// On every mesh we do exactly five Newton steps. We print the initial
// residual here and then start the iterations on this mesh.
//
// In every Newton step the system matrix and the right hand side have
// to be computed first, after which we store the norm of the right
// hand side as the residual to check against when deciding whether to
// stop the iterations. We then solve the linear system (the function
// also updates $u^{n+1}=u^n+\alpha^n\;\delta u^n$) and output the
// residual at the end of this Newton step:
pcout << " Initial residual: " << compute_residual(0) << std::endl;
for (unsigned int inner_iteration = 0; inner_iteration < 5;
++inner_iteration)
{
assemble_system();
assemble_multigrid();
previous_res = system_rhs.l2_norm();
solve();
pcout << " Residual: " << compute_residual(0) << std::endl;
}
// Every fifth iteration, i.e., just before we refine the mesh again,
// we output the solution as well as the Newton update. This happens
// as in all programs before:
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(present_solution, "solution");
Vector<float> subdomain(triangulation.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector(subdomain, "subdomain");
data_out.build_patches();
const std::string filename =
("solution-" + Utilities::int_to_string(refinement, 2) + "." +
Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
std::ofstream output(filename + ".vtu");
data_out.write_vtu(output);
if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{
std::vector<std::string> filenames;
for (unsigned int i = 0;
i < Utilities::MPI::n_mpi_processes(mpi_communicator);
++i)
filenames.push_back("solution-" + Utilities::int_to_string(refinement, 2) +
"." + Utilities::int_to_string(i, 4) + ".vtu");
std::ofstream master_output(
"solution-" + Utilities::int_to_string(refinement, 2) + ".pvtu");
data_out.write_pvtu_record(master_output, filenames);
}
}
}
} // namespace Step15
// @sect4{The main function}
// Finally the main function. This follows the scheme of all other main
// functions:
int main(int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace Step15;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
MinimalSurfaceProblem<2> laplace_problem_2d;
laplace_problem_2d.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;
}
