From 531222e76a1eeb32e5e4f8ac5b432f167450b78f Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 19 Jun 2022 09:50:04 +0200 Subject: [PATCH] Add test (MG+DG+hp) --- tests/matrix_free/poisson_dg_hp.cc | 363 ++++++++++++++++++ ...sson_dg_hp.with_p4est=true.mpirun=1.output | 3 + ...sson_dg_hp.with_p4est=true.mpirun=2.output | 3 + 3 files changed, 369 insertions(+) create mode 100644 tests/matrix_free/poisson_dg_hp.cc create mode 100644 tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=1.output create mode 100644 tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=2.output diff --git a/tests/matrix_free/poisson_dg_hp.cc b/tests/matrix_free/poisson_dg_hp.cc new file mode 100644 index 0000000000..7907f02386 --- /dev/null +++ b/tests/matrix_free/poisson_dg_hp.cc @@ -0,0 +1,363 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Solve Poisson problem problem on a with DG, MatrixFree and hp. + +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + + +double +get_penalty_parameter(const unsigned int i, + const unsigned int j, + const unsigned int degree) +{ + if (degree == 1) + { + if (i != j) + return 32.0; + if (i == 0) + return 32.0; + if (i == 1) + return 64.0; + } + else if (degree == 2) + { + if (i != j) + return 32.0; + if (i == 0) + return 32.0; + if (i == 1) + return 64.0; + } + + Assert(false, ExcNotImplemented()); + + return 0.0; +} + + + +template +class PoissonOperator +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + using number = double; + + using FECellIntegrator = FEEvaluation; + using FEFaceIntegrator = FEFaceEvaluation; + + PoissonOperator(const MatrixFree &matrix_free, + const unsigned int degree) + : matrix_free(matrix_free) + , degree(degree) + {} + + void + initialize_dof_vector(VectorType &vec) + { + matrix_free.initialize_dof_vector(vec); + } + + void + rhs(VectorType &vec) const + { + const int dummy = 0; + + matrix_free.template cell_loop( + [&](const auto &data, auto &dst, const auto &, const auto range) { + FECellIntegrator phi(matrix_free, range); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_value(1.0, q); + + phi.integrate_scatter(true, false, dst); + } + }, + vec, + dummy, + true); + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + matrix_free.template loop( + [&](const auto &data, auto &dst, const auto &src, const auto range) { + FECellIntegrator phi(matrix_free, range); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(src, EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate_scatter(EvaluationFlags::gradients, dst); + } + }, + [&](const auto &data, auto &dst, const auto &src, const auto range) { + FEFaceIntegrator fe_eval(data, range, true); + FEFaceIntegrator fe_eval_neighbor(data, range, false); + + for (unsigned int face = range.first; face < range.second; ++face) + { + fe_eval.reinit(face); + fe_eval_neighbor.reinit(face); + + fe_eval.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + fe_eval_neighbor.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArray sigmaF = + get_penalty_parameter(data.get_face_active_fe_index(range, true), + data.get_face_active_fe_index(range, false), + degree); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + VectorizedArray average_value = + (fe_eval.get_value(q) - fe_eval_neighbor.get_value(q)) * 0.5; + VectorizedArray average_valgrad = + fe_eval.get_normal_derivative(q) + + fe_eval_neighbor.get_normal_derivative(q); + average_valgrad = + average_value * 2. * sigmaF - average_valgrad * 0.5; + fe_eval.submit_normal_derivative(-average_value, q); + fe_eval_neighbor.submit_normal_derivative(-average_value, q); + fe_eval.submit_value(average_valgrad, q); + fe_eval_neighbor.submit_value(-average_valgrad, q); + } + fe_eval.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + fe_eval_neighbor.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }, + [&](const auto &data, auto &dst, const auto &src, const auto range) { + FEFaceIntegrator fe_eval(data, range, true); + + for (unsigned int face = range.first; face < range.second; ++face) + { + fe_eval.reinit(face); + fe_eval.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArray sigmaF = + get_penalty_parameter(data.get_face_active_fe_index(range), + data.get_face_active_fe_index(range), + degree); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + VectorizedArray average_value = fe_eval.get_value(q); + VectorizedArray average_valgrad = + -fe_eval.get_normal_derivative(q); + average_valgrad += average_value * sigmaF; + fe_eval.submit_normal_derivative(-average_value, q); + fe_eval.submit_value(average_valgrad, q); + } + + fe_eval.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }, + dst, + src, + true); + } + +private: + const MatrixFree &matrix_free; + const unsigned int degree; +}; + +template +void +test(const unsigned int degree) +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + + unsigned int subdivisions = degree == 1 ? 16 : 8; + + GridGenerator::subdivided_hyper_cube(tria, subdivisions); + + FE_DGQ fe1(degree); + FE_DGQ fe2(degree + 1); + hp::FECollection fes(fe1, fe2); + + QGauss quad(degree + 2); + MappingQ mapping(1); + + DoFHandler dof_handler(tria); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if (cell->center()[0] < 0.5) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); + } + + dof_handler.distribute_dofs(fes); + + AffineConstraints constraints; + constraints.close(); + + const auto solve_and_postprocess = + [&](const auto &poisson_operator, + auto & x, + auto & b) -> std::pair { + ReductionControl reduction_control(2000, 1e-7, 1e-2); + SolverCG::type> solver( + reduction_control); + + solver.solve(poisson_operator, x, b, PreconditionIdentity()); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + printf("Solved in %d iterations.\n", reduction_control.last_step()); + + constraints.distribute(x); + +#if 1 + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + x.update_ghost_values(); + data_out.add_data_vector(dof_handler, x, "solution"); + data_out.build_patches(mapping, degree + 1); + data_out.write_vtu_with_pvtu_record("./", + "result-" + std::to_string(dim) + "-" + + std::to_string(degree), + 0, + MPI_COMM_WORLD); +#endif + + Vector difference(tria.n_active_cells()); + + deallog << "dim=" << dim << ' '; + deallog << "degree=" << degree << ' '; + + VectorTools::integrate_difference(mapping, + dof_handler, + x, + Functions::ZeroFunction(), + difference, + quad, + VectorTools::NormType::L2_norm); + + const double error = + VectorTools::compute_global_error(tria, + difference, + VectorTools::NormType::L2_norm); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + printf("Error %f.\n", error); + + if (error < 0.042) + deallog << "OK" << std::endl; + else + deallog << "FAIL" << std::endl; + + return {reduction_control.last_step(), reduction_control.last_value()}; + }; + + const auto mf_algo = [&]() { + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags = update_gradients | update_values; + additional_data.mapping_update_flags_inner_faces = + update_gradients | update_values; + additional_data.mapping_update_flags_boundary_faces = + update_gradients | update_values; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, quad, additional_data); + + PoissonOperator poisson_operator(matrix_free, degree); + + LinearAlgebra::distributed::Vector x, b; + poisson_operator.initialize_dof_vector(x); + poisson_operator.initialize_dof_vector(b); + + poisson_operator.rhs(b); + + return solve_and_postprocess(poisson_operator, x, b); + }; + + mf_algo(); +} + + +int +main(int argc, char **argv) +{ + initlog(); + + deallog.depth_file(1); + + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + + test<2>(/*degree=*/1); + test<2>(/*degree=*/2); +} diff --git a/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=1.output b/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..9bb3640efd --- /dev/null +++ b/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=1.output @@ -0,0 +1,3 @@ + +DEAL::dim=2 degree=1 OK +DEAL::dim=2 degree=2 OK diff --git a/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=2.output b/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..9bb3640efd --- /dev/null +++ b/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=2.output @@ -0,0 +1,3 @@ + +DEAL::dim=2 degree=1 OK +DEAL::dim=2 degree=2 OK -- 2.39.5