ASPECT
stokes_matrix_free.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2018 - 2024 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 
22 #ifndef _aspect_stokes_matrix_free_h
23 #define _aspect_stokes_matrix_free_h
24 
25 #include <aspect/global.h>
26 
27 #include <aspect/simulator.h>
28 
29 #include <deal.II/matrix_free/matrix_free.h>
30 #include <deal.II/matrix_free/operators.h>
31 #include <deal.II/matrix_free/fe_evaluation.h>
32 
33 #include <deal.II/multigrid/mg_constrained_dofs.h>
34 #include <deal.II/multigrid/multigrid.h>
35 #include <deal.II/multigrid/mg_transfer_matrix_free.h>
36 #include <deal.II/multigrid/mg_transfer_global_coarsening.templates.h>
37 #include <deal.II/multigrid/mg_tools.h>
38 #include <deal.II/multigrid/mg_coarse.h>
39 #include <deal.II/multigrid/mg_smoother.h>
40 #include <deal.II/multigrid/mg_matrix.h>
41 
42 #include <deal.II/lac/vector.h>
43 #include <deal.II/lac/block_vector.h>
44 #include <deal.II/lac/la_parallel_vector.h>
45 #include <deal.II/lac/la_parallel_block_vector.h>
46 
50 using GMGNumberType = double;
51 
52 namespace aspect
53 {
54  namespace internal
55  {
61  namespace ChangeVectorTypes
62  {
63  void import(TrilinosWrappers::MPI::Vector &out,
64  const ::LinearAlgebra::ReadWriteVector<double> &rwv,
65  const VectorOperation::values operation);
66 
68  const ::LinearAlgebra::distributed::Vector<double> &in);
69 
70  void copy(::LinearAlgebra::distributed::Vector<double> &out,
72 
74  const ::LinearAlgebra::distributed::BlockVector<double> &in);
75 
76  void copy(::LinearAlgebra::distributed::BlockVector<double> &out,
78  }
79  }
80 
84  namespace MatrixFreeStokesOperators
85  {
86 
97  template <int dim, typename number>
99  {
104 
109 
114 
120 
125 
132  Table<2, VectorizedArray<number>> viscosity;
133 
138  Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>> strain_rate_table;
139 
145  Table<2, VectorizedArray<number>> newton_factor_wrt_pressure_table;
146 
154  Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
156 
162  Table<2, Tensor<1, dim, VectorizedArray<number>>> free_surface_stabilization_term_table;
163 
167  std::set<types::boundary_id> free_surface_boundary_indicators;
168 
173  std::size_t
174  memory_consumption() const;
175 
179  void clear();
180  };
181 
185  template <int dim, int degree_v, typename number>
187  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::BlockVector<number>>
188  {
189  public:
190 
194  StokesOperator ();
195 
199  void clear () override;
200 
204  void set_cell_data (const OperatorCellData<dim,number> &data);
205 
211  void compute_diagonal () override;
212 
213  private:
214 
219  void apply_add (::LinearAlgebra::distributed::BlockVector<number> &dst,
220  const ::LinearAlgebra::distributed::BlockVector<number> &src) const override;
221 
225  void local_apply (const ::MatrixFree<dim, number> &data,
226  ::LinearAlgebra::distributed::BlockVector<number> &dst,
227  const ::LinearAlgebra::distributed::BlockVector<number> &src,
228  const std::pair<unsigned int, unsigned int> &cell_range) const;
229 
233  void local_apply_face (const ::MatrixFree<dim, number> &data,
234  ::LinearAlgebra::distributed::BlockVector<number> &dst,
235  const ::LinearAlgebra::distributed::BlockVector<number> &src,
236  const std::pair<unsigned int, unsigned int> &face_range) const;
237 
241  void local_apply_boundary_face (const ::MatrixFree<dim, number> &data,
242  ::LinearAlgebra::distributed::BlockVector<number> &dst,
243  const ::LinearAlgebra::distributed::BlockVector<number> &src,
244  const std::pair<unsigned int, unsigned int> &face_range) const;
245 
250  };
251 
255  template <int dim, int degree_p, typename number>
257  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::Vector<number>>
258  {
259  public:
260 
265 
269  void clear () override;
270 
274  void set_cell_data (const OperatorCellData<dim,number> &data);
275 
281  void compute_diagonal () override;
282 
283  private:
284 
289  void apply_add (::LinearAlgebra::distributed::Vector<number> &dst,
290  const ::LinearAlgebra::distributed::Vector<number> &src) const override;
291 
295  void local_apply (const ::MatrixFree<dim, number> &data,
296  ::LinearAlgebra::distributed::Vector<number> &dst,
297  const ::LinearAlgebra::distributed::Vector<number> &src,
298  const std::pair<unsigned int, unsigned int> &cell_range) const;
299 
303  void inner_cell_operation(FEEvaluation<dim,
304  degree_p,
305  degree_p+2,
306  1,
307  number> &pressure) const;
308 
313  };
314 
319  template <int dim, int degree_v, typename number>
321  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::Vector<number>>
322  {
323  public:
324 
328  ABlockOperator ();
329 
333  void clear () override;
334 
338  void set_cell_data (const OperatorCellData<dim,number> &data);
339 
345  void compute_diagonal () override;
346 
352  void set_diagonal (const ::LinearAlgebra::distributed::Vector<number> &diag);
353 
354  private:
359  void inner_cell_operation(FEEvaluation<dim,
360  degree_v,
361  degree_v+1,
362  dim,
363  number> &velocity) const;
364 
369  void cell_operation(FEEvaluation<dim,
370  degree_v,
371  degree_v+1,
372  dim,
373  number> &velocity) const;
374 
380  void apply_add (::LinearAlgebra::distributed::Vector<number> &dst,
381  const ::LinearAlgebra::distributed::Vector<number> &src) const override;
382 
386  void local_apply (const ::MatrixFree<dim, number> &data,
387  ::LinearAlgebra::distributed::Vector<number> &dst,
388  const ::LinearAlgebra::distributed::Vector<number> &src,
389  const std::pair<unsigned int, unsigned int> &cell_range) const;
390 
395  };
396  }
397 
402  template <int dim>
404  {
405  public:
409  virtual ~StokesMatrixFreeHandler() = default;
410 
420  virtual std::pair<double,double> solve(LinearAlgebra::BlockVector &solution_vector) = 0;
421 
426  virtual void setup_dofs()=0;
427 
435  virtual void assemble()=0;
436 
442  virtual void build_preconditioner()=0;
443 
447  static
448  void declare_parameters (ParameterHandler &prm);
449 
454  virtual const DoFHandler<dim> &
455  get_dof_handler_v () const = 0;
456 
461  virtual const DoFHandler<dim> &
462  get_dof_handler_p () const = 0;
463 
468  virtual const DoFHandler<dim> &
469  get_dof_handler_projection () const = 0;
470 
475  virtual const AffineConstraints<double> &
476  get_constraints_v () const = 0;
477 
482  virtual const AffineConstraints<double> &
483  get_constraints_p () const = 0;
484 
489  virtual const MGTransferMF<dim,GMGNumberType> &
490  get_mg_transfer_A () const = 0;
491 
496  virtual const MGTransferMF<dim,GMGNumberType> &
497  get_mg_transfer_S () const = 0;
498 
503  virtual std::size_t get_cell_data_memory_consumption() const = 0;
504  };
505 
516  template <int dim, int velocity_degree>
518  {
519  public:
526  StokesMatrixFreeHandlerImplementation(Simulator<dim> &, ParameterHandler &prm);
527 
531  ~StokesMatrixFreeHandlerImplementation() override = default;
532 
542  std::pair<double,double> solve(LinearAlgebra::BlockVector &solution_vector) override;
543 
548  void setup_dofs() override;
549 
557  void assemble() override;
558 
563  void build_preconditioner() override;
564 
568  static
569  void declare_parameters (ParameterHandler &prm);
570 
575  const DoFHandler<dim> &
576  get_dof_handler_v () const override;
577 
582  const DoFHandler<dim> &
583  get_dof_handler_p () const override;
584 
589  const DoFHandler<dim> &
590  get_dof_handler_projection () const override;
591 
596  const AffineConstraints<double> &
597  get_constraints_v () const override;
598 
603  const AffineConstraints<double> &
604  get_constraints_p () const override;
605 
611  get_mg_transfer_A () const override;
612 
618  get_mg_transfer_S () const override;
619 
620 
625  std::size_t get_cell_data_memory_consumption() const override;
626 
627  private:
631  void parse_parameters (ParameterHandler &prm);
632 
638  void evaluate_material_model();
639 
644  void correct_stokes_rhs();
645 
646 
648 
650 
656 
662 
663  DoFHandler<dim> dof_handler_v;
664  DoFHandler<dim> dof_handler_p;
665  DoFHandler<dim> dof_handler_projection;
666 
667  FESystem<dim> fe_v;
668  FESystem<dim> fe_p;
669  FESystem<dim> fe_projection;
670 
675 
679  MGLevelObject<MatrixFreeStokesOperators::OperatorCellData<dim, GMGNumberType>> level_cell_data;
680 
684 
687 
691 
692  AffineConstraints<double> constraints_v;
693  AffineConstraints<double> constraints_p;
694 
695  MGLevelObject<GMGABlockMatrixType> mg_matrices_A_block;
696  MGLevelObject<GMGSchurComplementMatrixType> mg_matrices_Schur_complement;
697 
698  MGConstrainedDoFs mg_constrained_dofs_A_block;
701 
704 
705  std::vector<std::shared_ptr<MatrixFree<dim,double>>> matrix_free_objects;
706  };
707 }
708 
709 
710 #endif
const OperatorCellData< dim, number > * cell_data
std::set< types::boundary_id > free_surface_boundary_indicators
MGLevelObject< MatrixFreeStokesOperators::OperatorCellData< dim, GMGNumberType > > level_cell_data
::TrilinosWrappers::MPI::BlockVector BlockVector
Definition: global.h:268
MatrixFreeStokesOperators::OperatorCellData< dim, GMGNumberType > active_cell_data
MGLevelObject< GMGSchurComplementMatrixType > mg_matrices_Schur_complement
MGLevelObject< GMGABlockMatrixType > mg_matrices_A_block
const OperatorCellData< dim, number > * cell_data
Table< 2, Tensor< 1, dim, VectorizedArray< number > > > free_surface_stabilization_term_table
::TrilinosWrappers::MPI::Vector Vector
Definition: global.h:262
const OperatorCellData< dim, number > * cell_data
void declare_parameters(ParameterHandler &prm)
Table< 2, SymmetricTensor< 2, dim, VectorizedArray< number > > > strain_rate_table
MGTransferMF< dim, GMGNumberType > mg_transfer_A_block
MGTransferMatrixFree< dim, NumberType > MGTransferMF
Definition: compat.h:45
std::vector< std::shared_ptr< MatrixFree< dim, double > > > matrix_free_objects
Definition: compat.h:59
MGTransferMF< dim, GMGNumberType > mg_transfer_Schur_complement
Table< 2, VectorizedArray< number > > viscosity
Table< 2, SymmetricTensor< 2, dim, VectorizedArray< number > > > newton_factor_wrt_strain_rate_table
SchurComplementMatrixType Schur_complement_block_matrix
void copy(TrilinosWrappers::MPI::Vector &out, const ::LinearAlgebra::distributed::Vector< double > &in)
Table< 2, VectorizedArray< number > > newton_factor_wrt_pressure_table
double GMGNumberType