ExaDG
Loading...
Searching...
No Matches
inverse_mass_operator.h
1/* ______________________________________________________________________
2 *
3 * ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
4 *
5 * Copyright (C) 2021 by the ExaDG authors
6 *
7 * This program is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <https://www.gnu.org/licenses/>.
19 * ______________________________________________________________________
20 */
21
22#ifndef EXADG_OPERATORS_INVERSE_MASS_OPERATOR_H_
23#define EXADG_OPERATORS_INVERSE_MASS_OPERATOR_H_
24
25// deal.II
26#include <deal.II/fe/fe_data.h>
27#include <deal.II/lac/la_parallel_vector.h>
28#include <deal.II/matrix_free/operators.h>
29
30// ExaDG
31#include <exadg/grid/grid_data.h>
32#include <exadg/matrix_free/integrators.h>
33#include <exadg/operators/inverse_mass_parameters.h>
34#include <exadg/operators/mass_operator.h>
35#include <exadg/operators/variable_coefficients.h>
36#include <exadg/solvers_and_preconditioners/preconditioners/block_jacobi_preconditioner.h>
37#include <exadg/solvers_and_preconditioners/preconditioners/jacobi_preconditioner.h>
38
39namespace ExaDG
40{
41template<typename Number>
42struct InverseMassOperatorData
43{
44 InverseMassOperatorData()
45 : dof_index(0),
46 quad_index(0),
47 coefficient_is_variable(false),
48 consider_inverse_coefficient(false),
49 variable_coefficients(nullptr)
50 {
51 }
52
53 // Get optimal in the sense of (most likely) fastest implementation type of the inverse mass
54 // operator depending on the approximation space.
55 template<int dim>
56 static InverseMassType
57 get_optimal_inverse_mass_type(dealii::FiniteElement<dim> const & fe,
58 ElementType const element_type)
59 {
60 if(fe.conforms(dealii::FiniteElementData<dim>::L2))
61 {
62 if(element_type == ElementType::Hypercube)
63 {
64 return InverseMassType::MatrixfreeOperator;
65 }
66 else
67 {
68 return InverseMassType::ElementwiseKrylovSolver;
69 }
70 }
71 else
72 {
73 return InverseMassType::GlobalKrylovSolver;
74 }
75 }
76
77 // Parameters referring to dealii::MatrixFree
78 unsigned int dof_index;
79 unsigned int quad_index;
80
81 InverseMassParameters parameters;
82
83 // Enable variable coefficients.
84 bool coefficient_is_variable;
85
86 // Consider the regular form of the coefficient (1) or its inverse (2):
87 // (1) : (u_h , v_h * c)_Omega
88 // (2) : (u_h , v_h / c)_Omega
89 bool consider_inverse_coefficient;
90
91 VariableCoefficients<dealii::VectorizedArray<Number>> const * variable_coefficients;
92};
93
94template<int dim, int n_components, typename Number>
95class InverseMassOperator
96{
97private:
98 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
99
100 typedef InverseMassOperator<dim, n_components, Number> This;
101
102 typedef CellIntegrator<dim, n_components, Number> Integrator;
103
104 // use a template parameter of -1 to select the precompiled version of this operator
105 typedef dealii::MatrixFreeOperators::CellwiseInverseMassMatrix<dim, -1, n_components, Number>
106 InverseMassAsMatrixFreeOperator;
107
108 typedef std::pair<unsigned int, unsigned int> Range;
109
110public:
111 InverseMassOperator() : matrix_free(nullptr), dof_index(0), quad_index(0)
112 {
113 }
114
115 void
116 initialize(dealii::MatrixFree<dim, Number> const & matrix_free_in,
117 InverseMassOperatorData<Number> const inverse_mass_operator_data,
118 dealii::AffineConstraints<Number> const * constraints = nullptr)
119 {
120 this->matrix_free = &matrix_free_in;
121 dof_index = inverse_mass_operator_data.dof_index;
122 quad_index = inverse_mass_operator_data.quad_index;
123
124 data = inverse_mass_operator_data.parameters;
125
126 coefficient_is_variable = inverse_mass_operator_data.coefficient_is_variable;
127 consider_inverse_coefficient = inverse_mass_operator_data.consider_inverse_coefficient;
128 variable_coefficients = inverse_mass_operator_data.variable_coefficients;
129
130 // Variable coefficients only implemented for the matrix-free operator.
131 AssertThrow(not coefficient_is_variable or variable_coefficients != nullptr,
132 dealii::ExcMessage("Pointer to variable coefficients not set properly."));
133
134 dealii::FiniteElement<dim> const & fe = matrix_free->get_dof_handler(dof_index).get_fe();
135
136 // Some implementation variants of the inverse mass operator are based on assumptions on the
137 // discretization and more efficient choices are available for tensor-product elements or
138 // L2-conforming spaces.
139 if(data.implementation_type == InverseMassType::MatrixfreeOperator)
140 {
141 AssertThrow(
142 fe.base_element(0).dofs_per_cell == dealii::Utilities::pow(fe.degree + 1, dim),
143 dealii::ExcMessage(
144 "The matrix-free cell-wise inverse mass operator is currently only available for isotropic tensor-product elements."));
145
146 AssertThrow(
147 this->matrix_free->get_shape_info(0, quad_index).data[0].n_q_points_1d == fe.degree + 1,
148 dealii::ExcMessage(
149 "The matrix-free cell-wise inverse mass operator is currently only available if n_q_points_1d = n_nodes_1d."));
150
151 AssertThrow(
152 fe.conforms(dealii::FiniteElementData<dim>::L2),
153 dealii::ExcMessage(
154 "The matrix-free cell-wise inverse mass operator is only available for L2-conforming elements."));
155 }
156 else
157 {
158 // Setup MassOperator as underlying operator for cell-wise direct/iterative inverse or global
159 // solve.
160 MassOperatorData<dim, Number> mass_operator_data;
161 mass_operator_data.dof_index = dof_index;
162 mass_operator_data.quad_index = quad_index;
163 mass_operator_data.coefficient_is_variable = coefficient_is_variable;
164 mass_operator_data.variable_coefficients = variable_coefficients;
165 mass_operator_data.consider_inverse_coefficient = consider_inverse_coefficient;
166
167 if(data.implementation_type == InverseMassType::ElementwiseKrylovSolver)
168 {
169 mass_operator_data.implement_block_diagonal_preconditioner_matrix_free = true;
170 mass_operator_data.solver_block_diagonal = Elementwise::Solver::CG;
171 if(inverse_mass_operator_data.parameters.preconditioner == PreconditionerMass::None)
172 {
173 mass_operator_data.preconditioner_block_diagonal = Elementwise::Preconditioner::None;
174 }
175 else if(inverse_mass_operator_data.parameters.preconditioner ==
176 PreconditionerMass::PointJacobi)
177 {
178 mass_operator_data.preconditioner_block_diagonal =
179 Elementwise::Preconditioner::PointJacobi;
180 }
181 mass_operator_data.solver_data_block_diagonal =
182 inverse_mass_operator_data.parameters.solver_data;
183 }
184
185 // Use constraints if provided.
186 if(constraints == nullptr)
187 {
188 dealii::AffineConstraints<Number> dummy_constraints;
189 dummy_constraints.close();
190 mass_operator.initialize(*matrix_free, dummy_constraints, mass_operator_data);
191 }
192 else
193 {
194 mass_operator.initialize(*matrix_free, *constraints, mass_operator_data);
195 }
196
197 if(data.implementation_type == InverseMassType::ElementwiseKrylovSolver or
198 data.implementation_type == InverseMassType::BlockMatrices)
199 {
200 // Non-L2-conforming elements are asserted here because the cell-wise inverse neglecting
201 // cell-coupling terms is only an approximation of the inverse mass matrix.
202 AssertThrow(
203 fe.conforms(dealii::FiniteElementData<dim>::L2),
204 dealii::ExcMessage(
205 "The cell-wise inverse mass operator is only available for L2-conforming elements."));
206
207 block_jacobi_preconditioner =
208 std::make_shared<BlockJacobiPreconditioner<MassOperator<dim, n_components, Number>>>(
209 mass_operator, true /* initialize_preconditioner */);
210 }
211 else if(data.implementation_type == InverseMassType::GlobalKrylovSolver)
212 {
213 Krylov::SolverDataCG solver_data;
214 solver_data.max_iter = inverse_mass_operator_data.parameters.solver_data.max_iter;
215 solver_data.solver_tolerance_abs =
216 inverse_mass_operator_data.parameters.solver_data.abs_tol;
217 solver_data.solver_tolerance_rel =
218 inverse_mass_operator_data.parameters.solver_data.rel_tol;
219
220 solver_data.use_preconditioner =
221 inverse_mass_operator_data.parameters.preconditioner != PreconditionerMass::None;
222 if(inverse_mass_operator_data.parameters.preconditioner == PreconditionerMass::None)
223 {
224 // no setup required.
225 }
226 else if(inverse_mass_operator_data.parameters.preconditioner ==
227 PreconditionerMass::PointJacobi)
228 {
229 global_preconditioner =
230 std::make_shared<JacobiPreconditioner<MassOperator<dim, n_components, Number>>>(
231 mass_operator, true /* initialize_preconditioner */);
232 }
233 else if(inverse_mass_operator_data.parameters.preconditioner ==
234 PreconditionerMass::BlockJacobi)
235 {
236 global_preconditioner =
237 std::make_shared<BlockJacobiPreconditioner<MassOperator<dim, n_components, Number>>>(
238 mass_operator, true /* initialize_preconditioner */);
239 }
240 else
241 {
242 AssertThrow(false, dealii::ExcMessage("This `PreconditionerMass` is not implemented."));
243 }
244
245 global_solver = std::make_shared<Krylov::SolverCG<MassOperator<dim, n_components, Number>,
247 VectorType>>(mass_operator,
248 *global_preconditioner,
249 solver_data);
250 }
251 else
252 {
253 AssertThrow(false,
254 dealii::ExcMessage("The specified `InverseMassType` is not implemented."));
255 }
256 }
257 }
258
263 void
265 {
266 if(data.implementation_type == InverseMassType::MatrixfreeOperator)
267 {
268 // no updates needed as long as the MatrixFree object is up-to-date (which is not the
269 // responsibility of the present class).
270 }
271 else if(data.implementation_type == InverseMassType::ElementwiseKrylovSolver or
272 data.implementation_type == InverseMassType::BlockMatrices)
273 {
274 // the mass operator does not need to be updated as long as the MatrixFree object is
275 // up-to-date (which is not the responsibility of the present class).
276
277 // update the matrix-based components of the block-Jacobi preconditioner
278 block_jacobi_preconditioner->update();
279 }
280 else if(data.implementation_type == InverseMassType::GlobalKrylovSolver)
281 {
282 global_preconditioner->update();
283 }
284 else
285 {
286 AssertThrow(false, dealii::ExcMessage("The specified InverseMassType is not implemented."));
287 }
288 }
289
290 // dst = M^-1 * src
291 void
292 apply(VectorType & dst, VectorType const & src) const
293 {
294 if(data.implementation_type == InverseMassType::GlobalKrylovSolver)
295 {
296 AssertThrow(global_solver.get() != 0,
297 dealii::ExcMessage("Global mass solver has not been initialized."));
298 global_solver->solve(dst, src);
299 }
300 else
301 {
302 dst.zero_out_ghost_values();
303
304 if(data.implementation_type == InverseMassType::MatrixfreeOperator)
305 {
306 matrix_free->cell_loop(&This::cell_loop_matrix_free_operator, this, dst, src);
307 }
308 else // ElementwiseKrylovSolver or BlockMatrices
309 {
310 AssertThrow(block_jacobi_preconditioner.get() != 0,
311 dealii::ExcMessage(
312 "Cell-wise iterative/direct block-Jacobi solver has not been initialized."));
313 block_jacobi_preconditioner->vmult(dst, src);
314 }
315 }
316 }
317
318 // dst = scaling_factor * (M^-1 * src)
319 void
320 apply_scale(VectorType & dst, double const scaling_factor, VectorType const & src) const
321 {
322 if(data.implementation_type == InverseMassType::MatrixfreeOperator)
323 {
324 // In the InverseMassType::MatrixfreeOperator case we can avoid
325 // streaming the vector from memory twice.
326
327 // ghost have to be zeroed out before MatrixFree::cell_loop().
328 dst.zero_out_ghost_values();
329
330 matrix_free->cell_loop(
331 &This::cell_loop_matrix_free_operator,
332 this,
333 dst,
334 src,
335 /*operation before cell operation*/ {}, /*operation after cell operation*/
336 [&](const unsigned int start_range, const unsigned int end_range) {
337 for(unsigned int i = start_range; i < end_range; ++i)
338 dst.local_element(i) *= scaling_factor;
339 },
340 dof_index);
341 }
342 else
343 {
344 apply(dst, src);
345 dst *= scaling_factor;
346 }
347 }
348
349
350private:
351 void
352 cell_loop_matrix_free_operator(dealii::MatrixFree<dim, Number> const &,
353 VectorType & dst,
354 VectorType const & src,
355 Range const & cell_range) const
356 {
357 Integrator integrator(*matrix_free, dof_index, quad_index);
358 InverseMassAsMatrixFreeOperator inverse_mass(integrator);
359
360 if(coefficient_is_variable)
361 {
362 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
363 {
364 integrator.reinit(cell);
365 integrator.read_dof_values(src, 0);
366
367 dealii::AlignedVector<dealii::VectorizedArray<Number>> inverse_JxW_times_coefficient(
368 integrator.n_q_points);
369 inverse_mass.fill_inverse_JxW_values(inverse_JxW_times_coefficient);
370
371 if(consider_inverse_coefficient)
372 {
373 // Consider a mass matrix of the form
374 // (u_h , v_h / c)_Omega
375 // hence fill the vector with (J / c)^-1 = c/J
376 for(unsigned int q = 0; q < integrator.n_q_points; ++q)
377 {
378 inverse_JxW_times_coefficient[q] *=
379 this->variable_coefficients->get_coefficient_cell(cell, q);
380 }
381 }
382 else
383 {
384 // Consider a mass matrix of the form
385 // (u_h , v_h * c)_Omega
386 // hence fill the vector with inv(J * c) = 1/(J * c)
387 for(unsigned int q = 0; q < integrator.n_q_points; ++q)
388 {
389 inverse_JxW_times_coefficient[q] /=
390 this->variable_coefficients->get_coefficient_cell(cell, q);
391 }
392 }
393
394 inverse_mass.apply(inverse_JxW_times_coefficient,
395 n_components,
396 integrator.begin_dof_values(),
397 integrator.begin_dof_values());
398
399 integrator.set_dof_values(dst, 0);
400 }
401 }
402 else
403 {
404 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
405 {
406 integrator.reinit(cell);
407 integrator.read_dof_values(src, 0);
408
409 inverse_mass.apply(integrator.begin_dof_values(), integrator.begin_dof_values());
410
411 integrator.set_dof_values(dst, 0);
412 }
413 }
414 }
415
416 dealii::MatrixFree<dim, Number> const * matrix_free;
417
418 unsigned int dof_index, quad_index;
419
420 InverseMassParameters data;
421
422 // Variable coefficients not managed by this class.
423 bool coefficient_is_variable;
424 bool consider_inverse_coefficient;
425
426 VariableCoefficients<dealii::VectorizedArray<Number>> const * variable_coefficients;
427
428 // Solver and preconditioner for solving a global linear system of equations for all degrees of
429 // freedom.
430 std::shared_ptr<PreconditionerBase<Number>> global_preconditioner;
431 std::shared_ptr<Krylov::SolverBase<VectorType>> global_solver;
432
433 // Block-Jacobi preconditioner used as cell-wise inverse for L2-conforming spaces. In this case,
434 // the mass matrix is block-diagonal and a block-Jacobi preconditioner inverts the mass operator
435 // exactly (up to solver tolerances). The implementation of the block-Jacobi preconditioner can be
436 // matrix-based or matrix-free, depending on the parameters specified.
437 std::shared_ptr<BlockJacobiPreconditioner<MassOperator<dim, n_components, Number>>>
438 block_jacobi_preconditioner;
439
440 // MassOperator as underlying operator for the cell-wise or global iterative solves.
441 MassOperator<dim, n_components, Number> mass_operator;
442};
443
444} // namespace ExaDG
445
446#endif /* EXADG_OPERATORS_INVERSE_MASS_OPERATOR_H_ */
void update()
Definition inverse_mass_operator.h:264
Definition preconditioner_base.h:35
Definition variable_coefficients.h:40
Definition driver.cpp:33
Definition inverse_mass_operator.h:43
Definition inverse_mass_parameters.h:54
Definition iterative_solvers_dealii_wrapper.h:92
Definition mass_operator.h:35