* ---------------------------------------------------------------------
*
- * Author: Wolfgang Bangerth, Texas A&M University, 2005, 2006
+ * Authors: Wolfgang Bangerth, Texas A&M University, 2005, 2006;
+ * (port to LinearOperator:) Matthias Maier, 2019
*/
#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/precondition.h>
+// The only two new header files that deserve some attention are those for
+// the LinearOperator and PackagedOperation classes:
+#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>
// @sect3{Linear solvers and preconditioners}
- // The linear solvers and preconditioners we use in this example have been
- // discussed in significant detail already in the introduction. We will
- // therefore not discuss the rationale for these classes here any more, but
- // rather only comment on implementational aspects.
-
- // @sect4{The <code>InverseMatrix</code> class template}
-
- // There are a few places in this program where we will need either the
- // action of the inverse of the mass matrix or the action of the inverse of
- // the approximate Schur complement. Rather than explicitly calling
- // SolverCG::solve every time that we need to solve such a system, we will
- // wrap the action of either inverse in a simple class. The only things we
- // would like to note are that this class is derived from
- // <code>Subscriptor</code> and, as mentioned above, it stores a pointer to
- // the underlying matrix with a <code>SmartPointer</code> object. This class
- // also appears in step-21 and a more advanced version of it appears in
- // step-22.
- template <class MatrixType>
- class InverseMatrix : public Subscriptor
- {
- public:
- InverseMatrix(const MatrixType &m);
-
- void vmult(Vector<double> &dst, const Vector<double> &src) const;
-
- private:
- const SmartPointer<const MatrixType> matrix;
- };
-
-
- template <class MatrixType>
- InverseMatrix<MatrixType>::InverseMatrix(const MatrixType &m)
- : matrix(&m)
- {}
-
+ // The linear solvers and preconditioners we use in this example have
+ // been discussed in significant detail already in the introduction. We
+ // will therefore not discuss the rationale for our approach here any
+ // more, but rather only comment on some remaining implementational
+ // aspects.
- template <class MatrixType>
- void InverseMatrix<MatrixType>::vmult(Vector<double> & dst,
- const Vector<double> &src) const
- {
- // To make the control flow simpler, we recreate both the ReductionControl
- // and SolverCG objects every time this is called. This is not the most
- // efficient choice because SolverCG instances allocate memory whenever
- // they are created; this is just a tutorial so such inefficiencies are
- // acceptable for the sake of exposition.
- SolverControl solver_control(std::max<unsigned int>(src.size(), 200),
- 1e-8 * src.l2_norm());
- SolverCG<> cg(solver_control);
-
- dst = 0;
-
- cg.solve(*matrix, dst, src, PreconditionIdentity());
- }
-
-
- // @sect4{The <code>SchurComplement</code> class}
+ // @sect4{MixedLaplace::solve}
- // The next class is the Schur complement class. Its rationale has also been
- // discussed in length in the introduction. Like <code>InverseMatrix</code>,
- // this class is derived from Subscriptor and stores SmartPointer s
- // pointing to the system matrix and <code>InverseMatrix</code> wrapper.
- //
- // The <code>vmult</code> function requires two temporary vectors that we do
- // not want to re-allocate and free every time we call this function. Since
- // here, we have full control over the use of these vectors (unlike above,
- // where a class called by the <code>vmult</code> function required these
- // vectors, not the <code>vmult</code> function itself), we allocate them
- // directly, rather than going through the <code>VectorMemory</code>
- // mechanism. However, again, these member variables do not carry any state
- // between successive calls to the member functions of this class (i.e., we
- // never care what values they were set to the last time a member function
- // was called), we mark these vectors as <code>mutable</code>.
- //
- // The rest of the (short) implementation of this class is straightforward
- // if you know the order of matrix-vector multiplications performed by the
- // <code>vmult</code> function:
- class SchurComplement : public Subscriptor
+ // As already outlined in the introduction, the solve function consists
+ // essentially of two steps. First, we have to form the first equation
+ // involving the Schur complement and solve for the pressure (component 1
+ // of the solution). Then, we can reconstruct the velocities from the
+ // second equation (component 0 of the solution).
+ template <int dim>
+ void MixedLaplaceProblem<dim>::solve()
{
- public:
- SchurComplement(const BlockSparseMatrix<double> & A,
- const InverseMatrix<SparseMatrix<double>> &Minv);
-
- void vmult(Vector<double> &dst, const Vector<double> &src) const;
-
- private:
- const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
- const SmartPointer<const InverseMatrix<SparseMatrix<double>>> m_inverse;
-
- mutable Vector<double> tmp1, tmp2;
- };
+ // As a first step we declare references to all block components of the
+ // matrix, the right hand side and the solution vector that we will
+ // need.
+ const auto &M = system_matrix.block(0, 0);
+ const auto &B = system_matrix.block(0, 1);
+ const auto &F = system_rhs.block(0);
+ const auto &G = system_rhs.block(1);
- SchurComplement ::SchurComplement(
- const BlockSparseMatrix<double> & A,
- const InverseMatrix<SparseMatrix<double>> &Minv)
- : system_matrix(&A)
- , m_inverse(&Minv)
- , tmp1(A.block(0, 0).m())
- , tmp2(A.block(0, 0).m())
- {}
+ auto &U = solution.block(0);
+ auto &P = solution.block(1);
+ // Then, we will create corresponding LinearOperator objects and create
+ // the <code>op_M_inv</code> operator:
- void SchurComplement::vmult(Vector<double> & dst,
- const Vector<double> &src) const
- {
- system_matrix->block(0, 1).vmult(tmp1, src);
- m_inverse->vmult(tmp2, tmp1);
- system_matrix->block(1, 0).vmult(dst, tmp2);
- }
+ const auto op_M = linear_operator(M);
+ const auto op_B = linear_operator(B);
+ ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10);
+ SolverCG<> solver_M(reduction_control_M);
+ PreconditionJacobi<> preconditioner_M;
- // @sect4{The <code>ApproximateSchurComplement</code> class}
+ preconditioner_M.initialize(M);
- // The third component of our solver and preconditioner system is the class
- // that approximates the Schur complement with the method described in the
- // introduction. We will use this class to build a preconditioner for our
- // system matrix.
- class ApproximateSchurComplement : public Subscriptor
- {
- public:
- ApproximateSchurComplement(const BlockSparseMatrix<double> &A);
+ const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M);
- void vmult(Vector<double> &dst, const Vector<double> &src) const;
+ // This puts us in the position to be able to declare the Schur
+ // complement <code>op_S</code> and the approximate Schur complement
+ // <code>op_aS</code>:
- private:
- const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
+ const auto op_S = transpose_operator(op_B) * op_M_inv * op_B;
+ const auto op_aS =
+ transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B;
- mutable Vector<double> tmp1, tmp2;
- };
+ // We now create a preconditioner out of <code>op_aS</code> that
+ // applies a few number of CG iterations (until a very modest relative
+ // reduction of $10^{-16}$ is reached):
+ ReductionControl reduction_control_aS(2000, 1.e-18, 1.0e-6);
+ SolverCG<> solver_aS(reduction_control_aS);
+ PreconditionIdentity preconditioner_aS;
+ const auto preconditioner_S =
+ inverse_operator(op_aS, solver_aS, preconditioner_aS);
- ApproximateSchurComplement::ApproximateSchurComplement(
- const BlockSparseMatrix<double> &A)
- : system_matrix(&A)
- , tmp1(A.block(0, 0).m())
- , tmp2(A.block(0, 0).m())
- {}
+ // Now on to the first equation. The right hand side of it is
+ // $B^TM^{-1}F-G$, which is what we compute in the first few lines. We
+ // then solve the first equation with a CG solver and the
+ // preconditioner we just declared.
+ const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G;
- void ApproximateSchurComplement::vmult(Vector<double> & dst,
- const Vector<double> &src) const
- {
- system_matrix->block(0, 1).vmult(tmp1, src);
- system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
- system_matrix->block(1, 0).vmult(dst, tmp2);
- }
+ SolverControl solver_control_S(2000, 1.e-12);
+ SolverCG<> solver_S(solver_control_S);
- // @sect4{MixedLaplace::solve}
+ const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S);
- // After all these preparations, we can finally write the function that
- // actually solves the linear problem. We will go through the two parts it
- // has that each solve one of the two equations, the first one for the
- // pressure (component 1 of the solution), then the velocities (component 0
- // of the solution).
- template <int dim>
- void MixedLaplaceProblem<dim>::solve()
- {
- InverseMatrix<SparseMatrix<double>> inverse_mass(system_matrix.block(0, 0));
- Vector<double> tmp(solution.block(0).size());
+ P = op_S_inv * schur_rhs;
- // Now on to the first equation. The right hand side of it is
- // $B^TM^{-1}F-G$, which is what we compute in the first few lines:
- {
- SchurComplement schur_complement(system_matrix, inverse_mass);
- Vector<double> schur_rhs(solution.block(1).size());
- inverse_mass.vmult(tmp, system_rhs.block(0));
- system_matrix.block(1, 0).vmult(schur_rhs, tmp);
- schur_rhs -= system_rhs.block(1);
-
- // Now that we have the right hand side we can go ahead and solve for the
- // pressure, using our approximation of the inverse as a preconditioner:
- SolverControl solver_control(solution.block(1).size(),
- 1e-12 * schur_rhs.l2_norm());
- SolverCG<> cg(solver_control);
-
- ApproximateSchurComplement approximate_schur(system_matrix);
- InverseMatrix<ApproximateSchurComplement> approximate_inverse(
- approximate_schur);
- cg.solve(schur_complement,
- solution.block(1),
- schur_rhs,
- approximate_inverse);
-
- std::cout << solver_control.last_step()
- << " CG Schur complement iterations to obtain convergence."
- << std::endl;
- }
+ std::cout << solver_control_S.last_step()
+ << " CG Schur complement iterations to obtain convergence."
+ << std::endl;
// After we have the pressure, we can compute the velocity. The equation
// reads $MU=-BP+F$, and we solve it by first computing the right hand
// side, and then multiplying it with the object that represents the
// inverse of the mass matrix:
- {
- system_matrix.block(0, 1).vmult(tmp, solution.block(1));
- tmp *= -1;
- tmp += system_rhs.block(0);
- inverse_mass.vmult(solution.block(0), tmp);
- }
+ U = op_M_inv * (F - op_B * P);
}
projects / dealii.git / commitdiff