Skip to content

Commit

Permalink
matrix-free: simplify diagonal computation
Browse files Browse the repository at this point in the history 
This uses compute_diagonal() from deal.II instead of hand-writting the
logic. For this, I need to break out the inner-most cell operation from
local_apply.
While we are here, clean-up the evaluation and remove some unnecessary
assignments.
  • Loading branch information
@tjhei
tjhei committed Jan 28, 2025
1 parent aa2bf45 commit 9bb7c56
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 100 deletions.
12 changes: 6 additions & 6 deletions include/aspect/stokes_matrix_free.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -299,14 +299,14 @@ namespace aspect
const dealii::LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;


/**
* Computes the diagonal contribution from a cell matrix.
* This function contains the inner-most operation done on a single cell
*/
void local_compute_diagonal (const MatrixFree<dim,number> &data,
dealii::LinearAlgebra::distributed::Vector<number> &dst,
const unsigned int &dummy,
const std::pair<unsigned int,unsigned int> &cell_range) const;
void inner_cell_operation(FEEvaluation<dim,
degree_p,
degree_p+2,
1,
number> &pressure) const;

/**
* A pointer to the current cell data that contains viscosity and other required parameters per cell.
Expand Down
144 changes: 50 additions & 94 deletions source/simulator/stokes_matrix_free.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -622,42 +622,53 @@ namespace aspect
{
FEEvaluation<dim,degree_p,degree_p+2,1,number> pressure (data, /*dofh*/1);

const bool use_viscosity_at_quadrature_points
= (cell_data->viscosity.size(1) == pressure.n_q_points);

for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
{
VectorizedArray<number> one_over_viscosity = cell_data->viscosity(cell, 0);
pressure.reinit (cell);
pressure.gather_evaluate (src, EvaluationFlags::values);
this->inner_cell_operation(pressure);
pressure.integrate_scatter (EvaluationFlags::values, dst);
}

const unsigned int n_components_filled = this->get_matrix_free()->n_active_entries_per_cell_batch(cell);
}

// The /= operator for VectorizedArray results in a floating point operation
// (divide by 0) since the (*viscosity)(cell) array is not completely filled.
// Therefore, we need to divide each entry manually.
for (unsigned int c=0; c<n_components_filled; ++c)
one_over_viscosity[c] = cell_data->pressure_scaling*cell_data->pressure_scaling/one_over_viscosity[c];

pressure.reinit (cell);
pressure.gather_evaluate (src, EvaluationFlags::values);
template <int dim, int degree_p, typename number>
void
MatrixFreeStokesOperators::MassMatrixOperator<dim,degree_p,number>
::inner_cell_operation(FEEvaluation<dim,
degree_p,
degree_p+2,
1,
number> &pressure) const
{
const bool use_viscosity_at_quadrature_points
= (cell_data->viscosity.size(1) == pressure.n_q_points);

for (const unsigned int q : pressure.quadrature_point_indices())
{
// Only update the viscosity if a Q1 projection is used.
if (use_viscosity_at_quadrature_points)
{
one_over_viscosity = cell_data->viscosity(cell, q);
const unsigned int cell = pressure.get_current_cell_index();
const unsigned int n_components_filled = this->get_matrix_free()->n_active_entries_per_cell_batch(cell);

const unsigned int n_components_filled = this->get_matrix_free()->n_active_entries_per_cell_batch(cell);
VectorizedArray<number> prefactor;

for (unsigned int c=0; c<n_components_filled; ++c)
one_over_viscosity[c] = cell_data->pressure_scaling*cell_data->pressure_scaling/one_over_viscosity[c];
}
// The /= operator for VectorizedArray results in a floating point operation
// (divide by 0) since the (*viscosity)(cell) array is not completely filled.
// Therefore, we need to divide each entry manually.
if (!use_viscosity_at_quadrature_points)
{
for (unsigned int c=0; c<n_components_filled; ++c)
prefactor[c] = cell_data->pressure_scaling*cell_data->pressure_scaling / cell_data->viscosity(cell, 0)[c];
}

pressure.submit_value(one_over_viscosity*
pressure.get_value(q),q);
for (const unsigned int q : pressure.quadrature_point_indices())
{
// Only update the viscosity if a Q1 projection is used.
if (use_viscosity_at_quadrature_points)
{
for (unsigned int c=0; c<n_components_filled; ++c)
prefactor[c] = cell_data->pressure_scaling*cell_data->pressure_scaling / cell_data->viscosity(cell, q)[c];
}

pressure.integrate_scatter (EvaluationFlags::values, dst);
pressure.submit_value(prefactor*pressure.get_value(q), q);
}
}

Expand Down Expand Up @@ -690,12 +701,23 @@ namespace aspect
dealii::LinearAlgebra::distributed::Vector<number> &diagonal =
this->diagonal_entries->get_vector();

unsigned int dummy = 0;
this->data->initialize_dof_vector(inverse_diagonal, /*dofh*/1);
this->data->initialize_dof_vector(diagonal, /*dofh*/1);

this->data->cell_loop (&MassMatrixOperator::local_compute_diagonal, this,
diagonal, dummy);
MatrixFreeTools::compute_diagonal<dim,degree_p,degree_p+2,1,number,VectorizedArray<number>,dealii::LinearAlgebra::distributed::Vector<number>>(
*(this->get_matrix_free()),
diagonal,
[&](FEEvaluation<dim,
degree_p,
degree_p+2,
1,
number> &pressure)
{
pressure.evaluate(EvaluationFlags::values);
this->inner_cell_operation(pressure);
pressure.integrate(EvaluationFlags::values);
},
1 /* dofhandler */);

this->set_constrained_entries_to_one(diagonal);
inverse_diagonal = diagonal;
Expand All @@ -715,72 +737,6 @@ namespace aspect



template <int dim, int degree_p, typename number>
void
MatrixFreeStokesOperators::MassMatrixOperator<dim,degree_p,number>
::local_compute_diagonal (const MatrixFree<dim,number> &data,
dealii::LinearAlgebra::distributed::Vector<number> &dst,
const unsigned int &,
const std::pair<unsigned int,unsigned int> &cell_range) const
{
FEEvaluation<dim,degree_p,degree_p+2,1,number> pressure (data, 1);

AlignedVector<VectorizedArray<number>> diagonal(pressure.dofs_per_cell);

const bool use_viscosity_at_quadrature_points
= (cell_data->viscosity.size(1) == pressure.n_q_points);

for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
{
VectorizedArray<number> one_over_viscosity = cell_data->viscosity(cell, 0);

const unsigned int n_components_filled = this->get_matrix_free()->n_active_entries_per_cell_batch(cell);

// The /= operator for VectorizedArray results in a floating point operation
// (divide by 0) since the (*viscosity)(cell) array is not completely filled.
// Therefore, we need to divide each entry manually.
for (unsigned int c=0; c<n_components_filled; ++c)
one_over_viscosity[c] = cell_data->pressure_scaling*cell_data->pressure_scaling/one_over_viscosity[c];

pressure.reinit (cell);
for (unsigned int i=0; i<pressure.dofs_per_cell; ++i)
{
for (unsigned int j=0; j<pressure.dofs_per_cell; ++j)
pressure.begin_dof_values()[j] = VectorizedArray<number>();
pressure.begin_dof_values()[i] = make_vectorized_array<number> (1.);

pressure.evaluate (EvaluationFlags::values);

for (const unsigned int q : pressure.quadrature_point_indices())
{
// Only update the viscosity if a Q1 projection is used.
if (use_viscosity_at_quadrature_points)
{
one_over_viscosity = cell_data->viscosity(cell, q);

const unsigned int n_components_filled = this->get_matrix_free()->n_active_entries_per_cell_batch(cell);

for (unsigned int c=0; c<n_components_filled; ++c)
one_over_viscosity[c] = cell_data->pressure_scaling*cell_data->pressure_scaling/one_over_viscosity[c];
}

pressure.submit_value(one_over_viscosity*
pressure.get_value(q),q);
}

pressure.integrate (EvaluationFlags::values);

diagonal[i] = pressure.begin_dof_values()[i];
}

for (unsigned int i=0; i<pressure.dofs_per_cell; ++i)
pressure.begin_dof_values()[i] = diagonal[i];
pressure.distribute_local_to_global (dst);
}
}



/**
* Velocity block operator
*/
Expand Down

0 comments on commit 9bb7c56

Please sign in to comment.