From 646d7c0a041fd91015f17f07c207ff5a4fed9e11 Mon Sep 17 00:00:00 2001 From: Halvor M Nilsen Date: Fri, 31 Jan 2025 18:24:09 +0100 Subject: [PATCH 1/2] Output More Geo-Mechanical Quantities to Summary/Restart In particular, this commit defines multiple additional 3D stress quantities in the restart file as well as some cell level stress vectors in the summary file. --- .../flow/GenericOutputBlackoilModule.cpp | 54 +++++++- .../flow/GenericOutputBlackoilModule.hpp | 12 ++ opm/simulators/flow/OutputBlackoilModule.hpp | 124 +++++++++++++----- 3 files changed, 153 insertions(+), 37 deletions(-) diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index abd8f3cc3..4155d0b06 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -548,12 +548,24 @@ assignToSolution(data::Solution& sol) 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{"STRESSXX", UnitSystem::measure::pressure, stressXX_}, + DataEntry{"STRESSYY", UnitSystem::measure::pressure, stressYY_}, + DataEntry{"STRESSZZ", UnitSystem::measure::pressure, stressZZ_}, + DataEntry{"STRESSXY", UnitSystem::measure::pressure, stressXY_}, + DataEntry{"STRESSXZ", UnitSystem::measure::pressure, stressXZ_}, + DataEntry{"STRESSYZ", UnitSystem::measure::pressure, stressYZ_}, + DataEntry{"LINSTRXX", UnitSystem::measure::pressure, linstressXX_}, + DataEntry{"LINSTRYY", UnitSystem::measure::pressure, linstressYY_}, + DataEntry{"LINSTRZZ", UnitSystem::measure::pressure, linstressZZ_}, + DataEntry{"LINSTRXY", UnitSystem::measure::pressure, linstressXY_}, + DataEntry{"LINSTRXZ", UnitSystem::measure::pressure, linstressXZ_}, + DataEntry{"LINSTRYZ", UnitSystem::measure::pressure, linstressYZ_}, + DataEntry{"FRCSTRXX", UnitSystem::measure::pressure, fracstressXX_}, + DataEntry{"FRCSTRYY", UnitSystem::measure::pressure, fracstressYY_}, + DataEntry{"FRCSTRZZ", UnitSystem::measure::pressure, fracstressZZ_}, + DataEntry{"FRCSTRXY", UnitSystem::measure::pressure, fracstressXY_}, + DataEntry{"FRCSTRXZ", UnitSystem::measure::pressure, fracstressXZ_}, + DataEntry{"FRCSTRYZ", UnitSystem::measure::pressure, fracstressYZ_}, DataEntry{"TEMPPOTF", UnitSystem::measure::pressure, mechPotentialTempForce_}, DataEntry{"TMULT_RC", UnitSystem::measure::identity, rockCompTransMultiplier_}, DataEntry{"UREA", UnitSystem::measure::density, cUrea_}, @@ -1037,6 +1049,36 @@ doAllocBuffers(const unsigned bufferSize, rstKeywords["DELSTRXY"] = 0; this->delstressYZ_.resize(bufferSize,0.0); rstKeywords["DELSTRYZ"] = 0; + + this->fracstressXX_.resize(bufferSize,0.0); + rstKeywords["FRCSTRXX"] = 0; + this->fracstressYY_.resize(bufferSize,0.0); + rstKeywords["FRCSTRYY"] = 0; + this->fracstressZZ_.resize(bufferSize,0.0); + rstKeywords["FRCSTRZZ"] = 0; + this->fracstressXY_.resize(bufferSize,0.0); + rstKeywords["FRCSTRXY"] = 0; + this->fracstressXZ_.resize(bufferSize,0.0); + rstKeywords["FRCSTRXZ"] = 0; + this->fracstressXY_.resize(bufferSize,0.0); + rstKeywords["FRCSTRXY"] = 0; + this->fracstressYZ_.resize(bufferSize,0.0); + rstKeywords["FRCSTRYZ"] = 0; + + this->linstressXX_.resize(bufferSize,0.0); + rstKeywords["LINSTRXX"] = 0; + this->linstressYY_.resize(bufferSize,0.0); + rstKeywords["LINSTRYY"] = 0; + this->linstressZZ_.resize(bufferSize,0.0); + rstKeywords["LINSTRZZ"] = 0; + this->linstressXY_.resize(bufferSize,0.0); + rstKeywords["LINSTRXY"] = 0; + this->linstressXZ_.resize(bufferSize,0.0); + rstKeywords["LINSTRXZ"] = 0; + this->linstressXY_.resize(bufferSize,0.0); + rstKeywords["LINSTRXY"] = 0; + this->linstressYZ_.resize(bufferSize,0.0); + rstKeywords["LINSTRYZ"] = 0; } // If TEMP is set in RPTRST we output temperature even if THERMAL diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index f148b4080..d88bd3441 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -522,6 +522,18 @@ protected: ScalarBuffer delstressXY_; ScalarBuffer delstressXZ_; ScalarBuffer delstressYZ_; + ScalarBuffer linstressXX_; + ScalarBuffer linstressYY_; + ScalarBuffer linstressZZ_; + ScalarBuffer linstressXY_; + ScalarBuffer linstressXZ_; + ScalarBuffer linstressYZ_; + ScalarBuffer fracstressXX_; + ScalarBuffer fracstressYY_; + ScalarBuffer fracstressZZ_; + ScalarBuffer fracstressXY_; + ScalarBuffer fracstressXZ_; + ScalarBuffer fracstressYZ_; ScalarBuffer strainXX_; ScalarBuffer strainYY_; ScalarBuffer strainZZ_; diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 44134cd5e..875cb1f57 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -196,44 +196,70 @@ public: void processElementMech(const ElementContext& elemCtx) { if constexpr (getPropValue()) { + if (this->mechPotentialForce_.empty()) { + return; + } + + enum Voigt { XX = 0, YY = 1, ZZ = 2, YZ = 3, XZ = 4, XY = 5, }; + const auto& problem = elemCtx.simulator().problem(); const auto& model = problem.geoMechModel(); + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { - unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - if (!this->mechPotentialForce_.empty()) { - // assume all mechanical things should be written - this->mechPotentialForce_[globalDofIdx] = model.mechPotentialForce(globalDofIdx); - this->mechPotentialPressForce_[globalDofIdx] = model.mechPotentialPressForce(globalDofIdx); - this->mechPotentialTempForce_[globalDofIdx] = model.mechPotentialTempForce(globalDofIdx); + const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - 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); + // Assume all mechanical things should be written + this->mechPotentialForce_[globalDofIdx] = model.mechPotentialForce(globalDofIdx); + this->mechPotentialPressForce_[globalDofIdx] = model.mechPotentialPressForce(globalDofIdx); + this->mechPotentialTempForce_[globalDofIdx] = model.mechPotentialTempForce(globalDofIdx); - 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); + const auto disp = model.disp(globalDofIdx, /*include_fracture*/true); + this->dispX_[globalDofIdx] = disp[Voigt::XX]; + this->dispY_[globalDofIdx] = disp[Voigt::YY]; + this->dispZ_[globalDofIdx] = disp[Voigt::ZZ]; + // Total stress is not stored but calculated result is Voigt notation + const auto stress = model.stress(globalDofIdx, /*include_fracture*/true); + this->stressXX_[globalDofIdx] = stress[Voigt::XX]; + this->stressYY_[globalDofIdx] = stress[Voigt::YY]; + this->stressZZ_[globalDofIdx] = stress[Voigt::ZZ]; + this->stressXY_[globalDofIdx] = stress[Voigt::XY]; + this->stressXZ_[globalDofIdx] = stress[Voigt::XZ]; + this->stressYZ_[globalDofIdx] = stress[Voigt::YZ]; - 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); - } + const auto strain = model.strain(globalDofIdx, /*include_fracture*/true); + this->strainXX_[globalDofIdx] = strain[Voigt::XX]; + this->strainYY_[globalDofIdx] = strain[Voigt::YY]; + this->strainZZ_[globalDofIdx] = strain[Voigt::ZZ]; + this->strainXY_[globalDofIdx] = strain[Voigt::XY]; + this->strainXZ_[globalDofIdx] = strain[Voigt::XZ]; + this->strainYZ_[globalDofIdx] = strain[Voigt::YZ]; + + // Not including fracture + const auto delstress = model.delstress(globalDofIdx); + this->delstressXX_[globalDofIdx] = delstress[Voigt::XX]; + this->delstressYY_[globalDofIdx] = delstress[Voigt::YY]; + this->delstressZZ_[globalDofIdx] = delstress[Voigt::ZZ]; + this->delstressXY_[globalDofIdx] = delstress[Voigt::XY]; + this->delstressXZ_[globalDofIdx] = delstress[Voigt::XZ]; + this->delstressYZ_[globalDofIdx] = delstress[Voigt::YZ]; + + const auto linstress = model.linstress(globalDofIdx); + this->linstressXX_[globalDofIdx] = linstress[Voigt::XX]; + this->linstressYY_[globalDofIdx] = linstress[Voigt::YY]; + this->linstressZZ_[globalDofIdx] = linstress[Voigt::ZZ]; + this->linstressXY_[globalDofIdx] = linstress[Voigt::XY]; + this->linstressXZ_[globalDofIdx] = linstress[Voigt::XZ]; + this->linstressYZ_[globalDofIdx] = linstress[Voigt::YZ]; + + // is the tresagii stress which make rock fracture + const auto fracstress = model.fractureStress(globalDofIdx); + this->fracstressXX_[globalDofIdx] = fracstress[Voigt::XX]; + this->fracstressYY_[globalDofIdx] = fracstress[Voigt::YY]; + this->fracstressZZ_[globalDofIdx] = fracstress[Voigt::ZZ]; + this->fracstressXY_[globalDofIdx] = fracstress[Voigt::XY]; + this->fracstressXZ_[globalDofIdx] = fracstress[Voigt::XZ]; + this->fracstressYZ_[globalDofIdx] = fracstress[Voigt::YZ]; } } } @@ -812,6 +838,7 @@ public: return; const auto& problem = elemCtx.simulator().problem(); + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { // Adding block data const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); @@ -850,6 +877,30 @@ public: else if (FluidSystem::phaseIsActive(waterPhaseIdx)) val.second = getValue(fs.temperature(waterPhaseIdx)); } + else if ((key.first == "BSTRSSXX") || + (key.first == "BSTRSSYY") || + (key.first == "BSTRSSZZ") || + (key.first == "BSTRSSXY") || + (key.first == "BSTRSSXZ") || + (key.first == "BSTRSSYZ")) + { + if constexpr (HasGeoMech>::value) { + enum Voigt { XX = 0, YY = 1, ZZ = 2, YZ = 3, XZ = 4, XY = 5, }; + + const auto stress = problem.geoMechModel() + .stress(globalDofIdx, /*include_fracture*/ true); + + if (key.first == "BSTRSSXX") { val.second = stress[Voigt::XX]; } + else if (key.first == "BSTRSSYY") { val.second = stress[Voigt::YY]; } + else if (key.first == "BSTRSSZZ") { val.second = stress[Voigt::ZZ]; } + else if (key.first == "BSTRSSXY") { val.second = stress[Voigt::XY]; } + else if (key.first == "BSTRSSXZ") { val.second = stress[Voigt::XZ]; } + else /* BSTRSSYZ */{ val.second = stress[Voigt::YZ]; } + } + else { + val.second = 0.0; + } + } else if (key.first == "BWKR" || key.first == "BKRW") val.second = getValue(intQuants.relativePermeability(waterPhaseIdx)); else if (key.first == "BGKR" || key.first == "BKRG") @@ -1243,6 +1294,17 @@ public: } private: + template + using RemoveCVR = std::remove_cv_t>; + + template + struct HasGeoMech : public std::false_type {}; + + template + struct HasGeoMech< + Problem, std::void_t().geoMechModel())> + > : public std::true_type {}; + bool isDefunctParallelWell(std::string wname) const override { if (simulator_.gridView().comm().size() == 1) From 9e235def1c60aeb23ec6e950b5c4609c60c0c416 Mon Sep 17 00:00:00 2001 From: Halvor M Nilsen Date: Fri, 31 Jan 2025 18:31:56 +0100 Subject: [PATCH 2/2] Add Way of Introducing New Well Connections at Runtime This commit adds a new member function void WellInterfaceGeneric<>::addPerforation() which is, initially, intended as a way of introducing new well connections arising from hydraulic fracturing. We expect more work in this area. --- CMakeLists_files.cmake | 1 + opm/simulators/wells/RuntimePerforation.hpp | 42 +++++++++++++++++++ opm/simulators/wells/WellInterfaceGeneric.cpp | 36 ++++++++++++++++ opm/simulators/wells/WellInterfaceGeneric.hpp | 4 +- 4 files changed, 82 insertions(+), 1 deletion(-) create mode 100644 opm/simulators/wells/RuntimePerforation.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f366b92b1..ff1e40129 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -1029,6 +1029,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/RatioCalculator.hpp opm/simulators/wells/RegionAttributeHelpers.hpp opm/simulators/wells/RegionAverageCalculator.hpp + opm/simulators/wells/RuntimePerforation.hpp opm/simulators/wells/SingleWellState.hpp opm/simulators/wells/StandardWell.hpp opm/simulators/wells/StandardWell_impl.hpp diff --git a/opm/simulators/wells/RuntimePerforation.hpp b/opm/simulators/wells/RuntimePerforation.hpp new file mode 100644 index 000000000..18038b40e --- /dev/null +++ b/opm/simulators/wells/RuntimePerforation.hpp @@ -0,0 +1,42 @@ +/* + Copyright 2025 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_RUNTIME_PERFORATION_HPP_INCLUDED +#define OPM_RUNTIME_PERFORATION_HPP_INCLUDED + +namespace Opm { + +/// Simple model of a well connection created at runtime, possibly as a +/// result of a geo-mechanical fracturing process. +struct RuntimePerforation +{ + /// Active cell index, on current rank, that is dynamically perforated + /// by a well. + int cell{}; + + /// Connection's transmissibility factor. + double ctf{}; + + /// Depth at which the new connection is created. + double depth{}; +}; + +} // namespace Opm + +#endif // OPM_RUNTIME_PERFORATION_HPP_INCLUDED diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 26ffb2f36..b242c7ddb 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -20,6 +20,7 @@ */ #include + #include #include @@ -33,7 +34,9 @@ #include #include #include + #include + #include #include #include @@ -45,10 +48,13 @@ #include +#include #include #include #include #include +#include +#include namespace Opm { @@ -676,6 +682,36 @@ void WellInterfaceGeneric::resetWellOperability() this->operability_status_.resetOperability(); } +template +void WellInterfaceGeneric::addPerforations(const std::vector& perfs) +{ + for (const auto& perf : perfs) { + auto it = std::find(well_cells_.begin(), well_cells_.end(), perf.cell); + if (it != this->well_cells_.end()) { + // If perforation to cell already exists, just add contribution. + const auto ind = std::distance(this->well_cells_.begin(), it); + this->well_index_[ind] += static_cast(perf.ctf); + } + else { + this->well_cells_.push_back(perf.cell); + this->well_index_.push_back(static_cast(perf.ctf)); + this->perf_depth_.push_back(static_cast(perf.depth)); + + // Not strictly needed. + const double nan = std::nan("1"); + this->perf_rep_radius_.push_back(nan); + this->perf_length_.push_back(nan); + this->bore_diameters_.push_back(nan); + + // For now use the saturation table for the first cell. + this->saturation_table_number_ + .push_back(this->saturation_table_number_.front()); + + ++this->number_of_perforations_; + } + } +} + template Scalar WellInterfaceGeneric::wmicrobes_() const { diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index f511eb097..f583f30b2 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -26,6 +26,7 @@ #include #include +#include #include #include @@ -189,7 +190,6 @@ public: void resetWellOperability(); - virtual std::vector getPrimaryVars() const { return {}; @@ -203,6 +203,8 @@ public: virtual Scalar connectionDensity(const int globalConnIdx, const int openConnIdx) const = 0; + void addPerforations(const std::vector& perfs); + protected: bool getAllowCrossFlow() const;