From 9bb7c56a2c7dca0fa955d234da910f7f0edfec78 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Tue, 28 Jan 2025 15:16:44 -0500 Subject: [PATCH] matrix-free: simplify diagonal computation 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. --- include/aspect/stokes_matrix_free.h | 12 +-- source/simulator/stokes_matrix_free.cc | 144 +++++++++---------------- 2 files changed, 56 insertions(+), 100 deletions(-) diff --git a/include/aspect/stokes_matrix_free.h b/include/aspect/stokes_matrix_free.h index 5721636ce0d..e005fdec5ae 100644 --- a/include/aspect/stokes_matrix_free.h +++ b/include/aspect/stokes_matrix_free.h @@ -299,14 +299,14 @@ namespace aspect const dealii::LinearAlgebra::distributed::Vector &src, const std::pair &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 &data, - dealii::LinearAlgebra::distributed::Vector &dst, - const unsigned int &dummy, - const std::pair &cell_range) const; + void inner_cell_operation(FEEvaluation &pressure) const; /** * A pointer to the current cell data that contains viscosity and other required parameters per cell. diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc index cac32af4725..dd8d3023dcd 100644 --- a/source/simulator/stokes_matrix_free.cc +++ b/source/simulator/stokes_matrix_free.cc @@ -622,42 +622,53 @@ namespace aspect { FEEvaluation 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 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; cpressure_scaling*cell_data->pressure_scaling/one_over_viscosity[c]; - pressure.reinit (cell); - pressure.gather_evaluate (src, EvaluationFlags::values); + template + void + MatrixFreeStokesOperators::MassMatrixOperator + ::inner_cell_operation(FEEvaluation &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 prefactor; - for (unsigned int c=0; cpressure_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; cpressure_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; cpressure_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); } } @@ -690,12 +701,23 @@ namespace aspect dealii::LinearAlgebra::distributed::Vector &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,dealii::LinearAlgebra::distributed::Vector>( + *(this->get_matrix_free()), + diagonal, + [&](FEEvaluation &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; @@ -715,72 +737,6 @@ namespace aspect - template - void - MatrixFreeStokesOperators::MassMatrixOperator - ::local_compute_diagonal (const MatrixFree &data, - dealii::LinearAlgebra::distributed::Vector &dst, - const unsigned int &, - const std::pair &cell_range) const - { - FEEvaluation pressure (data, 1); - - AlignedVector> 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 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; cpressure_scaling*cell_data->pressure_scaling/one_over_viscosity[c]; - - pressure.reinit (cell); - for (unsigned int i=0; i(); - pressure.begin_dof_values()[i] = make_vectorized_array (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; cpressure_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