mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
86bc13912d
most warnings were in DUNE and ALUGrid, but these have been ignored. Also fixing some of these warnings (in particular the "parameter foo is unused" ones) would make the code harder to read and understand, so they have been ignored, too...
355 lines
11 KiB
C++
355 lines
11 KiB
C++
/*
|
|
Copyright (C) 2009-2013 by Andreas Lauser
|
|
|
|
This file is part of the Open Porous Media project (OPM).
|
|
|
|
OPM 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 2 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
/*!
|
|
* \file
|
|
* \brief A test for numerical integration using the vertex-centered finite
|
|
* volume geometries.
|
|
*/
|
|
#include "config.h"
|
|
|
|
#include <dune/common/exceptions.hh>
|
|
#include <dune/grid/yaspgrid.hh>
|
|
#include <dune/grid/io/file/vtk/vtkwriter.hh>
|
|
|
|
#if HAVE_ALUGRID
|
|
#include <dune/grid/alugrid.hh>
|
|
#endif // HAVE_ALUGRID
|
|
|
|
#include <dune/common/version.hh>
|
|
#include <dune/geometry/quadraturerules.hh>
|
|
|
|
#include <ewoms/disc/vcfv/vcfvstencil.hh>
|
|
|
|
const unsigned dim = 3;
|
|
typedef double Scalar;
|
|
typedef Ewoms::QuadrialteralQuadratureGeometry<Scalar, dim> QuadratureGeom;
|
|
typedef QuadratureGeom::LocalPosition LocalPosition;
|
|
typedef QuadratureGeom::GlobalPosition GlobalPosition;
|
|
|
|
// function prototypes
|
|
GlobalPosition::field_type f(const GlobalPosition &pos);
|
|
void testIdenityMapping();
|
|
template <class Grid>
|
|
void writeTetrahedronSubControlVolumes(const Grid &grid);
|
|
void testTetrahedron();
|
|
template <class Grid>
|
|
void writeCubeSubControlVolumes(const Grid &grid);
|
|
void testCube();
|
|
void testQuadrature();
|
|
|
|
GlobalPosition::field_type f(const GlobalPosition &pos)
|
|
{
|
|
GlobalPosition::field_type result = 1;
|
|
for (int i = 0; i < GlobalPosition::dimension; ++i)
|
|
result *= pos[i];
|
|
return result;
|
|
}
|
|
|
|
void testIdenityMapping()
|
|
{
|
|
QuadratureGeom foo;
|
|
|
|
Scalar corners[][3] = { { 0, 0, 0 },
|
|
{ 1, 0, 0 },
|
|
{ 0, 1, 0 },
|
|
{ 1, 1, 0 },
|
|
{ 0, 0, 1 },
|
|
{ 1, 0, 1 },
|
|
{ 0, 1, 1 },
|
|
{ 1, 1, 1 } };
|
|
foo.setCorners(corners, 8);
|
|
|
|
std::cout << "testing identity mapping...\n";
|
|
int n = 100;
|
|
for (int i = 0; i < n; ++i) {
|
|
for (int j = 0; j < n; ++j) {
|
|
for (int k = 0; k < n; ++k) {
|
|
LocalPosition localPos;
|
|
|
|
localPos[0] = Scalar(i) / (n - 1);
|
|
localPos[1] = Scalar(j) / (n - 1);
|
|
localPos[2] = Scalar(k) / (n - 1);
|
|
|
|
GlobalPosition globalPos = foo.global(localPos);
|
|
|
|
GlobalPosition diff(localPos);
|
|
diff -= globalPos;
|
|
assert(diff.two_norm() < 1e-10);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template <class Grid>
|
|
void writeTetrahedronSubControlVolumes(const Grid &grid)
|
|
{
|
|
#if HAVE_ALUGRID
|
|
typedef typename Grid::LeafGridView GridView;
|
|
|
|
typedef Dune::ALUGrid<dim, dim, Dune::cube, Dune::nonconforming> Grid2;
|
|
typedef typename Grid2::LeafGridView GridView2;
|
|
typedef Dune::GridFactory<Grid2> GridFactory2;
|
|
|
|
GridFactory2 gf2;
|
|
const auto &gridView = grid.leafView();
|
|
typedef Ewoms::VcfvStencil<Scalar, GridView> Stencil;
|
|
Stencil stencil(gridView);
|
|
|
|
auto eIt = gridView.template begin<0>();
|
|
const auto &eEndIt = gridView.template end<0>();
|
|
for (; eIt != eEndIt; ++eIt) {
|
|
stencil.update(*eIt);
|
|
for (int scvIdx = 0; scvIdx < stencil.numDof(); ++scvIdx) {
|
|
const auto &scvLocalGeom = stencil.subControlVolume(scvIdx).localGeometry();
|
|
|
|
for (int i = 0; i < scvLocalGeom.numCorners; ++i) {
|
|
GlobalPosition pos(
|
|
eIt->geometry().global(scvLocalGeom.corner(i)));
|
|
gf2.insertVertex(pos);
|
|
}
|
|
}
|
|
}
|
|
|
|
int cornerOffset = 0;
|
|
eIt = gridView.template begin<0>();
|
|
for (; eIt != eEndIt; ++eIt) {
|
|
stencil.update(*eIt);
|
|
for (int scvIdx = 0; scvIdx < stencil.numDof(); ++scvIdx) {
|
|
const auto &scvLocalGeom = stencil.subControlVolume(scvIdx).localGeometry();
|
|
|
|
std::vector<unsigned int> vertexIndices;
|
|
for (int i = 0; i < scvLocalGeom.numCorners; ++i) {
|
|
vertexIndices.push_back(cornerOffset);
|
|
++cornerOffset;
|
|
}
|
|
|
|
gf2.insertElement(Dune::GeometryType(Dune::GeometryType::cube, dim),
|
|
vertexIndices);
|
|
}
|
|
}
|
|
|
|
const auto &grid2 = *gf2.createGrid();
|
|
typedef Dune::VTKWriter<GridView2> VtkWriter;
|
|
VtkWriter writer(grid2.leafView(), Dune::VTK::conforming);
|
|
writer.write("tetrahedron-scvs", Dune::VTK::ascii);
|
|
#endif // HAVE_ALUGRID
|
|
}
|
|
|
|
void testTetrahedron()
|
|
{
|
|
#if HAVE_ALUGRID
|
|
typedef Dune::ALUGrid<dim, dim, Dune::simplex, Dune::nonconforming> Grid;
|
|
typedef Dune::GridFactory<Grid> GridFactory;
|
|
GridFactory gf;
|
|
Scalar corners[][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
|
|
|
|
for (unsigned i = 0; i < sizeof(corners) / sizeof(corners[0]); ++i) {
|
|
GlobalPosition pos;
|
|
for (unsigned j = 0; j < dim; ++j)
|
|
pos[j] = corners[i][j];
|
|
gf.insertVertex(pos);
|
|
}
|
|
std::vector<unsigned int> v = { 0, 1, 2, 3 };
|
|
gf.insertElement(Dune::GeometryType(Dune::GeometryType::simplex, dim), v);
|
|
auto *grid = gf.createGrid();
|
|
|
|
// write the sub-control volumes to a VTK file.
|
|
writeTetrahedronSubControlVolumes(*grid);
|
|
|
|
delete grid;
|
|
#endif // HAVE_ALUGRID
|
|
}
|
|
|
|
template <class Grid>
|
|
void writeCubeSubControlVolumes(const Grid &grid)
|
|
{
|
|
#if HAVE_ALUGRID
|
|
typedef typename Grid::LeafGridView GridView;
|
|
|
|
typedef Dune::ALUGrid<dim, dim, Dune::cube, Dune::nonconforming> Grid2;
|
|
typedef typename Grid2::LeafGridView GridView2;
|
|
typedef Dune::GridFactory<Grid2> GridFactory2;
|
|
typedef Ewoms::VcfvStencil<Scalar, GridView> Stencil;
|
|
|
|
GridFactory2 gf2;
|
|
const auto &gridView = grid.leafView();
|
|
Stencil stencil(gridView);
|
|
auto eIt = gridView.template begin<0>();
|
|
const auto &eEndIt = gridView.template end<0>();
|
|
for (; eIt != eEndIt; ++eIt) {
|
|
stencil.update(*eIt);
|
|
for (int scvIdx = 0; scvIdx < stencil.numDof(); ++scvIdx) {
|
|
const auto &scvLocalGeom = stencil.subControlVolume(scvIdx).localGeometry();
|
|
|
|
for (int i = 0; i < scvLocalGeom.numCorners; ++i) {
|
|
GlobalPosition pos(
|
|
eIt->geometry().global(scvLocalGeom.corner(i)));
|
|
gf2.insertVertex(pos);
|
|
}
|
|
}
|
|
}
|
|
|
|
int cornerOffset = 0;
|
|
eIt = gridView.template begin<0>();
|
|
for (; eIt != eEndIt; ++eIt) {
|
|
stencil.update(*eIt);
|
|
for (int scvIdx = 0; scvIdx < stencil.numDof(); ++scvIdx) {
|
|
const auto &scvLocalGeom = stencil.subControlVolume(scvIdx).localGeometry();
|
|
|
|
std::vector<unsigned int> vertexIndices;
|
|
for (int i = 0; i < scvLocalGeom.numCorners; ++i) {
|
|
vertexIndices.push_back(cornerOffset);
|
|
++cornerOffset;
|
|
}
|
|
|
|
gf2.insertElement(Dune::GeometryType(Dune::GeometryType::cube, dim),
|
|
vertexIndices);
|
|
}
|
|
}
|
|
|
|
const auto &grid2 = *gf2.createGrid();
|
|
typedef Dune::VTKWriter<GridView2> VtkWriter;
|
|
VtkWriter writer(grid2.leafView(), Dune::VTK::conforming);
|
|
writer.write("cube-scvs", Dune::VTK::ascii);
|
|
#endif // HAVE_ALUGRID
|
|
}
|
|
|
|
void testCube()
|
|
{
|
|
#if HAVE_ALUGRID
|
|
typedef Dune::ALUGrid<dim, dim, Dune::cube, Dune::nonconforming> Grid;
|
|
typedef Dune::GridFactory<Grid> GridFactory;
|
|
GridFactory gf;
|
|
Scalar corners[][3] = { { 0, 0, 0 },
|
|
{ 1, 0, 0 },
|
|
{ 0, 2, 0 },
|
|
{ 3, 3, 0 },
|
|
{ 0, 0, 4 },
|
|
{ 5, 0, 5 },
|
|
{ 0, 6, 6 },
|
|
{ 7, 7, 7 }, };
|
|
|
|
for (unsigned i = 0; i < sizeof(corners) / sizeof(corners[0]); ++i) {
|
|
GlobalPosition pos;
|
|
for (unsigned j = 0; j < dim; ++j)
|
|
pos[j] = corners[i][j];
|
|
gf.insertVertex(pos);
|
|
}
|
|
std::vector<unsigned int> v = { 0, 1, 2, 3, 4, 5, 6, 7 };
|
|
gf.insertElement(Dune::GeometryType(Dune::GeometryType::cube, dim), v);
|
|
auto *grid = gf.createGrid();
|
|
|
|
// write the sub-control volumes to a VTK file.
|
|
writeCubeSubControlVolumes(*grid);
|
|
|
|
delete grid;
|
|
#endif // HAVE_ALUGRID
|
|
}
|
|
|
|
void testQuadrature()
|
|
{
|
|
std::cout << "testing SCV quadrature...\n";
|
|
|
|
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
|
|
std::bitset<dim> isPeriodic(false);
|
|
std::array<int, dim> cellRes;
|
|
#else
|
|
Dune::FieldVector<bool, dim> isPeriodic(false);
|
|
Dune::FieldVector<int, dim> cellRes;
|
|
#endif
|
|
std::fill(cellRes.begin(), cellRes.end(), 10);
|
|
|
|
GlobalPosition upperRight(1.0);
|
|
|
|
typedef Dune::YaspGrid<dim> Grid;
|
|
typedef Grid::LeafGridView GridView;
|
|
Grid grid(
|
|
#ifdef HAVE_MPI
|
|
Dune::MPIHelper::getCommunicator(),
|
|
#endif
|
|
upperRight, // upper right
|
|
cellRes, // number of cells
|
|
isPeriodic, 0); // overlap
|
|
|
|
// compute approximate integral
|
|
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
|
|
auto gridView = grid.leafGridView();
|
|
#else
|
|
auto gridView = grid.leafView();
|
|
#endif
|
|
auto eIt = gridView.begin<0>();
|
|
const auto eEndIt = gridView.end<0>();
|
|
Scalar result = 0;
|
|
// instanciate a stencil
|
|
typedef Ewoms::VcfvStencil<Scalar, GridView> Stencil;
|
|
Stencil stencil(gridView);
|
|
for (; eIt != eEndIt; ++eIt) {
|
|
const auto &elemGeom = eIt->geometry();
|
|
|
|
stencil.update(*eIt);
|
|
|
|
// loop over all sub-control volumes
|
|
for (int scvIdx = 0; scvIdx < stencil.numDof(); ++scvIdx) {
|
|
const auto &scvLocalGeom = stencil.subControlVolume(scvIdx).localGeometry();
|
|
|
|
Dune::GeometryType geomType = scvLocalGeom.type();
|
|
static const int quadratureOrder = 2;
|
|
const auto &rule
|
|
= Dune::QuadratureRules<Scalar, dim>::rule(geomType,
|
|
quadratureOrder);
|
|
|
|
// integrate f over the sub-control volume
|
|
for (auto it = rule.begin(); it != rule.end(); ++it) {
|
|
auto posScvLocal = it->position();
|
|
auto posElemLocal = scvLocalGeom.global(posScvLocal);
|
|
auto posGlobal = elemGeom.global(posScvLocal);
|
|
|
|
Scalar fval = f(posGlobal);
|
|
Scalar weight = it->weight();
|
|
Scalar detjac = scvLocalGeom.integrationElement(posScvLocal)
|
|
* elemGeom.integrationElement(posElemLocal);
|
|
|
|
result += fval * weight * detjac;
|
|
}
|
|
}
|
|
}
|
|
|
|
std::cout << "result: " << result
|
|
<< " (expected value: " << 1.0 / (1 << dim) << ")\n";
|
|
}
|
|
|
|
int main(int argc, char **argv)
|
|
{
|
|
// initialize MPI, finalize is done automatically on exit
|
|
Dune::MPIHelper::instance(argc, argv);
|
|
|
|
testIdenityMapping();
|
|
// test the quadrature in a tetrahedron. since the CLang compiler
|
|
// prior to version 3.2 generates incorrect code here, we do not
|
|
// do it if the compiler is clang 3.1 or older.
|
|
#if !__clang__ || __clang_major__ > 3 \
|
|
|| (__clang_major__ == 3 && __clang_minor__ >= 3)
|
|
testTetrahedron();
|
|
testCube();
|
|
#endif
|
|
testQuadrature();
|
|
|
|
return 0;
|
|
}
|