Merge pull request #472 from totto82/freeBC

Add support for using free boundary conditions in Flow

ebos/eclequilinitializer.hh

@@ -146,11 +146,18 @@ public:
             if (enableTemperature || enableEnergy)
                 fluidState.setTemperature(initialState.temperature()[elemIdx]);
 
-            // set the phase pressures.
+            // set the phase pressures, invB factor and density
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                 if (!FluidSystem::phaseIsActive(phaseIdx))
                     continue;
                 fluidState.setPressure(phaseIdx, initialState.press()[phaseIdx][elemIdx]);
+
+                const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, regionIdx);
+                fluidState.setInvB(phaseIdx, b);
+
+                const auto& rho = FluidSystem::density(fluidState, phaseIdx, regionIdx);
+                fluidState.setDensity(phaseIdx, rho);
+
             }
         }
     }

ebos/eclfluxmodule.hh

@@ -361,11 +361,12 @@ protected:
                                      unsigned timeIdx,
                                      const FluidState& exFluidState)
     {
-        bool enableBoundaryMassFlux = false;
+        const auto& problem = elemCtx.problem();
+
+        bool enableBoundaryMassFlux = problem.hasFreeBoundaryConditions();
         if (!enableBoundaryMassFlux)
             return;
 
-        const auto& problem = elemCtx.problem();
         const auto& stencil = elemCtx.stencil(timeIdx);
         const auto& scvf = stencil.boundaryFace(scvfIdx);
 

ebos/eclproblem.hh

@@ -623,6 +623,8 @@ public:
         }
 
         tracerModel_.init();
+
+        readBoundaryConditions_();
     }
 
     void prefetch(const Element& elem) const
@@ -823,24 +825,6 @@ public:
         aquiferModel_.endTimeStep();
         tracerModel_.endTimeStep();
 
-
-        // we no longer need the initial soluiton
-        if (this->simulator().episodeIndex() == 0 && !initialFluidStates_.empty())  {
-            // we always need to provide a temperature and if energy is not conserved, we
-            // use the initial one. This means we have to "salvage" the temperature from
-            // the initial fluid states before deleting the array.
-            if (!enableEnergy) {
-                initialTemperature_.resize(initialFluidStates_.size());
-                for (unsigned i = 0; i < initialFluidStates_.size(); ++i) {
-                    const auto& fs = initialFluidStates_[i];
-                    initialTemperature_[i] = fs.temperature(/*phaseIdx=*/0);
-                }
-            }
-
-            initialFluidStates_.clear();
-        }
-
-
         updateCompositionChangeLimits_();
     }
 
@@ -1308,6 +1292,48 @@ public:
             values.setThermalFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
         }
 
+        if (hasFreeBoundaryConditions()) {
+
+            unsigned indexInInside  = context.intersection(spaceIdx).indexInInside();
+            unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
+            unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
+            switch (indexInInside) {
+            case 0:
+            {
+                if (freebcXMinus_[globalDofIdx])
+                    values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
+                break;
+            }
+            case 1:
+                if (freebcX_[globalDofIdx])
+                    values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
+
+                break;
+            case 2:
+                if (freebcYMinus_[globalDofIdx])
+                    values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
+
+                break;
+            case 3:
+                if (freebcY_[globalDofIdx])
+                    values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
+
+                break;
+            case 4:
+                if (freebcZMinus_[globalDofIdx])
+                    values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
+
+                break;
+            case 5:
+                if (freebcZ_[globalDofIdx])
+                    values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
+
+                break;
+            default:
+                throw std::logic_error("invalid face index for boundary condition");
+
+            }
+        }
     }
 
     /*!
@@ -1488,6 +1514,9 @@ public:
         return (oilVaporizationControl.getType() == Opm::OilVaporizationEnum::VAPPARS);
     }
 
+    bool hasFreeBoundaryConditions() const
+    { return hasFreeBoundaryConditions_; }
+
 private:
     bool drsdtActive_() const
     {
@@ -2264,6 +2293,35 @@ private:
         pffDofData_.update(distFn);
     }
 
+    void readBoundaryConditions_()
+    {
+        hasFreeBoundaryConditions_ = false;
+        readBoundaryConditionKeyword_("FREEBCX", freebcX_ );
+        readBoundaryConditionKeyword_("FREEBCX-", freebcXMinus_ );
+        readBoundaryConditionKeyword_("FREEBCY", freebcY_ );
+        readBoundaryConditionKeyword_("FREEBCY-", freebcYMinus_ );
+        readBoundaryConditionKeyword_("FREEBCZ", freebcZ_ );
+        readBoundaryConditionKeyword_("FREEBCZ-", freebcZMinus_ );
+    }
+
+    void readBoundaryConditionKeyword_(const std::string& name, std::vector<bool>& compressedData )
+    {
+        const auto& eclProps = this->simulator().vanguard().eclState().get3DProperties();
+        const auto& vanguard = this->simulator().vanguard();
+
+        unsigned numElems = vanguard.gridView().size(/*codim=*/0);
+        compressedData.resize(numElems, false);
+
+        if (eclProps.hasDeckDoubleGridProperty(name)) {
+            const std::vector<double>& data = eclProps.getDoubleGridProperty(name).getData();
+            for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
+                unsigned cartElemIdx = vanguard.cartesianIndex(elemIdx);
+                compressedData[elemIdx] = (data[cartElemIdx] > 0);
+            }
+            hasFreeBoundaryConditions_ = true;
+        }
+    }
+
     static std::string briefDescription_;
 
     std::vector<Scalar> porosity_;
@@ -2309,6 +2367,14 @@ private:
     PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
     TracerModel tracerModel_;
 
+    bool hasFreeBoundaryConditions_;
+    std::vector<bool> freebcX_;
+    std::vector<bool> freebcXMinus_;
+    std::vector<bool> freebcY_;
+    std::vector<bool> freebcYMinus_;
+    std::vector<bool> freebcZ_;
+    std::vector<bool> freebcZMinus_;
+
 
 };