From: Peter Munch Date: Sat, 28 Nov 2020 19:16:07 +0000 (+0100) Subject: Add Poisson test: MatrixFree + DG + simplex mesh X-Git-Tag: v9.3.0-rc1~827^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11235%2Fhead;p=dealii.git Add Poisson test: MatrixFree + DG + simplex mesh --- diff --git a/tests/simplex/matrix_free_03.cc b/tests/simplex/matrix_free_03.cc new file mode 100644 index 0000000000..5eee6625c6 --- /dev/null +++ b/tests/simplex/matrix_free_03.cc @@ -0,0 +1,611 @@ +// --------------------------------------------------------------------- +// +// 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 and Helmholtz problem on a simplex mesh with +// continuous elements and compare results between matrix-free and matrix-based +// implementations. + +#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; + +const double PENALTY = 8; + + +template +class SmoothSolution : public Function +{ +public: + SmoothSolution() + : Function() + {} + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const override; +}; + +template +void +SmoothSolution::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] = 0.0; +} + +template +class SmoothRightHandSide : public Function +{ +public: + SmoothRightHandSide() + : Function() + {} + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int /*component*/ = 0) const override; +}; + +template +void +SmoothRightHandSide::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] = 1.0; +} + + +template +class PoissonOperator +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + using number = double; + + PoissonOperator(const MatrixFree &matrix_free) + : matrix_free(matrix_free) + {} + + 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 cells) { + FEEvaluation phi(data); + for (unsigned int cell = cells.first; cell < cells.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); + } + + const int fe_degree = 5; /*TODO*/ + + + void + vmult(VectorType &dst, const VectorType &src) const + { + matrix_free.template loop( + [&](const auto &data, auto &dst, const auto &src, const auto cell_range) { + FEEvaluation phi(data); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; + ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(EvaluationFlags::gradients); + phi.set_dof_values(dst); + } + }, + [&](const auto &data, auto &dst, const auto &src, const auto face_range) { + FEFaceEvaluation fe_eval(data, true); + FEFaceEvaluation fe_eval_neighbor(data, false); + + for (unsigned int face = face_range.first; face < 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 = PENALTY; + // (std::abs((fe_eval.get_normal_vector(0) * + // fe_eval.inverse_jacobian(0))[dim - 1]) + + // std::abs((fe_eval.get_normal_vector(0) * + // fe_eval_neighbor.inverse_jacobian(0))[dim - 1])) * + // (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); + + 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 face_range) { + FEFaceEvaluation fe_eval(data, true); + for (unsigned int face = face_range.first; face < face_range.second; + face++) + { + fe_eval.reinit(face); + fe_eval.read_dof_values(src); + fe_eval.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArray sigmaF = PENALTY; + // std::abs((fe_eval.get_normal_vector(0) * + // fe_eval.inverse_jacobian(0))[dim - 1]) * + // number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; + + 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); + } + +private: + const MatrixFree &matrix_free; +}; + + +struct CopyDataFace +{ + FullMatrix cell_matrix; + std::vector joint_dof_indices; + std::array values; + std::array cell_indices; +}; + +struct CopyData +{ + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + double value; + unsigned int cell_index; + template + void + reinit(const Iterator &cell, unsigned int dofs_per_cell) + { + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + } +}; + +template +void +test(const unsigned int degree) +{ + Triangulation tria; + +#if true + unsigned int n_subdivisions = 0; + + if (dim == 2 && degree == 1) + n_subdivisions = 16; + + if (dim == 2 && degree == 2) + n_subdivisions = 8; + + if (dim == 3 && degree == 1) + n_subdivisions = 8; + + if (dim == 3 && degree == 2) + n_subdivisions = 4; + + GridGenerator::subdivided_hyper_cube_with_simplices(tria, n_subdivisions); + + Simplex::FE_DGP fe(degree); + Simplex::QGauss quadrature(degree + 1); + Simplex::QGauss face_quadrature(degree + 1); + MappingFE mapping(Simplex::FE_P(1)); +#else + GridGenerator::subdivided_hyper_cube(tria, dim == 2 ? 16 : 8); + + FE_DGQ fe(degree); + QGauss quadrature(degree + 1); + MappingFE mapping(FE_Q(1)); + +#endif + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + AffineConstraints constraints; + constraints.close(); + + const auto solve_and_postprocess = + [&](const auto &poisson_operator, + auto & x, + auto & b) -> std::pair { + ReductionControl reduction_control(1000, 1e-7, 1e-3); + SolverCG::type> solver( + reduction_control); + + try + { + solver.solve(poisson_operator, x, b, PreconditionIdentity()); + } + catch (const std::exception &e) + { + deallog << e.what() << std::endl; + } + + 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, 2); + data_out.write_vtu_with_pvtu_record("./", "result", 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, + quadrature, + VectorTools::NormType::L2_norm); + + deallog << VectorTools::compute_global_error(tria, + difference, + VectorTools::NormType::L2_norm) + << 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, quadrature, additional_data); + + PoissonOperator poisson_operator(matrix_free); + + 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); + }; + + const auto mb_algo = [&]() { + Vector solution, system_rhs; + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + SparsityPattern sparsity_pattern; + DoFTools::make_flux_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + + SparseMatrix system_matrix; + system_matrix.reinit(sparsity_pattern); + + const double diffusion_coefficient = 1.0; + + const auto exact_solution = std::make_shared>(); + const auto rhs_function = std::make_shared>(); + + const auto cell_worker = + [&](const auto &cell, auto &scratch_data, auto ©_data) { + const FEValues &fe_v = scratch_data.reinit(cell); + const unsigned int dofs_per_cell = fe_v.dofs_per_cell; + copy_data.reinit(cell, dofs_per_cell); + + const auto & q_points = scratch_data.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const std::vector &JxW = scratch_data.get_JxW_values(); + + std::vector rhs(n_q_points); + rhs_function->value_list(q_points, rhs); + + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + diffusion_coefficient * // nu + fe_v.shape_grad(i, point) * // grad v_h + fe_v.shape_grad(j, point) * // grad u_h + JxW[point]; // dx + + copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h + rhs[point] * // f + JxW[point]; // dx + } + }; + + const auto boundary_worker = [&](const auto & cell, + const unsigned int &face_no, + auto & scratch_data, + auto & copy_data) { + const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); + + const auto & q_points = scratch_data.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const unsigned int dofs_per_cell = fe_fv.dofs_per_cell; + + const std::vector & JxW = scratch_data.get_JxW_values(); + const std::vector> &normals = + scratch_data.get_normal_vectors(); + + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); + + const double penalty = PENALTY; + + for (unsigned int point = 0; point < n_q_points; ++point) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + (-diffusion_coefficient * // - nu + fe_fv.shape_value(i, point) * // v_h + (fe_fv.shape_grad(j, point) * // (grad u_h . + normals[point]) // n) + + - diffusion_coefficient * // - nu + (fe_fv.shape_grad(i, point) * // (grad v_h . + normals[point]) * // n) + fe_fv.shape_value(j, point) // u_h + + + diffusion_coefficient * penalty * // + nu sigma + fe_fv.shape_value(i, point) * // v_h + fe_fv.shape_value(j, point) // u_h + + ) * + JxW[point]; // dx + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + copy_data.cell_rhs(i) += + (-diffusion_coefficient * // - nu + (fe_fv.shape_grad(i, point) * // (grad v_h . + normals[point]) * // n) + g[point] // g + + + + diffusion_coefficient * penalty * // + nu sigma + fe_fv.shape_value(i, point) * g[point] // v_h g + + ) * + JxW[point]; // dx + } + }; + + const auto face_worker = [&](const auto & cell, + const unsigned int &f, + const unsigned int &sf, + const auto & ncell, + const unsigned int &nf, + const unsigned int &nsf, + auto & scratch_data, + auto & copy_data) { + const FEInterfaceValues &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + + const auto & q_points = fe_iv.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + + copy_data.face_data.emplace_back(); + CopyDataFace & copy_data_face = copy_data.face_data.back(); + const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs(); + copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices(); + copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face); + + const std::vector & JxW = fe_iv.get_JxW_values(); + const std::vector> &normals = fe_iv.get_normal_vectors(); + + const double penalty = PENALTY; + + for (unsigned int point = 0; point < n_q_points; ++point) + { + for (unsigned int i = 0; i < n_dofs_face; ++i) + for (unsigned int j = 0; j < n_dofs_face; ++j) + copy_data_face.cell_matrix(i, j) += + (-diffusion_coefficient * // - nu + fe_iv.jump(i, point) * // [v_h] + (fe_iv.average_gradient(j, point) * // ({grad u_h} . + normals[point]) // n) + + - diffusion_coefficient * // - nu + (fe_iv.average_gradient(i, point) * // (grad v_h . + normals[point]) * // n) + fe_iv.jump(j, point) // [u_h] + + + diffusion_coefficient * penalty * // + nu sigma + fe_iv.jump(i, point) * // [v_h] + fe_iv.jump(j, point) // [u_h] + + ) * + JxW[point]; // dx + } + }; + + AffineConstraints constraints; + constraints.close(); + const auto copier = [&](const auto &c) { + constraints.distribute_local_to_global(c.cell_matrix, + c.cell_rhs, + c.local_dof_indices, + system_matrix, + system_rhs); + + for (auto &cdf : c.face_data) + { + constraints.distribute_local_to_global(cdf.cell_matrix, + cdf.joint_dof_indices, + system_matrix); + } + }; + + UpdateFlags cell_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + UpdateFlags face_flags = update_values | update_gradients | + update_quadrature_points | update_normal_vectors | + update_JxW_values; + + MeshWorker::ScratchData scratch_data( + mapping, fe, quadrature, cell_flags, face_quadrature, face_flags); + CopyData cd; + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + cd, + MeshWorker::assemble_own_cells | + MeshWorker::assemble_boundary_faces | + MeshWorker::assemble_own_interior_faces_once, + boundary_worker, + face_worker); + + return solve_and_postprocess(system_matrix, solution, system_rhs); + }; + + mb_algo(); + 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); + test<3>(/*degree=*/1); + test<3>(/*degree=*/2); +} diff --git a/tests/simplex/matrix_free_03.with_simplex_support=on.output b/tests/simplex/matrix_free_03.with_simplex_support=on.output new file mode 100644 index 0000000000..730828db40 --- /dev/null +++ b/tests/simplex/matrix_free_03.with_simplex_support=on.output @@ -0,0 +1,9 @@ + +DEAL::dim=2 degree=1 0.0412453 +DEAL::dim=2 degree=1 0.0412453 +DEAL::dim=2 degree=2 0.0412685 +DEAL::dim=2 degree=2 0.0412685 +DEAL::dim=3 degree=1 0.0258555 +DEAL::dim=3 degree=1 0.0258555 +DEAL::dim=3 degree=2 0.0250053 +DEAL::dim=3 degree=2 0.0250053