From: Peter Munch Date: Wed, 2 Dec 2020 21:10:03 +0000 (+0100) Subject: Add Poisson test: MatrixFree + DG + mixed mesh X-Git-Tag: v9.3.0-rc1~801^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11273%2Fhead;p=dealii.git Add Poisson test: MatrixFree + DG + mixed mesh --- diff --git a/tests/simplex/matrix_free_04.cc b/tests/simplex/matrix_free_04.cc new file mode 100644 index 0000000000..6d46b58889 --- /dev/null +++ b/tests/simplex/matrix_free_04.cc @@ -0,0 +1,372 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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 mixed mesh with DG and MatrixFree. + +#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" + +#include "./simplex_grids.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 8.0; + if (i == 0) + return 8.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 version, const unsigned int degree) +{ + Triangulation tria; + + unsigned int subdivisions = degree == 1 ? 16 : 8; + + if (version == 0) + GridGenerator::subdivided_hyper_cube_with_simplices(tria, subdivisions); + else if (version == 1) + GridGenerator::subdivided_hyper_cube(tria, subdivisions); + else if (version == 2) + GridGenerator::subdivided_hyper_cube_with_simplices_mix(tria, subdivisions); + + Simplex::FE_DGP fe1(degree); + FE_DGQ fe2(degree); + hp::FECollection fes(fe1, fe2); + + Simplex::QGauss quad1(degree + 1); + QGauss quad2(degree + 1); + hp::QCollection quads(quad1, quad2); + + MappingFE mapping1(Simplex::FE_P(1)); + MappingQ mapping2(1); + hp::MappingCollection mappings(mapping1, mapping2); + + DoFHandler dof_handler(tria); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->reference_cell_type() == ReferenceCell::Type::Tri || + cell->reference_cell_type() == ReferenceCell::Type::Tet) + 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 0 + 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(mappings, degree); + data_out.write_vtu_with_pvtu_record("./", + "result-" + std::to_string(dim) + "-" + + std::to_string(degree) + "-" + + std::to_string(version), + 0, + MPI_COMM_WORLD); +#endif + + Vector difference(tria.n_active_cells()); + + deallog << "dim=" << dim << " "; + deallog << "degree=" << degree << " "; + + VectorTools::integrate_difference(mappings, + dof_handler, + x, + Functions::ZeroFunction(), + difference, + quads, + VectorTools::NormType::L2_norm); + + const double error = + VectorTools::compute_global_error(tria, + difference, + VectorTools::NormType::L2_norm); + + 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( + mappings, dof_handler, constraints, quads, 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); + + for (unsigned int i = 0; i < 3; ++i) + test<2>(i, /*degree=*/1); + + for (unsigned int i = 0; i < 3; ++i) + test<2>(i, /*degree=*/2); +} diff --git a/tests/simplex/matrix_free_04.with_simplex_support=on.output b/tests/simplex/matrix_free_04.with_simplex_support=on.output new file mode 100644 index 0000000000..6e86258005 --- /dev/null +++ b/tests/simplex/matrix_free_04.with_simplex_support=on.output @@ -0,0 +1,7 @@ + +DEAL::dim=2 degree=1 OK +DEAL::dim=2 degree=1 OK +DEAL::dim=2 degree=1 OK +DEAL::dim=2 degree=2 OK +DEAL::dim=2 degree=2 OK +DEAL::dim=2 degree=2 OK