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