Merge pull request #4791 from svenn-t/aquifer_h2store

Enable aquifers with H2STORE

ebos/equil/initstateequil.hh

@@ -715,7 +715,7 @@ private:
 
     void applyNumericalAquifers_(const GridView& gridView,
                                  const NumericalAquifers& aquifer,
-                                 const bool co2store);
+                                 const bool co2store_or_h2store);
 
     template<class RMap>
     void setRegionPvtIdx(const EclipseState& eclState, const RMap& reg);

ebos/equil/initstateequil_impl.hh

@@ -1507,7 +1507,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
     calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
 
     // modify the pressure and saturation for numerical aquifer cells
-    applyNumericalAquifers_(gridView, num_aquifers, eclipseState.runspec().co2Storage());
+    applyNumericalAquifers_(gridView, num_aquifers, eclipseState.runspec().co2Storage() || eclipseState.runspec().h2Storage());
 
     // Modify oil pressure in no-oil regions so that the pressures of present phases can
     // be recovered from the oil pressure and capillary relations.
@@ -1654,14 +1654,23 @@ void InitialStateComputer<FluidSystem,
                           CartesianIndexMapper>::
 applyNumericalAquifers_(const GridView& gridView,
                         const NumericalAquifers& aquifer,
-                        const bool co2store)
+                        const bool co2store_or_h2store)
 {
     const auto num_aqu_cells = aquifer.allAquiferCells();
     if (num_aqu_cells.empty()) return;
 
+    // Check if water phase is active, or in the case of CO2STORE and H2STORE, water is modelled as oil phase
+    bool oil_as_brine = co2store_or_h2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
+    const auto watPos =  oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx;
+    if (!FluidSystem::phaseIsActive(watPos)){
+        throw std::logic_error  { "Water phase has to be active for numerical aquifer case" };
+    }
+
     ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
     auto elemIt = gridView.template begin</*codim=*/0>();
     const auto& elemEndIt = gridView.template end</*codim=*/0>();
+    const auto oilPos = FluidSystem::oilPhaseIdx;
+    const auto gasPos = FluidSystem::gasPhaseIdx;
     for (; elemIt != elemEndIt; ++elemIt) {
         const Element& element = *elemIt;
         const unsigned int elemIdx = elemMapper.index(element);
@@ -1669,21 +1678,12 @@ applyNumericalAquifers_(const GridView& gridView,
         const auto search = num_aqu_cells.find(cartIx);
         if (search != num_aqu_cells.end()) {
             // numerical aquifer cells are filled with water initially
-            // for co2store the oilphase may be used for the brine
+            this->sat_[watPos][elemIdx] = 1.;
-            bool co2store_oil_as_brine = co2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
-            const auto watPos =  co2store_oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx;
-            if (FluidSystem::phaseIsActive(watPos)) {
-                this->sat_[watPos][elemIdx] = 1.;
-            } else {
-                throw std::logic_error  { "Water phase has to be active for numerical aquifer case" };
-            }
 
-            const auto oilPos = FluidSystem::oilPhaseIdx;
+            if (!co2store_or_h2store && FluidSystem::phaseIsActive(oilPos)) {
-            if (!co2store && FluidSystem::phaseIsActive(oilPos)) {
                 this->sat_[oilPos][elemIdx] = 0.;
             }
 
-            const auto gasPos = FluidSystem::gasPhaseIdx;
             if (FluidSystem::phaseIsActive(gasPos)) {
                 this->sat_[gasPos][elemIdx] = 0.;
             }

opm/simulators/aquifers/AquiferAnalytical.hpp

@@ -246,7 +246,7 @@ protected:
 
     int compIdx_() const
     {
-        if (this->co2store_())
+        if (this->co2store_or_h2store_())
             return FluidSystem::oilCompIdx;
 
         return FluidSystem::waterCompIdx;

opm/simulators/aquifers/AquiferCarterTracy.hpp

@@ -104,7 +104,7 @@ public:
         aquCT->dimensionless_pressure = this->dimensionless_pressure_;
         aquCT->influxConstant = this->aquct_data_.influxConstant();
 
-        if (!this->co2store_()) {
+        if (!this->co2store_or_h2store_()) {
             aquCT->timeConstant = this->aquct_data_.timeConstant();
             aquCT->waterDensity = this->aquct_data_.waterDensity();
             aquCT->waterViscosity = this->aquct_data_.waterViscosity();
@@ -217,7 +217,7 @@ protected:
 
     void calculateAquiferConstants() override
     {
-        this->Tc_ = this->co2store_()
+        this->Tc_ = this->co2store_or_h2store_()
             ? this->timeConstantCO2Store()
             : this->aquct_data_.timeConstant();
 
@@ -245,7 +245,7 @@ protected:
             this->Ta0_ = this->aquct_data_.initial_temperature.value();
         }
 
-        this->rhow_ = this->co2store_()
+        this->rhow_ = this->co2store_or_h2store_()
             ? this->waterDensityCO2Store()
             : this->aquct_data_.waterDensity();
     }

opm/simulators/aquifers/AquiferConstantFlux.hpp

@@ -214,7 +214,7 @@ private:
     // TODO: this is a function from AquiferAnalytical
     int compIdx_() const
     {
-        if (this->co2store_())
+        if (this->co2store_or_h2store_())
             return FluidSystem::oilCompIdx;
 
         return FluidSystem::waterCompIdx;

opm/simulators/aquifers/AquiferInterface.hpp

@@ -78,15 +78,16 @@ public:
     int aquiferID() const { return this->aquiferID_; }
 
 protected:
-    bool co2store_() const
+    bool co2store_or_h2store_() const
     {
-        return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
+        const auto& rspec = ebos_simulator_.vanguard().eclState().runspec();
+        return rspec.co2Storage() || rspec.h2Storage();
     }
 
     int phaseIdx_() const
     {
         // If OIL is used to model brine the aquifer should do the same
-        if (co2store_() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
+        if (co2store_or_h2store_() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
             return FluidSystem::oilPhaseIdx;
 
         return FluidSystem::waterPhaseIdx;