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

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 poromult = Opm::linearInterpolation(p_, poromult_, pressure);
-            const double dporomultdp = Opm::linearInterpolationDerivative(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;
         }
     }

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){

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)  )

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;

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");
+    LinearSolverAGMG::LinearSolverAGMG(const parameter::ParameterGroup& param)
-#endif  // HAVE_AGMG
+        : 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    ,

opm/core/linalg/LinearSolverAGMG.hpp

@@ -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();
 

opm/core/linalg/LinearSolverFactory.cpp

@@ -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 if (ls == "agmg") {
-        else {
+#if HAVE_AGMG
+            solver_.reset(new LinearSolverAGMG(param));
+#endif
+       }
+       else {
             THROW("Linear solver " << ls << " is unknown.");
         }
 

opm/core/utility/linearInterpolation.hpp

@@ -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