Skip to content

Commit

Permalink
Adding std::function to variableContainer
Browse files Browse the repository at this point in the history
  • Loading branch information
@landinjm
landinjm committed Jan 14, 2025
1 parent 7dfef89 commit 90eaee1
Show file tree
Hide file tree
Showing 3 changed files with 115 additions and 199 deletions.
71 changes: 18 additions & 53 deletions include/core/matrix_free_operator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,30 +105,16 @@ matrixFreeOperator<dim, degree, number>::local_apply(
// Constructor for FEEvaluation objects
variableContainer<dim, degree, number> variable_list(data);

for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
// 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)
{
// Initialize, read DOFs, and set evaulation flags for each variable
variable_list.reinit_and_eval(src, cell);

// Number of quadrature points
unsigned int n_q_points = variable_list.get_num_q_points();

for (unsigned int q = 0; q < n_q_points; ++q)
{
// Set the quadrature point
variable_list.q_point = q;

// Grab the quadrature point location
Point<dim, VectorizedArray<number>> q_point_loc =
variable_list.get_q_point_location();

// Calculate the residuals
nonexplicit_RHS(variable_list, q_point_loc);
}

// Integrate and add to global vector dst
variable_list.integrate_and_distribute(dst);
}
this->nonexplicit_RHS(var_list, q_point_loc);
},
dst,
src,
cell_range);
}

template <int dim, int degree, typename number>
Expand Down Expand Up @@ -193,37 +179,16 @@ matrixFreeOperator<dim, degree, number>::local_compute_diagonal(
// Constructor for FEEvaluation objects
variableContainer<dim, degree, number> variable_list(data);

variable_list.init_diagonal(field_index);

for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
// 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)
{
variable_list.reinit_cell(cell);

for (unsigned int i = 0; i < variable_list.n_dofs_per_cell; ++i)
{
variable_list.reinit_and_eval_diagonal(field_index, i);

// Number of quadrature points
unsigned int n_q_points = variable_list.get_num_q_points();

for (unsigned int q = 0; q < n_q_points; ++q)
{
// Set the quadrature point
variable_list.q_point = q;

// Grab the quadrature point location
Point<dim, VectorizedArray<number>> q_point_loc =
variable_list.get_q_point_location();

// Calculate the residuals
nonexplicit_RHS(variable_list, q_point_loc);
}

variable_list.integrate_and_assemble_local_diagonal(field_index, i);
}

variable_list.distribute_diagonal(dst, field_index);
}
this->nonexplicit_RHS(var_list, q_point_loc);
},
dst,
cell_range,
field_index);
}

#endif
122 changes: 97 additions & 25 deletions include/core/variableContainer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

#include <core/exceptions.h>
#include <core/model_variables.h>
#include <core/variableAttributes.h>

Expand Down Expand Up @@ -71,37 +72,108 @@ class variableContainer
integrate_and_distribute(dealii::LinearAlgebra::distributed::Vector<number> &dst);

/**
* \brief Initialize the `AlignedVector` for the diagonal based on the DoFs per cell.
* Additionally, assign `n_dofs_per_cell` based on the current field index.
* \brief TODO: Add comments
*/
void
init_diagonal(unsigned int global_var_index);
eval_local_operator(
const std::function<void(variableContainer &,
const dealii::Point<dim, dealii::VectorizedArray<number>> &)>
&func,
dealii::LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range)
{
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
{
// Initialize, read DOFs, and set evaulation flags for each variable
reinit_and_eval(src, cell);

// Number of quadrature points
unsigned int n_q_points = get_num_q_points();

for (unsigned int q = 0; q < n_q_points; ++q)
{
// Set the quadrature point
q_point = q;

// Grab the quadrature point location
Point<dim, VectorizedArray<number>> q_point_loc = get_q_point_location();

// Calculate the residuals
func(*this, q_point_loc);
}

// Integrate and add to global vector dst
integrate_and_distribute(dst);
}
};

/**
* \brief Initialize the operation ptr to the cell for all field indices.
* \brief TODO: Add comments
*/
void
reinit_cell(unsigned int cell);

/**
* \brief Submit the DoFs for the diagonal identity matrix and evaluate based on the
* evaluation flags.
*/
void
reinit_and_eval_diagonal(unsigned int global_var_index, unsigned int i);

/**
* \brief Integrate the residuals and assemble the local diagonal
*/
void
integrate_and_assemble_local_diagonal(unsigned int global_var_index, unsigned int i);

/**
* \brief Distribute the diagonal from local to global.
*/
void
distribute_diagonal(dealii::LinearAlgebra::distributed::Vector<number> &dst,
unsigned int global_var_index);
eval_local_diagonal(
const std::function<void(variableContainer &,
const dealii::Point<dim, dealii::VectorizedArray<number>> &)>
&func,
dealii::LinearAlgebra::distributed::Vector<number> &dst,
const std::pair<unsigned int, unsigned int> &cell_range,
const unsigned int &global_var_index)
{
const auto &var_info = varInfoList[global_var_index];

AssertThrow(var_info.var_type != fieldType::VECTOR,
FeatureNotImplemented("Vector multigrid"));

auto *scalar_FEEval_ptr = scalar_vars_map.at(global_var_index).get();

n_dofs_per_cell = scalar_FEEval_ptr->dofs_per_cell;
diagonal = std::make_unique<dealii::AlignedVector<dealii::VectorizedArray<number>>>(
n_dofs_per_cell);

for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
{
scalar_FEEval_ptr->reinit(cell);

for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < n_dofs_per_cell; ++j)
{
scalar_FEEval_ptr->submit_dof_value(dealii::VectorizedArray<number>(), j);
}
scalar_FEEval_ptr->submit_dof_value(dealii::make_vectorized_array<number>(
1.0),
i);

scalar_FEEval_ptr->evaluate(var_info.evaluation_flags);

// Number of quadrature points
unsigned int n_q_points = get_num_q_points();

for (unsigned int q = 0; q < n_q_points; ++q)
{
// Set the quadrature point
q_point = q;

// Grab the quadrature point location
dealii::Point<dim, dealii::VectorizedArray<number>> q_point_loc =
get_q_point_location();

// Calculate the residuals
func(*this, q_point_loc);
}

scalar_FEEval_ptr->integrate(var_info.residual_flags);
(*diagonal)[i] = scalar_FEEval_ptr->get_dof_value(i);
}

for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
{
scalar_FEEval_ptr->submit_dof_value((*diagonal)[i], i);
}
scalar_FEEval_ptr->distribute_local_to_global(dst);
}
};

// The quadrature point index, a method to get the number of quadrature points
// per cell, and a method to get the xyz coordinates for the quadrature point
Expand Down
121 changes: 0 additions & 121 deletions src/core/variableContainer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,127 +143,6 @@ variableContainer<dim, degree, number>::integrate_and_distribute(
}
}

template <int dim, int degree, typename number>
void
variableContainer<dim, degree, number>::init_diagonal(unsigned int global_var_index)
{
const auto &var_info = varInfoList[global_var_index];

if (var_info.var_type == fieldType::SCALAR)
{
n_dofs_per_cell = scalar_vars_map.at(global_var_index)->dofs_per_cell;
}
else if (var_info.var_type == fieldType::VECTOR)
{
n_dofs_per_cell = vector_vars_map.at(global_var_index)->dofs_per_cell;
}

diagonal = std::make_unique<dealii::AlignedVector<dealii::VectorizedArray<number>>>(
n_dofs_per_cell);
}

template <int dim, int degree, typename number>
void
variableContainer<dim, degree, number>::reinit_cell(unsigned int cell)
{
for (unsigned int i = 0; i < num_var; i++)
{
const auto &var_info = varInfoList[i];

const unsigned int var_index = var_info.global_var_index;

if (var_info.var_type == fieldType::SCALAR)
{
scalar_vars_map.at(var_index).get()->reinit(cell);
}
else if (var_info.var_type == fieldType::VECTOR)
{
vector_vars_map.at(var_index).get()->reinit(cell);
}
}
}

template <int dim, int degree, typename number>
void
variableContainer<dim, degree, number>::reinit_and_eval_diagonal(
unsigned int global_var_index,
unsigned int i)
{
const auto &var_info = varInfoList[global_var_index];

if (var_info.var_type == fieldType::SCALAR)
{
auto *scalar_FEEval_ptr = scalar_vars_map.at(global_var_index).get();
for (unsigned int j = 0; j < n_dofs_per_cell; ++j)
{
scalar_FEEval_ptr->submit_dof_value(dealii::VectorizedArray<number>(), j);
}
scalar_FEEval_ptr->submit_dof_value(dealii::make_vectorized_array<number>(1.0), i);
}
else if (var_info.var_type == fieldType::VECTOR)
{
AssertThrow(false, FeatureNotImplemented("Vector multigrid"));
}

for (unsigned int i = 0; i < num_var; i++)
{
const auto &var_info = varInfoList[i];
const unsigned int var_index = var_info.global_var_index;

if (var_info.var_type == fieldType::SCALAR)
{
scalar_vars_map.at(var_index)->evaluate(var_info.evaluation_flags);
}
else if (var_info.var_type == fieldType::VECTOR)
{
vector_vars_map.at(var_index)->evaluate(var_info.evaluation_flags);
}
}
}

template <int dim, int degree, typename number>
void
variableContainer<dim, degree, number>::integrate_and_assemble_local_diagonal(
unsigned int global_var_index,
unsigned int i)
{
const auto &var_info = varInfoList[global_var_index];

if (var_info.var_type == fieldType::SCALAR)
{
auto *scalar_FEEval_ptr = scalar_vars_map.at(global_var_index).get();
scalar_FEEval_ptr->integrate(var_info.residual_flags);
(*diagonal)[i] = scalar_FEEval_ptr->get_dof_value(i);
}
else if (var_info.var_type == fieldType::VECTOR)
{
AssertThrow(false, FeatureNotImplemented("Vector multigrid"));
}
}

template <int dim, int degree, typename number>
void
variableContainer<dim, degree, number>::distribute_diagonal(
dealii::LinearAlgebra::distributed::Vector<number> &dst,
unsigned int global_var_index)
{
const auto &var_info = varInfoList[global_var_index];

if (var_info.var_type == fieldType::SCALAR)
{
auto *scalar_FEEval_ptr = scalar_vars_map.at(global_var_index).get();
for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
{
scalar_FEEval_ptr->submit_dof_value((*diagonal)[i], i);
}
scalar_FEEval_ptr->distribute_local_to_global(dst);
}
else if (var_info.var_type == fieldType::VECTOR)
{
AssertThrow(false, FeatureNotImplemented("Vector multigrid"));
}
}

template <int dim, int degree, typename number>
dealii::VectorizedArray<number>
variableContainer<dim, degree, number>::get_scalar_value(
Expand Down

0 comments on commit 90eaee1

Please sign in to comment.