stokes_matrix_free.h

Go to the documentation of this file.

/*
  Copyright (C) 2011 - 2025 by the authors of the ASPECT code.

  This file is part of ASPECT.

  ASPECT is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2, or (at your option)
  any later version.

  ASPECT is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with ASPECT; see the file LICENSE.  If not see
  <http://www.gnu.org/licenses/>.
*/

#ifndef _aspect_simulator_solver_stokes_matrix_free_h
#define _aspect_simulator_solver_stokes_matrix_free_h

#include <aspect/simulator/solver/matrix_free_operators.h>
#include <aspect/global.h>
#include <aspect/parameters.h>
#include <aspect/simulator/solver/interface.h>
#include <aspect/utilities.h>

#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_transfer_global_coarsening.templates.h>

#include <deal.II/lac/solver_bicgstab.h>

namespace aspect
{
  using GMGNumberType = double;

  template <int dim>
  class StokesMatrixFreeHandler: public StokesSolver::Interface<dim>
  {
    public:
      virtual void setup_dofs()=0;

      virtual void assemble()=0;

      virtual void build_preconditioner()=0;

      static
      void declare_parameters (ParameterHandler &prm);

      virtual std::size_t get_dof_handler_memory_consumption() const = 0;

      virtual std::size_t get_mg_transfer_memory_consumption() const = 0;

      virtual std::size_t get_constraint_memory_consumption() const = 0;

      virtual std::size_t get_cell_data_memory_consumption() const = 0;
  };

  template <int dim>
  std::unique_ptr<StokesMatrixFreeHandler<dim>> create_matrix_free_solver(Simulator<dim> &simulator, const Parameters<dim> &parameters);


  namespace internal
  {

    template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType,
              class SchurComplementMatrixType, class ABlockPreconditionerType,
              class SchurComplementPreconditionerType>
    class BlockSchurGMGPreconditioner : public
#if DEAL_II_VERSION_GTE(9,7,0)
      EnableObserverPointer
#else
      Subscriptor
#endif
    {
      public:
        BlockSchurGMGPreconditioner (const StokesMatrixType                  &Stokes_matrix,
                                     const ABlockMatrixType                  &A_block,
                                     const BTBlockOperatorType               &BT_block,
                                     const SchurComplementMatrixType         &Schur_complement_block,
                                     const ABlockPreconditionerType          &A_block_preconditioner,
                                     const SchurComplementPreconditionerType &Schur_complement_preconditioner,
                                     const bool                               do_solve_A,
                                     const bool                               do_solve_Schur_complement,
                                     const bool                               A_block_is_symmetric,
                                     const double                             A_block_tolerance,
                                     const double                             Schur_complement_tolerance);

        void vmult (::LinearAlgebra::distributed::BlockVector<double>       &dst,
                    const ::LinearAlgebra::distributed::BlockVector<double> &src) const;

        unsigned int n_iterations_A_block() const;
        unsigned int n_iterations_Schur_complement() const;


      private:
        const StokesMatrixType                  &stokes_matrix;
        const ABlockMatrixType                  &A_block;
        const BTBlockOperatorType               &BT_block;
        const SchurComplementMatrixType         &Schur_complement_block;
        const ABlockPreconditionerType          &A_block_preconditioner;
        const SchurComplementPreconditionerType &Schur_complement_preconditioner;

        const bool                                                      do_solve_A;
        const bool                                                      do_solve_Schur_complement;
        const bool                                                      A_block_is_symmetric;
        mutable unsigned int                                            n_iterations_A_;
        mutable unsigned int                                            n_iterations_Schur_complement_;
        const double                                                    A_block_tolerance;
        const double                                                    Schur_complement_tolerance;
        mutable ::LinearAlgebra::distributed::BlockVector<double> utmp;
    };

    template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType, class SchurComplementMatrixType,
              class ABlockPreconditionerType, class SchurComplementPreconditionerType>
    BlockSchurGMGPreconditioner<StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType,
                                ABlockPreconditionerType, SchurComplementPreconditionerType>::
                                BlockSchurGMGPreconditioner (const StokesMatrixType                  &Stokes_matrix,
                                                             const ABlockMatrixType                  &A_block,
                                                             const BTBlockOperatorType               &BT_block,
                                                             const SchurComplementMatrixType         &Schur_complement_block,
                                                             const ABlockPreconditionerType          &A_block_preconditioner,
                                                             const SchurComplementPreconditionerType &Schur_complement_preconditioner,
                                                             const bool                               do_solve_A,
                                                             const bool                               do_solve_Schur_complement,
                                                             const bool                               A_block_symmetric,
                                                             const double                             A_block_tolerance,
                                                             const double                             Schur_complement_tolerance)
                                  :
                                  stokes_matrix                   (Stokes_matrix),
                                  A_block                         (A_block),
                                  BT_block                        (BT_block),
                                  Schur_complement_block          (Schur_complement_block),
                                  A_block_preconditioner          (A_block_preconditioner),
                                  Schur_complement_preconditioner (Schur_complement_preconditioner),
                                  do_solve_A                      (do_solve_A),
                                  do_solve_Schur_complement       (do_solve_Schur_complement),
                                  A_block_is_symmetric            (A_block_symmetric),
                                  n_iterations_A_                 (0),
                                  n_iterations_Schur_complement_  (0),
                                  A_block_tolerance               (A_block_tolerance),
                                  Schur_complement_tolerance      (Schur_complement_tolerance)
    {}

    template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType,
              class SchurComplementMatrixType, class ABlockPreconditionerType,
              class SchurComplementPreconditionerType>
    unsigned int
    BlockSchurGMGPreconditioner<StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType,
                                ABlockPreconditionerType, SchurComplementPreconditionerType>::
                                n_iterations_A_block() const
    {
      return n_iterations_A_;
    }

    template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType, class SchurComplementMatrixType,
              class ABlockPreconditionerType, class SchurComplementPreconditionerType>
    unsigned int
    BlockSchurGMGPreconditioner<StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType,
                                ABlockPreconditionerType, SchurComplementPreconditionerType>::
                                n_iterations_Schur_complement() const
    {
      return n_iterations_Schur_complement_;
    }

    template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType, class SchurComplementMatrixType,
              class ABlockPreconditionerType, class SchurComplementPreconditionerType>
    void
    BlockSchurGMGPreconditioner<StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType,
                                ABlockPreconditionerType, SchurComplementPreconditionerType>::
                                vmult (::LinearAlgebra::distributed::BlockVector<double>       &dst,
                                       const ::LinearAlgebra::distributed::BlockVector<double>  &src) const
    {
      if (utmp.size()==0)
        utmp.reinit(src);

      // This needs to be done explicitly, as GMRES does not
      // initialize the data of the vector dst before calling
      // us. Otherwise we might use random data as our initial guess.
      dst = 0.0;

      // either solve with the Schur complement matrix (if do_solve_Schur_complement==true)
      // or just apply one preconditioner sweep (for the first few
      // iterations of our two-stage outer GMRES iteration)
      if (do_solve_Schur_complement)
        {
          // first solve with the bottom right block, which we have built
          // as a mass matrix with the inverse of the viscosity
          SolverControl solver_control(100, src.block(1).l2_norm() * Schur_complement_tolerance,true);

          SolverCG<::LinearAlgebra::distributed::Vector<double>> solver(solver_control);
          // Trilinos reports a breakdown
          // in case src=dst=0, even
          // though it should return
          // convergence without
          // iterating. We simply skip
          // solving in this case.
          if (src.block(1).l2_norm() > 1e-50)
            {
              try
                {
                  solver.solve(Schur_complement_block,
                               dst.block(1), src.block(1),
                               Schur_complement_preconditioner);
                  n_iterations_Schur_complement_ += solver_control.last_step();
                }
              // if the solver fails, report the error from processor 0 with some additional
              // information about its location, and throw a quiet exception on all other
              // processors
              catch (const std::exception &exc)
                {
                  Utilities::throw_linear_solver_failure_exception("iterative (bottom right) solver",
                                                                   "BlockSchurGMGPreconditioner::vmult",
                                                                   std::vector<SolverControl> {solver_control},
                                                                   exc,
                                                                   src.block(0).get_mpi_communicator());
                }
            }
        }
      else
        {
          Schur_complement_preconditioner.vmult(dst.block(1),src.block(1));
          n_iterations_Schur_complement_ += 1;
        }

      dst.block(1) *= -1.0;

      {
        BT_block.vmult(utmp,dst);
        utmp.block(0) *= -1.0;
        utmp.block(0) += src.block(0);
      }

      // now either solve with the top left block (if do_solve_A==true)
      // or just apply one preconditioner sweep (for the first few
      // iterations of our two-stage outer GMRES iteration)
      if (do_solve_A == true)
        {
          SolverControl solver_control(1000, utmp.block(0).l2_norm() * A_block_tolerance);
          PrimitiveVectorMemory<::LinearAlgebra::distributed::Vector<double>> mem;

          try
            {
              if (A_block_is_symmetric)
                {
                  SolverCG<::LinearAlgebra::distributed::Vector<double>> solver(solver_control,mem);
                  solver.solve(A_block, dst.block(0), utmp.block(0),
                               A_block_preconditioner);
                }
              else
                {
                  // Use BiCGStab for non-symmetric matrices.
                  // BiCGStab can also solve indefinite systems if necessary.
                  // Do not compute the exact residual, as this
                  // is more expensive, and we only need an approximate solution.
                  SolverBicgstab<::LinearAlgebra::distributed::Vector<double>>
                  solver(solver_control,
                         mem,
                         SolverBicgstab<::LinearAlgebra::distributed::Vector<double>>::AdditionalData(/*exact_residual=*/ false));
                  solver.solve(A_block, dst.block(0), utmp.block(0),
                               A_block_preconditioner);
                }

              n_iterations_A_ += solver_control.last_step();
            }
          // if the solver fails, report the error from processor 0 with some additional
          // information about its location, and throw a quiet exception on all other
          // processors
          catch (const std::exception &exc)
            {
              Utilities::throw_linear_solver_failure_exception("iterative (top left) solver",
                                                               "BlockSchurGMGPreconditioner::vmult",
                                                               std::vector<SolverControl> {solver_control},
                                                               exc,
                                                               src.block(0).get_mpi_communicator());
            }
        }
      else
        {
          A_block_preconditioner.vmult (dst.block(0), utmp.block(0));
          n_iterations_A_ += 1;
        }
    }
  }
}

#endif