Skip to content

Commit

Permalink
user facing functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@landinjm
landinjm committed Jan 16, 2025
1 parent 47919cd commit 260cc37
Show file tree
Hide file tree
Showing 4 changed files with 210 additions and 191 deletions.
32 changes: 30 additions & 2 deletions applications/multigrid_test/custom_PDE.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,36 @@ using namespace dealii;
template <int dim, int degree, typename number>
class customPDE : public matrixFreeOperator<dim, degree, number>
{
public:
using scalarValue = VectorizedArray<number>;
using scalarGrad = Tensor<1, dim, VectorizedArray<number>>;
using scalarHess = Tensor<2, dim, VectorizedArray<number>>;
using vectorValue = Tensor<1, dim, VectorizedArray<number>>;
using vectorGrad = Tensor<2, dim, VectorizedArray<number>>;
using vectorHess = Tensor<3, dim, VectorizedArray<number>>;

/**
* \brief Constructor.
*/
customPDE();
};
customPDE() = default;

/**
* \brief User-implemented class for the RHS of nonexplicit equations.
*/
void
nonexplicit_RHS(variableContainer<dim, degree, number> &variable_list,
const Point<dim, VectorizedArray<number>> &q_point_loc) const override;
};

template <int dim, int degree, typename number>
void
customPDE<dim, degree, number>::nonexplicit_RHS(
variableContainer<dim, degree, number> &variable_list,
const Point<dim, VectorizedArray<number>> &q_point_loc) const
{
scalarGrad phi = variable_list.get_scalar_gradient(0);

scalarValue coefficient = 1.0 / (0.05 + 2.0 * q_point_loc.square());

variable_list.set_scalar_gradient_term(0, phi * coefficient);
}
28 changes: 15 additions & 13 deletions applications/multigrid_test/main.cc
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
#include "custom_PDE.h"

#include <core/temp_test.h>

int
Expand All @@ -12,23 +14,23 @@ main(int argc, char *argv[])
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << '\n'
<< '\n'
<< "----------------------------------------------------" << '\n';
std::cerr << "Exception on processing: " << '\n'
<< exc.what() << '\n'
<< "Aborting!" << '\n'
<< "----------------------------------------------------" << '\n';
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << '\n'
<< '\n'
<< "----------------------------------------------------" << '\n';
std::cerr << "Unknown exception!" << '\n'
<< "Aborting!" << '\n'
<< "----------------------------------------------------" << '\n';
return 1;
}

Expand Down
90 changes: 36 additions & 54 deletions include/core/matrix_free_operator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,6 @@

#include <core/variableContainer.h>

using namespace dealii;

/**
* \brief This is the abstract base class for the matrix free implementation of some PDE.
*
Expand All @@ -19,16 +17,10 @@ using namespace dealii;
*/
template <int dim, int degree, typename number>
class matrixFreeOperator
: public MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>
: public dealii::MatrixFreeOperators::
Base<dim, dealii::LinearAlgebra::distributed::Vector<number>>
{
public:
using scalarValue = VectorizedArray<number>;
using scalarGrad = Tensor<1, dim, VectorizedArray<number>>;
using scalarHess = Tensor<2, dim, VectorizedArray<number>>;
using vectorValue = Tensor<1, dim, VectorizedArray<number>>;
using vectorGrad = Tensor<2, dim, VectorizedArray<number>>;
using vectorHess = Tensor<3, dim, VectorizedArray<number>>;

/**
* \brief Constructor.
*/
Expand All @@ -49,66 +41,69 @@ class matrixFreeOperator

protected:
/**
* \brief User-implemented PDE.
* \brief User-implemented class for the RHS of nonexplicit equations.
*/
virtual void
nonexplicit_RHS(variableContainer<dim, degree, number> &variable_list,
const Point<dim, VectorizedArray<number>> &q_point_loc) const;
nonexplicit_RHS(
variableContainer<dim, degree, number> &variable_list,
const dealii::Point<dim, dealii::VectorizedArray<number>> &q_point_loc) const = 0;

private:
/**
* \brief Apply operator to src and add result in dst.
*/
void
apply_add(LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src) const override;
apply_add(dealii::LinearAlgebra::distributed::Vector<number> &dst,
const dealii::LinearAlgebra::distributed::Vector<number> &src) const override;

/**
* \brief Local application of the operator.
*/
void
local_apply(const MatrixFree<dim, number> &data,
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;
local_apply(const dealii::MatrixFree<dim, number> &data,
dealii::LinearAlgebra::distributed::Vector<number> &dst,
const dealii::LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;

/**
* \brief Local computation of the diagonal of the operator.
*/
void
local_compute_diagonal(const MatrixFree<dim, number> &data,
LinearAlgebra::distributed::Vector<number> &dst,
const unsigned int &dummy,
local_compute_diagonal(const dealii::MatrixFree<dim, number> &data,
dealii::LinearAlgebra::distributed::Vector<number> &dst,
const unsigned int &dummy,
const std::pair<unsigned int, unsigned int> &cell_range) const;
};

template <int dim, int degree, typename number>
matrixFreeOperator<dim, degree, number>::matrixFreeOperator()
: MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>()
: dealii::MatrixFreeOperators::
Base<dim, dealii::LinearAlgebra::distributed::Vector<number>>()
{}

template <int dim, int degree, typename number>
void
matrixFreeOperator<dim, degree, number>::clear()
{
MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>::clear();
dealii::MatrixFreeOperators::
Base<dim, dealii::LinearAlgebra::distributed::Vector<number>>::clear();
}

template <int dim, int degree, typename number>
void
matrixFreeOperator<dim, degree, number>::local_apply(
const MatrixFree<dim, number> &data,
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const
const dealii::MatrixFree<dim, number> &data,
dealii::LinearAlgebra::distributed::Vector<number> &dst,
const dealii::LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const
{
// Constructor for FEEvaluation objects
variableContainer<dim, degree, number> variable_list(data);

// Initialize, evaluate, and submit based on user function.
variable_list.eval_local_operator(
[this](variableContainer<dim, degree, number> &var_list,
const Point<dim, VectorizedArray<number>> &q_point_loc)
[this](variableContainer<dim, degree, number> &var_list,
const dealii::Point<dim, dealii::VectorizedArray<number>> &q_point_loc)
{
this->nonexplicit_RHS(var_list, q_point_loc);
},
Expand All @@ -117,24 +112,11 @@ matrixFreeOperator<dim, degree, number>::local_apply(
cell_range);
}

template <int dim, int degree, typename number>
void
matrixFreeOperator<dim, degree, number>::nonexplicit_RHS(
variableContainer<dim, degree, number> &variable_list,
const Point<dim, VectorizedArray<number>> &q_point_loc) const
{
scalarGrad phi = variable_list.get_scalar_gradient(0);

scalarValue coefficient = 1.0 / (0.05 + 2.0 * q_point_loc.square());

variable_list.set_scalar_gradient_term(0, phi * coefficient);
}

template <int dim, int degree, typename number>
void
matrixFreeOperator<dim, degree, number>::apply_add(
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src) const
dealii::LinearAlgebra::distributed::Vector<number> &dst,
const dealii::LinearAlgebra::distributed::Vector<number> &src) const
{
this->data->cell_loop(&matrixFreeOperator::local_apply, this, dst, src);
}
Expand All @@ -144,8 +126,8 @@ void
matrixFreeOperator<dim, degree, number>::compute_diagonal()
{
this->inverse_diagonal_entries.reset(
new DiagonalMatrix<LinearAlgebra::distributed::Vector<number>>());
LinearAlgebra::distributed::Vector<number> &inverse_diagonal =
new DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<number>>());
dealii::LinearAlgebra::distributed::Vector<number> &inverse_diagonal =
this->inverse_diagonal_entries->get_vector();
this->data->initialize_dof_vector(inverse_diagonal);
unsigned int dummy = 0;
Expand All @@ -159,7 +141,7 @@ matrixFreeOperator<dim, degree, number>::compute_diagonal()
for (unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i)
{
Assert(inverse_diagonal.local_element(i) > 0.0,
ExcMessage(
dealii::ExcMessage(
"No diagonal entry in a positive definite operator should be zero"));
inverse_diagonal.local_element(i) = 1.0 / inverse_diagonal.local_element(i);
}
Expand All @@ -168,10 +150,10 @@ matrixFreeOperator<dim, degree, number>::compute_diagonal()
template <int dim, int degree, typename number>
void
matrixFreeOperator<dim, degree, number>::local_compute_diagonal(
const MatrixFree<dim, number> &data,
LinearAlgebra::distributed::Vector<number> &dst,
[[maybe_unused]] const unsigned int &dummy,
const std::pair<unsigned int, unsigned int> &cell_range) const
const dealii::MatrixFree<dim, number> &data,
dealii::LinearAlgebra::distributed::Vector<number> &dst,
[[maybe_unused]] const unsigned int &dummy,
const std::pair<unsigned int, unsigned int> &cell_range) const
{
// Field index
unsigned int field_index = 0;
Expand All @@ -181,8 +163,8 @@ matrixFreeOperator<dim, degree, number>::local_compute_diagonal(

// Initialize, evaluate, and submit diagonal based on user function.
variable_list.eval_local_diagonal(
[this](variableContainer<dim, degree, number> &var_list,
const Point<dim, VectorizedArray<number>> &q_point_loc)
[this](variableContainer<dim, degree, number> &var_list,
const dealii::Point<dim, dealii::VectorizedArray<number>> &q_point_loc)
{
this->nonexplicit_RHS(var_list, q_point_loc);
},
Expand Down
Loading

0 comments on commit 260cc37

Please sign in to comment.