From 260cc37c8da71d49fa2c6c95f382a318129b42bc Mon Sep 17 00:00:00 2001 From: landinjm Date: Thu, 16 Jan 2025 16:54:37 -0500 Subject: [PATCH] user facing functions --- applications/multigrid_test/custom_PDE.h | 32 ++- applications/multigrid_test/main.cc | 28 +-- include/core/matrix_free_operator.h | 90 ++++---- include/core/temp_test.h | 251 ++++++++++++----------- 4 files changed, 210 insertions(+), 191 deletions(-) diff --git a/applications/multigrid_test/custom_PDE.h b/applications/multigrid_test/custom_PDE.h index 81ae379a..92d758a7 100644 --- a/applications/multigrid_test/custom_PDE.h +++ b/applications/multigrid_test/custom_PDE.h @@ -13,8 +13,36 @@ using namespace dealii; template class customPDE : public matrixFreeOperator { +public: + using scalarValue = VectorizedArray; + using scalarGrad = Tensor<1, dim, VectorizedArray>; + using scalarHess = Tensor<2, dim, VectorizedArray>; + using vectorValue = Tensor<1, dim, VectorizedArray>; + using vectorGrad = Tensor<2, dim, VectorizedArray>; + using vectorHess = Tensor<3, dim, VectorizedArray>; + /** * \brief Constructor. */ - customPDE(); -}; \ No newline at end of file + customPDE() = default; + + /** + * \brief User-implemented class for the RHS of nonexplicit equations. + */ + void + nonexplicit_RHS(variableContainer &variable_list, + const Point> &q_point_loc) const override; +}; + +template +void +customPDE::nonexplicit_RHS( + variableContainer &variable_list, + const Point> &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); +} \ No newline at end of file diff --git a/applications/multigrid_test/main.cc b/applications/multigrid_test/main.cc index bc73f53a..c4bfe737 100644 --- a/applications/multigrid_test/main.cc +++ b/applications/multigrid_test/main.cc @@ -1,3 +1,5 @@ +#include "custom_PDE.h" + #include int @@ -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; } diff --git a/include/core/matrix_free_operator.h b/include/core/matrix_free_operator.h index 0f2e5dc7..aa05f29b 100644 --- a/include/core/matrix_free_operator.h +++ b/include/core/matrix_free_operator.h @@ -7,8 +7,6 @@ #include -using namespace dealii; - /** * \brief This is the abstract base class for the matrix free implementation of some PDE. * @@ -19,16 +17,10 @@ using namespace dealii; */ template class matrixFreeOperator - : public MatrixFreeOperators::Base> + : public dealii::MatrixFreeOperators:: + Base> { public: - using scalarValue = VectorizedArray; - using scalarGrad = Tensor<1, dim, VectorizedArray>; - using scalarHess = Tensor<2, dim, VectorizedArray>; - using vectorValue = Tensor<1, dim, VectorizedArray>; - using vectorGrad = Tensor<2, dim, VectorizedArray>; - using vectorHess = Tensor<3, dim, VectorizedArray>; - /** * \brief Constructor. */ @@ -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 &variable_list, - const Point> &q_point_loc) const; + nonexplicit_RHS( + variableContainer &variable_list, + const dealii::Point> &q_point_loc) const = 0; private: /** * \brief Apply operator to src and add result in dst. */ void - apply_add(LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &src) const override; + apply_add(dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const override; /** * \brief Local application of the operator. */ void - local_apply(const MatrixFree &data, - LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const; + local_apply(const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const; /** * \brief Local computation of the diagonal of the operator. */ void - local_compute_diagonal(const MatrixFree &data, - LinearAlgebra::distributed::Vector &dst, - const unsigned int &dummy, + local_compute_diagonal(const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const unsigned int &dummy, const std::pair &cell_range) const; }; template matrixFreeOperator::matrixFreeOperator() - : MatrixFreeOperators::Base>() + : dealii::MatrixFreeOperators:: + Base>() {} template void matrixFreeOperator::clear() { - MatrixFreeOperators::Base>::clear(); + dealii::MatrixFreeOperators:: + Base>::clear(); } template void matrixFreeOperator::local_apply( - const MatrixFree &data, - LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const + const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const { // Constructor for FEEvaluation objects variableContainer variable_list(data); // Initialize, evaluate, and submit based on user function. variable_list.eval_local_operator( - [this](variableContainer &var_list, - const Point> &q_point_loc) + [this](variableContainer &var_list, + const dealii::Point> &q_point_loc) { this->nonexplicit_RHS(var_list, q_point_loc); }, @@ -117,24 +112,11 @@ matrixFreeOperator::local_apply( cell_range); } -template -void -matrixFreeOperator::nonexplicit_RHS( - variableContainer &variable_list, - const Point> &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 void matrixFreeOperator::apply_add( - LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &src) const + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const { this->data->cell_loop(&matrixFreeOperator::local_apply, this, dst, src); } @@ -144,8 +126,8 @@ void matrixFreeOperator::compute_diagonal() { this->inverse_diagonal_entries.reset( - new DiagonalMatrix>()); - LinearAlgebra::distributed::Vector &inverse_diagonal = + new DiagonalMatrix>()); + dealii::LinearAlgebra::distributed::Vector &inverse_diagonal = this->inverse_diagonal_entries->get_vector(); this->data->initialize_dof_vector(inverse_diagonal); unsigned int dummy = 0; @@ -159,7 +141,7 @@ matrixFreeOperator::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); } @@ -168,10 +150,10 @@ matrixFreeOperator::compute_diagonal() template void matrixFreeOperator::local_compute_diagonal( - const MatrixFree &data, - LinearAlgebra::distributed::Vector &dst, - [[maybe_unused]] const unsigned int &dummy, - const std::pair &cell_range) const + const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + [[maybe_unused]] const unsigned int &dummy, + const std::pair &cell_range) const { // Field index unsigned int field_index = 0; @@ -181,8 +163,8 @@ matrixFreeOperator::local_compute_diagonal( // Initialize, evaluate, and submit diagonal based on user function. variable_list.eval_local_diagonal( - [this](variableContainer &var_list, - const Point> &q_point_loc) + [this](variableContainer &var_list, + const dealii::Point> &q_point_loc) { this->nonexplicit_RHS(var_list, q_point_loc); }, diff --git a/include/core/temp_test.h b/include/core/temp_test.h index d0d87949..d4221f54 100644 --- a/include/core/temp_test.h +++ b/include/core/temp_test.h @@ -32,6 +32,12 @@ #include #include +/** + * Forward declaration for user-implemented PDE class. + */ +template +class customPDE; + const unsigned int degree = 2; const unsigned int dimension = 3; @@ -39,6 +45,9 @@ template class LaplaceProblem { public: + using SystemMatrixType = customPDE; + using LevelMatrixType = customPDE; + LaplaceProblem(); void run(); @@ -53,23 +62,23 @@ class LaplaceProblem void output_results(unsigned int cycle) const; - parallel::distributed::Triangulation triangulation; + dealii::parallel::distributed::Triangulation triangulation; - const FE_Q fe; - DoFHandler dof_handler; + const dealii::FE_Q fe; + dealii::DoFHandler dof_handler; - MappingQ1 mapping; + dealii::MappingQ1 mapping; + + dealii::AffineConstraints constraints; - AffineConstraints constraints; - using SystemMatrixType = matrixFreeOperator; SystemMatrixType system_matrix; - MGConstrainedDoFs mg_constrained_dofs; - using LevelMatrixType = matrixFreeOperator; - MGLevelObject mg_matrices; + dealii::MGConstrainedDoFs mg_constrained_dofs; + + dealii::MGLevelObject mg_matrices; - LinearAlgebra::distributed::Vector solution; - LinearAlgebra::distributed::Vector system_rhs; + dealii::LinearAlgebra::distributed::Vector solution; + dealii::LinearAlgebra::distributed::Vector system_rhs; double setup_time {}; }; @@ -78,8 +87,8 @@ template LaplaceProblem::LaplaceProblem() : triangulation( MPI_COMM_WORLD, - Triangulation::limit_level_difference_at_vertices, - parallel::distributed::Triangulation::construct_multigrid_hierarchy) + dealii::Triangulation::limit_level_difference_at_vertices, + dealii::parallel::distributed::Triangulation::construct_multigrid_hierarchy) , fe(degree) , dof_handler(triangulation) {} @@ -90,96 +99,92 @@ LaplaceProblem::setup_system() { Timer time; setup_time = 0; - { - system_matrix.clear(); - mg_matrices.clear_elements(); - - dof_handler.distribute_dofs(fe); - dof_handler.distribute_mg_dofs(); - - conditionalOStreams::pout_base - << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; - - constraints.clear(); - constraints.reinit(dof_handler.locally_owned_dofs(), - DoFTools::extract_locally_relevant_dofs(dof_handler)); - DoFTools::make_hanging_node_constraints(dof_handler, constraints); - VectorTools::interpolate_boundary_values(mapping, - dof_handler, - 0, - Functions::ZeroFunction(), - constraints); - constraints.close(); - } + + system_matrix.clear(); + mg_matrices.clear_elements(); + + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + conditionalOStreams::pout_base + << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; + + constraints.clear(); + constraints.reinit(dof_handler.locally_owned_dofs(), + dealii::DoFTools::extract_locally_relevant_dofs(dof_handler)); + dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints); + dealii::VectorTools::interpolate_boundary_values(mapping, + dof_handler, + 0, + dealii::Functions::ZeroFunction(), + constraints); + constraints.close(); + setup_time += time.wall_time(); conditionalOStreams::pout_base << "Distribute DoFs & B.C. (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << 's' << '\n'; time.restart(); - { - { - typename MatrixFree::AdditionalData additional_data; - additional_data.tasks_parallel_scheme = - MatrixFree::AdditionalData::none; - additional_data.mapping_update_flags = - (update_gradients | update_JxW_values | update_quadrature_points); - std::shared_ptr> system_mf_storage( - new MatrixFree()); - system_mf_storage->reinit(mapping, - dof_handler, - constraints, - QGauss<1>(fe.degree + 1), - additional_data); - system_matrix.initialize(system_mf_storage); - } - system_matrix.initialize_dof_vector(solution); - system_matrix.initialize_dof_vector(system_rhs); - } + typename dealii::MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + dealii::MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + std::shared_ptr> system_mf_storage( + new dealii::MatrixFree()); + system_mf_storage->reinit(mapping, + dof_handler, + constraints, + dealii::QGauss<1>(fe.degree + 1), + additional_data); + system_matrix.initialize(system_mf_storage); + + system_matrix.initialize_dof_vector(solution); + system_matrix.initialize_dof_vector(system_rhs); + setup_time += time.wall_time(); conditionalOStreams::pout_base << "Setup matrix-free system (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << 's' << '\n'; time.restart(); - { - const unsigned int nlevels = triangulation.n_global_levels(); - mg_matrices.resize(0, nlevels - 1); - - const std::set dirichlet_boundary_ids = {0}; - mg_constrained_dofs.initialize(dof_handler); - mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, - dirichlet_boundary_ids); - - for (unsigned int level = 0; level < nlevels; ++level) - { - AffineConstraints level_constraints( - dof_handler.locally_owned_mg_dofs(level), - DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); - for (const types::global_dof_index dof_index : - mg_constrained_dofs.get_boundary_indices(level)) - { - level_constraints.constrain_dof_to_zero(dof_index); - } - level_constraints.close(); - - typename MatrixFree::AdditionalData additional_data; - additional_data.tasks_parallel_scheme = - MatrixFree::AdditionalData::none; - additional_data.mapping_update_flags = - (update_gradients | update_JxW_values | update_quadrature_points); - additional_data.mg_level = level; - std::shared_ptr> mg_mf_storage_level = - std::make_shared>(); - mg_mf_storage_level->reinit(mapping, - dof_handler, - level_constraints, - QGauss<1>(fe.degree + 1), - additional_data); - - mg_matrices[level].initialize(mg_mf_storage_level, mg_constrained_dofs, level); - } - } + const unsigned int nlevels = triangulation.n_global_levels(); + mg_matrices.resize(0, nlevels - 1); + + const std::set dirichlet_boundary_ids = {0}; + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, dirichlet_boundary_ids); + + for (unsigned int level = 0; level < nlevels; ++level) + { + dealii::AffineConstraints level_constraints( + dof_handler.locally_owned_mg_dofs(level), + dealii::DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); + for (const dealii::types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + { + level_constraints.constrain_dof_to_zero(dof_index); + } + level_constraints.close(); + + typename dealii::MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + dealii::MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + additional_data.mg_level = level; + std::shared_ptr> mg_mf_storage_level = + std::make_shared>(); + mg_mf_storage_level->reinit(mapping, + dof_handler, + level_constraints, + dealii::QGauss<1>(fe.degree + 1), + additional_data); + + mg_matrices[level].initialize(mg_mf_storage_level, mg_constrained_dofs, level); + } + setup_time += time.wall_time(); conditionalOStreams::pout_base << "Setup matrix-free levels (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << 's' @@ -193,19 +198,19 @@ LaplaceProblem::assemble_rhs() Timer time; system_rhs = 0; - FEEvaluation phi(*system_matrix.get_matrix_free()); + dealii::FEEvaluation phi(*system_matrix.get_matrix_free()); for (unsigned int cell = 0; cell < system_matrix.get_matrix_free()->n_cell_batches(); ++cell) { phi.reinit(cell); for (const unsigned int q : phi.quadrature_point_indices()) { - phi.submit_value(make_vectorized_array(1.0), q); + phi.submit_value(dealii::make_vectorized_array(1.0), q); } - phi.integrate(EvaluationFlags::values); + phi.integrate(dealii::EvaluationFlags::EvaluationFlags::values); phi.distribute_local_to_global(system_rhs); } - system_rhs.compress(VectorOperation::add); + system_rhs.compress(dealii::VectorOperation::add); setup_time += time.wall_time(); conditionalOStreams::pout_base << "Assemble right hand side (CPU/wall) " @@ -217,8 +222,8 @@ template void LaplaceProblem::solve() { - Timer time; - MGTransferMatrixFree mg_transfer(mg_constrained_dofs); + Timer time; + dealii::MGTransferMatrixFree mg_transfer(mg_constrained_dofs); mg_transfer.build(dof_handler); setup_time += time.wall_time(); conditionalOStreams::pout_base << "MG build transfer time (CPU/wall) " @@ -226,10 +231,12 @@ LaplaceProblem::solve() time.restart(); using SmootherType = - PreconditionChebyshev>; - mg::SmootherRelaxation> - mg_smoother; - MGLevelObject smoother_data; + dealii::PreconditionChebyshev>; + dealii::mg::SmootherRelaxation> + mg_smoother; + dealii::MGLevelObject smoother_data; smoother_data.resize(0, triangulation.n_global_levels() - 1); for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level) { @@ -242,7 +249,7 @@ LaplaceProblem::solve() else { smoother_data[0].smoothing_range = 1e-3; - smoother_data[0].degree = numbers::invalid_unsigned_int; + smoother_data[0].degree = dealii::numbers::invalid_unsigned_int; smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m(); } mg_matrices[level].compute_diagonal(); @@ -251,19 +258,21 @@ LaplaceProblem::solve() } mg_smoother.initialize(mg_matrices, smoother_data); - MGCoarseGridApplySmoother> mg_coarse; + dealii::MGCoarseGridApplySmoother> + mg_coarse; mg_coarse.initialize(mg_smoother); - mg::Matrix> mg_matrix(mg_matrices); + dealii::mg::Matrix> mg_matrix( + mg_matrices); - MGLevelObject> + dealii::MGLevelObject> mg_interface_matrices; mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1); for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level) { mg_interface_matrices[level].initialize(mg_matrices[level]); } - mg::Matrix> mg_interface( + dealii::mg::Matrix> mg_interface( mg_interface_matrices); Multigrid> mg(mg_matrix, @@ -273,13 +282,13 @@ LaplaceProblem::solve() mg_smoother); mg.set_edge_matrices(mg_interface, mg_interface); - PreconditionMG, - MGTransferMatrixFree> + dealii::PreconditionMG, + dealii::MGTransferMatrixFree> preconditioner(dof_handler, mg, mg_transfer); - SolverControl solver_control(100, 1e-12 * system_rhs.l2_norm()); - SolverCG> cg(solver_control); + dealii::SolverControl solver_control(100, 1e-12 * system_rhs.l2_norm()); + dealii::SolverCG> cg(solver_control); setup_time += time.wall_time(); conditionalOStreams::pout_base << "MG build smoother time (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << "s\n"; @@ -309,15 +318,15 @@ LaplaceProblem::output_results(unsigned int cycle) const return; } - DataOut data_out; + dealii::DataOut data_out; solution.update_ghost_values(); data_out.attach_dof_handler(dof_handler); data_out.add_data_vector(solution, "solution"); data_out.build_patches(mapping); - DataOutBase::VtkFlags flags; - flags.compression_level = DataOutBase::CompressionLevel::best_speed; + dealii::DataOutBase::VtkFlags flags; + flags.compression_level = dealii::DataOutBase::CompressionLevel::best_speed; data_out.set_flags(flags); data_out.write_vtu_with_pvtu_record("./", "solution", cycle, MPI_COMM_WORLD, 3); @@ -329,16 +338,14 @@ template void LaplaceProblem::run() { - { - const unsigned int n_vect_doubles = VectorizedArray::size(); - const unsigned int n_vect_bits = 8 * sizeof(double) * n_vect_doubles; + const unsigned int n_vect_doubles = dealii::VectorizedArray::size(); + const unsigned int n_vect_bits = 8 * sizeof(double) * n_vect_doubles; - conditionalOStreams::pout_base - << "Vectorization over " << n_vect_doubles << " doubles = " << n_vect_bits - << " bits (" << Utilities::System::get_current_vectorization_level() << ')' << '\n'; - } + conditionalOStreams::pout_base + << "Vectorization over " << n_vect_doubles << " doubles = " << n_vect_bits + << " bits (" << Utilities::System::get_current_vectorization_level() << ')' << '\n'; - GridGenerator::hyper_cube(triangulation, 0.0, 1.0); + dealii::GridGenerator::hyper_cube(triangulation, 0.0, 1.0); triangulation.refine_global(6); setup_system();