OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Changed to extrapolation outside tables for pormult. Added facetags to the grid structure. Changed default fluid to Linear.

Browse Source
This commit is contained in:
Halvor Møll Nilsen 2012-08-14 09:36:19 +02:00
parent 3d45ca3512
commit 064184bd0c
8 changed files with 89 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
opm/core/fluid/RockCompressibility.cpp
View File

@ -69,7 +69,8 @@ namespace Opm
const double cpnorm = rock_comp_*(pressure - pref_);
return (1.0 + cpnorm + 0.5*cpnorm*cpnorm);
} else {
return Opm::linearInterpolation(p_, poromult_, pressure);
// return Opm::linearInterpolation(p_, poromult_, pressure);
return Opm::linearInterpolationExtrap(p_, poromult_, pressure);
}
}
@ -78,8 +79,11 @@ namespace Opm
if (p_.empty()) {
return rock_comp_;
} else {
const double poromult = Opm::linearInterpolation(p_, poromult_, pressure);
const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure);
//const double poromult = Opm::linearInterpolation(p_, poromult_, pressure);
//const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure);
const double poromult = Opm::linearInterpolationExtrap(p_, poromult_, pressure);
const double dporomultdp = Opm::linearInterpolationDerivativeExtrap(p_, poromult_, pressure);
return dporomultdp/poromult;
}
}

3
opm/core/fluid/SaturationPropsBasic.cpp
View File

@ -114,7 +114,8 @@ namespace Opm
THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases);
}
num_phases_ = num_phases;
std::string rpf = param.getDefault("relperm_func", std::string("Unset"));
//std::string rpf = param.getDefault("relperm_func", std::string("Unset"));
std::string rpf = param.getDefault("relperm_func", std::string("Linear"));
if (rpf == "Constant") {
relperm_func_ = Constant;
if(num_phases!=1){

3
opm/core/grid.c
View File

@ -104,6 +104,8 @@ allocate_grid(size_t ndims ,
nel = ncellfaces;
G->cell_faces = malloc(nel * sizeof *G->cell_faces);
G->cell_facetag = malloc(nel * sizeof *G->cell_facetag);
nel = ncells + 1;
G->cell_facepos = malloc(nel * sizeof *G->cell_facepos);
@ -121,6 +123,7 @@ allocate_grid(size_t ndims ,
(G->face_normals == NULL) ||
(G->face_areas == NULL) ||
(G->cell_faces == NULL) ||
(G->cell_facetag == NULL) ||
(G->cell_facepos == NULL) ||
(G->cell_centroids == NULL) ||
(G->cell_volumes == NULL) )

13
opm/core/grid/cart_grid.c
View File

@ -262,6 +262,13 @@ fill_cart_topology_3d(struct UnstructuredGrid *G)
}
}
}
for (k=0; k< nx*ny*nz;++k){
for (i=0; i < 6; ++i){
G->cell_facetag[k*6+i]=i;
}
}
fnodes = G->face_nodes;
@ -565,6 +572,12 @@ fill_cart_topology_2d(struct UnstructuredGrid *G)
}
}
for (j=0; j< nx*ny;++j){
G->cell_facetag[j*4+0]=0;
G->cell_facetag[j*4+1]=2;
G->cell_facetag[j*4+2]=1;
G->cell_facetag[j*4+3]=3;
}
fnodes = G->face_nodes;

15
opm/core/linalg/LinearSolverAGMG.cpp
View File

@ -65,12 +65,21 @@ namespace Opm
rtol_ (rtol) ,
is_spd_(is_spd)
{
#if !HAVE_AGMG
THROW("AGMG support is not enabled in this library");
#endif // HAVE_AGMG
}
LinearSolverAGMG::LinearSolverAGMG(const parameter::ParameterGroup& param)
: max_it_(100) ,
rtol_ (1.0e-6),
is_spd_(false)
{
max_it_ = param.getDefault("max_it", max_it_);
rtol_ = param.getDefault("rtol" , rtol_ );
is_spd_ = param.getDefault("is_spd", is_spd_);
}
LinearSolverAGMG::~LinearSolverAGMG() {}
LinearSolverInterface::LinearSolverReport
LinearSolverAGMG::solve(const int size ,

9
opm/core/linalg/LinearSolverAGMG.hpp
View File

@ -43,6 +43,8 @@
*/
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <string>
namespace Opm
{
@ -61,7 +63,12 @@ namespace Opm
LinearSolverAGMG(const int max_it = 100 ,
const double rtol = 1.0e-6,
const bool is_spd = false);
/**
* Constructor.
* \param[in] param ParameterGroup object contianing the fields
* max_it,rtol,is_spd as used in the constructor
*/
LinearSolverAGMG(const parameter::ParameterGroup& param);
/// Destructor.
virtual ~LinearSolverAGMG();

12
opm/core/linalg/LinearSolverFactory.cpp
View File

@ -30,6 +30,10 @@
#if HAVE_DUNE_ISTL
#include <opm/core/linalg/LinearSolverIstl.hpp>
#endif
#if HAVE_AGMG
#include <opm/core/linalg/LinearSolverAGMG.hpp>
#endif
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
@ -69,8 +73,12 @@ namespace Opm
solver_.reset(new LinearSolverIstl(param));
#endif
}
else {
else if (ls == "agmg") {
#if HAVE_AGMG
solver_.reset(new LinearSolverAGMG(param));
#endif
}
else {
THROW("Linear solver " << ls << " is unknown.");
}

34
opm/core/utility/linearInterpolation.hpp
View File

@ -85,6 +85,40 @@ namespace Opm
}
}
template <typename T>
T linearInterpolationDerivativeExtrap(const std::vector<double>& xv,
const std::vector<T>& yv,
double x)
{
std::vector<double>::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x);
int lb_ix = lb - xv.begin();
int nend = int(xv.size());
if (lb_ix == 0) {
return (yv[1]-yv[0])/(xv[1]-xv[0]);
} else if (lb_ix == int(xv.size())) {
return (yv[nend-1]-yv[nend-2])/(xv[nend-1]-xv[nend-2]);
} else {
return (yv[lb_ix] - yv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]);
}
}
template <typename T>
T linearInterpolationExtrap(const std::vector<double>& xv,
const std::vector<T>& yv,
double x)
{
std::vector<double>::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x);
int lb_ix = lb - xv.begin();
if (lb_ix == 0) {
lb_ix=1;
} else if (lb_ix == int(xv.size())){
lb_ix = int(xv.size())-1;
}
double w = (x - xv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]);
return (1.0 - w)*yv[lb_ix - 1] + w*yv[lb_ix];
}
} // namespace Opm