Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
trilinos_precondition_ml.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2022 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
35namespace 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
84 if (higher_order_elements)
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,
163 ExcDimensionMismatch(constant_modes[d].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 static_cast<double>(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
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1697
std::shared_ptr< SparseMatrix > trilinos_matrix
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
Teuchos::RCP< Epetra_Operator > preconditioner
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
TrilinosWrappers::types::int_type n_global_rows(const Epetra_CrsGraph &graph)
TrilinosWrappers::types::int_type global_length(const Epetra_MultiVector &vector)
AdditionalData(const bool elliptic=true, const bool higher_order_elements=false, const unsigned int n_cycles=1, const bool w_cycle=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_operator_null_space(Teuchos::ParameterList &parameter_list, std::unique_ptr< Epetra_MultiVector > &distributed_constant_modes, const Epetra_RowMatrix &matrix) const
void set_parameters(Teuchos::ParameterList &parameter_list, std::unique_ptr< Epetra_MultiVector > &distributed_constant_modes, const Epetra_RowMatrix &matrix) const