Dear all,
I'm trying to use a block-diagonal preconditioner for the mixed Poisson
problem from step-20, but haven't made it yet.
The diagonal preconditioner of interest is from Powell & Silvester,2003
<https://www.google.com/url?q=https%3A%2F%2Fepubs.siam.org%2Fdoi%2Fabs%2F10.1137%2FS0895479802404428&sa=D&sntz=1&usg=AFQjCNGCFHHh92fZkVu4-ueHwkaDbH_SiA>.
Specifically, for the matrix system after discretization with the RT0-DG0
element pair
[A BT
B 0],
they suggested a block diagonal preconditioner
P =
[DA 0
0 PS],
where DA is diag(A), and PS is AMG applied to the approximated Schur
complement SD = B(DA^-1)BT.
As the above idea is quite similar to that in step-55, so I borrowed quite
a few code blocks from there, e.g. InverseMatrix and
BlockDiagonalPreconditioner. The resulting solve() with a SolverMinRes solver
is as follows. (The rest of the code is mainly from step-20)
template <int dim>
void MixedLaplaceProblem<dim>::solve()
{
TrilinosWrappers::PreconditionJacobi diagA;
diagA.initialize(system_matrix.block(0, 0));//DA
ApproximateSchurComplement approximate_schur(system_matrix, solution);
//SD
InverseMatrix<ApproximateSchurComplement> approximate_inverse(
approximate_schur);//PS
const BlockDiagonalPreconditioner<TrilinosWrappers::PreconditionJacobi,
InverseMatrix<ApproximateSchurComplement> >
preconditioner(diagA, approximate_inverse);
SolverControl solver_control(system_matrix.m(),
1e-6 * system_rhs.l2_norm());
SolverMinRes<TrilinosWrappers::MPI::BlockVector> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs,
preconditioner
);
std::cout<<"last step = "<<solver_control.last_step()<<std::endl;
}
However, the code can't compile, and gives the following error message:
/home/dealii/test.cc:336:5: error: no matching function for call to ‘
dealii::TrilinosWrappers::PreconditionAMG::initialize(const Step20::
ApproximateSchurComplement&, dealii::TrilinosWrappers::PreconditionAMG::
AdditionalData&)’
preconditioner.initialize(*matrix, data);
^
In file included from /home/dealii/dealii-v9.0.0/include/deal.II/lac/
generic_linear_algebra.h:141:0,
from /home/dealii/test.cc:31:
/home/dealii/dealii-v9.0.0/include/deal.II/lac/trilinos_precondition.h:1516
:10: note: candidate: void dealii::TrilinosWrappers::PreconditionAMG::
initialize(const dealii::TrilinosWrappers::SparseMatrix&, const dealii::
TrilinosWrappers::PreconditionAMG::AdditionalData&)
void initialize (const SparseMatrix &matrix,
The error is related to the use of PreconditionAMG in my implementation of
InverseMatrix<MatrixType>::vmult:
template <class MatrixType>
void InverseMatrix<MatrixType>::vmult (LinearAlgebraTrilinos::MPI::Vector
&dst,
const LinearAlgebraTrilinos::MPI::
Vector &src) const
{
SolverControl solver_control(src.size(),
1e-6 * src.l2_norm());
SolverCG<TrilinosWrappers::MPI::Vector> cg (solver_control);
dst = 0;
LinearAlgebraTrilinos::MPI::PreconditionAMG preconditioner;
LinearAlgebraTrilinos::MPI::PreconditionAMG::AdditionalData data;
preconditioner.initialize(*matrix, data);
// PreconditionIdentity preconditioner;//A Identity preconditioner
works fine here
cg.solve (*matrix, dst, src, preconditioner);
std::cout<<"inner last step = "<<solver_control.last_step()<<std::endl;
}
I find a conflict here is that::
1) TrilinosWrappers::PreconditionAMG can only be initialized with a matrix,
but not a ApproximateSchurComplement;
2) step-20 uses ApproximateSchurComplement because we actually don't want
to/can't assemble the Schur matrix or its inverse.
Could you please give me a hint on how to solve this conflict in deal.II?
Is it actually possible to assemble this B(DA^-1)BT directly?
Attached is a minimal (but not working) code.
Best regards,
Charlie
--
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/73183270-5edb-4c8f-a998-17da5cb12221%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/packaged_operation.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/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
namespace Step20
{
using namespace dealii;
template <int dim>
class MixedLaplaceProblem
{
public:
MixedLaplaceProblem(const unsigned int degree);
void run();
private:
void make_grid_and_dofs();
void assemble_system();
void solve();
void compute_errors() const;
void output_results() const;
const unsigned int degree;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
BlockSparsityPattern sparsity_pattern;
std::vector<IndexSet> stokes_partitioning;
TrilinosWrappers::BlockSparseMatrix system_matrix;
TrilinosWrappers::MPI::BlockVector solution;
TrilinosWrappers::MPI::BlockVector system_rhs;
bool verbose;
};
template <class PreconditionerA, class PreconditionerS>
class BlockDiagonalPreconditioner : public Subscriptor
{
public:
BlockDiagonalPreconditioner(const PreconditionerA &preconditioner_A,
const PreconditionerS &preconditioner_S);
void vmult(LinearAlgebraTrilinos::MPI::BlockVector & dst,
const LinearAlgebraTrilinos::MPI::BlockVector &src) const;
private:
const PreconditionerA &preconditioner_A;
const PreconditionerS &preconditioner_S;
};
template <class PreconditionerA, class PreconditionerS>
BlockDiagonalPreconditioner<PreconditionerA, PreconditionerS>::
BlockDiagonalPreconditioner(const PreconditionerA &preconditioner_A,
const PreconditionerS &preconditioner_S)
: preconditioner_A(preconditioner_A)
, preconditioner_S(preconditioner_S)
{}
template <class PreconditionerA, class PreconditionerS>
void BlockDiagonalPreconditioner<PreconditionerA, PreconditionerS>::vmult(
LinearAlgebraTrilinos::MPI::BlockVector & dst,
const LinearAlgebraTrilinos::MPI::BlockVector &src) const
{
preconditioner_A.vmult(dst.block(0), src.block(0));
preconditioner_S.vmult(dst.block(1), src.block(1));
}
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide()
: Function<dim>(1)
{}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
template <int dim>
class PressureBoundaryValues : public Function<dim>
{
public:
PressureBoundaryValues()
: Function<dim>(1)
{}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
template <int dim>
class ExactSolution : public Function<dim>
{
public:
ExactSolution()
: Function<dim>(dim + 1)
{}
virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override;
};
template <int dim>
double RightHandSide<dim>::value(const Point<dim> & /*p*/,
const unsigned int /*component*/) const
{
return 0;
}
template <int dim>
double
PressureBoundaryValues<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
const double alpha = 0.3;
const double beta = 1;
return -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
alpha * p[0] * p[0] * p[0] / 6);
}
template <int dim>
void ExactSolution<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
{
Assert(values.size() == dim + 1,
ExcDimensionMismatch(values.size(), dim + 1));
const double alpha = 0.3;
const double beta = 1;
values(0) = alpha * p[1] * p[1] / 2 + beta - alpha * p[0] * p[0] / 2;
values(1) = alpha * p[0] * p[1];
values(2) = -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
alpha * p[0] * p[0] * p[0] / 6);
}
template <int dim>
class KInverse : public TensorFunction<2, dim>
{
public:
KInverse()
: TensorFunction<2, dim>()
{}
virtual void value_list(const std::vector<Point<dim>> &points,
std::vector<Tensor<2, dim>> &values) const override;
};
template <int dim>
void KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
std::vector<Tensor<2, dim>> & values) const
{
(void)points;
AssertDimension(points.size(), values.size());
for (auto &value : values)
{
value.clear();
for (unsigned int d = 0; d < dim; ++d)
value[d][d] = 1.;
}
}
template <int dim>
MixedLaplaceProblem<dim>::MixedLaplaceProblem(const unsigned int degree)
: degree(degree)
, fe(FE_RaviartThomas<dim>(degree), 1, FE_DGQ<dim>(degree), 1)
, dof_handler(triangulation)
{
verbose = false;
}
template <int dim>
void MixedLaplaceProblem<dim>::make_grid_and_dofs()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(5);
dof_handler.distribute_dofs(fe);
DoFRenumbering::component_wise(dof_handler);
std::vector<types::global_dof_index> dofs_per_component(dim + 1);
DoFTools::count_dofs_per_component(dof_handler, dofs_per_component);
const unsigned int n_u = dofs_per_component[0],
n_p = dofs_per_component[dim];
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< "Total number of cells: " << triangulation.n_cells()
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< " (" << n_u << '+' << n_p << ')' << std::endl;
BlockDynamicSparsityPattern dsp(2, 2);
dsp.block(0, 0).reinit(n_u, n_u);
dsp.block(1, 0).reinit(n_p, n_u);
dsp.block(0, 1).reinit(n_u, n_p);
dsp.block(1, 1).reinit(n_p, n_p);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern(dof_handler, dsp);
system_matrix.reinit(dsp);
stokes_partitioning.resize(2);
stokes_partitioning[0] = complete_index_set(n_u);
stokes_partitioning[1] = complete_index_set(n_p);
solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
system_rhs.reinit(stokes_partitioning, MPI_COMM_WORLD);
}
template <int dim>
void MixedLaplaceProblem<dim>::assemble_system()
{
QGauss<dim> quadrature_formula(degree + 2);
QGauss<dim - 1> face_quadrature_formula(degree + 2);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> fe_face_values(fe,
face_quadrature_formula,
update_values | update_normal_vectors |
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();
const unsigned int n_face_q_points = face_quadrature_formula.size();
FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const RightHandSide<dim> right_hand_side;
const PressureBoundaryValues<dim> pressure_boundary_values;
const KInverse<dim> k_inverse;
std::vector<double> rhs_values(n_q_points);
std::vector<double> boundary_values(n_face_q_points);
std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(dim);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
local_matrix = 0;
local_rhs = 0;
right_hand_side.value_list(fe_values.get_quadrature_points(),
rhs_values);
k_inverse.value_list(fe_values.get_quadrature_points(),
k_inverse_values);
for (unsigned int q = 0; q < n_q_points; ++q)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
const double div_phi_i_u = fe_values[velocities].divergence(i, q);
const double phi_i_p = fe_values[pressure].value(i, q);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const Tensor<1, dim> phi_j_u =
fe_values[velocities].value(j, q);
const double div_phi_j_u =
fe_values[velocities].divergence(j, q);
const double phi_j_p = fe_values[pressure].value(j, q);
local_matrix(i, j) +=
(phi_i_u * k_inverse_values[q] * phi_j_u
- phi_i_p * div_phi_j_u
- div_phi_i_u * phi_j_p)
* fe_values.JxW(q);
}
local_rhs(i) += -phi_i_p * rhs_values[q] * fe_values.JxW(q);
}
for (unsigned int face_n = 0;
face_n < GeometryInfo<dim>::faces_per_cell;
++face_n)
if (cell->at_boundary(face_n))
{
fe_face_values.reinit(cell, face_n);
pressure_boundary_values.value_list(
fe_face_values.get_quadrature_points(), boundary_values);
for (unsigned int q = 0; q < n_face_q_points; ++q)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
local_rhs(i) += -(fe_face_values[velocities].value(i, q) *
fe_face_values.normal_vector(q) *
boundary_values[q] *
fe_face_values.JxW(q));
}
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
local_matrix(i, j));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += local_rhs(i);
}
if (verbose)
{
std::cout<<" *** \nafter constrants.condense\n";
system_matrix.print(std::cout);
system_rhs.print(std::cout);
}
}
template <class MatrixType>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix(const MatrixType &m);
void vmult(LinearAlgebraTrilinos::MPI::Vector &dst,
const LinearAlgebraTrilinos::MPI::Vector &src) const;
private:
const SmartPointer<const MatrixType> matrix;
};
template <class MatrixType>
InverseMatrix<MatrixType>::InverseMatrix (const MatrixType &m)
:
matrix (&m)
{}
template <class MatrixType>
void InverseMatrix<MatrixType>::vmult (LinearAlgebraTrilinos::MPI::Vector &dst,
const LinearAlgebraTrilinos::MPI::Vector &src) const
{
SolverControl solver_control(src.size(),
1e-6 * src.l2_norm());
SolverCG<TrilinosWrappers::MPI::Vector> cg (solver_control);
dst = 0;
LinearAlgebraTrilinos::MPI::PreconditionAMG preconditioner;
LinearAlgebraTrilinos::MPI::PreconditionAMG::AdditionalData data;
preconditioner.initialize(*matrix, data);
// PreconditionIdentity preconditioner;//a Identity preconditioner works fine
cg.solve (*matrix, dst, src, preconditioner);
std::cout<<"inner last step = "<<solver_control.last_step()<<std::endl;
}
class ApproximateSchurComplement : public Subscriptor
{
public:
ApproximateSchurComplement (const LinearAlgebraTrilinos::MPI::BlockSparseMatrix &A, const LinearAlgebraTrilinos::MPI::BlockVector &V);
void vmult (LinearAlgebraTrilinos::MPI::Vector &dst,
const LinearAlgebraTrilinos::MPI::Vector &src) const;
private:
const SmartPointer<const LinearAlgebraTrilinos::MPI::BlockSparseMatrix > system_matrix;
mutable LinearAlgebraTrilinos::MPI::Vector tmp1, tmp2;
};
ApproximateSchurComplement::ApproximateSchurComplement
(const LinearAlgebraTrilinos::MPI::BlockSparseMatrix &A, const LinearAlgebraTrilinos::MPI::BlockVector &V) :
system_matrix (&A),
tmp1(V.block(0)),
tmp2(V.block(0))
{}
void
ApproximateSchurComplement::vmult
(LinearAlgebraTrilinos::MPI::Vector &dst,
const LinearAlgebraTrilinos::MPI::Vector &src) const
{
tmp1 = 0;
tmp2 = 0;
system_matrix->block(0,1).vmult (tmp1, src);
TrilinosWrappers::PreconditionJacobi diagA;
diagA.initialize(system_matrix->block(0, 0));
diagA.vmult(tmp2, tmp1);
system_matrix->block(1,0).vmult (dst, tmp2);
}
template <int dim>
void MixedLaplaceProblem<dim>::solve()
{
TrilinosWrappers::PreconditionJacobi diagA;
diagA.initialize(system_matrix.block(0, 0));
ApproximateSchurComplement approximate_schur(system_matrix, solution);
InverseMatrix<ApproximateSchurComplement> approximate_inverse(approximate_schur);
const BlockDiagonalPreconditioner<TrilinosWrappers::PreconditionJacobi,
InverseMatrix<ApproximateSchurComplement> >
preconditioner(diagA, approximate_inverse);
SolverControl solver_control(system_matrix.m(),
1e-6 * system_rhs.l2_norm());
SolverMinRes<TrilinosWrappers::MPI::BlockVector> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs,
preconditioner
);
std::cout<<"last step = "<<solver_control.last_step()<<std::endl;
}
template <int dim>
void MixedLaplaceProblem<dim>::compute_errors() const
{
const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
dim + 1);
ExactSolution<dim> exact_solution;
Vector<double> cellwise_errors(triangulation.n_active_cells());
QTrapez<1> q_trapez;
QIterated<dim> quadrature(q_trapez, degree + 2);
VectorTools::integrate_difference(dof_handler,
solution,
exact_solution,
cellwise_errors,
quadrature,
VectorTools::L2_norm,
&pressure_mask);
const double p_l2_error =
VectorTools::compute_global_error(triangulation,
cellwise_errors,
VectorTools::L2_norm);
VectorTools::integrate_difference(dof_handler,
solution,
exact_solution,
cellwise_errors,
quadrature,
VectorTools::L2_norm,
&velocity_mask);
const double u_l2_error =
VectorTools::compute_global_error(triangulation,
cellwise_errors,
VectorTools::L2_norm);
std::cout << "Errors: ||e_p||_L2 = " << p_l2_error
<< ", ||e_u||_L2 = " << u_l2_error << std::endl;
}
template <int dim>
void MixedLaplaceProblem<dim>::output_results() const
{
std::vector<std::string> solution_names(dim, "u");
solution_names.emplace_back("p");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
interpretation(dim,
DataComponentInterpretation::component_is_part_of_vector);
interpretation.push_back(DataComponentInterpretation::component_is_scalar);
DataOut<dim> data_out;
data_out.add_data_vector(dof_handler,
solution,
solution_names,
interpretation);
data_out.build_patches(degree + 1);
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
template <int dim>
void MixedLaplaceProblem<dim>::run()
{
make_grid_and_dofs();
assemble_system();
solve();
compute_errors();
output_results();
}
}
int main(int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace Step20;
Utilities::MPI::MPI_InitFinalize mpi_initialization(
argc, argv, numbers::invalid_unsigned_int);
AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
ExcMessage(
"This program can only be run in serial."));
MixedLaplaceProblem<2> mixed_laplace_problem(0);
mixed_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;
}##
# CMake script for the step-21 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "test")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
FIND_PACKAGE(deal.II 9.0.0 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
