From: Martin Kronbichler Date: Wed, 4 Jul 2018 13:36:15 +0000 (+0200) Subject: Add test cases. X-Git-Tag: v9.1.0-rc1~880^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4124b79ead8d704e6817f6ae955a953d20f704f5;p=dealii.git Add test cases. --- diff --git a/tests/matrix_free/advect_1d.cc b/tests/matrix_free/advect_1d.cc new file mode 100644 index 0000000000..1fb04f7dcc --- /dev/null +++ b/tests/matrix_free/advect_1d.cc @@ -0,0 +1,327 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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. +// +// --------------------------------------------------------------------- + + + +// similar to matrix_vector_faces_15 and the code in matrix_vector_common but +// for 1D + +#include +#include + +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +std::ofstream logfile("output"); + + +template , + int n_components = 1> +class MatrixFreeAdvectionBasic +{ +public: + MatrixFreeAdvectionBasic(const MatrixFree &data, + const bool zero_within_loop = true, + const unsigned int start_vector_component = 0) + : data(data) + , zero_within_loop(zero_within_loop) + , start_vector_component(start_vector_component) + { + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + if (!zero_within_loop) + dst = 0; + data.loop(&MatrixFreeAdvectionBasic::local_apply, + &MatrixFreeAdvectionBasic::local_apply_face, + &MatrixFreeAdvectionBasic::local_apply_boundary_face, + this, + dst, + src, + zero_within_loop, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + + FEEvaluation phi(data); + + const unsigned int dofs_per_cell = phi.dofs_per_cell; + + AlignedVector> coefficients(phi.dofs_per_cell); + MatrixFreeOperators:: + CellwiseInverseMassMatrix + inverse(phi); + + for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell) + { + phi.reinit(cell); + phi.read_dof_values(dst); + + inverse.fill_inverse_JxW_values(coefficients); + inverse.apply(coefficients, + n_components, + phi.begin_dof_values(), + phi.begin_dof_values()); + + phi.set_dof_values(dst); + } + } + +private: + void + local_apply(const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation phi( + data, 0, 0, start_vector_component); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(advection * phi.get_value(q), q); + phi.integrate(false, true); + phi.distribute_local_to_global(dst); + } + } + + void + local_apply_face( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &face_range) const + { + FEFaceEvaluation phi_m( + data, true, 0, 0, start_vector_component); + FEFaceEvaluation phi_p( + data, false, 0, 0, start_vector_component); + typedef typename FEFaceEvaluation::value_type value_type; + + for (unsigned int face = face_range.first; face < face_range.second; face++) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + phi_m.evaluate(true, false); + phi_p.reinit(face); + phi_p.read_dof_values(src); + phi_p.evaluate(true, false); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + value_type u_minus = phi_m.get_value(q), + u_plus = phi_p.get_value(q); + const VectorizedArray normal_times_advection = + advection * phi_m.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + phi_m.submit_value(-flux_times_normal, q); + phi_p.submit_value(flux_times_normal, q); + } + + phi_m.integrate(true, false); + phi_m.distribute_local_to_global(dst); + phi_p.integrate(true, false); + phi_p.distribute_local_to_global(dst); + } + } + + void + local_apply_boundary_face( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &face_range) const + { + FEFaceEvaluation + fe_eval(data, true, 0, 0, start_vector_component); + typedef typename FEFaceEvaluation::value_type value_type; + + 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(true, false); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + value_type u_minus = fe_eval.get_value(q); + value_type u_plus = -u_minus; + const VectorizedArray normal_times_advection = + advection * fe_eval.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + fe_eval.submit_value(-flux_times_normal, q); + } + + fe_eval.integrate(true, false); + fe_eval.distribute_local_to_global(dst); + } + } + + const MatrixFree & data; + const bool zero_within_loop; + const unsigned int start_vector_component; + Tensor<1, dim, VectorizedArray> advection; +}; + + + +template +class AnalyticFunction : public Function +{ +public: + static_assert(dim == 1, "Only 1D implemented"); + AnalyticFunction() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int) const override + { + return std::sin(3 * numbers::PI * p[0] / 0.8); + } +}; + + + +template +class AnalyticDerivative : public Function +{ +public: + static_assert(dim == 1, "Only 1D implemented"); + AnalyticDerivative() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int) const override + { + Tensor<1, dim> advection; + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + + return -std::cos(3 * numbers::PI * p[0] / 0.8) * advection[0] * 3 * + numbers::PI / 0.8; + } +}; + + + +template +void +test(const unsigned int n_refine) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria, 0, 0.8); + tria.refine_global(n_refine); + + FE_DGQ fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + if (n_refine == 3) + { + deallog << "Testing " << dof.get_fe().get_name(); + deallog << std::endl; + } + + LinearAlgebra::distributed::Vector in, out; + + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + MatrixFree mf_data; + mf_data.reinit(dof, constraints, quad, data); + + mf_data.initialize_dof_vector(in); + mf_data.initialize_dof_vector(out); + + VectorTools::interpolate(dof, AnalyticFunction(), in); + + MatrixFreeAdvectionBasic> + mf2(mf_data); + mf2.vmult(out, in); + + VectorTools::interpolate(dof, AnalyticDerivative(), in); + out -= in; + + double diff_norm = out.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << " " << std::endl; +} + + + +int +main() +{ + initlog(); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 2>(r); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 4>(r); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 5>(r); +} diff --git a/tests/matrix_free/advect_1d.output b/tests/matrix_free/advect_1d.output new file mode 100644 index 0000000000..b213520047 --- /dev/null +++ b/tests/matrix_free/advect_1d.output @@ -0,0 +1,22 @@ + +DEAL::Testing FE_DGQ<1>(2) +DEAL::Norm of difference: 0.531539 +DEAL::Norm of difference: 0.135407 +DEAL::Norm of difference: 0.0340112 +DEAL::Norm of difference: 0.00851279 +DEAL::Norm of difference: 0.00212882 +DEAL::Norm of difference: 0.000532245 +DEAL::Testing FE_DGQ<1>(4) +DEAL::Norm of difference: 0.00529591 +DEAL::Norm of difference: 0.000336014 +DEAL::Norm of difference: 2.10801e-05 +DEAL::Norm of difference: 1.31874e-06 +DEAL::Norm of difference: 8.24413e-08 +DEAL::Norm of difference: 5.15331e-09 +DEAL::Testing FE_DGQ<1>(5) +DEAL::Norm of difference: 0.000346887 +DEAL::Norm of difference: 1.09982e-05 +DEAL::Norm of difference: 3.44938e-07 +DEAL::Norm of difference: 1.07881e-08 +DEAL::Norm of difference: 3.35085e-10 +DEAL::Norm of difference: 9.72689e-12 diff --git a/tests/matrix_free/advect_1d_system.cc b/tests/matrix_free/advect_1d_system.cc new file mode 100644 index 0000000000..0e89637ba7 --- /dev/null +++ b/tests/matrix_free/advect_1d_system.cc @@ -0,0 +1,336 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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. +// +// --------------------------------------------------------------------- + + + +// similar to advect_1d but for a system of 3 components + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +std::ofstream logfile("output"); + + +template , + int n_components = 1> +class MatrixFreeAdvectionBasic +{ +public: + MatrixFreeAdvectionBasic(const MatrixFree &data, + const bool zero_within_loop = true, + const unsigned int start_vector_component = 0) + : data(data) + , zero_within_loop(zero_within_loop) + , start_vector_component(start_vector_component) + { + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + if (!zero_within_loop) + dst = 0; + data.loop(&MatrixFreeAdvectionBasic::local_apply, + &MatrixFreeAdvectionBasic::local_apply_face, + &MatrixFreeAdvectionBasic::local_apply_boundary_face, + this, + dst, + src, + zero_within_loop, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + + FEEvaluation phi(data); + + const unsigned int dofs_per_cell = phi.dofs_per_cell; + + AlignedVector> coefficients(phi.dofs_per_cell); + MatrixFreeOperators:: + CellwiseInverseMassMatrix + inverse(phi); + + for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell) + { + phi.reinit(cell); + phi.read_dof_values(dst); + + inverse.fill_inverse_JxW_values(coefficients); + inverse.apply(coefficients, + n_components, + phi.begin_dof_values(), + phi.begin_dof_values()); + + phi.set_dof_values(dst); + } + } + +private: + void + local_apply(const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation phi( + data, 0, 0, start_vector_component); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + Tensor<1, n_components, Tensor<1, dim, VectorizedArray>> + grad; + for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int d = 0; d < dim; ++d) + grad[c][d] = advection[d] * phi.get_value(q)[c]; + phi.submit_gradient(grad, q); + } + phi.integrate(false, true); + phi.distribute_local_to_global(dst); + } + } + + void + local_apply_face( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &face_range) const + { + FEFaceEvaluation phi_m( + data, true, 0, 0, start_vector_component); + FEFaceEvaluation phi_p( + data, false, 0, 0, start_vector_component); + typedef typename FEFaceEvaluation::value_type value_type; + + for (unsigned int face = face_range.first; face < face_range.second; face++) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + phi_m.evaluate(true, false); + phi_p.reinit(face); + phi_p.read_dof_values(src); + phi_p.evaluate(true, false); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + value_type u_minus = phi_m.get_value(q), + u_plus = phi_p.get_value(q); + const VectorizedArray normal_times_advection = + advection * phi_m.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + phi_m.submit_value(-flux_times_normal, q); + phi_p.submit_value(flux_times_normal, q); + } + + phi_m.integrate(true, false); + phi_m.distribute_local_to_global(dst); + phi_p.integrate(true, false); + phi_p.distribute_local_to_global(dst); + } + } + + void + local_apply_boundary_face( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &face_range) const + { + FEFaceEvaluation + fe_eval(data, true, 0, 0, start_vector_component); + typedef typename FEFaceEvaluation::value_type value_type; + + 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(true, false); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + value_type u_minus = fe_eval.get_value(q); + value_type u_plus = -u_minus; + const VectorizedArray normal_times_advection = + advection * fe_eval.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + fe_eval.submit_value(-flux_times_normal, q); + } + + fe_eval.integrate(true, false); + fe_eval.distribute_local_to_global(dst); + } + } + + const MatrixFree & data; + const bool zero_within_loop; + const unsigned int start_vector_component; + Tensor<1, dim, VectorizedArray> advection; +}; + + + +template +class AnalyticFunction : public Function +{ +public: + static_assert(dim == 1, "Only 1D implemented"); + AnalyticFunction() + : Function(3) + {} + + virtual double + value(const Point &p, const unsigned int component) const override + { + return std::sin(3 * numbers::PI * p[0] / 0.8) * (component + 1.2); + } +}; + + + +template +class AnalyticDerivative : public Function +{ +public: + static_assert(dim == 1, "Only 1D implemented"); + AnalyticDerivative() + : Function(3) + {} + + virtual double + value(const Point &p, const unsigned int component) const override + { + Tensor<1, dim> advection; + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + + return -std::cos(3 * numbers::PI * p[0] / 0.8) * advection[0] * 3 * + numbers::PI / 0.8 * (component + 1.2); + } +}; + + + +template +void +test(const unsigned int n_refine) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria, 0, 0.8); + tria.refine_global(n_refine); + + FE_DGQ fe_dgq(fe_degree); + FESystem fe(fe_dgq, 3); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + if (n_refine == 3) + { + deallog << "Testing " << dof.get_fe().get_name(); + deallog << std::endl; + } + + LinearAlgebra::distributed::Vector in, out; + + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + MatrixFree mf_data; + mf_data.reinit(dof, constraints, quad, data); + + mf_data.initialize_dof_vector(in); + mf_data.initialize_dof_vector(out); + + VectorTools::interpolate(dof, AnalyticFunction(), in); + + MatrixFreeAdvectionBasic, + 3> + mf2(mf_data); + mf2.vmult(out, in); + + VectorTools::interpolate(dof, AnalyticDerivative(), in); + out -= in; + + double diff_norm = out.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << " " << std::endl; +} + + + +int +main() +{ + initlog(); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 2>(r); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 4>(r); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 5>(r); +} diff --git a/tests/matrix_free/advect_1d_system.output b/tests/matrix_free/advect_1d_system.output new file mode 100644 index 0000000000..6f3e54286b --- /dev/null +++ b/tests/matrix_free/advect_1d_system.output @@ -0,0 +1,22 @@ + +DEAL::Testing FESystem<1>[FE_DGQ<1>(2)^3] +DEAL::Norm of difference: 1.70093 +DEAL::Norm of difference: 0.433302 +DEAL::Norm of difference: 0.108836 +DEAL::Norm of difference: 0.0272409 +DEAL::Norm of difference: 0.00681224 +DEAL::Norm of difference: 0.00170318 +DEAL::Testing FESystem<1>[FE_DGQ<1>(4)^3] +DEAL::Norm of difference: 0.0169469 +DEAL::Norm of difference: 0.00107525 +DEAL::Norm of difference: 6.74562e-05 +DEAL::Norm of difference: 4.21998e-06 +DEAL::Norm of difference: 2.63812e-07 +DEAL::Norm of difference: 1.64905e-08 +DEAL::Testing FESystem<1>[FE_DGQ<1>(5)^3] +DEAL::Norm of difference: 0.00111004 +DEAL::Norm of difference: 3.51943e-05 +DEAL::Norm of difference: 1.10380e-06 +DEAL::Norm of difference: 3.45206e-08 +DEAL::Norm of difference: 1.07406e-09 +DEAL::Norm of difference: 2.59508e-11