Finalize integration of foam module.

opm/core/props/BlackoilPhases.hpp

@@ -33,11 +33,11 @@ namespace Opm
         // sense that they can be active or not and canonical indices can be translated
         // to and from active ones. That said, they are not considered by num_phases or
         // MaxNumPhases. The crypto phases which are currently implemented are solvent,
-        // polymer, energy and polymer molecular weight.
+        // polymer, energy, polymer molecular weight and foam.
-        static const int NumCryptoPhases = 4;
+        static const int NumCryptoPhases = 5;
 
         // enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 };
-        enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Solvent = 3, Polymer = 4, Energy = 5, PolymerMW = 6 };
+        enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Solvent = 3, Polymer = 4, Energy = 5, PolymerMW = 6, Foam = 7 };
     };
 
     struct PhaseUsage : public BlackoilPhases
@@ -50,6 +50,7 @@ namespace Opm
         bool has_energy;
         // polymer molecular weight
         bool has_polymermw;
+        bool has_foam;
     };
 
     /// Check or assign presence of a formed, free phase.  Limited to

opm/core/props/phaseUsageFromDeck.hpp

@@ -119,6 +119,15 @@ namespace Opm
         else
             pu.phase_pos[BlackoilPhases::PolymerMW] = -1;
 
+        // Add foam info
+        pu.has_foam = phases.active(Phase::FOAM);
+        if (pu.has_foam) {
+            pu.phase_pos[BlackoilPhases::Foam] = numActivePhases;
+            ++ numActivePhases;
+        }
+        else
+            pu.phase_pos[BlackoilPhases::Foam] = -1;
+
         return pu;
     }
 

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -83,6 +83,7 @@ SET_BOOL_PROP(EclFlowProblem, EnablePolymer, false);
 SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false);
 SET_BOOL_PROP(EclFlowProblem, EnableTemperature, true);
 SET_BOOL_PROP(EclFlowProblem, EnableEnergy, false);
+SET_BOOL_PROP(EclFlowProblem, EnableFoam, false);
 
 SET_TYPE_PROP(EclFlowProblem, EclWellModel, Opm::BlackoilWellModel<TypeTag>);
 SET_TAG_PROP(EclFlowProblem, LinearSolverSplice, FlowIstlSolver);
@@ -122,10 +123,12 @@ namespace Opm {
         static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
         static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
         static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
+        static const int contiFoamEqIdx = Indices::contiFoamEqIdx;
         static const int solventSaturationIdx = Indices::solventSaturationIdx;
         static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
         static const int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
         static const int temperatureIdx = Indices::temperatureIdx;
+        static const int foamConcentrationIdx = Indices::foamConcentrationIdx;
 
         typedef Dune::FieldVector<Scalar, numEq >        VectorBlockType;
         typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
@@ -160,6 +163,7 @@ namespace Opm {
         , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
         , has_polymermw_(GET_PROP_VALUE(TypeTag, EnablePolymerMW))
         , has_energy_(GET_PROP_VALUE(TypeTag, EnableEnergy))
+        , has_foam_(GET_PROP_VALUE(TypeTag, EnableFoam))
         , param_( param )
         , well_model_ (well_model)
         , terminal_output_ (terminal_output)
@@ -620,6 +624,12 @@ namespace Opm {
                     R_sum[ contiPolymerEqIdx ] += R2;
                     maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
                 }
+                if (has_foam_ ) {
+                    B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
+                    const auto R2 = ebosResid[cell_idx][contiFoamEqIdx];
+                    R_sum[ contiFoamEqIdx ] += R2;
+                    maxCoeff[ contiFoamEqIdx ] = std::max( maxCoeff[ contiFoamEqIdx ], std::abs( R2 ) / pvValue );
+                }
 
                 if (has_polymermw_) {
                     assert(has_polymer_);
@@ -706,6 +716,9 @@ namespace Opm {
                 if (has_energy_) {
                     compNames[temperatureIdx] = "Energy";
                 }
+                if (has_foam_) {
+                    compNames[foamConcentrationIdx] = "Foam";
+                }
             }
 
             // Create convergence report.
@@ -857,6 +870,7 @@ namespace Opm {
         const bool has_polymer_;
         const bool has_polymermw_;
         const bool has_energy_;
+        const bool has_foam_;
 
         ModelParameters                 param_;
         SimulatorReport failureReport_;

opm/simulators/wells/StandardWell.hpp

@@ -28,6 +28,10 @@
 #include <opm/simulators/wells/WellInterface.hpp>
 #include <opm/simulators/linalg/ISTLSolverEbos.hpp>
 
+#include <ewoms/models/blackoil/blackoilpolymermodules.hh>
+#include <ewoms/models/blackoil/blackoilsolventmodules.hh>
+#include <ewoms/models/blackoil/blackoilfoammodules.hh>
+
 #include <opm/material/densead/DynamicEvaluation.hpp>
 
 #include <dune/common/dynvector.hh>
@@ -53,8 +57,6 @@ namespace Opm
         using typename Base::MaterialLaw;
         using typename Base::ModelParameters;
         using typename Base::Indices;
-        using typename Base::PolymerModule;
-        using typename Base::FoamModule;
         using typename Base::RateConverterType;
 
         using Base::numEq;
@@ -64,6 +66,9 @@ namespace Opm
         using Base::has_foam;
         using Base::has_energy;
 
+        using PolymerModule =  Ewoms::BlackOilPolymerModule<TypeTag>;
+        using FoamModule = Ewoms::BlackOilFoamModule<TypeTag>;
+
         // polymer concentration and temperature are already known by the well, so
         // polymer and energy conservation do not need to be considered explicitly
         static const int numPolymerEq = Indices::numPolymers;
@@ -129,6 +134,7 @@ namespace Opm
 
         using Base::contiSolventEqIdx;
         using Base::contiPolymerEqIdx;
+        using Base::contiFoamEqIdx;
         static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
 
         StandardWell(const Well2& well, const int time_step, const Wells* wells,
@@ -202,6 +208,7 @@ namespace Opm
         using Base::ebosCompIdxToFlowCompIdx;
         using Base::wsolvent;
         using Base::wpolymer;
+        using Base::wfoam;
         using Base::wellHasTHPConstraints;
         using Base::mostStrictBhpFromBhpLimits;
         using Base::scalingFactor;

opm/simulators/wells/StandardWell_impl.hpp

@@ -643,6 +643,18 @@ namespace Opm
                 }
             }
 
+            if (has_foam) {
+                // TODO: the application of well efficiency factor has not been tested with an example yet
+                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+                EvalWell cq_s_foam = cq_s[gasCompIdx] * well_efficiency_factor_;
+                if (well_type_ == INJECTOR) {
+                    cq_s_foam *= wfoam();
+                } else {
+                    cq_s_foam *= extendEval(intQuants.foamConcentration());
+                }
+                connectionRates_[perf][contiFoamEqIdx] = Base::restrictEval(cq_s_foam);
+            }
+
             // Store the perforation pressure for later usage.
             well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf];
 
@@ -1973,7 +1985,7 @@ namespace Opm
                        Opm::DeferredLogger& deferred_logger) const
     {
         // the following implementation assume that the polymer is always after the w-o-g phases
-        // For the polymer case and the energy case, there is one more mass balance equations of reservoir than wells
+        // For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells
         assert((int(B_avg.size()) == num_components_) || has_polymer || has_energy || has_foam);
 
         const double tol_wells = param_.tolerance_wells_;

opm/simulators/wells/WellInterface.hpp

@@ -47,10 +47,6 @@
 #include <opm/simulators/timestepping/ConvergenceReport.hpp>
 #include <opm/simulators/utils/DeferredLogger.hpp>
 
-#include <ewoms/models/blackoil/blackoilpolymermodules.hh>
-#include <ewoms/models/blackoil/blackoilsolventmodules.hh>
-#include <ewoms/models/blackoil/blackoilfoammodules.hh>
-
 #include<dune/common/fmatrix.hh>
 #include<dune/istl/bcrsmatrix.hh>
 #include<dune/istl/matrixmatrix.hh>
@@ -98,22 +94,19 @@ namespace Opm
         typedef Dune::BlockVector<VectorBlockType> BVector;
         typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
 
-        typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
-        typedef Ewoms::BlackOilFoamModule<TypeTag> FoamModule;
-
         static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
         static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
         static const bool has_energy = GET_PROP_VALUE(TypeTag, EnableEnergy);
         // flag for polymer molecular weight related
         static const bool has_polymermw = GET_PROP_VALUE(TypeTag, EnablePolymerMW);
+        static const bool has_foam = GET_PROP_VALUE(TypeTag, EnableFoam);
         static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
         static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
         // index for the polymer molecular weight continuity equation
         static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
+        static const int contiFoamEqIdx = Indices::contiFoamEqIdx;
 
-        static const bool has_foam = GET_PROP_VALUE(TypeTag, EnableFoam);
+        // For the conversion between the surface volume rate and reservoir voidage rate
-
-        // For the conversion between the surface volume rate and resrevoir voidage rate
         using RateConverterType = RateConverter::
         SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
 
@@ -370,6 +363,8 @@ namespace Opm
 
         double wpolymer() const;
 
+        double wfoam() const;
+
         bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
                                  const WellState& well_state,
                                  Opm::DeferredLogger& deferred_logger) const;

opm/simulators/wells/WellInterface_impl.hpp

@@ -350,6 +350,30 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    double
+    WellInterface<TypeTag>::
+    wfoam() const
+    {
+        if (!has_foam) {
+            return 0.0;
+        }
+
+        auto injectorType = well_ecl_.injectorType();
+
+        if (injectorType == WellInjector::GAS) {
+            WellFoamProperties fprop = well_ecl_.getFoamProperties();
+            return fprop.m_foamConcentration;
+        } else {
+            // Not a gas injection well => no foam.
+            return 0.0;
+        }
+    }
+
+
+
+
+
     template<typename TypeTag>
     double
     WellInterface<TypeTag>::