ASPECT
stokes_matrix_free.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2018 - 2022 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_tools.h>
37 #include <deal.II/multigrid/mg_coarse.h>
38 #include <deal.II/multigrid/mg_smoother.h>
39 #include <deal.II/multigrid/mg_matrix.h>
40 
41 #include <deal.II/lac/vector.h>
42 #include <deal.II/lac/block_vector.h>
43 #include <deal.II/lac/la_parallel_vector.h>
44 #include <deal.II/lac/la_parallel_block_vector.h>
45 
49 using GMGNumberType = double;
50 
51 namespace aspect
52 {
53  using namespace dealii;
54 
55  namespace internal
56  {
57  namespace TangentialBoundaryFunctions
58  {
59 
60  template <int dim>
61  void compute_no_normal_flux_constraints_box (const DoFHandler<dim> &dof,
62  const types::boundary_id bid,
63  const unsigned int first_vector_component,
64  MGConstrainedDoFs &mg_constrained_dofs);
65  template <int dim>
66  void
67  add_constraint(const std::array<types::global_dof_index,dim> &dof_indices,
68  const Tensor<1, dim> &constraining_vector,
69  AffineConstraints<double> &constraints,
70  const double inhomogeneity = 0);
71 
72  template <int dim, int spacedim>
73  void compute_no_normal_flux_constraints_shell(const DoFHandler<dim,spacedim> &dof_handler,
74  const MGConstrainedDoFs &mg_constrained_dofs,
75  const Mapping<dim> &mapping,
76  const unsigned int level,
77  const unsigned int first_vector_component,
78  const std::set<types::boundary_id> &boundary_ids,
79  AffineConstraints<double> &constraints);
80  }
81 
87  namespace ChangeVectorTypes
88  {
89  void import(TrilinosWrappers::MPI::Vector &out,
90  const ::LinearAlgebra::ReadWriteVector<double> &rwv,
91  const VectorOperation::values operation);
92 
94  const ::LinearAlgebra::distributed::Vector<double> &in);
95 
96  void copy(::LinearAlgebra::distributed::Vector<double> &out,
98 
100  const ::LinearAlgebra::distributed::BlockVector<double> &in);
101 
102  void copy(::LinearAlgebra::distributed::BlockVector<double> &out,
104  }
105  }
106 
110  namespace MatrixFreeStokesOperators
111  {
112 
123  template <int dim, typename number>
125  {
130 
135 
140 
146 
151 
158  Table<2, VectorizedArray<number>> viscosity;
159 
164  Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>> strain_rate_table;
165 
171  Table<2, VectorizedArray<number>> newton_factor_wrt_pressure_table;
172 
180  Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
182 
188  Table<2, Tensor<1, dim, VectorizedArray<number>>> free_surface_stabilization_term_table;
189 
193  std::set<types::boundary_id> free_surface_boundary_indicators;
194 
199  std::size_t
200  memory_consumption() const;
201 
205  void clear();
206  };
207 
211  template <int dim, int degree_v, typename number>
213  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::BlockVector<number>>
214  {
215  public:
216 
220  StokesOperator ();
221 
225  void clear () override;
226 
230  void set_cell_data (const OperatorCellData<dim,number> &data);
231 
237  void compute_diagonal () override;
238 
239  private:
240 
245  void apply_add (::LinearAlgebra::distributed::BlockVector<number> &dst,
246  const ::LinearAlgebra::distributed::BlockVector<number> &src) const override;
247 
251  void local_apply (const ::MatrixFree<dim, number> &data,
252  ::LinearAlgebra::distributed::BlockVector<number> &dst,
253  const ::LinearAlgebra::distributed::BlockVector<number> &src,
254  const std::pair<unsigned int, unsigned int> &cell_range) const;
255 
259  void local_apply_face (const ::MatrixFree<dim, number> &data,
260  ::LinearAlgebra::distributed::BlockVector<number> &dst,
261  const ::LinearAlgebra::distributed::BlockVector<number> &src,
262  const std::pair<unsigned int, unsigned int> &face_range) const;
263 
267  void local_apply_boundary_face (const ::MatrixFree<dim, number> &data,
268  ::LinearAlgebra::distributed::BlockVector<number> &dst,
269  const ::LinearAlgebra::distributed::BlockVector<number> &src,
270  const std::pair<unsigned int, unsigned int> &face_range) const;
271 
276  };
277 
281  template <int dim, int degree_p, typename number>
283  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::Vector<number>>
284  {
285  public:
286 
291 
295  void clear () override;
296 
300  void set_cell_data (const OperatorCellData<dim,number> &data);
301 
307  void compute_diagonal () override;
308 
309  private:
310 
315  void apply_add (::LinearAlgebra::distributed::Vector<number> &dst,
316  const ::LinearAlgebra::distributed::Vector<number> &src) const override;
317 
321  void local_apply (const ::MatrixFree<dim, number> &data,
322  ::LinearAlgebra::distributed::Vector<number> &dst,
323  const ::LinearAlgebra::distributed::Vector<number> &src,
324  const std::pair<unsigned int, unsigned int> &cell_range) const;
325 
326 
330  void local_compute_diagonal (const MatrixFree<dim,number> &data,
331  ::LinearAlgebra::distributed::Vector<number> &dst,
332  const unsigned int &dummy,
333  const std::pair<unsigned int,unsigned int> &cell_range) const;
334 
339  };
340 
345  template <int dim, int degree_v, typename number>
347  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::Vector<number>>
348  {
349  public:
350 
354  ABlockOperator ();
355 
359  void clear () override;
360 
364  void set_cell_data (const OperatorCellData<dim,number> &data);
365 
371  void compute_diagonal () override;
372 
378  void set_diagonal (const ::LinearAlgebra::distributed::Vector<number> &diag);
379 
380  private:
381 
386  void apply_add (::LinearAlgebra::distributed::Vector<number> &dst,
387  const ::LinearAlgebra::distributed::Vector<number> &src) const override;
388 
392  void local_apply (const ::MatrixFree<dim, number> &data,
393  ::LinearAlgebra::distributed::Vector<number> &dst,
394  const ::LinearAlgebra::distributed::Vector<number> &src,
395  const std::pair<unsigned int, unsigned int> &cell_range) const;
396 
400  void local_compute_diagonal (const MatrixFree<dim,number> &data,
401  ::LinearAlgebra::distributed::Vector<number> &dst,
402  const unsigned int &dummy,
403  const std::pair<unsigned int,unsigned int> &cell_range) const;
404 
409  };
410  }
411 
416  template<int dim>
418  {
419  public:
423  virtual ~StokesMatrixFreeHandler() = default;
424 
429  virtual std::pair<double,double> solve()=0;
430 
435  virtual void setup_dofs()=0;
436 
444  virtual void assemble()=0;
445 
451  virtual void build_preconditioner()=0;
452 
456  static
457  void declare_parameters (ParameterHandler &prm);
458 
463  virtual const DoFHandler<dim> &
464  get_dof_handler_v () const = 0;
465 
470  virtual const DoFHandler<dim> &
471  get_dof_handler_p () const = 0;
472 
477  virtual const DoFHandler<dim> &
478  get_dof_handler_projection () const = 0;
479 
484  virtual const AffineConstraints<double> &
485  get_constraints_v () const = 0;
486 
491  virtual const AffineConstraints<double> &
492  get_constraints_p () const = 0;
493 
498  virtual const MGTransferMatrixFree<dim,GMGNumberType> &
499  get_mg_transfer_A () const = 0;
500 
505  virtual const MGTransferMatrixFree<dim,GMGNumberType> &
506  get_mg_transfer_S () const = 0;
507 
512  virtual std::size_t get_cell_data_memory_consumption() const = 0;
513  };
514 
525  template<int dim, int velocity_degree>
527  {
528  public:
535  StokesMatrixFreeHandlerImplementation(Simulator<dim> &, ParameterHandler &prm);
536 
540  ~StokesMatrixFreeHandlerImplementation() override = default;
541 
546  std::pair<double,double> solve() override;
547 
552  void setup_dofs() override;
553 
561  virtual void assemble() override;
562 
567  void build_preconditioner() override;
568 
572  static
573  void declare_parameters (ParameterHandler &prm);
574 
579  const DoFHandler<dim> &
580  get_dof_handler_v () const override;
581 
586  const DoFHandler<dim> &
587  get_dof_handler_p () const override;
588 
593  const DoFHandler<dim> &
594  get_dof_handler_projection () const override;
595 
600  const AffineConstraints<double> &
601  get_constraints_v () const override;
602 
607  const AffineConstraints<double> &
608  get_constraints_p () const override;
609 
614  const MGTransferMatrixFree<dim,GMGNumberType> &
615  get_mg_transfer_A () const override;
616 
621  const MGTransferMatrixFree<dim,GMGNumberType> &
622  get_mg_transfer_S () const override;
623 
624 
629  std::size_t get_cell_data_memory_consumption() const override;
630 
631  private:
635  void parse_parameters (ParameterHandler &prm);
636 
642  void evaluate_material_model();
643 
648  void correct_stokes_rhs();
649 
650 
652 
654 
660 
666 
667  DoFHandler<dim> dof_handler_v;
668  DoFHandler<dim> dof_handler_p;
669  DoFHandler<dim> dof_handler_projection;
670 
671  FESystem<dim> fe_v;
672  FESystem<dim> fe_p;
673  FESystem<dim> fe_projection;
674 
679 
683  MGLevelObject<MatrixFreeStokesOperators::OperatorCellData<dim, GMGNumberType>> level_cell_data;
684 
685  // This variable is needed only in the setup in both evaluate_material_model()
686  // and build_preconditioner(). It will be deleted after the last use.
687  MGLevelObject<::LinearAlgebra::distributed::Vector<GMGNumberType>> level_viscosity_vector;
688 
692 
695 
699 
700  AffineConstraints<double> constraints_v;
701  AffineConstraints<double> constraints_p;
702 
703  MGLevelObject<GMGABlockMatrixType> mg_matrices_A_block;
704  MGLevelObject<GMGSchurComplementMatrixType> mg_matrices_Schur_complement;
705 
706  MGConstrainedDoFs mg_constrained_dofs_A_block;
709 
710  MGTransferMatrixFree<dim,GMGNumberType> mg_transfer_A_block;
711  MGTransferMatrixFree<dim,GMGNumberType> mg_transfer_Schur_complement;
712 
713  std::vector<std::shared_ptr<MatrixFree<dim,double>>> matrix_free_objects;
714  };
715 }
716 
717 
718 #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:229
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
void compute_no_normal_flux_constraints_box(const DoFHandler< dim > &dof, const types::boundary_id bid, const unsigned int first_vector_component, MGConstrainedDoFs &mg_constrained_dofs)
Table< 2, Tensor< 1, dim, VectorizedArray< number > > > free_surface_stabilization_term_table
void copy(::LinearAlgebra::distributed::BlockVector< double > &out, const TrilinosWrappers::MPI::BlockVector &in)
::TrilinosWrappers::MPI::Vector Vector
Definition: global.h:223
MGTransferMatrixFree< dim, GMGNumberType > mg_transfer_Schur_complement
const OperatorCellData< dim, number > * cell_data
void compute_no_normal_flux_constraints_shell(const DoFHandler< dim, spacedim > &dof_handler, const MGConstrainedDoFs &mg_constrained_dofs, const Mapping< dim > &mapping, const unsigned int level, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, AffineConstraints< double > &constraints)
MGTransferMatrixFree< dim, GMGNumberType > mg_transfer_A_block
void declare_parameters(ParameterHandler &prm)
Table< 2, SymmetricTensor< 2, dim, VectorizedArray< number > > > strain_rate_table
void add_constraint(const std::array< types::global_dof_index, dim > &dof_indices, const Tensor< 1, dim > &constraining_vector, AffineConstraints< double > &constraints, const double inhomogeneity=0)
std::vector< std::shared_ptr< MatrixFree< dim, double > > > matrix_free_objects
Table< 2, VectorizedArray< number > > viscosity
Table< 2, SymmetricTensor< 2, dim, VectorizedArray< number > > > newton_factor_wrt_strain_rate_table
SchurComplementMatrixType Schur_complement_block_matrix
Table< 2, VectorizedArray< number > > newton_factor_wrt_pressure_table
MGLevelObject<::LinearAlgebra::distributed::Vector< GMGNumberType > > level_viscosity_vector
double GMGNumberType