diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 2c2c0912265..0fef4368264 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3080-5511-8210a1a + baseline: integratedTests/baseline_integratedTests-pr3140-5514-fe26f80 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 8e8eed2bde5..eff45f70b87 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3140 (2024-06-11) +====================== +Fixed derivative in EzrokhiBrineDensity + PR #3080 (2024-06-07) ===================== Rebaseline after adding viscoelastic wave propagator. diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp index a61f15bbe52..907ca230af5 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp @@ -194,11 +194,13 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::checkTablesParameters( real64 const template< typename PHASE1, typename PHASE2, typename FLASH > void CO2BrineFluid< PHASE1, PHASE2, FLASH >::initializePreSubGroups() { +#if defined(GEOS_DEVICE_COMPILE) GEOS_THROW_IF( this->getCatalogName() == CO2BrineEzrokhiThermalFluid::catalogName(), GEOS_FMT( "The `{}` model is disabled for now. Please use the other thermal CO2-brine model instead: `{}`", CO2BrineEzrokhiThermalFluid::catalogName(), CO2BrinePhillipsThermalFluid::catalogName() ), InputError ); +#endif } template< typename PHASE1, typename PHASE2, typename FLASH > diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/EzrokhiBrineDensity.hpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/EzrokhiBrineDensity.hpp index 0f66d19ba4e..bffb350a4ac 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/EzrokhiBrineDensity.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/EzrokhiBrineDensity.hpp @@ -171,6 +171,7 @@ void EzrokhiBrineDensityUpdate::compute( real64 const & pressure, arraySlice1d< real64, USD3 > const & dValue, bool useMass ) const { + constexpr integer numDof = 4; using Deriv = multifluid::DerivativeOffset; real64 waterSatDensity_dTemperature = 0.0; @@ -188,32 +189,39 @@ void EzrokhiBrineDensityUpdate::compute( real64 const & pressure, real64 const coefPhaseComposition = m_coef0 + temperature * ( m_coef1 + m_coef2 * temperature ); + real64 const mw_co2 = m_componentMolarWeight[m_CO2Index]; + real64 const mw_h2o = m_componentMolarWeight[m_waterIndex]; + // we have to convert molar component phase fraction (phaseComposition[m_CO2Index]) to mass fraction - real64 const waterMWInv = 1.0 / (phaseComposition[m_CO2Index] * m_componentMolarWeight[m_CO2Index] - + phaseComposition[m_waterIndex] * m_componentMolarWeight[m_waterIndex]); + real64 const waterMWInv = 1.0 / (phaseComposition[m_CO2Index] * mw_co2 + phaseComposition[m_waterIndex] * mw_h2o); + real64 dWaterMWInv[numDof]{}; + for( integer dof = 0; dof < numDof; ++dof ) + { + dWaterMWInv[dof] = -(dPhaseComposition[m_CO2Index][dof] * mw_co2 + dPhaseComposition[m_waterIndex][dof] * mw_h2o)*waterMWInv*waterMWInv; + } - real64 const massPhaseCompositionCO2 = phaseComposition[m_CO2Index] * m_componentMolarWeight[m_CO2Index] * waterMWInv; + real64 const massPhaseCompositionCO2 = phaseComposition[m_CO2Index] * mw_co2 * waterMWInv; real64 const exponent = coefPhaseComposition * massPhaseCompositionCO2; + real64 dExponent[numDof]{}; + for( integer dof = 0; dof < numDof; ++dof ) + { + dExponent[dof] = coefPhaseComposition * mw_co2 * + (dPhaseComposition[m_CO2Index][dof] * waterMWInv + phaseComposition[m_CO2Index] * dWaterMWInv[dof]); + } + dExponent[Deriv::dT] += ( m_coef1 + 2.0 * m_coef2 * temperature) * massPhaseCompositionCO2; - real64 const exponent_dPressure = coefPhaseComposition * dPhaseComposition[m_CO2Index][Deriv::dP]; - real64 const exponent_dTemperature = coefPhaseComposition * dPhaseComposition[m_CO2Index][Deriv::dT] + - ( m_coef1 + 2 * m_coef2 * temperature) * massPhaseCompositionCO2; - // compute only common part of derivatives w.r.t. CO2 and water phase compositions - // later to be multiplied by (phaseComposition[m_waterIndex]) and ( -phaseComposition[m_CO2Index] ) respectively - real64 const exponent_dPhaseComp = coefPhaseComposition * m_componentMolarWeight[m_CO2Index] * m_componentMolarWeight[m_waterIndex] * waterMWInv * waterMWInv; - real64 exponentPowered = pow( 10, exponent ); + real64 exponentPowered = pow( 10.0, exponent ); value = waterDensity * exponentPowered; - real64 const dValueCoef = LvArray::math::log( 10 ) * value; - dValue[Deriv::dP] = dValueCoef * exponent_dPressure + waterDensity_dPressure * exponentPowered; - dValue[Deriv::dT] = dValueCoef * exponent_dTemperature + waterDensity_dTemperature * exponentPowered; + real64 const dValueCoef = LvArray::math::log( 10.0 ) * value; + + dValue[Deriv::dP] = dValueCoef * dExponent[Deriv::dP] + waterDensity_dPressure * exponentPowered; + dValue[Deriv::dT] = dValueCoef * dExponent[Deriv::dT] + waterDensity_dTemperature * exponentPowered; - // here, we multiply common part of derivatives by specific coefficients - real64 const dValue_dPhaseComp = dValueCoef * exponent_dPhaseComp; - dValue[Deriv::dC+m_CO2Index] = dValue_dPhaseComp * phaseComposition[m_waterIndex] * dPhaseComposition[m_CO2Index][Deriv::dC+m_CO2Index]; - dValue[Deriv::dC+m_waterIndex] = dValue_dPhaseComp * ( -phaseComposition[m_CO2Index] ) * dPhaseComposition[m_waterIndex][Deriv::dC+m_waterIndex]; + dValue[Deriv::dC+m_CO2Index] = dValueCoef * dExponent[Deriv::dC+m_CO2Index]; + dValue[Deriv::dC+m_waterIndex] = dValueCoef * dExponent[Deriv::dC+m_waterIndex]; if( !useMass ) { diff --git a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp index 83607986b96..79a1682cc1a 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp @@ -47,6 +47,8 @@ void constitutiveUpdatePassThru( MultiFluidBase const & fluid, CO2BrinePhillipsFluid, CO2BrineEzrokhiFluid, CO2BrinePhillipsThermalFluid, +// Including these in a CUDA build will lead to compiler segfault. +// Need to split compilation units for all the options #if !defined(GEOS_DEVICE_COMPILE) CO2BrineEzrokhiThermalFluid, CompositionalTwoPhasePengRobinsonLBCViscosity, @@ -69,6 +71,8 @@ void constitutiveUpdatePassThru( MultiFluidBase & fluid, CO2BrinePhillipsFluid, CO2BrineEzrokhiFluid, CO2BrinePhillipsThermalFluid, +// Including these in a CUDA build will lead to compiler segfault. +// Need to split compilation units for all the options" #if !defined(GEOS_DEVICE_COMPILE) CO2BrineEzrokhiThermalFluid, CompositionalTwoPhasePengRobinsonLBCViscosity, diff --git a/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt b/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt index 5f03cf2eaaa..44890ac660a 100644 --- a/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt @@ -4,6 +4,9 @@ set( gtest_geosx_tests testCO2BrinePVTModels.cpp testCO2SpycherPruessModels.cpp testDamage.cpp + testMultiFluidCO2Brine.cpp + testMultiFluidDeadOil.cpp + testMultiFluidLiveOil.cpp testRelPerm.cpp testRelPermHysteresis.cpp ) @@ -39,8 +42,7 @@ endif() if( ENABLE_PVTPackage ) list( APPEND gtest_geosx_tests - testMultiFluid.cpp ) - + testMultiFluidCompositionalMultiphasePVTPackage.cpp ) list( APPEND dependencyList PVTPackage ) endif() diff --git a/src/coreComponents/unitTests/constitutiveTests/MultiFluidTest.hpp b/src/coreComponents/unitTests/constitutiveTests/MultiFluidTest.hpp new file mode 100644 index 00000000000..b6f8d7dde47 --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/MultiFluidTest.hpp @@ -0,0 +1,529 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultiFluidTest.hpp + */ + +#ifndef GEOS_UNITTESTS_CONSTITUTIVETESTS_MULTIFLUIDTEST_HPP_ +#define GEOS_UNITTESTS_CONSTITUTIVETESTS_MULTIFLUIDTEST_HPP_ + +#include "constitutiveTestHelpers.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "constitutive/fluid/multifluid/MultiFluidFields.hpp" +#include "constitutive/fluid/multifluid/MultiFluidSelector.hpp" + +using namespace geos::dataRepository; +using namespace geos::constitutive; +using namespace geos::constitutive::multifluid; + +namespace geos +{ +namespace testing +{ + +template< typename FLUID_TYPE > +struct UsePVTPackage +{ + static constexpr bool value = false; +}; + +template< integer N > +using Feed = std::array< real64, N >; + +template< integer N, typename ARRAY1, typename ARRAY2 > +void copy( ARRAY1 const & inputArray, ARRAY2 & outputArray ) +{ + for( integer i = 0; i < N; ++i ) + { + outputArray[i] = inputArray[i]; + } +} + +template< integer N, typename T > +void fill( array1d< T > & outputArray, std::array< T const, N > const inputArray ) +{ + for( integer i = 0; i < N; ++i ) + { + outputArray.emplace_back( inputArray[i] ); + } +} + +template< integer N, integer U, typename T > +void fill( arraySlice1d< T, U > const & outputArray, std::array< T const, N > const inputArray ) +{ + for( integer i = 0; i < N; ++i ) + { + outputArray[i] = inputArray[i]; + } +} + +template< integer NUM_PHASE, integer NUM_COMP > +struct MultiFluidTestData +{ + MultiFluidTestData( real64 const pressure_, + real64 const temperature_, + Feed< NUM_COMP > && composition_ ); + + MultiFluidTestData( real64 const pressure_, + real64 const temperature_, + arraySlice1d< real64 const > const & composition_ ); + + MultiFluidTestData( real64 const pressure_, + real64 const temperature_, + Feed< NUM_COMP > && composition_, + real64 const totalDensity_, + Feed< NUM_PHASE > && phaseFraction_, + Feed< NUM_PHASE > && phaseDensity_, + Feed< NUM_PHASE > && phaseMassDensity_, + Feed< NUM_PHASE > && phaseViscosity_, + Feed< NUM_PHASE > && phaseEnthalpy_, + Feed< NUM_PHASE > && phaseInternalEnergy_ ); + + real64 const pressure{0.0}; + real64 const temperature{0.0}; + Feed< NUM_COMP > composition{}; + real64 const totalDensity{0.0}; + Feed< NUM_PHASE > phaseFraction{}; + Feed< NUM_PHASE > phaseDensity{}; + Feed< NUM_PHASE > phaseMassDensity{}; + Feed< NUM_PHASE > phaseViscosity{}; + Feed< NUM_PHASE > phaseEnthalpy{}; + Feed< NUM_PHASE > phaseInternalEnergy{}; +}; + +template< typename FLUID_TYPE, integer NUM_PHASE, integer NUM_COMP > +class MultiFluidTest : public ConstitutiveTestBase< MultiFluidBase > +{ +public: + static constexpr integer numPhase = NUM_PHASE; + static constexpr integer numComp = NUM_COMP; + using TestData = MultiFluidTestData< numPhase, numComp >; + +public: + MultiFluidTest() = default; + ~MultiFluidTest() override = default; + + MultiFluidBase & getFluid() const { return *m_model; } + + Group & getParent() { return m_parent; } + + void testValuesAgainstPreviousImplementation( typename FLUID_TYPE::KernelWrapper const & wrapper, + MultiFluidTestData< NUM_PHASE, NUM_COMP > const & testData, + real64 const relTol ) const; + + void testNumericalDerivatives( MultiFluidBase & fluid, + Group * parent, + MultiFluidTestData< NUM_PHASE, NUM_COMP > const & testData, + real64 const perturbParameter, + real64 const relTol, + real64 const absTol = std::numeric_limits< real64 >::max() ); + + +protected: + virtual void resetFluid( MultiFluidBase & fluid ) const + { + GEOS_UNUSED_VAR( fluid ); + } + + static void writeTableToFile( string const & fileName, char const * content ) + { + std::ofstream os( fileName ); + ASSERT_TRUE( os.is_open() ); + os << content; + os.close(); + } + + static void removeFile( string const & fileName ) + { + int const ret = std::remove( fileName.c_str() ); + ASSERT_TRUE( ret == 0 ); + } +}; + +template< typename FLUID_TYPE, integer NUM_PHASE, integer NUM_COMP > +void MultiFluidTest< FLUID_TYPE, NUM_PHASE, NUM_COMP >:: +testNumericalDerivatives( MultiFluidBase & fluid, + Group * parent, + MultiFluidTestData< NUM_PHASE, NUM_COMP > const & testData, + real64 const perturbParameter, + real64 const relTol, + real64 const absTol ) +{ + using Deriv = multifluid::DerivativeOffset; + + integer const NC = fluid.numFluidComponents(); + integer const NP = fluid.numFluidPhases(); + integer const NDOF = NC+2; + + bool const isThermal = fluid.isThermal(); + + // Copy input values into an array with expected layout + array2d< real64, compflow::LAYOUT_COMP > compositionValues( 1, NC ); + for( integer i = 0; i < NC; ++i ) + { + compositionValues[0][i] = testData.composition[i]; + } + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const composition = compositionValues[0]; + + auto const & components = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); + auto const & phases = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); + + // create a clone of the fluid to run updates on + string const fluidCopyName = fluid.getName() + "Copy"; + std::unique_ptr< ConstitutiveBase > fluidCopyPtr = fluid.deliverClone( fluidCopyName, parent ); + MultiFluidBase & fluidCopy = dynamicCast< MultiFluidBase & >( *fluidCopyPtr ); + + fluid.allocateConstitutiveData( fluid.getParent(), 1 ); + fluidCopy.allocateConstitutiveData( fluid.getParent(), 1 ); + + // extract data views from both fluids + #define GET_FLUID_DATA( FLUID, TRAIT ) \ + FLUID.getReference< TRAIT::type >( TRAIT::key() )[0][0] + + MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > phaseFrac { + GET_FLUID_DATA( fluid, fields::multifluid::phaseFraction ), + GET_FLUID_DATA( fluid, fields::multifluid::dPhaseFraction ) + }; + + MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > phaseDens { + GET_FLUID_DATA( fluid, fields::multifluid::phaseDensity ), + GET_FLUID_DATA( fluid, fields::multifluid::dPhaseDensity ) + }; + + MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > phaseVisc { + GET_FLUID_DATA( fluid, fields::multifluid::phaseViscosity ), + GET_FLUID_DATA( fluid, fields::multifluid::dPhaseViscosity ) + }; + + MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > phaseEnthalpy { + GET_FLUID_DATA( fluid, fields::multifluid::phaseEnthalpy ), + GET_FLUID_DATA( fluid, fields::multifluid::dPhaseEnthalpy ) + }; + + MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > phaseInternalEnergy { + GET_FLUID_DATA( fluid, fields::multifluid::phaseInternalEnergy ), + GET_FLUID_DATA( fluid, fields::multifluid::dPhaseInternalEnergy ) + }; + + MultiFluidVarSlice< real64, 2, USD_PHASE_COMP - 2, USD_PHASE_COMP_DC - 2 > phaseCompFrac { + GET_FLUID_DATA( fluid, fields::multifluid::phaseCompFraction ), + GET_FLUID_DATA( fluid, fields::multifluid::dPhaseCompFraction ) + }; + + MultiFluidVarSlice< real64, 0, USD_FLUID - 2, USD_FLUID_DC - 2 > totalDens { + GET_FLUID_DATA( fluid, fields::multifluid::totalDensity ), + GET_FLUID_DATA( fluid, fields::multifluid::dTotalDensity ) + }; + + auto const & phaseFracCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseFraction ); + auto const & phaseDensCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseDensity ); + auto const & phaseViscCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseViscosity ); + auto const & phaseEnthCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseEnthalpy ); + auto const & phaseEnergyCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseInternalEnergy ); + auto const & phaseCompFracCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseCompFraction ); + auto const & totalDensCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::totalDensity ); + +#undef GET_FLUID_DATA + + // Enthalpy and internal energy values can be quite large and prone to round-off errors when + // performing finite difference derivatives. We will therefore scale the values using this + // reference enthalpy value to bring them to a more manageable scale. + real64 constexpr referenceEnthalpy = 5.0584e5; + auto const scaleEnthalpy = []( auto & arraySlice ) + { + LvArray::forValuesInSlice( arraySlice, []( real64 & value ){ value /= referenceEnthalpy; } ); + }; + + real64 const pressure = testData.pressure; + real64 const temperature = testData.temperature; + + // set the original fluid state to current + constitutive::constitutiveUpdatePassThru( fluid, [&] ( auto & castedFluid ) + { + typename TYPEOFREF( castedFluid ) ::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); + fluidWrapper.update( 0, 0, pressure, temperature, composition ); + } ); + + if( isThermal ) + { + scaleEnthalpy( phaseEnthalpy.value ); + scaleEnthalpy( phaseEnthalpy.derivs ); + scaleEnthalpy( phaseInternalEnergy.value ); + scaleEnthalpy( phaseInternalEnergy.derivs ); + } + + // now perturb variables and update the copied fluid's state + constitutive::constitutiveUpdatePassThru( fluidCopy, [&] ( auto & castedFluid ) + { + typename TYPEOFREF( castedFluid ) ::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); + + // to be able to use the checkDerivative utility function, we have to invert the layout + auto dPhaseFrac = invertLayout( phaseFrac.derivs.toSliceConst(), NP, NDOF ); + auto dPhaseDens = invertLayout( phaseDens.derivs.toSliceConst(), NP, NDOF ); + auto dPhaseVisc = invertLayout( phaseVisc.derivs.toSliceConst(), NP, NDOF ); + auto dPhaseEnth = invertLayout( phaseEnthalpy.derivs.toSliceConst(), NP, NDOF ); + auto dPhaseEnergy = invertLayout( phaseInternalEnergy.derivs.toSliceConst(), NP, NDOF ); + auto dTotalDens = invertLayout( totalDens.derivs.toSliceConst(), NDOF ); + auto dPhaseCompFrac = invertLayout( phaseCompFrac.derivs.toSliceConst(), NP, NC, NDOF ); + + // update pressure and check derivatives + { + real64 const dP = perturbParameter * (pressure + perturbParameter); + resetFluid( fluidCopy ); + fluidWrapper.update( 0, 0, pressure + dP, temperature, composition ); + + checkDerivative( phaseFracCopy.toSliceConst(), phaseFrac.value.toSliceConst(), dPhaseFrac[Deriv::dP].toSliceConst(), + dP, relTol, absTol, "phaseFrac", "Pres", phases ); + checkDerivative( phaseDensCopy.toSliceConst(), phaseDens.value.toSliceConst(), dPhaseDens[Deriv::dP].toSliceConst(), + dP, relTol, absTol, "phaseDens", "Pres", phases ); + checkDerivative( phaseViscCopy.toSliceConst(), phaseVisc.value.toSliceConst(), dPhaseVisc[Deriv::dP].toSliceConst(), + dP, relTol, absTol, "phaseVisc", "Pres", phases ); + checkDerivative( totalDensCopy, totalDens.value, dTotalDens[Deriv::dP], + dP, relTol, absTol, "totalDens", "Pres" ); + checkDerivative( phaseCompFracCopy.toSliceConst(), phaseCompFrac.value.toSliceConst(), dPhaseCompFrac[Deriv::dP].toSliceConst(), + dP, relTol, absTol, "phaseCompFrac", "Pres", phases, components ); + if( isThermal ) + { + scaleEnthalpy( phaseEnthCopy ); + scaleEnthalpy( phaseEnergyCopy ); + checkDerivative( phaseEnthCopy.toSliceConst(), phaseEnthalpy.value.toSliceConst(), dPhaseEnth[Deriv::dP].toSliceConst(), + dP, relTol, absTol, "phaseEnth", "Pres", phases ); + checkDerivative( phaseEnergyCopy.toSliceConst(), phaseInternalEnergy.value.toSliceConst(), dPhaseEnergy[Deriv::dP].toSliceConst(), + dP, relTol, absTol, "phaseEnergy", "Pres", phases ); + } + } + + // update temperature and check derivatives + { + real64 const dT = perturbParameter * (temperature + perturbParameter); + resetFluid( fluidCopy ); + fluidWrapper.update( 0, 0, pressure, temperature + dT, composition ); + + checkDerivative( phaseFracCopy.toSliceConst(), phaseFrac.value.toSliceConst(), dPhaseFrac[Deriv::dT].toSliceConst(), + dT, relTol, absTol, "phaseFrac", "Temp", phases ); + checkDerivative( phaseDensCopy.toSliceConst(), phaseDens.value.toSliceConst(), dPhaseDens[Deriv::dT].toSliceConst(), + dT, relTol, absTol, "phaseDens", "Temp", phases ); + checkDerivative( phaseViscCopy.toSliceConst(), phaseVisc.value.toSliceConst(), dPhaseVisc[Deriv::dT].toSliceConst(), + dT, relTol, absTol, "phaseVisc", "Temp", phases ); + checkDerivative( totalDensCopy, totalDens.value, dTotalDens[Deriv::dT], + dT, relTol, absTol, "totalDens", "Temp" ); + checkDerivative( phaseCompFracCopy.toSliceConst(), phaseCompFrac.value.toSliceConst(), dPhaseCompFrac[Deriv::dT].toSliceConst(), + dT, relTol, absTol, "phaseCompFrac", "Temp", phases, components ); + if( isThermal ) + { + scaleEnthalpy( phaseEnthCopy ); + scaleEnthalpy( phaseEnergyCopy ); + checkDerivative( phaseEnthCopy.toSliceConst(), phaseEnthalpy.value.toSliceConst(), dPhaseEnth[Deriv::dT].toSliceConst(), + dT, relTol, absTol, "phaseEnth", "Temp", phases ); + checkDerivative( phaseEnergyCopy.toSliceConst(), phaseInternalEnergy.value.toSliceConst(), dPhaseEnergy[Deriv::dT].toSliceConst(), + dT, relTol, absTol, "phaseEnergy", "Temp", phases ); + } + } + + array2d< real64, compflow::LAYOUT_COMP > compNew( 1, NC ); + for( integer jc = 0; jc < NC; ++jc ) + { + real64 const dC = LvArray::math::max( 1.0e-7, perturbParameter * ( composition[jc] + perturbParameter ) ); + for( integer ic = 0; ic < NC; ++ic ) + { + compNew[0][ic] = composition[ic]; + } + compNew[0][jc] += dC; + + // Note: in PVTPackage, derivatives are obtained with finite-difference approx **with normalization of the comp fraction** + // The component fraction is perturbed (just as above), and then all the component fractions are normalized (as below) + // But, in the native DO model and in CO2BrinePhillips, derivatives are computed analytically, which results in different + // derivatives wrt component fractions--although the derivatives wrt component densities obtained with the chain rule + // in the solver will be very similar (see discussion on PR #1325 on GitHub). + // + // Since both approaches--FD approximation of derivatives with normalization, and analytical derivatives--are correct, + // we have to support both when we check the intermediate derivatives wrt component fractions below. Therefore, if the + // PVTPackage is used, then we normalize the perturbed component fractions before taking the FD approx. If the native + // DO or CO2-brine models are used, we skip the normalization below. + if constexpr ( UsePVTPackage< FLUID_TYPE >::value ) + { + // renormalize + real64 sum = 0.0; + for( integer ic = 0; ic < NC; ++ic ) + { + sum += compNew[0][ic]; + } + for( integer ic = 0; ic < NC; ++ic ) + { + compNew[0][ic] /= sum; + } + } + + resetFluid( fluidCopy ); + fluidWrapper.update( 0, 0, pressure, temperature, compNew[0] ); + + string const var = "compFrac[" + components[jc] + "]"; + checkDerivative( phaseFracCopy.toSliceConst(), phaseFrac.value.toSliceConst(), dPhaseFrac[Deriv::dC+jc].toSliceConst(), + dC, relTol, absTol, "phaseFrac", var, phases ); + checkDerivative( phaseDensCopy.toSliceConst(), phaseDens.value.toSliceConst(), dPhaseDens[Deriv::dC+jc].toSliceConst(), + dC, relTol, absTol, "phaseDens", var, phases ); + checkDerivative( phaseViscCopy.toSliceConst(), phaseVisc.value.toSliceConst(), dPhaseVisc[Deriv::dC+jc].toSliceConst(), + dC, relTol, absTol, "phaseVisc", var, phases ); + checkDerivative( totalDensCopy, totalDens.value, dTotalDens[Deriv::dC+jc], + dC, relTol, absTol, "totalDens", var ); + checkDerivative( phaseCompFracCopy.toSliceConst(), phaseCompFrac.value.toSliceConst(), dPhaseCompFrac[Deriv::dC+jc].toSliceConst(), + dC, relTol, absTol, "phaseCompFrac", var, phases, components ); + if( isThermal ) + { + scaleEnthalpy( phaseEnthCopy ); + scaleEnthalpy( phaseEnergyCopy ); + checkDerivative( phaseEnthCopy.toSliceConst(), phaseEnthalpy.value.toSliceConst(), dPhaseEnth[Deriv::dC+jc].toSliceConst(), + dC, relTol, absTol, "phaseEnth", var, phases ); + checkDerivative( phaseEnergyCopy.toSliceConst(), phaseInternalEnergy.value.toSliceConst(), dPhaseEnergy[Deriv::dC+jc].toSliceConst(), + dC, relTol, absTol, "phaseEnergy", var, phases ); + } + } + } ); +} + +template< typename FLUID_TYPE, integer NUM_PHASE, integer NUM_COMP > +void MultiFluidTest< FLUID_TYPE, NUM_PHASE, NUM_COMP >::testValuesAgainstPreviousImplementation( typename FLUID_TYPE::KernelWrapper const & wrapper, + MultiFluidTestData< NUM_PHASE, NUM_COMP > const & testData, + real64 const relTol ) const +{ + integer constexpr numDof = numComp + 2; + + // Copy input values into an array with expected layout + array2d< real64, compflow::LAYOUT_COMP > compositionValues( 1, numComp ); + for( integer i = 0; i < numComp; ++i ) + { + compositionValues[0][i] = testData.composition[i]; + } + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const composition = compositionValues[0]; + + StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseFraction( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseFraction( 1, 1, numPhase, numDof ); + StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseDensity( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseDensity( 1, 1, numPhase, numDof ); + StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseMassDensity( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseMassDensity( 1, 1, numPhase, numDof ); + StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseViscosity( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseViscosity( 1, 1, numPhase, numDof ); + StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseEnthalpy( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseEnthalpy( 1, 1, numPhase, numDof ); + StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseInternalEnergy( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseInternalEnergy( 1, 1, numPhase, numDof ); + StackArray< real64, 4, numComp *numPhase, LAYOUT_PHASE_COMP > phaseCompFraction( 1, 1, numPhase, numComp ); + StackArray< real64, 5, numDof *numComp *numPhase, LAYOUT_PHASE_COMP_DC > dPhaseCompFraction( 1, 1, numPhase, numComp, numDof ); + StackArray< real64, 2, 1, LAYOUT_FLUID > totalDensity( 1, 1 ); + StackArray< real64, 3, numDof, LAYOUT_FLUID_DC > dTotalDensity( 1, 1, numDof ); + + wrapper.compute( testData.pressure, + testData.temperature, + composition, + { phaseFraction[0][0], dPhaseFraction[0][0] }, + { phaseDensity[0][0], dPhaseDensity[0][0] }, + { phaseMassDensity[0][0], dPhaseMassDensity[0][0] }, + { phaseViscosity[0][0], dPhaseViscosity[0][0] }, + { phaseEnthalpy[0][0], dPhaseEnthalpy[0][0] }, + { phaseInternalEnergy[0][0], dPhaseInternalEnergy[0][0] }, + { phaseCompFraction[0][0], dPhaseCompFraction[0][0] }, + { totalDensity[0][0], dTotalDensity[0][0] } ); + + checkRelativeError( totalDensity[0][0], testData.totalDensity, relTol ); + for( integer ip = 0; ip < numPhase; ++ip ) + { + checkRelativeError( phaseFraction[0][0][ip], testData.phaseFraction[ip], relTol ); + checkRelativeError( phaseDensity[0][0][ip], testData.phaseDensity[ip], relTol ); + checkRelativeError( phaseMassDensity[0][0][ip], testData.phaseMassDensity[ip], relTol ); + checkRelativeError( phaseViscosity[0][0][ip], testData.phaseViscosity[ip], relTol ); + checkRelativeError( phaseEnthalpy[0][0][ip], testData.phaseEnthalpy[ip], relTol ); + checkRelativeError( phaseInternalEnergy[0][0][ip], testData.phaseInternalEnergy[ip], relTol ); + } +} + +template< integer NUM_PHASE, integer NUM_COMP > +MultiFluidTestData< NUM_PHASE, NUM_COMP >::MultiFluidTestData( real64 const pressure_, + real64 const temperature_, + Feed< NUM_COMP > && composition_ ): + pressure( pressure_ ), + temperature( temperature_ ), + composition( composition_ ) +{} + +template< integer NUM_PHASE, integer NUM_COMP > +MultiFluidTestData< NUM_PHASE, NUM_COMP >::MultiFluidTestData( real64 const pressure_, + real64 const temperature_, + arraySlice1d< real64 const > const & composition_ ): + pressure( pressure_ ), + temperature( temperature_ ) +{ + copy< NUM_COMP >( composition_, composition ); +} + +template< integer NUM_PHASE, integer NUM_COMP > +MultiFluidTestData< NUM_PHASE, NUM_COMP >::MultiFluidTestData( real64 const pressure_, + real64 const temperature_, + Feed< NUM_COMP > && composition_, + real64 const totalDensity_, + Feed< NUM_PHASE > && phaseFraction_, + Feed< NUM_PHASE > && phaseDensity_, + Feed< NUM_PHASE > && phaseMassDensity_, + Feed< NUM_PHASE > && phaseViscosity_, + Feed< NUM_PHASE > && phaseEnthalpy_, + Feed< NUM_PHASE > && phaseInternalEnergy_ ): + pressure( pressure_ ), + temperature( temperature_ ), + composition( composition_ ), + totalDensity( totalDensity_ ), + phaseFraction( phaseFraction_ ), + phaseDensity( phaseDensity_ ), + phaseMassDensity( phaseMassDensity_ ), + phaseViscosity( phaseViscosity_ ), + phaseEnthalpy( phaseEnthalpy_ ), + phaseInternalEnergy( phaseInternalEnergy_ ) +{} + +template< integer N > +inline ::std::ostream & PrintTo( string const & name, Feed< N > const & data, ::std::ostream & os ) +{ + os << " \"" << name << "\": [ " << data[0]; + for( integer i = 1; i < N; i++ ) + { + os << ", " << data[i]; + } + os << " ]"; + return os; +} + +template< integer NUM_PHASE, integer NUM_COMP > +inline void PrintTo( MultiFluidTestData< NUM_PHASE, NUM_COMP > const & data, ::std::ostream *s ) +{ + std::ostream & os = *s; + + os << "{\n"; + os << " \"pressure\": " << data.pressure << ",\n"; + os << " \"temperature\": " << data.temperature << ",\n"; + PrintTo< NUM_COMP >( "composition", data.composition, os ) << ",\n"; + os << " \"totalDensity\": " << data.totalDensity << ",\n"; + PrintTo< NUM_PHASE >( "phaseFraction", data.phaseFraction, os ) << ",\n"; + PrintTo< NUM_PHASE >( "phaseDensity", data.phaseDensity, os ) << ",\n"; + PrintTo< NUM_PHASE >( "phaseMassDensity", data.phaseMassDensity, os ) << ",\n"; + PrintTo< NUM_PHASE >( "phaseViscosity", data.phaseViscosity, os ) << ",\n"; + PrintTo< NUM_PHASE >( "phaseEnthalpy", data.phaseEnthalpy, os ) << ",\n"; + PrintTo< NUM_PHASE >( "phaseInternalEnergy", data.phaseInternalEnergy, os ) << "\n"; + os << "}"; +} + +} // namespace testing + +} // namespace geos + +#endif //GEOS_UNITTESTS_CONSTITUTIVETESTS_MULTIFLUIDTEST_HPP_ diff --git a/src/coreComponents/unitTests/constitutiveTests/constitutiveTestHelpers.hpp b/src/coreComponents/unitTests/constitutiveTests/constitutiveTestHelpers.hpp index 2874dda757f..4b90ca2840a 100644 --- a/src/coreComponents/unitTests/constitutiveTests/constitutiveTestHelpers.hpp +++ b/src/coreComponents/unitTests/constitutiveTests/constitutiveTestHelpers.hpp @@ -21,7 +21,7 @@ #include "constitutive/capillaryPressure/capillaryPressureSelector.hpp" #include "functions/FunctionManager.hpp" #include "functions/TableFunction.hpp" -#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp" +#include "unitTests/fluidFlowTests/testFlowUtils.hpp" // TPL includes #include diff --git a/src/coreComponents/unitTests/constitutiveTests/testMultiFluid.cpp b/src/coreComponents/unitTests/constitutiveTests/testMultiFluid.cpp deleted file mode 100644 index 2297aaba3f3..00000000000 --- a/src/coreComponents/unitTests/constitutiveTests/testMultiFluid.cpp +++ /dev/null @@ -1,1173 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2018-2020 TotalEnergies - * Copyright (c) 2019- GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -// Source includes -#include "common/DataTypes.hpp" -#include "common/TimingMacros.hpp" -#include "constitutive/fluid/multifluid/MultiFluidFields.hpp" -#include "constitutive/fluid/multifluid/MultiFluidSelector.hpp" -#include "constitutive/fluid/multifluid/MultiFluidUtils.hpp" -#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp" -#include "mainInterface/initialization.hpp" -#include "functions/FunctionManager.hpp" -#include "mainInterface/GeosxState.hpp" - -// TPL includes -#include -#include - -using namespace geos; -using namespace geos::testing; -using namespace geos::constitutive; -using namespace geos::constitutive::multifluid; -using namespace geos::dataRepository; -using namespace geos::constitutive::PVTProps; - -/// Black-oil tables written into temporary files during testing - -static const char * pvtoTableContent = "# Rs[sm3/sm3]\tPbub[Pa]\tBo[m3/sm3]\tVisc(Pa.s)\n" - "\n" - " 2\t 2000000\t 1.02\t 0.000975\n" - " 5\t 5000000\t 1.03\t 0.00091\n" - " 10\t 10000000\t1.04\t 0.00083\n" - " 15\t 20000000\t1.05\t 0.000695\n" - " 90000000\t1.03\t 0.000985 -- some line comment\n" - " 30\t 30000000\t1.07\t 0.000594\n" - " 40\t 40000000\t1.08\t 0.00051\n" - " 50000000\t1.07\t 0.000549 -- another one\n" - " 90000000\t1.06\t 0.00074\n" - " 50\t 50000000.7\t1.09\t 0.000449\n" - " 90000000.7\t1.08\t 0.000605"; - -static const char * pvtwTableContent = "#\tPref[Pa]\tBw[m3/sm3]\tCp[1/Pa]\t Visc[Pa.s]\n" - "\t30600000.1\t1.03\t\t0.00000000041\t0.0003"; - -/// Dead-oil tables written into temporary files during testing - -static const char * pvdgTableContent = "# Pg(Pa) Bg(m3/sm3) Visc(Pa.s)\n" - "3000000 0.04234 0.00001344\n" - "6000000 0.02046 0.0000142\n" - "9000000 0.01328 0.00001526\n" - "12000000 0.00977 0.0000166\n" - "15000000 0.00773 0.00001818\n" - "18000000 0.006426 0.00001994\n" - "21000000 0.005541 0.00002181\n" - "24000000 0.004919 0.0000237\n" - "27000000 0.004471 0.00002559\n" - "29500000 0.004194 0.00002714\n" - "31000000 0.004031 0.00002806\n" - "33000000 0.00391 0.00002832\n" - "53000000 0.003868 0.00002935"; - -static const char * pvdoTableContent = "#P[Pa] Bo[m3/sm3] Visc(Pa.s)\n" - "10000000.0 1.23331 0.00015674\n" - "12500000.0 1.21987 0.00016570\n" - "15000000.0 1.20802 0.00017445\n" - "20000000.0 1.18791 0.00019143\n" - "25000000.0 1.17137 0.00020779\n" - "30000000.0 1.15742 0.00022361\n" - "33200000.3 1.14946 0.00023359\n" - "35000000.0 1.14543 0.00023894\n" - "40000000.0 1.13498 0.00025383\n" - "50000000.0 1.11753 0.00028237\n" - "60000000.0 1.10346 0.00030941\n" - "70000000.0 1.09180 0.00033506\n" - "80000000.0 1.08194 0.00035945\n" - "90000000.0 1.07347 0.00038266\n" - "95000000.0 1.06966 0.00039384\n" - "100000000.0 1.06610 0.00040476\n" - "110000000.0 1.05961 0.00042584\n" - "112500000.0 1.05811 0.00043096\n" - "115000000.0 1.05665 0.00043602\n" - "117500000.0 1.05523 0.00044102\n" - "120000000.0 1.05385 0.00044596\n"; - -static const char * pvdwTableContent = "# Pref[Pa] Bw[m3/sm3] Cp[1/Pa] Visc[Pa.s]\n" - " 30600000.1 1.03 0.00000000041 0.0003"; - -// CO2-brine model - -static const char * pvtLiquidPhillipsTableContent = "DensityFun PhillipsBrineDensity 1e6 1.5e7 5e4 367.15 369.15 1 0.2\n" - "ViscosityFun PhillipsBrineViscosity 0.1"; - -// the last are set relatively high (1e-4) to increase derivative value and check it properly -static const char * pvtLiquidEzrokhiTableContent = "DensityFun EzrokhiBrineDensity 2.01e-6 -6.34e-7 1e-4\n" - "ViscosityFun EzrokhiBrineViscosity 2.42e-7 0 1e-4"; - -static const char * pvtGasTableContent = "DensityFun SpanWagnerCO2Density 1e6 1.5e7 5e4 367.15 369.15 1\n" - "ViscosityFun FenghourCO2Viscosity 1e6 1.5e7 5e4 367.15 369.15 1"; - -static const char * co2FlashTableContent = "FlashModel CO2Solubility 1e6 1.5e7 5e4 367.15 369.15 1 0.15"; - -void testNumericalDerivatives( MultiFluidBase & fluid, - Group & parent, - real64 const P, - real64 const T, - arraySlice1d< real64 > const & compositionInput, - real64 const perturbParameter, - bool usePVTPackage, - real64 const relTol, - real64 const absTol = std::numeric_limits< real64 >::max() ) -{ - using Deriv = multifluid::DerivativeOffset; - - integer const NC = fluid.numFluidComponents(); - integer const NP = fluid.numFluidPhases(); - integer const NDOF = NC+2; - - // Copy input values into an array with expected layout - array2d< real64, compflow::LAYOUT_COMP > compositionValues( 1, NC ); - for( integer i = 0; i < NC; ++i ) - { - compositionValues[0][i] = compositionInput[i]; - } - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const composition = compositionValues[0]; - - auto const & components = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); - auto const & phases = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); - - // create a clone of the fluid to run updates on - std::unique_ptr< ConstitutiveBase > fluidCopyPtr = fluid.deliverClone( "fluidCopy", &parent ); - MultiFluidBase & fluidCopy = dynamicCast< MultiFluidBase & >( *fluidCopyPtr ); - - fluid.allocateConstitutiveData( fluid.getParent(), 1 ); - fluidCopy.allocateConstitutiveData( fluid.getParent(), 1 ); - - // extract data views from both fluids - #define GET_FLUID_DATA( FLUID, TRAIT ) \ - FLUID.getReference< TRAIT::type >( TRAIT::key() )[0][0] - - MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > phaseFrac { - GET_FLUID_DATA( fluid, fields::multifluid::phaseFraction ), - GET_FLUID_DATA( fluid, fields::multifluid::dPhaseFraction ) - }; - - MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > phaseDens { - GET_FLUID_DATA( fluid, fields::multifluid::phaseDensity ), - GET_FLUID_DATA( fluid, fields::multifluid::dPhaseDensity ) - }; - - MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > phaseVisc { - GET_FLUID_DATA( fluid, fields::multifluid::phaseViscosity ), - GET_FLUID_DATA( fluid, fields::multifluid::dPhaseViscosity ) - }; - - MultiFluidVarSlice< real64, 2, USD_PHASE_COMP - 2, USD_PHASE_COMP_DC - 2 > phaseCompFrac { - GET_FLUID_DATA( fluid, fields::multifluid::phaseCompFraction ), - GET_FLUID_DATA( fluid, fields::multifluid::dPhaseCompFraction ) - }; - - MultiFluidVarSlice< real64, 0, USD_FLUID - 2, USD_FLUID_DC - 2 > totalDens { - GET_FLUID_DATA( fluid, fields::multifluid::totalDensity ), - GET_FLUID_DATA( fluid, fields::multifluid::dTotalDensity ) - }; - - auto const & phaseFracCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseFraction ); - auto const & phaseDensCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseDensity ); - auto const & phaseViscCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseViscosity ); - auto const & phaseCompFracCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::phaseCompFraction ); - auto const & totalDensCopy = GET_FLUID_DATA( fluidCopy, fields::multifluid::totalDensity ); - -#undef GET_FLUID_DATA - - // set the original fluid state to current - constitutive::constitutiveUpdatePassThru( fluid, [&] ( auto & castedFluid ) - { - typename TYPEOFREF( castedFluid ) ::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); - fluidWrapper.update( 0, 0, P, T, composition ); - } ); - - // now perturb variables and update the copied fluid's state - constitutive::constitutiveUpdatePassThru( fluidCopy, [&] ( auto & castedFluid ) - { - typename TYPEOFREF( castedFluid ) ::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); - - // to be able to use the checkDerivative utility function, we have to invert the layout - auto dPhaseFrac = invertLayout( phaseFrac.derivs.toSliceConst(), NP, NDOF ); - auto dPhaseDens = invertLayout( phaseDens.derivs.toSliceConst(), NP, NDOF ); - auto dPhaseVisc = invertLayout( phaseVisc.derivs.toSliceConst(), NP, NDOF ); - auto dTotalDens = invertLayout( totalDens.derivs.toSliceConst(), NDOF ); - auto dPhaseCompFrac = invertLayout( phaseCompFrac.derivs.toSliceConst(), NP, NC, NDOF ); - - // update pressure and check derivatives - { - real64 const dP = perturbParameter * (P + perturbParameter); - fluidWrapper.update( 0, 0, P + dP, T, composition ); - - checkDerivative( phaseFracCopy.toSliceConst(), phaseFrac.value.toSliceConst(), dPhaseFrac[Deriv::dP].toSliceConst(), - dP, relTol, absTol, "phaseFrac", "Pres", phases ); - checkDerivative( phaseDensCopy.toSliceConst(), phaseDens.value.toSliceConst(), dPhaseDens[Deriv::dP].toSliceConst(), - dP, relTol, absTol, "phaseDens", "Pres", phases ); - checkDerivative( phaseViscCopy.toSliceConst(), phaseVisc.value.toSliceConst(), dPhaseVisc[Deriv::dP].toSliceConst(), - dP, relTol, absTol, "phaseVisc", "Pres", phases ); - checkDerivative( totalDensCopy, totalDens.value, dTotalDens[Deriv::dP], - dP, relTol, absTol, "totalDens", "Pres" ); - checkDerivative( phaseCompFracCopy.toSliceConst(), phaseCompFrac.value.toSliceConst(), dPhaseCompFrac[Deriv::dP].toSliceConst(), - dP, relTol, absTol, "phaseCompFrac", "Pres", phases, components ); - } - - // update temperature and check derivatives - { - real64 const dT = perturbParameter * (T + perturbParameter); - fluidWrapper.update( 0, 0, P, T + dT, composition ); - - checkDerivative( phaseFracCopy.toSliceConst(), phaseFrac.value.toSliceConst(), dPhaseFrac[Deriv::dT].toSliceConst(), - dT, relTol, absTol, "phaseFrac", "Temp", phases ); - checkDerivative( phaseDensCopy.toSliceConst(), phaseDens.value.toSliceConst(), dPhaseDens[Deriv::dT].toSliceConst(), - dT, relTol, absTol, "phaseDens", "Temp", phases ); - checkDerivative( phaseViscCopy.toSliceConst(), phaseVisc.value.toSliceConst(), dPhaseVisc[Deriv::dT].toSliceConst(), - dT, relTol, absTol, "phaseVisc", "Temp", phases ); - checkDerivative( totalDensCopy, totalDens.value, dTotalDens[Deriv::dT], - dT, relTol, absTol, "totalDens", "Temp" ); - checkDerivative( phaseCompFracCopy.toSliceConst(), phaseCompFrac.value.toSliceConst(), dPhaseCompFrac[Deriv::dT].toSliceConst(), - dT, relTol, absTol, "phaseCompFrac", "Temp", phases, components ); - } - - array2d< real64, compflow::LAYOUT_COMP > compNew( 1, NC ); - for( integer jc = 0; jc < NC; ++jc ) - { - real64 const dC = perturbParameter * ( composition[jc] + perturbParameter ); - for( integer ic = 0; ic < NC; ++ic ) - { - compNew[0][ic] = composition[ic]; - } - compNew[0][jc] += dC; - - // Note: in PVTPackage, derivatives are obtained with finite-difference approx **with normalization of the comp fraction** - // The component fraction is perturbed (just as above), and then all the component fractions are normalized (as below) - // But, in the native DO model and in CO2BrinePhillips, derivatives are computed analytically, which results in different - // derivatives wrt component fractions--although the derivatives wrt component densities obtained with the chain rule - // in the solver will be very similar (see discussion on PR #1325 on GitHub). - // - // Since both approaches--FD approximation of derivatives with normalization, and analytical derivatives--are correct, - // we have to support both when we check the intermediate derivatives wrt component fractions below. Therefore, if the - // PVTPackage is used, then we normalize the perturbed component fractions before taking the FD approx. If the native - // DO or CO2-brine models are used, we skip the normalization below. - if( usePVTPackage ) - { - // renormalize - real64 sum = 0.0; - for( integer ic = 0; ic < NC; ++ic ) - { - sum += compNew[0][ic]; - } - for( integer ic = 0; ic < NC; ++ic ) - { - compNew[0][ic] /= sum; - } - } - - fluidWrapper.update( 0, 0, P, T, compNew[0] ); - - string const var = "compFrac[" + components[jc] + "]"; - checkDerivative( phaseFracCopy.toSliceConst(), phaseFrac.value.toSliceConst(), dPhaseFrac[Deriv::dC+jc].toSliceConst(), - dC, relTol, absTol, "phaseFrac", var, phases ); - checkDerivative( phaseDensCopy.toSliceConst(), phaseDens.value.toSliceConst(), dPhaseDens[Deriv::dC+jc].toSliceConst(), - dC, relTol, absTol, "phaseDens", var, phases ); - checkDerivative( phaseViscCopy.toSliceConst(), phaseVisc.value.toSliceConst(), dPhaseVisc[Deriv::dC+jc].toSliceConst(), - dC, relTol, absTol, "phaseVisc", var, phases ); - checkDerivative( totalDensCopy, totalDens.value, dTotalDens[Deriv::dC+jc], - dC, relTol, absTol, "totalDens", var ); - checkDerivative( phaseCompFracCopy.toSliceConst(), phaseCompFrac.value.toSliceConst(), dPhaseCompFrac[Deriv::dC+jc].toSliceConst(), - dC, relTol, absTol, "phaseCompFrac", var, phases, components ); - } - } ); -} - -void testValuesAgainstPreviousImplementation( CO2BrinePhillipsFluid::KernelWrapper const & wrapper, - real64 const P, - real64 const T, - arraySlice1d< real64 const > const & compositionInput, - real64 const & savedTotalDensity, - real64 const & savedGasPhaseFrac, - real64 const & savedWaterDens, - real64 const & savedGasDens, - real64 const & savedWaterMassDens, - real64 const & savedGasMassDens, - real64 const & savedWaterVisc, - real64 const & savedGasVisc, - real64 const & savedWaterPhaseGasComp, - real64 const & savedWaterPhaseWaterComp, - real64 const relTol ) -{ - integer constexpr numPhase = 2; - integer constexpr numComp = 2; - integer constexpr numDof = numComp + 2; - - // Copy input values into an array with expected layout - array2d< real64, compflow::LAYOUT_COMP > compositionValues( 1, numComp ); - for( integer i = 0; i < numComp; ++i ) - { - compositionValues[0][i] = compositionInput[i]; - } - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const composition = compositionValues[0]; - - StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseFraction( 1, 1, numPhase ); - StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseFraction( 1, 1, numPhase, numDof ); - StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseDensity( 1, 1, numPhase ); - StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseDensity( 1, 1, numPhase, numDof ); - StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseMassDensity( 1, 1, numPhase ); - StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseMassDensity( 1, 1, numPhase, numDof ); - StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseViscosity( 1, 1, numPhase ); - StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseViscosity( 1, 1, numPhase, numDof ); - StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseEnthalpy( 1, 1, numPhase ); - StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseEnthalpy( 1, 1, numPhase, numDof ); - StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseInternalEnergy( 1, 1, numPhase ); - StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseInternalEnergy( 1, 1, numPhase, numDof ); - StackArray< real64, 4, numComp *numPhase, LAYOUT_PHASE_COMP > phaseCompFraction( 1, 1, numPhase, numComp ); - StackArray< real64, 5, numDof *numComp *numPhase, LAYOUT_PHASE_COMP_DC > dPhaseCompFraction( 1, 1, numPhase, numComp, numDof ); - StackArray< real64, 2, 1, LAYOUT_FLUID > totalDensity( 1, 1 ); - StackArray< real64, 3, numDof, LAYOUT_FLUID_DC > dTotalDensity( 1, 1, numDof ); - - wrapper.compute( P, T, composition, - { - phaseFraction[0][0], - dPhaseFraction[0][0] - }, - { - phaseDensity[0][0], - dPhaseDensity[0][0] - }, - { - phaseMassDensity[0][0], - dPhaseMassDensity[0][0] - }, - { - phaseViscosity[0][0], - dPhaseViscosity[0][0] - }, - { - phaseEnthalpy[0][0], - dPhaseEnthalpy[0][0] - }, - { - phaseInternalEnergy[0][0], - dPhaseInternalEnergy[0][0] - }, - { - phaseCompFraction[0][0], - dPhaseCompFraction[0][0] - }, - { - totalDensity[0][0], - dTotalDensity[0][0] - } ); - - checkRelativeError( totalDensity[0][0], savedTotalDensity, relTol ); - for( integer ip = 0; ip < numPhase; ++ip ) - { - real64 const savedPhaseFrac = ( ip == 0 ) ? savedGasPhaseFrac : 1 - savedGasPhaseFrac; - checkRelativeError( phaseFraction[0][0][ip], savedPhaseFrac, relTol ); - real64 const savedPhaseDens = ( ip == 0 ) ? savedGasDens : savedWaterDens; - checkRelativeError( phaseDensity[0][0][ip], savedPhaseDens, relTol ); - real64 const savedPhaseMassDens = ( ip == 0 ) ? savedGasMassDens : savedWaterMassDens; - checkRelativeError( phaseMassDensity[0][0][ip], savedPhaseMassDens, relTol ); - real64 const savedPhaseVisc = ( ip == 0 ) ? savedGasVisc : savedWaterVisc; - checkRelativeError( phaseViscosity[0][0][ip], savedPhaseVisc, relTol ); - for( integer ic = 0; ic < numComp; ++ic ) - { - real64 savedCompFrac = 0.0; - if( ip == 0 ) - { - savedCompFrac = ( ic == 0 ) ? 1 : 0; - } - else - { - savedCompFrac = ( ic == 0 ) ? savedWaterPhaseGasComp : savedWaterPhaseWaterComp; - } - checkRelativeError( phaseCompFraction[0][0][ip][ic], savedCompFrac, relTol ); - } - } -} - -MultiFluidBase & makeCompositionalFluid( string const & name, Group & parent ) -{ - CompositionalMultiphaseFluidPVTPackage & fluid = parent.registerGroup< CompositionalMultiphaseFluidPVTPackage >( name ); - - // TODO we should actually create a fake XML node with data, but this seemed easier... - - auto & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); - compNames.resize( 4 ); - compNames[0] = "N2"; compNames[1] = "C10"; compNames[2] = "C20"; compNames[3] = "H20"; - - auto & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); - molarWgt.resize( 4 ); - molarWgt[0] = 28e-3; molarWgt[1] = 134e-3; molarWgt[2] = 275e-3; molarWgt[3] = 18e-3; - - auto & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); - phaseNames.resize( 2 ); - phaseNames[0] = "oil"; phaseNames[1] = "gas"; - - auto & eqnOfState = fluid.getReference< string_array >( CompositionalMultiphaseFluidPVTPackage::viewKeyStruct::equationsOfStateString() ); - eqnOfState.resize( 2 ); - eqnOfState[0] = "PR"; eqnOfState[1] = "PR"; - - auto & critPres = fluid.getReference< array1d< real64 > >( CompositionalMultiphaseFluidPVTPackage::viewKeyStruct::componentCriticalPressureString() ); - critPres.resize( 4 ); - critPres[0] = 34e5; critPres[1] = 25.3e5; critPres[2] = 14.6e5; critPres[3] = 220.5e5; - - auto & critTemp = fluid.getReference< array1d< real64 > >( CompositionalMultiphaseFluidPVTPackage::viewKeyStruct::componentCriticalTemperatureString() ); - critTemp.resize( 4 ); - critTemp[0] = 126.2; critTemp[1] = 622.0; critTemp[2] = 782.0; critTemp[3] = 647.0; - - auto & acFactor = fluid.getReference< array1d< real64 > >( CompositionalMultiphaseFluidPVTPackage::viewKeyStruct::componentAcentricFactorString() ); - acFactor.resize( 4 ); - acFactor[0] = 0.04; acFactor[1] = 0.443; acFactor[2] = 0.816; acFactor[3] = 0.344; - - fluid.postProcessInputRecursive(); - return fluid; -} - -class CompositionalFluidTestBase : public ::testing::Test -{ -public: - CompositionalFluidTestBase(): - node(), - parent( "parent", node ) - {} - -protected: - conduit::Node node; - Group parent; - MultiFluidBase * fluid; -}; - -class CompositionalFluidTest : public CompositionalFluidTestBase -{ -public: - CompositionalFluidTest() - { - parent.resize( 1 ); - fluid = &makeCompositionalFluid( "fluid", parent ); - - parent.initialize(); - parent.initializePostInitialConditions(); - } -}; - -TEST_F( CompositionalFluidTest, numericalDerivativesMolar ) -{ - fluid->setMassFlag( false ); - - // TODO test over a range of values - real64 const P = 5e6; - real64 const T = 297.15; - array1d< real64 > comp( 4 ); - comp[0] = 0.099; comp[1] = 0.3; comp[2] = 0.6; comp[3] = 0.001; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-4; - - testNumericalDerivatives( *fluid, parent, P, T, comp, eps, true, relTol ); -} - -TEST_F( CompositionalFluidTest, numericalDerivativesMass ) -{ - fluid->setMassFlag( true ); - - // TODO test over a range of values - real64 const P = 5e6; - real64 const T = 297.15; - array1d< real64 > comp( 4 ); - comp[0] = 0.099; comp[1] = 0.3; comp[2] = 0.6; comp[3] = 0.001; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-2; - - testNumericalDerivatives( *fluid, parent, P, T, comp, eps, true, relTol ); -} - -MultiFluidBase & makeLiveOilFluid( string const & name, Group * parent ) -{ - BlackOilFluid & fluid = parent->registerGroup< BlackOilFluid >( name ); - - string_array & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); - compNames.resize( 3 ); - compNames[0] = "oil"; compNames[1] = "gas"; compNames[2] = "water"; - - array1d< real64 > & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); - molarWgt.resize( 3 ); - molarWgt[0] = 114e-3; molarWgt[1] = 16e-3; molarWgt[2] = 18e-3; - - string_array & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); - phaseNames.resize( 3 ); - phaseNames[0] = "oil"; phaseNames[1] = "gas"; phaseNames[2] = "water"; - - array1d< real64 > & surfaceDens = fluid.getReference< array1d< real64 > >( BlackOilFluidBase::viewKeyStruct::surfacePhaseMassDensitiesString() ); - surfaceDens.resize( 3 ); - surfaceDens[0] = 800.0; surfaceDens[1] = 0.9907; surfaceDens[2] = 1022.0; - - path_array & tableNames = fluid.getReference< path_array >( BlackOilFluidBase::viewKeyStruct::tableFilesString() ); - tableNames.resize( 3 ); - tableNames[0] = "pvto.txt"; tableNames[1] = "pvdg.txt"; tableNames[2] = "pvtw.txt"; - - fluid.postProcessInputRecursive(); - return fluid; -} - -MultiFluidBase & makeDeadOilFluid( string const & name, Group * parent ) -{ - DeadOilFluid & fluid = parent->registerGroup< DeadOilFluid >( name ); - - string_array & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); - compNames.resize( 3 ); - compNames[0] = "oil"; compNames[1] = "water"; compNames[2] = "gas"; - - array1d< real64 > & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); - molarWgt.resize( 3 ); - molarWgt[0] = 114e-3; molarWgt[1] = 16e-3; molarWgt[2] = 18e-3; - - string_array & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); - phaseNames.resize( 3 ); - phaseNames[0] = "oil"; phaseNames[1] = "water"; phaseNames[2] = "gas"; - - array1d< real64 > & surfaceDens = fluid.getReference< array1d< real64 > >( BlackOilFluidBase::viewKeyStruct::surfacePhaseMassDensitiesString() ); - surfaceDens.resize( 3 ); - surfaceDens[0] = 800.0; surfaceDens[1] = 1022.0; surfaceDens[2] = 0.9907; - - path_array & tableNames = fluid.getReference< path_array >( BlackOilFluidBase::viewKeyStruct::tableFilesString() ); - tableNames.resize( 3 ); - tableNames[0] = "pvdo.txt"; tableNames[1] = "pvdw.txt"; tableNames[2] = "pvdg.txt"; - - fluid.postProcessInputRecursive(); - return fluid; -} - -MultiFluidBase & makeDeadOilFluidFromTable( string const & name, Group * parent ) -{ - FunctionManager & functionManager = FunctionManager::getInstance(); - - // 1) First, define the tables (PVDO, PVDG) - - // 1D table with linear interpolation - integer const NaxisPVDO = 21; - integer const NaxisPVDG = 13; - - array1d< real64_array > coordinatesPVDO; - real64_array valuesPVDO_Bo( NaxisPVDO ); - real64_array valuesPVDO_visc( NaxisPVDO ); - coordinatesPVDO.resize( 1 ); - coordinatesPVDO[0].resize( NaxisPVDO ); - coordinatesPVDO[0][0] = 10000000.0; valuesPVDO_Bo[0] = 1.23331; valuesPVDO_visc[0] = 0.00015674; - coordinatesPVDO[0][1] = 12500000.0; valuesPVDO_Bo[1] = 1.21987; valuesPVDO_visc[1] = 0.00016570; - coordinatesPVDO[0][2] = 15000000.0; valuesPVDO_Bo[2] = 1.20802; valuesPVDO_visc[2] = 0.00017445; - coordinatesPVDO[0][3] = 20000000.0; valuesPVDO_Bo[3] = 1.18791; valuesPVDO_visc[3] = 0.00019143; - coordinatesPVDO[0][4] = 25000000.0; valuesPVDO_Bo[4] = 1.17137; valuesPVDO_visc[4] = 0.00020779; - coordinatesPVDO[0][5] = 30000000.0; valuesPVDO_Bo[5] = 1.15742; valuesPVDO_visc[5] = 0.00022361; - coordinatesPVDO[0][6] = 33200000.3; valuesPVDO_Bo[6] = 1.14946; valuesPVDO_visc[6] = 0.00023359; - coordinatesPVDO[0][7] = 35000000.0; valuesPVDO_Bo[7] = 1.14543; valuesPVDO_visc[7] = 0.00023894; - coordinatesPVDO[0][8] = 40000000.0; valuesPVDO_Bo[8] = 1.13498; valuesPVDO_visc[8] = 0.00025383; - coordinatesPVDO[0][9] = 50000000.0; valuesPVDO_Bo[9] = 1.11753; valuesPVDO_visc[9] = 0.00028237; - coordinatesPVDO[0][10] = 60000000.0; valuesPVDO_Bo[10] = 1.10346; valuesPVDO_visc[10] = 0.00030941; - coordinatesPVDO[0][11] = 70000000.0; valuesPVDO_Bo[11] = 1.09180; valuesPVDO_visc[11] = 0.00033506; - coordinatesPVDO[0][12] = 80000000.0; valuesPVDO_Bo[12] = 1.08194; valuesPVDO_visc[12] = 0.00035945; - coordinatesPVDO[0][13] = 90000000.0; valuesPVDO_Bo[13] = 1.07347; valuesPVDO_visc[13] = 0.00038266; - coordinatesPVDO[0][14] = 95000000.0; valuesPVDO_Bo[14] = 1.06966; valuesPVDO_visc[14] = 0.00039384; - coordinatesPVDO[0][15] = 100000000.0; valuesPVDO_Bo[15] = 1.06610; valuesPVDO_visc[15] = 0.00040476; - coordinatesPVDO[0][16] = 110000000.0; valuesPVDO_Bo[16] = 1.05961; valuesPVDO_visc[16] = 0.00042584; - coordinatesPVDO[0][17] = 112500000.0; valuesPVDO_Bo[17] = 1.05811; valuesPVDO_visc[17] = 0.00043096; - coordinatesPVDO[0][18] = 115000000.0; valuesPVDO_Bo[18] = 1.05665; valuesPVDO_visc[18] = 0.00043602; - coordinatesPVDO[0][19] = 117500000.0; valuesPVDO_Bo[19] = 1.05523; valuesPVDO_visc[19] = 0.00044102; - coordinatesPVDO[0][20] = 120000000.0; valuesPVDO_Bo[20] = 1.05385; valuesPVDO_visc[20] = 0.00044596; - - array1d< real64_array > coordinatesPVDG; - real64_array valuesPVDG_Bg( NaxisPVDG ); - real64_array valuesPVDG_visc( NaxisPVDG ); - coordinatesPVDG.resize( 1 ); - coordinatesPVDG[0].resize( NaxisPVDG ); - coordinatesPVDG[0][0] = 3000000; valuesPVDG_Bg[0] = 0.04234; valuesPVDG_visc[0] = 0.00001344; - coordinatesPVDG[0][1] = 6000000; valuesPVDG_Bg[1] = 0.02046; valuesPVDG_visc[1] = 0.0000142; - coordinatesPVDG[0][2] = 9000000; valuesPVDG_Bg[2] = 0.01328; valuesPVDG_visc[2] = 0.00001526; - coordinatesPVDG[0][3] = 12000000; valuesPVDG_Bg[3] = 0.00977; valuesPVDG_visc[3] = 0.0000166; - coordinatesPVDG[0][4] = 15000000; valuesPVDG_Bg[4] = 0.00773; valuesPVDG_visc[4] = 0.00001818; - coordinatesPVDG[0][5] = 18000000; valuesPVDG_Bg[5] = 0.006426; valuesPVDG_visc[5] = 0.00001994; - coordinatesPVDG[0][6] = 21000000; valuesPVDG_Bg[6] = 0.005541; valuesPVDG_visc[6] = 0.00002181; - coordinatesPVDG[0][7] = 24000000; valuesPVDG_Bg[7] = 0.004919; valuesPVDG_visc[7] = 0.0000237; - coordinatesPVDG[0][8] = 27000000; valuesPVDG_Bg[8] = 0.004471; valuesPVDG_visc[8] = 0.00002559; - coordinatesPVDG[0][9] = 29500000; valuesPVDG_Bg[9] = 0.004194; valuesPVDG_visc[9] = 0.00002714; - coordinatesPVDG[0][10] = 31000000; valuesPVDG_Bg[10] = 0.004031; valuesPVDG_visc[10] = 0.00002806; - coordinatesPVDG[0][11] = 33000000; valuesPVDG_Bg[11] = 0.00391; valuesPVDG_visc[11] = 0.00002832; - coordinatesPVDG[0][12] = 53000000; valuesPVDG_Bg[12] = 0.003868; valuesPVDG_visc[12] = 0.00002935; - - TableFunction & tablePVDO_Bo = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDO_Bo" ) ); - tablePVDO_Bo.setTableCoordinates( coordinatesPVDO, { units::Pressure } ); - tablePVDO_Bo.setTableValues( valuesPVDO_Bo, units::Dimensionless ); - tablePVDO_Bo.reInitializeFunction(); - tablePVDO_Bo.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - - TableFunction & tablePVDO_visc = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDO_visc" ) ); - tablePVDO_visc.setTableCoordinates( coordinatesPVDO, { units::Pressure } ); - tablePVDO_visc.setTableValues( valuesPVDO_visc, units::Viscosity ); - tablePVDO_visc.reInitializeFunction(); - tablePVDO_visc.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - - TableFunction & tablePVDG_Bg = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDG_Bg" ) ); - tablePVDG_Bg.setTableCoordinates( coordinatesPVDG, { units::Pressure } ); - tablePVDG_Bg.setTableValues( valuesPVDG_Bg, units::Dimensionless ); - tablePVDG_Bg.reInitializeFunction(); - tablePVDG_Bg.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - - TableFunction & tablePVDG_visc = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDG_visc" ) ); - tablePVDG_visc.setTableCoordinates( coordinatesPVDG, { units::Pressure } ); - tablePVDG_visc.setTableValues( valuesPVDG_visc, units::Viscosity ); - tablePVDG_visc.reInitializeFunction(); - tablePVDG_visc.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - - // 2) Then, define the Dead-Oil constitutive model - - DeadOilFluid & fluid = parent->registerGroup< DeadOilFluid >( name ); - - string_array & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); - compNames.resize( 3 ); - compNames[0] = "gas"; compNames[1] = "water"; compNames[2] = "oil"; - - array1d< real64 > & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); - molarWgt.resize( 3 ); - molarWgt[0] = 18e-3; molarWgt[1] = 16e-3; molarWgt[2] = 114e-3; - - string_array & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); - phaseNames.resize( 3 ); - phaseNames[0] = "gas"; phaseNames[1] = "water"; phaseNames[2] = "oil"; - - array1d< real64 > & surfaceDens = fluid.getReference< array1d< real64 > >( DeadOilFluid::viewKeyStruct::surfacePhaseMassDensitiesString() ); - surfaceDens.resize( 3 ); - surfaceDens[0] = 0.9907; surfaceDens[1] = 1022.0; surfaceDens[2] = 800.0; - - string_array & FVFTableNames = fluid.getReference< string_array >( DeadOilFluid::viewKeyStruct::formationVolumeFactorTableNamesString() ); - FVFTableNames.resize( 2 ); - FVFTableNames[0] = "PVDG_Bg"; FVFTableNames[1] = "PVDO_Bo"; - - string_array & viscosityTableNames = fluid.getReference< string_array >( DeadOilFluid::viewKeyStruct::viscosityTableNamesString() ); - viscosityTableNames.resize( 2 ); - viscosityTableNames[0] = "PVDG_visc"; viscosityTableNames[1] = "PVDO_visc"; - - real64 & waterRefPressure = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterRefPressureString() ); - waterRefPressure = 30600000.1; - real64 & waterFormationVolumeFactor = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterFormationVolumeFactorString() ); - waterFormationVolumeFactor = 1.03; - real64 & waterCompressibility = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterCompressibilityString() ); - waterCompressibility = 0.00000000041; - real64 & waterViscosity = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterViscosityString() ); - waterViscosity = 0.0003; - - fluid.postProcessInputRecursive(); - return fluid; -} - -void writeTableToFile( string const & filename, char const * str ) -{ - std::ofstream os( filename ); - ASSERT_TRUE( os.is_open() ); - os << str; - os.close(); -} - -void removeFile( string const & filename ) -{ - int const ret = std::remove( filename.c_str() ); - ASSERT_TRUE( ret == 0 ); -} - -class LiveOilFluidTest : public CompositionalFluidTestBase -{ -public: - LiveOilFluidTest() - { - writeTableToFile( "pvto.txt", pvtoTableContent ); - writeTableToFile( "pvdg.txt", pvdgTableContent ); - writeTableToFile( "pvtw.txt", pvtwTableContent ); - - parent.resize( 1 ); - fluid = &makeLiveOilFluid( "fluid", &parent ); - - parent.initialize(); - parent.initializePostInitialConditions(); - } - - ~LiveOilFluidTest() - { - removeFile( "pvto.txt" ); - removeFile( "pvdg.txt" ); - removeFile( "pvtw.txt" ); - } -}; - -TEST_F( LiveOilFluidTest, numericalDerivativesMolar ) -{ - fluid->setMassFlag( false ); - - real64 const P[3] = { 5.4e6, 1.24e7, 3.21e7 }; - real64 const T = 297.15; - array1d< real64 > comp( 3 ); - comp[0] = 0.79999; comp[1] = 0.2; comp[2] = 0.00001; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-12; - - for( integer i = 0; i < 3; ++i ) - { - testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, false, relTol ); - } -} - -TEST_F( LiveOilFluidTest, numericalDerivativesMass ) -{ - fluid->setMassFlag( true ); - - real64 const P[3] = { 5.4e6, 1.24e7, 3.21e7 }; - real64 const T = 297.15; - array1d< real64 > comp( 3 ); - comp[0] = 0.79999; comp[1] = 0.2; comp[2] = 0.00001; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-12; - - for( integer i = 0; i < 3; ++i ) - { - testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, false, relTol ); - } -} - -class DeadOilFluidTest : public CompositionalFluidTestBase -{ -public: - - DeadOilFluidTest() - { - writeTableToFile( "pvdo.txt", pvdoTableContent ); - writeTableToFile( "pvdg.txt", pvdgTableContent ); - writeTableToFile( "pvdw.txt", pvdwTableContent ); - - parent.resize( 1 ); - fluid = &makeDeadOilFluid( "fluid", &parent ); - - parent.initialize(); - parent.initializePostInitialConditions(); - } - - ~DeadOilFluidTest() - { - removeFile( "pvdo.txt" ); - removeFile( "pvdg.txt" ); - removeFile( "pvdw.txt" ); - } -}; - -TEST_F( DeadOilFluidTest, numericalDerivativesMolar ) -{ - fluid->setMassFlag( false ); - - real64 const P[3] = { 1.24e7, 3.21e7, 5.01e7 }; - real64 const T = 297.15; - array1d< real64 > comp( 3 ); - comp[0] = 0.1; comp[1] = 0.3; comp[2] = 0.6; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-4; - - for( integer i = 0; i < 3; ++i ) - { - testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, false, relTol ); - } -} - -TEST_F( DeadOilFluidTest, numericalDerivativesMass ) -{ - fluid->setMassFlag( true ); - - real64 const P[3] = { 5.4e6, 1.24e7, 3.21e7 }; - real64 const T = 297.15; - array1d< real64 > comp( 3 ); - comp[0] = 0.1; comp[1] = 0.3; comp[2] = 0.6; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-4; - real64 const absTol = 1e-14; - - for( integer i = 0; i < 3; ++i ) - { - testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, false, relTol, absTol ); - } -} - -class DeadOilFluidFromTableTest : public CompositionalFluidTestBase -{ -public: - - DeadOilFluidFromTableTest() - { - parent.resize( 1 ); - fluid = &makeDeadOilFluidFromTable( "fluid", &parent ); - - parent.initialize(); - parent.initializePostInitialConditions(); - } -}; - -TEST_F( DeadOilFluidFromTableTest, numericalDerivativesMolar ) -{ - fluid->setMassFlag( false ); - - real64 const P[3] = { 5.4e6, 1.24e7, 3.21e7 }; - real64 const T = 297.15; - array1d< real64 > comp( 3 ); - comp[0] = 0.1; comp[1] = 0.3; comp[2] = 0.6; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-4; - - for( integer i = 0; i < 3; ++i ) - { - testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, false, relTol ); - } -} - -MultiFluidBase & makeCO2BrinePhillipsFluid( string const & name, Group * parent ) -{ - CO2BrinePhillipsFluid & fluid = parent->registerGroup< CO2BrinePhillipsFluid >( name ); - - auto & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); - compNames.resize( 2 ); - compNames[0] = "co2"; compNames[1] = "water"; - - auto & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); - molarWgt.resize( 2 ); - molarWgt[0] = 44e-3; molarWgt[1] = 18e-3; - - auto & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); - phaseNames.resize( 2 ); - phaseNames[0] = "gas"; phaseNames[1] = "liquid"; - - auto & phasePVTParaFileNames = fluid.getReference< path_array >( CO2BrinePhillipsFluid::viewKeyStruct::phasePVTParaFilesString() ); - phasePVTParaFileNames.resize( 2 ); - phasePVTParaFileNames[0] = "pvtgas.txt"; phasePVTParaFileNames[1] = "pvtliquid.txt"; - - auto & flashModelParaFileName = fluid.getReference< Path >( CO2BrinePhillipsFluid::viewKeyStruct::flashModelParaFileString() ); - flashModelParaFileName = "co2flash.txt"; - - fluid.postProcessInputRecursive(); - return fluid; -} - -class CO2BrinePhillipsFluidTest : public CompositionalFluidTestBase -{ -protected: - - CO2BrinePhillipsFluidTest() - { - writeTableToFile( "pvtliquid.txt", pvtLiquidPhillipsTableContent ); - writeTableToFile( "pvtgas.txt", pvtGasTableContent ); - writeTableToFile( "co2flash.txt", co2FlashTableContent ); - - parent.resize( 1 ); - fluid = &makeCO2BrinePhillipsFluid( "fluid", &parent ); - - parent.initialize(); - parent.initializePostInitialConditions(); - } - - ~CO2BrinePhillipsFluidTest() - { - removeFile( "pvtliquid.txt" ); - removeFile( "pvtgas.txt" ); - removeFile( "co2flash.txt" ); - } - -}; - - -TEST_F( CO2BrinePhillipsFluidTest, checkAgainstPreviousImplementationMolar ) -{ - fluid->setMassFlag( false ); - - real64 const P[3] = { 5e6, 7.5e6, 1.2e7 }; - real64 const T[3] = { 367.65, 368.15, 368.75 }; - array1d< real64 > comp( 2 ); - comp[0] = 0.3; comp[1] = 0.7; - - real64 const relTol = 1e-10; - - fluid->allocateConstitutiveData( fluid->getParent(), 1 ); - - CO2BrinePhillipsFluid::KernelWrapper wrapper = - dynamicCast< CO2BrinePhillipsFluid * >( fluid )->createKernelWrapper(); - - real64 const savedTotalDens[] = - { 5881.9010529428224, 5869.6094131788523, 5855.0332090690354, 9181.3523596865525, 9157.6071613646127, 9129.5728206336826, 15757.685798517123, 15698.877814368472, - 15630.149353340639 }; - real64 const savedGasPhaseFrac[] = - { 0.29413690046142371148, 0.29415754810481165027, 0.29418169867697463449, 0.29194010802017489326, 0.29196434961986583723, 0.29199266189550621142, 0.2890641335638892695, 0.28908718137828937067, - 0.28911404840933618843 }; - real64 const savedWaterDens[] = - { 53296.719183517576, 53274.578175308554, 53247.856162690216, 53248.577831698305, 53226.801031868054, 53200.505577694363, 53232.959859840405, 53211.345942175059, - 53185.244751356993 }; - real64 const savedGasDens[] = - { 1876.2436091302606656, 1872.184636376355229, 1867.3711104617746059, 3053.1548401973859654, 3044.5748249030266379, 3034.4507978134674886, 5769.0622621289458039, 5742.8476745352018042, - 5712.2837704249559465 }; - real64 const savedWaterMassDens[] = - { 970.85108546544745423, 970.4075834766143771, 969.87385780866463847, 974.23383396044232541, 973.78856424100911227, 973.25280170872576946, 979.48333010951580491, 979.04147229150635212, - 978.50977403260912979 }; - real64 const savedGasMassDens[] = - { 82.554718801731468147, 82.376124000559627802, 82.164328860318079251, 134.33881296868497657, 133.96129229573315911, 133.51583510379256836, 253.83873953367358922, 252.68529767954885301, - 251.34048589869803436 }; - real64 const savedWaterVisc[] = - { 0.0003032144206279845924, 0.00030157070452334377216, 0.00029959815370251820189, 0.0003032144206279845924, 0.00030157070452334377216, 0.00029959815370251820189, 0.0003032144206279845924, - 0.00030157070452334377216, 0.00029959815370251820189 }; - - real64 const savedGasVisc[] = - { 1.9042384704865343673e-05, 1.9062615947696152414e-05, 1.9086923154230274463e-05, 2.0061713844617985449e-05, 2.0075955757102255573e-05, 2.0093249989250199265e-05, 2.3889596884008691474e-05, - 2.3865756080512667728e-05, 2.3839170076324036522e-05 }; - real64 const savedWaterPhaseGasComp[] = - { 0.0083062842389820552153, 0.008277274736736653718, 0.0082433415400525456018, 0.011383065290266058955, 0.011349217198060387521, 0.011309682362800700661, 0.015382352969377973903, - 0.015350431636424789056, 0.015313218057419366105 }; - real64 const savedWaterPhaseWaterComp[] = - { 0.99169371576101794652, 0.9917227252632633272, 0.99175665845994742664, 0.98861693470973388553, 0.98865078280193963156, 0.98869031763719927852, 0.98461764703062204518, 0.98464956836357520054, - 0.98468678194258063563 }; - - integer counter = 0; - for( integer i = 0; i < 3; ++i ) - { - for( integer j = 0; j < 3; ++j ) - { - testValuesAgainstPreviousImplementation( wrapper, - P[i], T[j], comp, - savedTotalDens[counter], savedGasPhaseFrac[counter], - savedWaterDens[counter], savedGasDens[counter], - savedWaterMassDens[counter], savedGasMassDens[counter], - savedWaterVisc[counter], savedGasVisc[counter], - savedWaterPhaseGasComp[counter], savedWaterPhaseWaterComp[counter], - relTol ); - counter++; - } - } -} - -TEST_F( CO2BrinePhillipsFluidTest, checkAgainstPreviousImplementationMass ) -{ - fluid->setMassFlag( true ); - - real64 const P[3] = { 5e6, 7.5e6, 1.2e7 }; - real64 const T[3] = { 367.65, 368.15, 368.75 }; - array1d< real64 > comp( 2 ); - comp[0] = 0.3; comp[1] = 0.7; - - real64 const relTol = 1e-10; - - fluid->allocateConstitutiveData( fluid->getParent(), 1 ); - - CO2BrinePhillipsFluid::KernelWrapper wrapper = - dynamicCast< CO2BrinePhillipsFluid * >( fluid )->createKernelWrapper(); - - real64 const savedTotalDens[] = - { 238.31504112633914, 237.83897306400553, 237.27445306546298, 353.95258514794097, 353.12773295711992, 352.1532278769692, 549.90502586392017, 548.2725957521294, - 546.35992000222234 }; - real64 const savedGasPhaseFrac[] = - { 0.28566797890570228, 0.28571845092287312, 0.28577748565482669, 0.28029804182709406, 0.28035729907078311, 0.2804265068556816, 0.27326788204506247, 0.27332422114692956, - 0.27338989611171055 }; - real64 const savedWaterDens[] = - { 970.85108546544745423, 970.4075834766143771, 969.87385780866463847, 974.23383396044232541, 973.78856424100911227, 973.25280170872576946, 979.48333010951580491, 979.04147229150635212, - 978.50977403260912979 }; - real64 const savedGasDens[] = - { 82.554718801731468147, 82.376124000559627802, 82.164328860318079251, 134.33881296868497657, 133.96129229573315911, 133.51583510379256836, 253.83873953367358922, 252.68529767954885301, - 251.34048589869803436 }; - real64 const savedWaterMassDens[] = - { 970.85108546544745423, 970.4075834766143771, 969.87385780866463847, 974.23383396044232541, 973.78856424100911227, 973.25280170872576946, 979.48333010951580491, 979.04147229150635212, - 978.50977403260912979 }; - real64 const savedGasMassDens[] = - { 82.554718801731468147, 82.376124000559627802, 82.164328860318079251, 134.33881296868497657, 133.96129229573315911, 133.51583510379256836, 253.83873953367358922, 252.68529767954885301, - 251.34048589869803436 }; - real64 const savedWaterVisc[] = - { 0.0003032144206279845924, 0.00030157070452334377216, 0.00029959815370251820189, 0.0003032144206279845924, 0.00030157070452334377216, 0.00029959815370251820189, 0.0003032144206279845924, - 0.00030157070452334377216, 0.00029959815370251820189 }; - real64 const savedGasVisc[] = - { 1.9042384704865343673e-05, 1.9062615947696152414e-05, 1.9086923154230274463e-05, 2.0061713844617985449e-05, 2.0075955757102255573e-05, 2.0093249989250199265e-05, 2.3889596884008691474e-05, - 2.3865756080512667728e-05, 2.3839170076324036522e-05 }; - real64 const savedWaterPhaseGasComp[] = - { 0.020063528822832473, 0.019994285300494085, 0.019913282008777046, 0.027375162661670061, 0.027295074213708581, 0.02720152052681666, 0.036784005129927563, - 0.036709327088310859, 0.036622259649145526 }; - real64 const savedWaterPhaseWaterComp[] = - { 0.97993647117716742, 0.98000571469950604, 0.98008671799122282, 0.97262483733832983, 0.97270492578629131, 0.97279847947318321, 0.96321599487007259, 0.96329067291168902, - 0.96337774035085455 }; - - integer counter = 0; - for( integer i = 0; i < 3; ++i ) - { - for( integer j = 0; j < 3; ++j ) - { - testValuesAgainstPreviousImplementation( wrapper, - P[i], T[j], comp, - savedTotalDens[counter], savedGasPhaseFrac[counter], - savedWaterDens[counter], savedGasDens[counter], - savedWaterMassDens[counter], savedGasMassDens[counter], - savedWaterVisc[counter], savedGasVisc[counter], - savedWaterPhaseGasComp[counter], savedWaterPhaseWaterComp[counter], - relTol ); - counter++; - } - } -} - -TEST_F( CO2BrinePhillipsFluidTest, numericalDerivativesMolar ) -{ - fluid->setMassFlag( false ); - - // TODO test over a range of values - real64 const P[3] = { 5e6, 7.5e6, 1.2e7 }; - real64 const T[3] = { 367.65, 368.15, 368.75 }; - array1d< real64 > comp( 2 ); - comp[0] = 0.3; comp[1] = 0.7; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-4; - - for( integer i = 0; i < 3; ++i ) - { - for( integer j = 0; j < 3; ++j ) - { - testNumericalDerivatives( *fluid, parent, P[i], T[j], comp, eps, false, relTol ); - } - } -} - -TEST_F( CO2BrinePhillipsFluidTest, numericalDerivativesMass ) -{ - fluid->setMassFlag( true ); - - // TODO test over a range of values - real64 const P[3] = { 5e6, 7.5e6, 1.2e7 }; - real64 const T[3] = { 367.65, 368.15, 368.75 }; - array1d< real64 > comp( 2 ); - comp[0] = 0.3; comp[1] = 0.7; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-8; - - for( integer i = 0; i < 3; ++i ) - { - for( integer j = 0; j < 3; ++j ) - { - testNumericalDerivatives( *fluid, parent, P[i], T[j], comp, eps, false, relTol ); - } - } -} - -MultiFluidBase & makeCO2BrineEzrokhiFluid( string const & name, Group * parent ) -{ - CO2BrineEzrokhiFluid & fluid = parent->registerGroup< CO2BrineEzrokhiFluid >( name ); - - auto & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); - compNames.resize( 2 ); - compNames[0] = "co2"; compNames[1] = "water"; - - auto & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); - molarWgt.resize( 2 ); - molarWgt[0] = 44e-3; molarWgt[1] = 18e-3; - - auto & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); - phaseNames.resize( 2 ); - phaseNames[0] = "gas"; phaseNames[1] = "liquid"; - - auto & phasePVTParaFileNames = fluid.getReference< path_array >( CO2BrineEzrokhiFluid::viewKeyStruct::phasePVTParaFilesString() ); - phasePVTParaFileNames.resize( 2 ); - phasePVTParaFileNames[0] = "pvtgas.txt"; phasePVTParaFileNames[1] = "pvtliquid.txt"; - - auto & flashModelParaFileName = fluid.getReference< Path >( CO2BrineEzrokhiFluid::viewKeyStruct::flashModelParaFileString() ); - flashModelParaFileName = "co2flash.txt"; - - fluid.postProcessInputRecursive(); - return fluid; -} - -class CO2BrineEzrokhiFluidTest : public CompositionalFluidTestBase -{ -protected: - - CO2BrineEzrokhiFluidTest() - { - writeTableToFile( "pvtliquid.txt", pvtLiquidEzrokhiTableContent ); - writeTableToFile( "pvtgas.txt", pvtGasTableContent ); - writeTableToFile( "co2flash.txt", co2FlashTableContent ); - - parent.resize( 1 ); - fluid = &makeCO2BrineEzrokhiFluid( "fluid", &parent ); - - parent.initialize(); - parent.initializePostInitialConditions(); - } - - ~CO2BrineEzrokhiFluidTest() - { - removeFile( "pvtliquid.txt" ); - removeFile( "pvtgas.txt" ); - removeFile( "co2flash.txt" ); - } - -}; - -TEST_F( CO2BrineEzrokhiFluidTest, numericalDerivativesMolar ) -{ - fluid->setMassFlag( false ); - - // TODO test over a range of values - real64 const P[3] = { 5e6, 7.5e6, 1.2e7 }; - real64 const T[3] = { 367.65, 368.15, 368.75 }; - array1d< real64 > comp( 2 ); - comp[0] = 0.3; comp[1] = 0.7; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-4; - - for( integer i = 0; i < 3; ++i ) - { - for( integer j = 0; j < 3; ++j ) - { - testNumericalDerivatives( *fluid, parent, P[i], T[j], comp, eps, false, relTol ); - } - } -} - -TEST_F( CO2BrineEzrokhiFluidTest, numericalDerivativesMass ) -{ - fluid->setMassFlag( true ); - - // TODO test over a range of values - real64 const P[3] = { 5e6, 7.5e6, 1.2e7 }; - real64 const T[3] = { 367.65, 368.15, 368.75 }; - array1d< real64 > comp( 2 ); - comp[0] = 0.3; comp[1] = 0.7; - - real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); - real64 const relTol = 1e-8; - - for( integer i = 0; i < 3; ++i ) - { - for( integer j = 0; j < 3; ++j ) - { - testNumericalDerivatives( *fluid, parent, P[i], T[j], comp, eps, false, relTol ); - } - } -} - -int main( int argc, char * * argv ) -{ - ::testing::InitGoogleTest( &argc, argv ); - - geos::GeosxState state( geos::basicSetup( argc, argv ) ); - - int const result = RUN_ALL_TESTS(); - - geos::basicCleanup(); - - return result; -} diff --git a/src/coreComponents/unitTests/constitutiveTests/testMultiFluidCO2Brine.cpp b/src/coreComponents/unitTests/constitutiveTests/testMultiFluidCO2Brine.cpp new file mode 100644 index 00000000000..d1ddabca7c7 --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/testMultiFluidCO2Brine.cpp @@ -0,0 +1,373 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file testMultiFluidCO2Brine.cpp + */ + +#include "MultiFluidTest.hpp" +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/initialization.hpp" +#include "common/GeosxConfig.hpp" + +using namespace geos; +using namespace geos::testing; +using namespace geos::constitutive; + +enum class BrineModelType : int {Phillips, Ezrokhi}; +enum class FlashType : int {DuanSun, SpycherPruess}; + +ENUM_STRINGS( BrineModelType, "Phillips", "Ezrokhi" ); +ENUM_STRINGS( FlashType, "DuanSun", "SpycherPruess" ); + +template< BrineModelType BRINE, bool THERMAL > +struct FluidType {}; + +template<> +struct FluidType< BrineModelType::Phillips, false > +{ + using type = CO2BrinePhillipsFluid; + static constexpr const char * brineContent = "DensityFun PhillipsBrineDensity 1e6 1.5e7 5e4 367.15 369.15 1 0.2\n" + "ViscosityFun PhillipsBrineViscosity 0.1"; + static constexpr const char * gasContent = "DensityFun SpanWagnerCO2Density 1e6 1.5e7 5e4 367.15 369.15 1\n" + "ViscosityFun FenghourCO2Viscosity 1e6 1.5e7 5e4 367.15 369.15 1"; +}; +template<> +struct FluidType< BrineModelType::Phillips, true > +{ + using type = CO2BrinePhillipsThermalFluid; + static constexpr const char * brineContent = "DensityFun PhillipsBrineDensity 1e6 1.5e7 5e4 367.15 369.15 1 0.2\n" + "ViscosityFun PhillipsBrineViscosity 0.1\n" + "EnthalpyFun BrineEnthalpy 1e6 7.5e7 5e5 299.15 369.15 10 0"; + static constexpr const char * gasContent = "DensityFun SpanWagnerCO2Density 1e6 1.5e7 5e4 367.15 369.15 1\n" + "ViscosityFun FenghourCO2Viscosity 1e6 1.5e7 5e4 367.15 369.15 1\n" + "EnthalpyFun CO2Enthalpy 1e6 1.5e7 5e4 367.15 369.15 1"; +}; +template<> +struct FluidType< BrineModelType::Ezrokhi, false > +{ + using type = CO2BrineEzrokhiFluid; + static constexpr const char * brineContent = "DensityFun EzrokhiBrineDensity 2.01e-6 -6.34e-7 1e-4\n" + "ViscosityFun EzrokhiBrineViscosity 2.42e-7 0 1e-4"; + static constexpr const char * gasContent = "DensityFun SpanWagnerCO2Density 1e6 1.5e7 5e4 367.15 369.15 1\n" + "ViscosityFun FenghourCO2Viscosity 1e6 1.5e7 5e4 367.15 369.15 1"; +}; +template<> +struct FluidType< BrineModelType::Ezrokhi, true > +{ + using type = CO2BrineEzrokhiThermalFluid; + static constexpr const char * brineContent = "DensityFun EzrokhiBrineDensity 2.01e-6 -6.34e-7 1e-4\n" + "ViscosityFun EzrokhiBrineViscosity 2.42e-7 0 1e-4\n" + "EnthalpyFun BrineEnthalpy 1e6 7.5e7 5e5 299.15 369.15 10 0"; + static constexpr const char * gasContent = "DensityFun SpanWagnerCO2Density 1e6 1.5e7 5e4 367.15 369.15 1\n" + "ViscosityFun FenghourCO2Viscosity 1e6 1.5e7 5e4 367.15 369.15 1\n" + "EnthalpyFun CO2Enthalpy 1e6 1.5e7 5e4 367.15 369.15 1"; +}; + +template< FlashType FLASH > +struct FlashModel {}; + +template<> +struct FlashModel< FlashType::DuanSun > +{ + static constexpr const char * flashContent = "FlashModel CO2Solubility 1e6 1.5e7 5e4 367.15 369.15 1 0.15"; +}; +template<> +struct FlashModel< FlashType::SpycherPruess > +{ + static constexpr const char * flashContent = "FlashModel CO2Solubility 1e6 1.5e7 5e4 367.15 369.15 1 0.15 1.0e-10 SpycherPruess"; +}; + +template< BrineModelType BRINE, FlashType FLASH, bool THERMAL > +class MultiFluidCO2BrineTest : public MultiFluidTest< typename FluidType< BRINE, THERMAL >::type, 2, 2 >, + public ::testing::WithParamInterface< typename MultiFluidTest< typename FluidType< BRINE, THERMAL >::type, 2, 2 >::TestData > +{ +public: + using CO2BrineFluid = typename FluidType< BRINE, THERMAL >::type; + using Base = MultiFluidTest< typename FluidType< BRINE, THERMAL >::type, 2, 2 >; + static constexpr real64 relTol = 1.0e-4; + static constexpr real64 absTol = 1.0e-4; + +public: + MultiFluidCO2BrineTest() + { + Base::writeTableToFile( pvtGasFileName, FluidType< BRINE, THERMAL >::gasContent ); + Base::writeTableToFile( pvtLiquidFileName, FluidType< BRINE, THERMAL >::brineContent ); + Base::writeTableToFile( pvtFlashFileName, FlashModel< FLASH >::flashContent ); + + auto & parent = this->m_parent; + parent.resize( 1 ); + string const fluidName = GEOS_FMT( "fluid{}{}{}", + EnumStrings< BrineModelType >::toString( BRINE ), + EnumStrings< FlashType >::toString( FLASH ), + (THERMAL ? "Thermal" : "")); + this->m_model = makeCO2BrineFluid( fluidName, &parent ); + + parent.initialize(); + parent.initializePostInitialConditions(); + } + + ~MultiFluidCO2BrineTest() override + { + Base::removeFile( pvtGasFileName ); + Base::removeFile( pvtLiquidFileName ); + Base::removeFile( pvtFlashFileName ); + } + + // Test numerical derivatives at selected data points + void testNumericalDerivatives( bool const useMass ) + { + auto * parent = &(this->getParent()); + auto & fluid = this->getFluid(); + fluid.setMassFlag( useMass ); + + real64 constexpr eps = 1.0e-5; + + // Some of the functions are simply table lookups. We need to keep the test points away from + // the table nodes because the kink in the linear interpolation might cause numerical derivative + // mismatches. Some of these values have been manually inspected and the differences, although + // not meeting the tolerance here, are small as expected. + constexpr real64 temperatures[] = { 367.65, 368.00, 368.75 }; + constexpr real64 pressures[] = { 5.001e6, 7.501e6, 1.201e7 }; + + for( real64 const pressure : pressures ) + { + for( real64 const temperature : temperatures ) + { + typename Base::TestData data ( pressure, temperature, { 0.3, 0.7 } ); + Base::testNumericalDerivatives( fluid, parent, data, eps, relTol, absTol ); + } + } + } + +private: + static CO2BrineFluid * makeCO2BrineFluid( string const & name, Group * parent ); + static constexpr const char * pvtGasFileName = "pvtgas.txt"; + static constexpr const char * pvtLiquidFileName = "pvtliquid.txt"; + static constexpr const char * pvtFlashFileName = "co2flash.txt"; +}; + +template< BrineModelType BRINE, FlashType FLASH, bool THERMAL > +typename MultiFluidCO2BrineTest< BRINE, FLASH, THERMAL >::CO2BrineFluid * +MultiFluidCO2BrineTest< BRINE, FLASH, THERMAL >::makeCO2BrineFluid( string const & name, Group * parent ) +{ + CO2BrineFluid & co2BrineFluid = parent->registerGroup< CO2BrineFluid >( name ); + + Group & fluid = co2BrineFluid; + + auto & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); + fill< 2 >( phaseNames, {"gas", "liquid"} ); + + auto & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); + fill< 2 >( compNames, {"co2", "water"} ); + + auto & molarWeight = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); + fill< 2 >( molarWeight, {44e-3, 18e-3} ); + + auto & phasePVTParaFileNames = fluid.getReference< path_array >( CO2BrineFluid::viewKeyStruct::phasePVTParaFilesString() ); + fill< 2 >( phasePVTParaFileNames, {pvtGasFileName, pvtLiquidFileName} ); + + auto & flashModelParaFileName = fluid.getReference< Path >( CO2BrineFluid::viewKeyStruct::flashModelParaFileString() ); + flashModelParaFileName = pvtFlashFileName; + + co2BrineFluid.postProcessInputRecursive(); + + return &co2BrineFluid; +} + +using MultiFluidCO2BrineBrinePhillipsTest = MultiFluidCO2BrineTest< BrineModelType::Phillips, + FlashType::DuanSun, + false >; +using MultiFluidCO2BrineBrineEzrokhiTest = MultiFluidCO2BrineTest< BrineModelType::Ezrokhi, + FlashType::DuanSun, + false >; +using MultiFluidCO2BrineBrinePhillipsThermalTest = MultiFluidCO2BrineTest< BrineModelType::Phillips, + FlashType::DuanSun, + true >; +using MultiFluidCO2BrineBrinePhillipsSpycherPruessTest = MultiFluidCO2BrineTest< BrineModelType::Phillips, + FlashType::SpycherPruess, + false >; +#if !defined(GEOS_DEVICE_COMPILE) +using MultiFluidCO2BrineBrineEzrokhiThermalTest = MultiFluidCO2BrineTest< BrineModelType::Ezrokhi, + FlashType::DuanSun, + true >; +#endif + +TEST_F( MultiFluidCO2BrineBrinePhillipsTest, numericalDerivativesMolar ) +{ + testNumericalDerivatives( false ); +} +TEST_F( MultiFluidCO2BrineBrinePhillipsTest, numericalDerivativesMass ) +{ + testNumericalDerivatives( true ); +} + +TEST_F( MultiFluidCO2BrineBrineEzrokhiTest, numericalDerivativesMolar ) +{ + testNumericalDerivatives( false ); +} +TEST_F( MultiFluidCO2BrineBrineEzrokhiTest, numericalDerivativesMass ) +{ + testNumericalDerivatives( true ); +} + +TEST_F( MultiFluidCO2BrineBrinePhillipsThermalTest, numericalDerivativesMolar ) +{ + testNumericalDerivatives( false ); +} +TEST_F( MultiFluidCO2BrineBrinePhillipsThermalTest, numericalDerivativesMass ) +{ + testNumericalDerivatives( true ); +} + +#if !defined(GEOS_DEVICE_COMPILE) +TEST_F( MultiFluidCO2BrineBrineEzrokhiThermalTest, numericalDerivativesMolar ) +{ + testNumericalDerivatives( false ); +} +TEST_F( MultiFluidCO2BrineBrineEzrokhiThermalTest, numericalDerivativesMass ) +{ + testNumericalDerivatives( true ); +} +#endif + +TEST_F( MultiFluidCO2BrineBrinePhillipsSpycherPruessTest, numericalDerivativesMolar ) +{ + testNumericalDerivatives( false ); +} +TEST_F( MultiFluidCO2BrineBrinePhillipsSpycherPruessTest, numericalDerivativesMass ) +{ + testNumericalDerivatives( true ); +} + +TEST_P( MultiFluidCO2BrineBrinePhillipsTest, checkAgainstPreviousImplementationMass ) +{ + auto & fluid = dynamicCast< CO2BrineFluid & >( this->getFluid()); + auto fluidWrapper = fluid.createKernelWrapper(); + testValuesAgainstPreviousImplementation( fluidWrapper, GetParam(), relTol ); +} + +TEST_P( MultiFluidCO2BrineBrineEzrokhiTest, checkAgainstPreviousImplementationMass ) +{ + auto & fluid = dynamicCast< CO2BrineFluid & >( this->getFluid()); + auto fluidWrapper = fluid.createKernelWrapper(); + testValuesAgainstPreviousImplementation( fluidWrapper, GetParam(), relTol ); +} + +TEST_P( MultiFluidCO2BrineBrinePhillipsThermalTest, checkAgainstPreviousImplementationMass ) +{ + auto & fluid = dynamicCast< CO2BrineFluid & >( this->getFluid()); + auto fluidWrapper = fluid.createKernelWrapper(); + testValuesAgainstPreviousImplementation( fluidWrapper, GetParam(), relTol ); +} + +TEST_P( MultiFluidCO2BrineBrinePhillipsSpycherPruessTest, checkAgainstPreviousImplementationMass ) +{ + auto & fluid = dynamicCast< CO2BrineFluid & >( this->getFluid()); + auto fluidWrapper = fluid.createKernelWrapper(); + testValuesAgainstPreviousImplementation( fluidWrapper, GetParam(), relTol ); +} + +//------------------------------------------------------------------------------- +// Data +//------------------------------------------------------------------------------- + +/* UNCRUSTIFY-OFF */ + +#define D( ... ) MultiFluidCO2BrineBrinePhillipsTest::TestData{ __VA_ARGS__ } +INSTANTIATE_TEST_SUITE_P( + MultiFluidCO2Brine, MultiFluidCO2BrineBrinePhillipsTest, + ::testing::Values( + //| pressure | temp | composition | density | phase fraction | phase density | phase mass density | phase viscosity | phase enthalpy | phase internal energy | + D( 5.00100e+06, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 5.88317e+03, {2.94136e-01, 7.05864e-01}, {1.87668e+03, 5.32967e+04}, {8.25738e+01, 9.70853e+02}, {1.90427e-05, 3.03214e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 5.00100e+06, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 5.87456e+03, {2.94150e-01, 7.05850e-01}, {1.87383e+03, 5.32812e+04}, {8.24487e+01, 9.70542e+02}, {1.90569e-05, 3.02064e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 5.00100e+06, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 5.85630e+03, {2.94181e-01, 7.05819e-01}, {1.86780e+03, 5.32478e+04}, {8.21833e+01, 9.69875e+02}, {1.90872e-05, 2.99598e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 9.18273e+03, {2.91939e-01, 7.08061e-01}, {3.05367e+03, 5.32486e+04}, {1.34361e+02, 9.74235e+02}, {2.00622e-05, 3.03214e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 9.16610e+03, {2.91956e-01, 7.08044e-01}, {3.04766e+03, 5.32333e+04}, {1.34097e+02, 9.73923e+02}, {2.00722e-05, 3.02064e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 9.13094e+03, {2.91992e-01, 7.08008e-01}, {3.03496e+03, 5.32005e+04}, {1.33538e+02, 9.73254e+02}, {2.00938e-05, 2.99598e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 1.57730e+04, {2.89059e-01, 7.10941e-01}, {5.77607e+03, 5.32330e+04}, {2.54147e+02, 9.79494e+02}, {2.39022e-05, 3.03214e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 1.57318e+04, {2.89075e-01, 7.10925e-01}, {5.75768e+03, 5.32179e+04}, {2.53338e+02, 9.79185e+02}, {2.38854e-05, 3.02064e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 1.56452e+04, {2.89109e-01, 7.10891e-01}, {5.71917e+03, 5.31853e+04}, {2.51643e+02, 9.78520e+02}, {2.38514e-05, 2.99598e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ) + ) +); +#undef D + +#define D( ... ) MultiFluidCO2BrineBrineEzrokhiTest::TestData{ __VA_ARGS__ } +INSTANTIATE_TEST_SUITE_P( + MultiFluidCO2Brine, MultiFluidCO2BrineBrineEzrokhiTest, + ::testing::Values( + //| pressure | temp | composition | density | phase fraction | phase density | phase mass density | phase viscosity | phase enthalpy | phase internal energy | + D( 5.00100e+06, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 5.89876e+03, {2.94136e-01, 7.05864e-01}, {1.87668e+03, 5.51674e+04}, {8.25738e+01, 1.00493e+03}, {1.90427e-05, 3.12148e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 5.00100e+06, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 5.89023e+03, {2.94150e-01, 7.05850e-01}, {1.87383e+03, 5.51663e+04}, {8.24487e+01, 1.00488e+03}, {1.90569e-05, 3.11023e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 5.00100e+06, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 5.87213e+03, {2.94181e-01, 7.05819e-01}, {1.86780e+03, 5.51642e+04}, {8.21833e+01, 1.00478e+03}, {1.90872e-05, 3.08611e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 9.23469e+03, {2.91939e-01, 7.08061e-01}, {3.05367e+03, 5.58209e+04}, {1.34361e+02, 1.02130e+03}, {2.00622e-05, 3.16876e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 9.21829e+03, {2.91956e-01, 7.08044e-01}, {3.04766e+03, 5.58254e+04}, {1.34097e+02, 1.02135e+03}, {2.00722e-05, 3.15764e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 9.18360e+03, {2.91992e-01, 7.08008e-01}, {3.03496e+03, 5.58353e+04}, {1.33538e+02, 1.02146e+03}, {2.00938e-05, 3.13381e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 1.59791e+04, {2.89059e-01, 7.10941e-01}, {5.77607e+03, 5.67057e+04}, {2.54147e+02, 1.04339e+03}, {2.39022e-05, 3.23075e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 1.59385e+04, {2.89075e-01, 7.10925e-01}, {5.75768e+03, 5.67188e+04}, {2.53338e+02, 1.04360e+03}, {2.38854e-05, 3.21991e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 1.58533e+04, {2.89109e-01, 7.10891e-01}, {5.71917e+03, 5.67472e+04}, {2.51643e+02, 1.04405e+03}, {2.38514e-05, 3.19665e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ) + ) +); +#undef D + +#define D( ... ) MultiFluidCO2BrineBrinePhillipsThermalTest::TestData{ __VA_ARGS__ } +INSTANTIATE_TEST_SUITE_P( + MultiFluidCO2Brine, MultiFluidCO2BrineBrinePhillipsThermalTest, + ::testing::Values( + //| pressure | temp | composition | density | phase fraction | phase density | phase mass density | phase viscosity | phase enthalpy | phase internal energy | + D( 5.00100e+06, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 5.88317e+03, {2.94136e-01, 7.05864e-01}, {1.87668e+03, 5.32967e+04}, {8.25738e+01, 9.70853e+02}, {1.90427e-05, 3.03214e-04}, {1.21447e+07, 2.18796e+07}, {1.20841e+07, 2.18744e+07} ), + D( 5.00100e+06, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 5.87456e+03, {2.94150e-01, 7.05850e-01}, {1.87383e+03, 5.32812e+04}, {8.24487e+01, 9.70542e+02}, {1.90569e-05, 3.02064e-04}, {1.21537e+07, 2.19628e+07}, {1.20931e+07, 2.19576e+07} ), + D( 5.00100e+06, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 5.85630e+03, {2.94181e-01, 7.05819e-01}, {1.86780e+03, 5.32478e+04}, {8.21833e+01, 9.69875e+02}, {1.90872e-05, 2.99598e-04}, {1.21731e+07, 2.21410e+07}, {1.21123e+07, 2.21359e+07} ), + D( 7.50100e+06, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 9.18273e+03, {2.91939e-01, 7.08061e-01}, {3.05367e+03, 5.32486e+04}, {1.34361e+02, 9.74235e+02}, {2.00622e-05, 3.03214e-04}, {1.17163e+07, 2.18445e+07}, {1.16605e+07, 2.18368e+07} ), + D( 7.50100e+06, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 9.16610e+03, {2.91956e-01, 7.08044e-01}, {3.04766e+03, 5.32333e+04}, {1.34097e+02, 9.73923e+02}, {2.00722e-05, 3.02064e-04}, {1.17269e+07, 2.19275e+07}, {1.16709e+07, 2.19198e+07} ), + D( 7.50100e+06, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 9.13094e+03, {2.91992e-01, 7.08008e-01}, {3.03496e+03, 5.32005e+04}, {1.33538e+02, 9.73254e+02}, {2.00938e-05, 2.99598e-04}, {1.17494e+07, 2.21054e+07}, {1.16932e+07, 2.20977e+07} ), + D( 1.20100e+07, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 1.57730e+04, {2.89059e-01, 7.10941e-01}, {5.77607e+03, 5.32330e+04}, {2.54147e+02, 9.79494e+02}, {2.39022e-05, 3.03214e-04}, {1.08306e+07, 2.17898e+07}, {1.07833e+07, 2.17775e+07} ), + D( 1.20100e+07, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 1.57318e+04, {2.89075e-01, 7.10925e-01}, {5.75768e+03, 5.32179e+04}, {2.53338e+02, 9.79185e+02}, {2.38854e-05, 3.02064e-04}, {1.08453e+07, 2.18726e+07}, {1.07979e+07, 2.18603e+07} ), + D( 1.20100e+07, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 1.56452e+04, {2.89109e-01, 7.10891e-01}, {5.71917e+03, 5.31853e+04}, {2.51643e+02, 9.78520e+02}, {2.38514e-05, 2.99598e-04}, {1.08767e+07, 2.20500e+07}, {1.08290e+07, 2.20378e+07} ) + ) +); +#undef D + +#define D( ... ) MultiFluidCO2BrineBrinePhillipsSpycherPruessTest::TestData{ __VA_ARGS__ } +INSTANTIATE_TEST_SUITE_P( + MultiFluidCO2Brine, MultiFluidCO2BrineBrinePhillipsSpycherPruessTest, + ::testing::Values( + //| pressure | temp | composition | density | phase fraction | phase density | phase mass density | phase viscosity | phase enthalpy | phase internal energy | + D( 5.00100e+06, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 5.77217e+03, {3.00488e-01, 6.99512e-01}, {1.87668e+03, 5.32842e+04}, {8.25738e+01, 9.70985e+02}, {1.90427e-05, 3.03214e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 5.00100e+06, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 5.76239e+03, {3.00580e-01, 6.99420e-01}, {1.87383e+03, 5.32686e+04}, {8.24487e+01, 9.70677e+02}, {1.90569e-05, 3.02064e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 5.00100e+06, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 5.74154e+03, {3.00781e-01, 6.99219e-01}, {1.86780e+03, 5.32348e+04}, {8.21833e+01, 9.70013e+02}, {1.90872e-05, 2.99598e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 9.05928e+03, {2.96731e-01, 7.03269e-01}, {3.05367e+03, 5.32295e+04}, {1.34361e+02, 9.74440e+02}, {2.00622e-05, 3.03214e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 9.04143e+03, {2.96804e-01, 7.03196e-01}, {3.04766e+03, 5.32141e+04}, {1.34097e+02, 9.74130e+02}, {2.00722e-05, 3.02064e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 7.50100e+06, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 9.00358e+03, {2.96962e-01, 7.03038e-01}, {3.03496e+03, 5.31808e+04}, {1.33538e+02, 9.73465e+02}, {2.00938e-05, 2.99598e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.67650e+02, {3.00000e-01, 7.00000e-01}, 1.56195e+04, {2.93050e-01, 7.06950e-01}, {5.77607e+03, 5.32048e+04}, {2.54147e+02, 9.79798e+02}, {2.39022e-05, 3.03214e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.68000e+02, {3.00000e-01, 7.00000e-01}, 1.55771e+04, {2.93105e-01, 7.06895e-01}, {5.75768e+03, 5.31894e+04}, {2.53338e+02, 9.79492e+02}, {2.38854e-05, 3.02064e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ), + D( 1.20100e+07, 3.68750e+02, {3.00000e-01, 7.00000e-01}, 1.54878e+04, {2.93225e-01, 7.06775e-01}, {5.71917e+03, 5.31561e+04}, {2.51643e+02, 9.78834e+02}, {2.38514e-05, 2.99598e-04}, {0.00000e+00, 0.00000e+00}, {0.00000e+00, 0.00000e+00} ) + ) +); +#undef D + +/* UNCRUSTIFY-ON */ + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + + geos::GeosxState state( geos::basicSetup( argc, argv ) ); + + int const result = RUN_ALL_TESTS(); + + geos::basicCleanup(); + + return result; +} diff --git a/src/coreComponents/unitTests/constitutiveTests/testMultiFluidCompositionalMultiphasePVTPackage.cpp b/src/coreComponents/unitTests/constitutiveTests/testMultiFluidCompositionalMultiphasePVTPackage.cpp new file mode 100644 index 00000000000..45a7651f812 --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/testMultiFluidCompositionalMultiphasePVTPackage.cpp @@ -0,0 +1,198 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file testMultiFluidCompositionalMultiphasePVTPackage.cpp + */ + +#include "MultiFluidTest.hpp" +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/initialization.hpp" + +using namespace geos; +using namespace geos::testing; +using namespace geos::constitutive; + +namespace geos +{ +namespace testing +{ + +template< > +struct UsePVTPackage< CompositionalMultiphaseFluidPVTPackage > +{ + static constexpr bool value = true; +}; + +} // namespace testing + +} // namespace geos + +enum class EOS_TYPE : int { PR, SRK }; +ENUM_STRINGS( EOS_TYPE, "PR", "SRK" ); + +template< integer NUM_COMP > +struct Fluid +{}; + +template< EOS_TYPE EOS, integer NUM_COMP > +class MultiFluidCompositionalMultiphasePVTPackageTest : public MultiFluidTest< CompositionalMultiphaseFluidPVTPackage, 2, NUM_COMP > +{ +public: + using Base = MultiFluidTest< CompositionalMultiphaseFluidPVTPackage, 2, NUM_COMP >; + static constexpr real64 relTol = 1.0e-4; +public: + MultiFluidCompositionalMultiphasePVTPackageTest() + { + auto & parent = this->m_parent; + parent.resize( 1 ); + + string fluidName = GEOS_FMT( "fluid{}{}", EnumStrings< EOS_TYPE >::toString( EOS ), NUM_COMP ); + this->m_model = makeFluid( fluidName, &parent ); + + parent.initialize(); + parent.initializePostInitialConditions(); + } + + ~MultiFluidCompositionalMultiphasePVTPackageTest() override = default; + + void testNumericatDerivatives( const bool useMass ) + { + auto & fluid = this->getFluid(); + fluid.setMassFlag( useMass ); + + auto * parent = &(this->getParent()); + + array2d< real64 > samples; + Fluid< Base::numComp >::getSamples( samples ); + integer const sampleCount = samples.size( 0 ); + + real64 constexpr eps = 1.0e-7; + + constexpr real64 pressures[] = { 1.0e5, 50.0e5, 100.0e5, 600.0e5 }; + constexpr real64 temperatures[] = { 15.5, 24.0, 40.0, 80.0 }; + + for( integer sampleIndex = 0; sampleIndex < sampleCount; ++sampleIndex ) + { + for( real64 const pressure : pressures ) + { + for( real64 const temperature : temperatures ) + { + typename Base::TestData data ( pressure, units::convertCToK( temperature ), samples[sampleIndex].toSliceConst() ); + Base::testNumericalDerivatives( fluid, parent, data, eps, relTol ); + } + } + } + } + +private: + static CompositionalMultiphaseFluidPVTPackage * makeFluid( string const & name, Group * parent ); +}; + +template< integer NUM_COMP > +static void fillBinaryCoeffs( array2d< real64 > & binaryCoeff, std::array< real64 const, NUM_COMP *(NUM_COMP-1)/2 > const data ) +{ + auto bic = data.begin(); + binaryCoeff.resize( NUM_COMP, NUM_COMP ); + for( integer i = 0; i < NUM_COMP; ++i ) + { + binaryCoeff( i, i ) = 0.0; + for( integer j = i+1; j < NUM_COMP; ++j ) + { + binaryCoeff( i, j ) = *bic++; + binaryCoeff( j, i ) = binaryCoeff( i, j ); + } + } +} + +template< > +struct Fluid< 4 > +{ + static void fillProperties( Group & fluid ) + { + string_array & componentNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); + fill< 4 >( componentNames, {"N2", "C10", "C20", "H20"} ); + + array1d< real64 > & molarWeight = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); + fill< 4 >( molarWeight, {28e-3, 134e-3, 275e-3, 18e-3} ); + + array1d< real64 > & criticalPressure = fluid.getReference< array1d< real64 > >( CompositionalMultiphaseFluidPVTPackage::viewKeyStruct::componentCriticalPressureString() ); + fill< 4 >( criticalPressure, {34e5, 25.3e5, 14.6e5, 220.5e5} ); + array1d< real64 > & criticalTemperature = fluid.getReference< array1d< real64 > >( CompositionalMultiphaseFluidPVTPackage::viewKeyStruct::componentCriticalTemperatureString() ); + fill< 4 >( criticalTemperature, {126.2, 622.0, 782.0, 647.0} ); + array1d< real64 > & acentricFactor = fluid.getReference< array1d< real64 > >( CompositionalMultiphaseFluidPVTPackage::viewKeyStruct::componentAcentricFactorString() ); + fill< 4 >( acentricFactor, {0.04, 0.443, 0.816, 0.344} ); + } + + static void getSamples( array2d< real64 > & samples ) + { + samples.resize( 2, 4 ); + fill< 4 >( samples[0], {0.099, 0.300, 0.600, 0.001} ); + fill< 4 >( samples[1], {0.000, 0.000, 0.000, 1.000} ); + } +}; + +template< EOS_TYPE EOS, integer NUM_COMP > +CompositionalMultiphaseFluidPVTPackage * +MultiFluidCompositionalMultiphasePVTPackageTest< EOS, NUM_COMP >:: +makeFluid( string const & name, Group * parent ) +{ + CompositionalMultiphaseFluidPVTPackage & fluid = parent->registerGroup< CompositionalMultiphaseFluidPVTPackage >( name ); + + string_array & phaseNames = fluid.getReference< string_array >( string( MultiFluidBase::viewKeyStruct::phaseNamesString()) ); + fill< 2 >( phaseNames, {"oil", "gas"} ); + + string const eosName = EnumStrings< EOS_TYPE >::toString( EOS ); + string_array & equationOfState = fluid.getReference< string_array >( CompositionalMultiphaseFluidPVTPackage::viewKeyStruct::equationsOfStateString() ); + fill< 2 >( equationOfState, {eosName, eosName} ); + + Fluid< NUM_COMP >::fillProperties( fluid ); + + fluid.postProcessInputRecursive(); + return &fluid; +} + +using PengRobinson4Test = MultiFluidCompositionalMultiphasePVTPackageTest< EOS_TYPE::PR, 4 >; +using SoaveRedlichKwong4Test = MultiFluidCompositionalMultiphasePVTPackageTest< EOS_TYPE::SRK, 4 >; + +TEST_F( PengRobinson4Test, numericalDerivativesMolar ) +{ + testNumericatDerivatives( false ); +} +TEST_F( PengRobinson4Test, numericalDerivativesMass ) +{ + testNumericatDerivatives( true ); +} + +TEST_F( SoaveRedlichKwong4Test, numericalDerivativesMolar ) +{ + testNumericatDerivatives( false ); +} +TEST_F( SoaveRedlichKwong4Test, numericalDerivativesMass ) +{ + testNumericatDerivatives( true ); +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + + geos::GeosxState state( geos::basicSetup( argc, argv ) ); + + int const result = RUN_ALL_TESTS(); + + geos::basicCleanup(); + + return result; +} diff --git a/src/coreComponents/unitTests/constitutiveTests/testMultiFluidDeadOil.cpp b/src/coreComponents/unitTests/constitutiveTests/testMultiFluidDeadOil.cpp new file mode 100644 index 00000000000..d8356ae9a81 --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/testMultiFluidDeadOil.cpp @@ -0,0 +1,275 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file testMultiFluidDeadOil.cpp + */ + +#include "MultiFluidTest.hpp" +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/initialization.hpp" + +using namespace geos; +using namespace geos::testing; +using namespace geos::constitutive; + +static constexpr char const * pvdgTableContent = "# Pg(Pa) Bg(m3/sm3) Visc(Pa.s)\n" + "3000000 0.04234 0.00001344\n" + "6000000 0.02046 0.0000142\n" + "9000000 0.01328 0.00001526\n" + "12000000 0.00977 0.0000166\n" + "15000000 0.00773 0.00001818\n" + "18000000 0.006426 0.00001994\n" + "21000000 0.005541 0.00002181\n" + "24000000 0.004919 0.0000237\n" + "27000000 0.004471 0.00002559 -- this is a comment\n" + "29500000 0.004194 0.00002714\n" + "31000000 0.004031 0.00002806\n" + "33000000 0.00391 0.00002832\n" + "53000000 0.003868 0.00002935"; + +static constexpr char const * pvdoTableContent = "#P[Pa] Bo[m3/sm3] Visc(Pa.s)\n" + "10000000.0 1.23331 0.00015674\n" + "12500000.0 1.21987 0.00016570\n" + "15000000.0 1.20802 0.00017445\n" + "20000000.0 1.18791 0.00019143\n" + "25000000.0 1.17137 0.00020779\n" + "30000000.0 1.15742 0.00022361\n" + "33200000.3 1.14946 0.00023359\n" + "35000000.0 1.14543 0.00023894\n" + "40000000.0 1.13498 0.00025383 -- this is a comment\n" + "50000000.0 1.11753 0.00028237\n" + "60000000.0 1.10346 0.00030941\n" + "70000000.0 1.09180 0.00033506\n" + "80000000.0 1.08194 0.00035945\n" + "90000000.0 1.07347 0.00038266\n" + "95000000.0 1.06966 0.00039384\n" + "100000000.0 1.06610 0.00040476\n" + "110000000.0 1.05961 0.00042584\n" + "112500000.0 1.05811 0.00043096\n" + "115000000.0 1.05665 0.00043602\n" + "117500000.0 1.05523 0.00044102\n" + "120000000.0 1.05385 0.00044596\n"; + +static constexpr char const * pvdwTableContent = "# Pref[Pa] Bw[m3/sm3] Cp[1/Pa] Visc[Pa.s]\n" + " 30600000.1 1.03 0.00000000041 0.0003"; + +template< bool FROM_TABLE > +class MultiFluidDeadOilTest : public MultiFluidTest< DeadOilFluid, 3, 3 > +{ +public: + static constexpr real64 relTol = 1.0e-4; + static constexpr real64 absTol = 1.0e-4; +public: + MultiFluidDeadOilTest() + { + if constexpr (!FROM_TABLE) + { + writeTableToFile( "pvdo.txt", pvdoTableContent ); + writeTableToFile( "pvdg.txt", pvdgTableContent ); + writeTableToFile( "pvdw.txt", pvdwTableContent ); + } + + m_parent.resize( 1 ); + string const fluidName = GEOS_FMT( "fluid{}", (FROM_TABLE ? "Tables" : "Files")); + m_model = makeDeadOilFluid( fluidName, &m_parent ); + + m_parent.initialize(); + m_parent.initializePostInitialConditions(); + } + + ~MultiFluidDeadOilTest() override + { + if constexpr (!FROM_TABLE) + { + removeFile( "pvdo.txt" ); + removeFile( "pvdg.txt" ); + removeFile( "pvdw.txt" ); + } + } + +private: + static DeadOilFluid * makeDeadOilFluid( string const & name, Group * parent ); + static void fillPhysicalProperties( DeadOilFluid & fluid ); +}; + +template< bool FROM_TABLE > +void MultiFluidDeadOilTest< FROM_TABLE >::fillPhysicalProperties( DeadOilFluid & fluid ) +{ + string_array & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); + fill< 3 >( phaseNames, {"oil", "water", "gas"} ); + + string_array & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); + fill< 3 >( compNames, {"oil", "water", "gas"} ); + + array1d< real64 > & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); + fill< 3 >( molarWgt, {114e-3, 18e-3, 16e-3} ); + + array1d< real64 > & surfaceDens = fluid.getReference< array1d< real64 > >( BlackOilFluidBase::viewKeyStruct::surfacePhaseMassDensitiesString() ); + fill< 3 >( surfaceDens, {800.0, 1022.0, 0.9907} ); +} + +template<> +DeadOilFluid * MultiFluidDeadOilTest< false >::makeDeadOilFluid( string const & name, Group * parent ) +{ + DeadOilFluid & fluid = parent->registerGroup< DeadOilFluid >( name ); + + fillPhysicalProperties( fluid ); + + path_array & tableNames = fluid.getReference< path_array >( BlackOilFluidBase::viewKeyStruct::tableFilesString() ); + fill< 3 >( tableNames, {"pvdo.txt", "pvdw.txt", "pvdg.txt"} ); + + fluid.postProcessInputRecursive(); + return &fluid; +} + +template<> +DeadOilFluid * MultiFluidDeadOilTest< true >::makeDeadOilFluid( string const & name, Group * parent ) +{ + // 1) First, define the tables (PVDO, PVDG) + + // 1D table with linear interpolation + integer constexpr NaxisPVDO = 21; + integer constexpr NaxisPVDG = 13; + + array1d< real64_array > coordinatesPVDO( 1 ); + real64_array valuesPVDO_Bo; + real64_array valuesPVDO_visc; + fill< NaxisPVDO >( coordinatesPVDO[0], { + 1.000e+07, 1.250e+07, 1.500e+07, 2.000e+07, 2.500e+07, 3.000e+07, 3.320e+07, 3.500e+07, + 4.000e+07, 5.000e+07, 6.000e+07, 7.000e+07, 8.000e+07, 9.000e+07, 9.500e+07, 1.000e+08, + 1.100e+08, 1.125e+08, 1.150e+08, 1.175e+08, 1.200e+08 } ); + fill< NaxisPVDO >( valuesPVDO_Bo, { + 1.23331, 1.21987, 1.20802, 1.18791, 1.17137, 1.15742, 1.14946, 1.14543, + 1.13498, 1.11753, 1.10346, 1.09180, 1.08194, 1.07347, 1.06966, 1.06610, + 1.05961, 1.05811, 1.05665, 1.05523, 1.05385 } ); + fill< NaxisPVDO >( valuesPVDO_visc, { + 1.56740e-04, 1.65700e-04, 1.74450e-04, 1.91430e-04, 2.07790e-04, 2.23610e-04, 2.33590e-04, 2.38940e-04, + 2.53830e-04, 2.82370e-04, 3.09410e-04, 3.35060e-04, 3.59450e-04, 3.82660e-04, 3.93840e-04, 4.04760e-04, + 4.25840e-04, 4.30960e-04, 4.36020e-04, 4.41020e-04, 4.45960e-04 } ); + + array1d< real64_array > coordinatesPVDG( 1 ); + real64_array valuesPVDG_Bg; + real64_array valuesPVDG_visc; + fill< NaxisPVDG >( coordinatesPVDG[0], { + 3.000e+06, 6.000e+06, 9.000e+06, 1.200e+07, 1.500e+07, 1.800e+07, 2.100e+07, 2.400e+07, + 2.700e+07, 2.950e+07, 3.100e+07, 3.300e+07, 5.300e+07 } ); + fill< NaxisPVDG >( valuesPVDG_Bg, { + 4.23400e-02, 2.04600e-02, 1.32800e-02, 9.77000e-03, 7.73000e-03, 6.42600e-03, 5.54100e-03, 4.91900e-03, + 4.47100e-03, 4.19400e-03, 4.03100e-03, 3.91000e-03, 3.86800e-03 } ); + fill< NaxisPVDG >( valuesPVDG_visc, { + 1.34400e-05, 1.42000e-05, 1.52600e-05, 1.66000e-05, 1.81800e-05, 1.99400e-05, 2.18100e-05, 2.37000e-05, + 2.55900e-05, 2.71400e-05, 2.80600e-05, 2.83200e-05, 2.93500e-05 } ); + + initializeTable( "PVDO_Bo", coordinatesPVDO, valuesPVDO_Bo ); + initializeTable( "PVDO_visc", coordinatesPVDO, valuesPVDO_visc ); + initializeTable( "PVDG_Bg", coordinatesPVDG, valuesPVDG_Bg ); + initializeTable( "PVDG_visc", coordinatesPVDG, valuesPVDG_visc ); + + // 2) Then, define the Dead-Oil constitutive model + + DeadOilFluid & fluid = parent->registerGroup< DeadOilFluid >( name ); + + fillPhysicalProperties( fluid ); + + string_array & FVFTableNames = fluid.getReference< string_array >( DeadOilFluid::viewKeyStruct::formationVolumeFactorTableNamesString() ); + fill< 2 >( FVFTableNames, {"PVDG_Bg", "PVDO_Bo"} ); + + string_array & viscosityTableNames = fluid.getReference< string_array >( DeadOilFluid::viewKeyStruct::viscosityTableNamesString() ); + fill< 2 >( viscosityTableNames, {"PVDG_visc", "PVDO_visc"} ); + + real64 & waterRefPressure = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterRefPressureString() ); + waterRefPressure = 30600000.1; + real64 & waterFormationVolumeFactor = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterFormationVolumeFactorString() ); + waterFormationVolumeFactor = 1.03; + real64 & waterCompressibility = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterCompressibilityString() ); + waterCompressibility = 0.00000000041; + real64 & waterViscosity = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterViscosityString() ); + waterViscosity = 0.0003; + + fluid.postProcessInputRecursive(); + return &fluid; +} + +using MultiFluidDeadOilTestFromFiles = MultiFluidDeadOilTest< false >; +using MultiFluidDeadOilTestFromTables = MultiFluidDeadOilTest< true >; + +TEST_F( MultiFluidDeadOilTestFromFiles, numericalDerivativesMolar ) +{ + auto & fluid = getFluid(); + fluid.setMassFlag( false ); + + real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); + + for( real64 const pressure : { 1.24e7, 3.21e7, 5.01e7 } ) + { + TestData data ( pressure, 297.15, { 0.1, 0.3, 0.6 } ); + testNumericalDerivatives( fluid, &getParent(), data, eps, relTol, absTol ); + } +} + +TEST_F( MultiFluidDeadOilTestFromFiles, numericalDerivativesMass ) +{ + auto & fluid = getFluid(); + fluid.setMassFlag( true ); + + real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); + + for( real64 const pressure : { 1.24e7, 3.21e7, 5.01e7 } ) + { + TestData data ( pressure, 297.15, { 0.1, 0.3, 0.6 } ); + testNumericalDerivatives( fluid, &getParent(), data, eps, relTol, absTol ); + } +} + +TEST_F( MultiFluidDeadOilTestFromTables, numericalDerivativesMolar ) +{ + auto & fluid = getFluid(); + fluid.setMassFlag( false ); + + real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); + + for( real64 const pressure : { 1.24e7, 3.21e7, 5.01e7 } ) + { + TestData data ( pressure, 297.15, { 0.1, 0.3, 0.6 } ); + testNumericalDerivatives( fluid, &getParent(), data, eps, relTol, absTol ); + } +} + +TEST_F( MultiFluidDeadOilTestFromTables, numericalDerivativesMass ) +{ + auto & fluid = getFluid(); + fluid.setMassFlag( true ); + + real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); + + for( real64 const pressure : { 1.24e7, 3.21e7, 5.01e7 } ) + { + TestData data ( pressure, 297.15, { 0.1, 0.3, 0.6 } ); + testNumericalDerivatives( fluid, &getParent(), data, eps, relTol, absTol ); + } +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + + geos::GeosxState state( geos::basicSetup( argc, argv ) ); + + int const result = RUN_ALL_TESTS(); + + geos::basicCleanup(); + + return result; +} diff --git a/src/coreComponents/unitTests/constitutiveTests/testMultiFluidLiveOil.cpp b/src/coreComponents/unitTests/constitutiveTests/testMultiFluidLiveOil.cpp new file mode 100644 index 00000000000..4ea16870680 --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/testMultiFluidLiveOil.cpp @@ -0,0 +1,168 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file testMultiFluidLiveOil.cpp + */ + +#include "MultiFluidTest.hpp" +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/initialization.hpp" + +using namespace geos; +using namespace geos::testing; +using namespace geos::constitutive; + +static constexpr char const * pvdgTableContent = "# Pg(Pa) Bg(m3/sm3) Visc(Pa.s)\n" + "3000000 0.04234 0.00001344\n" + "6000000 0.02046 0.0000142\n" + "9000000 0.01328 0.00001526\n" + "12000000 0.00977 0.0000166\n" + "15000000 0.00773 0.00001818\n" + "18000000 0.006426 0.00001994\n" + "21000000 0.005541 0.00002181\n" + "24000000 0.004919 0.0000237\n" + "27000000 0.004471 0.00002559 -- this is a comment\n" + "29500000 0.004194 0.00002714\n" + "31000000 0.004031 0.00002806\n" + "33000000 0.00391 0.00002832\n" + "53000000 0.003868 0.00002935"; + +static const char * pvtoTableContent = "# Rs[sm3/sm3]\tPbub[Pa]\tBo[m3/sm3]\tVisc(Pa.s)\n" + "\n" + " 2\t 2000000\t 1.02\t 0.000975\n" + " 5\t 5000000\t 1.03\t 0.00091\n" + " 10\t 10000000\t1.04\t 0.00083\n" + " 15\t 20000000\t1.05\t 0.000695\n" + " 90000000\t1.03\t 0.000985 -- some line comment\n" + " 30\t 30000000\t1.07\t 0.000594\n" + " 40\t 40000000\t1.08\t 0.00051\n" + " 50000000\t1.07\t 0.000549 -- another one\n" + " 90000000\t1.06\t 0.00074\n" + " 50\t 50000000.7\t1.09\t 0.000449\n" + " 90000000.7\t1.08\t 0.000605"; + +static const char * pvtwTableContent = "#\tPref[Pa]\tBw[m3/sm3]\tCp[1/Pa]\t Visc[Pa.s]\n" + "\t30600000.1\t1.03\t\t0.00000000041\t0.0003"; + +class MultiFluidLiveOilTest : public MultiFluidTest< BlackOilFluid, 3, 3 > +{ +public: + static constexpr real64 relTol = 1.0e-4; + static constexpr real64 absTol = 1.0e-4; +public: + MultiFluidLiveOilTest() + { + writeTableToFile( pvtoFileName, pvtoTableContent ); + writeTableToFile( pvdgFileName, pvdgTableContent ); + writeTableToFile( pvtwFileName, pvtwTableContent ); + + m_parent.resize( 1 ); + m_model = makeLiveOilFluid( "fluid", &m_parent ); + + m_parent.initialize(); + m_parent.initializePostInitialConditions(); + } + + ~MultiFluidLiveOilTest() override + { + removeFile( pvtoFileName ); + removeFile( pvdgFileName ); + removeFile( pvtwFileName ); + } + +private: + static BlackOilFluid * makeLiveOilFluid( string const & name, Group * parent ); + static constexpr const char * pvtoFileName = "pvto.txt"; + static constexpr const char * pvdgFileName = "pvdg.txt"; + static constexpr const char * pvtwFileName = "pvtw.txt"; +}; + +BlackOilFluid * MultiFluidLiveOilTest::makeLiveOilFluid( string const & name, Group * parent ) +{ + BlackOilFluid & fluid = parent->registerGroup< BlackOilFluid >( name ); + + string_array & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); + fill< 3 >( phaseNames, {"oil", "gas", "water"} ); + + string_array & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); + fill< 3 >( compNames, {"oil", "gas", "water"} ); + + array1d< real64 > & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); + fill< 3 >( molarWgt, {114e-3, 16e-3, 18e-3} ); + + array1d< real64 > & surfaceDens = fluid.getReference< array1d< real64 > >( BlackOilFluidBase::viewKeyStruct::surfacePhaseMassDensitiesString() ); + fill< 3 >( surfaceDens, {800.0, 0.9907, 1022.0} ); + + path_array & tableNames = fluid.getReference< path_array >( BlackOilFluidBase::viewKeyStruct::tableFilesString() ); + fill< 3 >( tableNames, {pvtoFileName, pvdgFileName, pvtwFileName} ); + + fluid.postProcessInputRecursive(); + return &fluid; +} + +TEST_F( MultiFluidLiveOilTest, numericalDerivativesMolar ) +{ + auto & fluid = getFluid(); + fluid.setMassFlag( false ); + + real64 constexpr eps = 1.0e-6; + + array2d< real64 > samples( 2, 3 ); + fill< 3 >( samples[0], { 0.10000, 0.3, 0.60000 } ); + fill< 3 >( samples[1], { 0.79999, 0.2, 0.00001 } ); + + for( real64 const pressure : { 1.24e7, 3.21e7, 5.01e7 } ) + { + for( integer sampleIndex = 0; sampleIndex < samples.size( 0 ); sampleIndex++ ) + { + TestData data ( pressure, 297.15, samples[sampleIndex] ); + testNumericalDerivatives( fluid, &getParent(), data, eps, relTol, absTol ); + } + } +} + +TEST_F( MultiFluidLiveOilTest, numericalDerivativesMass ) +{ + auto & fluid = getFluid(); + fluid.setMassFlag( true ); + + real64 constexpr eps = 1.0e-6; + + array2d< real64 > samples( 2, 3 ); + fill< 3 >( samples[0], { 0.10000, 0.3, 0.60000 } ); + fill< 3 >( samples[1], { 0.79999, 0.2, 0.00001 } ); + + for( real64 const pressure : { 1.24e7, 3.21e7, 5.01e7 } ) + { + for( integer sampleIndex = 0; sampleIndex < samples.size( 0 ); sampleIndex++ ) + { + TestData data ( pressure, 297.15, samples[sampleIndex] ); + testNumericalDerivatives( fluid, &getParent(), data, eps, relTol, absTol ); + } + } +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + + geos::GeosxState state( geos::basicSetup( argc, argv ) ); + + int const result = RUN_ALL_TESTS(); + + geos::basicCleanup(); + + return result; +} diff --git a/src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp b/src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp index 5b336a0481a..07447f822cd 100644 --- a/src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp +++ b/src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp @@ -24,6 +24,7 @@ #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "testFlowUtils.hpp" namespace geos { @@ -34,136 +35,6 @@ namespace testing using namespace geos::constitutive; using namespace geos::constitutive::multifluid; -void checkDerivative( real64 const valueEps, - real64 const value, - real64 const deriv, - real64 const eps, - real64 const relTol, - real64 const absTol, - string const & name, - string const & var ) -{ - real64 const numDeriv = (valueEps - value) / eps; - checkRelativeError( deriv, numDeriv, relTol, absTol, "d(" + name + ")/d(" + var + ")" ); -} - -void checkDerivative( real64 const valueEps, - real64 const value, - real64 const deriv, - real64 const eps, - real64 const relTol, - string const & name, - string const & var ) -{ return checkDerivative( valueEps, value, deriv, eps, relTol, DEFAULT_ABS_TOL, name, var ); } - -template< int USD1, int USD2, int USD3 > -void checkDerivative( arraySlice1d< real64 const, USD1 > const & valueEps, - arraySlice1d< real64 const, USD2 > const & value, - arraySlice1d< real64 const, USD3 > const & deriv, - real64 const eps, - real64 const relTol, - real64 const absTol, - string const & name, - string const & var, - arrayView1d< string const > const & labels ) -{ - localIndex const size = labels.size( 0 ); - - for( localIndex i = 0; i < size; ++i ) - { - checkDerivative( valueEps[i], value[i], deriv[i], eps, relTol, absTol, - name + "[" + labels[i] + "]", var ); - } -} - -template< int DIM, int USD1, int USD2, int USD3, typename ... Args > -void checkDerivative( ArraySlice< real64 const, DIM, USD1 > const & valueEps, - ArraySlice< real64 const, DIM, USD2 > const & value, - ArraySlice< real64 const, DIM, USD3 > const & deriv, - real64 const eps, - real64 const relTol, - real64 const absTol, - string const & name, - string const & var, - arrayView1d< string const > const & labels, - Args ... label_lists ) -{ - localIndex const size = labels.size( 0 ); - - for( localIndex i = 0; i < size; ++i ) - { - checkDerivative( valueEps[i], value[i], deriv[i], eps, relTol, absTol, - name + "[" + labels[i] + "]", var, label_lists ... ); - } -} - -template< int DIM, int USD1, int USD2, int USD3, typename ... Args > -void checkDerivative( ArraySlice< real64 const, DIM, USD1 > const & valueEps, - ArraySlice< real64 const, DIM, USD2 > const & value, - ArraySlice< real64 const, DIM, USD3 > const & deriv, - real64 const eps, - real64 const relTol, - string const & name, - string const & var, - arrayView1d< string const > const & labels, - Args ... label_lists ) -{ return checkDerivative( valueEps, value, deriv, eps, relTol, DEFAULT_ABS_TOL, name, var, labels, label_lists ... ); } - -// invert compositional derivative array layout to move innermost slice on the top -// (this is needed so we can use checkDerivative() to check derivative w.r.t. for each compositional var) -template< int USD > -array1d< real64 > invertLayout( arraySlice1d< real64 const, USD > const & input, - localIndex N ) -{ - array1d< real64 > output( N ); - for( int i = 0; i < N; ++i ) - { - output[i] = input[i]; - } - - return output; -} - -template< int USD > -array2d< real64 > invertLayout( arraySlice2d< real64 const, USD > const & input, - localIndex N1, - localIndex N2 ) -{ - array2d< real64 > output( N2, N1 ); - - for( localIndex i = 0; i < N1; ++i ) - { - for( localIndex j = 0; j < N2; ++j ) - { - output( j, i ) = input( i, j ); - } - } - - return output; -} - -template< int USD > -array3d< real64 > invertLayout( arraySlice3d< real64 const, USD > const & input, - localIndex N1, - localIndex N2, - localIndex N3 ) -{ - array3d< real64 > output( N3, N1, N2 ); - - for( localIndex i = 0; i < N1; ++i ) - { - for( localIndex j = 0; j < N2; ++j ) - { - for( localIndex k = 0; k < N3; ++k ) - { - output( k, i, j ) = input( i, j, k ); - } - } - } - - return output; -} - void fillNumericalJacobian( arrayView1d< real64 const > const & residual, arrayView1d< real64 const > const & residualOrig, globalIndex const dofIndex, diff --git a/src/coreComponents/unitTests/fluidFlowTests/testFlowUtils.hpp b/src/coreComponents/unitTests/fluidFlowTests/testFlowUtils.hpp new file mode 100644 index 00000000000..6dd42bda36d --- /dev/null +++ b/src/coreComponents/unitTests/fluidFlowTests/testFlowUtils.hpp @@ -0,0 +1,160 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_UNITTESTS_FLUIDFLOWTESTS_TESTFLOWUTILS_HPP +#define GEOS_UNITTESTS_FLUIDFLOWTESTS_TESTFLOWUTILS_HPP + +#include "codingUtilities/UnitTestUtilities.hpp" + +namespace geos +{ + +namespace testing +{ + +void checkDerivative( real64 const valueEps, + real64 const value, + real64 const deriv, + real64 const eps, + real64 const relTol, + real64 const absTol, + string const & name, + string const & var ) +{ + real64 const numDeriv = (valueEps - value) / eps; + checkRelativeError( deriv, numDeriv, relTol, absTol, "d(" + name + ")/d(" + var + ")" ); +} + +void checkDerivative( real64 const valueEps, + real64 const value, + real64 const deriv, + real64 const eps, + real64 const relTol, + string const & name, + string const & var ) +{ return checkDerivative( valueEps, value, deriv, eps, relTol, DEFAULT_ABS_TOL, name, var ); } + +template< int USD1, int USD2, int USD3 > +void checkDerivative( arraySlice1d< real64 const, USD1 > const & valueEps, + arraySlice1d< real64 const, USD2 > const & value, + arraySlice1d< real64 const, USD3 > const & deriv, + real64 const eps, + real64 const relTol, + real64 const absTol, + string const & name, + string const & var, + arrayView1d< string const > const & labels ) +{ + localIndex const size = labels.size( 0 ); + + for( localIndex i = 0; i < size; ++i ) + { + checkDerivative( valueEps[i], value[i], deriv[i], eps, relTol, absTol, + name + "[" + labels[i] + "]", var ); + } +} + +template< int DIM, int USD1, int USD2, int USD3, typename ... Args > +void checkDerivative( ArraySlice< real64 const, DIM, USD1 > const & valueEps, + ArraySlice< real64 const, DIM, USD2 > const & value, + ArraySlice< real64 const, DIM, USD3 > const & deriv, + real64 const eps, + real64 const relTol, + real64 const absTol, + string const & name, + string const & var, + arrayView1d< string const > const & labels, + Args ... label_lists ) +{ + localIndex const size = labels.size( 0 ); + + for( localIndex i = 0; i < size; ++i ) + { + checkDerivative( valueEps[i], value[i], deriv[i], eps, relTol, absTol, + name + "[" + labels[i] + "]", var, label_lists ... ); + } +} + +template< int DIM, int USD1, int USD2, int USD3, typename ... Args > +void checkDerivative( ArraySlice< real64 const, DIM, USD1 > const & valueEps, + ArraySlice< real64 const, DIM, USD2 > const & value, + ArraySlice< real64 const, DIM, USD3 > const & deriv, + real64 const eps, + real64 const relTol, + string const & name, + string const & var, + arrayView1d< string const > const & labels, + Args ... label_lists ) +{ return checkDerivative( valueEps, value, deriv, eps, relTol, DEFAULT_ABS_TOL, name, var, labels, label_lists ... ); } + +// invert compositional derivative array layout to move innermost slice on the top +// (this is needed so we can use checkDerivative() to check derivative w.r.t. for each compositional var) +template< int USD > +array1d< real64 > invertLayout( arraySlice1d< real64 const, USD > const & input, + localIndex N ) +{ + array1d< real64 > output( N ); + for( int i = 0; i < N; ++i ) + { + output[i] = input[i]; + } + + return output; +} + +template< int USD > +array2d< real64 > invertLayout( arraySlice2d< real64 const, USD > const & input, + localIndex N1, + localIndex N2 ) +{ + array2d< real64 > output( N2, N1 ); + + for( localIndex i = 0; i < N1; ++i ) + { + for( localIndex j = 0; j < N2; ++j ) + { + output( j, i ) = input( i, j ); + } + } + + return output; +} + +template< int USD > +array3d< real64 > invertLayout( arraySlice3d< real64 const, USD > const & input, + localIndex N1, + localIndex N2, + localIndex N3 ) +{ + array3d< real64 > output( N3, N1, N2 ); + + for( localIndex i = 0; i < N1; ++i ) + { + for( localIndex j = 0; j < N2; ++j ) + { + for( localIndex k = 0; k < N3; ++k ) + { + output( k, i, j ) = input( i, j, k ); + } + } + } + + return output; +} + +} // namespace testing + +} // namespace geos + +#endif //GEOS_UNITTESTS_FLUIDFLOWTESTS_TESTFLOWUTILS_HPP diff --git a/src/coreComponents/unitTests/fluidFlowTests/testSingleFlowUtils.hpp b/src/coreComponents/unitTests/fluidFlowTests/testSingleFlowUtils.hpp index 2d9a1ee6ea1..0477b968d6a 100644 --- a/src/coreComponents/unitTests/fluidFlowTests/testSingleFlowUtils.hpp +++ b/src/coreComponents/unitTests/fluidFlowTests/testSingleFlowUtils.hpp @@ -24,6 +24,7 @@ #include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseFVM.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "testFlowUtils.hpp" namespace geos { @@ -33,29 +34,6 @@ namespace testing using namespace geos::constitutive; -void checkDerivative( real64 const valueEps, - real64 const value, - real64 const deriv, - real64 const eps, - real64 const relTol, - real64 const absTol, - string const & name, - string const & var ) -{ - real64 const numDeriv = (valueEps - value) / eps; - checkRelativeError( deriv, numDeriv, relTol, absTol, "d(" + name + ")/d(" + var + ")" ); -} - -void checkDerivative( real64 const valueEps, - real64 const value, - real64 const deriv, - real64 const eps, - real64 const relTol, - string const & name, - string const & var ) -{ return checkDerivative( valueEps, value, deriv, eps, relTol, DEFAULT_ABS_TOL, name, var ); } - - void fillNumericalJacobian( arrayView1d< real64 const > const & residual, arrayView1d< real64 const > const & residualOrig, globalIndex const dofIndex,