From: Martin Kronbichler Date: Tue, 26 Mar 2024 16:40:05 +0000 (+0100) Subject: New test case X-Git-Tag: v9.6.0-rc1~446^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=23971dd1060af4695bd81a27721041660659cdce;p=dealii.git New test case --- diff --git a/tests/matrix_free/matrix_vector_stokes_1d.cc b/tests/matrix_free/matrix_vector_stokes_1d.cc new file mode 100644 index 0000000000..de6274e94c --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_1d.cc @@ -0,0 +1,343 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// this function tests the correctness of the implementation of matrix free +// matrix-vector products with vector-valued evaluators and their use in 1d in +// particular using a variant of the stokes equations. + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + +#include "create_mesh.h" + + + +template +class MatrixFreeTest +{ +public: + using CellIterator = typename DoFHandler::active_cell_iterator; + using Number = double; + + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in){}; + + void + local_apply(const MatrixFree &data, + VectorType &dst, + const VectorType &src, + const std::pair &cell_range) const + { + using vector_t = VectorizedArray; + FEEvaluation velocity(data, + 0); + FEEvaluation pressure(data, 1); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + velocity.reinit(cell); + velocity.read_dof_values(src, 0); + velocity.evaluate(EvaluationFlags::gradients | EvaluationFlags::values); + pressure.reinit(cell); + pressure.read_dof_values(src, dim); + pressure.evaluate(EvaluationFlags::values); + + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + Tensor<2, dim, vector_t> grad_u; + if constexpr (dim == 1) + grad_u[0] = velocity.get_gradient(q); + else + grad_u = velocity.get_gradient(q); + + Tensor<1, dim, vector_t> val_u; + if constexpr (dim == 1) + val_u[0] = velocity.get_value(q); + else + val_u = velocity.get_value(q); + + vector_t pres = pressure.get_value(q); + vector_t div = -trace(grad_u); + pressure.submit_value(div, q); + + // subtract p * I + for (unsigned int d = 0; d < dim; ++d) + grad_u[d][d] -= pres; + + velocity.submit_gradient(grad_u, q); + velocity.submit_value(val_u, q); + } + + velocity.integrate(EvaluationFlags::gradients | + EvaluationFlags::values); + velocity.distribute_local_to_global(dst, 0); + pressure.integrate(EvaluationFlags::values); + pressure.distribute_local_to_global(dst, dim); + } + } + + + void + vmult(VectorType &dst, const VectorType &src) const + { + AssertDimension(dst.size(), dim + 1); + for (unsigned int d = 0; d < dim + 1; ++d) + dst[d] = 0; + data.cell_loop(&MatrixFreeTest::local_apply, + this, + dst, + src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + Triangulation triangulation; + if constexpr (dim == 1) + GridGenerator::hyper_cube(triangulation); + else + create_mesh(triangulation); + + if (fe_degree == 1) + triangulation.refine_global(4 - dim); + else + triangulation.refine_global(3 - dim); + + FE_Q fe_u(fe_degree + 1); + FE_Q fe_p(fe_degree); + FESystem fe(fe_u, dim, fe_p, 1); + DoFHandler dof_handler_u(triangulation); + DoFHandler dof_handler_p(triangulation); + DoFHandler dof_handler(triangulation); + + MatrixFree mf_data; + + AffineConstraints constraints; + constraints.close(); + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + std::vector> vec1, vec2; + + dof_handler.distribute_dofs(fe); + dof_handler_u.distribute_dofs(fe_u); + dof_handler_p.distribute_dofs(fe_p); + DoFRenumbering::component_wise(dof_handler); + + const std::vector dofs_per_block = + DoFTools::count_dofs_per_fe_component(dof_handler); + + { + BlockDynamicSparsityPattern csp(dim + 1, dim + 1); + + for (unsigned int d = 0; d < dim + 1; ++d) + for (unsigned int e = 0; e < dim + 1; ++e) + csp.block(d, e).reinit(dofs_per_block[d], dofs_per_block[e]); + + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false); + sparsity_pattern.copy_from(csp); + } + + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dim + 1); + for (unsigned int i = 0; i < dim + 1; ++i) + solution.block(i).reinit(dofs_per_block[i]); + solution.collect_sizes(); + + system_rhs.reinit(solution); + + vec1.resize(dim + 1); + vec2.resize(dim + 1); + vec1[0].reinit(dofs_per_block[0]); + vec2[0].reinit(vec1[0]); + for (unsigned int i = 1; i < dim; ++i) + { + vec1[i].reinit(vec1[0]); + vec2[i].reinit(vec1[0]); + } + vec1[dim].reinit(dofs_per_block[dim]); + vec2[dim].reinit(vec1[dim]); + + // this is from step-22 + { + QGauss quadrature_formula(fe_degree + 2); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + std::vector> phi_grad_u(dofs_per_cell); + std::vector> phi_u(dofs_per_cell); + std::vector div_phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + local_matrix = 0; + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + phi_grad_u[k] = fe_values[velocities].gradient(k, q); + phi_u[k] = fe_values[velocities].value(k, q); + div_phi_u[k] = fe_values[velocities].divergence(k, q); + phi_p[k] = fe_values[pressure].value(k, q); + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j <= i; ++j) + { + local_matrix(i, j) += + (phi_u[i] * phi_u[j] + + scalar_product(phi_grad_u[i], phi_grad_u[j]) - + div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) * + fe_values.JxW(q); + } + } + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + local_matrix(i, j) = local_matrix(j, i); + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(local_matrix, + local_dof_indices, + system_matrix); + } + } + + // first system_rhs with random numbers + for (unsigned int i = 0; i < dim + 1; ++i) + for (unsigned int j = 0; j < system_rhs.block(i).size(); ++j) + { + const double val = -1. + 2. * random_value(); + system_rhs.block(i)(j) = val; + vec1[i](j) = val; + } + + // setup matrix-free structure + { + std::vector *> dofs; + dofs.push_back(&dof_handler_u); + dofs.push_back(&dof_handler_p); + AffineConstraints dummy_constraints; + dummy_constraints.close(); + std::vector *> constraints; + constraints.push_back(&dummy_constraints); + constraints.push_back(&dummy_constraints); + QGauss<1> quad(fe_degree + 2); + mf_data.reinit(MappingQ1{}, + dofs, + constraints, + quad, + typename MatrixFree::AdditionalData( + MatrixFree::AdditionalData::none)); + } + + system_matrix.vmult(solution, system_rhs); + + using VectorType = std::vector>; + MatrixFreeTest mf(mf_data); + mf.vmult(vec2, vec1); + + // Verification + double error = 0.; + for (unsigned int i = 0; i < dim + 1; ++i) + for (unsigned int j = 0; j < system_rhs.block(i).size(); ++j) + error += std::fabs(solution.block(i)(j) - vec2[i](j)); + double relative = solution.block(0).l1_norm(); + deallog << " Verification fe degree " << fe_degree << ": " + << error / relative << std::endl + << std::endl; +} + + + +int +main() +{ + initlog(); + deallog << std::setprecision(3); + + { + deallog << std::endl << "Test with doubles" << std::endl << std::endl; + deallog.push("1d"); + test<1, 1>(); + deallog.pop(); + deallog.push("2d"); + test<2, 1>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_stokes_1d.output b/tests/matrix_free/matrix_vector_stokes_1d.output new file mode 100644 index 0000000000..37c9d6ea38 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_1d.output @@ -0,0 +1,10 @@ + +DEAL:: +DEAL::Test with doubles +DEAL:: +DEAL:1d:: Verification fe degree 1: 2.89e-16 +DEAL:1d:: +DEAL:2d:: Verification fe degree 1: 1.05e-15 +DEAL:2d:: +DEAL:3d:: Verification fe degree 1: 1.31e-15 +DEAL:3d::