deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20:00+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17
18#ifdef DEAL_II_WITH_TRILINOS
19
22# include <deal.II/lac/vector.h>
23
24# include <Epetra_MultiVector.h>
25# include <Ifpack.h>
26# include <Ifpack_Chebyshev.h>
27# include <Teuchos_ParameterList.hpp>
28# include <Teuchos_RCP.hpp>
29# include <ml_MultiLevelPreconditioner.h>
30# include <ml_include.h>
31
33
34namespace TrilinosWrappers
35{
36 /* -------------------------- PreconditionAMG -------------------------- */
37
39 const bool elliptic,
40 const bool higher_order_elements,
41 const unsigned int n_cycles,
42 const bool w_cycle,
43 const double aggregation_threshold,
44 const std::vector<std::vector<bool>> &constant_modes,
45 const unsigned int smoother_sweeps,
46 const unsigned int smoother_overlap,
47 const bool output_details,
48 const char *smoother_type,
49 const char *coarse_type)
50 : elliptic(elliptic)
51 , higher_order_elements(higher_order_elements)
52 , n_cycles(n_cycles)
53 , w_cycle(w_cycle)
54 , aggregation_threshold(aggregation_threshold)
55 , constant_modes(constant_modes)
56 , smoother_sweeps(smoother_sweeps)
57 , smoother_overlap(smoother_overlap)
58 , output_details(output_details)
59 , smoother_type(smoother_type)
60 , coarse_type(coarse_type)
61 {}
62
63
64
65 void
67 Teuchos::ParameterList &parameter_list,
68 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
69 const Epetra_RowMatrix &matrix) const
70 {
71 if (elliptic == true)
72 {
73 ML_Epetra::SetDefaults("SA", parameter_list);
74
75 // uncoupled mode can give a lot of warnings or even fail when there
76 // are too many entries per row and aggregation gets complicated, but
77 // MIS does not work if too few elements are located on one
78 // processor. work around these warnings by choosing the different
79 // strategies in different situations: for low order, always use the
80 // standard choice uncoupled. if higher order, right now we also just
81 // use Uncoupled, but we should be aware that maybe MIS might be
82 // needed
83 if (higher_order_elements)
84 parameter_list.set("aggregation: type", "Uncoupled");
85 }
86 else
87 {
88 ML_Epetra::SetDefaults("NSSA", parameter_list);
89 parameter_list.set("aggregation: type", "Uncoupled");
90 parameter_list.set("aggregation: block scaling", true);
91 }
92
93 parameter_list.set("smoother: type", smoother_type);
94 parameter_list.set("coarse: type", coarse_type);
95
96 // Force re-initialization of the random seed to make ML deterministic
97 // (only supported in trilinos >12.2):
98# if DEAL_II_TRILINOS_VERSION_GTE(12, 4, 0)
99 parameter_list.set("initialize random seed", true);
100# endif
101
102 parameter_list.set("smoother: sweeps", static_cast<int>(smoother_sweeps));
103 parameter_list.set("cycle applications", static_cast<int>(n_cycles));
104 if (w_cycle == true)
105 parameter_list.set("prec type", "MGW");
106 else
107 parameter_list.set("prec type", "MGV");
108
109 parameter_list.set("smoother: Chebyshev alpha", 10.);
110 parameter_list.set("smoother: ifpack overlap",
111 static_cast<int>(smoother_overlap));
112 parameter_list.set("aggregation: threshold", aggregation_threshold);
113 parameter_list.set("coarse: max size", 2000);
114
115 if (output_details)
116 parameter_list.set("ML output", 10);
117 else
118 parameter_list.set("ML output", 0);
119
120 set_operator_null_space(parameter_list, distributed_constant_modes, matrix);
121 }
122
123
124
125 void
127 Teuchos::ParameterList &parameter_list,
128 std::unique_ptr<Epetra_MultiVector> &ptr_distributed_constant_modes,
129 const Epetra_RowMatrix &matrix) const
130 {
131 const auto run = [&](const auto &constant_modes) {
132 const Epetra_Map &domain_map = matrix.OperatorDomainMap();
133
134 const size_type constant_modes_dimension = constant_modes.size();
135 ptr_distributed_constant_modes =
136 std::make_unique<Epetra_MultiVector>(domain_map,
137 constant_modes_dimension > 0 ?
138 constant_modes_dimension :
139 1);
140 Assert(ptr_distributed_constant_modes, ExcNotInitialized());
141 Epetra_MultiVector &distributed_constant_modes =
142 *ptr_distributed_constant_modes;
143
144 if (constant_modes_dimension > 0)
145 {
146 const size_type global_size = TrilinosWrappers::n_global_rows(matrix);
147 (void)global_length; // work around compiler warning about unused
148 // function in release mode
149 Assert(global_size ==
151 distributed_constant_modes)),
152 ExcDimensionMismatch(global_size,
154 distributed_constant_modes)));
155 const bool constant_modes_are_global =
156 constant_modes[0].size() == global_size;
157 const size_type my_size = domain_map.NumMyElements();
158
159 // Reshape null space as a contiguous vector of doubles so that
160 // Trilinos can read from it.
161 const size_type expected_mode_size =
162 constant_modes_are_global ? global_size : my_size;
163 for (size_type d = 0; d < constant_modes_dimension; ++d)
164 {
165 Assert(constant_modes[d].size() == expected_mode_size,
166 ExcDimensionMismatch(constant_modes[d].size(),
167 expected_mode_size));
168 for (size_type row = 0; row < my_size; ++row)
169 {
170 const TrilinosWrappers::types::int_type mode_index =
171 constant_modes_are_global ?
172 TrilinosWrappers::global_index(domain_map, row) :
173 row;
174 distributed_constant_modes[d][row] =
175 static_cast<double>(constant_modes[d][mode_index]);
176 }
177 }
178 (void)expected_mode_size;
179
180 parameter_list.set("null space: type", "pre-computed");
181 parameter_list.set("null space: dimension",
182 distributed_constant_modes.NumVectors());
183 parameter_list.set("null space: vectors",
184 distributed_constant_modes.Values());
185 }
186 };
187
188 if (!constant_modes_values.empty())
189 {
190 AssertDimension(constant_modes.size(), 0);
191 run(constant_modes_values);
192 }
193 else
194 run(constant_modes);
195 }
196
197
198
199 void
201 Teuchos::ParameterList &parameter_list,
202 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
203 const SparseMatrix &matrix) const
204 {
205 return set_parameters(parameter_list,
206 distributed_constant_modes,
207 matrix.trilinos_matrix());
208 }
209
210
211
212 void
214 Teuchos::ParameterList &parameter_list,
215 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
216 const SparseMatrix &matrix) const
217 {
218 return set_operator_null_space(parameter_list,
219 distributed_constant_modes,
220 matrix.trilinos_matrix());
221 }
222
223
224
226 {
227 preconditioner.reset();
228 trilinos_matrix.reset();
229 }
230
231
232
233 void
235 const AdditionalData &additional_data)
236 {
237 initialize(matrix.trilinos_matrix(), additional_data);
238 }
239
240
241
242 void
243 PreconditionAMG::initialize(const Epetra_RowMatrix &matrix,
244 const AdditionalData &additional_data)
245 {
246 // Build the AMG preconditioner.
247 Teuchos::ParameterList ml_parameters;
248 std::unique_ptr<Epetra_MultiVector> distributed_constant_modes;
249 additional_data.set_parameters(ml_parameters,
250 distributed_constant_modes,
251 matrix);
252
253 initialize(matrix, ml_parameters);
254
255 if (additional_data.output_details)
256 {
257 ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
258 dynamic_cast<ML_Epetra::MultiLevelPreconditioner *>(
259 preconditioner.get());
260 Assert(multilevel_operator != nullptr,
261 ExcMessage("Preconditioner setup failed."));
262 multilevel_operator->PrintUnused(0);
263 }
264 }
265
266
267
268 void
270 const Teuchos::ParameterList &ml_parameters)
271 {
272 initialize(matrix.trilinos_matrix(), ml_parameters);
273 }
274
275
276
277 void
278 PreconditionAMG::initialize(const Epetra_RowMatrix &matrix,
279 const Teuchos::ParameterList &ml_parameters)
280 {
281 preconditioner.reset(
282 new ML_Epetra::MultiLevelPreconditioner(matrix, ml_parameters));
283 }
284
285
286
287 template <typename number>
288 void
290 const ::SparseMatrix<number> &deal_ii_sparse_matrix,
291 const AdditionalData &additional_data,
292 const double drop_tolerance,
293 const ::SparsityPattern *use_this_sparsity)
294 {
295 preconditioner.reset();
296 const size_type n_rows = deal_ii_sparse_matrix.m();
297
298 // Init Epetra Matrix using an equidistributed map; avoid storing the
299 // nonzero elements.
300 IndexSet distributor(n_rows);
301 const unsigned int n_mpi_processes = communicator.NumProc();
302 const unsigned int my_id = communicator.MyPID();
303 distributor.add_range(my_id * n_rows / n_mpi_processes,
304 (my_id + 1) * n_rows / n_mpi_processes);
305
306 if (trilinos_matrix.get() == nullptr)
307 trilinos_matrix = std::make_shared<SparseMatrix>();
308
309 trilinos_matrix->reinit(distributor,
310 distributor,
311 deal_ii_sparse_matrix,
312 communicator.Comm(),
313 drop_tolerance,
314 true,
315 use_this_sparsity);
316
317 initialize(*trilinos_matrix, additional_data);
318 }
319
320
321
322 void
324 {
325 ML_Epetra::MultiLevelPreconditioner *multilevel_operator =
326 dynamic_cast<ML_Epetra::MultiLevelPreconditioner *>(preconditioner.get());
327 multilevel_operator->ReComputePreconditioner();
328 }
329
330
331
332 void
338
339
340
343 {
344 unsigned int memory = sizeof(*this);
345
346 // todo: find a way to read out ML's data
347 // sizes
348 if (trilinos_matrix.get() != nullptr)
349 memory += trilinos_matrix->memory_consumption();
350 return memory;
351 }
352
353
354
355# ifndef DOXYGEN
356 // explicit instantiations
357 template void
358 PreconditionAMG::initialize(const ::SparseMatrix<double> &,
359 const AdditionalData &,
360 const double,
361 const ::SparsityPattern *);
362 template void
363 PreconditionAMG::initialize(const ::SparseMatrix<float> &,
364 const AdditionalData &,
365 const double,
366 const ::SparsityPattern *);
367# endif
368
369
370
371} // namespace TrilinosWrappers
372
374
375#endif // DEAL_II_WITH_TRILINOS
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1803
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
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