Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
trilinos_precondition_muelu.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 # if DEAL_II_TRILINOS_VERSION_GTE(11, 14, 0)
21 # include <deal.II/lac/sparse_matrix.h>
22 # include <deal.II/lac/trilinos_sparse_matrix.h>
23 
24 # include <MueLu_CreateEpetraPreconditioner.hpp>
25 # include <ml_MultiLevelPreconditioner.h>
26 
27 DEAL_II_NAMESPACE_OPEN
28 
29 namespace TrilinosWrappers
30 {
32  const bool elliptic,
33  const unsigned int n_cycles,
34  const bool w_cycle,
35  const double aggregation_threshold,
36  const std::vector<std::vector<bool>> &constant_modes,
37  const unsigned int smoother_sweeps,
38  const unsigned int smoother_overlap,
39  const bool output_details,
40  const char * smoother_type,
41  const char * coarse_type)
42  : elliptic(elliptic)
43  , n_cycles(n_cycles)
44  , w_cycle(w_cycle)
45  , aggregation_threshold(aggregation_threshold)
46  , constant_modes(constant_modes)
47  , smoother_sweeps(smoother_sweeps)
48  , smoother_overlap(smoother_overlap)
49  , output_details(output_details)
50  , smoother_type(smoother_type)
51  , coarse_type(coarse_type)
52  {}
53 
54 
55 
57  {
58  // clang-tidy wants to default the constructor if we disable the check
59  // in case we compile without 64-bit indices
60 # ifdef DEAL_II_WITH_64BIT_INDICES
61  constexpr bool enabled = false;
62 # else
63  constexpr bool enabled = true;
64 # endif
65  AssertThrow(enabled,
66  ExcMessage(
67  "PreconditionAMGMueLu does not support 64bit-indices!"));
68  }
69 
70 
71 
72  void
74  const AdditionalData &additional_data)
75  {
76  initialize(matrix.trilinos_matrix(), additional_data);
77  }
78 
79 
80 
81  void
82  PreconditionAMGMueLu::initialize(const Epetra_CrsMatrix &matrix,
83  const AdditionalData & additional_data)
84  {
85  // Build the AMG preconditioner.
86  Teuchos::ParameterList parameter_list;
87 
88  parameter_list.set("parameterlist: syntax", "ml");
89 
90  if (additional_data.elliptic == true)
91  ML_Epetra::SetDefaults("SA", parameter_list);
92  else
93  {
94  ML_Epetra::SetDefaults("NSSA", parameter_list);
95  parameter_list.set("aggregation: block scaling", true);
96  }
97  // MIS does not exist anymore, only choice are uncoupled and coupled. When
98  // using uncoupled, aggregates cannot span multiple processes. When using
99  // coupled aggregates can span multiple processes.
100  parameter_list.set("aggregation: type", "Uncoupled");
101 
102  parameter_list.set("smoother: type", additional_data.smoother_type);
103  parameter_list.set("coarse: type", additional_data.coarse_type);
104 
105  parameter_list.set("smoother: sweeps",
106  static_cast<int>(additional_data.smoother_sweeps));
107  parameter_list.set("cycle applications",
108  static_cast<int>(additional_data.n_cycles));
109  if (additional_data.w_cycle == true)
110  parameter_list.set("prec type", "MGW");
111  else
112  parameter_list.set("prec type", "MGV");
113 
114  parameter_list.set("smoother: Chebyshev alpha", 10.);
115  parameter_list.set("smoother: ifpack overlap",
116  static_cast<int>(additional_data.smoother_overlap));
117  parameter_list.set("aggregation: threshold",
118  additional_data.aggregation_threshold);
119  parameter_list.set("coarse: max size", 2000);
120 
121  if (additional_data.output_details)
122  parameter_list.set("ML output", 10);
123  else
124  parameter_list.set("ML output", 0);
125 
126  const Epetra_Map &domain_map = matrix.OperatorDomainMap();
127 
128  const size_type constant_modes_dimension =
129  additional_data.constant_modes.size();
130  Epetra_MultiVector distributed_constant_modes(
131  domain_map, constant_modes_dimension > 0 ? constant_modes_dimension : 1);
132  std::vector<double> dummy(constant_modes_dimension);
133 
134  if (constant_modes_dimension > 0)
135  {
136  const size_type n_rows = TrilinosWrappers::n_global_rows(matrix);
137  const bool constant_modes_are_global =
138  additional_data.constant_modes[0].size() == n_rows;
139  const size_type n_relevant_rows =
140  constant_modes_are_global ? n_rows :
141  additional_data.constant_modes[0].size();
142  const size_type my_size = domain_map.NumMyElements();
143  if (constant_modes_are_global == false)
144  Assert(n_relevant_rows == my_size,
145  ExcDimensionMismatch(n_relevant_rows, my_size));
146  Assert(n_rows == static_cast<size_type>(TrilinosWrappers::global_length(
147  distributed_constant_modes)),
148  ExcDimensionMismatch(n_rows,
150  distributed_constant_modes)));
151 
152  (void)n_relevant_rows;
153  (void)global_length;
154 
155  // Reshape null space as a contiguous vector of doubles so that
156  // Trilinos can read from it.
157  for (size_type d = 0; d < constant_modes_dimension; ++d)
158  for (size_type row = 0; row < my_size; ++row)
159  {
160  TrilinosWrappers::types::int_type global_row_id =
161  constant_modes_are_global ?
162  TrilinosWrappers::global_index(domain_map, row) :
163  row;
164  distributed_constant_modes[d][row] =
165  additional_data.constant_modes[d][global_row_id];
166  }
167 
168  parameter_list.set("null space: type", "pre-computed");
169  parameter_list.set("null space: dimension",
170  distributed_constant_modes.NumVectors());
171  if (my_size > 0)
172  parameter_list.set("null space: vectors",
173  distributed_constant_modes.Values());
174  // We need to set a valid pointer to data even if there is no data on
175  // the current processor. Therefore, pass a dummy in that case
176  else
177  parameter_list.set("null space: vectors", dummy.data());
178  }
179 
180  initialize(matrix, parameter_list);
181  }
182 
183 
184 
185  void
187  Teuchos::ParameterList &muelu_parameters)
188  {
189  initialize(matrix.trilinos_matrix(), muelu_parameters);
190  }
191 
192 
193 
194  void
195  PreconditionAMGMueLu::initialize(const Epetra_CrsMatrix &matrix,
196  Teuchos::ParameterList &muelu_parameters)
197  {
198  const auto teuchos_wrapped_matrix =
199  Teuchos::rcp(const_cast<Epetra_CrsMatrix *>(&matrix), false);
200  preconditioner = MueLu::CreateEpetraPreconditioner(teuchos_wrapped_matrix,
201  muelu_parameters);
202  }
203 
204 
205 
206  template <typename number>
207  void
209  const ::SparseMatrix<number> &deal_ii_sparse_matrix,
210  const AdditionalData & additional_data,
211  const double drop_tolerance,
212  const ::SparsityPattern * use_this_sparsity)
213  {
214  preconditioner.reset();
215  const size_type n_rows = deal_ii_sparse_matrix.m();
216 
217  // Init Epetra Matrix using an
218  // equidistributed map; avoid
219  // storing the nonzero
220  // elements.
221  vector_distributor = std::make_shared<Epetra_Map>(
222  static_cast<TrilinosWrappers::types::int_type>(n_rows), 0, communicator);
223 
224  if (trilinos_matrix.get() == nullptr)
225  trilinos_matrix = std::make_shared<SparseMatrix>();
226 
229  deal_ii_sparse_matrix,
230  drop_tolerance,
231  true,
232  use_this_sparsity);
233 
234  initialize(*trilinos_matrix, additional_data);
235  }
236 
237 
238 
239  void
241  {
243  trilinos_matrix.reset();
244  }
245 
246 
247 
250  {
251  unsigned int memory = sizeof(*this);
252 
253  // todo: find a way to read out ML's data
254  // sizes
255  if (trilinos_matrix.get() != nullptr)
256  memory += trilinos_matrix->memory_consumption();
257  return memory;
258  }
259 
260 
261 
262  // explicit instantiations
263  template void
264  PreconditionAMGMueLu::initialize(const ::SparseMatrix<double> &,
265  const AdditionalData &,
266  const double,
267  const ::SparsityPattern *);
268  template void
269  PreconditionAMGMueLu::initialize(const ::SparseMatrix<float> &,
270  const AdditionalData &,
271  const double,
272  const ::SparsityPattern *);
273 
274 } // namespace TrilinosWrappers
275 
276 DEAL_II_NAMESPACE_CLOSE
277 
278 # endif // DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
279 #endif // DEAL_II_WITH_TRILINOS
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
TrilinosWrappers::types::int_type global_length(const Epetra_MultiVector &vector)
std::shared_ptr< Epetra_Map > vector_distributor
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
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)
AdditionalData(const bool elliptic=true, 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")
std::shared_ptr< SparseMatrix > trilinos_matrix
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
Teuchos::RCP< Epetra_Operator > preconditioner