blackoil: implement the immiscible solvent extension

Conceptually this is IMO pretty questionable, since it adds a second
"gas phase" that does not mix with "ordinary" gas. I suppose the
reason why this extension was conceived by E100 is that if all you
have is hammer, everything looks like a nail...

Functionality-wise, this patch is still not fully complete because
miscibility of the solvent "phase" is not yet implemented. As far as I
can see, the API changes required by miscibility are quite limited,
though.

ebos/eclfluxmodule.hh

@@ -32,6 +32,7 @@
 #define EWOMS_ECL_FLUX_MODULE_HH
 
 #include <ewoms/disc/common/fvbaseproperties.hh>
+#include <ewoms/models/blackoil/blackoilproperties.hh>
 #include <ewoms/common/signum.hh>
 
 #include <opm/common/Valgrind.hpp>
@@ -102,6 +103,8 @@ protected:
 template <class TypeTag>
 class EclTransExtensiveQuantities
 {
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
+
     typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
     typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
@@ -110,7 +113,9 @@ class EclTransExtensiveQuantities
     typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
 
     enum { dimWorld = GridView::dimensionworld };
-    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { numPhases = FluidSystem::numPhases };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
 
     typedef Opm::MathToolbox<Evaluation> Toolbox;
     typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
@@ -201,6 +206,9 @@ protected:
         return dnIdx_[phaseIdx];
     }
 
+    void updateSolvent(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
+    { asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx); }
+
     /*!
      * \brief Update the required gradients for interior faces
      */
@@ -347,6 +355,13 @@ protected:
     void calculateFluxes_(const ElementContext& elemCtx OPM_UNUSED, unsigned scvfIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
     { }
 
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+
     // the volumetric flux of all phases [m^3/s]
     Evaluation volumeFlux_[numPhases];
 

ebos/eclpeacemanwell.hh

@@ -28,6 +28,7 @@
 #ifndef EWOMS_ECL_PEACEMAN_WELL_HH
 #define EWOMS_ECL_PEACEMAN_WELL_HH
 
+#include <ewoms/models/blackoil/blackoilproperties.hh>
 #include <ewoms/aux/baseauxiliarymodule.hh>
 #include <ewoms/common/propertysystem.hh>
 #include <ewoms/common/alignedallocator.hh>
@@ -46,17 +47,6 @@
 #include <map>
 
 namespace Ewoms {
-namespace Properties {
-NEW_PROP_TAG(Scalar);
-NEW_PROP_TAG(Discretization);
-NEW_PROP_TAG(FluidSystem);
-NEW_PROP_TAG(Simulator);
-NEW_PROP_TAG(ElementContext);
-NEW_PROP_TAG(RateVector);
-NEW_PROP_TAG(GridView);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(NumComponents);
-}
 
 template <class TypeTag>
 class EcfvDiscretization;
@@ -124,6 +114,7 @@ class EclPeacemanWell : public BaseAuxiliaryModule<TypeTag>
     static const unsigned gasCompIdx = FluidSystem::gasCompIdx;
 
     static const unsigned numModelEq = GET_PROP_VALUE(TypeTag, NumEq);
+    static const unsigned conti0EqIdx = GET_PROP_TYPE(TypeTag, Indices)::conti0EqIdx;
 
     typedef Opm::CompositionalFluidState<Scalar, FluidSystem, /*storeEnthalpy=*/false> FluidState;
     typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
@@ -446,10 +437,9 @@ public:
 
             // now we put this derivative into the right place in the Jacobian
             // matrix. This is a bit hacky because it assumes that the model uses a mass
-            // rate for each component as its first conservation equations, but we
+            // rate for each component as its first conservation equation, but we require
-            // require the black-oil model for now anyway, so this should not be too much
+            // the black-oil model for now anyway, so this should not be too much of a
-            // of a problem...
+            // problem...
-            assert(numModelEq == numComponents);
             Opm::Valgrind::CheckDefined(q);
             auto& matrixEntry = matrix[gridDofIdx][wellGlobalDofIdx];
             matrixEntry = 0.0;
@@ -1127,8 +1117,8 @@ public:
                 continue;
 
             modelRate.setVolumetricRate(intQuants.fluidState(), phaseIdx, volumetricRates[phaseIdx]);
-            for (unsigned eqIdx = 0; eqIdx < q.size(); ++eqIdx)
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
-                q[eqIdx] += modelRate[eqIdx];
+                q[conti0EqIdx + compIdx] += modelRate[conti0EqIdx + compIdx];
         }
 
         Opm::Valgrind::CheckDefined(q);

ebos/eclproblem.hh

@@ -248,6 +248,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
     enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
     enum { numPhases = FluidSystem::numPhases };
     enum { numComponents = FluidSystem::numComponents };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
     enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
     enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
@@ -267,6 +268,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
     typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
     typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
 
+    typedef BlackOilSolventModule<TypeTag> SolventModule;
     typedef Opm::CompositionalFluidState<Scalar, FluidSystem> ScalarFluidState;
     typedef Opm::MathToolbox<Evaluation> Toolbox;
     typedef Ewoms::EclSummaryWriter<TypeTag> EclSummaryWriter;
@@ -315,6 +317,10 @@ public:
     {
         // add the output module for the Ecl binary output
         simulator.model().addOutputModule(new Ewoms::EclOutputBlackOilModule<TypeTag>(simulator));
+
+        // Tell the solvent module to initialize its internal data structures
+        const auto& gridManager = simulator.gridManager();
+        SolventModule::initFromDeck(gridManager.deck(), gridManager.eclState());
     }
 
     /*!
@@ -754,6 +760,24 @@ public:
         return pvtnum_[elemIdx];
     }
 
+    /*!
+     * \brief Returns the index of the relevant region for thermodynmic properties
+     */
+    template <class Context>
+    unsigned satnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
+    { return satnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
+
+    /*!
+     * \brief Returns the index the relevant saturation function region given a cell index
+     */
+    unsigned satnumRegionIndex(unsigned elemIdx) const
+    {
+        if (satnum_.empty())
+            return 0;
+
+        return satnum_[elemIdx];
+    }
+
     /*!
      * \copydoc FvBaseProblem::name
      */
@@ -803,6 +827,9 @@ public:
         }
         else
             values.assignNaive(initialFluidStates_[globalDofIdx]);
+
+        //if (enableSolvent)
+        //    values[Indices::solventSaturationIdx] = whatever;
     }
 
     /*!
@@ -1004,8 +1031,9 @@ private:
         const auto& deck = gridManager.deck();
         const auto& eclState = gridManager.eclState();
 
-        // the PVT region number
+        // the PVT and saturation region numbers
         updatePvtnum_();
+        updateSatnum_();
 
         ////////////////////////////////
         // porosity
@@ -1357,6 +1385,25 @@ private:
         }
     }
 
+    void updateSatnum_()
+    {
+        const auto& eclState = this->simulator().gridManager().eclState();
+        const auto& eclProps = eclState.get3DProperties();
+
+        if (!eclProps.hasDeckIntGridProperty("SATNUM"))
+            return;
+
+        const auto& satnumData = eclProps.getIntGridProperty("SATNUM").getData();
+        const auto& gridManager = this->simulator().gridManager();
+
+        unsigned numElems = gridManager.gridView().size(/*codim=*/0);
+        satnum_.resize(numElems);
+        for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
+            unsigned cartElemIdx = gridManager.cartesianIndex(elemIdx);
+            satnum_[elemIdx] = satnumData[cartElemIdx] - 1;
+        }
+    }
+
     struct PffDofData_
     {
         Scalar transmissibility;
@@ -1401,6 +1448,7 @@ private:
     EclThresholdPressure<TypeTag> thresholdPressures_;
 
     std::vector<unsigned short> pvtnum_;
+    std::vector<unsigned short> satnum_;
     std::vector<unsigned short> rockTableIdx_;
     std::vector<RockParams> rockParams_;