Only compute psat when needed #4515

totto82 · 2025-03-10T08:11:53Z

Note In Spycher Pruess 2010 they say 31C which is inconsistent with the critical temperature (30.95C) we use in OPM. For consistancy we change to critical temperature.

totto82 · 2025-03-10T08:12:05Z

jenkins build this please

totto82 · 2025-03-10T08:48:07Z

@svenn-t Can you review this?

svenn-t

Looks good to me, just a small comment.

svenn-t · 2025-03-10T09:16:55Z

opm/material/binarycoefficients/Brine_CO2.hpp

+            if (model1) {
+                // Computing the vapor pressure is not trivial and is also not defined for T > criticalTemperature
+                Evaluation psat = CO2::vaporPressure(temperature);
+                model1 = model1 && pg > psat;


This can be written model1 = pg > psat since model1 is always true inside this if-clause.

thanks. updated.

totto82 · 2025-03-10T11:58:42Z

jenkins build this please

blattms · 2025-03-10T21:17:14Z

This prevents the erroneous std::power call with negative base for the default value of ACTCO2S.
Unfortunately, for the other two possible values it is still called somewhere during the set up phase. I will try to find out where exactly.

blattms · 2025-03-10T22:00:38Z

Here is the relevant part of the stack trace when running ACTCO2S being 1 or 2 with this PR and my PR that throws:


Error: [opm-common/opm/material/components/CO2.hpp:150] Non-integer exponent for pow function of negative value (because temperature of CO2 is larger than critical temperature.

Thread 1 "flow" hit Catchpoint 3 (exception thrown), 0x00007ffff1ea90a1 in __cxa_throw () from /lib/x86_64-linux-gnu/libstdc++.so.6
(gdb) bt
#0  0x00007ffff1ea90a1 in __cxa_throw ()
   from /lib/x86_64-linux-gnu/libstdc++.so.6
#1  0x0000555555810a85 in Opm::CO2<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > >::vaporPressure<double> (T=<optimized out>)
    at opm-common/opm/material/components/CO2.hpp:150
#2  0x0000555555810cb3 in Opm::BinaryCoeff::Brine_CO2<double, Opm::SimpleHuDuanH2O<double>, Opm::CO2<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > >, true>::equilibriumConstantCO2_<double> (
    temperature=@0x7fffffffca80: 321.54618432792739, 
    pg=@0x7fffffffca88: 9969819.0459214486, highTemp=@0x7fffffffc8ef: false, 
    spycherPruess2005=spycherPruess2005@entry=false)
    at opm-common/opm/material/binarycoefficients/Brine_CO2.hpp:933
#3  0x0000555555810dab in Opm::BinaryCoeff::Brine_CO2<double, Opm::SimpleHuDuanH2O<double>, Opm::CO2<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > >, true>::computeB_<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > > (params=..., 
    temperature=@0x7fffffffca80: 321.54618432792739, 
    pg=@0x7fffffffca88: 9969819.0459214486, 
    yH2O=@0x7fffffffc900: 0.0011435699163673098, 
    xCO2=@0x7fffffffc8f8: 0.0089999999999999993, 
    highTemp=@0x7fffffffc8ef: false, extrapolate=true, spycherPruess2005=false)
    at opm-common/opm/material/binarycoefficients/Brine_CO2.hpp:772
#4  0x000055555581112c in Opm::BinaryCoeff::Brine_CO2<double, Opm::SimpleHuDuanH2O<double>, Opm::CO2<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > >, true>::mutualSolubility_<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > > (params=..., 
    temperature=@0x7fffffffca80: 321.54618432792739, 
    pg=@0x7fffffffca88: 9969819.0459214486, 
    xCO2=@0x7fffffffc8f8: 0.0089999999999999993, 
    yH2O=@0x7fffffffc900: 0.0011435699163673098, 
    m_NaCl=@0x7fffffffc978: 0.97999999999999998, 
    gammaNaCl=@0x7fffffffc8f0: 1.2535994027722603, 
    highTemp=@0x7fffffffc8ef: false, iterate=@0x7fffffffc8ee: true, 
    extrapolate=true)
    at opm-common/opm/material/binarycoefficients/Brine_CO2.hpp:655
#5  0x00005555558112a8 in Opm::BinaryCoeff::Brine_CO2<double, Opm::SimpleHuDuanH2O<double>, Opm::CO2<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > >, true>::fixPointIterSolubility_<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > > (params=..., 
    temperature=@0x7fffffffca80: 321.54618432792739, 
    pg=@0x7fffffffca88: 9969819.0459214486, 
    m_NaCl=@0x7fffffffc978: 0.97999999999999998, 
    activityModel=@0x555557704d24: 1, extrapolate=true)
    at opm-common/opm/material/binarycoefficients/Brine_CO2.hpp:591
#6  0x00005555558114f5 in Opm::BinaryCoeff::Brine_CO2<double, Opm::SimpleHuDuanH2O<double>, Opm::CO2<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > >, true>::calculateMoleFractions<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > > > (params=..., 
    temperature=@0x7fffffffca80: 321.54618432792739, 
    pg=@0x7fffffffca88: 9969819.0459214486, 
    salinity=@0x7fffffffca48: 0.054168883064250682, 
    knownPhaseIdx=knownPhaseIdx@entry=-1, xlCO2=@0x7fffffffc9e8: 0, 
    ygH2O=@0x7fffffffc9f0: 0, activityModel=@0x555557704d24: 1, 
    extrapolate=true)
    at opm-common/opm/material/binarycoefficients/Brine_CO2.hpp:137
#7  0x00005555558116bf in Opm::BrineCo2Pvt<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > >, std::vector<double, std::allocator<double> > >::rsSat<double> (this=0x555557704c90, regionIdx=regionIdx@entry=0, 
    temperature=@0x7fffffffca80: 321.54618432792739, 
    pressure=@0x7fffffffca88: 9969819.0459214486, 
    salinity=@0x7fffffffca48: 0.054168883064250682)
    at opm-common/opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp:536
#8  0x00005555558117a0 in Opm::BrineCo2Pvt<double, Opm::CO2Tables<double, std::vector<double, std::allocator<double> > >, std::vector<double, std::allocator<double> > >::saturatedGasDissolutionFactor<double> (this=<optimized out>, 
    regionIdx=regionIdx@entry=0, 
    temperature=@0x7fffffffca80: 321.54618432792739, 
    pressure=@0x7fffffffca88: 9969819.0459214486)
    at opm-common/opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp:436
#9  0x000055555581193d in Opm::OilPvtMultiplexer<double, true>::saturatedGasDissolutionFactor<double> (this=<optimized out>, regionIdx=0, 
    temperature=@0x7fffffffca80: 321.54618432792739, 
    pressure=@0x7fffffffca88: 9969819.0459214486)
    at opm-common/opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp:210
#10 0x00007ffff6515918 in Opm::EQUIL::Miscibility::RsVD<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr> >::satRs (this=this@entry=0x55556dfefb20, press=<optimized out>, 
    temp=<optimized out>)
    at opm-simulators/opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp:92
#11 0x00007ffff651593f in Opm::EQUIL::Miscibility::RsVD<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr> >::operator() (this=0x55556dfefb20, depth=<optimized out>, 
    press=<optimized out>, temp=<optimized out>, satGas=0)
    at opm-simulators/opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp:74
#12 0x00007ffff651f837 in Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> >::density (
    this=this@entry=0x7fffffffccf0, depth=<optimized out>, 
    press=<optimized out>)
    at opm-simulators/opm/simulators/flow/equil/InitStateEquil_impl.hpp:319
#13 0x00007ffff651f8eb in Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> >::operator() (this=this@entry=0x7fffffffccf0, depth=<optimized out>, press=<optimized out>)
    at opm-simulators/opm/simulators/flow/equil/InitStateEquil_impl.hpp:307
#14 0x00007ffff651f9c7 in Opm::EQUIL::Details::RK4IVP<double, Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> > >::RK4IVP (this=0x55556e6db380, f=..., span=..., 
    y0=<optimized out>, N=2000) at /usr/include/c++/12/array:206
#15 0x00007ffff651fbc4 in std::make_unique<Opm::EQUIL::Details::RK4IVP<double, Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> > >, Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> > const&, std::array<double, 2ul>, double const&, int const&> ()
    at /usr/include/c++/12/bits/unique_ptr.h:1065
#16 0x00007ffff651fc57 in Opm::EQUIL::Details::PressureTable<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::EquilReg<double> >::PressureFunction<Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> > >::PressureFunction (this=0x55561ef59370, 
    ode=..., ic=..., nsample=<optimized out>, span=...)
    at opm-simulators/opm/simulators/flow/equil/InitStateEquil_impl.hpp:443
#17 0x00007ffff651fd1b in std::make_unique<Opm::EQUIL::Details::PressureTable<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::EquilReg<double> >::PressureFunction<Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> > >, Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> > const&, Opm::EQUIL::Details::PressureTable<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::EquilReg<double> >::PressureFunction<Opm::EQUIL::Details::PhasePressODE::Oil<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::Miscibility::RsFunction<double> > >::InitCond const&, int&, std::array<double, 2ul> const&> () at /usr/include/c++/12/bits/unique_ptr.h:1065
#18 0x00007ffff651fdbb in Opm::EQUIL::Details::PressureTable<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::EquilReg<double> >::makeOilPressure (
    this=this@entry=0x7fffffffcf70, ic=..., reg=..., span=...)
    at opm-simulators/opm/simulators/flow/equil/InitStateEquil_impl.hpp:1256
#19 0x00007ffff6520d5e in Opm::EQUIL::Details::PressureTable<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::EquilReg<double> >::equil_OWG (
    this=0x7fffffffcf70, reg=..., span=...)
    at opm-simulators/opm/simulators/flow/equil/InitStateEquil_impl.hpp:1222
#20 0x00007ffff65166fc in Opm::EQUIL::Details::PressureTable<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Opm::EQUIL::EquilReg<double> >::equilibrate (
    this=this@entry=0x7fffffffcf70, reg=..., span=...)
    at opm-simulators/opm/simulators/flow/equil/InitStateEquil_impl.hpp:1049
#21 0x00007ffff651c603 in Opm::EQUIL::DeckDependent::InitialStateComputer<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Dune::CpGrid, Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid> >, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid> > >, Dune::CartesianIndexMapper<Dune::CpGrid> >::calcPressSatRsRv<Opm::RegionMapping<std::vector<int, std::allocator<int> > >, Opm::EclMaterialLawManager<Opm::ThreePhaseMaterialTraits<double, 0, 1, 2> >, Dune::Communication<ompi_communicator_t*> > (
    this=this@entry=0x7fffffffd620, reg=..., 
    rec=std::vector of length 1, capacity 1 = {...}, materialLawManager=..., 
    comm=..., grav=grav@entry=9.8066499999999994)
    at opm-simulators/opm/simulators/flow/equil/InitStateEquil_impl.hpp:1870
#22 0x00007ffff651ef7b in Opm::EQUIL::DeckDependent::InitialStateComputer<Opm::BlackOilFluidSystem<double, Opm::BlackOilDefaultIndexTraits, Opm::VectorWithDefaultAllocator, std::shared_ptr>, Dune::CpGrid, Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid> >, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid> > >, Dune::CartesianIndexMapper<Dune::CpGrid> >::InitialStateComputer<Opm::EclMaterialLawManager<Opm::ThreePhaseMaterialTraits<double, 0, 1, 2> > > (this=this@entry=0x7fffffffd620, 
    materialLawManager=..., eclipseState=..., grid=..., gridView=..., 
    cartMapper=..., grav=grav@entry=9.8066499999999994, 
    num_pressure_points=2000, applySwatInit=true)
    at opm-simulators/opm/simulators/flow/equil/InitStateEquil_impl.hpp:1569
#23 0x0000555555b9253b in Opm::EquilInitializer<Opm::Properties::TTag::FlowGasOilProblem>::EquilInitializer<Opm::EclMaterialLawManager<Opm::ThreePhaseMaterialTraits<double, 0, 1, 2> > > (this=0x7fffffffd840, simulator=..., 
    materialLawManager=...)
    at opm-simulators/opm/simulators/flow/EquilInitializer.hpp:119

totto82 mentioned this pull request Mar 10, 2025

CO2Store triggers division by zero for high temperatures (>304 K) #4512

Open

svenn-t reviewed Mar 10, 2025

View reviewed changes

Only compute psat when needed

48de096

totto82 force-pushed the guard_c branch from 7122f4a to 48de096 Compare March 10, 2025 11:58

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Only compute psat when needed #4515

Only compute psat when needed #4515

totto82 commented Mar 10, 2025

totto82 commented Mar 10, 2025

totto82 commented Mar 10, 2025

svenn-t left a comment

svenn-t Mar 10, 2025

totto82 Mar 10, 2025

totto82 commented Mar 10, 2025

blattms commented Mar 10, 2025

blattms commented Mar 10, 2025

Only compute psat when needed #4515

Are you sure you want to change the base?

Only compute psat when needed #4515

Conversation

totto82 commented Mar 10, 2025

totto82 commented Mar 10, 2025

totto82 commented Mar 10, 2025

svenn-t left a comment

Choose a reason for hiding this comment

svenn-t Mar 10, 2025

Choose a reason for hiding this comment

totto82 Mar 10, 2025

Choose a reason for hiding this comment

totto82 commented Mar 10, 2025

blattms commented Mar 10, 2025

blattms commented Mar 10, 2025