From 3457d46a00b48066d2183d3f4abcd1e8c2da0957 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 8 Jan 2021 10:44:40 +0100 Subject: [PATCH 1/5] add flow_co2 gas-oil simulator with diffusion for testing --- CMakeLists.txt | 10 ++++++ flow/flow_co2_diffusion.cpp | 64 +++++++++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+) create mode 100644 flow/flow_co2_diffusion.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 1e7f592a9..a27fded43 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -339,6 +339,16 @@ opm_add_test(flow_onephase DEPENDS opmsimulators LIBRARIES opmsimulators) +opm_add_test(flow_co2 + ONLY_COMPILE + DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} + SOURCES + flow/flow_co2_diffusion.cpp + $ + EXE_NAME flow_co2 + DEPENDS opmsimulators + LIBRARIES opmsimulators) + opm_add_test(flow_onephase_energy ONLY_COMPILE DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} diff --git a/flow/flow_co2_diffusion.cpp b/flow/flow_co2_diffusion.cpp new file mode 100644 index 000000000..548108b63 --- /dev/null +++ b/flow/flow_co2_diffusion.cpp @@ -0,0 +1,64 @@ +/* + 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 . +*/ +#include "config.h" +#include +#include + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlowCo2Problem { + using InheritsFrom = std::tuple; +}; +} + +template +struct EnableDiffusion { + static constexpr bool value = true; +}; + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::EclFlowProblem; + using FluidSystem = GetPropType; + +public: + typedef Opm::BlackOilTwoPhaseIndices(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*disabledCompIdx=*/FluidSystem::waterCompIdx> type; +}; +}} + + +// ----------------- Main program ----------------- +int main(int argc, char** argv) +{ + using TypeTag = Opm::Properties::TTag::EclFlowCo2Problem; + auto mainObject = Opm::Main(argc, argv); + return mainObject.runStatic(); +} From 1abbd9c77642e70e1cc31518289cb7f8952a3b5d Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 14 Jan 2021 15:33:28 +0100 Subject: [PATCH 2/5] add diffusivity to eclTranmissibility --- ebos/eclproblem.hh | 17 ++++++++ ebos/ecltransmissibility.hh | 82 +++++++++++++++++++++++++++++++++++++ 2 files changed, 99 insertions(+) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index cfc4abfff..a5258d301 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -622,6 +622,7 @@ class EclProblem : public GetPropType enum { enableExtbo = getPropValue() }; enum { enableTemperature = getPropValue() }; enum { enableEnergy = getPropValue() }; + enum { enableDiffusion = getPropValue() }; enum { enableThermalFluxBoundaries = getPropValue() }; enum { enableApiTracking = getPropValue() }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; @@ -1372,6 +1373,18 @@ public: return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility; } + /*! + * \copydoc EclTransmissiblity::diffusivity + */ + template + Scalar diffusivity(const Context& context, + unsigned OPM_OPTIM_UNUSED fromDofLocalIdx, + unsigned toDofLocalIdx) const + { + assert(fromDofLocalIdx == 0); + return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity; + } + /*! * \copydoc EclTransmissiblity::transmissibilityBoundary */ @@ -3072,6 +3085,8 @@ private: struct PffDofData_ { Opm::ConditionalStorage thermalHalfTrans; + Opm::ConditionalStorage diffusivity; + //Scalar diffusivity; Scalar transmissibility; }; @@ -3093,6 +3108,8 @@ private: if (enableEnergy) *dofData.thermalHalfTrans = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx); + if (enableDiffusion) + *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx); } }; diff --git a/ebos/ecltransmissibility.hh b/ebos/ecltransmissibility.hh index 3218da5c0..9d2bcaf3b 100644 --- a/ebos/ecltransmissibility.hh +++ b/ebos/ecltransmissibility.hh @@ -72,6 +72,7 @@ class EclTransmissibility using Intersection = typename GridView::Intersection; static const bool enableEnergy = getPropValue(); + static const bool enableDiffusion = getPropValue(); // Grid and world dimension enum { dimWorld = GridView::dimensionworld }; @@ -181,6 +182,13 @@ public: thermalHalfTransBoundary_.clear(); } + // if diffusion is enabled, let's do the same for the "diffusivity" + if (enableDiffusion) { + diffusivity_->clear(); + diffusivity_->reserve(numElements*3*1.05); + extractPorosity_(); + } + // The MULTZ needs special case if the option is ALL // Then the smallest multiplier is applied. // Default is to apply the top and bottom multiplier @@ -386,6 +394,42 @@ public: faceDir); trans_[isId_(elemIdx, outsideElemIdx)] = trans; + + // update the "thermal half transmissibility" for the intersection + if (enableDiffusion) { + + Scalar halfDiffusivity1; + Scalar halfDiffusivity2; + + computeHalfDiffusivity_(halfDiffusivity1, + faceAreaNormal, + distanceVector_(faceCenterInside, + intersection.indexInInside(), + elemIdx, + axisCentroids), + porosity_[elemIdx]); + computeHalfDiffusivity_(halfDiffusivity2, + faceAreaNormal, + distanceVector_(faceCenterOutside, + intersection.indexInOutside(), + outsideElemIdx, + axisCentroids), + porosity_[outsideElemIdx]); + + applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg); + applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg); + + //TODO Add support for multipliers + Scalar diffusivity; + if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30) + // avoid division by zero + diffusivity = 0.0; + else + diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2); + + + (*diffusivity_)[isId_(elemIdx, outsideElemIdx)] = diffusivity; + } } } @@ -450,6 +494,12 @@ public: Scalar thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const { return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx)); } + /*! + * \brief Return the diffusivity for the intersection between two elements. + */ + Scalar diffusivity(unsigned elemIdx1, unsigned elemIdx2) const + { return diffusivity_->at(isId_(elemIdx1, elemIdx2)); } + private: void removeSmallNonCartesianTransmissibilities_() @@ -881,6 +931,21 @@ private: "(The PERM{X,Y,Z} keywords are missing)"); } + void extractPorosity_() + { + // read the intrinsic porosity from the eclState. Note that all arrays + // provided by eclState are one-per-cell of "uncompressed" grid, whereas the + // simulation grid might remove a few elements. (e.g. because it is distributed + // over several processes.) + const auto& fp = vanguard_.eclState().fieldProps(); + if (fp.has_double("PORO")) { + porosity_ = fp.get_double("PORO"); + } + else + throw std::logic_error("Can't read the porosityfrom the ecl state. " + "(The PORO keywords are missing)"); + } + std::uint64_t isId_(std::uint32_t elemIdx1, std::uint32_t elemIdx2) const { std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2); @@ -924,6 +989,20 @@ private: halfTrans /= distance.two_norm2(); } + void computeHalfDiffusivity_(Scalar& halfDiff, + const DimVector& areaNormal, + const DimVector& distance, + const Scalar& poro) const + { + halfDiff = poro; + Scalar val = 0; + for (unsigned i = 0; i < areaNormal.size(); ++i) + val += areaNormal[i]*distance[i]; + + halfDiff *= std::abs(val); + halfDiff /= distance.two_norm2(); + } + DimVector distanceVector_(const DimVector& center, int faceIdx, // in the reference element that contains the intersection unsigned elemIdx, @@ -1000,11 +1079,14 @@ private: const Vanguard& vanguard_; Scalar transmissibilityThreshold_; std::vector permeability_; + std::vector porosity_; std::unordered_map trans_; std::map, Scalar> transBoundary_; std::map, Scalar> thermalHalfTransBoundary_; Opm::ConditionalStorage > thermalHalfTrans_; + Opm::ConditionalStorage > diffusivity_; }; } // namespace Opm From 361de71c960625b60225d66a22cd528b79b1fa5b Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 20 Jan 2021 09:02:44 +0100 Subject: [PATCH 3/5] add optional flow_bo_diffusion simulator for testing --- CMakeLists.txt | 10 +++++++++ flow/flow_bo_diffusion.cpp | 44 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+) create mode 100644 flow/flow_bo_diffusion.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index a27fded43..f802662a3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -349,6 +349,16 @@ opm_add_test(flow_co2 DEPENDS opmsimulators LIBRARIES opmsimulators) +opm_add_test(flow_bo_diffusion + ONLY_COMPILE + DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} + SOURCES + flow/flow_bo_diffusion.cpp + $ + EXE_NAME flow_bo_diffusion + DEPENDS opmsimulators + LIBRARIES opmsimulators) + opm_add_test(flow_onephase_energy ONLY_COMPILE DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} diff --git a/flow/flow_bo_diffusion.cpp b/flow/flow_bo_diffusion.cpp new file mode 100644 index 000000000..cdeb74e55 --- /dev/null +++ b/flow/flow_bo_diffusion.cpp @@ -0,0 +1,44 @@ +/* + 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 . +*/ +#include "config.h" +#include +#include + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlowCo2Problem { + using InheritsFrom = std::tuple; +}; +} + +template +struct EnableDiffusion { + static constexpr bool value = true; +}; + + +}} + + +// ----------------- Main program ----------------- +int main(int argc, char** argv) +{ + using TypeTag = Opm::Properties::TTag::EclFlowCo2Problem; + auto mainObject = Opm::Main(argc, argv); + return mainObject.runStatic(); +} From 5deb410bd5ef943d2b97b02647063324fa347402 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 21 Jan 2021 12:06:26 +0100 Subject: [PATCH 4/5] clean-up in flow_bo_diffusion.cpp and flow_co2.cpp --- CMakeLists.txt | 4 ++-- flow/flow_bo_diffusion.cpp | 7 +++---- flow/flow_co2_diffusion.cpp | 2 +- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f802662a3..678fb41ca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -339,13 +339,13 @@ opm_add_test(flow_onephase DEPENDS opmsimulators LIBRARIES opmsimulators) -opm_add_test(flow_co2 +opm_add_test(flow_co2_diffusion ONLY_COMPILE DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} SOURCES flow/flow_co2_diffusion.cpp $ - EXE_NAME flow_co2 + EXE_NAME flow_co2_diffusion DEPENDS opmsimulators LIBRARIES opmsimulators) diff --git a/flow/flow_bo_diffusion.cpp b/flow/flow_bo_diffusion.cpp index cdeb74e55..f93e638de 100644 --- a/flow/flow_bo_diffusion.cpp +++ b/flow/flow_bo_diffusion.cpp @@ -16,18 +16,17 @@ */ #include "config.h" #include -#include namespace Opm { namespace Properties { namespace TTag { -struct EclFlowCo2Problem { +struct EclFlowDiffusionProblem { using InheritsFrom = std::tuple; }; } template -struct EnableDiffusion { +struct EnableDiffusion { static constexpr bool value = true; }; @@ -38,7 +37,7 @@ struct EnableDiffusion { // ----------------- Main program ----------------- int main(int argc, char** argv) { - using TypeTag = Opm::Properties::TTag::EclFlowCo2Problem; + using TypeTag = Opm::Properties::TTag::EclFlowDiffusionProblem; auto mainObject = Opm::Main(argc, argv); return mainObject.runStatic(); } diff --git a/flow/flow_co2_diffusion.cpp b/flow/flow_co2_diffusion.cpp index 548108b63..8b4fe9a95 100644 --- a/flow/flow_co2_diffusion.cpp +++ b/flow/flow_co2_diffusion.cpp @@ -27,7 +27,7 @@ struct EclFlowCo2Problem { } template -struct EnableDiffusion { +struct EnableDiffusion { static constexpr bool value = true; }; From 4ca3c2af7240e65199ac6d91c9348781e70b6263 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 18 Feb 2021 13:49:35 +0100 Subject: [PATCH 5/5] clean-up commented code --- ebos/eclproblem.hh | 1 - 1 file changed, 1 deletion(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index a5258d301..dfc62d689 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -3086,7 +3086,6 @@ private: { Opm::ConditionalStorage thermalHalfTrans; Opm::ConditionalStorage diffusivity; - //Scalar diffusivity; Scalar transmissibility; };