Files renamed, forgot to add new files. Here they are
This commit is contained in:
Torbjrn Skille
parent
5374faa9b8
commit
bee5699d79
24
opm/common/utility/numeric/CalculateCellvol.hpp
Executable file
24
opm/common/utility/numeric/CalculateCellvol.hpp
Executable file
@ -0,0 +1,24 @@
|
||||
/*
|
||||
Copyright 2018 Statoil ASA.
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <vector>
|
||||
#include <math.h>
|
||||
|
||||
|
||||
double CalculateCellvol(std::vector<double>& X, std::vector<double>& Y, std::vector<double>& Z);
|
||||
|
||||
|
||||
153
src/opm/common/utility/numeric/CalculateCellvol.cpp
Executable file
153
src/opm/common/utility/numeric/CalculateCellvol.cpp
Executable file
@ -0,0 +1,153 @@
|
||||
/*
|
||||
Copyright 2018 Statoil ASA.
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
|
||||
#include <algorithm>
|
||||
#include <opm/common/utility/numeric/CalculateCellvol.hpp>
|
||||
|
||||
|
||||
/*
|
||||
Cell volume calculation based on following publication:
|
||||
|
||||
D. K Pointing, Corner Point Geometry in Reservoir Simulation ,
|
||||
ECMOR I - 1st European Conference on the Mathematics of Oil Recovery,
|
||||
1989
|
||||
*/
|
||||
|
||||
|
||||
/*
|
||||
The expressions {C(0,0,0),..C(1,1,1)} have a nice interpretation in
|
||||
terms of a type of multipole expansion - the last four terms are differences
|
||||
in the lengths of cell face diagonals and of the diagonals across the
|
||||
cell. For a cubical block only the first four terms would exist.
|
||||
|
||||
*/
|
||||
|
||||
double C(double* r, int i1, int i2, int i3){
|
||||
|
||||
double temp;
|
||||
|
||||
temp=0.0;
|
||||
|
||||
if ((i1==0) && (i2==0) && (i3==0)){
|
||||
temp=r[0];
|
||||
}
|
||||
|
||||
if ((i1==1) && (i2==0) && (i3==0)){
|
||||
temp=r[1]-r[0];
|
||||
}
|
||||
|
||||
if ((i1==0) && (i2==1) && (i3==0)){
|
||||
temp=r[2]-r[0];
|
||||
}
|
||||
|
||||
if ((i1==0) && (i2==0) && (i3==1)){
|
||||
temp=r[4]-r[0];
|
||||
}
|
||||
|
||||
if ((i1==1) && (i2==1) && (i3==0)){
|
||||
temp=r[3]+r[0]-r[2]-r[1];
|
||||
}
|
||||
|
||||
if ((i1==0) && (i2==1) && (i3==1)){
|
||||
temp=r[6]+r[0]-r[4]-r[2];
|
||||
}
|
||||
|
||||
if ((i1==1) && (i2==0) && (i3==1)){
|
||||
temp=r[5]+r[0]-r[4]-r[1];
|
||||
}
|
||||
|
||||
if ((i1==1) && (i2==1) && (i3==1)){
|
||||
temp=r[7]+r[4]+r[2]+r[1]-r[6]-r[5]-r[3]-r[0];
|
||||
}
|
||||
|
||||
return temp;
|
||||
|
||||
}
|
||||
|
||||
|
||||
double perm123sign(int i1, int i2, int i3){
|
||||
|
||||
double temp;
|
||||
|
||||
if ((i1==1) && (i2==2) && (i3==3)){
|
||||
temp=1.0;
|
||||
}
|
||||
|
||||
if ((i1==1) && (i2==3) && (i3==2)){
|
||||
temp=-1.0;
|
||||
}
|
||||
|
||||
if ((i1==2) && (i2==1) && (i3==3)){
|
||||
temp=-1.0;
|
||||
}
|
||||
|
||||
if ((i1==2) && (i2==3) && (i3==1)){
|
||||
temp=1.0;
|
||||
}
|
||||
|
||||
if ((i1==3) && (i2==1) && (i3==2)){
|
||||
temp=1.0;
|
||||
}
|
||||
|
||||
if ((i1==3) && (i2==2) && (i3==1)){
|
||||
temp=-1.0;
|
||||
}
|
||||
|
||||
return temp;
|
||||
}
|
||||
|
||||
double CalculateCellvol(std::vector<double>& X, std::vector<double>& Y, std::vector<double>& Z){
|
||||
|
||||
double volume=0.0;
|
||||
double* vect[3];
|
||||
int permutation[]={1, 2, 3};
|
||||
|
||||
do {
|
||||
for (int j=0;j<3;j++){
|
||||
if (permutation[j]==1){
|
||||
vect[j]=&X[0];
|
||||
} else if (permutation[j]==2){
|
||||
vect[j]=&Y[0];
|
||||
} else if (permutation[j]==3){
|
||||
vect[j]=&Z[0];
|
||||
}
|
||||
}
|
||||
|
||||
for (int pb=0;pb<2;pb++){
|
||||
for (int pg=0;pg<2;pg++){
|
||||
for (int qa=0;qa<2;qa++){
|
||||
for (int qg=0;qg<2;qg++){
|
||||
for (int ra=0;ra<2;ra++){
|
||||
for (int rb=0;rb<2;rb++){
|
||||
const double cprod =C(vect[0], 1, pb, pg)*C(vect[1], qa, 1, qg)*C(vect[2], ra, rb, 1);
|
||||
const double denom =(qa+ra+1) * (pb+rb+1) * (pg+qg+1);
|
||||
volume += perm123sign(permutation[0], permutation[1], permutation[2]) * cprod / denom;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} while (std::next_permutation(std::begin(permutation), std::end(permutation)));
|
||||
|
||||
return std::fabs(volume);
|
||||
}
|
||||
|
||||
|
||||
|
||||
63
tests/test_CalculateCellvol.cpp
Executable file
63
tests/test_CalculateCellvol.cpp
Executable file
@ -0,0 +1,63 @@
|
||||
|
||||
/*
|
||||
Copyright 2018 Statoil ASA.
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
|
||||
#define NVERBOSE // Suppress own messages when throw()ing
|
||||
|
||||
#define BOOST_TEST_MODULE CubicTest
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <boost/test/floating_point_comparison.hpp>
|
||||
|
||||
/* --- our own headers --- */
|
||||
#include <opm/common/utility/numeric/CalculateCellvol.hpp>
|
||||
|
||||
|
||||
|
||||
//using namespace Opm;
|
||||
|
||||
BOOST_AUTO_TEST_SUITE ()
|
||||
|
||||
BOOST_AUTO_TEST_CASE (calc_cellvol)
|
||||
{
|
||||
|
||||
/* This is four example cells with different geometries. */
|
||||
|
||||
std::vector<double> x1 {488100.140035, 488196.664549, 488085.584866, 488182.365605, 488099.065709, 488195.880889, 488084.559409, 488181.633495};
|
||||
std::vector<double> y1 {6692539.945578, 6692550.834909, 6692638.574346, 6692650.086244, 6692538.810649, 6692550.080826, 6692637.628127, 6692649.429649};
|
||||
std::vector<double> z1 {2841.856000, 2840.138000, 2842.042000, 2839.816000, 2846.142000, 2844.252000, 2846.244000, 2843.868000};
|
||||
|
||||
std::vector<double> x2 {489539.050892, 489638.562319, 489527.025803, 489627.914272, 489537.173839, 489638.562319, 489524.119410, 489627.914272};
|
||||
std::vector<double> y2 {6695907.426615, 6695923.399538, 6696003.688945, 6696020.073168, 6695908.526577, 6695923.399538, 6696005.533276, 6696020.073168};
|
||||
std::vector<double> z2 {2652.859000, 2651.765000, 2652.381000, 2652.608000, 2655.276000, 2651.765000, 2656.381000, 2652.608000};
|
||||
|
||||
std::vector<double> x3 {486819.009056, 486904.877179, 486812.975440, 486900.527416, 486818.450563, 486903.978383, 486812.975440, 486900.527416};
|
||||
std::vector<double> y3 {6694966.253697, 6694998.360834, 6695061.104896, 6695086.717018, 6694967.135567, 6695000.197284, 6695061.104896, 6695086.717018};
|
||||
std::vector<double> z3 {2795.193000, 2799.997000, 2792.240000, 2793.972000, 2795.671000, 2800.927000, 2792.240000, 2793.972000};
|
||||
|
||||
std::vector<double> x4 {489194.711544, 489259.232449, 489206.220219, 489294.099099, 489194.711544, 489254.725623, 489201.723535, 489289.423088};
|
||||
std::vector<double> y4 {6694510.504054, 6694525.862296, 6694608.410134, 6694621.029382, 6694510.504054, 6694526.272988, 6694608.819829, 6694621.467114};
|
||||
std::vector<double> z4 {2740.7340, 2754.7810, 2730.0150, 2726.1080, 2740.7340, 2758.3530, 2733.9490, 2730.1080};
|
||||
|
||||
BOOST_REQUIRE_CLOSE (CalculateCellvol(x1,y1,z1), 40368.7852157, 1e-9);
|
||||
BOOST_REQUIRE_CLOSE (CalculateCellvol(x2,y2,z2), 15766.9187847524, 1e-9);
|
||||
BOOST_REQUIRE_CLOSE (CalculateCellvol(x3,y3,z3), 3268.8819007839, 1e-9);
|
||||
BOOST_REQUIRE_CLOSE (CalculateCellvol(x4,y4,z4), 23391.4917234564, 1e-9);
|
||||
}
|
||||
|
||||
BOOST_AUTO_TEST_SUITE_END()
|
||||
Loading…
Reference in New Issue
Block a user