4#ifndef matrix_free_operator_h
5#define matrix_free_operator_h
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>
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>
19PRISMS_PF_BEGIN_NAMESPACE
29template <
int dim,
int degree,
typename number>
33 using VectorType = dealii::LinearAlgebra::distributed::Vector<number>;
34 using value_type = number;
35 using size_type = dealii::VectorizedArray<number>;
46 const std::map<unsigned int, variableAttributes> &_attributes_list);
52 const unsigned int &_current_index,
53 const std::map<unsigned int, variableAttributes> &_attributes_list);
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>());
66 dealii::types::global_dof_index
75 el(
const unsigned int &row,
const unsigned int &col)
const;
100 std::shared_ptr<const dealii::MatrixFree<dim, number, size_type>>
106 const std::shared_ptr<dealii::DiagonalMatrix<VectorType>> &
114 const std::unordered_map<std::pair<unsigned int, dependencyType>,
116 pairHash> &_global_to_local_solution);
123 const std::vector<VectorType *> &_src_solution_subset = std::vector<VectorType *>());
129 vmult(VectorType &dst,
const VectorType &src)
const;
135 Tvmult(VectorType &dst,
const VectorType &src)
const;
142 const std::vector<VectorType *> &src)
const;
149 const std::vector<VectorType *> &src)
const;
156 const std::vector<VectorType *> &src)
const;
176 const dealii::Point<dim, size_type> &q_point_loc)
const = 0;
183 const dealii::Point<dim, size_type> &q_point_loc)
const = 0;
190 const dealii::Point<dim, size_type> &q_point_loc)
const = 0;
198 const dealii::Point<dim, size_type> &q_point_loc)
const = 0;
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;
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;
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;
245 compute_local_residual(
const dealii::MatrixFree<dim, number, size_type> &data,
247 const VectorType &src,
248 const std::pair<unsigned int, unsigned int> &cell_range)
const;
254 compute_local_newton_update(
255 const dealii::MatrixFree<dim, number, size_type> &data,
257 const VectorType &src,
258 const std::pair<unsigned int, unsigned int> &cell_range)
const;
264 local_compute_diagonal(
const dealii::MatrixFree<dim, number, size_type> &data,
266 const unsigned int &dummy,
267 const std::pair<unsigned int, unsigned int> &cell_range)
const;
272 const std::map<unsigned int, variableAttributes> &attributes_list;
277 std::shared_ptr<const dealii::MatrixFree<dim, number, size_type>> data;
282 std::vector<unsigned int> selected_fields;
287 std::vector<std::vector<unsigned int>> edge_constrained_indices;
292 std::unordered_map<std::pair<unsigned int, dependencyType>,
unsigned int,
pairHash>
293 global_to_local_solution;
298 std::vector<VectorType *> src_solution_subset;
303 std::shared_ptr<dealii::DiagonalMatrix<VectorType>> diagonal_entries;
308 std::shared_ptr<dealii::DiagonalMatrix<VectorType>> inverse_diagonal_entries;
311template <
int dim,
int degree,
typename number>
314 const std::map<unsigned int, variableAttributes> &_attributes_list)
316 , user_inputs(_user_inputs)
317 , attributes_list(_attributes_list)
320template <
int dim,
int degree,
typename number>
323 const unsigned int &_current_index,
324 const std::map<unsigned int, variableAttributes> &_attributes_list)
326 , user_inputs(_user_inputs)
327 , current_index(_current_index)
328 , attributes_list(_attributes_list)
331template <
int dim,
int degree,
typename number>
334 std::shared_ptr<
const dealii::MatrixFree<dim, number, size_type>> _data,
335 const std::vector<unsigned int> &selected_field_indexes)
339 selected_fields.clear();
340 if (selected_field_indexes.empty())
342 for (
unsigned int i = 0; i < _data->n_components(); ++i)
344 selected_fields.push_back(i);
349 for (
unsigned int i = 0; i < selected_field_indexes.size(); ++i)
351 AssertIndexRange(selected_field_indexes[i], _data->n_components());
352 for (
unsigned int j = 0; j < selected_field_indexes.size(); ++j)
356 Assert(selected_field_indexes[j] != selected_field_indexes[i],
357 dealii::ExcMessage(
"Given row indices must be unique"));
360 selected_fields.push_back(selected_field_indexes[i]);
364 edge_constrained_indices.clear();
365 edge_constrained_indices.resize(selected_fields.size());
368template <
int dim,
int degree,
typename number>
369dealii::types::global_dof_index
372 Assert(data.get() !=
nullptr, dealii::ExcNotInitialized());
373 unsigned int total_size = 0;
374 for (
const unsigned int field : selected_fields)
376 total_size += data->get_vector_partitioner(field)->size();
381template <
int dim,
int degree,
typename number>
384 [[maybe_unused]]
const unsigned int &row,
385 [[maybe_unused]]
const unsigned int &col)
const
387 AssertThrow(
false, FeatureNotImplemented(
"el()"));
391template <
int dim,
int degree,
typename number>
396 inverse_diagonal_entries.reset();
397 global_to_local_solution.clear();
400template <
int dim,
int degree,
typename number>
404 unsigned int dof_handler_index)
const
406 data->initialize_dof_vector(dst, dof_handler_index);
409template <
int dim,
int degree,
typename number>
412 VectorType &dst)
const
414 for (
unsigned int j = 0; j < dealii::MatrixFreeOperators::BlockHelper::n_blocks(dst);
417 const std::vector<unsigned int> &constrained_dofs =
418 data->get_constrained_dofs(selected_fields[j]);
419 for (
const auto constrained_dof : constrained_dofs)
421 dealii::MatrixFreeOperators::BlockHelper::subblock(dst, j).local_element(
422 constrained_dof) = 1.0;
424 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
426 dealii::MatrixFreeOperators::BlockHelper::subblock(dst, j).local_element(
427 edge_constrained_indices[j][i]) = 1.0;
432template <
int dim,
int degree,
typename number>
435 const std::unordered_map<std::pair<unsigned int, dependencyType>,
437 pairHash> &_global_to_local_solution)
439 global_to_local_solution = _global_to_local_solution;
442template <
int dim,
int degree,
typename number>
443std::shared_ptr<const dealii::MatrixFree<dim, number>>
449template <
int dim,
int degree,
typename number>
450const std::shared_ptr<
451 dealii::DiagonalMatrix<typename matrixFreeOperator<dim, degree, number>::VectorType>> &
454 Assert(inverse_diagonal_entries.get() !=
nullptr && inverse_diagonal_entries->m() > 0,
455 dealii::ExcNotInitialized());
456 return inverse_diagonal_entries;
459template <
int dim,
int degree,
typename number>
462 const std::vector<VectorType *> &_src_solution_subset)
464 src_solution_subset = _src_solution_subset;
467template <
int dim,
int degree,
typename number>
470 std::vector<VectorType *> &dst,
471 const std::vector<VectorType *> &src)
const
473 Assert(!global_to_local_solution.empty(),
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"));
480 this->data->cell_loop(&matrixFreeOperator::compute_local_explicit_update,
487template <
int dim,
int degree,
typename number>
490 std::vector<VectorType *> &dst,
491 const std::vector<VectorType *> &src)
const
493 Assert(!global_to_local_solution.empty(),
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"));
500 this->data->cell_loop(&matrixFreeOperator::compute_local_postprocess_explicit_update,
507template <
int dim,
int degree,
typename number>
510 std::vector<VectorType *> &dst,
511 const std::vector<VectorType *> &src)
const
513 Assert(!global_to_local_solution.empty(),
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"));
520 this->data->cell_loop(&matrixFreeOperator::compute_local_nonexplicit_auxiliary_update,
527template <
int dim,
int degree,
typename number>
530 const VectorType &src)
const
532 Assert(!global_to_local_solution.empty(),
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"));
543 this->data->cell_loop(&matrixFreeOperator::compute_local_residual,
550template <
int dim,
int degree,
typename number>
553 const VectorType &src)
const
555 Assert(!global_to_local_solution.empty(),
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"));
564 this->data->cell_loop(&matrixFreeOperator::compute_local_newton_update,
571template <
int dim,
int degree,
typename number>
574 const VectorType &src)
const
576 this->vmult(dst, src);
579template <
int dim,
int degree,
typename number>
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
590 global_to_local_solution,
591 solveType::EXPLICIT_RHS);
594 variable_list.eval_local_operator(
596 const dealii::Point<dim, size_type> &q_point_loc)
598 this->compute_explicit_RHS(var_list, q_point_loc);
605template <
int dim,
int degree,
typename number>
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
616 global_to_local_solution,
617 solveType::POSTPROCESS);
620 variable_list.eval_local_operator(
622 const dealii::Point<dim, size_type> &q_point_loc)
624 this->compute_postprocess_explicit_RHS(var_list, q_point_loc);
631template <
int dim,
int degree,
typename number>
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
642 global_to_local_solution,
643 solveType::NONEXPLICIT_RHS);
646 variable_list.eval_local_operator(
648 const dealii::Point<dim, size_type> &q_point_loc)
650 this->compute_nonexplicit_RHS(var_list, q_point_loc);
657template <
int dim,
int degree,
typename number>
660 const dealii::MatrixFree<dim, number> &data,
662 [[maybe_unused]]
const VectorType &src,
663 const std::pair<unsigned int, unsigned int> &cell_range)
const
668 global_to_local_solution,
669 solveType::NONEXPLICIT_RHS);
672 variable_list.eval_local_operator(
674 const dealii::Point<dim, size_type> &q_point_loc)
676 this->compute_nonexplicit_RHS(var_list, q_point_loc);
683template <
int dim,
int degree,
typename number>
686 const dealii::MatrixFree<dim, number> &data,
688 const VectorType &src,
689 const std::pair<unsigned int, unsigned int> &cell_range)
const
694 global_to_local_solution,
695 solveType::NONEXPLICIT_LHS);
699 variable_list.eval_local_operator(
701 const dealii::Point<dim, size_type> &q_point_loc)
703 this->compute_nonexplicit_LHS(var_list, q_point_loc);
711template <
int dim,
int degree,
typename number>
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,
724 set_constrained_entries_to_one(inverse_diagonal);
726 for (
unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i)
728 Assert(inverse_diagonal.local_element(i) > 0.0,
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);
735template <
int dim,
int degree,
typename number>
738 const dealii::MatrixFree<dim, number> &data,
740 [[maybe_unused]]
const unsigned int &dummy,
741 const std::pair<unsigned int, unsigned int> &cell_range)
const
746 global_to_local_solution,
747 solveType::NONEXPLICIT_LHS);
750 variable_list.eval_local_diagonal(
752 const dealii::Point<dim, size_type> &q_point_loc)
754 this->compute_nonexplicit_LHS(var_list, q_point_loc);
761PRISMS_PF_END_NAMESPACE
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
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