Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_precondition_ml.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
18 
19 #ifdef DEAL_II_WITH_TRILINOS
20 
23 # include <deal.II/lac/vector.h>
24 
25 # include <Epetra_MultiVector.h>
26 # include <Ifpack.h>
27 # include <Ifpack_Chebyshev.h>
28 # include <Teuchos_ParameterList.hpp>
29 # include <Teuchos_RCP.hpp>
30 # include <ml_MultiLevelPreconditioner.h>
31 # include <ml_include.h>
32 
34 
35 namespace TrilinosWrappers
36 {
37  /* -------------------------- PreconditionAMG -------------------------- */
38 
40  const bool elliptic,
41  const bool higher_order_elements,
42  const unsigned int n_cycles,
43  const bool w_cycle,
44  const double aggregation_threshold,
45  const std::vector<std::vector<bool>> &constant_modes,
46  const unsigned int smoother_sweeps,
47  const unsigned int smoother_overlap,
48  const bool output_details,
49  const char * smoother_type,
50  const char * coarse_type)
51  : elliptic(elliptic)
52  , higher_order_elements(higher_order_elements)
53  , n_cycles(n_cycles)
54  , w_cycle(w_cycle)
55  , aggregation_threshold(aggregation_threshold)
56  , constant_modes(constant_modes)
57  , smoother_sweeps(smoother_sweeps)
58  , smoother_overlap(smoother_overlap)
59  , output_details(output_details)
60  , smoother_type(smoother_type)
61  , coarse_type(coarse_type)
62  {}
63 
64 
65 
66  void
68  Teuchos::ParameterList & parameter_list,
69  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
70  const Epetra_RowMatrix & matrix) const
71  {
72  if (elliptic == true)
73  {
74  ML_Epetra::SetDefaults("SA", parameter_list);
75 
76  // uncoupled mode can give a lot of warnings or even fail when there
77  // are too many entries per row and aggreggation gets complicated, but
78  // MIS does not work if too few elements are located on one
79  // processor. work around these warnings by choosing the different
80  // strategies in different situations: for low order, always use the
81  // standard choice uncoupled. if higher order, right now we also just
82  // use Uncoupled, but we should be aware that maybe MIS might be
83  // needed
85  parameter_list.set("aggregation: type", "Uncoupled");
86  }
87  else
88  {
89  ML_Epetra::SetDefaults("NSSA", parameter_list);
90  parameter_list.set("aggregation: type", "Uncoupled");
91  parameter_list.set("aggregation: block scaling", true);
92  }
93 
94  parameter_list.set("smoother: type", smoother_type);
95  parameter_list.set("coarse: type", coarse_type);
96 
97  // Force re-initialization of the random seed to make ML deterministic
98  // (only supported in trilinos >12.2):
99 # if DEAL_II_TRILINOS_VERSION_GTE(12, 4, 0)
100  parameter_list.set("initialize random seed", true);
101 # endif
102 
103  parameter_list.set("smoother: sweeps", static_cast<int>(smoother_sweeps));
104  parameter_list.set("cycle applications", static_cast<int>(n_cycles));
105  if (w_cycle == true)
106  parameter_list.set("prec type", "MGW");
107  else
108  parameter_list.set("prec type", "MGV");
109 
110  parameter_list.set("smoother: Chebyshev alpha", 10.);
111  parameter_list.set("smoother: ifpack overlap",
112  static_cast<int>(smoother_overlap));
113  parameter_list.set("aggregation: threshold", aggregation_threshold);
114  parameter_list.set("coarse: max size", 2000);
115 
116  if (output_details)
117  parameter_list.set("ML output", 10);
118  else
119  parameter_list.set("ML output", 0);
120 
121  set_operator_null_space(parameter_list, distributed_constant_modes, matrix);
122  }
123 
124 
125 
126  void
128  Teuchos::ParameterList & parameter_list,
129  std::unique_ptr<Epetra_MultiVector> &ptr_distributed_constant_modes,
130  const Epetra_RowMatrix & matrix) const
131  {
132  const Epetra_Map &domain_map = matrix.OperatorDomainMap();
133 
134  const size_type constant_modes_dimension = constant_modes.size();
135  ptr_distributed_constant_modes = std::make_unique<Epetra_MultiVector>(
136  domain_map, constant_modes_dimension > 0 ? constant_modes_dimension : 1);
137  Assert(ptr_distributed_constant_modes, ExcNotInitialized());
138  Epetra_MultiVector &distributed_constant_modes =
139  *ptr_distributed_constant_modes;
140 
141  if (constant_modes_dimension > 0)
142  {
143  const size_type global_size = TrilinosWrappers::n_global_rows(matrix);
144  (void)global_length; // work around compiler warning about unused
145  // function in release mode
146  Assert(global_size ==
147  static_cast<size_type>(
148  TrilinosWrappers::global_length(distributed_constant_modes)),
149  ExcDimensionMismatch(global_size,
151  distributed_constant_modes)));
152  const bool constant_modes_are_global =
153  constant_modes[0].size() == global_size;
154  const size_type my_size = domain_map.NumMyElements();
155 
156  // Reshape null space as a contiguous vector of doubles so that
157  // Trilinos can read from it.
158  const size_type expected_mode_size =
159  constant_modes_are_global ? global_size : my_size;
160  for (size_type d = 0; d < constant_modes_dimension; ++d)
161  {
162  Assert(constant_modes[d].size() == expected_mode_size,
164  expected_mode_size));
165  for (size_type row = 0; row < my_size; ++row)
166  {
167  const TrilinosWrappers::types::int_type mode_index =
168  constant_modes_are_global ?
169  TrilinosWrappers::global_index(domain_map, row) :
170  row;
171  distributed_constant_modes[d][row] =
172  constant_modes[d][mode_index];
173  }
174  }
175  (void)expected_mode_size;
176 
177  parameter_list.set("null space: type", "pre-computed");
178  parameter_list.set("null space: dimension",
179  distributed_constant_modes.NumVectors());
180  parameter_list.set("null space: vectors",
181  distributed_constant_modes.Values());
182  }
183  }
184 
185 
186 
187  void
189  Teuchos::ParameterList & parameter_list,
190  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
191  const SparseMatrix & matrix) const
192  {
193  return set_parameters(parameter_list,
194  distributed_constant_modes,
195  matrix.trilinos_matrix());
196  }
197 
198 
199 
200  void
202  Teuchos::ParameterList & parameter_list,
203  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
204  const SparseMatrix & matrix) const
205  {
206  return set_operator_null_space(parameter_list,
207  distributed_constant_modes,
208  matrix.trilinos_matrix());
209  }
210 
211 
212 
214  {
215  preconditioner.reset();
216  trilinos_matrix.reset();
217  }
218 
219 
220 
221  void
223  const AdditionalData &additional_data)
224  {
225  initialize(matrix.trilinos_matrix(), additional_data);
226  }
227 
228 
229 
230  void
231  PreconditionAMG::initialize(const Epetra_RowMatrix &matrix,
232  const AdditionalData & additional_data)
233  {
234  // Build the AMG preconditioner.
235  Teuchos::ParameterList ml_parameters;
236  std::unique_ptr<Epetra_MultiVector> distributed_constant_modes;
237  additional_data.set_parameters(ml_parameters,
238  distributed_constant_modes,
239  matrix);
240 
241  initialize(matrix, ml_parameters);
242 
243  if (additional_data.output_details)
244  {
245  ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
246  dynamic_cast<ML_Epetra::MultiLevelPreconditioner *>(
247  preconditioner.get());
248  Assert(multilevel_operator != nullptr,
249  ExcMessage("Preconditioner setup failed."));
250  multilevel_operator->PrintUnused(0);
251  }
252  }
253 
254 
255 
256  void
258  const Teuchos::ParameterList &ml_parameters)
259  {
260  initialize(matrix.trilinos_matrix(), ml_parameters);
261  }
262 
263 
264 
265  void
266  PreconditionAMG::initialize(const Epetra_RowMatrix & matrix,
267  const Teuchos::ParameterList &ml_parameters)
268  {
269  preconditioner.reset(
270  new ML_Epetra::MultiLevelPreconditioner(matrix, ml_parameters));
271  }
272 
273 
274 
275  template <typename number>
276  void
278  const ::SparseMatrix<number> &deal_ii_sparse_matrix,
279  const AdditionalData & additional_data,
280  const double drop_tolerance,
281  const ::SparsityPattern * use_this_sparsity)
282  {
283  preconditioner.reset();
284  const size_type n_rows = deal_ii_sparse_matrix.m();
285 
286  // Init Epetra Matrix using an equidistributed map; avoid storing the
287  // nonzero elements.
288  IndexSet distributor(n_rows);
289  const unsigned int n_mpi_processes = communicator.NumProc();
290  const unsigned int my_id = communicator.MyPID();
291  distributor.add_range(my_id * n_rows / n_mpi_processes,
292  (my_id + 1) * n_rows / n_mpi_processes);
293 
294  if (trilinos_matrix.get() == nullptr)
295  trilinos_matrix = std::make_shared<SparseMatrix>();
296 
297  trilinos_matrix->reinit(distributor,
298  distributor,
299  deal_ii_sparse_matrix,
300  communicator.Comm(),
301  drop_tolerance,
302  true,
303  use_this_sparsity);
304 
305  initialize(*trilinos_matrix, additional_data);
306  }
307 
308 
309 
310  void
312  {
313  ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
314  dynamic_cast<ML_Epetra::MultiLevelPreconditioner *>(preconditioner.get());
315  multilevel_operator->ReComputePreconditioner();
316  }
317 
318 
319 
320  void
322  {
324  trilinos_matrix.reset();
325  }
326 
327 
328 
331  {
332  unsigned int memory = sizeof(*this);
333 
334  // todo: find a way to read out ML's data
335  // sizes
336  if (trilinos_matrix.get() != nullptr)
337  memory += trilinos_matrix->memory_consumption();
338  return memory;
339  }
340 
341 
342 
343 # ifndef DOXYGEN
344  // explicit instantiations
345  template void
346  PreconditionAMG::initialize(const ::SparseMatrix<double> &,
347  const AdditionalData &,
348  const double,
349  const ::SparsityPattern *);
350  template void
351  PreconditionAMG::initialize(const ::SparseMatrix<float> &,
352  const AdditionalData &,
353  const double,
354  const ::SparsityPattern *);
355 # endif
356 
357 
358 
359 } // namespace TrilinosWrappers
360 
362 
363 #endif // DEAL_II_WITH_TRILINOS
AdditionalData(const bool elliptic=true, const bool higher_order_elements=false, const unsigned int n_cycles=1, const bool w_cyle=false, const double aggregation_threshold=1e-4, const std::vector< std::vector< bool >> &constant_modes=std::vector< std::vector< bool >>(0), const unsigned int smoother_sweeps=2, const unsigned int smoother_overlap=0, const bool output_details=false, const char *smoother_type="Chebyshev", const char *coarse_type="Amesos-KLU")
Contents is actually a matrix.
void set_parameters(Teuchos::ParameterList &parameter_list, std::unique_ptr< Epetra_MultiVector > &distributed_constant_modes, const Epetra_RowMatrix &matrix) const
TrilinosWrappers::types::int_type global_length(const Epetra_MultiVector &vector)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::types::int_type n_global_rows(const Epetra_CrsGraph &graph)
#define Assert(cond, exc)
Definition: exceptions.h:1461
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
const Epetra_CrsMatrix & trilinos_matrix() const
std::shared_ptr< SparseMatrix > trilinos_matrix
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1678
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
void set_operator_null_space(Teuchos::ParameterList &parameter_list, std::unique_ptr< Epetra_MultiVector > &distributed_constant_modes, const Epetra_RowMatrix &matrix) const
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
Teuchos::RCP< Epetra_Operator > preconditioner