Reference documentation for deal.II version Git 193422c69f 2020-07-08 17:07:46 +0200
\(\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_muelu.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 # ifdef DEAL_II_TRILINOS_WITH_MUELU
23 
24 # include <MueLu_CreateEpetraPreconditioner.hpp>
25 # include <ml_MultiLevelPreconditioner.h>
26 
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 equidistributed map; avoid storing the
218  // nonzero elements.
219  IndexSet distributor(n_rows);
220  const unsigned int n_mpi_processes = communicator.NumProc();
221  const unsigned int my_id = communicator.MyPID();
222  distributor.add_range(my_id * n_rows / n_mpi_processes,
223  (my_id + 1) * n_rows / n_mpi_processes);
224 
225  if (trilinos_matrix.get() == nullptr)
226  trilinos_matrix = std::make_shared<SparseMatrix>();
227 
228  trilinos_matrix->reinit(distributor,
229  distributor,
230  deal_ii_sparse_matrix,
231  communicator.Comm(),
232  drop_tolerance,
233  true,
234  use_this_sparsity);
235 
236  initialize(*trilinos_matrix, additional_data);
237  }
238 
239 
240 
241  void
243  {
245  trilinos_matrix.reset();
246  }
247 
248 
249 
252  {
253  unsigned int memory = sizeof(*this);
254 
255  // todo: find a way to read out ML's data
256  // sizes
257  if (trilinos_matrix.get() != nullptr)
258  memory += trilinos_matrix->memory_consumption();
259  return memory;
260  }
261 
262 
263 
264 # ifndef DOXYGEN
265  // explicit instantiations
266  template void
267  PreconditionAMGMueLu::initialize(const ::SparseMatrix<double> &,
268  const AdditionalData &,
269  const double,
270  const ::SparsityPattern *);
271  template void
272  PreconditionAMGMueLu::initialize(const ::SparseMatrix<float> &,
273  const AdditionalData &,
274  const double,
275  const ::SparsityPattern *);
276 # endif
277 
278 } // namespace TrilinosWrappers
279 
281 
282 # endif // DEAL_II_TRILINOS_WITH_MUELU
283 #endif // DEAL_II_WITH_TRILINOS
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
Contents is actually a matrix.
TrilinosWrappers::types::int_type global_length(const Epetra_MultiVector &vector)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::types::int_type n_global_rows(const Epetra_CrsGraph &graph)
#define Assert(cond, exc)
Definition: exceptions.h:1403
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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")
const Epetra_CrsMatrix & trilinos_matrix() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::shared_ptr< SparseMatrix > trilinos_matrix
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:1671
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
Teuchos::RCP< Epetra_Operator > preconditioner