From 064184bd0c249dc7a2a35f1b805112f865662402 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Tue, 14 Aug 2012 09:36:19 +0200 Subject: [PATCH] Changed to extrapolation outside tables for pormult. Added facetags to the grid structure. Changed default fluid to Linear. --- opm/core/fluid/RockCompressibility.cpp | 10 ++++--- opm/core/fluid/SaturationPropsBasic.cpp | 3 ++- opm/core/grid.c | 3 +++ opm/core/grid/cart_grid.c | 13 +++++++++ opm/core/linalg/LinearSolverAGMG.cpp | 15 ++++++++--- opm/core/linalg/LinearSolverAGMG.hpp | 9 ++++++- opm/core/linalg/LinearSolverFactory.cpp | 12 +++++++-- opm/core/utility/linearInterpolation.hpp | 34 ++++++++++++++++++++++++ 8 files changed, 89 insertions(+), 10 deletions(-) diff --git a/opm/core/fluid/RockCompressibility.cpp b/opm/core/fluid/RockCompressibility.cpp index 17d86beb..10486101 100644 --- a/opm/core/fluid/RockCompressibility.cpp +++ b/opm/core/fluid/RockCompressibility.cpp @@ -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; } } diff --git a/opm/core/fluid/SaturationPropsBasic.cpp b/opm/core/fluid/SaturationPropsBasic.cpp index 8249a31f..760f934d 100644 --- a/opm/core/fluid/SaturationPropsBasic.cpp +++ b/opm/core/fluid/SaturationPropsBasic.cpp @@ -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){ diff --git a/opm/core/grid.c b/opm/core/grid.c index 60093e83..07250ea3 100644 --- a/opm/core/grid.c +++ b/opm/core/grid.c @@ -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) ) diff --git a/opm/core/grid/cart_grid.c b/opm/core/grid/cart_grid.c index c715aef3..25c04d24 100644 --- a/opm/core/grid/cart_grid.c +++ b/opm/core/grid/cart_grid.c @@ -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; diff --git a/opm/core/linalg/LinearSolverAGMG.cpp b/opm/core/linalg/LinearSolverAGMG.cpp index ae76bd4e..5ce7dfca 100644 --- a/opm/core/linalg/LinearSolverAGMG.cpp +++ b/opm/core/linalg/LinearSolverAGMG.cpp @@ -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 , diff --git a/opm/core/linalg/LinearSolverAGMG.hpp b/opm/core/linalg/LinearSolverAGMG.hpp index 813d57b8..f8c95ad7 100644 --- a/opm/core/linalg/LinearSolverAGMG.hpp +++ b/opm/core/linalg/LinearSolverAGMG.hpp @@ -43,6 +43,8 @@ */ #include +#include +#include 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(); diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp index 9b6d8909..98552c84 100644 --- a/opm/core/linalg/LinearSolverFactory.cpp +++ b/opm/core/linalg/LinearSolverFactory.cpp @@ -30,6 +30,10 @@ #if HAVE_DUNE_ISTL #include #endif +#if HAVE_AGMG +#include +#endif + #include #include @@ -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."); } diff --git a/opm/core/utility/linearInterpolation.hpp b/opm/core/utility/linearInterpolation.hpp index 0ed8b57b..743f396c 100644 --- a/opm/core/utility/linearInterpolation.hpp +++ b/opm/core/utility/linearInterpolation.hpp @@ -85,6 +85,40 @@ namespace Opm } } + template + T linearInterpolationDerivativeExtrap(const std::vector& xv, + const std::vector& yv, + double x) + { + std::vector::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 + T linearInterpolationExtrap(const std::vector& xv, + const std::vector& yv, + double x) + { + std::vector::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