diff --git a/test/discretization/CMakeLists.txt b/test/discretization/CMakeLists.txt index f127b603f79c94575badc24a7eb39c9328c88892..25c1bdd297e74a772cd7d985038aa7cde682c038 100644 --- a/test/discretization/CMakeLists.txt +++ b/test/discretization/CMakeLists.txt @@ -2,3 +2,4 @@ add_subdirectory(cellcentered) add_subdirectory(staggered) add_subdirectory(box) add_subdirectory(projection) +add_subdirectory(rotationsymmetry) diff --git a/test/discretization/rotationsymmetry/CMakeLists.txt b/test/discretization/rotationsymmetry/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..a8237b62853ccb9adcd4b836bc631644bef27442 --- /dev/null +++ b/test/discretization/rotationsymmetry/CMakeLists.txt @@ -0,0 +1,2 @@ +dune_add_test(SOURCES test_rotationsymmetric_gridgeometry.cc + LABELS unit discretization) diff --git a/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc new file mode 100644 index 0000000000000000000000000000000000000000..15072fb5227dfca3f54eb95dc4881cb5d35ca793 --- /dev/null +++ b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc @@ -0,0 +1,208 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \brief Test for rotation symmetric grid geometry + */ +#include <config.h> + +#include <iostream> + +#include <dune/common/parallel/mpihelper.hh> +#include <dune/common/exceptions.hh> +#include <dune/common/float_cmp.hh> + +#include <dune/grid/utility/structuredgridfactory.hh> +#include <dune/grid/yaspgrid.hh> +#include <dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh> +#include <dumux/discretization/rotationpolicy.hh> +#include <dumux/discretization/rotationsymmetricscv.hh> +#include <dumux/discretization/rotationsymmetricscvf.hh> + +namespace Dumux { + +template<class BaseTraits, RotationPolicy rotPolicy> +struct GGRotSymTraits : public BaseTraits +{ + using SubControlVolume = RotationSymmetricSubControlVolume<typename BaseTraits::SubControlVolume, rotPolicy>; + using SubControlVolumeFace = RotationSymmetricSubControlVolumeFace<typename BaseTraits::SubControlVolumeFace, rotPolicy>; +}; + +template<class GG> +void runTest(const GG& gg, const double refVolume, const double refSurface) +{ + double volume = 0.0; + double surface = 0.0; + for (const auto& element : elements(gg.gridView())) + { + auto fvGeometry = localView(gg); + fvGeometry.bind(element); + + for (const auto& scv : scvs(fvGeometry)) + volume += scv.volume(); + + for (const auto& scvf : scvfs(fvGeometry)) + if (scvf.boundary()) + surface += scvf.area(); + } + + // compare to reference + if (!Dune::FloatCmp::eq(refVolume, volume)) + DUNE_THROW(Dune::Exception, "Volume not correct! Reference: " << refVolume << ", computed: " << volume); + if (!Dune::FloatCmp::eq(refSurface, surface)) + DUNE_THROW(Dune::Exception, "Surface not correct! Reference: " << refSurface << ", computed: " << surface); +} + +} // end namespace Dumux + +int main (int argc, char *argv[]) try +{ + using namespace Dumux; + + // maybe initialize mpi + Dune::MPIHelper::instance(argc, argv); + + // test the disc policy + { + using Grid = Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>; + using GGTraits = GGRotSymTraits<CCTpfaDefaultGridGeometryTraits<typename Grid::LeafGridView>, RotationPolicy::disc>; + using GridGeometry = CCTpfaFVGridGeometry<typename Grid::LeafGridView, /*caching=*/false, GGTraits>; + using GlobalPosition = typename GridGeometry::SubControlVolume::GlobalPosition; + + // make a grid + const double innerRadius = 0.1; + const double outerRadius = 1.0; + GlobalPosition lower(innerRadius); + GlobalPosition upper(outerRadius); + std::array<unsigned int, Grid::dimension> els{{10}}; + std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els); + + // obtain leaf and make GridGeometry + auto leafGridView = grid->leafGridView(); + GridGeometry gg(leafGridView); + gg.update(); + + // compute the annulus area and the surface + const double refVolume = M_PI*(outerRadius*outerRadius - innerRadius*innerRadius); + const double refSurface = 2.0*M_PI*(innerRadius + outerRadius); + runTest(gg, refVolume, refSurface); + + std::cout << "Successfully tested disc (annulus) policy." << std::endl; + + } // end disc policy + + // test the ball policy + { + using Grid = Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>; + using GGTraits = GGRotSymTraits<CCTpfaDefaultGridGeometryTraits<typename Grid::LeafGridView>, RotationPolicy::ball>; + using GridGeometry = CCTpfaFVGridGeometry<typename Grid::LeafGridView, /*caching=*/false, GGTraits>; + using GlobalPosition = typename GridGeometry::SubControlVolume::GlobalPosition; + + // make a grid + const double innerRadius = 0.1; + const double outerRadius = 1.0; + GlobalPosition lower(innerRadius); + GlobalPosition upper(outerRadius); + std::array<unsigned int, Grid::dimension> els{{10}}; + std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els); + + // obtain leaf and make GridGeometry + auto leafGridView = grid->leafGridView(); + GridGeometry gg(leafGridView); + gg.update(); + + // compute the annulus area and the surface + const double refVolume = 4.0/3.0*M_PI*(outerRadius*outerRadius*outerRadius - innerRadius*innerRadius*innerRadius); + const double refSurface = 4.0*M_PI*(innerRadius*innerRadius + outerRadius*outerRadius); + runTest(gg, refVolume, refSurface); + + std::cout << "Successfully tested ball (shell) policy." << std::endl; + + } // end ball policy + + // test the toroid policy + { + using Grid = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>; + using GGTraits = GGRotSymTraits<CCTpfaDefaultGridGeometryTraits<typename Grid::LeafGridView>, RotationPolicy::toroid>; + using GridGeometry = CCTpfaFVGridGeometry<typename Grid::LeafGridView, /*caching=*/false, GGTraits>; + using GlobalPosition = typename GridGeometry::SubControlVolume::GlobalPosition; + + // make a grid + const double innerRadius = 0.1; + const double outerRadius = 1.0; + GlobalPosition lower({innerRadius, innerRadius}); + GlobalPosition upper({outerRadius, outerRadius}); + std::array<unsigned int, Grid::dimension> els({20, 20}); + std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els); + + // obtain leaf and make GridGeometry + auto leafGridView = grid->leafGridView(); + GridGeometry gg(leafGridView); + gg.update(); + + // compute the annulus area and the surface + const auto centroidRadius = 0.5*(innerRadius + outerRadius); + const auto side = (outerRadius - innerRadius); + const double refVolume = side*side*2.0*M_PI*centroidRadius; + const double refSurface = 4.0*side*2.0*M_PI*centroidRadius; + runTest(gg, refVolume, refSurface); + + std::cout << "Successfully tested toroid policy." << std::endl; + + } // end toroid policy + + // test the toroid policy for perfect cylinder + { + using Grid = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>; + using GGTraits = GGRotSymTraits<CCTpfaDefaultGridGeometryTraits<typename Grid::LeafGridView>, RotationPolicy::toroid>; + using GridGeometry = CCTpfaFVGridGeometry<typename Grid::LeafGridView, /*caching=*/false, GGTraits>; + using GlobalPosition = typename GridGeometry::SubControlVolume::GlobalPosition; + + // make a grid + const double innerRadius = 0.0; + const double outerRadius = 1.0; + GlobalPosition lower({innerRadius, innerRadius}); + GlobalPosition upper({outerRadius, outerRadius}); + std::array<unsigned int, Grid::dimension> els({20, 20}); + std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els); + + // obtain leaf and make GridGeometry + auto leafGridView = grid->leafGridView(); + GridGeometry gg(leafGridView); + gg.update(); + + // compute the annulus area and the surface + const double refVolume = outerRadius*M_PI*outerRadius*outerRadius; + const double refSurface = 2.0*M_PI*outerRadius + 2.0*M_PI*outerRadius*outerRadius; + runTest(gg, refVolume, refSurface); + + std::cout << "Successfully tested toroid policy for perfect cylinder." << std::endl; + + } // end toroid policy + + return 0; +} +// ////////////////////////////////// +// Error handler +// ///////////////////////////////// +catch (Dune::Exception &e) { + + std::cout << e << std::endl; + return 1; +}