added output writing for geomechanical module

ebos/eclgenericoutputblackoilmodule.cc

@@ -151,6 +151,7 @@ EclGenericOutputBlackoilModule(const EclipseState& eclState,
                            const SummaryState& summaryState,
                            bool enableEnergy,
                            bool enableTemperature,
+                           bool enableMech,
                            bool enableSolvent,
                            bool enablePolymer,
                            bool enableFoam,
@@ -166,6 +167,7 @@ EclGenericOutputBlackoilModule(const EclipseState& eclState,
     , logOutput_(eclState, schedule, summaryState)
     , enableEnergy_(enableEnergy)
     , enableTemperature_(enableTemperature)
+    , enableMech_(enableMech)
     , enableSolvent_(enableSolvent)
     , enablePolymer_(enablePolymer)
     , enableFoam_(enableFoam)
@@ -423,15 +425,26 @@ assignToSolution(data::Solution& sol)
     const auto extendedSolutionArrays = std::array {
         DataEntry{"BIOFILM",  UnitSystem::measure::identity,           cBiofilm_},
         DataEntry{"CALCITE",  UnitSystem::measure::identity,           cCalcite_},
+        DataEntry{"DELSTRXX", UnitSystem::measure::pressure,           delstressXX_},
+        DataEntry{"DELSTRYY", UnitSystem::measure::pressure,           delstressYY_},
+        DataEntry{"DELSTRZZ", UnitSystem::measure::pressure,           delstressZZ_},
+        DataEntry{"DELSTRXY", UnitSystem::measure::pressure,           delstressXY_},
+        DataEntry{"DELSTRXZ", UnitSystem::measure::pressure,           delstressXZ_},
+        DataEntry{"DELSTRYZ", UnitSystem::measure::pressure,           delstressYZ_},
+        DataEntry{"DISPX",    UnitSystem::measure::length,             dispX_},
+        DataEntry{"DISPY",    UnitSystem::measure::length,             dispY_},
+        DataEntry{"DISPZ",    UnitSystem::measure::length,             dispZ_},
         DataEntry{"DRSDTCON", UnitSystem::measure::gas_oil_ratio_rate, drsdtcon_},
         DataEntry{"KRNSW_GO", UnitSystem::measure::identity,           krnSwMdcGo_},
         DataEntry{"KRNSW_OW", UnitSystem::measure::identity,           krnSwMdcOw_},
+        DataEntry{"MECHPOTF", UnitSystem::measure::pressure,           mechPotentialForce_},
         DataEntry{"MICROBES", UnitSystem::measure::density,            cMicrobes_},
         DataEntry{"OXYGEN",   UnitSystem::measure::density,            cOxygen_},
         DataEntry{"PCSWM_GO", UnitSystem::measure::identity,           pcSwMdcGo_},
         DataEntry{"PCSWM_OW", UnitSystem::measure::identity,           pcSwMdcOw_},
         DataEntry{"PERMFACT", UnitSystem::measure::identity,           permFact_},
         DataEntry{"PORV_RC",  UnitSystem::measure::identity,           rockCompPorvMultiplier_},
+        DataEntry{"PRESPOTF", UnitSystem::measure::pressure,           mechPotentialPressForce_},
         DataEntry{"PRES_OVB", UnitSystem::measure::pressure,           overburdenPressure_},
         DataEntry{"RSW",      UnitSystem::measure::gas_oil_ratio,      rsw_},
         DataEntry{"RVW",      UnitSystem::measure::oil_gas_ratio,      rvw_},
@@ -442,6 +455,19 @@ assignToSolution(data::Solution& sol)
         DataEntry{"STD_CO2",  UnitSystem::measure::identity,           mFracCo2_},
         DataEntry{"STD_GAS",  UnitSystem::measure::identity,           mFracGas_},
         DataEntry{"STD_OIL",  UnitSystem::measure::identity,           mFracOil_},
+        DataEntry{"STRAINXX", UnitSystem::measure::identity,           strainXX_},
+        DataEntry{"STRAINYY", UnitSystem::measure::identity,           strainYY_},
+        DataEntry{"STRAINZZ", UnitSystem::measure::identity,           strainZZ_},
+        DataEntry{"STRAINXY", UnitSystem::measure::identity,           strainXY_},
+        DataEntry{"STRAINXZ", UnitSystem::measure::identity,           strainXZ_},
+        DataEntry{"STRAINYZ", UnitSystem::measure::identity,           strainYZ_},
+        DataEntry{"STRESSXX", UnitSystem::measure::length,             stressXX_},
+        DataEntry{"STRESSYY", UnitSystem::measure::length,             stressYY_},
+        DataEntry{"STRESSZZ", UnitSystem::measure::length,             stressZZ_},
+        DataEntry{"STRESSXY", UnitSystem::measure::length,             stressXY_},
+        DataEntry{"STRESSXZ", UnitSystem::measure::length,             stressXZ_},
+        DataEntry{"STRESSYZ", UnitSystem::measure::length,             stressYZ_},
+        DataEntry{"TEMPPOTF", UnitSystem::measure::pressure,           mechPotentialTempForce_},
         DataEntry{"TMULT_RC", UnitSystem::measure::identity,           rockCompTransMultiplier_},
         DataEntry{"UREA",     UnitSystem::measure::density,            cUrea_},
     };
@@ -760,6 +786,66 @@ doAllocBuffers(unsigned bufferSize,
     rstKeywords["PRES"] = 0;
     rstKeywords["PRESSURE"] = 0;
 
+    if (enableMech_ && eclState_.runspec().mech()) {
+        this->mechPotentialForce_.resize(bufferSize,0.0);
+        rstKeywords["MECHPOTF"] = 0;
+        this->mechPotentialTempForce_.resize(bufferSize,0.0);
+        rstKeywords["TEMPPOTF"] = 0;
+        this->mechPotentialPressForce_.resize(bufferSize,0.0);
+        rstKeywords["PRESPOTF"] = 0;
+
+        this->dispX_.resize(bufferSize,0.0);
+        rstKeywords["DISPX"] = 0;
+        this->dispY_.resize(bufferSize,0.0);
+        rstKeywords["DISPY"] = 0;
+        this->dispZ_.resize(bufferSize,0.0);
+        rstKeywords["DISPZ"] = 0;
+        this->stressXX_.resize(bufferSize,0.0);
+        rstKeywords["STRESSXX"] = 0;
+        this->stressYY_.resize(bufferSize,0.0);
+        rstKeywords["STRESSYY"] = 0;
+        this->stressZZ_.resize(bufferSize,0.0);
+        rstKeywords["STRESSZZ"] = 0;
+        this->stressXY_.resize(bufferSize,0.0);
+        rstKeywords["STRESSXY"] = 0;
+        this->stressXZ_.resize(bufferSize,0.0);
+        rstKeywords["STRESSXZ"] = 0;
+        this->stressXY_.resize(bufferSize,0.0);
+        rstKeywords["STRESSXY"] = 0;
+        this->stressYZ_.resize(bufferSize,0.0);
+        rstKeywords["STRESSYZ"] = 0;
+
+        this->strainXX_.resize(bufferSize,0.0);
+        rstKeywords["STRAINXX"] = 0;
+        this->strainYY_.resize(bufferSize,0.0);
+        rstKeywords["STRAINYY"] = 0;
+        this->strainZZ_.resize(bufferSize,0.0);
+        rstKeywords["STRAINZZ"] = 0;
+        this->strainXY_.resize(bufferSize,0.0);
+        rstKeywords["STRAINXY"] = 0;
+        this->strainXZ_.resize(bufferSize,0.0);
+        rstKeywords["STRAINXZ"] = 0;
+        this->strainXY_.resize(bufferSize,0.0);
+        rstKeywords["STRAINXY"] = 0;
+        this->strainYZ_.resize(bufferSize,0.0);
+        rstKeywords["STRAINYZ"] = 0;
+
+        this->delstressXX_.resize(bufferSize,0.0);
+        rstKeywords["DELSTRXX"] = 0;
+        this->delstressYY_.resize(bufferSize,0.0);
+        rstKeywords["DELSTRYY"] = 0;
+        this->delstressZZ_.resize(bufferSize,0.0);
+        rstKeywords["DELSTRZZ"] = 0;
+        this->delstressXY_.resize(bufferSize,0.0);
+        rstKeywords["DELSTRXY"] = 0;
+        this->delstressXZ_.resize(bufferSize,0.0);
+        rstKeywords["DELSTRXZ"] = 0;
+        this->delstressXY_.resize(bufferSize,0.0);
+        rstKeywords["DELSTRXY"] = 0;
+        this->delstressYZ_.resize(bufferSize,0.0);
+        rstKeywords["DELSTRYZ"] = 0;
+    }
+
     // If TEMP is set in RPTRST we output temperature even if THERMAL
     // is not activated
     if (enableEnergy_ || rstKeywords["TEMP"] > 0) {

ebos/eclgenericoutputblackoilmodule.hh

@@ -56,11 +56,11 @@ public:
      Scalar* getPRESSURE_ptr(void) {
         return (this->fluidPressure_.data()) ;
     };
-    
+
     int  getPRESSURE_size( void ) {
         return (this->fluidPressure_.size()) ;
     };
-    
+
     // write cumulative production and injection reports to output
     void outputCumLog(std::size_t reportStepNum,
                       const bool substep,
@@ -286,6 +286,7 @@ protected:
                                    const SummaryState& summaryState,
                                    bool enableEnergy,
                                    bool enableTemperature,
+                                   bool enableMech,
                                    bool enableSolvent,
                                    bool enablePolymer,
                                    bool enableFoam,
@@ -347,6 +348,7 @@ protected:
 
     bool enableEnergy_;
     bool enableTemperature_;
+    bool enableMech_;
 
     bool enableSolvent_;
     bool enablePolymer_;
@@ -427,6 +429,33 @@ protected:
     ScalarBuffer pcow_;
     ScalarBuffer pcog_;
 
+    // buffers for mechanical output
+    ScalarBuffer mechPotentialForce_;
+    ScalarBuffer mechPotentialPressForce_;
+    ScalarBuffer mechPotentialTempForce_;
+
+    ScalarBuffer dispX_;
+    ScalarBuffer dispY_;
+    ScalarBuffer dispZ_;
+    ScalarBuffer stressXX_;
+    ScalarBuffer stressYY_;
+    ScalarBuffer stressZZ_;
+    ScalarBuffer stressXY_;
+    ScalarBuffer stressXZ_;
+    ScalarBuffer stressYZ_;
+    ScalarBuffer delstressXX_;
+    ScalarBuffer delstressYY_;
+    ScalarBuffer delstressZZ_;
+    ScalarBuffer delstressXY_;
+    ScalarBuffer delstressXZ_;
+    ScalarBuffer delstressYZ_;
+    ScalarBuffer strainXX_;
+    ScalarBuffer strainYY_;
+    ScalarBuffer strainZZ_;
+    ScalarBuffer strainXY_;
+    ScalarBuffer strainXZ_;
+    ScalarBuffer strainYZ_;
+
     std::array<ScalarBuffer, numPhases> saturation_;
     std::array<ScalarBuffer, numPhases> invB_;
     std::array<ScalarBuffer, numPhases> density_;

ebos/ecloutputblackoilmodule.hh

@@ -145,6 +145,7 @@ public:
                    simulator.vanguard().summaryState(),
                    getPropValue<TypeTag, Properties::EnableEnergy>(),
                    getPropValue<TypeTag, Properties::EnableTemperature>(),
+                   getPropValue<TypeTag, Properties::EnableMech>(),
                    getPropValue<TypeTag, Properties::EnableSolvent>(),
                    getPropValue<TypeTag, Properties::EnablePolymer>(),
                    getPropValue<TypeTag, Properties::EnableFoam>(),
@@ -216,6 +217,53 @@ public:
                              problem.eclWriter()->getOutputNnc().size());
     }
 
+    void processElementMech(const ElementContext& elemCtx)
+    {
+        if constexpr (getPropValue<TypeTag, Properties::EnableMech>()) {
+            const auto& problem = elemCtx.simulator().problem();
+            const auto& model = problem.geoMechModel();
+            for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
+                const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
+                const auto& fs = intQuants.fluidState();
+                unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
+                if (!this->mechPotentialForce_.empty()) {
+                    //assume all mechancal tings should be written
+                    this->mechPotentialForce_[globalDofIdx] = model.mechPotentialForce(globalDofIdx);
+                    this->mechPotentialPressForce_[globalDofIdx] = model.mechPotentialPressForce(globalDofIdx);
+                    this->mechPotentialTempForce_[globalDofIdx] = model.mechPotentialTempForce(globalDofIdx);
+
+                    this->dispX_[globalDofIdx] = model.disp(globalDofIdx, 0);
+                    this->dispY_[globalDofIdx] = model.disp(globalDofIdx, 1);
+                    this->dispZ_[globalDofIdx] = model.disp(globalDofIdx, 2);
+                    this->stressXX_[globalDofIdx] = model.stress(globalDofIdx, 0);
+                    this->stressYY_[globalDofIdx] = model.stress(globalDofIdx, 1);
+                    this->stressZZ_[globalDofIdx] = model.stress(globalDofIdx, 2);
+                    //voight notation
+                    this->stressXY_[globalDofIdx] = model.stress(globalDofIdx, 5);
+                    this->stressXZ_[globalDofIdx] = model.stress(globalDofIdx, 4);
+                    this->stressYZ_[globalDofIdx] = model.stress(globalDofIdx, 3);
+
+                    this->strainXX_[globalDofIdx] = model.strain(globalDofIdx, 0);
+                    this->strainYY_[globalDofIdx] = model.strain(globalDofIdx, 1);
+                    this->strainZZ_[globalDofIdx] = model.strain(globalDofIdx, 2);
+                    //voight notation
+                    this->strainXY_[globalDofIdx] = model.strain(globalDofIdx, 5);
+                    this->strainXZ_[globalDofIdx] = model.strain(globalDofIdx, 4);
+                    this->strainYZ_[globalDofIdx] = model.strain(globalDofIdx, 3);
+
+
+                    this->delstressXX_[globalDofIdx] = model.delstress(globalDofIdx, 0);
+                    this->delstressYY_[globalDofIdx] = model.delstress(globalDofIdx, 1);
+                    this->delstressZZ_[globalDofIdx] = model.delstress(globalDofIdx, 2);
+                    //voight notation
+                    this->delstressXY_[globalDofIdx] = model.delstress(globalDofIdx, 5);
+                    this->delstressXZ_[globalDofIdx] = model.delstress(globalDofIdx, 4);
+                    this->delstressYZ_[globalDofIdx] = model.delstress(globalDofIdx, 3);
+                }
+            }
+        }
+    }
+
     /*!
      * \brief Modify the internal buffers according to the intensive
      *        quanties relevant for an element

ebos/eclwriter.hh

@@ -123,6 +123,7 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
     using BaseType = EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>;
 
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
+    enum { enableMech = getPropValue<TypeTag, Properties::EnableMech>() };
     enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
     enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
 
@@ -574,6 +575,17 @@ private:
             }
         }
 
+        if constexpr (enableMech) {
+            if (simulator_.vanguard().eclState().runspec().mech()) {
+                OPM_TIMEBLOCK(prepareMechData);
+                for (const auto& elem : elements(gridView)) {
+                    elemCtx.updatePrimaryStencil(elem);
+                    elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+                    eclOutputModule_->processElementMech(elemCtx);
+                }
+            }
+        }
+
         if (! this->simulator_.model().linearizer().getFlowsInfo().empty()) {
             OPM_TIMEBLOCK(prepareFlowsData);
             for (const auto& elem : elements(gridView)) {
@@ -660,6 +672,7 @@ private:
     bool damarisUpdate_ = false;  ///< Whenever this is true writeOutput() will set up Damaris offsets of model fields
 #endif
 };
+
 } // namespace Opm
 
 #endif