From 75a6b603c5196b5bc7e2824827abd2d547a68eae Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 12 Dec 2022 13:32:44 +0100 Subject: [PATCH 1/2] EclThermalLawManager: move code to separate compile unit --- CMakeLists_files.cmake | 1 + opm/material/thermal/EclThermalLawManager.hpp | 234 +-------------- .../material/thermal/EclThermalLawManager.cpp | 275 ++++++++++++++++++ 3 files changed, 291 insertions(+), 219 deletions(-) create mode 100644 src/opm/material/thermal/EclThermalLawManager.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index d8a25cfff..cca6bb40e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -269,6 +269,7 @@ if(ENABLE_ECL_INPUT) src/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.cpp src/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.cpp src/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.cpp + src/opm/material/thermal/EclThermalLawManager.cpp ) list( APPEND PYTHON_CXX_SOURCE_FILES diff --git a/opm/material/thermal/EclThermalLawManager.hpp b/opm/material/thermal/EclThermalLawManager.hpp index 4d5f4199a..1f0d3fa9f 100644 --- a/opm/material/thermal/EclThermalLawManager.hpp +++ b/opm/material/thermal/EclThermalLawManager.hpp @@ -37,15 +37,12 @@ #include "EclThermalConductionLawMultiplexer.hpp" #include "EclThermalConductionLawMultiplexerParams.hpp" -#include -#include - -#include -#include #include namespace Opm { +class EclipseState; + /*! * \ingroup fluidmatrixinteractions * @@ -64,249 +61,48 @@ public: using ThermalConductionLaw = EclThermalConductionLawMultiplexer; using ThermalConductionLawParams = typename ThermalConductionLaw::Params; - EclThermalLawManager() - { - solidEnergyApproach_ = SolidEnergyLawParams::undefinedApproach; - thermalConductivityApproach_ = ThermalConductionLawParams::undefinedApproach; - } + void initParamsForElements(const EclipseState& eclState, size_t numElems); - void initParamsForElements(const EclipseState& eclState, size_t numElems) - { - const auto& fp = eclState.fieldProps(); - const auto& tableManager = eclState.getTableManager(); - bool has_heatcr = fp.has_double("HEATCR"); - bool has_thconr = fp.has_double("THCONR"); - bool has_thc = fp.has_double("THCROCK") || fp.has_double("THCOIL") || fp.has_double("THCGAS") || fp.has_double("THCWATER"); + const SolidEnergyLawParams& solidEnergyLawParams(unsigned elemIdx) const; - if (has_heatcr) - initHeatcr_(eclState, numElems); - else if (tableManager.hasTables("SPECROCK")) - initSpecrock_(eclState, numElems); - else - initNullRockEnergy_(); - - if (has_thconr) - initThconr_(eclState, numElems); - else if (has_thc) - initThc_(eclState, numElems); - else - initNullCond_(); - } - - const SolidEnergyLawParams& solidEnergyLawParams(unsigned elemIdx) const - { - switch (solidEnergyApproach_) { - case SolidEnergyLawParams::heatcrApproach: - assert(elemIdx < solidEnergyLawParams_.size()); - return solidEnergyLawParams_[elemIdx]; - - case SolidEnergyLawParams::specrockApproach: - { - assert(elemIdx < elemToSatnumIdx_.size()); - unsigned satnumIdx = elemToSatnumIdx_[elemIdx]; - assert(satnumIdx < solidEnergyLawParams_.size()); - return solidEnergyLawParams_[satnumIdx]; - } - - case SolidEnergyLawParams::nullApproach: - return solidEnergyLawParams_[0]; - - default: - throw std::runtime_error("Attempting to retrieve solid energy storage parameters " - "without a known approach being defined by the deck."); - } - } - - const ThermalConductionLawParams& thermalConductionLawParams(unsigned elemIdx) const - { - switch (thermalConductivityApproach_) { - case ThermalConductionLawParams::thconrApproach: - case ThermalConductionLawParams::thcApproach: - assert(elemIdx < thermalConductionLawParams_.size()); - return thermalConductionLawParams_[elemIdx]; - - case ThermalConductionLawParams::nullApproach: - return thermalConductionLawParams_[0]; - - default: - throw std::runtime_error("Attempting to retrieve thermal conduction parameters without " - "a known approach being defined by the deck."); - } - } + const ThermalConductionLawParams& thermalConductionLawParams(unsigned elemIdx) const; private: /*! * \brief Initialize the parameters for the solid energy law using using HEATCR and friends. */ - void initHeatcr_(const EclipseState& eclState, - size_t numElems) - { - solidEnergyApproach_ = SolidEnergyLawParams::heatcrApproach; - // actually the value of the reference temperature does not matter for energy - // conservation. We set it anyway to faciliate comparisons with ECL - HeatcrLawParams::setReferenceTemperature(FluidSystem::surfaceTemperature); - - const auto& fp = eclState.fieldProps(); - 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(); - - heatcrElemParams.setReferenceRockHeatCapacity(heatcrData[elemIdx]); - heatcrElemParams.setDRockHeatCapacity_dT(heatcrtData[elemIdx]); - heatcrElemParams.finalize(); - elemParam.finalize(); - } - } + void initHeatcr_(const EclipseState& eclState, size_t numElems); /*! * \brief Initialize the parameters for the solid energy law using using SPECROCK and friends. */ - void initSpecrock_(const EclipseState& eclState, - size_t numElems) - { - solidEnergyApproach_ = SolidEnergyLawParams::specrockApproach; - - // initialize the element index -> SATNUM index mapping - const auto& fp = eclState.fieldProps(); - 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[elemIdx] - 1; - } - // internalize the SPECROCK table - unsigned numSatRegions = eclState.runspec().tabdims().getNumSatTables(); - const auto& tableManager = eclState.getTableManager(); - solidEnergyLawParams_.resize(numSatRegions); - for (unsigned satnumIdx = 0; satnumIdx < numSatRegions; ++satnumIdx) { - const auto& specrockTable = tableManager.getSpecrockTables()[satnumIdx]; - - auto& multiplexerParams = solidEnergyLawParams_[satnumIdx]; - - multiplexerParams.setSolidEnergyApproach(SolidEnergyLawParams::specrockApproach); - - auto& specrockParams = multiplexerParams.template getRealParams(); - const auto& temperatureColumn = specrockTable.getColumn("TEMPERATURE"); - const auto& cvRockColumn = specrockTable.getColumn("CV_ROCK"); - specrockParams.setHeatCapacities(temperatureColumn, cvRockColumn); - specrockParams.finalize(); - - multiplexerParams.finalize(); - } - } + void initSpecrock_(const EclipseState& eclState, size_t numElems); /*! * \brief Specify the solid energy law by setting heat capacity of rock to 0 */ - void initNullRockEnergy_() - { - solidEnergyApproach_ = SolidEnergyLawParams::nullApproach; - - solidEnergyLawParams_.resize(1); - solidEnergyLawParams_[0].finalize(); - } + void initNullRockEnergy_(); /*! * \brief Initialize the parameters for the thermal conduction law using THCONR and friends. */ - void initThconr_(const EclipseState& eclState, - size_t numElems) - { - thermalConductivityApproach_ = ThermalConductionLawParams::thconrApproach; - - const auto& fp = eclState.fieldProps(); - std::vector thconrData; - std::vector thconsfData; - if (fp.has_double("THCONR")) - thconrData = fp.get_double("THCONR"); - - if (fp.has_double("THCONSF")) - thconsfData = fp.get_double("THCONSF"); - - thermalConductionLawParams_.resize(numElems); - for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { - auto& elemParams = thermalConductionLawParams_[elemIdx]; - elemParams.setThermalConductionApproach(ThermalConductionLawParams::thconrApproach); - auto& thconrElemParams = elemParams.template getRealParams(); - - double thconr = thconrData.empty() ? 0.0 : thconrData[elemIdx]; - double thconsf = thconsfData.empty() ? 0.0 : thconsfData[elemIdx]; - thconrElemParams.setReferenceTotalThermalConductivity(thconr); - thconrElemParams.setDTotalThermalConductivity_dSg(thconsf); - - thconrElemParams.finalize(); - elemParams.finalize(); - } - } + void initThconr_(const EclipseState& eclState, size_t numElems); /*! * \brief Initialize the parameters for the thermal conduction law using THCROCK and friends. */ - void initThc_(const EclipseState& eclState, - size_t numElems) - { - thermalConductivityApproach_ = ThermalConductionLawParams::thcApproach; - - const auto& fp = eclState.fieldProps(); - std::vector thcrockData; - std::vector thcoilData; - std::vector thcgasData; - std::vector thcwaterData = fp.get_double("THCWATER"); - - if (fp.has_double("THCROCK")) - thcrockData = fp.get_double("THCROCK"); - - if (fp.has_double("THCOIL")) - thcoilData = fp.get_double("THCOIL"); - - if (fp.has_double("THCGAS")) - thcgasData = fp.get_double("THCGAS"); - - if (fp.has_double("THCWATER")) - thcwaterData = fp.get_double("THCWATER"); - - const std::vector& poroData = fp.get_double("PORO"); - - thermalConductionLawParams_.resize(numElems); - for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { - auto& elemParams = thermalConductionLawParams_[elemIdx]; - elemParams.setThermalConductionApproach(ThermalConductionLawParams::thcApproach); - auto& thcElemParams = elemParams.template getRealParams(); - - 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); - thcElemParams.setThcwater(thcwater); - - thcElemParams.finalize(); - elemParams.finalize(); - } - } + void initThc_(const EclipseState& eclState, size_t numElems); /*! * \brief Disable thermal conductivity */ - void initNullCond_() - { - thermalConductivityApproach_ = ThermalConductionLawParams::nullApproach; - - thermalConductionLawParams_.resize(1); - thermalConductionLawParams_[0].finalize(); - } + void initNullCond_(); private: - typename ThermalConductionLawParams::ThermalConductionApproach thermalConductivityApproach_; - typename SolidEnergyLawParams::SolidEnergyApproach solidEnergyApproach_; + using ConductionApproach = typename ThermalConductionLawParams::ThermalConductionApproach; + ConductionApproach thermalConductivityApproach_ = ThermalConductionLawParams::undefinedApproach; + using SolidEnergyApproach = typename SolidEnergyLawParams::SolidEnergyApproach; + SolidEnergyApproach solidEnergyApproach_ = SolidEnergyLawParams::undefinedApproach; std::vector elemToSatnumIdx_; diff --git a/src/opm/material/thermal/EclThermalLawManager.cpp b/src/opm/material/thermal/EclThermalLawManager.cpp new file mode 100644 index 000000000..3dd573d1b --- /dev/null +++ b/src/opm/material/thermal/EclThermalLawManager.cpp @@ -0,0 +1,275 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#include +#include + +#include +#include + +#include +#include + +#include +#include + +namespace Opm { + +template +void EclThermalLawManager:: +initParamsForElements(const EclipseState& eclState, size_t numElems) +{ + const auto& fp = eclState.fieldProps(); + const auto& tableManager = eclState.getTableManager(); + bool has_heatcr = fp.has_double("HEATCR"); + bool has_thconr = fp.has_double("THCONR"); + bool has_thc = fp.has_double("THCROCK") || + fp.has_double("THCOIL") || + fp.has_double("THCGAS") || + fp.has_double("THCWATER"); + + if (has_heatcr) + initHeatcr_(eclState, numElems); + else if (tableManager.hasTables("SPECROCK")) + initSpecrock_(eclState, numElems); + else + initNullRockEnergy_(); + + if (has_thconr) + initThconr_(eclState, numElems); + else if (has_thc) + initThc_(eclState, numElems); + else + initNullCond_(); +} + +template +const typename EclThermalLawManager::SolidEnergyLawParams& +EclThermalLawManager:: +solidEnergyLawParams(unsigned elemIdx) const +{ + switch (solidEnergyApproach_) { + case SolidEnergyLawParams::heatcrApproach: + assert(elemIdx < solidEnergyLawParams_.size()); + return solidEnergyLawParams_[elemIdx]; + + case SolidEnergyLawParams::specrockApproach: + { + assert(elemIdx < elemToSatnumIdx_.size()); + unsigned satnumIdx = elemToSatnumIdx_[elemIdx]; + assert(satnumIdx < solidEnergyLawParams_.size()); + return solidEnergyLawParams_[satnumIdx]; + } + + case SolidEnergyLawParams::nullApproach: + return solidEnergyLawParams_[0]; + + default: + throw std::runtime_error("Attempting to retrieve solid energy storage parameters " + "without a known approach being defined by the deck."); + } +} + +template +const typename EclThermalLawManager::ThermalConductionLawParams& +EclThermalLawManager:: +thermalConductionLawParams(unsigned elemIdx) const +{ + switch (thermalConductivityApproach_) { + case ThermalConductionLawParams::thconrApproach: + case ThermalConductionLawParams::thcApproach: + assert(elemIdx < thermalConductionLawParams_.size()); + return thermalConductionLawParams_[elemIdx]; + + case ThermalConductionLawParams::nullApproach: + return thermalConductionLawParams_[0]; + + default: + throw std::runtime_error("Attempting to retrieve thermal conduction parameters without " + "a known approach being defined by the deck."); + } +} + +template +void EclThermalLawManager:: +initHeatcr_(const EclipseState& eclState, size_t numElems) +{ + solidEnergyApproach_ = SolidEnergyLawParams::heatcrApproach; + // actually the value of the reference temperature does not matter for energy + // conservation. We set it anyway to faciliate comparisons with ECL + HeatcrLawParams::setReferenceTemperature(FluidSystem::surfaceTemperature); + + const auto& fp = eclState.fieldProps(); + 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(); + + heatcrElemParams.setReferenceRockHeatCapacity(heatcrData[elemIdx]); + heatcrElemParams.setDRockHeatCapacity_dT(heatcrtData[elemIdx]); + heatcrElemParams.finalize(); + elemParam.finalize(); + } +} + +template +void EclThermalLawManager:: +initSpecrock_(const EclipseState& eclState, size_t numElems) +{ + solidEnergyApproach_ = SolidEnergyLawParams::specrockApproach; + + // initialize the element index -> SATNUM index mapping + const auto& fp = eclState.fieldProps(); + 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[elemIdx] - 1; + } + // internalize the SPECROCK table + unsigned numSatRegions = eclState.runspec().tabdims().getNumSatTables(); + const auto& tableManager = eclState.getTableManager(); + solidEnergyLawParams_.resize(numSatRegions); + for (unsigned satnumIdx = 0; satnumIdx < numSatRegions; ++satnumIdx) { + const auto& specrockTable = tableManager.getSpecrockTables()[satnumIdx]; + + auto& multiplexerParams = solidEnergyLawParams_[satnumIdx]; + + multiplexerParams.setSolidEnergyApproach(SolidEnergyLawParams::specrockApproach); + + auto& specrockParams = multiplexerParams.template getRealParams(); + const auto& temperatureColumn = specrockTable.getColumn("TEMPERATURE"); + const auto& cvRockColumn = specrockTable.getColumn("CV_ROCK"); + specrockParams.setHeatCapacities(temperatureColumn, cvRockColumn); + specrockParams.finalize(); + + multiplexerParams.finalize(); + } +} + +template +void EclThermalLawManager:: +initNullRockEnergy_() +{ + solidEnergyApproach_ = SolidEnergyLawParams::nullApproach; + + solidEnergyLawParams_.resize(1); + solidEnergyLawParams_[0].finalize(); +} + +template +void EclThermalLawManager:: +initThconr_(const EclipseState& eclState, size_t numElems) +{ + thermalConductivityApproach_ = ThermalConductionLawParams::thconrApproach; + + const auto& fp = eclState.fieldProps(); + std::vector thconrData; + std::vector thconsfData; + if (fp.has_double("THCONR")) + thconrData = fp.get_double("THCONR"); + + if (fp.has_double("THCONSF")) + thconsfData = fp.get_double("THCONSF"); + + thermalConductionLawParams_.resize(numElems); + for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { + auto& elemParams = thermalConductionLawParams_[elemIdx]; + elemParams.setThermalConductionApproach(ThermalConductionLawParams::thconrApproach); + auto& thconrElemParams = elemParams.template getRealParams(); + + double thconr = thconrData.empty() ? 0.0 : thconrData[elemIdx]; + double thconsf = thconsfData.empty() ? 0.0 : thconsfData[elemIdx]; + thconrElemParams.setReferenceTotalThermalConductivity(thconr); + thconrElemParams.setDTotalThermalConductivity_dSg(thconsf); + + thconrElemParams.finalize(); + elemParams.finalize(); + } +} + +template +void EclThermalLawManager:: +initThc_(const EclipseState& eclState, size_t numElems) +{ + thermalConductivityApproach_ = ThermalConductionLawParams::thcApproach; + + const auto& fp = eclState.fieldProps(); + std::vector thcrockData; + std::vector thcoilData; + std::vector thcgasData; + std::vector thcwaterData = fp.get_double("THCWATER"); + + if (fp.has_double("THCROCK")) + thcrockData = fp.get_double("THCROCK"); + + if (fp.has_double("THCOIL")) + thcoilData = fp.get_double("THCOIL"); + + if (fp.has_double("THCGAS")) + thcgasData = fp.get_double("THCGAS"); + + if (fp.has_double("THCWATER")) + thcwaterData = fp.get_double("THCWATER"); + + const std::vector& poroData = fp.get_double("PORO"); + + thermalConductionLawParams_.resize(numElems); + for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { + auto& elemParams = thermalConductionLawParams_[elemIdx]; + elemParams.setThermalConductionApproach(ThermalConductionLawParams::thcApproach); + auto& thcElemParams = elemParams.template getRealParams(); + + 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); + thcElemParams.setThcwater(thcwater); + + thcElemParams.finalize(); + elemParams.finalize(); + } +} + +template +void EclThermalLawManager:: +initNullCond_() +{ + thermalConductivityApproach_ = ThermalConductionLawParams::nullApproach; + + thermalConductionLawParams_.resize(1); + thermalConductionLawParams_[0].finalize(); +} + +using FS = BlackOilFluidSystem; +template class EclThermalLawManager; + +} // namespace Opm From ee5b4240be92a8721c8a74c29825b939e76b0df6 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 22 Dec 2022 18:15:49 +0100 Subject: [PATCH 2/2] EclThermalLawManager: replace naked throws with OPM_THROW --- src/opm/material/thermal/EclThermalLawManager.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/opm/material/thermal/EclThermalLawManager.cpp b/src/opm/material/thermal/EclThermalLawManager.cpp index 3dd573d1b..430a80c9e 100644 --- a/src/opm/material/thermal/EclThermalLawManager.cpp +++ b/src/opm/material/thermal/EclThermalLawManager.cpp @@ -24,6 +24,8 @@ #include #include +#include + #include #include @@ -85,8 +87,9 @@ solidEnergyLawParams(unsigned elemIdx) const return solidEnergyLawParams_[0]; default: - throw std::runtime_error("Attempting to retrieve solid energy storage parameters " - "without a known approach being defined by the deck."); + OPM_THROW(std::runtime_error, + "Attempting to retrieve solid energy storage parameters " + "without a known approach being defined by the deck."); } } @@ -105,8 +108,9 @@ thermalConductionLawParams(unsigned elemIdx) const return thermalConductionLawParams_[0]; default: - throw std::runtime_error("Attempting to retrieve thermal conduction parameters without " - "a known approach being defined by the deck."); + OPM_THROW(std::runtime_error, + "Attempting to retrieve thermal conduction parameters without " + "a known approach being defined by the deck."); } }