From edd6945f582fae7f96284be83ef39df6e4cbe514 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Thu, 5 Sep 2024 14:29:16 +0200 Subject: [PATCH] Add utility to compute errors --- examples/diffusion_reaction.cc | 18 +- examples/poisson.cc | 131 +++- include/poly_utils.h | 540 +++++++++++----- test/polydeal/exact_solutions_dgp.cc | 749 +++++++++++++++++++++++ test/polydeal/exact_solutions_dgp.output | 4 + 5 files changed, 1275 insertions(+), 167 deletions(-) create mode 100644 test/polydeal/exact_solutions_dgp.cc create mode 100644 test/polydeal/exact_solutions_dgp.output diff --git a/examples/diffusion_reaction.cc b/examples/diffusion_reaction.cc index 0f7e34e3..8bdb5a6a 100644 --- a/examples/diffusion_reaction.cc +++ b/examples/diffusion_reaction.cc @@ -718,6 +718,8 @@ DiffusionReactionProblem::solve() interpolated_solution, completely_distributed_solution); + // Use other method to compute the error +#ifdef FALSE Vector cellwise_error(completely_distributed_solution.size()); VectorTools::integrate_difference(ah->output_dh, interpolated_solution, @@ -740,9 +742,19 @@ DiffusionReactionProblem::solve() VectorTools::compute_global_error(tria_pft, cellwise_error, VectorTools::NormType::H1_seminorm); - - pcout << "L2 error (exponential solution): " << error << std::endl; - pcout << "Semi H1 error (exponential solution): " << semiH1error << std::endl; +#endif + + std::vector global_errors; + PolyUtils::compute_global_error(*ah, + completely_distributed_solution, + Solution(), + {VectorTools::L2_norm, + VectorTools::H1_seminorm}, + global_errors); + + pcout << "L2 error (exponential solution): " << global_errors[0] << std::endl; + pcout << "Semi H1 error (exponential solution): " << global_errors[1] + << std::endl; } diff --git a/examples/poisson.cc b/examples/poisson.cc index eae66191..29f6e2e4 100644 --- a/examples/poisson.cc +++ b/examples/poisson.cc @@ -11,6 +11,8 @@ // ----------------------------------------------------------------------------- +#include + #include #include #include @@ -34,6 +36,43 @@ #include + +struct ConvergenceInfo +{ + ConvergenceInfo() = default; + void + add(const std::pair> + &dofs_and_errs) + { + vec_data.push_back(dofs_and_errs); + } + + + + void + print() + { + Assert(vec_data.size() > 0, ExcInternalError()); + for (const auto &dof_and_errs : vec_data) + std::cout << std::left << std::setw(24) << std::scientific + << "N DoFs: " << dof_and_errs.first << std::endl; + + for (const auto &dof_and_errs : vec_data) + std::cout << std::left << std::setw(24) << std::scientific + << "L2 error: " << dof_and_errs.second.first << std::endl; + for (const auto &dof_and_errs : vec_data) + std::cout << std::left << std::setw(24) << std::scientific + << "H1 error: " << dof_and_errs.second.second << std::endl; + } + + + + std::vector>> + vec_data; +}; + + + template constexpr T constexpr_pow(T num, unsigned int pow) @@ -261,6 +300,11 @@ class SolutionProduct : public Function virtual Tensor<1, dim> gradient(const Point &p, const unsigned int component = 0) const override; + + virtual void + gradient_list(const std::vector> &points, + std::vector> &gradients, + const unsigned int /*component*/) const override; }; template @@ -293,6 +337,18 @@ SolutionProduct::value_list(const std::vector> &points, +template +void +SolutionProduct::gradient_list(const std::vector> &points, + std::vector> &gradients, + const unsigned int /*component*/) const +{ + for (unsigned int i = 0; i < gradients.size(); ++i) + gradients[i] = this->gradient(points[i]); +} + + + template class SolutionProductSine : public Function { @@ -389,12 +445,20 @@ class Poisson void run(); + types::global_dof_index + get_n_dofs() const; + + std::pair + get_error() const; + GridType grid_type; PartitionerType partitioner_type; SolutionType solution_type; unsigned int extraction_level; unsigned int n_subdomains; double penalty_constant = 60.; // 10*(p+1)(p+d) for p = 1 and d = 2 => 60 + double l2_err; + double semih1_err; }; @@ -461,7 +525,7 @@ Poisson::make_grid() else { GridGenerator::hyper_cube(tria, 0., 1.); - tria.refine_global(9); + tria.refine_global(8); } std::cout << "Size of tria: " << tria.n_active_cells() << std::endl; cached_tria = std::make_unique>(tria, mapping); @@ -1041,7 +1105,24 @@ Poisson::output_results() "agglo_idx", DataOut::type_cell_data); + data_out.build_patches(mapping); + data_out.write_vtu(output); + // Compute L2 and semiH1 norm of error + std::vector errors; + PolyUtils::compute_global_error(*ah, + solution, + *analytical_solution, + {VectorTools::L2_norm, + VectorTools::H1_seminorm}, + errors); + l2_err = errors[0]; + semih1_err = errors[1]; + std::cout << "Error (L2): " << l2_err << std::endl; + std::cout << "Error (H1): " << semih1_err << std::endl; + +// Old version +#ifdef FALSE { // L2 error Vector difference_per_cell(tria.n_active_cells()); @@ -1085,12 +1166,30 @@ Poisson::output_results() std::cout << "H1 seminorm:" << H1_seminorm << std::endl; } - - data_out.build_patches(mapping); - data_out.write_vtu(output); +#endif } } + + +template +inline types::global_dof_index +Poisson::get_n_dofs() const +{ + return ah->n_dofs(); +} + + + +template +inline std::pair +Poisson::get_error() const +{ + return std::make_pair(l2_err, semih1_err); +} + + + template void Poisson::run() @@ -1107,32 +1206,28 @@ Poisson::run() int main() { - std::cout << "Benchmarking with Rtree:" << std::endl; - // Testing p-convergence + ConvergenceInfo convergence_info; std::cout << "Testing p-convergence" << std::endl; { - // for (unsigned int fe_degree : {1, 2, 3, 4, 5}) - for (unsigned int fe_degree : {1}) + for (unsigned int fe_degree : {1, 2, 3, 4, 5, 6}) { std::cout << "Fe degree: " << fe_degree << std::endl; - Poisson<2> poisson_problem{GridType::grid_generator, - PartitionerType::rtree, + Poisson<2> poisson_problem{GridType::unstructured, + PartitionerType::metis, SolutionType::product_sine, - 5 /*extaction_level*/, - 0, + 0 /*extraction_level*/, + 364, fe_degree}; poisson_problem.run(); + convergence_info.add( + std::make_pair>( + poisson_problem.get_n_dofs(), poisson_problem.get_error())); } } + convergence_info.print(); std::cout << std::endl; return 0; - // for (const unsigned int target_partitions : - // {16, 64, 256, 1024, 4096}) // ball - // // {6, 23, 91, 364, 1456, 5824}) // - // // unstructured {16, 64, 256, 1024, - // // 4096}) - // // // structured square } diff --git a/include/poly_utils.h b/include/poly_utils.h index f96ef74d..ae9a7f7e 100644 --- a/include/poly_utils.h +++ b/include/poly_utils.h @@ -41,6 +41,8 @@ #include #include +#include + #include #include #include @@ -835,19 +837,218 @@ namespace dealii::PolyUtils + namespace internal + { + /** + * Same as the public free function with the same name, but storing + * explicitly the interpolation matrix and performing interpolation through + * matrix-vector product. + */ + template + void + interpolate_to_fine_grid( + const AgglomerationHandler &agglomeration_handler, + VectorType &dst, + const VectorType &src) + { + Assert((dim == spacedim), ExcNotImplemented()); + Assert( + dst.size() == 0, + ExcMessage( + "The destination vector must the empt upon calling this function.")); + + using NumberType = typename VectorType::value_type; + constexpr bool is_trilinos_vector = + std::is_same_v; + using MatrixType = std::conditional_t>; + + MatrixType interpolation_matrix; + + [[maybe_unused]] + typename std::conditional_t + sp; + + // Get some info from the handler + const DoFHandler &agglo_dh = + agglomeration_handler.agglo_dh; + + DoFHandler *output_dh = + const_cast *>( + &agglomeration_handler.output_dh); + const FiniteElement &fe = agglomeration_handler.get_fe(); + const Triangulation &tria = + agglomeration_handler.get_triangulation(); + const auto &bboxes = agglomeration_handler.get_local_bboxes(); + + // Setup an auxiliary DoFHandler for output purposes + output_dh->reinit(tria); + output_dh->distribute_dofs(fe); + + const IndexSet &locally_owned_dofs = output_dh->locally_owned_dofs(); + const IndexSet locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(*output_dh); + + const IndexSet &locally_owned_dofs_agglo = agglo_dh.locally_owned_dofs(); + + + DynamicSparsityPattern dsp(output_dh->n_dofs(), + agglo_dh.n_dofs(), + locally_relevant_dofs); + + std::vector agglo_dof_indices(fe.dofs_per_cell); + std::vector standard_dof_indices( + fe.dofs_per_cell); + std::vector output_dof_indices(fe.dofs_per_cell); + + Quadrature quad(fe.get_unit_support_points()); + FEValues output_fe_values(fe, + quad, + update_quadrature_points); + + for (const auto &cell : agglo_dh.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if (agglomeration_handler.is_master_cell(cell)) + { + auto slaves = agglomeration_handler.get_slaves_of_idx( + cell->active_cell_index()); + slaves.emplace_back(cell); + + cell->get_dof_indices(agglo_dof_indices); + + for (const auto &slave : slaves) + { + // addd master-slave relationship + const auto slave_output = + slave->as_dof_handler_iterator(*output_dh); + slave_output->get_dof_indices(output_dof_indices); + for (const auto row : output_dof_indices) + dsp.add_entries(row, + agglo_dof_indices.begin(), + agglo_dof_indices.end()); + } + } + } + + + const auto assemble_interpolation_matrix = [&]() { + FullMatrix local_matrix(fe.dofs_per_cell, fe.dofs_per_cell); + std::vector> reference_q_points(fe.dofs_per_cell); + + // Dummy AffineConstraints, only needed for loc2glb + AffineConstraints c; + c.close(); + + for (const auto &cell : agglo_dh.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if (agglomeration_handler.is_master_cell(cell)) + { + auto slaves = agglomeration_handler.get_slaves_of_idx( + cell->active_cell_index()); + slaves.emplace_back(cell); + + cell->get_dof_indices(agglo_dof_indices); + + const types::global_cell_index polytope_index = + agglomeration_handler.cell_to_polytope_index(cell); + + // Get the box of this agglomerate. + const BoundingBox &box = bboxes[polytope_index]; + + for (const auto &slave : slaves) + { + // add master-slave relationship + const auto slave_output = + slave->as_dof_handler_iterator(*output_dh); + + slave_output->get_dof_indices(output_dof_indices); + output_fe_values.reinit(slave_output); + + local_matrix = 0.; + + const auto &q_points = + output_fe_values.get_quadrature_points(); + for (const auto i : output_fe_values.dof_indices()) + { + const auto &p = box.real_to_unit(q_points[i]); + for (const auto j : output_fe_values.dof_indices()) + { + local_matrix(i, j) = fe.shape_value(j, p); + } + } + c.distribute_local_to_global(local_matrix, + output_dof_indices, + agglo_dof_indices, + interpolation_matrix); + } + } + } + }; + + + if constexpr (std::is_same_v) + { + const MPI_Comm &communicator = tria.get_communicator(); + SparsityTools::distribute_sparsity_pattern(dsp, + locally_owned_dofs, + communicator, + locally_relevant_dofs); + + interpolation_matrix.reinit(locally_owned_dofs, + locally_owned_dofs_agglo, + dsp, + communicator); + dst.reinit(locally_owned_dofs); + assemble_interpolation_matrix(); + } + else if constexpr (std::is_same_v>) + { + sp.copy_from(dsp); + interpolation_matrix.reinit(sp); + dst.reinit(output_dh->n_dofs()); + assemble_interpolation_matrix(); + } + else + { + // PETSc, LA::d::v options not implemented. + (void)agglomeration_handler; + (void)dst; + (void)src; + AssertThrow(false, ExcNotImplemented()); + } + + // If tria is distributed + if (dynamic_cast *>( + &tria) != nullptr) + interpolation_matrix.compress(VectorOperation::add); + + // Finally, perform the interpolation. + interpolation_matrix.vmult(dst, src); + } + } // namespace internal + + + /** * Given a vector @p src, typically the solution stemming after the * agglomerate problem has been solved, this function interpolates @p src - * onto the finer grid and stores the result in vector @p dst. + * onto the finer grid and stores the result in vector @p dst. The last + * argument @p on_the_fly does not build any interpolation matrix and allows + * computing the entries in @p dst in a matrix-free fashion. * - * @note Supported parallel types are TrilinosWrappers::SparseMatrix and TrilinosWrappers::MPI::Vector + * @note Supported parallel types are TrilinosWrappers::SparseMatrix and + * TrilinosWrappers::MPI::Vector. */ template void interpolate_to_fine_grid( const AgglomerationHandler &agglomeration_handler, VectorType &dst, - const VectorType &src) + const VectorType &src, + const bool on_the_fly = true) { Assert((dim == spacedim), ExcNotImplemented()); Assert( @@ -856,179 +1057,110 @@ namespace dealii::PolyUtils "The destination vector must the empt upon calling this function.")); using NumberType = typename VectorType::value_type; - constexpr bool is_trilinos_vector = + static constexpr bool is_trilinos_vector = std::is_same_v; - using MatrixType = std::conditional_t>; - MatrixType interpolation_matrix; - - [[maybe_unused]] - typename std::conditional_t - sp; + static constexpr bool is_supported_vector = + std::is_same_v> || is_trilinos_vector; + static_assert(is_supported_vector); - // Get some info from the handler - const DoFHandler &agglo_dh = agglomeration_handler.agglo_dh; - - DoFHandler *output_dh = - const_cast *>(&agglomeration_handler.output_dh); - const FiniteElement &fe = agglomeration_handler.get_fe(); - const Triangulation &tria = - agglomeration_handler.get_triangulation(); - const auto &bboxes = agglomeration_handler.get_local_bboxes(); - - // Setup an auxiliary DoFHandler for output purposes - output_dh->reinit(tria); - output_dh->distribute_dofs(fe); - - const IndexSet &locally_owned_dofs = output_dh->locally_owned_dofs(); - const IndexSet locally_relevant_dofs = - DoFTools::extract_locally_relevant_dofs(*output_dh); + // First, check for an easy return + if (on_the_fly == false) + { + return internal::interpolate_to_fine_grid(agglomeration_handler, + dst, + src); + } + else + { + // otherwise, do not create any matrix + const Triangulation &tria = + agglomeration_handler.get_triangulation(); + const FiniteElement &original_fe = + agglomeration_handler.get_fe(); - const IndexSet &locally_owned_dofs_agglo = agglo_dh.locally_owned_dofs(); + // We use DGQ nodal elements of the same degree as the ones in the + // agglomeration handler to generate the output also in the case in + // which different elements are used, such as DGP. + FE_DGQ output_fe(original_fe.degree); + DoFHandler &output_dh = + const_cast &>(agglomeration_handler.output_dh); - DynamicSparsityPattern dsp(output_dh->n_dofs(), - agglo_dh.n_dofs(), - locally_relevant_dofs); + output_dh.reinit(tria); + output_dh.distribute_dofs(output_fe); - std::vector agglo_dof_indices(fe.dofs_per_cell); - std::vector standard_dof_indices(fe.dofs_per_cell); - std::vector output_dof_indices(fe.dofs_per_cell); + if constexpr (std::is_same_v) + { + const IndexSet &locally_owned_dofs = output_dh.locally_owned_dofs(); + dst.reinit(locally_owned_dofs); + } + else if constexpr (std::is_same_v>) + { + dst.reinit(output_dh.n_dofs()); + } + else + { + // PETSc, LA::d::v options not implemented. + (void)agglomeration_handler; + (void)dst; + (void)src; + AssertThrow(false, ExcNotImplemented()); + } - Quadrature quad(fe.get_unit_support_points()); - FEValues output_fe_values(fe, - quad, - update_quadrature_points); - for (const auto &cell : agglo_dh.active_cell_iterators()) - if (cell->is_locally_owned()) - { - if (agglomeration_handler.is_master_cell(cell)) - { - auto slaves = agglomeration_handler.get_slaves_of_idx( - cell->active_cell_index()); - slaves.emplace_back(cell); - cell->get_dof_indices(agglo_dof_indices); + const unsigned int dofs_per_cell = + agglomeration_handler.n_dofs_per_cell(); + const unsigned int output_dofs_per_cell = output_fe.n_dofs_per_cell(); + Quadrature quad(output_fe.get_unit_support_points()); + FEValues output_fe_values(output_fe, + quad, + update_quadrature_points); - for (const auto &slave : slaves) - { - // addd master-slave relationship - const auto slave_output = - slave->as_dof_handler_iterator(*output_dh); - slave_output->get_dof_indices(output_dof_indices); - for (const auto row : output_dof_indices) - dsp.add_entries(row, - agglo_dof_indices.begin(), - agglo_dof_indices.end()); - } - } - } + std::vector local_dof_indices(dofs_per_cell); + std::vector local_dof_indices_output( + output_dofs_per_cell); - - const auto assemble_interpolation_matrix = [&]() { - FullMatrix local_matrix(fe.dofs_per_cell, fe.dofs_per_cell); - std::vector> reference_q_points(fe.dofs_per_cell); - - // Dummy AffineConstraints, only needed for loc2glb - AffineConstraints c; - c.close(); - - for (const auto &cell : agglo_dh.active_cell_iterators()) - if (cell->is_locally_owned()) + const auto &bboxes = agglomeration_handler.get_local_bboxes(); + for (const auto &polytope : agglomeration_handler.polytope_iterators()) { - if (agglomeration_handler.is_master_cell(cell)) + if (polytope->is_locally_owned()) { - auto slaves = agglomeration_handler.get_slaves_of_idx( - cell->active_cell_index()); - slaves.emplace_back(cell); - - cell->get_dof_indices(agglo_dof_indices); + polytope->get_dof_indices(local_dof_indices); + const BoundingBox &box = bboxes[polytope->index()]; - const types::global_cell_index polytope_index = - agglomeration_handler.cell_to_polytope_index(cell); - - // Get the box of this agglomerate. - const BoundingBox &box = bboxes[polytope_index]; - - for (const auto &slave : slaves) + const auto &deal_cells = + polytope->get_agglomerate(); // fine deal.II cells + for (const auto &cell : deal_cells) { - // add master-slave relationship - const auto slave_output = - slave->as_dof_handler_iterator(*output_dh); - - slave_output->get_dof_indices(output_dof_indices); + const auto slave_output = cell->as_dof_handler_iterator( + agglomeration_handler.output_dh); + slave_output->get_dof_indices(local_dof_indices_output); output_fe_values.reinit(slave_output); - local_matrix = 0.; - - const auto &q_points = + const auto &qpoints = output_fe_values.get_quadrature_points(); - for (const auto i : output_fe_values.dof_indices()) + + for (unsigned int j = 0; j < output_dofs_per_cell; ++j) { - const auto &p = box.real_to_unit(q_points[i]); - for (const auto j : output_fe_values.dof_indices()) - { - local_matrix(i, j) = fe.shape_value(j, p); - } + const auto &ref_qpoint = box.real_to_unit(qpoints[j]); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + dst(local_dof_indices_output[j]) += + src(local_dof_indices[i]) * + original_fe.shape_value(i, ref_qpoint); } - c.distribute_local_to_global(local_matrix, - output_dof_indices, - agglo_dof_indices, - interpolation_matrix); } } } - }; - - - if constexpr (std::is_same_v) - { - const MPI_Comm &communicator = tria.get_communicator(); - SparsityTools::distribute_sparsity_pattern(dsp, - locally_owned_dofs, - communicator, - locally_relevant_dofs); - - interpolation_matrix.reinit(locally_owned_dofs, - locally_owned_dofs_agglo, - dsp, - communicator); - dst.reinit(locally_owned_dofs); - assemble_interpolation_matrix(); - } - else if constexpr (std::is_same_v>) - { - sp.copy_from(dsp); - interpolation_matrix.reinit(sp); - dst.reinit(output_dh->n_dofs()); - assemble_interpolation_matrix(); } - else - { - // PETSc, LA::d::v options not implemented. - (void)agglomeration_handler; - (void)dst; - (void)src; - AssertThrow(false, ExcNotImplemented()); - } - - // If tria is distributed - if (dynamic_cast *>( - &tria) != nullptr) - interpolation_matrix.compress(VectorOperation::add); - - // Finally, perform the interpolation. - interpolation_matrix.vmult(dst, src); } /** - * Construct the interpolation matrix from the DG space defined the polytopic - * elements + * Construct the interpolation matrix from the DG space defined the + * polytopic elements * defined in @p agglomeration_handler to the DG space defined on the DoFHandler associated * to standard shapes. The interpolation matrix is assumed to be * default-constructed and is filled inside this function. @@ -1200,6 +1332,122 @@ namespace dealii::PolyUtils + /** + * Similar to VectorTools::compute_global_error(), but customized for + * polytopic elements. Aside from the solution vector and a reference + * function, this function takes in addition a vector @p norms with types + * VectorTools::NormType to be computed and later stored in the last + * argument @p global_errors. + * In case of a parallel vector, the local errors are collected over each + * processor and later a classical reduction operation is performed. + */ + template + void + compute_global_error(const AgglomerationHandler &agglomeration_handler, + const VectorType &solution, + const Function &exact_solution, + const std::vector &norms, + std::vector &global_errors) + { + Assert(solution.size() > 0, + ExcNotImplemented( + "Solution vector must be non-empty upon calling this function.")); + Assert(std::any_of(norms.cbegin(), + norms.cend(), + [](VectorTools::NormType norm_type) { + return (norm_type == + VectorTools::NormType::H1_seminorm || + norm_type == VectorTools::NormType::L2_norm); + }), + ExcMessage("Norm type not supported")); + global_errors.resize(norms.size()); + std::fill(global_errors.begin(), global_errors.end(), 0.); + + // Vector storing errors local to the current processor. + std::vector local_errors(norms.size()); + std::fill(local_errors.begin(), local_errors.end(), 0.); + + // Get some info from the handler + const unsigned int dofs_per_cell = agglomeration_handler.n_dofs_per_cell(); + + const bool compute_semi_H1 = + std::any_of(norms.cbegin(), + norms.cend(), + [](VectorTools::NormType norm_type) { + return norm_type == VectorTools::NormType::H1_seminorm; + }); + + std::vector local_dof_indices(dofs_per_cell); + for (const auto &polytope : agglomeration_handler.polytope_iterators()) + { + if (polytope->is_locally_owned()) + { + const auto &agglo_values = agglomeration_handler.reinit(polytope); + polytope->get_dof_indices(local_dof_indices); + + const auto &q_points = agglo_values.get_quadrature_points(); + const unsigned int n_qpoints = q_points.size(); + std::vector analyical_sol_at_qpoints(n_qpoints); + exact_solution.value_list(q_points, analyical_sol_at_qpoints); + std::vector> grad_analyical_sol_at_qpoints( + n_qpoints); + + if (compute_semi_H1) + exact_solution.gradient_list(q_points, + grad_analyical_sol_at_qpoints); + + for (unsigned int q_index : agglo_values.quadrature_point_indices()) + { + double solution_at_qpoint = 0.; + Tensor<1, dim> grad_solution_at_qpoint; + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + solution_at_qpoint += solution(local_dof_indices[i]) * + agglo_values.shape_value(i, q_index); + + if (compute_semi_H1) + grad_solution_at_qpoint += + solution(local_dof_indices[i]) * + agglo_values.shape_grad(i, q_index); + } + // L2 + local_errors[0] += std::pow((analyical_sol_at_qpoints[q_index] - + solution_at_qpoint), + 2) * + agglo_values.JxW(q_index); + + // H1 seminorm + if (compute_semi_H1) + for (unsigned int d = 0; d < dim; ++d) + local_errors[1] += + std::pow((grad_analyical_sol_at_qpoints[q_index][d] - + grad_solution_at_qpoint[d]), + 2) * + agglo_values.JxW(q_index); + } + } + } + + // Perform reduction and take sqrt of each error + global_errors[0] = Utilities::MPI::reduce( + local_errors[0], + agglomeration_handler.get_triangulation().get_mpi_communicator(), + [](const double a, const double b) { return a + b; }); + + global_errors[0] = std::sqrt(global_errors[0]); + + if (compute_semi_H1) + { + global_errors[1] = Utilities::MPI::reduce( + local_errors[1], + agglomeration_handler.get_triangulation().get_mpi_communicator(), + [](const double a, const double b) { return a + b; }); + global_errors[1] = std::sqrt(global_errors[1]); + } + } + + + /** * Utility function that builds the multilevel hierarchy from the tree level * @p starting_level. This function fills the vector of diff --git a/test/polydeal/exact_solutions_dgp.cc b/test/polydeal/exact_solutions_dgp.cc new file mode 100644 index 00000000..11c82372 --- /dev/null +++ b/test/polydeal/exact_solutions_dgp.cc @@ -0,0 +1,749 @@ +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + + +// Same as exact_solutions.cc, but with a DGP space. + +static constexpr double TOL = 1e-14; + +enum class SolutionType +{ + LinearSolution, + QuadraticSolution +}; + + +template +class LinearFunction : public Function +{ +public: + LinearFunction(const std::vector &coeffs) + { + Assert(coeffs.size() <= dim, ExcMessage("Wrong size!")); + coefficients.resize(coeffs.size()); + for (size_t i = 0; i < coeffs.size(); i++) + coefficients[i] = coeffs[i]; + } + virtual double + value(const Point &p, const unsigned int component = 0) const override; + + virtual void + value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const override; + + std::vector coefficients; +}; + +template +double +LinearFunction::value(const Point &p, const unsigned int) const +{ + double value = 0.; + for (size_t i = 0; i < coefficients.size(); i++) + value += coefficients[i] * p[i]; + return value; +} + + + +template +void +LinearFunction::value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const +{ + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = value(points[i]); +} + + +template +class RightHandSide : public Function +{ +public: + RightHandSide(const SolutionType &solution_type) + : Function() + { + sol_type = solution_type; + } + + virtual void + value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const override; + + SolutionType sol_type; +}; + + +template +void +RightHandSide::value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const +{ + (void)points; + if (sol_type == SolutionType::LinearSolution) + { + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = 0.; //-4.; // ball, radial solution + } + else + { + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = -4.; // quadratic solution + } +} + +template +class SolutionLinear : public Function +{ +public: + SolutionLinear() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const override; + + virtual void + value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const override; + + virtual Tensor<1, dim> + gradient(const Point &p, + const unsigned int component = 0) const override; +}; + +template +double +SolutionLinear::value(const Point &p, const unsigned int) const +{ + return p[0] + p[1] - 1; // linear +} + +template +Tensor<1, dim> +SolutionLinear::gradient(const Point &p, const unsigned int) const +{ + Assert(dim == 2, ExcMessage("This test only works in 2D.")); + (void)p; + Tensor<1, dim> return_value; + return_value[0] = 1.; + return_value[1] = 1.; + return return_value; +} + + +template +void +SolutionLinear::value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const +{ + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = this->value(points[i]); +} + + + +template +class SolutionQuadratic : public Function +{ +public: + SolutionQuadratic() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const override; + + virtual void + value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const override; + + virtual Tensor<1, dim> + gradient(const Point &p, + const unsigned int component = 0) const override; +}; + +template +double +SolutionQuadratic::value(const Point &p, const unsigned int) const +{ + return p[0] * p[0] + p[1] * p[1] - 1; +} + +template +Tensor<1, dim> +SolutionQuadratic::gradient(const Point &p, const unsigned int) const +{ + Assert(dim == 2, ExcMessage("This test only works in 2D.")); + (void)p; + Tensor<1, dim> return_value; + return_value[0] = 2 * p[0]; + return_value[1] = 2 * p[1]; + return return_value; +} + + + +template +void +SolutionQuadratic::value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const +{ + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = this->value(points[i]); +} + + + +template +class Poisson +{ +private: + void + make_grid(); + void + setup_agglomeration(); + void + assemble_system(); + void + solve(); + void + output_results(); + + + Triangulation tria; + MappingQ mapping; + FE_DGP *dg_fe; + std::unique_ptr> ah; + // no hanging node in DG discretization, we define an AffineConstraints + // object + // so we can use the distribute_local_to_global() directly. + AffineConstraints constraints; + SparsityPattern sparsity; + DynamicSparsityPattern dsp; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + std::unique_ptr> cached_tria; + std::unique_ptr> rhs_function; + Function *analytical_solution; + +public: + Poisson(const SolutionType &solution_type); + ~Poisson(); + + void + run(); + + double penalty_constant = 10; + SolutionType sol_type; +}; + + + +template +Poisson::Poisson(const SolutionType &solution_type) + : mapping(1) + , sol_type(solution_type) +{ + if (sol_type == SolutionType::LinearSolution) + { + dg_fe = new FE_DGP{1}; + analytical_solution = new SolutionLinear(); + } + else + { + dg_fe = new FE_DGP{2}; + analytical_solution = new SolutionQuadratic(); + } +} + + + +template +Poisson::~Poisson() +{ + delete dg_fe; + delete analytical_solution; +} + + + +template +void +Poisson::make_grid() +{ + GridIn grid_in; + GridGenerator::hyper_cube(tria, 0, 1); + tria.refine_global(2); + GridTools::distort_random(0.25, tria); + + std::cout << "Size of tria: " << tria.n_active_cells() << std::endl; + cached_tria = std::make_unique>(tria, mapping); + rhs_function = std::make_unique>(sol_type); + + constraints.close(); +} + +template +void +Poisson::setup_agglomeration() +{ + ah = std::make_unique>(*cached_tria); + + std::vector idxs_to_be_agglomerated = {0, 1, 2, 3}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated, + cells_to_be_agglomerated); + + std::vector idxs_to_be_agglomerated2 = {4, 5, 6, 7}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated2; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated2, + cells_to_be_agglomerated2); + + std::vector idxs_to_be_agglomerated3 = {8, + 9, + 10, + 11}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated3; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated3, + cells_to_be_agglomerated3); + + std::vector idxs_to_be_agglomerated4 = {12, + 13, + 14, + 15}; + + std::vector::active_cell_iterator> + cells_to_be_agglomerated4; + PolyUtils::collect_cells_for_agglomeration(tria, + idxs_to_be_agglomerated4, + cells_to_be_agglomerated4); + + // Agglomerate the cells just stored + ah->define_agglomerate(cells_to_be_agglomerated); + ah->define_agglomerate(cells_to_be_agglomerated2); + ah->define_agglomerate(cells_to_be_agglomerated3); + ah->define_agglomerate(cells_to_be_agglomerated4); + + ah->distribute_agglomerated_dofs(*dg_fe); + ah->create_agglomeration_sparsity_pattern(dsp); + sparsity.copy_from(dsp); +} + + + +template +void +Poisson::assemble_system() +{ + system_matrix.reinit(sparsity); + solution.reinit(ah->n_dofs()); + system_rhs.reinit(ah->n_dofs()); + + const unsigned int quadrature_degree = 2 * dg_fe->get_degree() + 1; + const unsigned int face_quadrature_degree = 2 * dg_fe->get_degree() + 1; + ah->initialize_fe_values(QGauss(quadrature_degree), + update_gradients | update_JxW_values | + update_quadrature_points | update_JxW_values | + update_values, + QGauss(face_quadrature_degree)); + + const unsigned int dofs_per_cell = ah->n_dofs_per_cell(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + // Next, we define the four dofsxdofs matrices needed to assemble jumps and + // averages. + FullMatrix M11(dofs_per_cell, dofs_per_cell); + FullMatrix M12(dofs_per_cell, dofs_per_cell); + FullMatrix M21(dofs_per_cell, dofs_per_cell); + FullMatrix M22(dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + + LinearFunction linear_func{{1, 1}}; + double test_integral = 0.; + double test_bdary = 0.; + double test_volume = 0.; + for (const auto &polytope : ah->polytope_iterators()) + { + // local_volume = 0.; + cell_matrix = 0; + cell_rhs = 0; + const auto &agglo_values = ah->reinit(polytope); + + const auto &q_points = agglo_values.get_quadrature_points(); + + const unsigned int n_qpoints = q_points.size(); + std::vector rhs(n_qpoints); + rhs_function->value_list(q_points, rhs); + std::vector linear_values(n_qpoints); + linear_func.value_list(q_points, linear_values, 1); + + for (unsigned int q_index : agglo_values.quadrature_point_indices()) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + cell_matrix(i, j) += agglo_values.shape_grad(i, q_index) * + agglo_values.shape_grad(j, q_index) * + agglo_values.JxW(q_index); + } + cell_rhs(i) += agglo_values.shape_value(i, q_index) * + rhs[q_index] * agglo_values.JxW(q_index); + } + test_integral += linear_values[q_index] * agglo_values.JxW(q_index); + test_volume += agglo_values.JxW(q_index); + } + + polytope->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + + // Face terms + const unsigned int n_faces = polytope->n_faces(); + AssertThrow(n_faces > 0, ExcMessage("Invalid element!")); + + + + // auto polygon_boundary_vertices = ah->polytope_boundary(cell); + for (unsigned int f = 0; f < n_faces; ++f) + { + if (polytope->at_boundary(f)) + { + const auto &fe_face = ah->reinit(polytope, f); + + const unsigned int dofs_per_cell = fe_face.dofs_per_cell; + std::vector local_dof_indices_bdary_cell( + dofs_per_cell); + + const auto &face_q_points = fe_face.get_quadrature_points(); + std::vector analytical_solution_values( + face_q_points.size()); + analytical_solution->value_list(face_q_points, + analytical_solution_values, + 1); + + // Get normal vectors seen from each agglomeration. + const auto &normals = fe_face.get_normal_vectors(); + + const double penalty = + penalty_constant / std::fabs(polytope->diameter()); + + cell_matrix = 0.; + cell_rhs = 0.; + for (unsigned int q_index : fe_face.quadrature_point_indices()) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + cell_matrix(i, j) += + (-fe_face.shape_value(i, q_index) * + fe_face.shape_grad(j, q_index) * + normals[q_index] - + fe_face.shape_grad(i, q_index) * normals[q_index] * + fe_face.shape_value(j, q_index) + + (penalty)*fe_face.shape_value(i, q_index) * + fe_face.shape_value(j, q_index)) * + fe_face.JxW(q_index); + } + cell_rhs(i) += + (penalty * analytical_solution_values[q_index] * + fe_face.shape_value(i, q_index) - + fe_face.shape_grad(i, q_index) * normals[q_index] * + analytical_solution_values[q_index]) * + fe_face.JxW(q_index); + } + + test_bdary += fe_face.JxW(q_index); + } + + // distribute DoFs + polytope->get_dof_indices(local_dof_indices_bdary_cell); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + else + { + const auto &neigh_polytope = polytope->neighbor(f); + + // This is necessary to loop over internal faces only once. + if (polytope->index() < neigh_polytope->index()) + { + unsigned int nofn = + polytope->neighbor_of_agglomerated_neighbor(f); + + const auto &fe_faces = + ah->reinit_interface(polytope, neigh_polytope, f, nofn); + + const auto &fe_faces0 = fe_faces.first; + const auto &fe_faces1 = fe_faces.second; + +#ifdef AGGLO_DEBUG + const auto &points0 = fe_faces0.get_quadrature_points(); + const auto &points1 = fe_faces1.get_quadrature_points(); + for (size_t i = 0; + i < fe_faces1.get_quadrature_points().size(); + ++i) + { + double d = (points0[i] - points1[i]).norm(); + Assert(d < 1e-15, + ExcMessage( + "Face qpoints at the interface do not match!")); + } + +#endif + + std::vector + local_dof_indices_neighbor(dofs_per_cell); + + M11 = 0.; + M12 = 0.; + M21 = 0.; + M22 = 0.; + + const auto &normals = fe_faces0.get_normal_vectors(); + const double penalty = + penalty_constant / std::fabs(polytope->diameter()); + + // M11 + for (unsigned int q_index : + fe_faces0.quadrature_point_indices()) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + M11(i, j) += + (-0.5 * fe_faces0.shape_grad(i, q_index) * + normals[q_index] * + fe_faces0.shape_value(j, q_index) - + 0.5 * fe_faces0.shape_grad(j, q_index) * + normals[q_index] * + fe_faces0.shape_value(i, q_index) + + (penalty)*fe_faces0.shape_value(i, q_index) * + fe_faces0.shape_value(j, q_index)) * + fe_faces0.JxW(q_index); + + M12(i, j) += + (0.5 * fe_faces0.shape_grad(i, q_index) * + normals[q_index] * + fe_faces1.shape_value(j, q_index) - + 0.5 * fe_faces1.shape_grad(j, q_index) * + normals[q_index] * + fe_faces0.shape_value(i, q_index) - + (penalty)*fe_faces0.shape_value(i, q_index) * + fe_faces1.shape_value(j, q_index)) * + fe_faces1.JxW(q_index); + + // A10 + M21(i, j) += + (-0.5 * fe_faces1.shape_grad(i, q_index) * + normals[q_index] * + fe_faces0.shape_value(j, q_index) + + 0.5 * fe_faces0.shape_grad(j, q_index) * + normals[q_index] * + fe_faces1.shape_value(i, q_index) - + (penalty)*fe_faces1.shape_value(i, q_index) * + fe_faces0.shape_value(j, q_index)) * + fe_faces1.JxW(q_index); + + // A11 + M22(i, j) += + (0.5 * fe_faces1.shape_grad(i, q_index) * + normals[q_index] * + fe_faces1.shape_value(j, q_index) + + 0.5 * fe_faces1.shape_grad(j, q_index) * + normals[q_index] * + fe_faces1.shape_value(i, q_index) + + (penalty)*fe_faces1.shape_value(i, q_index) * + fe_faces1.shape_value(j, q_index)) * + fe_faces1.JxW(q_index); + } + } + } + + // distribute DoFs accordingly + // Retrieve DoFs info from the cell iterator. + neigh_polytope->get_dof_indices(local_dof_indices_neighbor); + + constraints.distribute_local_to_global(M11, + local_dof_indices, + system_matrix); + constraints.distribute_local_to_global( + M12, + local_dof_indices, + local_dof_indices_neighbor, + system_matrix); + constraints.distribute_local_to_global( + M21, + local_dof_indices_neighbor, + local_dof_indices, + system_matrix); + constraints.distribute_local_to_global( + M22, local_dof_indices_neighbor, system_matrix); + } // Loop only once trough internal faces + } + } // Loop over faces of current cell + } // Loop over cells + + AssertThrow( + std::fabs(test_integral - 1.) < TOL, + ExcMessage( + "Value for integral of linear function on this domain is not correct.")); + AssertThrow(std::fabs(test_volume - 1.) < TOL, + ExcMessage("Value of measure of domain is not correct.")); + AssertThrow(std::fabs(test_bdary - 4.) < TOL, + ExcMessage( + "Value for the measure of the boundary is not correct.")); +} + + + +template +void +Poisson::solve() +{ + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + A_direct.vmult(solution, system_rhs); +} + + + +template +void +Poisson::output_results() +{ + // Compute errors. + + // Prepare interpolation matrix onto the finer grid. + Vector interpolated_solution; + PolyUtils::interpolate_to_fine_grid(*ah, interpolated_solution, solution); + + // L2 error + Vector difference_per_cell(tria.n_active_cells()); + + VectorTools::integrate_difference(mapping, + ah->output_dh, + interpolated_solution, + *analytical_solution, + difference_per_cell, + QGauss(dg_fe->degree + 1), + VectorTools::L2_norm); + + const double L2_error = + VectorTools::compute_global_error(tria, + difference_per_cell, + VectorTools::L2_norm); + + // std::cout << "L2 error:" << L2_error << std::endl; + AssertThrow(L2_error < TOL, ExcMessage("L2 error too large.")); + + + + // H1 seminorm + Vector difference_per_cell_H1_semi(tria.n_active_cells()); + + VectorTools::integrate_difference(mapping, + ah->output_dh, + interpolated_solution, + *analytical_solution, + difference_per_cell_H1_semi, + QGauss(dg_fe->degree + 1), + VectorTools::H1_seminorm); + + const double H1_seminorm = + VectorTools::compute_global_error(tria, + difference_per_cell_H1_semi, + VectorTools::H1_seminorm); + AssertThrow(H1_seminorm < TOL, ExcMessage("H1 seminorm too large.")); + + // std::cout << "H1 seminorm:" << H1_seminorm << std::endl; +} + + +template +void +Poisson::run() +{ + make_grid(); + setup_agglomeration(); + assemble_system(); + solve(); + output_results(); +} + + + +int +main() +{ + deallog.depth_console(1); + + try + { + Poisson<2> poisson_problem_linear_sol{SolutionType::LinearSolution}; + poisson_problem_linear_sol.run(); + std::cout << "Linear: OK" << std::endl; + + Poisson<2> poisson_problem_quadratic_sol{SolutionType::QuadraticSolution}; + poisson_problem_quadratic_sol.run(); + std::cout << "Quadratic: OK" << std::endl; + } + catch (const std::exception &exc) + { + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + } + + return 0; +} diff --git a/test/polydeal/exact_solutions_dgp.output b/test/polydeal/exact_solutions_dgp.output new file mode 100644 index 00000000..a25026aa --- /dev/null +++ b/test/polydeal/exact_solutions_dgp.output @@ -0,0 +1,4 @@ +Size of tria: 16 +Linear: OK +Size of tria: 16 +Quadratic: OK