ASPECT
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
matrix_free_operators.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 #ifndef _aspect_simulator_stokes_matrix_free_operators_h
22 #define _aspect_simulator_stokes_matrix_free_operators_h
23 
24 #include <aspect/global.h>
26 #include <aspect/simulator.h>
27 
28 #include <deal.II/matrix_free/matrix_free.h>
29 #include <deal.II/matrix_free/operators.h>
30 #include <deal.II/matrix_free/fe_evaluation.h>
31 
32 #include <deal.II/multigrid/mg_constrained_dofs.h>
33 #include <deal.II/multigrid/multigrid.h>
34 #include <deal.II/multigrid/mg_transfer_matrix_free.h>
35 #include <deal.II/multigrid/mg_transfer_global_coarsening.templates.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 
46 namespace aspect
47 {
48  namespace internal
49  {
55  namespace ChangeVectorTypes
56  {
57  void import(TrilinosWrappers::MPI::Vector &out,
58  const ::LinearAlgebra::ReadWriteVector<double> &rwv,
59  const VectorOperation::values operation);
60 
62  const ::LinearAlgebra::distributed::Vector<double> &in);
63 
64  void copy(::LinearAlgebra::distributed::Vector<double> &out,
66 
68  const ::LinearAlgebra::distributed::BlockVector<double> &in);
69 
70  void copy(::LinearAlgebra::distributed::BlockVector<double> &out,
72  }
73  }
74 
78  namespace MatrixFreeStokesOperators
79  {
80 
91  template <int dim, typename number>
93  {
98 
103 
108 
113 
119 
124 
131  Table<2, VectorizedArray<number>> viscosity;
132 
137  Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>> strain_rate_table;
138 
144  Table<2, VectorizedArray<number>> newton_factor_wrt_pressure_table;
145 
153  Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
155 
161  Table<2, VectorizedArray<number>> dilation_lhs_term_table;
162 
170  Table<2, VectorizedArray<number>>
172 
177  Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
179 
185  Table<2, Tensor<1, dim, VectorizedArray<number>>> free_surface_stabilization_term_table;
186 
190  std::set<types::boundary_id> free_surface_boundary_indicators;
191 
196  std::size_t
197  memory_consumption() const;
198 
202  void clear();
203  };
204 
208  template <int dim, int degree_v, typename number>
210  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::BlockVector<number>>
211  {
212  public:
213 
217  StokesOperator ();
218 
222  void clear () override;
223 
227  void set_cell_data (const OperatorCellData<dim,number> &data);
228 
234  void compute_diagonal () override;
235 
236  private:
237 
242  void apply_add (::LinearAlgebra::distributed::BlockVector<number> &dst,
243  const ::LinearAlgebra::distributed::BlockVector<number> &src) const override;
244 
248  void local_apply (const ::MatrixFree<dim, number> &data,
249  ::LinearAlgebra::distributed::BlockVector<number> &dst,
250  const ::LinearAlgebra::distributed::BlockVector<number> &src,
251  const std::pair<unsigned int, unsigned int> &cell_range) const;
252 
256  void local_apply_face (const ::MatrixFree<dim, number> &data,
257  ::LinearAlgebra::distributed::BlockVector<number> &dst,
258  const ::LinearAlgebra::distributed::BlockVector<number> &src,
259  const std::pair<unsigned int, unsigned int> &face_range) const;
260 
264  void local_apply_boundary_face (const ::MatrixFree<dim, number> &data,
265  ::LinearAlgebra::distributed::BlockVector<number> &dst,
266  const ::LinearAlgebra::distributed::BlockVector<number> &src,
267  const std::pair<unsigned int, unsigned int> &face_range) const;
268 
273  };
274 
275 
276 
280  template <int dim, int degree_v, typename number>
282  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::BlockVector<number>>
283  {
284  public:
285 
289  BTBlockOperator ();
290 
294  void clear () override;
295 
299  void set_cell_data (const OperatorCellData<dim,number> &data);
300 
306  void compute_diagonal () override;
307 
308  private:
309 
314  void apply_add (::LinearAlgebra::distributed::BlockVector<number> &dst,
315  const ::LinearAlgebra::distributed::BlockVector<number> &src) const override;
316 
320  void local_apply (const ::MatrixFree<dim, number> &data,
321  ::LinearAlgebra::distributed::BlockVector<number> &dst,
322  const ::LinearAlgebra::distributed::BlockVector<number> &src,
323  const std::pair<unsigned int, unsigned int> &cell_range) const;
324 
328  void local_apply_face (const ::MatrixFree<dim, number> &data,
329  ::LinearAlgebra::distributed::BlockVector<number> &dst,
330  const ::LinearAlgebra::distributed::BlockVector<number> &src,
331  const std::pair<unsigned int, unsigned int> &face_range) const;
332 
337  };
338 
339 
340 
344  template <int dim, int degree_p, typename number>
346  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::Vector<number>>
347  {
348  public:
349 
354 
358  void clear () override;
359 
363  void set_cell_data (const OperatorCellData<dim,number> &data);
364 
370  void compute_diagonal () override;
371 
372  private:
373 
378  void apply_add (::LinearAlgebra::distributed::Vector<number> &dst,
379  const ::LinearAlgebra::distributed::Vector<number> &src) const override;
380 
384  void local_apply (const ::MatrixFree<dim, number> &data,
385  ::LinearAlgebra::distributed::Vector<number> &dst,
386  const ::LinearAlgebra::distributed::Vector<number> &src,
387  const std::pair<unsigned int, unsigned int> &cell_range) const;
388 
392  void inner_cell_operation(FEEvaluation<dim,
393  degree_p,
394  degree_p+2,
395  1,
396  number> &pressure) const;
397 
402  };
403 
408  template <int dim, int degree_v, typename number>
410  : public MatrixFreeOperators::Base<dim, ::LinearAlgebra::distributed::Vector<number>>
411  {
412  public:
413 
417  ABlockOperator ();
418 
422  void clear () override;
423 
427  void set_cell_data (const OperatorCellData<dim,number> &data);
428 
434  void compute_diagonal () override;
435 
441  void set_diagonal (const ::LinearAlgebra::distributed::Vector<number> &diag);
442 
443  private:
448  void inner_cell_operation(FEEvaluation<dim,
449  degree_v,
450  degree_v+1,
451  dim,
452  number> &velocity) const;
453 
458  void cell_operation(FEEvaluation<dim,
459  degree_v,
460  degree_v+1,
461  dim,
462  number> &velocity) const;
463 
469  void apply_add (::LinearAlgebra::distributed::Vector<number> &dst,
470  const ::LinearAlgebra::distributed::Vector<number> &src) const override;
471 
475  void local_apply (const ::MatrixFree<dim, number> &data,
476  ::LinearAlgebra::distributed::Vector<number> &dst,
477  const ::LinearAlgebra::distributed::Vector<number> &src,
478  const std::pair<unsigned int, unsigned int> &cell_range) const;
479 
484  };
485  }
486 
487 }
488 
489 #endif
Table< 2, VectorizedArray< number > > dilation_derivative_wrt_pressure_table
const OperatorCellData< dim, number > * cell_data
std::set< types::boundary_id > free_surface_boundary_indicators
::TrilinosWrappers::MPI::BlockVector BlockVector
Definition: global.h:269
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:263
const OperatorCellData< dim, number > * cell_data
Table< 2, SymmetricTensor< 2, dim, VectorizedArray< number > > > strain_rate_table
const OperatorCellData< dim, number > * cell_data
Table< 2, SymmetricTensor< 2, dim, VectorizedArray< number > > > dilation_derivative_wrt_strain_rate_table
Table< 2, VectorizedArray< number > > dilation_lhs_term_table
Table< 2, SymmetricTensor< 2, dim, VectorizedArray< number > > > newton_factor_wrt_strain_rate_table
void copy(TrilinosWrappers::MPI::Vector &out, const ::LinearAlgebra::distributed::Vector< double > &in)
Table< 2, VectorizedArray< number > > newton_factor_wrt_pressure_table