diff --git a/opm/core/utility/VelocityInterpolation.cpp b/opm/core/utility/VelocityInterpolation.cpp index b235fc77..c8678898 100644 --- a/opm/core/utility/VelocityInterpolation.cpp +++ b/opm/core/utility/VelocityInterpolation.cpp @@ -86,6 +86,15 @@ namespace Opm namespace { + /// Calculates the determinant of a 2 x 2 matrix, represented as + /// two two-dimensional arrays. + double determinantOf(const double* a0, + const double* a1) + { + return + a0[0] * a1[1] - a0[1] * a1[0]; + } + /// Calculates the determinant of a 3 x 3 matrix, represented as /// three three-dimensional arrays. double determinantOf(const double* a0, @@ -97,6 +106,16 @@ namespace Opm a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); } + + /// Calculates the volume of the parallelepiped given by + /// the vectors n[i] for i = 0..(dim-1), each n[i] is of size dim. + double cornerVolume(double** n, const int dim) + { + ASSERT(dim == 2 || dim == 3); + double det = (dim == 2) ? determinantOf(n[0], n[1]) : determinantOf(n[0], n[1], n[2]); + return std::fabs(det); + } + } // anonymous namespace @@ -139,8 +158,7 @@ namespace Opm CornerInfo ci; ci.corner_id = corner_id_count++;; ci.vertex = *it; - typedef double* DblPtr; - DblPtr fnorm[3] = { 0, 0, 0 }; + double* fnorm[Maxdim] = { 0 }; typedef std::multimap::const_iterator MMIt; std::pair frange = vertex_adj_faces.equal_range(ci.vertex); int fi = 0; @@ -155,7 +173,7 @@ namespace Opm } ASSERT(fi == dim); adj_faces_.insert(adj_faces_.end(), vert_adj_faces.begin(), vert_adj_faces.end()); - const double corner_vol = std::fabs(determinantOf(fnorm[0], fnorm[1], fnorm[2])); + const double corner_vol = cornerVolume(fnorm, dim); ci.volume = corner_vol; cell_corner_info.push_back(ci); std::sort(vert_adj_faces.begin(), vert_adj_faces.end()); diff --git a/tests/Makefile.am b/tests/Makefile.am index 1b6457c3..6cd55301 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -24,6 +24,7 @@ test_read_grid \ test_read_vag \ test_readpolymer \ test_sf2p \ +test_velocityinterpolation \ test_wells \ test_writeVtkData \ unit_test @@ -59,6 +60,8 @@ test_read_vag_SOURCES = test_read_vag.cpp test_sf2p_SOURCES = test_sf2p.cpp +test_velocityinterpolation_SOURCES = test_velocityinterpolation.cpp + test_writeVtkData_SOURCES = test_writeVtkData.cpp test_wells_SOURCES = test_wells.cpp diff --git a/tests/test_velocityinterpolation.cpp b/tests/test_velocityinterpolation.cpp new file mode 100644 index 00000000..372868ea --- /dev/null +++ b/tests/test_velocityinterpolation.cpp @@ -0,0 +1,482 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + 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 3 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 . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE VelocityInterpolationTest +#include + +#include +#include +#include +#include + +using namespace Opm; + +namespace +{ + + // Compute flux corresponding to a constant velocity vector v. + void computeFlux(const UnstructuredGrid& grid, const std::vector& v, std::vector& flux) + { + const int dim = v.size(); + ASSERT(dim == grid.dimensions); + flux.resize(grid.number_of_faces); + for (int face = 0; face < grid.number_of_faces; ++face) { + flux[face] = std::inner_product(v.begin(), v.end(), grid.face_normals + face*dim, 0.0); + } + } + + // Compute a linear vector function v = v0 + x*v1. + void computeLinearVec(const std::vector& v0, + const std::vector& v1, + const std::vector& x, + std::vector& v) + { + ASSERT(v0.size() == v1.size() && v0.size() == x.size()); + const int dim = v0.size(); + v.resize(dim); + for (int dd = 0; dd < dim; ++dd) { + v[dd] = v0[dd] + x[dd]*v1[dd]; + } + } + + + // Compute flux corresponding to a velocity vector v = v0 + x*v1. + void computeFluxLinear(const UnstructuredGrid& grid, + const std::vector& v0, + const std::vector& v1, + std::vector& flux) + { + const int dim = v0.size(); + ASSERT(dim == grid.dimensions); + flux.resize(grid.number_of_faces); + std::vector x(dim); + std::vector v(dim); + for (int face = 0; face < grid.number_of_faces; ++face) { + const double* fc = grid.face_centroids + face*dim; + std::copy(fc, fc + dim, x.begin()); + computeLinearVec(v0, v1, x, v); + flux[face] = std::inner_product(v.begin(), v.end(), grid.face_normals + face*dim, 0.0); + } + } + + + double vectorDiff2(const std::vector& v1, const std::vector& v2) + { + ASSERT(v1.size() == v2.size()); + const int sz = v1.size(); + double vdiff = 0.0; + for (int i = 0; i < sz; ++i) { + vdiff += (v1[i] - v2[i])*(v1[i] - v2[i]); + } + return vdiff; + } + +} // anonymous namespace + + +template +void testConstantVelRepro2d() +{ + // Set up 2d 1-cell cartesian case. + GridManager g(1, 1); + const UnstructuredGrid& grid = *g.c_grid(); + std::vector v(2); + v[0] = 0.12345; + v[1] = -0.6789; + std::vector flux; + computeFlux(grid, v, flux); + VelInterp vic(grid); + vic.setupFluxes(&flux[0]); + + // Test a few points + std::vector x(2); + x[0] = 0.23456; + x[1] = 0.87654; + std::vector v_interp(2); + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.5; + x[1] = 0.5; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 1.0; + x[1] = 0.5; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 1.0; + x[1] = 1.0; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); +} + + +namespace +{ + + // Data for a pyramid. Node numbering goes + // lexicographic on bottom, then top. + // Face numbering goes xmin, xmax, ymin, ymax, bottom. + namespace Pyramid + { + static int face_nodes[] = { 0, 4, 2, 3, 4, 1, 0, 1, 4, 4, 3, 2, 0, 2, 3, 1, }; + static int face_nodepos[] = { 0, 3, 6, 9, 12, 16 }; + static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; + static int cell_faces[] = { 0, 1, 2, 3, 4 }; + static int cell_facepos[] = { 0, 5 }; + static double node_coordinates[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0 }; + static double face_centroids[] = { 0, 1.0/3.0, 1.0/3.0, + 2.0/3.0, 1.0/3.0, 1.0/3.0, + 1.0/3.0, 0, 1.0/3.0, + 1.0/3.0, 2.0/3.0, 1.0/3.0, + 0.5, 0.5, 0 }; + static double face_areas[] = { 0.5, std::sqrt(2.0), 0.5, std::sqrt(2.0), 1.0 }; + static double face_normals[] = { -0.5000, 0, 0, + 0.5000, 0, 0.5000, + 0, -0.5000, 0, + 0, 0.5000, 0.5000, + 0, 0, -1.0000 }; + static double cell_centroids[] = { 0.375, 0.375, 0.25 }; + static double cell_volumes[] = { 1.0/3.0 }; + + } // namespace Pyramid + + UnstructuredGrid makePyramid() + { + // Make a 3d 1-cell grid, where the single cell is a pyramid. + UnstructuredGrid grid; + grid.dimensions = 3; + grid.number_of_cells = 1; + grid.number_of_faces = 5; + grid.number_of_nodes = 5; + grid.face_nodes = Pyramid::face_nodes; + grid.face_nodepos = Pyramid::face_nodepos; + grid.face_cells = Pyramid::face_cells; + grid.cell_faces = Pyramid::cell_faces; + grid.cell_facepos = Pyramid::cell_facepos; + grid.node_coordinates = Pyramid::node_coordinates; + grid.face_centroids = Pyramid::face_centroids; + grid.face_areas = Pyramid::face_areas; + grid.face_normals = Pyramid::face_normals; + grid.cell_centroids = Pyramid::cell_centroids; + grid.cell_volumes = Pyramid::cell_volumes; + return grid; + } + +} // anonymous namespace + + +template +void testConstantVelReproPyramid() +{ + // Set up a 3d 1-cell non-cartesian case (a pyramid). + UnstructuredGrid grid = makePyramid(); + std::vector v(3); + v[0] = 0.12345; + v[1] = -0.6789; + v[2] = 0.3456; + std::vector flux; + computeFlux(grid, v, flux); + VelInterp vic(grid); + vic.setupFluxes(&flux[0]); + + // Test a few points + std::vector x(3); + x[0] = 0.123; + x[1] = 0.0123; + x[2] = 0.213; + std::vector v_interp(3); + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.0; + x[1] = 0.0; + x[2] = 1.0; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.5; + x[1] = 0.5; + x[2] = 0.0; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.5; + x[1] = 0.5; + x[2] = 0.1; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); +} + + +namespace +{ + // Data for an irregular 2d polygon. + namespace Irreg2d + { + static int face_nodes[] = { 0, 1, 1, 2, 2, 3, 3, 4, 4, 0 }; + static int face_nodepos[] = { 0, 2, 4, 6, 8, 10 }; + static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; + static int cell_faces[] = { 0, 1, 2, 3, 4 }; + static int cell_facepos[] = { 0, 5 }; + static double node_coordinates[] = { 0, 0, 3, 0, 3, 2, 1, 3, 0, 2 }; + static double face_centroids[] = { 1.5, 0, 3, 1, 2, 2.5, 0.5, 2.5, 0, 1 }; + static double face_areas[] = { 3, 2, std::sqrt(5.0), std::sqrt(2.0), 2 }; + static double face_normals[] = { 0, -3, 2, 0, 1, 2, -1, 1, -2, 0 }; + static double cell_centroids[] = { 22.0/15.0, 19.0/15.0 }; + static double cell_volumes[] = { 7.5 }; + } // namespace Irreg2d + + UnstructuredGrid makeIrreg2d() + { + // Make a 2d 1-cell grid, where the single cell is a polygon. + UnstructuredGrid grid; + grid.dimensions = 2; + grid.number_of_cells = 1; + grid.number_of_faces = 5; + grid.number_of_nodes = 5; + grid.face_nodes = Irreg2d::face_nodes; + grid.face_nodepos = Irreg2d::face_nodepos; + grid.face_cells = Irreg2d::face_cells; + grid.cell_faces = Irreg2d::cell_faces; + grid.cell_facepos = Irreg2d::cell_facepos; + grid.node_coordinates = Irreg2d::node_coordinates; + grid.face_centroids = Irreg2d::face_centroids; + grid.face_areas = Irreg2d::face_areas; + grid.face_normals = Irreg2d::face_normals; + grid.cell_centroids = Irreg2d::cell_centroids; + grid.cell_volumes = Irreg2d::cell_volumes; + return grid; + } + +} // anonymous namespace + + +template +void testConstantVelReproIrreg2d() +{ + // Set up a 2d 1-cell non-cartesian case (a pyramid). + UnstructuredGrid grid = makeIrreg2d(); + std::vector v(2); + v[0] = 0.12345; + v[1] = -0.6789; + std::vector flux; + computeFlux(grid, v, flux); + VelInterp vic(grid); + vic.setupFluxes(&flux[0]); + + // Test a few points + std::vector x(2); + x[0] = 1.2345; + x[1] = 2.0123; + std::vector v_interp(2); + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.0; + x[1] = 0.0; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 1.0; + x[1] = 3.0; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 3.0; + x[1] = 1.0; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); +} + + + +namespace +{ + // Data for an irregular 3d prism. + namespace IrregPrism + { + static int face_nodes[] = { 0, 4, 2, 1, 3, 5, 0, 1, 5, 4, 2, 4, 5, 3, 2, 3, 0, 1}; + static int face_nodepos[] = { 0, 3, 6, 10, 14, 18 }; + static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; + static int cell_faces[] = { 0, 1, 2, 3, 4 }; + static int cell_facepos[] = { 0, 5 }; + static double node_coordinates[] = { 0, 0, 0, + 2, 0, 0, + 0, 1, 0, + 2, 1, 0, + 0, 0, 1, + 1, 0, 1 }; + static double face_centroids[] = { 0, 1.0/3.0, 1.0/3.0, + 5.0/3.0, 1.0/3.0, 1.0/3.0, + 7.0/9.0, 0, 4.0/9.0, + 7.0/9.0, 5.0/9.0, 4.0/9.0, + 1, 0.5, 0 }; + static double face_areas[] = { 0.500000000000000, + 0.707106781186548, + 1.500000000000000, + 2.121320343559642, + 2.000000000000000 }; + static double face_normals[] = { -0.500000000000000, 0, 0, + 0.500000000000000, 0.000000000000000, 0.500000000000000, + 0, -1.500000000000000, 0, + 0, 1.500000000000000, 1.500000000000000, + 0, 0, -2.000000000000000 }; + static double cell_centroids[] = { 0.85, 0.35, 0.3 }; + static double cell_volumes[] = { 5.0/6.0 }; + } // namespace IrregPrism + + UnstructuredGrid makeIrregPrism() + { + // Make a 3d 1-cell grid, where the single cell is a prism. + UnstructuredGrid grid; + grid.dimensions = 3; + grid.number_of_cells = 1; + grid.number_of_faces = 5; + grid.number_of_nodes = 6; + grid.face_nodes = IrregPrism::face_nodes; + grid.face_nodepos = IrregPrism::face_nodepos; + grid.face_cells = IrregPrism::face_cells; + grid.cell_faces = IrregPrism::cell_faces; + grid.cell_facepos = IrregPrism::cell_facepos; + grid.node_coordinates = IrregPrism::node_coordinates; + grid.face_centroids = IrregPrism::face_centroids; + grid.face_areas = IrregPrism::face_areas; + grid.face_normals = IrregPrism::face_normals; + grid.cell_centroids = IrregPrism::cell_centroids; + grid.cell_volumes = IrregPrism::cell_volumes; + return grid; + } + +} // anonymous namespace + +template +void testConstantVelReproIrregPrism() +{ + // Set up a 3d 1-cell non-cartesian case (a pyramid). + UnstructuredGrid grid = makeIrregPrism(); + std::vector v(3); + v[0] = 0.12345; + v[1] = -0.6789; + v[2] = 0.3456; + std::vector flux; + computeFlux(grid, v, flux); + VelInterp vic(grid); + vic.setupFluxes(&flux[0]); + + // Test a few points + std::vector x(3); + x[0] = 0.123; + x[1] = 0.0123; + x[2] = 0.213; + std::vector v_interp(3); + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.0; + x[1] = 0.0; + x[2] = 1.0; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 1.0; + x[1] = 0.0; + x[2] = 1.0; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.5; + x[1] = 0.5; + x[2] = 0.5; + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); +} + +template +void testLinearVelReproIrregPrism() +{ + // Set up a 3d 1-cell non-cartesian case (a pyramid). + UnstructuredGrid grid = makeIrregPrism(); + std::vector v0(3); + // v0[0] = 0.12345; + // v0[1] = -0.6789; + // v0[2] = 0.423; + v0[0] = 0.0; + v0[1] = 0.0; + v0[2] = 0.0; + std::vector v1(3); + // v1[0] = -0.1; + // v1[1] = 0.454; + // v1[2] = 0.21; + v1[0] = 0.0; + v1[1] = 0.0; + v1[2] = 1.0; + std::vector flux; + computeFluxLinear(grid, v0, v1, flux); + VelInterp vic(grid); + vic.setupFluxes(&flux[0]); + + // Test a few points + std::vector v(3); + std::vector x(3); + x[0] = 0.123; + x[1] = 0.0123; + x[2] = 0.213; + computeLinearVec(v0, v1, x, v); + std::vector v_interp(3); + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.0; + x[1] = 0.0; + x[2] = 1.0; + computeLinearVec(v0, v1, x, v); + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 1.0; + x[1] = 0.0; + x[2] = 1.0; + computeLinearVec(v0, v1, x, v); + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); + x[0] = 0.5; + x[1] = 0.5; + x[2] = 0.5; + computeLinearVec(v0, v1, x, v); + vic.interpolate(0, &x[0], &v_interp[0]); + BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); +} + + +BOOST_AUTO_TEST_CASE(test_VelocityInterpolationConstant) +{ + testConstantVelRepro2d(); + testConstantVelReproPyramid(); + testConstantVelReproIrreg2d(); + // testConstantVelReproIrregPrism(); // We should verify if this is expected to fail or not. +} + +BOOST_AUTO_TEST_CASE(test_VelocityInterpolationECVI) +{ + testConstantVelRepro2d(); + BOOST_CHECK_THROW(testConstantVelReproPyramid(), std::exception); + testConstantVelReproIrreg2d(); + testConstantVelReproIrregPrism(); + // Though the interpolation has linear precision, the corner velocity + // construction does not, so the below test cannot be expected to succeed. + // testLinearVelReproIrregPrism(); +} + +