Merge pull request #6181 from daavid00/summaryMassWat

Summary vectors for water mass

opm/simulators/flow/FIPContainer.cpp

@@ -422,6 +422,26 @@ assignCalciteMass(const unsigned globalDofIdx,
     this->fip_[Inplace::Phase::CalciteMass][globalDofIdx] = calciteMass;
 }
 
+template<class FluidSystem>
+bool
+FIPContainer<FluidSystem>::
+hasWaterMass() const
+{
+    return has(Inplace::Phase::WaterMass);
+}
+
+template<class FluidSystem>
+void
+FIPContainer<FluidSystem>::
+assignWaterMass(const unsigned globalDofIdx,
+                const std::array<Scalar, numPhases>& fip,
+                const Scalar   rhoW)
+{
+    if (this->has(Inplace::Phase::WaterMass)) {
+        this->fip_[Inplace::Phase::WaterMass][globalDofIdx] = fip[waterPhaseIdx] * rhoW;
+    }
+}
+
 template<class FluidSystem>
 void
 FIPContainer<FluidSystem>::

opm/simulators/flow/FIPContainer.hpp

@@ -106,6 +106,11 @@ public:
     void assignCalciteMass(const unsigned globalDofIdx,
                            const Scalar   calciteMass);
 
+    bool hasWaterMass() const;
+    void assignWaterMass(const unsigned globalDofIdx,
+                         const std::array<Scalar, numPhases>& fip,
+                         const Scalar   rhoW);
+
     void assignGasWater(const unsigned  globalDofIdx,
                         const std::array<Scalar, numPhases>& fip,
                         const Scalar    gasInPlaceWater,

opm/simulators/flow/LogOutputHelper.cpp

@@ -1504,6 +1504,7 @@ fipUnitConvert_(std::unordered_map<Inplace::Phase, Scalar>& fip) const
         {Inplace::Phase::CO2InGasPhaseMobKrg,       M::moles},
         {Inplace::Phase::WaterInWaterPhase,         M::liquid_surface_volume},
         {Inplace::Phase::WaterInGasPhase,           M::liquid_surface_volume},
+        {Inplace::Phase::WaterMass,                 M::mass},
         {Inplace::Phase::CO2Mass,                   M::mass},
         {Inplace::Phase::CO2MassInWaterPhase,       M::mass},
         {Inplace::Phase::CO2MassInGasPhase,         M::mass},

opm/simulators/flow/OutputBlackoilModule.hpp

@@ -893,6 +893,11 @@ private:
                 this->updateCalciteMass(globalDofIdx, intQuants, totVolume);
             }
         }
+
+        if (this->fipC_.hasWaterMass() && FluidSystem::phaseIsActive(waterPhaseIdx))
+        {
+            this->updateWaterMass(globalDofIdx, fs, fip);
+        }
     }
 
     template <typename FluidState, typename FIPArray>
@@ -1011,6 +1016,17 @@ private:
         return xoG * pv * rhoo * so / mM;
     }
 
+    template <typename FluidState, typename FIPArray>
+    void updateWaterMass(const unsigned    globalDofIdx,
+                         const FluidState& fs,
+                         const FIPArray&   fip
+                         )
+    {
+        const Scalar rhoW = FluidSystem::referenceDensity(waterPhaseIdx, fs.pvtRegionIndex());
+
+        this->fipC_.assignWaterMass(globalDofIdx, fip, rhoW);
+    }   
+
     template <typename IntensiveQuantities>
     void updateMicrobialMass(const unsigned             globalDofIdx,
                              const IntensiveQuantities& intQuants,
@@ -1907,8 +1923,8 @@ private:
                                           getValue(ectx.fs.Rv()) *
                                           getValue(ectx.fs.invB(gasPhaseIdx)) *
                                           getValue(ectx.fs.saturation(gasPhaseIdx))) *
-                                         model.dofTotalVolume(ectx.globalDofIdx) *
-                                         getValue(ectx.intQuants.porosity());
+                                          model.dofTotalVolume(ectx.globalDofIdx) *
+                                          getValue(ectx.intQuants.porosity());
                               }
                   }
             },
@@ -2003,8 +2019,8 @@ private:
                               &problem = this->simulator_.problem(),
                               &regions = this->regions_](const unsigned phaseIdx, const Context& ectx)
                              {
-                                 auto phase = RegionPhasePoreVolAverage::Phase{};
-                                 phase.ix = phaseIdx;
+                                auto phase = RegionPhasePoreVolAverage::Phase{};
+                                phase.ix = phaseIdx;
 
                                 // Note different region handling here.  FIPNUM is
                                 // one-based, but we need zero-based lookup in
@@ -2028,6 +2044,19 @@ private:
                             }
                   }
             },
+            Entry{ScalarEntry{"BAMIP",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  const Scalar rhoW = FluidSystem::referenceDensity(waterPhaseIdx,
+                                                                   ectx.intQuants.pvtRegionIndex());
+                                  return getValue(ectx.fs.invB(waterPhaseIdx)) *
+                                         getValue(ectx.fs.saturation(waterPhaseIdx)) *
+                                         rhoW *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity());
+                              }
+                  }
+            },
             Entry{ScalarEntry{"BMMIP",
                               [&model = this->simulator_.model()](const Context& ectx)
                               {

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1842,7 +1842,6 @@ assignMassGasRate(data::Wells& wsrpt,
     }
 }
 
-
 template<class Scalar>
 void BlackoilWellModelGeneric<Scalar>::
 assignWellTracerRates(data::Wells& wsrpt,

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1764,6 +1764,11 @@ namespace Opm
         const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
 
         const int nseg = this->numberOfSegments();
+        
+        const Scalar rhow = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ?
+            FluidSystem::referenceDensity( FluidSystem::waterPhaseIdx, Base::pvtRegionIdx() ) : 0.0;
+        const unsigned watCompIdx = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ?
+            Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) : 0;
 
         for (int seg = 0; seg < nseg; ++seg) {
             // calculating the perforation rate for each perforation that belongs to this segment
@@ -1802,6 +1807,11 @@ namespace Opm
                 }
                 perf_press_state[local_perf_index] = perf_press.value();
 
+                // mass rates, for now only water
+                if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                    perf_data.wat_mass_rates[local_perf_index] = cq_s[watCompIdx].value() * rhow;
+                }
+
                 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
                     // the cq_s entering mass balance equations need to consider the efficiency factors.
                     const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;

opm/simulators/wells/PerfData.cpp

@@ -52,6 +52,7 @@ PerfData<Scalar>::PerfData(const std::size_t num_perf,
     , satnum_id(num_perf)
     , ecl_index(num_perf)
     , gas_mass_rates(num_perf)
+    , wat_mass_rates(num_perf)
 {
     if (injector) {
         prepareInjectorContainers();
@@ -97,6 +98,7 @@ PerfData<Scalar> PerfData<Scalar>::serializationTestObject()
     result.filtrate_data = ConnFiltrateData<Scalar>::serializationTestObject();
     result.connFracStatistics.assign(3, ConnFracStatistics<Scalar>::serializationTestObject());
     result.gas_mass_rates = {31.0};
+    result.wat_mass_rates = {32.0};
 
     return result;
 }
@@ -142,6 +144,7 @@ bool PerfData<Scalar>::try_assign(const PerfData& other)
     this->filtrate_data = other.filtrate_data;
     this->connFracStatistics = other.connFracStatistics;
     this->gas_mass_rates = other.gas_mass_rates;
+    this->wat_mass_rates = other.wat_mass_rates;
 
     return true;
 }
@@ -174,6 +177,7 @@ bool PerfData<Scalar>::operator==(const PerfData& rhs) const
         && (this->filtrate_data == rhs.filtrate_data)
         && (this->connFracStatistics == rhs.connFracStatistics)
         && (this->gas_mass_rates == rhs.gas_mass_rates)
+        && (this->wat_mass_rates == rhs.wat_mass_rates)
         ;
 }
 

opm/simulators/wells/PerfData.hpp

@@ -81,6 +81,7 @@ public:
         serializer(filtrate_data);
         serializer(connFracStatistics);
         serializer(gas_mass_rates);
+        serializer(wat_mass_rates);
     }
 
     bool operator==(const PerfData&) const;
@@ -110,6 +111,7 @@ public:
     std::vector<int> satnum_id{};
     std::vector<std::size_t> ecl_index{};
     std::vector<Scalar> gas_mass_rates{};
+    std::vector<Scalar> wat_mass_rates{};
 
     // The water_throughput, skin_pressure and water_velocity variables are
     // only used for injectors to check the injectivity.

opm/simulators/wells/SingleWellState.cpp

@@ -211,6 +211,12 @@ Scalar SingleWellState<Scalar>::sum_urea_rates() const
     return this->sum_connection_rates(this->perf_data.urea_rates);
 }
 
+template<class Scalar>
+Scalar SingleWellState<Scalar>::sum_wat_mass_rates() const
+{
+    return this->sum_connection_rates(this->perf_data.wat_mass_rates);
+}
+
 template<class Scalar>
 Scalar SingleWellState<Scalar>::sum_filtrate_rate() const
 {

opm/simulators/wells/SingleWellState.hpp

@@ -147,6 +147,7 @@ public:
     Scalar sum_microbial_rates() const;
     Scalar sum_oxygen_rates() const;
     Scalar sum_urea_rates() const;
+    Scalar sum_wat_mass_rates() const;
 
     Scalar sum_filtrate_rate() const;
     Scalar sum_filtrate_total() const;

opm/simulators/wells/StandardWell_impl.hpp

@@ -644,6 +644,13 @@ namespace Opm
             const Scalar rho = FluidSystem::referenceDensity( FluidSystem::gasPhaseIdx, Base::pvtRegionIdx() );
             perf_data.gas_mass_rates[perf] = cq_s[gas_comp_idx].value() * rho;
         }
+
+        // Store the perforation water mass rate.
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            const unsigned wat_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+            const Scalar rho = FluidSystem::referenceDensity( FluidSystem::waterPhaseIdx, Base::pvtRegionIdx() );
+            perf_data.wat_mass_rates[perf] = cq_s[wat_comp_idx].value() * rho;
+        }
     }
 
 

opm/simulators/wells/WellState.cpp

@@ -549,6 +549,7 @@ WellState<Scalar>::report(const int* globalCellIdxMap,
             well.rates.set(rt::reservoir_water, reservoir_rates[pu.phase_pos[BlackoilPhases::Aqua]]);
             well.rates.set(rt::productivity_index_water, wpi[pu.phase_pos[BlackoilPhases::Aqua]]);
             well.rates.set(rt::well_potential_water, well_potentials[pu.phase_pos[BlackoilPhases::Aqua]]);
+            well.rates.set(rt::mass_wat, ws.sum_wat_mass_rates());
         }
 
         if (pu.phase_used[BlackoilPhases::Liquid]) {
@@ -1159,6 +1160,12 @@ reportConnectionPressuresAndRates(const std::size_t well_index,
             connections[i].rates.set(rt::mass_gas, perf_data.gas_mass_rates[i]);
         }
     }
+
+    if (pu.phase_used[Water]) {
+        for (auto i = 0*num_perf_well; i < num_perf_well; ++i) {
+            connections[i].rates.set(rt::mass_wat, perf_data.wat_mass_rates[i]);
+        }
+    }
 }
 
 template <class Scalar>