From 562da9ef70614fcf527fd13abf409f57628c72d8 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 23 Sep 2024 11:51:53 +0200 Subject: [PATCH] support ZMF explicit initialization --- opm/simulators/flow/FlowProblemComp.hpp | 97 ++++++++++++++++++------- 1 file changed, 69 insertions(+), 28 deletions(-) diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index 13edcefce..1faff38a7 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -350,55 +350,66 @@ protected: const auto& vanguard = simulator.vanguard(); const auto& eclState = vanguard.eclState(); const auto& fp = eclState.fieldProps(); - // const bool has_xmf = fp.has_double("XMF"); - assert(fp.has_double("XMF")); - // const bool has_ymf = fp.has_double("YMF"); - assert(fp.has_double("YMF")); - const bool has_temp = fp.has_double("TEMPI"); const bool has_pressure = fp.has_double("PRESSURE"); + if (!has_pressure) + throw std::runtime_error("The ECL input file requires the presence of the PRESSURE " + "keyword if the model is initialized explicitly"); + + const bool has_xmf = fp.has_double("XMF"); + const bool has_ymf = fp.has_double("YMF"); + const bool has_zmf = fp.has_double("ZMF"); + if (! (has_zmf || (has_xmf && has_ymf))) { + throw std::runtime_error("The ECL input file requires the presence of ZMF or XMF and YMF " + "keyword if the model is initialized explicitly"); + } + + if (has_zmf && (has_xmf || has_ymf)) { + throw std::runtime_error("The ECL input file can not handle explicit initialization " + "with both ZMF and XMF or YMF"); + } + + if (has_xmf != has_ymf) { + throw std::runtime_error("The ECL input file needs XMF and YMF combined to do the explicit " + "initializtion when using XMF or YMF"); + } + + const bool has_temp = fp.has_double("TEMPI"); // const bool has_gas = fp.has_double("SGAS"); assert(fp.has_double("SGAS")); - if (!has_pressure) - throw std::runtime_error("The ECL input file requires the presence of the PRESSURE " - "keyword if the model is initialized explicitly"); std::size_t numDof = this->model().numGridDof(); initialFluidStates_.resize(numDof); - std::vector xmfData; - std::vector ymfData; std::vector waterSaturationData; std::vector gasSaturationData; std::vector soilData; std::vector pressureData; std::vector tempiData; - if (waterPhaseIdx > 0 && Indices::numPhases > 1) + const bool water_active = FluidSystem::phaseIsActive(waterPhaseIdx); + const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx); + const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx); + + if (water_active && Indices::numPhases > 1) waterSaturationData = fp.get_double("SWAT"); else waterSaturationData.resize(numDof); pressureData = fp.get_double("PRESSURE"); - assert(waterPhaseIdx < 0); - - xmfData = fp.get_double("XMF"); - ymfData = fp.get_double("YMF"); if (has_temp) { tempiData = fp.get_double("TEMPI"); } else { ; // TODO: throw? } - if (gasPhaseIdx > 0) // && FluidSystem::phaseIsActive(oilPhaseIdx)) + if (gas_active) // && FluidSystem::phaseIsActive(oilPhaseIdx)) gasSaturationData = fp.get_double("SGAS"); else gasSaturationData.resize(numDof); - // constexpr std::size_t num_comps = 3; - for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { auto& dofFluidState = initialFluidStates_[dofIdx]; // dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx)); @@ -407,11 +418,11 @@ protected: assert(std::isfinite(temperatureLoc) && temperatureLoc > 0); dofFluidState.setTemperature(temperatureLoc); - if (gasPhaseIdx > 0) { - dofFluidState.setSaturation(FluidSystem::gasPhaseIdx, - gasSaturationData[dofIdx]); + if (gas_active) { + dofFluidState.setSaturation(FluidSystem::gasPhaseIdx, + gasSaturationData[dofIdx]); } - if (oilPhaseIdx > 0) { + if (oil_active) { dofFluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - waterSaturationData[dofIdx] @@ -438,13 +449,43 @@ protected: dofFluidState.setPressure(phaseIdx, pressure); } - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - const std::size_t data_idx = compIdx * numDof + dofIdx; - const Scalar xmf = xmfData[data_idx]; - const Scalar ymf = ymfData[data_idx]; + if (has_xmf && has_ymf) { + const auto& xmfData = fp.get_double("XMF"); + const auto& ymfData = fp.get_double("YMF"); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + const std::size_t data_idx = compIdx * numDof + dofIdx; + const Scalar xmf = xmfData[data_idx]; + const Scalar ymf = ymfData[data_idx]; - dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf); - dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf); + dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf); + dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf); + } + } + + if (has_zmf) { + const auto& zmfData = fp.get_double("ZMF"); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + const std::size_t data_idx = compIdx * numDof + dofIdx; + const Scalar zmf = zmfData[data_idx]; + // TODO: I know we do this for co2_ptflash example, but maybe we should not do it this way? + if (gas_active) { + const auto& sat_gas = dofFluidState.saturation(gasPhaseIdx); + if (sat_gas > 0.) { + dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, zmf); + } else { + dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, 0.); + } + } + if (oil_active) { + dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, zmf); + const auto& sat_oil = dofFluidState.saturation(oilPhaseIdx); + if (sat_oil > 0.) { + dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, zmf); + } else { + dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, 0.); + } + } + } } } }