From 9f7a1d4e7741612578c457fba7517ce067603cdc Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 4 Oct 2018 23:35:59 -0600 Subject: [PATCH] Add a test for the parallel hp implementation. --- tests/mpi/hp_integrate_difference.cc | 254 ++++++++++++++++++ ...erence.with_trilinos=true.mpirun=10.output | 5 + ...ference.with_trilinos=true.mpirun=4.output | 5 + 3 files changed, 264 insertions(+) create mode 100644 tests/mpi/hp_integrate_difference.cc create mode 100644 tests/mpi/hp_integrate_difference.with_trilinos=true.mpirun=10.output create mode 100644 tests/mpi/hp_integrate_difference.with_trilinos=true.mpirun=4.output diff --git a/tests/mpi/hp_integrate_difference.cc b/tests/mpi/hp_integrate_difference.cc new file mode 100644 index 0000000000..acf45be4cf --- /dev/null +++ b/tests/mpi/hp_integrate_difference.cc @@ -0,0 +1,254 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2017 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. +// +// --------------------------------------------------------------------- + + + +// Test VectorTools::integrate_difference for parallel computations +// with the hp::DoFHandler. This includes applying hanging node +// constraints and consequently verifies that we compute them +// correctly. +// +// The way this test works is this: We create a domain [-1,1]^d +// and use the following function f(x,y) (or f(x,y,z)): +// - in the lower left quadrant, f(x,y)=xy +// - in the lower right quadrant, f(x,y)=xy+(xy)^2 +// - in the upper left quadrant, f(x,y)=xy+(xy)^3 +// - in the upper right quadrant, f(x,y)=xy+(xy)^2+(xy)^3+(xy)*4 +// +// We interpolate this function onto a finite element space that is +// chosen as follows: +// - in the lower left quadrant, Q1 +// - in the lower right quadrant, Q2 +// - in the upper left quadrant, Q3 +// - in the upper right quadrant, Q4 +// In other words, the function f(...) is in the space. +// +// We can then run two tests with it: +// - Interpolate the function onto the finite element space and compute +// its L2 norm. This can be done analytically. In particular, the +// area under the square of the functions above is, for the four +// quadrants: +// . 1/9 +// . 47/1800 +// . 2332/11025 +// . 1816349/3175200 +// This makes the total sum under the square of the function equal to +// 2923673/3175200 and the L2 norm under the function equal to +// sqrt(5847346)/2520, which is about 0.9595748472. +// (In 3d, we also integrate the same function over z=-1..1, so the +// volume under f(...)^2 increases by a factor of 2, and the L2 norm +// by a factor of sqrt(2); the numerical value is then 1.357043763.) +// - Interpolate the function onto the finite element space and compute +// the L2 norm of the difference between the interpolated function +// and the original function. This should be zero. +// This test does both. + +#include +#include + +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + + +template +class CheckFunction : public Function +{ +public: + double + value(const Point &p, const unsigned int) const + { + const double x = p[0]; + const double y = p[1]; + + // function is bilinear everywhere + double f = x * y; + + // on the right half of the domain, add something biquadratic that's + // zero at x=0 + if (x >= 0) + f += x * x * y * y; + + + // in the top half of the domain, add something bicubic that's + // zero at y=0 + if (y >= 0) + f += y * y * y * x * x * x; + + // in the top right quadrant, add something biquartic that's + // zero at x=0 and y=0 + if (x >= 0 && y >= 0) + f += x * x * x * x * y * y * y * y; + + return f; + } +}; + + +template +void +test() +{ + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::hyper_cube(tr, -1, 1); + tr.refine_global(3); + + const hp::FECollection fe(FE_Q(1), + FE_Q(2), + FE_Q(3), + FE_Q(4)); + hp::DoFHandler dof_handler(tr); + + // set DoF indices as described at the top of the file + for (auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + const double x = cell->center()[0]; + const double y = cell->center()[1]; + + if (x < 0 && y < 0) + cell->set_active_fe_index(0); + else if (x > 0 && y < 0) + cell->set_active_fe_index(1); + else if (x < 0 && y > 0) + cell->set_active_fe_index(2); + else if (x > 0 && y > 0) + cell->set_active_fe_index(3); + } + + dof_handler.distribute_dofs(fe); + + // interpolate the function above onto the finite element space + TrilinosWrappers::MPI::Vector interpolated(dof_handler.locally_owned_dofs(), + MPI_COMM_WORLD); + VectorTools::interpolate(dof_handler, CheckFunction(), interpolated); + + // then also apply constraints + AffineConstraints hanging_node_constraints; + DoFTools::make_hanging_node_constraints(dof_handler, + hanging_node_constraints); + hanging_node_constraints.close(); + hanging_node_constraints.distribute(interpolated); + + // extract a vector that has ghost elements + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs(dof_handler, relevant_set); + TrilinosWrappers::MPI::Vector x_rel(relevant_set, MPI_COMM_WORLD); + x_rel = interpolated; + + // Create a sufficiently high order quadrature formula + hp::QCollection quadrature(QGauss(3), + QGauss(4), + QGauss(5), + QGauss(6)); + + { + // integrate the difference between the function above and + // the zero function. for this case, we can compute the exact values + // by hand. the ones printed in the output are correct + Vector results(tr.n_active_cells()); + VectorTools::integrate_difference(dof_handler, + x_rel, + Functions::ZeroFunction(), + results, + quadrature, + VectorTools::L2_norm); + const double global = + VectorTools::compute_global_error(tr, results, VectorTools::L2_norm); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + deallog << "L2 norm = " << global << std::endl; + + Assert(std::fabs(global - std::sqrt(5847346.) / 2520. * + (dim == 3 ? std::sqrt(2) : 1)) < 1e-7, + ExcInternalError()); + } + + + { + // Now also integrate the difference between the function above and + // the its interpolant. This should then of course be zero + Vector results(tr.n_active_cells()); + VectorTools::integrate_difference(dof_handler, + x_rel, + CheckFunction(), + results, + quadrature, + VectorTools::L2_norm); + const double global = + VectorTools::compute_global_error(tr, results, VectorTools::L2_norm); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + deallog << "L2 error = " << global << std::endl; + + Assert(std::fabs(global) < 1e-15, ExcInternalError()); + } +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + initlog(); + + deallog.push("2d"); + test<2>(); + deallog.pop(); + + deallog.push("3d"); + test<3>(); + deallog.pop(); + } + else + { + deallog.push("2d"); + test<2>(); + deallog.pop(); + + deallog.push("3d"); + test<3>(); + deallog.pop(); + } +} diff --git a/tests/mpi/hp_integrate_difference.with_trilinos=true.mpirun=10.output b/tests/mpi/hp_integrate_difference.with_trilinos=true.mpirun=10.output new file mode 100644 index 0000000000..f37ee44b7f --- /dev/null +++ b/tests/mpi/hp_integrate_difference.with_trilinos=true.mpirun=10.output @@ -0,0 +1,5 @@ + +DEAL:0:2d::L2 norm = 0.959575 +DEAL:0:2d::L2 error = 3.84786e-16 +DEAL:0:3d::L2 norm = 1.35704 +DEAL:0:3d::L2 error = 6.75122e-16 diff --git a/tests/mpi/hp_integrate_difference.with_trilinos=true.mpirun=4.output b/tests/mpi/hp_integrate_difference.with_trilinos=true.mpirun=4.output new file mode 100644 index 0000000000..f37ee44b7f --- /dev/null +++ b/tests/mpi/hp_integrate_difference.with_trilinos=true.mpirun=4.output @@ -0,0 +1,5 @@ + +DEAL:0:2d::L2 norm = 0.959575 +DEAL:0:2d::L2 error = 3.84786e-16 +DEAL:0:3d::L2 norm = 1.35704 +DEAL:0:3d::L2 error = 6.75122e-16 -- 2.39.5