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_muelu.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# ifdef DEAL_II_TRILINOS_WITH_MUELU
23
24# include <MueLu_CreateEpetraPreconditioner.hpp>
25# include <ml_MultiLevelPreconditioner.h>
26
28
29namespace 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,
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)),
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 {
161 constant_modes_are_global ?
162 TrilinosWrappers::global_index(domain_map, row) :
163 row;
164 distributed_constant_modes[d][row] = static_cast<double>(
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 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 & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
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 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")