PRISMS-PF Manual v3.0-pre
All Classes Functions Variables Enumerations Pages
matrix_free_operator.h
1// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
2// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
3
4#ifndef matrix_free_operator_h
5#define matrix_free_operator_h
6
7#include <deal.II/base/subscriptor.h>
8#include <deal.II/base/vectorization.h>
9#include <deal.II/lac/la_parallel_vector.h>
10#include <deal.II/matrix_free/matrix_free.h>
11#include <deal.II/matrix_free/operators.h>
12
13#include <prismspf/config.h>
14#include <prismspf/core/type_enums.h>
15#include <prismspf/core/variable_attributes.h>
16#include <prismspf/core/variable_container.h>
17#include <prismspf/user_inputs/user_input_parameters.h>
18
19PRISMS_PF_BEGIN_NAMESPACE
20
29template <int dim, int degree, typename number>
30class matrixFreeOperator : public dealii::Subscriptor
31{
32public:
33 using VectorType = dealii::LinearAlgebra::distributed::Vector<number>;
34 using value_type = number;
35 using size_type = dealii::VectorizedArray<number>;
36
40 matrixFreeOperator() = default;
41
46 const std::map<unsigned int, variableAttributes> &_attributes_list);
47
52 const unsigned int &_current_index,
53 const std::map<unsigned int, variableAttributes> &_attributes_list);
54
58 void
59 initialize(std::shared_ptr<const dealii::MatrixFree<dim, number, size_type>> _data,
60 const std::vector<unsigned int> &selected_field_indexes =
61 std::vector<unsigned int>());
62
66 dealii::types::global_dof_index
67 m() const;
68
74 number
75 el(const unsigned int &row, const unsigned int &col) const;
76
81 void
82 clear();
83
88 void
89 initialize_dof_vector(VectorType &dst, unsigned int dof_handler_index = 0) const;
90
94 void
95 set_constrained_entries_to_one(VectorType &dst) const;
96
100 std::shared_ptr<const dealii::MatrixFree<dim, number, size_type>>
101 get_matrix_free() const;
102
106 const std::shared_ptr<dealii::DiagonalMatrix<VectorType>> &
108
112 void
114 const std::unordered_map<std::pair<unsigned int, dependencyType>,
115 unsigned int,
116 pairHash> &_global_to_local_solution);
117
121 void
123 const std::vector<VectorType *> &_src_solution_subset = std::vector<VectorType *>());
124
128 void
129 vmult(VectorType &dst, const VectorType &src) const;
130
134 void
135 Tvmult(VectorType &dst, const VectorType &src) const;
136
140 void
141 compute_explicit_update(std::vector<VectorType *> &dst,
142 const std::vector<VectorType *> &src) const;
143
147 void
148 compute_postprocess_explicit_update(std::vector<VectorType *> &dst,
149 const std::vector<VectorType *> &src) const;
150
154 void
155 compute_nonexplicit_auxiliary_update(std::vector<VectorType *> &dst,
156 const std::vector<VectorType *> &src) const;
157
161 void
162 compute_residual(VectorType &dst, const VectorType &src) const;
163
167 void
168 compute_diagonal(unsigned int field_index);
169
170protected:
174 virtual void
176 const dealii::Point<dim, size_type> &q_point_loc) const = 0;
177
181 virtual void
183 const dealii::Point<dim, size_type> &q_point_loc) const = 0;
184
188 virtual void
190 const dealii::Point<dim, size_type> &q_point_loc) const = 0;
191
195 virtual void
198 const dealii::Point<dim, size_type> &q_point_loc) const = 0;
199
204
208 const unsigned int current_index = numbers::invalid_index;
209
210private:
214 void
215 compute_local_explicit_update(
216 const dealii::MatrixFree<dim, number, size_type> &data,
217 std::vector<VectorType *> &dst,
218 const std::vector<VectorType *> &src,
219 const std::pair<unsigned int, unsigned int> &cell_range) const;
220
224 void
225 compute_local_postprocess_explicit_update(
226 const dealii::MatrixFree<dim, number, size_type> &data,
227 std::vector<VectorType *> &dst,
228 const std::vector<VectorType *> &src,
229 const std::pair<unsigned int, unsigned int> &cell_range) const;
230
234 void
235 compute_local_nonexplicit_auxiliary_update(
236 const dealii::MatrixFree<dim, number, size_type> &data,
237 std::vector<VectorType *> &dst,
238 const std::vector<VectorType *> &src,
239 const std::pair<unsigned int, unsigned int> &cell_range) const;
240
244 void
245 compute_local_residual(const dealii::MatrixFree<dim, number, size_type> &data,
246 VectorType &dst,
247 const VectorType &src,
248 const std::pair<unsigned int, unsigned int> &cell_range) const;
249
253 void
254 compute_local_newton_update(
255 const dealii::MatrixFree<dim, number, size_type> &data,
256 VectorType &dst,
257 const VectorType &src,
258 const std::pair<unsigned int, unsigned int> &cell_range) const;
259
263 void
264 local_compute_diagonal(const dealii::MatrixFree<dim, number, size_type> &data,
265 VectorType &dst,
266 const unsigned int &dummy,
267 const std::pair<unsigned int, unsigned int> &cell_range) const;
268
272 const std::map<unsigned int, variableAttributes> &attributes_list;
273
277 std::shared_ptr<const dealii::MatrixFree<dim, number, size_type>> data;
278
282 std::vector<unsigned int> selected_fields;
283
287 std::vector<std::vector<unsigned int>> edge_constrained_indices;
288
292 std::unordered_map<std::pair<unsigned int, dependencyType>, unsigned int, pairHash>
293 global_to_local_solution;
294
298 std::vector<VectorType *> src_solution_subset;
299
303 std::shared_ptr<dealii::DiagonalMatrix<VectorType>> diagonal_entries;
304
308 std::shared_ptr<dealii::DiagonalMatrix<VectorType>> inverse_diagonal_entries;
309};
310
311template <int dim, int degree, typename number>
313 const userInputParameters<dim> &_user_inputs,
314 const std::map<unsigned int, variableAttributes> &_attributes_list)
315 : Subscriptor()
316 , user_inputs(_user_inputs)
317 , attributes_list(_attributes_list)
318{}
319
320template <int dim, int degree, typename number>
322 const userInputParameters<dim> &_user_inputs,
323 const unsigned int &_current_index,
324 const std::map<unsigned int, variableAttributes> &_attributes_list)
325 : Subscriptor()
326 , user_inputs(_user_inputs)
327 , current_index(_current_index)
328 , attributes_list(_attributes_list)
329{}
330
331template <int dim, int degree, typename number>
332void
334 std::shared_ptr<const dealii::MatrixFree<dim, number, size_type>> _data,
335 const std::vector<unsigned int> &selected_field_indexes)
336{
337 data = _data;
338
339 selected_fields.clear();
340 if (selected_field_indexes.empty())
341 {
342 for (unsigned int i = 0; i < _data->n_components(); ++i)
343 {
344 selected_fields.push_back(i);
345 }
346 }
347 else
348 {
349 for (unsigned int i = 0; i < selected_field_indexes.size(); ++i)
350 {
351 AssertIndexRange(selected_field_indexes[i], _data->n_components());
352 for (unsigned int j = 0; j < selected_field_indexes.size(); ++j)
353 {
354 if (j != i)
355 {
356 Assert(selected_field_indexes[j] != selected_field_indexes[i],
357 dealii::ExcMessage("Given row indices must be unique"));
358 }
359 }
360 selected_fields.push_back(selected_field_indexes[i]);
361 }
362 }
363
364 edge_constrained_indices.clear();
365 edge_constrained_indices.resize(selected_fields.size());
366}
367
368template <int dim, int degree, typename number>
369dealii::types::global_dof_index
371{
372 Assert(data.get() != nullptr, dealii::ExcNotInitialized());
373 unsigned int total_size = 0;
374 for (const unsigned int field : selected_fields)
375 {
376 total_size += data->get_vector_partitioner(field)->size();
377 }
378 return total_size;
379}
380
381template <int dim, int degree, typename number>
382number
384 [[maybe_unused]] const unsigned int &row,
385 [[maybe_unused]] const unsigned int &col) const
386{
387 AssertThrow(false, FeatureNotImplemented("el()"));
388 return 0.0;
389}
390
391template <int dim, int degree, typename number>
392void
394{
395 data.reset();
396 inverse_diagonal_entries.reset();
397 global_to_local_solution.clear();
398}
399
400template <int dim, int degree, typename number>
401void
403 VectorType &dst,
404 unsigned int dof_handler_index) const
405{
406 data->initialize_dof_vector(dst, dof_handler_index);
407}
408
409template <int dim, int degree, typename number>
410void
412 VectorType &dst) const
413{
414 for (unsigned int j = 0; j < dealii::MatrixFreeOperators::BlockHelper::n_blocks(dst);
415 ++j)
416 {
417 const std::vector<unsigned int> &constrained_dofs =
418 data->get_constrained_dofs(selected_fields[j]);
419 for (const auto constrained_dof : constrained_dofs)
420 {
421 dealii::MatrixFreeOperators::BlockHelper::subblock(dst, j).local_element(
422 constrained_dof) = 1.0;
423 }
424 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
425 {
426 dealii::MatrixFreeOperators::BlockHelper::subblock(dst, j).local_element(
427 edge_constrained_indices[j][i]) = 1.0;
428 }
429 }
430}
431
432template <int dim, int degree, typename number>
433void
435 const std::unordered_map<std::pair<unsigned int, dependencyType>,
436 unsigned int,
437 pairHash> &_global_to_local_solution)
438{
439 global_to_local_solution = _global_to_local_solution;
440}
441
442template <int dim, int degree, typename number>
443std::shared_ptr<const dealii::MatrixFree<dim, number>>
448
449template <int dim, int degree, typename number>
450const std::shared_ptr<
451 dealii::DiagonalMatrix<typename matrixFreeOperator<dim, degree, number>::VectorType>> &
453{
454 Assert(inverse_diagonal_entries.get() != nullptr && inverse_diagonal_entries->m() > 0,
455 dealii::ExcNotInitialized());
456 return inverse_diagonal_entries;
457}
458
459template <int dim, int degree, typename number>
460void
462 const std::vector<VectorType *> &_src_solution_subset)
463{
464 src_solution_subset = _src_solution_subset;
465}
466
467template <int dim, int degree, typename number>
468void
470 std::vector<VectorType *> &dst,
471 const std::vector<VectorType *> &src) const
472{
473 Assert(!global_to_local_solution.empty(),
474 dealii::ExcMessage(
475 "The global to local solution mapping must not be empty. Make sure to call "
476 "add_global_to_local_mapping() prior to any computations."));
477 Assert(!dst.empty(), dealii::ExcMessage("The dst vector must not be empty"));
478 Assert(!src.empty(), dealii::ExcMessage("The src vector must not be empty"));
479
480 this->data->cell_loop(&matrixFreeOperator::compute_local_explicit_update,
481 this,
482 dst,
483 src,
484 true);
485}
486
487template <int dim, int degree, typename number>
488void
490 std::vector<VectorType *> &dst,
491 const std::vector<VectorType *> &src) const
492{
493 Assert(!global_to_local_solution.empty(),
494 dealii::ExcMessage(
495 "The global to local solution mapping must not be empty. Make sure to call "
496 "add_global_to_local_mapping() prior to any computations."));
497 Assert(!dst.empty(), dealii::ExcMessage("The dst vector must not be empty"));
498 Assert(!src.empty(), dealii::ExcMessage("The src vector must not be empty"));
499
500 this->data->cell_loop(&matrixFreeOperator::compute_local_postprocess_explicit_update,
501 this,
502 dst,
503 src,
504 true);
505}
506
507template <int dim, int degree, typename number>
508void
510 std::vector<VectorType *> &dst,
511 const std::vector<VectorType *> &src) const
512{
513 Assert(!global_to_local_solution.empty(),
514 dealii::ExcMessage(
515 "The global to local solution mapping must not be empty. Make sure to call "
516 "add_global_to_local_mapping() prior to any computations."));
517 Assert(!dst.empty(), dealii::ExcMessage("The dst vector must not be empty"));
518 Assert(!src.empty(), dealii::ExcMessage("The src vector must not be empty"));
519
520 this->data->cell_loop(&matrixFreeOperator::compute_local_nonexplicit_auxiliary_update,
521 this,
522 dst,
523 src,
524 true);
525}
526
527template <int dim, int degree, typename number>
528void
530 const VectorType &src) const
531{
532 Assert(!global_to_local_solution.empty(),
533 dealii::ExcMessage(
534 "The global to local solution mapping must not be empty. Make sure to call "
535 "add_global_to_local_mapping() prior to any computations."));
536 Assert(!src_solution_subset.empty(),
537 dealii::ExcMessage("The src_solution_subset vector must not be empty"));
538 Assert(dst.size() != 0,
539 dealii::ExcMessage("The dst vector should not have size equal to 0"));
540 Assert(src.size() != 0,
541 dealii::ExcMessage("The src vector should not have size equal to 0"));
542
543 this->data->cell_loop(&matrixFreeOperator::compute_local_residual,
544 this,
545 dst,
546 src,
547 true);
548}
549
550template <int dim, int degree, typename number>
551void
553 const VectorType &src) const
554{
555 Assert(!global_to_local_solution.empty(),
556 dealii::ExcMessage(
557 "The global to local solution mapping must not be empty. Make sure to call "
558 "add_global_to_local_mapping() prior to any computations."));
559 Assert(dst.size() != 0,
560 dealii::ExcMessage("The dst vector should not have size equal to 0"));
561 Assert(src.size() != 0,
562 dealii::ExcMessage("The src vector should not have size equal to 0"));
563
564 this->data->cell_loop(&matrixFreeOperator::compute_local_newton_update,
565 this,
566 dst,
567 src,
568 true);
569}
570
571template <int dim, int degree, typename number>
572void
574 const VectorType &src) const
575{
576 this->vmult(dst, src);
577}
578
579template <int dim, int degree, typename number>
580void
582 const dealii::MatrixFree<dim, number> &data,
583 std::vector<VectorType *> &dst,
584 const std::vector<VectorType *> &src,
585 const std::pair<unsigned int, unsigned int> &cell_range) const
586{
587 // Constructor for FEEvaluation objects
589 attributes_list,
590 global_to_local_solution,
591 solveType::EXPLICIT_RHS);
592
593 // Initialize, evaluate, and submit based on user function.
594 variable_list.eval_local_operator(
596 const dealii::Point<dim, size_type> &q_point_loc)
597 {
598 this->compute_explicit_RHS(var_list, q_point_loc);
599 },
600 dst,
601 src,
602 cell_range);
603}
604
605template <int dim, int degree, typename number>
606void
608 const dealii::MatrixFree<dim, number> &data,
609 std::vector<VectorType *> &dst,
610 const std::vector<VectorType *> &src,
611 const std::pair<unsigned int, unsigned int> &cell_range) const
612{
613 // Constructor for FEEvaluation objects
615 attributes_list,
616 global_to_local_solution,
617 solveType::POSTPROCESS);
618
619 // Initialize, evaluate, and submit based on user function.
620 variable_list.eval_local_operator(
622 const dealii::Point<dim, size_type> &q_point_loc)
623 {
624 this->compute_postprocess_explicit_RHS(var_list, q_point_loc);
625 },
626 dst,
627 src,
628 cell_range);
629}
630
631template <int dim, int degree, typename number>
632void
634 const dealii::MatrixFree<dim, number> &data,
635 std::vector<VectorType *> &dst,
636 const std::vector<VectorType *> &src,
637 const std::pair<unsigned int, unsigned int> &cell_range) const
638{
639 // Constructor for FEEvaluation objects
641 attributes_list,
642 global_to_local_solution,
643 solveType::NONEXPLICIT_RHS);
644
645 // Initialize, evaluate, and submit based on user function.
646 variable_list.eval_local_operator(
648 const dealii::Point<dim, size_type> &q_point_loc)
649 {
650 this->compute_nonexplicit_RHS(var_list, q_point_loc);
651 },
652 dst,
653 src,
654 cell_range);
655}
656
657template <int dim, int degree, typename number>
658void
660 const dealii::MatrixFree<dim, number> &data,
661 VectorType &dst,
662 [[maybe_unused]] const VectorType &src,
663 const std::pair<unsigned int, unsigned int> &cell_range) const
664{
665 // Constructor for FEEvaluation objects
667 attributes_list,
668 global_to_local_solution,
669 solveType::NONEXPLICIT_RHS);
670
671 // Initialize, evaluate, and submit based on user function.
672 variable_list.eval_local_operator(
674 const dealii::Point<dim, size_type> &q_point_loc)
675 {
676 this->compute_nonexplicit_RHS(var_list, q_point_loc);
677 },
678 dst,
679 src_solution_subset,
680 cell_range);
681}
682
683template <int dim, int degree, typename number>
684void
686 const dealii::MatrixFree<dim, number> &data,
687 VectorType &dst,
688 const VectorType &src,
689 const std::pair<unsigned int, unsigned int> &cell_range) const
690{
691 // Constructor for FEEvaluation objects
693 attributes_list,
694 global_to_local_solution,
695 solveType::NONEXPLICIT_LHS);
696
697 // Initialize, evaluate, and submit based on user function. Note that the src solution
698 // subset must not include the src vector.
699 variable_list.eval_local_operator(
701 const dealii::Point<dim, size_type> &q_point_loc)
702 {
703 this->compute_nonexplicit_LHS(var_list, q_point_loc);
704 },
705 dst,
706 src,
707 src_solution_subset,
708 cell_range);
709}
710
711template <int dim, int degree, typename number>
712void
714{
715 inverse_diagonal_entries.reset(new dealii::DiagonalMatrix<VectorType>());
716 VectorType &inverse_diagonal = inverse_diagonal_entries->get_vector();
717 data->initialize_dof_vector(inverse_diagonal, field_index);
718 unsigned int dummy = 0;
719 data->cell_loop(&matrixFreeOperator::local_compute_diagonal,
720 this,
721 inverse_diagonal,
722 dummy);
723
724 set_constrained_entries_to_one(inverse_diagonal);
725
726 for (unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i)
727 {
728 Assert(inverse_diagonal.local_element(i) > 0.0,
729 dealii::ExcMessage(
730 "No diagonal entry in a positive definite operator should be zero"));
731 inverse_diagonal.local_element(i) = 1.0 / inverse_diagonal.local_element(i);
732 }
733}
734
735template <int dim, int degree, typename number>
736void
738 const dealii::MatrixFree<dim, number> &data,
739 VectorType &dst,
740 [[maybe_unused]] const unsigned int &dummy,
741 const std::pair<unsigned int, unsigned int> &cell_range) const
742{
743 // Constructor for FEEvaluation objects
745 attributes_list,
746 global_to_local_solution,
747 solveType::NONEXPLICIT_LHS);
748
749 // Initialize, evaluate, and submit diagonal based on user function.
750 variable_list.eval_local_diagonal(
752 const dealii::Point<dim, size_type> &q_point_loc)
753 {
754 this->compute_nonexplicit_LHS(var_list, q_point_loc);
755 },
756 dst,
757 src_solution_subset,
758 cell_range);
759}
760
761PRISMS_PF_END_NAMESPACE
762
763#endif
This is the abstract base class for the matrix-free implementation of some PDE.
Definition matrix_free_operator.h:31
void compute_postprocess_explicit_update(std::vector< VectorType * > &dst, const std::vector< VectorType * > &src) const
Compute the explicit update for postprocessed fields.
Definition matrix_free_operator.h:489
void compute_explicit_update(std::vector< VectorType * > &dst, const std::vector< VectorType * > &src) const
Compute the explicit update.
Definition matrix_free_operator.h:469
void compute_residual(VectorType &dst, const VectorType &src) const
Compute the residual of this operator. This is the b in Ax=b.
Definition matrix_free_operator.h:529
virtual void compute_postprocess_explicit_RHS(variableContainer< dim, degree, number > &variable_list, const dealii::Point< dim, size_type > &q_point_loc) const =0
User-implemented class for the RHS of postprocessed explicit equations.
matrixFreeOperator()=default
Default constructor.
void compute_nonexplicit_auxiliary_update(std::vector< VectorType * > &dst, const std::vector< VectorType * > &src) const
Compute a nonexplicit auxiliary update.
Definition matrix_free_operator.h:509
const userInputParameters< dim > & user_inputs
The user-inputs.
Definition matrix_free_operator.h:203
void initialize(std::shared_ptr< const dealii::MatrixFree< dim, number, size_type > > _data, const std::vector< unsigned int > &selected_field_indexes=std::vector< unsigned int >())
Initialize operator.
Definition matrix_free_operator.h:333
virtual void compute_nonexplicit_RHS(variableContainer< dim, degree, number > &variable_list, const dealii::Point< dim, size_type > &q_point_loc) const =0
User-implemented class for the RHS of nonexplicit equations.
number el(const unsigned int &row, const unsigned int &col) const
Return the value of the matrix entry. This function is only valid when row == col and when the diagon...
Definition matrix_free_operator.h:383
const unsigned int current_index
The current index that is being solved.
Definition matrix_free_operator.h:208
void set_constrained_entries_to_one(VectorType &dst) const
Set constrained entries to one.
Definition matrix_free_operator.h:411
void compute_diagonal(unsigned int field_index)
Compute the diagonal of this operator.
Definition matrix_free_operator.h:713
dealii::types::global_dof_index m() const
Return the number of DoFs.
Definition matrix_free_operator.h:370
void add_src_solution_subset(const std::vector< VectorType * > &_src_solution_subset=std::vector< VectorType * >())
Add the solution subset for src vector.
Definition matrix_free_operator.h:461
void clear()
Release all memory and return to state like having called the default constructor.
Definition matrix_free_operator.h:393
void Tvmult(VectorType &dst, const VectorType &src) const
Transpose matrix-vector multiplication.
Definition matrix_free_operator.h:573
std::shared_ptr< const dealii::MatrixFree< dim, number, size_type > > get_matrix_free() const
Get read access to the MatrixFree object stored with this operator.
Definition matrix_free_operator.h:444
virtual void compute_nonexplicit_LHS(variableContainer< dim, degree, number > &variable_list, const dealii::Point< dim, size_type > &q_point_loc) const =0
User-implemented class for the LHS of nonexplicit equations.
const std::shared_ptr< dealii::DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Get read access to the inverse diagonal of this operator.
Definition matrix_free_operator.h:452
void add_global_to_local_mapping(const std::unordered_map< std::pair< unsigned int, dependencyType >, unsigned int, pairHash > &_global_to_local_solution)
Add the mappings from global to local solution vectors.
Definition matrix_free_operator.h:434
void initialize_dof_vector(VectorType &dst, unsigned int dof_handler_index=0) const
Initialize a given vector with the MatrixFree object that this object contains.
Definition matrix_free_operator.h:402
virtual void compute_explicit_RHS(variableContainer< dim, degree, number > &variable_list, const dealii::Point< dim, size_type > &q_point_loc) const =0
User-implemented class for the RHS of explicit equations.
void vmult(VectorType &dst, const VectorType &src) const
Matrix-vector multiplication.
Definition matrix_free_operator.h:552
Definition user_input_parameters.h:22
This class permits the access of a subset of indexed fields and gives an error if any non-allowed fie...
Definition variable_container.h:30
Simple hash function for pairs.
Definition variable_attributes.h:24