ASPECT
stokes_matrix_free.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2025 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef _aspect_simulator_solver_stokes_matrix_free_h
22 #define _aspect_simulator_solver_stokes_matrix_free_h
23 
25 #include <aspect/global.h>
26 #include <aspect/parameters.h>
28 #include <aspect/utilities.h>
29 
30 #include <deal.II/multigrid/mg_transfer_matrix_free.h>
31 #include <deal.II/multigrid/mg_transfer_global_coarsening.templates.h>
32 
33 #include <deal.II/lac/solver_bicgstab.h>
34 
35 namespace aspect
36 {
40  using GMGNumberType = double;
41 
46  template <int dim>
48  {
49  public:
54  virtual void setup_dofs()=0;
55 
63  virtual void assemble()=0;
64 
70  virtual void build_preconditioner()=0;
71 
75  static
76  void declare_parameters (ParameterHandler &prm);
77 
81  virtual std::size_t get_dof_handler_memory_consumption() const = 0;
82 
86  virtual std::size_t get_mg_transfer_memory_consumption() const = 0;
87 
91  virtual std::size_t get_constraint_memory_consumption() const = 0;
92 
97  virtual std::size_t get_cell_data_memory_consumption() const = 0;
98  };
99 
103  template <int dim>
104  std::unique_ptr<StokesMatrixFreeHandler<dim>> create_matrix_free_solver(Simulator<dim> &simulator, const Parameters<dim> &parameters);
105 
106 
107  namespace internal
108  {
109 
113  template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType,
114  class SchurComplementMatrixType, class ABlockPreconditionerType,
115  class SchurComplementPreconditionerType>
117  {
118  public:
138  BlockSchurGMGPreconditioner (const StokesMatrixType &Stokes_matrix,
139  const ABlockMatrixType &A_block,
140  const BTBlockOperatorType &BT_block,
141  const SchurComplementMatrixType &Schur_complement_block,
142  const ABlockPreconditionerType &A_block_preconditioner,
143  const SchurComplementPreconditionerType &Schur_complement_preconditioner,
144  const bool do_solve_A,
145  const bool do_solve_Schur_complement,
146  const bool A_block_is_symmetric,
147  const double A_block_tolerance,
148  const double Schur_complement_tolerance);
149 
153  void vmult (::LinearAlgebra::distributed::BlockVector<double> &dst,
154  const ::LinearAlgebra::distributed::BlockVector<double> &src) const;
155 
156  unsigned int n_iterations_A_block() const;
157  unsigned int n_iterations_Schur_complement() const;
158 
159 
160  private:
164  const StokesMatrixType &stokes_matrix;
165  const ABlockMatrixType &A_block;
166  const BTBlockOperatorType &BT_block;
167  const SchurComplementMatrixType &Schur_complement_block;
168  const ABlockPreconditionerType &A_block_preconditioner;
169  const SchurComplementPreconditionerType &Schur_complement_preconditioner;
170 
175  const bool do_solve_A;
178  mutable unsigned int n_iterations_A_;
179  mutable unsigned int n_iterations_Schur_complement_;
180  const double A_block_tolerance;
182  mutable ::LinearAlgebra::distributed::BlockVector<double> utmp;
183  };
184 
185  template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType, class SchurComplementMatrixType,
186  class ABlockPreconditionerType, class SchurComplementPreconditionerType>
187  BlockSchurGMGPreconditioner<StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType,
188  ABlockPreconditionerType, SchurComplementPreconditionerType>::
189  BlockSchurGMGPreconditioner (const StokesMatrixType &Stokes_matrix,
190  const ABlockMatrixType &A_block,
191  const BTBlockOperatorType &BT_block,
192  const SchurComplementMatrixType &Schur_complement_block,
193  const ABlockPreconditionerType &A_block_preconditioner,
194  const SchurComplementPreconditionerType &Schur_complement_preconditioner,
195  const bool do_solve_A,
196  const bool do_solve_Schur_complement,
197  const bool A_block_symmetric,
198  const double A_block_tolerance,
199  const double Schur_complement_tolerance)
200  :
201  stokes_matrix (Stokes_matrix),
202  A_block (A_block),
203  BT_block (BT_block),
204  Schur_complement_block (Schur_complement_block),
205  A_block_preconditioner (A_block_preconditioner),
206  Schur_complement_preconditioner (Schur_complement_preconditioner),
207  do_solve_A (do_solve_A),
208  do_solve_Schur_complement (do_solve_Schur_complement),
209  A_block_is_symmetric (A_block_symmetric),
210  n_iterations_A_ (0),
211  n_iterations_Schur_complement_ (0),
212  A_block_tolerance (A_block_tolerance),
213  Schur_complement_tolerance (Schur_complement_tolerance)
214  {}
215 
216  template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType,
217  class SchurComplementMatrixType, class ABlockPreconditionerType,
218  class SchurComplementPreconditionerType>
219  unsigned int
220  BlockSchurGMGPreconditioner<StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType,
221  ABlockPreconditionerType, SchurComplementPreconditionerType>::
223  {
224  return n_iterations_A_;
225  }
226 
227  template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType, class SchurComplementMatrixType,
228  class ABlockPreconditionerType, class SchurComplementPreconditionerType>
229  unsigned int
230  BlockSchurGMGPreconditioner<StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType,
231  ABlockPreconditionerType, SchurComplementPreconditionerType>::
233  {
235  }
236 
237  template <class StokesMatrixType, class ABlockMatrixType, class BTBlockOperatorType, class SchurComplementMatrixType,
238  class ABlockPreconditionerType, class SchurComplementPreconditionerType>
239  void
240  BlockSchurGMGPreconditioner<StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType,
241  ABlockPreconditionerType, SchurComplementPreconditionerType>::
242  vmult (::LinearAlgebra::distributed::BlockVector<double> &dst,
243  const ::LinearAlgebra::distributed::BlockVector<double> &src) const
244  {
245  if (utmp.size()==0)
246  utmp.reinit(src);
247 
248  // This needs to be done explicitly, as GMRES does not
249  // initialize the data of the vector dst before calling
250  // us. Otherwise we might use random data as our initial guess.
251  dst = 0.0;
252 
253  // either solve with the Schur complement matrix (if do_solve_Schur_complement==true)
254  // or just apply one preconditioner sweep (for the first few
255  // iterations of our two-stage outer GMRES iteration)
257  {
258  // first solve with the bottom right block, which we have built
259  // as a mass matrix with the inverse of the viscosity
260  SolverControl solver_control(100, src.block(1).l2_norm() * Schur_complement_tolerance,true);
261 
262  SolverCG<::LinearAlgebra::distributed::Vector<double>> solver(solver_control);
263  // Trilinos reports a breakdown
264  // in case src=dst=0, even
265  // though it should return
266  // convergence without
267  // iterating. We simply skip
268  // solving in this case.
269  if (src.block(1).l2_norm() > 1e-50)
270  {
271  try
272  {
273  solver.solve(Schur_complement_block,
274  dst.block(1), src.block(1),
276  n_iterations_Schur_complement_ += solver_control.last_step();
277  }
278  // if the solver fails, report the error from processor 0 with some additional
279  // information about its location, and throw a quiet exception on all other
280  // processors
281  catch (const std::exception &exc)
282  {
283  Utilities::throw_linear_solver_failure_exception("iterative (bottom right) solver",
284  "BlockSchurGMGPreconditioner::vmult",
285  std::vector<SolverControl> {solver_control},
286  exc,
287  src.block(0).get_mpi_communicator());
288  }
289  }
290  }
291  else
292  {
293  Schur_complement_preconditioner.vmult(dst.block(1),src.block(1));
295  }
296 
297  dst.block(1) *= -1.0;
298 
299  {
300  BT_block.vmult(utmp,dst);
301  utmp.block(0) *= -1.0;
302  utmp.block(0) += src.block(0);
303  }
304 
305  // now either solve with the top left block (if do_solve_A==true)
306  // or just apply one preconditioner sweep (for the first few
307  // iterations of our two-stage outer GMRES iteration)
308  if (do_solve_A == true)
309  {
310  SolverControl solver_control(1000, utmp.block(0).l2_norm() * A_block_tolerance);
311  PrimitiveVectorMemory<::LinearAlgebra::distributed::Vector<double>> mem;
312 
313  try
314  {
316  {
317  SolverCG<::LinearAlgebra::distributed::Vector<double>> solver(solver_control,mem);
318  solver.solve(A_block, dst.block(0), utmp.block(0),
320  }
321  else
322  {
323  // Use BiCGStab for non-symmetric matrices.
324  // BiCGStab can also solve indefinite systems if necessary.
325  // Do not compute the exact residual, as this
326  // is more expensive, and we only need an approximate solution.
327  SolverBicgstab<::LinearAlgebra::distributed::Vector<double>>
328  solver(solver_control,
329  mem,
330  SolverBicgstab<::LinearAlgebra::distributed::Vector<double>>::AdditionalData(/*exact_residual=*/ false));
331  solver.solve(A_block, dst.block(0), utmp.block(0),
333  }
334 
335  n_iterations_A_ += solver_control.last_step();
336  }
337  // if the solver fails, report the error from processor 0 with some additional
338  // information about its location, and throw a quiet exception on all other
339  // processors
340  catch (const std::exception &exc)
341  {
342  Utilities::throw_linear_solver_failure_exception("iterative (top left) solver",
343  "BlockSchurGMGPreconditioner::vmult",
344  std::vector<SolverControl> {solver_control},
345  exc,
346  src.block(0).get_mpi_communicator());
347  }
348  }
349  else
350  {
351  A_block_preconditioner.vmult (dst.block(0), utmp.block(0));
352  n_iterations_A_ += 1;
353  }
354  }
355  }
356 }
357 
358 #endif
virtual std::size_t get_dof_handler_memory_consumption() const =0
const SchurComplementMatrixType & Schur_complement_block
const Simulator< dim > * simulator
void throw_linear_solver_failure_exception(const std::string &solver_name, const std::string &function_name, const std::vector< SolverControl > &solver_controls, const std::exception &exc, const MPI_Comm mpi_communicator, const std::string &output_filename="")
const SchurComplementPreconditionerType & Schur_complement_preconditioner
virtual std::size_t get_mg_transfer_memory_consumption() const =0
void vmult(::LinearAlgebra::distributed::BlockVector< double > &dst, const ::LinearAlgebra::distributed::BlockVector< double > &src) const
virtual std::size_t get_cell_data_memory_consumption() const =0
std::unique_ptr< StokesMatrixFreeHandler< dim > > create_matrix_free_solver(Simulator< dim > &simulator, const Parameters< dim > &parameters)
double GMGNumberType
mutable ::LinearAlgebra::distributed::BlockVector< double > utmp
virtual void build_preconditioner()=0
const ABlockPreconditionerType & A_block_preconditioner
static void declare_parameters(ParameterHandler &prm)
virtual std::size_t get_constraint_memory_consumption() const =0