Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
trilinos_precondition_ml.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2019 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 
16 #include <deal.II/lac/trilinos_index_access.h>
17 #include <deal.II/lac/trilinos_precondition.h>
18 
19 #ifdef DEAL_II_WITH_TRILINOS
20 
21 # include <deal.II/lac/sparse_matrix.h>
22 # include <deal.II/lac/trilinos_index_access.h>
23 # include <deal.II/lac/trilinos_sparse_matrix.h>
24 # include <deal.II/lac/vector.h>
25 
26 # include <Epetra_MultiVector.h>
27 # include <Ifpack.h>
28 # include <Ifpack_Chebyshev.h>
29 # include <Teuchos_ParameterList.hpp>
30 # include <Teuchos_RCP.hpp>
31 # include <ml_MultiLevelPreconditioner.h>
32 # include <ml_include.h>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 namespace TrilinosWrappers
37 {
38  /* -------------------------- PreconditionAMG -------------------------- */
39 
41  const bool elliptic,
42  const bool higher_order_elements,
43  const unsigned int n_cycles,
44  const bool w_cycle,
45  const double aggregation_threshold,
46  const std::vector<std::vector<bool>> &constant_modes,
47  const unsigned int smoother_sweeps,
48  const unsigned int smoother_overlap,
49  const bool output_details,
50  const char * smoother_type,
51  const char * coarse_type)
52  : elliptic(elliptic)
53  , higher_order_elements(higher_order_elements)
54  , n_cycles(n_cycles)
55  , w_cycle(w_cycle)
56  , aggregation_threshold(aggregation_threshold)
57  , constant_modes(constant_modes)
58  , smoother_sweeps(smoother_sweeps)
59  , smoother_overlap(smoother_overlap)
60  , output_details(output_details)
61  , smoother_type(smoother_type)
62  , coarse_type(coarse_type)
63  {}
64 
65 
66 
67  void
69  Teuchos::ParameterList & parameter_list,
70  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
71  const Epetra_RowMatrix & matrix) const
72  {
73  if (elliptic == true)
74  {
75  ML_Epetra::SetDefaults("SA", parameter_list);
76 
77  // uncoupled mode can give a lot of warnings or even fail when there
78  // are too many entries per row and aggreggation gets complicated, but
79  // MIS does not work if too few elements are located on one
80  // processor. work around these warnings by choosing the different
81  // strategies in different situations: for low order, always use the
82  // standard choice uncoupled. if higher order, right now we also just
83  // use Uncoupled, but we should be aware that maybe MIS might be
84  // needed
85  if (higher_order_elements)
86  parameter_list.set("aggregation: type", "Uncoupled");
87  }
88  else
89  {
90  ML_Epetra::SetDefaults("NSSA", parameter_list);
91  parameter_list.set("aggregation: type", "Uncoupled");
92  parameter_list.set("aggregation: block scaling", true);
93  }
94 
95  parameter_list.set("smoother: type", smoother_type);
96  parameter_list.set("coarse: type", coarse_type);
97 
98  // Force re-initialization of the random seed to make ML deterministic
99  // (only supported in trilinos >12.2):
100 # if DEAL_II_TRILINOS_VERSION_GTE(12, 4, 0)
101  parameter_list.set("initialize random seed", true);
102 # endif
103 
104  parameter_list.set("smoother: sweeps", static_cast<int>(smoother_sweeps));
105  parameter_list.set("cycle applications", static_cast<int>(n_cycles));
106  if (w_cycle == true)
107  parameter_list.set("prec type", "MGW");
108  else
109  parameter_list.set("prec type", "MGV");
110 
111  parameter_list.set("smoother: Chebyshev alpha", 10.);
112  parameter_list.set("smoother: ifpack overlap",
113  static_cast<int>(smoother_overlap));
114  parameter_list.set("aggregation: threshold", aggregation_threshold);
115  parameter_list.set("coarse: max size", 2000);
116 
117  if (output_details)
118  parameter_list.set("ML output", 10);
119  else
120  parameter_list.set("ML output", 0);
121 
122  set_operator_null_space(parameter_list, distributed_constant_modes, matrix);
123  }
124 
125 
126 
127  void
129  Teuchos::ParameterList & parameter_list,
130  std::unique_ptr<Epetra_MultiVector> &ptr_distributed_constant_modes,
131  const Epetra_RowMatrix & matrix) const
132  {
133  const Epetra_Map &domain_map = matrix.OperatorDomainMap();
134 
135  const size_type constant_modes_dimension = constant_modes.size();
136  ptr_distributed_constant_modes.reset(new Epetra_MultiVector(
137  domain_map, constant_modes_dimension > 0 ? constant_modes_dimension : 1));
138  Assert(ptr_distributed_constant_modes, ExcNotInitialized());
139  Epetra_MultiVector &distributed_constant_modes =
140  *ptr_distributed_constant_modes;
141 
142  if (constant_modes_dimension > 0)
143  {
144  const size_type global_size = TrilinosWrappers::n_global_rows(matrix);
145  (void)global_length; // work around compiler warning about unused
146  // function in release mode
147  Assert(global_size ==
148  static_cast<size_type>(
149  TrilinosWrappers::global_length(distributed_constant_modes)),
150  ExcDimensionMismatch(global_size,
152  distributed_constant_modes)));
153  const bool constant_modes_are_global =
154  constant_modes[0].size() == global_size;
155  const size_type my_size = domain_map.NumMyElements();
156 
157  // Reshape null space as a contiguous vector of doubles so that
158  // Trilinos can read from it.
159  const size_type expected_mode_size =
160  constant_modes_are_global ? global_size : my_size;
161  for (size_type d = 0; d < constant_modes_dimension; ++d)
162  {
163  Assert(constant_modes[d].size() == expected_mode_size,
164  ExcDimensionMismatch(constant_modes[d].size(),
165  expected_mode_size));
166  for (size_type row = 0; row < my_size; ++row)
167  {
168  const TrilinosWrappers::types::int_type mode_index =
169  constant_modes_are_global ?
170  TrilinosWrappers::global_index(domain_map, row) :
171  row;
172  distributed_constant_modes[d][row] =
173  constant_modes[d][mode_index];
174  }
175  }
176  (void)expected_mode_size;
177 
178  parameter_list.set("null space: type", "pre-computed");
179  parameter_list.set("null space: dimension",
180  distributed_constant_modes.NumVectors());
181  parameter_list.set("null space: vectors",
182  distributed_constant_modes.Values());
183  }
184  }
185 
186 
187 
188  void
190  Teuchos::ParameterList & parameter_list,
191  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
192  const SparseMatrix & matrix) const
193  {
194  return set_parameters(parameter_list,
195  distributed_constant_modes,
196  matrix.trilinos_matrix());
197  }
198 
199 
200 
201  void
203  Teuchos::ParameterList & parameter_list,
204  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
205  const SparseMatrix & matrix) const
206  {
207  return set_operator_null_space(parameter_list,
208  distributed_constant_modes,
209  matrix.trilinos_matrix());
210  }
211 
212 
213 
215  {
216  preconditioner.reset();
217  trilinos_matrix.reset();
218  }
219 
220 
221 
222  void
224  const AdditionalData &additional_data)
225  {
226  initialize(matrix.trilinos_matrix(), additional_data);
227  }
228 
229 
230 
231  void
232  PreconditionAMG::initialize(const Epetra_RowMatrix &matrix,
233  const AdditionalData & additional_data)
234  {
235  // Build the AMG preconditioner.
236  Teuchos::ParameterList ml_parameters;
237  std::unique_ptr<Epetra_MultiVector> distributed_constant_modes;
238  additional_data.set_parameters(ml_parameters,
239  distributed_constant_modes,
240  matrix);
241 
242  initialize(matrix, ml_parameters);
243 
244  if (additional_data.output_details)
245  {
246  ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
247  dynamic_cast<ML_Epetra::MultiLevelPreconditioner *>(
248  preconditioner.get());
249  Assert(multilevel_operator != nullptr,
250  ExcMessage("Preconditioner setup failed."));
251  multilevel_operator->PrintUnused(0);
252  }
253  }
254 
255 
256 
257  void
259  const Teuchos::ParameterList &ml_parameters)
260  {
261  initialize(matrix.trilinos_matrix(), ml_parameters);
262  }
263 
264 
265 
266  void
267  PreconditionAMG::initialize(const Epetra_RowMatrix & matrix,
268  const Teuchos::ParameterList &ml_parameters)
269  {
270  preconditioner.reset(
271  new ML_Epetra::MultiLevelPreconditioner(matrix, ml_parameters));
272  }
273 
274 
275 
276  template <typename number>
277  void
279  const ::SparseMatrix<number> &deal_ii_sparse_matrix,
280  const AdditionalData & additional_data,
281  const double drop_tolerance,
282  const ::SparsityPattern * use_this_sparsity)
283  {
284  preconditioner.reset();
285  const size_type n_rows = deal_ii_sparse_matrix.m();
286 
287  // Init Epetra Matrix using an
288  // equidistributed map; avoid
289  // storing the nonzero
290  // elements.
291  vector_distributor = std::make_shared<Epetra_Map>(
292  static_cast<TrilinosWrappers::types::int_type>(n_rows), 0, communicator);
293 
294  if (trilinos_matrix.get() == nullptr)
295  trilinos_matrix = std::make_shared<SparseMatrix>();
296 
299  deal_ii_sparse_matrix,
300  drop_tolerance,
301  true,
302  use_this_sparsity);
303 
304  initialize(*trilinos_matrix, additional_data);
305  }
306 
307 
308 
309  void
311  {
312  ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
313  dynamic_cast<ML_Epetra::MultiLevelPreconditioner *>(preconditioner.get());
314  multilevel_operator->ReComputePreconditioner();
315  }
316 
317 
318 
319  void
321  {
323  trilinos_matrix.reset();
324  }
325 
326 
327 
330  {
331  unsigned int memory = sizeof(*this);
332 
333  // todo: find a way to read out ML's data
334  // sizes
335  if (trilinos_matrix.get() != nullptr)
336  memory += trilinos_matrix->memory_consumption();
337  return memory;
338  }
339 
340 
341 
342  // explicit instantiations
343  template void
344  PreconditionAMG::initialize(const ::SparseMatrix<double> &,
345  const AdditionalData &,
346  const double,
347  const ::SparsityPattern *);
348  template void
349  PreconditionAMG::initialize(const ::SparseMatrix<float> &,
350  const AdditionalData &,
351  const double,
352  const ::SparsityPattern *);
353 
354 
355 
356 } // namespace TrilinosWrappers
357 
358 DEAL_II_NAMESPACE_CLOSE
359 
360 #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")
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)
std::shared_ptr< Epetra_Map > vector_distributor
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:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::shared_ptr< SparseMatrix > trilinos_matrix
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
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