From efdb6a4dfc5672eb39a0c7a564d87b04afbaa238 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 28 Jan 2020 11:16:52 +0100 Subject: [PATCH 1/2] use compressed field properties in EclMaterialLawManager --- .../EclEpsGridProperties.hpp | 60 +++++++------------ .../EclMaterialLawManager.hpp | 18 +++--- tests/test_eclmateriallawmanager.cpp | 14 ++--- 3 files changed, 35 insertions(+), 57 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclEpsGridProperties.hpp b/opm/material/fluidmatrixinteractions/EclEpsGridProperties.hpp index cb8cb0286..4ade4c7fe 100644 --- a/opm/material/fluidmatrixinteractions/EclEpsGridProperties.hpp +++ b/opm/material/fluidmatrixinteractions/EclEpsGridProperties.hpp @@ -59,23 +59,10 @@ namespace Opm { */ namespace { -template -std::vector compressed_copy(const std::vector& global_vector, const std::vector& compressedToCartesianElemIdx) { - std::vector compressed = std::vector(compressedToCartesianElemIdx.size()); - for (std::size_t active_index = 0; active_index < compressedToCartesianElemIdx.size(); active_index++) { - auto global_index = compressedToCartesianElemIdx[active_index]; - compressed[active_index] = global_vector[global_index]; - } - - return compressed; -} - - - -std::vector try_get(const FieldPropsManager& fp, const std::string& keyword, const std::vector& compressedToCartesianElemIdx) { +std::vector try_get(const FieldPropsManager& fp, const std::string& keyword) { if (fp.has_double(keyword)) - return compressed_copy(fp.get_global_double(keyword), compressedToCartesianElemIdx); + return fp.get_double(keyword); return {}; } @@ -91,48 +78,47 @@ public: #if HAVE_ECL_INPUT EclEpsGridProperties(const Opm::EclipseState& eclState, - bool useImbibition, - const std::vector& compressedToCartesianElemIdx) + bool useImbibition) { std::string kwPrefix = useImbibition?"I":""; const auto& fp = eclState.fieldProps(); if (useImbibition) - compressed_satnum = compressed_copy(fp.get_global_int("IMBNUM"), compressedToCartesianElemIdx); + compressed_satnum = fp.get_int("IMBNUM"); else - compressed_satnum = compressed_copy(fp.get_global_int("SATNUM"), compressedToCartesianElemIdx); + compressed_satnum = fp.get_int("SATNUM"); - this->compressed_swl = try_get( fp, kwPrefix+"SWL", compressedToCartesianElemIdx); - this->compressed_sgl = try_get( fp, kwPrefix+"SGL", compressedToCartesianElemIdx); - this->compressed_swcr = try_get( fp, kwPrefix+"SWCR", compressedToCartesianElemIdx); - this->compressed_sgcr = try_get( fp, kwPrefix+"SGCR", compressedToCartesianElemIdx); - this->compressed_sowcr = try_get( fp, kwPrefix+"SOWCR", compressedToCartesianElemIdx); - this->compressed_sogcr = try_get( fp, kwPrefix+"SOGCR", compressedToCartesianElemIdx); - this->compressed_swu = try_get( fp, kwPrefix+"SWU", compressedToCartesianElemIdx); - this->compressed_sgu = try_get( fp, kwPrefix+"SGU", compressedToCartesianElemIdx); - this->compressed_pcw = try_get( fp, kwPrefix+"PCW", compressedToCartesianElemIdx); - this->compressed_pcg = try_get( fp, kwPrefix+"PCG", compressedToCartesianElemIdx); - this->compressed_krw = try_get( fp, kwPrefix+"KRW", compressedToCartesianElemIdx); - this->compressed_kro = try_get( fp, kwPrefix+"KRO", compressedToCartesianElemIdx); - this->compressed_krg = try_get( fp, kwPrefix+"KRG", compressedToCartesianElemIdx); + this->compressed_swl = try_get( fp, kwPrefix+"SWL"); + this->compressed_sgl = try_get( fp, kwPrefix+"SGL"); + this->compressed_swcr = try_get( fp, kwPrefix+"SWCR"); + this->compressed_sgcr = try_get( fp, kwPrefix+"SGCR"); + this->compressed_sowcr = try_get( fp, kwPrefix+"SOWCR"); + this->compressed_sogcr = try_get( fp, kwPrefix+"SOGCR"); + this->compressed_swu = try_get( fp, kwPrefix+"SWU"); + this->compressed_sgu = try_get( fp, kwPrefix+"SGU"); + this->compressed_pcw = try_get( fp, kwPrefix+"PCW"); + this->compressed_pcg = try_get( fp, kwPrefix+"PCG"); + this->compressed_krw = try_get( fp, kwPrefix+"KRW"); + this->compressed_kro = try_get( fp, kwPrefix+"KRO"); + this->compressed_krg = try_get( fp, kwPrefix+"KRG"); // _may_ be needed to calculate the Leverett capillary pressure scaling factor if (fp.has_double("PORO")) - this->compressed_poro = compressed_copy(fp.get_global_double("PORO"), compressedToCartesianElemIdx); + this->compressed_poro = fp.get_double("PORO"); if (fp.has_double("PERMX")) - this->compressed_permx = compressed_copy(fp.get_global_double("PERMX"), compressedToCartesianElemIdx); + this->compressed_permx = fp.get_double("PERMX"); else this->compressed_permx = std::vector(this->compressed_satnum.size()); if (fp.has_double("PERMY")) - this->compressed_permy= compressed_copy(fp.get_global_double("PERMY"), compressedToCartesianElemIdx); + this->compressed_permy = fp.get_double("PERMY"); else - this->compressed_permy= this->compressed_permx; + this->compressed_permy = this->compressed_permx; if (fp.has_double("PERMZ")) - this->compressed_permz= compressed_copy(fp.get_global_double("PERMZ"), compressedToCartesianElemIdx); + this->compressed_permz = fp.get_double("PERMZ"); else this->compressed_permz= this->compressed_permx; } diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp index 7e313bb4e..216b71a62 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp @@ -148,8 +148,7 @@ public: } } - void initParamsForElements(const EclipseState& eclState, - const std::vector& compressedToCartesianElemIdx) + void initParamsForElements(const EclipseState& eclState, size_t numCompressedElems) { // get the number of saturation regions const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables(); @@ -169,15 +168,13 @@ public: readOilWaterEffectiveParameters_(oilWaterEffectiveParamVector_, eclState, satRegionIdx); } - unsigned numCompressedElems = static_cast(compressedToCartesianElemIdx.size()); // copy the SATNUM grid property. in some cases this is not necessary, but it // should not require much memory anyway... satnumRegionArray_.resize(numCompressedElems); if (eclState.fieldProps().has_int("SATNUM")) { - const auto& satnumRawData = eclState.fieldProps().get_global_int("SATNUM"); + const auto& satnumRawData = eclState.fieldProps().get_int("SATNUM"); for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { - unsigned cartesianElemIdx = static_cast(compressedToCartesianElemIdx[elemIdx]); - satnumRegionArray_[elemIdx] = satnumRawData[cartesianElemIdx] - 1; + satnumRegionArray_[elemIdx] = satnumRawData[elemIdx] - 1; } } else @@ -187,10 +184,9 @@ public: // the same as the saturation region (SATNUM) imbnumRegionArray_ = satnumRegionArray_; if (eclState.fieldProps().has_int("IMBNUM")) { - const auto& imbnumRawData = eclState.fieldProps().get_global_int("IMBNUM"); + const auto& imbnumRawData = eclState.fieldProps().get_int("IMBNUM"); for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { - int cartesianElemIdx = compressedToCartesianElemIdx[elemIdx]; - imbnumRegionArray_[elemIdx] = imbnumRawData[cartesianElemIdx] - 1; + imbnumRegionArray_[elemIdx] = imbnumRawData[elemIdx] - 1; } } @@ -213,7 +209,7 @@ public: oilWaterScaledImbPointsVector.resize(numCompressedElems); } - EclEpsGridProperties epsGridProperties(eclState, false, compressedToCartesianElemIdx); + EclEpsGridProperties epsGridProperties(eclState, false); for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { readGasOilScaledPoints_(gasOilScaledInfoVector, @@ -233,7 +229,7 @@ public: } if (enableHysteresis()) { - EclEpsGridProperties epsImbGridProperties(eclState, true, compressedToCartesianElemIdx); + EclEpsGridProperties epsImbGridProperties(eclState, true); for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { readGasOilScaledPoints_(gasOilScaledImbInfoVector, gasOilScaledImbPointsVector, diff --git a/tests/test_eclmateriallawmanager.cpp b/tests/test_eclmateriallawmanager.cpp index 9e62f5e2e..2ad2ebf1a 100644 --- a/tests/test_eclmateriallawmanager.cpp +++ b/tests/test_eclmateriallawmanager.cpp @@ -447,14 +447,10 @@ inline void testAll() const auto& eclGrid = eclState.getInputGrid(); size_t n = eclGrid.getCartesianSize(); - std::vector compressedToCartesianIdx(n); - - for (size_t i = 0; i < n; ++ i) - compressedToCartesianIdx[i] = static_cast(i); MaterialLawManager materialLawManager; materialLawManager.initFromDeck(deck, eclState); - materialLawManager.initParamsForElements(eclState, compressedToCartesianIdx); + materialLawManager.initParamsForElements(eclState, n); if (materialLawManager.enableEndPointScaling()) throw std::logic_error("Discrepancy between the deck and the EclMaterialLawManager"); @@ -468,7 +464,7 @@ inline void testAll() Opm::EclMaterialLawManager fam2MaterialLawManager; fam2MaterialLawManager.initFromDeck(fam2Deck, fam2EclState); - fam2MaterialLawManager.initParamsForElements(fam2EclState, compressedToCartesianIdx); + fam2MaterialLawManager.initParamsForElements(fam2EclState, n); if (fam2MaterialLawManager.enableEndPointScaling()) throw std::logic_error("Discrepancy between the deck and the EclMaterialLawManager"); @@ -481,7 +477,7 @@ inline void testAll() Opm::EclMaterialLawManager hysterMaterialLawManager; hysterMaterialLawManager.initFromDeck(hysterDeck, hysterEclState); - hysterMaterialLawManager.initParamsForElements(hysterEclState, compressedToCartesianIdx); + hysterMaterialLawManager.initParamsForElements(hysterEclState, n); if (hysterMaterialLawManager.enableEndPointScaling()) throw std::logic_error("Discrepancy between the deck and the EclMaterialLawManager"); @@ -574,14 +570,14 @@ inline void testAll() MaterialLawManager fam1materialLawManager; fam1materialLawManager.initFromDeck(fam1Deck, fam1EclState); - fam1materialLawManager.initParamsForElements(fam1EclState, compressedToCartesianIdx); + fam1materialLawManager.initParamsForElements(fam1EclState, n); const auto fam2Deck = parser.parseString(fam2DeckStringGasOil); const Opm::EclipseState fam2EclState(fam2Deck); Opm::EclMaterialLawManager fam2MaterialLawManager; fam2MaterialLawManager.initFromDeck(fam2Deck, fam2EclState); - fam2MaterialLawManager.initParamsForElements(fam2EclState, compressedToCartesianIdx); + fam2MaterialLawManager.initParamsForElements(fam2EclState, n); for (unsigned elemIdx = 0; elemIdx < n; ++ elemIdx) { for (int i = 0; i < 100; ++ i) { From ac9207cba14bbe2b1f75959f1adac0862d72973c Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 28 Jan 2020 11:16:52 +0100 Subject: [PATCH 2/2] use compressed field properties in EclThermalLawManager --- opm/material/thermal/EclThermalLawManager.hpp | 73 ++++++++----------- 1 file changed, 32 insertions(+), 41 deletions(-) diff --git a/opm/material/thermal/EclThermalLawManager.hpp b/opm/material/thermal/EclThermalLawManager.hpp index 26299c165..6b290e05a 100644 --- a/opm/material/thermal/EclThermalLawManager.hpp +++ b/opm/material/thermal/EclThermalLawManager.hpp @@ -69,8 +69,7 @@ public: thermalConductivityApproach_ = ThermalConductionLawParams::undefinedApproach; } - void initParamsForElements(const Opm::EclipseState& eclState, - const std::vector& compressedToCartesianElemIdx) + void initParamsForElements(const Opm::EclipseState& eclState, size_t numElems) { const auto& fp = eclState.fieldProps(); const auto& tableManager = eclState.getTableManager(); @@ -79,16 +78,16 @@ public: bool has_thc = fp.has_double("THCROCK") || fp.has_double("THCOIL") || fp.has_double("THCGAS") || fp.has_double("THCWATER"); if (has_heatcr) - initHeatcr_(eclState, compressedToCartesianElemIdx); + initHeatcr_(eclState, numElems); else if (tableManager.hasTables("SPECROCK")) - initSpecrock_(eclState, compressedToCartesianElemIdx); + initSpecrock_(eclState, numElems); else initNullRockEnergy_(); if (has_thconr) - initThconr_(eclState, compressedToCartesianElemIdx); + initThconr_(eclState, numElems); else if (has_thc) - initThc_(eclState, compressedToCartesianElemIdx); + initThc_(eclState, numElems); else initNullCond_(); } @@ -139,7 +138,7 @@ private: * \brief Initialize the parameters for the solid energy law using using HEATCR and friends. */ void initHeatcr_(const Opm::EclipseState& eclState, - const std::vector& compressedToCartesianElemIdx) + size_t numElems) { solidEnergyApproach_ = SolidEnergyLawParams::heatcrApproach; // actually the value of the reference temperature does not matter for energy @@ -147,18 +146,16 @@ private: HeatcrLawParams::setReferenceTemperature(FluidSystem::surfaceTemperature); const auto& fp = eclState.fieldProps(); - const std::vector& heatcrData = fp.get_global_double("HEATCR"); - const std::vector& heatcrtData = fp.get_global_double("HEATCRT"); - unsigned numElems = compressedToCartesianElemIdx.size(); + const std::vector& heatcrData = fp.get_double("HEATCR"); + const std::vector& heatcrtData = fp.get_double("HEATCRT"); solidEnergyLawParams_.resize(numElems); for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { auto& elemParam = solidEnergyLawParams_[elemIdx]; elemParam.setSolidEnergyApproach(SolidEnergyLawParams::heatcrApproach); auto& heatcrElemParams = elemParam.template getRealParams(); - unsigned cartesianElemIdx = compressedToCartesianElemIdx[elemIdx]; - heatcrElemParams.setReferenceRockHeatCapacity(heatcrData[cartesianElemIdx]); - heatcrElemParams.setDRockHeatCapacity_dT(heatcrtData[cartesianElemIdx]); + heatcrElemParams.setReferenceRockHeatCapacity(heatcrData[elemIdx]); + heatcrElemParams.setDRockHeatCapacity_dT(heatcrtData[elemIdx]); heatcrElemParams.finalize(); elemParam.finalize(); } @@ -168,20 +165,18 @@ private: * \brief Initialize the parameters for the solid energy law using using SPECROCK and friends. */ void initSpecrock_(const Opm::EclipseState& eclState, - const std::vector& compressedToCartesianElemIdx) + size_t numElems) { solidEnergyApproach_ = SolidEnergyLawParams::specrockApproach; // initialize the element index -> SATNUM index mapping const auto& fp = eclState.fieldProps(); - const std::vector& satnumData = fp.get_global_int("SATNUM"); - elemToSatnumIdx_.resize(compressedToCartesianElemIdx.size()); - for (unsigned elemIdx = 0; elemIdx < compressedToCartesianElemIdx.size(); ++ elemIdx) { - unsigned cartesianElemIdx = compressedToCartesianElemIdx[elemIdx]; - + const std::vector& satnumData = fp.get_int("SATNUM"); + elemToSatnumIdx_.resize(numElems); + for (unsigned elemIdx = 0; elemIdx < numElems; ++ elemIdx) { // satnumData contains Fortran-style indices, i.e., they start with 1 instead // of 0! - elemToSatnumIdx_[elemIdx] = satnumData[cartesianElemIdx] - 1; + elemToSatnumIdx_[elemIdx] = satnumData[elemIdx] - 1; } // internalize the SPECROCK table unsigned numSatRegions = eclState.runspec().tabdims().getNumSatTables(); @@ -219,7 +214,7 @@ private: * \brief Initialize the parameters for the thermal conduction law using THCONR and friends. */ void initThconr_(const Opm::EclipseState& eclState, - const std::vector& compressedToCartesianElemIdx) + size_t numElems) { thermalConductivityApproach_ = ThermalConductionLawParams::thconrApproach; @@ -227,21 +222,19 @@ private: std::vector thconrData; std::vector thconsfData; if (fp.has_double("THCONR")) - thconrData = fp.get_global_double("THCONR"); + thconrData = fp.get_double("THCONR"); if (fp.has_double("THCONSF")) - thconsfData = fp.get_global_double("THCONSF"); + thconsfData = fp.get_double("THCONSF"); - unsigned numElems = compressedToCartesianElemIdx.size(); thermalConductionLawParams_.resize(numElems); for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { auto& elemParams = thermalConductionLawParams_[elemIdx]; elemParams.setThermalConductionApproach(ThermalConductionLawParams::thconrApproach); auto& thconrElemParams = elemParams.template getRealParams(); - int cartElemIdx = compressedToCartesianElemIdx[elemIdx]; - double thconr = thconrData.empty() ? 0.0 : thconrData[cartElemIdx]; - double thconsf = thconsfData.empty() ? 0.0 : thconsfData[cartElemIdx]; + double thconr = thconrData.empty() ? 0.0 : thconrData[elemIdx]; + double thconsf = thconsfData.empty() ? 0.0 : thconsfData[elemIdx]; thconrElemParams.setReferenceTotalThermalConductivity(thconr); thconrElemParams.setDTotalThermalConductivity_dSg(thconsf); @@ -254,7 +247,7 @@ private: * \brief Initialize the parameters for the thermal conduction law using THCROCK and friends. */ void initThc_(const Opm::EclipseState& eclState, - const std::vector& compressedToCartesianElemIdx) + size_t numElems) { thermalConductivityApproach_ = ThermalConductionLawParams::thcApproach; @@ -262,35 +255,33 @@ private: std::vector thcrockData; std::vector thcoilData; std::vector thcgasData; - std::vector thcwaterData = fp.get_global_double("THCWATER"); + std::vector thcwaterData = fp.get_double("THCWATER"); if (fp.has_double("THCROCK")) - thcrockData = fp.get_global_double("THCROCK"); + thcrockData = fp.get_double("THCROCK"); if (fp.has_double("THCOIL")) - thcoilData = fp.get_global_double("THCOIL"); + thcoilData = fp.get_double("THCOIL"); if (fp.has_double("THCGAS")) - thcgasData = fp.get_global_double("THCGAS"); + thcgasData = fp.get_double("THCGAS"); if (fp.has_double("THCWATER")) - thcwaterData = fp.get_global_double("THCWATER"); + thcwaterData = fp.get_double("THCWATER"); - const std::vector& poroData = fp.get_global_double("PORO"); + const std::vector& poroData = fp.get_double("PORO"); - unsigned numElems = compressedToCartesianElemIdx.size(); thermalConductionLawParams_.resize(numElems); for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { auto& elemParams = thermalConductionLawParams_[elemIdx]; elemParams.setThermalConductionApproach(ThermalConductionLawParams::thcApproach); auto& thcElemParams = elemParams.template getRealParams(); - int cartElemIdx = compressedToCartesianElemIdx[elemIdx]; - thcElemParams.setPorosity(poroData[cartElemIdx]); - double thcrock = thcrockData.empty() ? 0.0 : thcrockData[cartElemIdx]; - double thcoil = thcoilData.empty() ? 0.0 : thcoilData[cartElemIdx]; - double thcgas = thcgasData.empty() ? 0.0 : thcgasData[cartElemIdx]; - double thcwater = thcwaterData.empty() ? 0.0 : thcwaterData[cartElemIdx]; + thcElemParams.setPorosity(poroData[elemIdx]); + double thcrock = thcrockData.empty() ? 0.0 : thcrockData[elemIdx]; + double thcoil = thcoilData.empty() ? 0.0 : thcoilData[elemIdx]; + double thcgas = thcgasData.empty() ? 0.0 : thcgasData[elemIdx]; + double thcwater = thcwaterData.empty() ? 0.0 : thcwaterData[elemIdx]; thcElemParams.setThcrock(thcrock); thcElemParams.setThcoil(thcoil); thcElemParams.setThcgas(thcgas);