compute rock fraction

ebos/eclgenericproblem.cc

@@ -347,6 +347,20 @@ porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
     return this->referencePorosity_[timeIdx][globalSpaceIdx];
 }
 
+template<class GridView, class FluidSystem, class Scalar>
+Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
+rockFraction(unsigned elementIdx, unsigned timeIdx) const
+{
+    const auto& fp = eclState_.fieldProps();
+    const std::vector<double>& poro  = fp.get_double("PORO");
+    // the reference porosity is defined as the accumulated pore volume divided by the
+    // geometric volume of the element. Note that it can
+    // be larger than 1.0 if porevolume multipliers are used
+    // to for instance implement larger boundary cells
+    Scalar porosity = poro[elementIdx];
+    return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
+}
+
 template<class GridView, class FluidSystem, class Scalar>
 template<class T>
 void EclGenericProblem<GridView,FluidSystem,Scalar>::

ebos/eclgenericproblem.hh

@@ -138,6 +138,17 @@ public:
     Scalar referencePorosity(unsigned elementIdx, unsigned timeIdx) const
     { return referencePorosity_[timeIdx][elementIdx]; }
 
+
+    /*!
+     * \brief Returns the rockFraction of an element
+     *
+     * Usually 1 - porosity, but if pvmult is used to modify porosity
+     * we will apply the same multiplier to the rock fraction
+     * i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume
+     * multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
+     */
+    Scalar rockFraction(unsigned elementIdx, unsigned timeIdx) const;
+
     /*!
      * \brief Sets the porosity of an element
      *