From 46a6f65ee3792fa991506bde7c4828908846d956 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 23 Aug 2024 16:35:17 -0400 Subject: [PATCH] Update step-64 --- examples/step-64/step-64.cc | 121 +++++++++++++++++++++++++++++++++++- 1 file changed, 119 insertions(+), 2 deletions(-) diff --git a/examples/step-64/step-64.cc b/examples/step-64/step-64.cc index 5e3624681a..3962daf896 100644 --- a/examples/step-64/step-64.cc +++ b/examples/step-64/step-64.cc @@ -136,9 +136,22 @@ namespace Step64 Portable::FEEvaluation *fe_eval, const int q_point) const; + DEAL_II_HOST_DEVICE void + set_matrix_free_data(const Portable::MatrixFree::Data &data) + { + gpu_data = &data; + } + + DEAL_II_HOST_DEVICE void set_cell(int new_cell) + { + cell = new_cell; + } + static const unsigned int n_q_points = dealii::Utilities::pow(fe_degree + 1, dim); + static const unsigned int n_local_dofs = n_q_points; + private: const typename Portable::MatrixFree::Data *gpu_data; double *coef; @@ -233,7 +246,7 @@ namespace Step64 // needs to have a `vmult()` function that performs the action of // the linear operator on a source vector. template - class HelmholtzOperator + class HelmholtzOperator : public Subscriptor { public: HelmholtzOperator(const DoFHandler &dof_handler, @@ -248,9 +261,25 @@ namespace Step64 LinearAlgebra::distributed::Vector &vec) const; + void compute_diagonal(); + + std::shared_ptr>> + get_matrix_diagonal_inverse() const; + + types::global_dof_index m() const; + + types::global_dof_index n() const; + + double el(const types::global_dof_index row, + const types::global_dof_index col) const; + private: Portable::MatrixFree mf_data; LinearAlgebra::distributed::Vector coef; + std::shared_ptr>> + inverse_diagonal_entries; }; @@ -324,6 +353,81 @@ namespace Step64 } + + template + void HelmholtzOperator::compute_diagonal() + { + this->inverse_diagonal_entries.reset( + new DiagonalMatrix< + LinearAlgebra::distributed::Vector>()); + LinearAlgebra::distributed::Vector + &inverse_diagonal = inverse_diagonal_entries->get_vector(); + initialize_dof_vector(inverse_diagonal); + + HelmholtzOperatorQuad helmholtz_operator_quad( + nullptr, coef.get_values(), -1); + + MatrixFreeTools::compute_diagonal( + mf_data, + inverse_diagonal, + helmholtz_operator_quad, + EvaluationFlags::values | EvaluationFlags::gradients, + EvaluationFlags::values | EvaluationFlags::gradients); + + double *raw_diagonal = inverse_diagonal.get_values(); + + Kokkos::parallel_for( + inverse_diagonal.locally_owned_size(), KOKKOS_LAMBDA(int i) { + Assert(raw_diagonal[i] > 0., + ExcMessage("No diagonal entry in a positive definite operator " + "should be zero")); + raw_diagonal[i] = 1. / raw_diagonal[i]; + }); + } + + + + template + std::shared_ptr>> + HelmholtzOperator::get_matrix_diagonal_inverse() const + { + return inverse_diagonal_entries; + } + + + + template + types::global_dof_index HelmholtzOperator::m() const + { + return mf_data.get_vector_partitioner()->size(); + } + + + + template + types::global_dof_index HelmholtzOperator::n() const + { + return mf_data.get_vector_partitioner()->size(); + } + + + + template + double + HelmholtzOperator::el(const types::global_dof_index row, + const types::global_dof_index col) const + { + (void)col; + Assert(row == col, ExcNotImplemented()); + Assert(inverse_diagonal_entries.get() != nullptr && + inverse_diagonal_entries->m() > 0, + ExcNotInitialized()); + return 1.0 / (*inverse_diagonal_entries)(row, row); + } + + + // @sect3{Class HelmholtzProblem} // This is the main class of this program. It defines the usual @@ -509,7 +613,20 @@ namespace Step64 template void HelmholtzProblem::solve() { - PreconditionIdentity preconditioner; + system_matrix_dev->compute_diagonal(); + + using PreconditionerType = PreconditionChebyshev< + HelmholtzOperator, + LinearAlgebra::distributed::Vector>; + typename PreconditionerType::AdditionalData additional_data; + additional_data.smoothing_range = 15.; + additional_data.degree = 5; + additional_data.eig_cg_n_iterations = 10; + additional_data.preconditioner = + system_matrix_dev->get_matrix_diagonal_inverse(); + + PreconditionerType preconditioner; + preconditioner.initialize(*system_matrix_dev, additional_data); SolverControl solver_control(system_rhs_dev.size(), 1e-12 * system_rhs_dev.l2_norm()); -- 2.39.5