From 20c95de33474c3d56eb31462ad98a0ae92d12c67 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 25 Feb 2020 19:51:46 +0100 Subject: [PATCH] Add test case --- tests/matrix_free/matrix_vector_faces_28.cc | 174 ++++++++++++++++++ .../matrix_free/matrix_vector_faces_28.output | 17 ++ 2 files changed, 191 insertions(+) create mode 100644 tests/matrix_free/matrix_vector_faces_28.cc create mode 100644 tests/matrix_free/matrix_vector_faces_28.output diff --git a/tests/matrix_free/matrix_vector_faces_28.cc b/tests/matrix_free/matrix_vector_faces_28.cc new file mode 100644 index 0000000000..8cfe4fb36d --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_28.cc @@ -0,0 +1,174 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// check that MatrixFree::update_mapping() works as expected. The test is +// otherwise similar to matrix_vector_faces_03. + +#include + +#include +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +#include "create_mesh.h" + +std::ofstream logfile("output"); + +#include "matrix_vector_faces_common.h" + + +template +void +test() +{ + Triangulation tria; + create_mesh(tria); + + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + typename Triangulation::active_cell_iterator cell, endc; + cell = tria.begin_active(); + endc = tria.end(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 0.5) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.begin(tria.n_levels() - 1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + // if (fe_degree == 1) + // tria.refine_global(1); + cell = tria.begin_active(); + for (unsigned int i = 0; i < 9 - 3 * dim; ++i) + { + cell = tria.begin_active(); + endc = tria.end(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + deallog << "Degree of element: " << fe_degree << std::endl; + + FE_DGQ fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + FESystem fe_grid(FE_Q(dof.get_fe().degree), dim); + DoFHandler dof_grid(tria); + dof_grid.distribute_dofs(fe_grid); + Vector euler; + euler.reinit(dof_grid.n_dofs()); + const ComponentMask mask(dim, true); + VectorTools::get_position_vector(dof_grid, euler, mask); + + MappingFEField mapping(dof_grid, euler, mask); + + Vector in(dof.n_dofs()), out(dof.n_dofs()); + Vector out_dist(out); + + // Set random seed for reproducibility + Testing::srand(42); + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + if (constraints.is_constrained(i)) + continue; + const double entry = Testing::rand() / (double)RAND_MAX; + in(i) = entry; + } + + auto create_reference_result = [&]() { + // assemble sparse matrix with MeshWorker + SparsityPattern sparsity; + SparseMatrix matrix; + { + DynamicSparsityPattern d_sparsity(dof.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dof, d_sparsity); + sparsity.copy_from(d_sparsity); + } + const unsigned int n_q_points_1d = dof.get_fe().degree + 1; + matrix.reinit(sparsity); + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = + update_values | update_gradients | update_jacobians; + info_box.add_update_flags_all(update_flags); + info_box.initialize_gauss_quadrature(n_q_points_1d, + n_q_points_1d, + n_q_points_1d); + info_box.initialize(dof.get_fe(), mapping); + + MeshWorker::DoFInfo dof_info(dof); + + MeshWorker::Assembler::MatrixSimple> assembler; + assembler.initialize(matrix); + + MatrixIntegrator integrator; + MeshWorker::integration_loop( + dof.begin_active(), dof.end(), dof_info, info_box, integrator, assembler); + + matrix.vmult(out, in); + }; + + MatrixFree mf_data; + const QGauss<1> quad(dof.get_fe().degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.tasks_block_size = 3; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + mf_data.reinit(mapping, dof, constraints, quad, data); + + MatrixFreeTest, 1> mf( + mf_data); + mf.vmult(out_dist, in); + create_reference_result(); + + out_dist -= out; + double diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << std::endl; + + // deform the geometry slightly + euler(0) += 0.001; + euler(euler.size() - 1) += 0.001; + + mf.vmult(out_dist, in); + create_reference_result(); + + out_dist -= out; + diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference no update: " << diff_norm << std::endl; + + mf_data.update_mapping(mapping); + mf.vmult(out_dist, in); + out_dist -= out; + diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference w update: " << diff_norm << std::endl; +} diff --git a/tests/matrix_free/matrix_vector_faces_28.output b/tests/matrix_free/matrix_vector_faces_28.output new file mode 100644 index 0000000000..11ef5b8048 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_28.output @@ -0,0 +1,17 @@ + +DEAL:2d::Degree of element: 1 +DEAL:2d::Norm of difference: 4.87e-16 +DEAL:2d::Norm of difference no update: 0.0179 +DEAL:2d::Norm of difference w update: 4.87e-16 +DEAL:2d::Degree of element: 2 +DEAL:2d::Norm of difference: 1.92e-15 +DEAL:2d::Norm of difference no update: 0.00672 +DEAL:2d::Norm of difference w update: 2.50e-15 +DEAL:3d::Degree of element: 1 +DEAL:3d::Norm of difference: 1.88e-16 +DEAL:3d::Norm of difference no update: 5.11e-05 +DEAL:3d::Norm of difference w update: 1.88e-16 +DEAL:3d::Degree of element: 2 +DEAL:3d::Norm of difference: 1.06e-16 +DEAL:3d::Norm of difference no update: 5.82e-05 +DEAL:3d::Norm of difference w update: 1.06e-16 -- 2.39.5