include/meltpooldg/utilities/matrix_free_util.hpp Source File

Developer Documentation: include/meltpooldg/utilities/matrix_free_util.hpp Source File
Developer Documentation
matrix_free_util.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <deal.II/base/tensor.h>
4
5#include <deal.II/grid/grid_tools.h>
6#include <deal.II/grid/tria.h>
7
8#include <deal.II/matrix_free/fe_evaluation.h>
9#include <deal.II/matrix_free/matrix_free.h>
10#include <deal.II/matrix_free/tools.h>
11
12#include <meltpooldg/utilities/fe_integrator.hpp>
13
14#include <vector>
15
16namespace MeltPoolDG
17{
18
19 template <int dim, typename Number>
20 struct MatrixFreeContext
21 {
22 const dealii::MatrixFree<dim, Number> &mf;
23 unsigned int dof_idx;
24 unsigned int quad_idx;
25 };
26
43 template <typename FeEval>
44 dealii::Tensor<1, FeEval::n_components, dealii::VectorizedArray<typename FeEval::number_type>>
45 fe_evaluation_tensor_value_at_q(const FeEval &fe_eval, const unsigned q_index)
46 {
47 using ValueType = typename FeEval::value_type;
48 using VectorizedArrayType = dealii::VectorizedArray<typename FeEval::number_type>;
49
50 if constexpr (std::is_same_v<ValueType, VectorizedArrayType>)
51 {
52 dealii::Tensor<1, 1, VectorizedArrayType> t;
53 t[0] = fe_eval.get_value(q_index);
54 return t;
55 }
56 else
57 return fe_eval.get_value(q_index);
58 }
59
67 template <int dim, typename number>
68 std::vector<dealii::TriaIterator<dealii::CellAccessor<dim>>>
69 cells_in_cell_batch(const dealii::MatrixFree<dim, number> &mf, const unsigned int cell_batch_id)
70 {
71 unsigned int n_active_lanes = mf.n_active_entries_per_cell_batch(cell_batch_id);
72
73 std::vector<dealii::TriaIterator<dealii::CellAccessor<dim>>> cells;
74 cells.reserve(n_active_lanes);
75
76 for (unsigned int lane = 0; lane < n_active_lanes; ++lane)
77 cells.push_back(mf.get_cell_iterator(cell_batch_id, lane));
78
79 return cells;
80 }
81
89
95 template <typename OperatorType, int dim, typename number>
96 concept HasJacobianQuadraturePointKernels =
97 requires(OperatorType op,
98 dealii::FEEvaluation<dim, -1, 0, dim + 2, number> &fe_evaluator,
99 const dealii::FEEvaluation<dim, -1, 0, dim + 2, number> &const_fe_evaluator,
100 dealii::FEFaceEvaluation<dim, -1, 0, dim + 2, number> &fe_face_evaluator,
101 const dealii::FEFaceEvaluation<dim, -1, 0, dim + 2, number> &const_fe_face_evaluator,
102 unsigned int q_index) {
107 op.local_cell_jacobian_kernel(fe_evaluator, const_fe_evaluator, q_index);
108
113 op.local_face_jacobian_kernel(fe_face_evaluator,
114 fe_face_evaluator,
115 const_fe_face_evaluator,
116 const_fe_face_evaluator,
117 q_index);
118
123 op.local_boundary_face_jacobian_kernel(fe_face_evaluator, const_fe_face_evaluator, q_index);
124 };
125
141 template <int dim,
142 int n_components,
143 typename number,
144 HasJacobianQuadraturePointKernels<dim, number> OperatorType>
145 void
146 compute_jacobian_matrix_representation(
147 const OperatorType &implicit_operator,
148 std::variant<dealii::LinearAlgebra::distributed::Vector<number> *,
149 dealii::TrilinosWrappers::SparseMatrix *> dst,
150 const MatrixRepresentationType type,
151 const dealii::LinearAlgebra::distributed::Vector<number> &current_solution,
152 const dealii::MatrixFree<dim, number> &matrix_free,
153 const unsigned int dof_idx,
154 const unsigned int quad_idx)
155 {
156 FECellIntegrator<dim, n_components, number> phi(matrix_free, dof_idx, quad_idx);
157
158 // ## Cell ## //
159 std::function<void(FECellIntegrator<dim, n_components, number> &)> cell_operation =
160 [&](FECellIntegrator<dim, n_components, number> &delta_phi) -> void {
161 phi.reinit(delta_phi.get_cell_or_face_batch_id());
162 phi.gather_evaluate(current_solution,
163 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
164 delta_phi.evaluate(dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
165
166 for (const unsigned int q_index : delta_phi.quadrature_point_indices())
167 {
168 implicit_operator.local_cell_jacobian_kernel(delta_phi, phi, q_index);
169 }
170
171 delta_phi.integrate(dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
172 };
173
174 // ## Face ##//
175 FEFaceIntegrator<dim, n_components, number> phi_m(matrix_free,
176 true /*is_interior_face*/,
177 dof_idx,
178 quad_idx);
179 FEFaceIntegrator<dim, n_components, number> phi_p(matrix_free,
180 false /*is_interior_face*/,
181 dof_idx,
182 quad_idx);
183 std::function<void(FEFaceIntegrator<dim, n_components, number> &,
184 FEFaceIntegrator<dim, n_components, number> &)>
185 face_operation = [&](FEFaceIntegrator<dim, n_components, number> &delta_phi_m,
186 FEFaceIntegrator<dim, n_components, number> &delta_phi_p) -> void {
187 phi_m.reinit(delta_phi_m.get_cell_or_face_batch_id());
188 phi_p.reinit(delta_phi_p.get_cell_or_face_batch_id());
189 phi_m.gather_evaluate(current_solution,
190 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
191 phi_p.gather_evaluate(current_solution,
192 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
193
194 delta_phi_p.evaluate(dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
195 delta_phi_m.evaluate(dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
196
197
198 for (const unsigned int q_index : delta_phi_m.quadrature_point_indices())
199 {
200 implicit_operator.local_face_jacobian_kernel(
201 delta_phi_m, delta_phi_p, phi_m, phi_p, q_index);
202 }
203 delta_phi_m.integrate(dealii::EvaluationFlags::values);
204 delta_phi_p.integrate(dealii::EvaluationFlags::values);
205 };
206
207 // ## Boundary face ##//
208 FEFaceIntegrator<dim, n_components, number> phi_boundary(matrix_free,
209 true /*is_interior_face*/,
210 dof_idx,
211 quad_idx);
212 std::function<void(FEFaceIntegrator<dim, n_components, number> &)> boundary_operation =
213 [&](FEFaceIntegrator<dim, n_components, number> &delta_phi_boundary) -> void {
214 phi_boundary.reinit(delta_phi_boundary.get_cell_or_face_batch_id());
215 phi_boundary.gather_evaluate(current_solution,
216 dealii::EvaluationFlags::values |
217 dealii::EvaluationFlags::gradients);
218 delta_phi_boundary.evaluate(dealii::EvaluationFlags::values |
219 dealii::EvaluationFlags::gradients);
220
221
222 for (const unsigned int q_index : delta_phi_boundary.quadrature_point_indices())
223 {
224 implicit_operator.local_boundary_face_jacobian_kernel(delta_phi_boundary,
225 phi_boundary,
226 q_index);
227 }
228 delta_phi_boundary.integrate(dealii::EvaluationFlags::values);
229 };
230
231 switch (type)
232 {
233 case MatrixRepresentationType::DiagonalMatrix: {
234 dealii::MatrixFreeTools::compute_diagonal(
235 matrix_free,
236 *std::get<dealii::LinearAlgebra::distributed::Vector<number> *>(dst),
237 cell_operation,
238 face_operation,
239 boundary_operation,
240 dof_idx,
241 quad_idx);
242 break;
243 }
244 case MatrixRepresentationType::SystemMatrix: {
245 dealii::MatrixFreeTools::compute_matrix(
246 matrix_free,
247 dealii::AffineConstraints<number>(),
248 *std::get<dealii::TrilinosWrappers::SparseMatrix *>(dst),
249 cell_operation,
250 face_operation,
251 boundary_operation,
252 dof_idx,
253 quad_idx);
254 break;
255 }
256 default:
257 DEAL_II_ASSERT_UNREACHABLE();
258 }
259 }
260
276 template <int dim, typename number, int n_components>
277 std::vector<std::pair<dealii::Tensor<1, n_components, number>,
278 dealii::Tensor<1, n_components, dealii::Tensor<1, dim, number>>>>
279 compute_cell_average_quantities(
280 const MatrixFreeContext<dim, number> &mf_context,
281 const dealii::LinearAlgebra::distributed::Vector<number> &solution)
282 {
283 using VectorizedArrayType = dealii::VectorizedArray<number>;
284 using ActiveCellIterator = typename dealii::Triangulation<dim>::active_cell_iterator;
285
286 using AveragePairType =
287 std::pair<dealii::Tensor<1, n_components, number>,
288 dealii::Tensor<1, n_components, dealii::Tensor<1, dim, number>>>;
289
290 // A lambda used to exchange the cell average values of ghost cells after they have been
291 // computed for the locally owned cells.
292 const auto exchange_cell_data_to_ghosts =
293 [](const dealii::Triangulation<dim> &triangulation,
294 std::vector<AveragePairType> &cell_average_values) {
295 std::function<std::optional<AveragePairType>(const ActiveCellIterator &)> pack =
296 [&](const ActiveCellIterator &cell) -> std::optional<AveragePairType> {
297 return cell_average_values[cell->active_cell_index()];
298 };
299 std::function<void(const ActiveCellIterator &, const AveragePairType &)> unpack =
300 [&](const ActiveCellIterator &cell, const AveragePairType &value) {
301 cell_average_values[cell->active_cell_index()] = value;
302 };
303 dealii::GridTools::exchange_cell_data_to_ghosts(triangulation, pack, unpack);
304 };
305
306 // A lambda that computes the cell average values for a given range of cells used in the cell
307 // loop of the matrix free object.
308 const std::function<void(const dealii::MatrixFree<dim, number, VectorizedArrayType> &,
309 std::vector<AveragePairType> &,
310 const dealii::LinearAlgebra::distributed::Vector<number> &,
311 const std::pair<unsigned int, unsigned int> &)>
312 cell_loop = [&](const dealii::MatrixFree<dim, number, VectorizedArrayType> &matrix_free,
313 std::vector<AveragePairType> &dst_cell_average_values,
314 const dealii::LinearAlgebra::distributed::Vector<number> &src,
315 const std::pair<unsigned int, unsigned int> &cell_range) {
316 FECellIntegrator<dim, n_components, number> fe_cell_integrator(matrix_free,
317 mf_context.dof_idx,
318 mf_context.quad_idx);
319
320 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
321 {
322 fe_cell_integrator.reinit(cell);
323 fe_cell_integrator.gather_evaluate(src,
324 dealii::EvaluationFlags::values |
325 dealii::EvaluationFlags::gradients);
326
327 typename FECellIntegrator<dim, n_components, number>::value_type cell_average_value;
328 typename FECellIntegrator<dim, n_components, number>::gradient_type
329 cell_average_gradient;
330 VectorizedArrayType cell_volume = 0.;
331 for (const unsigned int q : fe_cell_integrator.quadrature_point_indices())
332 {
333 cell_average_value += fe_cell_integrator.get_value(q) * fe_cell_integrator.JxW(q);
334 cell_average_gradient +=
335 fe_cell_integrator.get_gradient(q) * fe_cell_integrator.JxW(q);
336 cell_volume += fe_cell_integrator.JxW(q);
337 }
338 cell_average_value /= cell_volume;
339 for (unsigned int c = 0; c < n_components; ++c)
340 cell_average_gradient[c] /= cell_volume;
341
342 const auto cells = cells_in_cell_batch(mf_context.mf, cell);
343
344 for (unsigned int lane = 0; lane < mf_context.mf.n_active_entries_per_cell_batch(cell);
345 ++lane)
346 {
347 for (unsigned int c = 0; c < n_components; ++c)
348 {
349 dst_cell_average_values[cells[lane]->active_cell_index()].first[c] =
350 cell_average_value[c][lane];
351 for (unsigned int d = 0; d < dim; ++d)
352 {
353 dst_cell_average_values[cells[lane]->active_cell_index()].second[c][d] =
354 cell_average_gradient[c][d][lane];
355 }
356 }
357 }
358 }
359 };
360
361 std::vector<AveragePairType> cell_average_values(
362 mf_context.mf.get_dof_handler(mf_context.dof_idx).get_triangulation().n_active_cells());
363 mf_context.mf.cell_loop(cell_loop, cell_average_values, solution);
364
365 exchange_cell_data_to_ghosts(
366 mf_context.mf.get_dof_handler(mf_context.dof_idx).get_triangulation(), cell_average_values);
367
368 return cell_average_values;
369 }
370} // namespace MeltPoolDG
MeltPoolDG::HasJacobianQuadraturePointKernels
Definition matrix_free_util.hpp:96
fe_integrator.hpp
MeltPoolDG
Interface for a general preconditioner.
Definition boundary_condition_functions.hpp:17
MeltPoolDG::fe_evaluation_tensor_value_at_q
dealii::Tensor< 1, FeEval::n_components, dealii::VectorizedArray< typename FeEval::number_type > > fe_evaluation_tensor_value_at_q(const FeEval &fe_eval, const unsigned q_index)
Definition matrix_free_util.hpp:45
MeltPoolDG::FEFaceIntegrator
dealii::FEFaceEvaluation< dim, -1, 0, n_components, number, VectorizedArrayType > FEFaceIntegrator
Definition fe_integrator.hpp:22
MeltPoolDG::compute_cell_average_quantities
std::vector< std::pair< dealii::Tensor< 1, n_components, number >, dealii::Tensor< 1, n_components, dealii::Tensor< 1, dim, number > > > > compute_cell_average_quantities(const MatrixFreeContext< dim, number > &mf_context, const dealii::LinearAlgebra::distributed::Vector< number > &solution)
Definition matrix_free_util.hpp:279
MeltPoolDG::cells_in_cell_batch
std::vector< dealii::TriaIterator< dealii::CellAccessor< dim > > > cells_in_cell_batch(const dealii::MatrixFree< dim, number > &mf, const unsigned int cell_batch_id)
Definition matrix_free_util.hpp:69
MeltPoolDG::compute_jacobian_matrix_representation
void compute_jacobian_matrix_representation(const OperatorType &implicit_operator, std::variant< dealii::LinearAlgebra::distributed::Vector< number > *, dealii::TrilinosWrappers::SparseMatrix * > dst, const MatrixRepresentationType type, const dealii::LinearAlgebra::distributed::Vector< number > &current_solution, const dealii::MatrixFree< dim, number > &matrix_free, const unsigned int dof_idx, const unsigned int quad_idx)
Definition matrix_free_util.hpp:146
MeltPoolDG::MatrixRepresentationType
MatrixRepresentationType
Definition matrix_free_util.hpp:85
MeltPoolDG::DiagonalMatrix
@ DiagonalMatrix
Definition matrix_free_util.hpp:86
MeltPoolDG::SystemMatrix
@ SystemMatrix
Definition matrix_free_util.hpp:87
MeltPoolDG::FECellIntegrator
dealii::FEEvaluation< dim, -1, 0, n_components, number, VectorizedArrayType > FECellIntegrator
Definition fe_integrator.hpp:14
MeltPoolDG::MatrixFreeContext
Definition matrix_free_util.hpp:21
MeltPoolDG::MatrixFreeContext::quad_idx
unsigned int quad_idx
Definition matrix_free_util.hpp:24
MeltPoolDG::MatrixFreeContext::mf
const dealii::MatrixFree< dim, Number > & mf
Definition matrix_free_util.hpp:22
MeltPoolDG::MatrixFreeContext::dof_idx
unsigned int dof_idx
Definition matrix_free_util.hpp:23
Generated by doxygen 1.9.8