From cfef232ff723f7bde26672ce5da67720e43c20b5 Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Tue, 9 Feb 2021 10:34:23 -0800 Subject: [PATCH 1/6] implemented DO model --- .../constitutive/CMakeLists.txt | 2 + .../contact/ContactRelationBase.cpp | 16 +- .../constitutive/fluid/DeadOilFluid.cpp | 229 ++++++ .../constitutive/fluid/DeadOilFluid.hpp | 714 ++++++++++++++++++ .../constitutive/fluid/multiFluidSelector.hpp | 7 +- .../constitutive/unitTests/testMultiFluid.cpp | 171 +++++ .../fileIO/schema/docs/Constitutive.rst | 1 + .../fileIO/schema/docs/Constitutive_other.rst | 1 + .../fileIO/schema/docs/DeadOilFluid.rst | 21 + .../fileIO/schema/docs/DeadOilFluid_other.rst | 33 + src/coreComponents/fileIO/schema/schema.xsd | 27 + .../fileIO/schema/schema.xsd.other | 53 ++ .../managers/Functions/TableFunction.cpp | 229 +++--- .../managers/Functions/TableFunction.hpp | 314 +++++++- .../managers/unitTests/testFunctions.cpp | 234 ++++-- src/docs/sphinx/CompleteXMLSchema.rst | 14 + 16 files changed, 1866 insertions(+), 200 deletions(-) create mode 100644 src/coreComponents/constitutive/fluid/DeadOilFluid.cpp create mode 100644 src/coreComponents/constitutive/fluid/DeadOilFluid.hpp create mode 100644 src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst create mode 100644 src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 9c648804f49..f7d33c12f69 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -13,6 +13,7 @@ set( constitutive_headers capillaryPressure/VanGenuchtenCapillaryPressure.hpp contact/ContactRelationBase.hpp contact/Coulomb.hpp + fluid/DeadOilFluid.hpp fluid/MultiPhaseMultiComponentFluid.hpp fluid/PVTFunctions/BrineCO2DensityFunction.hpp fluid/PVTFunctions/BrineViscosityFunction.hpp @@ -68,6 +69,7 @@ set( constitutive_sources capillaryPressure/VanGenuchtenCapillaryPressure.cpp contact/ContactRelationBase.cpp contact/Coulomb.cpp + fluid/DeadOilFluid.cpp fluid/CompressibleSinglePhaseFluid.cpp fluid/MultiPhaseMultiComponentFluid.cpp fluid/PVTFunctions/BrineCO2DensityFunction.cpp diff --git a/src/coreComponents/constitutive/contact/ContactRelationBase.cpp b/src/coreComponents/constitutive/contact/ContactRelationBase.cpp index f6548c9271b..26386acf29e 100644 --- a/src/coreComponents/constitutive/contact/ContactRelationBase.cpp +++ b/src/coreComponents/constitutive/contact/ContactRelationBase.cpp @@ -107,15 +107,16 @@ void ContactRelationBase::initializePreSubGroups( Group * const ) TableFunction * const apertureTable = dynamic_cast< TableFunction * >(m_apertureFunction); if( apertureTable!=nullptr ) { - array1d< array1d< real64 > > & xvals0 = apertureTable->getCoordinates(); + ArrayOfArraysView< real64 > xvals0 = apertureTable->getCoordinates(); array1d< real64 > & yvals = apertureTable->getValues(); GEOSX_ERROR_IF( xvals0.size() > 1, "Aperture limiter table cannot be greater than a 1d table." ); - array1d< real64 > & xvals = xvals0[0]; + arraySlice1d< real64 > xvals = xvals0[0]; + localIndex const size = xvals.size(); - GEOSX_ERROR_IF( xvals.back() > 0.0 || xvals.back() < 0.0, + GEOSX_ERROR_IF( xvals0( 0, size-1 ) > 0.0 || xvals0( 0, size-1 ) < 0.0, "Invalid aperture limiter table. Last coordinate must be zero!!" ); GEOSX_ERROR_IF( xvals.size() < 2, @@ -129,20 +130,15 @@ void ContactRelationBase::initializePreSubGroups( Group * const ) real64 m_apertureTransition = (yvals[n] - slope * xvals[n] ) / ( 1.0 - slope ); - xvals.emplace_back( m_apertureTransition ); + xvals0.emplaceBack( 0, m_apertureTransition ); yvals.emplace_back( m_apertureTransition ); - xvals.emplace_back( m_apertureTransition*10e9 ); + xvals0.emplaceBack( 0, m_apertureTransition*10e9 ); yvals.emplace_back( m_apertureTransition*10e9 ); apertureTable->reInitializeFunction(); } -// for( int i=0 ; i<200 ; ++i ) -// { -// real64 coord = 0.01*i-1.0; -// std::cout<Evaluate( &coord )< const phaseDict = +{ + { "gas", DeadOilFluid::PhaseType::GAS }, + { "oil", DeadOilFluid::PhaseType::OIL }, + { "water", DeadOilFluid::PhaseType::WATER } +}; + +} + +DeadOilFluid::DeadOilFluid( string const & name, + Group * const parent ) + : + MultiFluidBase( name, parent ), + m_waterRefPressure( 0.0 ), + m_waterFormationVolFactor( 0.0 ), + m_waterCompressibility( 0.0 ), + m_waterViscosity( 0.0 ) +{ + getWrapperBase( viewKeyStruct::componentMolarWeightString )->setInputFlag( InputFlags::REQUIRED ); + getWrapperBase( viewKeyStruct::phaseNamesString )->setInputFlag( InputFlags::REQUIRED ); + + registerWrapper( viewKeyStruct::surfacePhaseMassDensitiesString, &m_surfacePhaseMassDensity )-> + setInputFlag( InputFlags::REQUIRED )-> + setDescription( "List of surface mass densities for each phase" ); + + registerWrapper( viewKeyStruct::formationVolumeFactorTableNamesString, &m_formationVolFactorTableNames )-> + setInputFlag( InputFlags::REQUIRED )-> + setDescription( "Names of the formation volume factor tables (one per hydrocarbon phase, in the order provided in \"phaseNames\"). \n" + "For instance, if \"oil\" is before \"gas\" in \"phaseNames\", the table order should be: oilTableName, gasTableName" ); + + registerWrapper( viewKeyStruct::viscosityTableNamesString, &m_viscosityTableNames )-> + setInputFlag( InputFlags::REQUIRED )-> + setDescription( "Names of the viscosity tables (one per hydrocarbon phase, in the order provided in \"phaseNames\") \n" + "For instance, if \"oil\" is before \"gas\" in \"phaseNames\", the table order should be: oilTableName, gasTableName" ); + + registerWrapper( viewKeyStruct::waterRefPressureString, &m_waterRefPressure )-> + setInputFlag( InputFlags::OPTIONAL )-> + setDescription( "Water reference pressure" ); + + registerWrapper( viewKeyStruct::waterFormationVolumeFactorString, &m_waterFormationVolFactor )-> + setInputFlag( InputFlags::OPTIONAL )-> + setDescription( "Water formation volume factor" ); + + registerWrapper( viewKeyStruct::waterCompressibilityString, &m_waterCompressibility )-> + setInputFlag( InputFlags::OPTIONAL )-> + setDescription( "Water compressibility" ); + + registerWrapper( viewKeyStruct::waterViscosityString, &m_waterViscosity )-> + setInputFlag( InputFlags::OPTIONAL )-> + setDescription( "Water viscosity" ); + +} + +DeadOilFluid::~DeadOilFluid() +{ } + + +void DeadOilFluid::postProcessInput() +{ + m_componentNames = m_phaseNames; + + MultiFluidBase::postProcessInput(); + + m_phaseTypes.resize( numFluidPhases() ); + m_phaseOrder.resize( PhaseType::MAX_NUM_PHASES ); + m_phaseOrder.setValues< serialPolicy >( -1 ); + + for( localIndex ip = 0; ip < numFluidPhases(); ++ip ) + { + auto it = phaseDict.find( m_phaseNames[ip] ); + GEOSX_ERROR_IF( it == phaseDict.end(), "DeadOilFluid: phase not supported: " << m_phaseNames[ip] ); + integer const phaseIndex = it->second; + GEOSX_ERROR_IF( phaseIndex >= PhaseType::MAX_NUM_PHASES, "DeadOilFluid: invalid phase index " << phaseIndex ); + + m_phaseTypes[ip] = phaseIndex; + m_phaseOrder[phaseIndex] = LvArray::integerConversion< integer >( ip ); + if( phaseIndex == PhaseType::OIL || phaseIndex == PhaseType::GAS ) + { + m_hydrocarbonPhaseOrder.emplace_back( LvArray::integerConversion< integer >( ip ) ); + } + } + + // check the water phase parameters + integer const ipWater = m_phaseOrder[PhaseType::WATER]; + integer const ipGas = m_phaseOrder[PhaseType::GAS]; + if( ipWater >= 0 ) // if water is present + { + GEOSX_ERROR_IF( m_waterRefPressure <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " << viewKeyStruct::waterRefPressureString ); + GEOSX_ERROR_IF( m_waterFormationVolFactor <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " << viewKeyStruct::waterFormationVolumeFactorString ); + GEOSX_ERROR_IF( m_waterViscosity <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " << viewKeyStruct::waterViscosityString ); + GEOSX_ERROR_IF( m_waterCompressibility <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " << viewKeyStruct::waterCompressibilityString ); + } + else + { + GEOSX_ERROR_IF( m_waterRefPressure > 0.0, + "DeadOilFluid: if water is absent, this keyword is not needed " << viewKeyStruct::waterRefPressureString ); + GEOSX_ERROR_IF( m_waterFormationVolFactor > 0.0, + "DeadOilFluid: if water is absent, this keyword is not needed " << viewKeyStruct::waterFormationVolumeFactorString ); + GEOSX_ERROR_IF( m_waterViscosity > 0.0, + "DeadOilFluid: if water is absent, this keyword is not needed " << viewKeyStruct::waterViscosityString ); + GEOSX_ERROR_IF( m_waterCompressibility > 0.0, + "DeadOilFluid: if water is absent, this keyword is not needed " << viewKeyStruct::waterCompressibilityString ); + } + + // check the table names + localIndex const numExpectedTables = (ipGas >= 0) ? 2 : 1; + GEOSX_ERROR_IF( m_formationVolFactorTableNames.size() != numExpectedTables, + "DeadOilFluid: one formation volume factor table must be provided for each hydrocarbon phase" ); + GEOSX_ERROR_IF( m_viscosityTableNames.size() != numExpectedTables, + "DeadOilFluid: one viscosity table must be provided for each hydrocarbon phase" ); + + // check the size of the additional parameters + GEOSX_ERROR_IF( m_surfacePhaseMassDensity.size() != m_phaseNames.size(), + "DeadOilFluid: the number of surfacePhaseMassDensities is inconsistent with the phase names" ); + GEOSX_ERROR_IF( m_componentMolarWeight.size() != m_phaseNames.size(), + "DeadOilFluid: the number of componentMolarWeights is inconsistent with the phase names" ); +} + +void DeadOilFluid::initializePostSubGroups( Group * const group ) +{ + MultiFluidBase::initializePostSubGroups( group ); + createAllKernelWrappers(); +} + +void DeadOilFluid::createAllKernelWrappers() +{ + FunctionManager const & functionManager = FunctionManager::instance(); + + GEOSX_ERROR_IF( m_hydrocarbonPhaseOrder.size() != 1 && m_hydrocarbonPhaseOrder.size() != 2, + "DeadOilFluid: the number of hydrocarbon phases should be equal to 1 (oil) or 2 (oil+gas)" ); + + if( m_formationVolFactorTables.size() == 0 && m_viscosityTables.size() == 0 ) + { + + // loop over the hydrocarbon phases + for( localIndex iph = 0; iph < m_hydrocarbonPhaseOrder.size(); ++iph ) + { + localIndex const ip = m_hydrocarbonPhaseOrder[iph]; + + // grab the tables by name from the function manager + TableFunction const & FVFTable = *functionManager.getGroup< TableFunction const >( m_formationVolFactorTableNames[ip] ); + TableFunction const & viscosityTable = *functionManager.getGroup< TableFunction const >( m_viscosityTableNames[ip] ); + validateTable( FVFTable ); + validateTable( viscosityTable ); + + // create the table wrapper for the oil and (if present) the gas phases + m_formationVolFactorTables.emplace_back( FVFTable.createKernelWrapper() ); + m_viscosityTables.emplace_back( viscosityTable.createKernelWrapper() ); + } + } + +} + +void DeadOilFluid::validateTable( TableFunction const & table ) const +{ + array1d< real64 > const & property = table.getValues(); + + GEOSX_ERROR_IF( table.getInterpolationMethod() != TableFunction::InterpolationType::Linear, + "DeadOilFluid: the interpolation method for the PVT tables must be linear" ); + GEOSX_ERROR_IF( property.size() <= 1, + "DeadOilFluid: each PVT table must contain at least 2 values" ); + + for( localIndex i = 2; i < property.size(); ++i ) + { + GEOSX_ERROR_IF( (property[i] - property[i-1]) * (property[i-1] - property[i-2]) < 0, + "DeadOilFluid: the values in each PVT table must monotone" ); + } +} + +std::unique_ptr< ConstitutiveBase > +DeadOilFluid::deliverClone( string const & name, + Group * const parent ) const +{ + std::unique_ptr< ConstitutiveBase > clone = MultiFluidBase::deliverClone( name, parent ); + + DeadOilFluid & model = dynamicCast< DeadOilFluid & >( *clone ); + + model.m_phaseTypes = m_phaseTypes; + model.m_phaseOrder = m_phaseOrder; + model.m_hydrocarbonPhaseOrder = m_hydrocarbonPhaseOrder; + + model.createAllKernelWrappers(); + + return clone; +} + + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, DeadOilFluid, string const &, Group * const ) +} //namespace constitutive + +} //namespace geosx diff --git a/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp new file mode 100644 index 00000000000..f682e1ffce1 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp @@ -0,0 +1,714 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file DeadOilFluid.hpp + */ + +#ifndef GEOSX_CONSTITUTIVE_FLUID_DEADOILFLUID_HPP_ +#define GEOSX_CONSTITUTIVE_FLUID_DEADOILFLUID_HPP_ + +#include "constitutive/fluid/MultiFluidBase.hpp" +#include "managers/Functions/TableFunction.hpp" + +namespace geosx +{ + +namespace constitutive +{ + +/** + * @brief Kernel wrapper class for DeadOilFluild + */ +class DeadOilFluidUpdate final : public MultiFluidBaseUpdate +{ +public: + + DeadOilFluidUpdate( arrayView1d< integer const > const & phaseTypes, + arrayView1d< integer const > const & phaseOrder, + arrayView1d< integer const > const & hydrocarbonPhaseOrder, + arrayView1d< real64 const > const & surfacePhaseMassDensity, + arrayView1d< TableFunction::KernelWrapper const > const & formationVolFactorTables, + arrayView1d< TableFunction::KernelWrapper const > const & viscosityTables, + real64 const waterRefPressure, + real64 const waterFormationVolFactor, + real64 const waterCompressibility, + real64 const waterViscosity, + arrayView1d< real64 const > const & componentMolarWeight, + bool useMass, + arrayView3d< real64 > const & phaseFraction, + arrayView3d< real64 > const & dPhaseFraction_dPressure, + arrayView3d< real64 > const & dPhaseFraction_dTemperature, + arrayView4d< real64 > const & dPhaseFraction_dGlobalCompFraction, + arrayView3d< real64 > const & phaseDensity, + arrayView3d< real64 > const & dPhaseDensity_dPressure, + arrayView3d< real64 > const & dPhaseDensity_dTemperature, + arrayView4d< real64 > const & dPhaseDensity_dGlobalCompFraction, + arrayView3d< real64 > const & phaseMassDensity, + arrayView3d< real64 > const & dPhaseMassDensity_dPressure, + arrayView3d< real64 > const & dPhaseMassDensity_dTemperature, + arrayView4d< real64 > const & dPhaseMassDensity_dGlobalCompFraction, + arrayView3d< real64 > const & phaseViscosity, + arrayView3d< real64 > const & dPhaseViscosity_dPressure, + arrayView3d< real64 > const & dPhaseViscosity_dTemperature, + arrayView4d< real64 > const & dPhaseViscosity_dGlobalCompFraction, + arrayView4d< real64 > const & phaseCompFraction, + arrayView4d< real64 > const & dPhaseCompFraction_dPressure, + arrayView4d< real64 > const & dPhaseCompFraction_dTemperature, + arrayView5d< real64 > const & dPhaseCompFraction_dGlobalCompFraction, + arrayView2d< real64 > const & totalDensity, + arrayView2d< real64 > const & dTotalDensity_dPressure, + arrayView2d< real64 > const & dTotalDensity_dTemperature, + arrayView3d< real64 > const & dTotalDensity_dGlobalCompFraction ) + : MultiFluidBaseUpdate( componentMolarWeight, + useMass, + phaseFraction, + dPhaseFraction_dPressure, + dPhaseFraction_dTemperature, + dPhaseFraction_dGlobalCompFraction, + phaseDensity, + dPhaseDensity_dPressure, + dPhaseDensity_dTemperature, + dPhaseDensity_dGlobalCompFraction, + phaseMassDensity, + dPhaseMassDensity_dPressure, + dPhaseMassDensity_dTemperature, + dPhaseMassDensity_dGlobalCompFraction, + phaseViscosity, + dPhaseViscosity_dPressure, + dPhaseViscosity_dTemperature, + dPhaseViscosity_dGlobalCompFraction, + phaseCompFraction, + dPhaseCompFraction_dPressure, + dPhaseCompFraction_dTemperature, + dPhaseCompFraction_dGlobalCompFraction, + totalDensity, + dTotalDensity_dPressure, + dTotalDensity_dTemperature, + dTotalDensity_dGlobalCompFraction ), + m_phaseTypes( phaseTypes ), + m_phaseOrder( phaseOrder ), + m_hydrocarbonPhaseOrder( hydrocarbonPhaseOrder ), + m_surfacePhaseMassDensity( surfacePhaseMassDensity ), + m_formationVolFactorTables( formationVolFactorTables ), + m_viscosityTables( viscosityTables ), + m_waterRefPressure( waterRefPressure ), + m_waterFormationVolFactor( waterFormationVolFactor ), + m_waterCompressibility( waterCompressibility ), + m_waterViscosity( waterViscosity ) + {} + + /// Default copy constructor + DeadOilFluidUpdate( DeadOilFluidUpdate const & ) = default; + + /// Default move constructor + DeadOilFluidUpdate( DeadOilFluidUpdate && ) = default; + + /// Deleted copy assignment operator + DeadOilFluidUpdate & operator=( DeadOilFluidUpdate const & ) = delete; + + /// Deleted move assignment operator + DeadOilFluidUpdate & operator=( DeadOilFluidUpdate && ) = delete; + + GEOSX_HOST_DEVICE + GEOSX_FORCE_INLINE + virtual void compute( real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition, + arraySlice1d< real64 > const & phaseFraction, + arraySlice1d< real64 > const & phaseDensity, + arraySlice1d< real64 > const & phaseMassDensity, + arraySlice1d< real64 > const & phaseViscosity, + arraySlice2d< real64 > const & phaseCompFraction, + real64 & totalDensity ) const override; + + GEOSX_HOST_DEVICE + GEOSX_FORCE_INLINE + virtual void compute( real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition, + arraySlice1d< real64 > const & phaseFraction, + arraySlice1d< real64 > const & dPhaseFraction_dPressure, + arraySlice1d< real64 > const & dPhaseFraction_dTemperature, + arraySlice2d< real64 > const & dPhaseFraction_dGlobalCompFraction, + arraySlice1d< real64 > const & phaseDensity, + arraySlice1d< real64 > const & dPhaseDensity_dPressure, + arraySlice1d< real64 > const & dPhaseDensity_dTemperature, + arraySlice2d< real64 > const & dPhaseDensity_dGlobalCompFraction, + arraySlice1d< real64 > const & phaseMassDensity, + arraySlice1d< real64 > const & dPhaseMassDensity_dPressure, + arraySlice1d< real64 > const & dPhaseMassDensity_dTemperature, + arraySlice2d< real64 > const & dPhaseMassDensity_dGlobalCompFraction, + arraySlice1d< real64 > const & phaseViscosity, + arraySlice1d< real64 > const & dPhaseViscosity_dPressure, + arraySlice1d< real64 > const & dPhaseViscosity_dTemperature, + arraySlice2d< real64 > const & dPhaseViscosity_dGlobalCompFraction, + arraySlice2d< real64 > const & phaseCompFraction, + arraySlice2d< real64 > const & dPhaseCompFraction_dPressure, + arraySlice2d< real64 > const & dPhaseCompFraction_dTemperature, + arraySlice3d< real64 > const & dPhaseCompFraction_dGlobalCompFraction, + real64 & totalDensity, + real64 & dTotalDensity_dPressure, + real64 & dTotalDensity_dTemperature, + arraySlice1d< real64 > const & dTotalDensity_dGlobalCompFraction ) const override; + + GEOSX_HOST_DEVICE + GEOSX_FORCE_INLINE + virtual void update( localIndex const k, + localIndex const q, + real64 const pressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition ) const override + { + compute( pressure, + temperature, + composition, + m_phaseFraction[k][q], + m_dPhaseFraction_dPressure[k][q], + m_dPhaseFraction_dTemperature[k][q], + m_dPhaseFraction_dGlobalCompFraction[k][q], + m_phaseDensity[k][q], + m_dPhaseDensity_dPressure[k][q], + m_dPhaseDensity_dTemperature[k][q], + m_dPhaseDensity_dGlobalCompFraction[k][q], + m_phaseMassDensity[k][q], + m_dPhaseMassDensity_dPressure[k][q], + m_dPhaseMassDensity_dTemperature[k][q], + m_dPhaseMassDensity_dGlobalCompFraction[k][q], + m_phaseViscosity[k][q], + m_dPhaseViscosity_dPressure[k][q], + m_dPhaseViscosity_dTemperature[k][q], + m_dPhaseViscosity_dGlobalCompFraction[k][q], + m_phaseCompFraction[k][q], + m_dPhaseCompFraction_dPressure[k][q], + m_dPhaseCompFraction_dTemperature[k][q], + m_dPhaseCompFraction_dGlobalCompFraction[k][q], + m_totalDensity[k][q], + m_dTotalDensity_dPressure[k][q], + m_dTotalDensity_dTemperature[k][q], + m_dTotalDensity_dGlobalCompFraction[k][q] ); + } + +private: + + GEOSX_HOST_DEVICE + GEOSX_FORCE_INLINE + void computeDensities( real64 pressure, + arraySlice1d< real64, 0 > const & phaseMassDens ) const; + + GEOSX_HOST_DEVICE + GEOSX_FORCE_INLINE + void computeDensities( real64 pressure, + arraySlice1d< real64 > const & phaseMassDens, + arraySlice1d< real64 > const & dPhaseMassDens_dPres, + arraySlice2d< real64 > const & dPhaseMassDens_dGlobalCompFrac ) const; + + GEOSX_HOST_DEVICE + GEOSX_FORCE_INLINE + void computeViscosities( real64 pressure, + arraySlice1d< real64, 0 > const & phaseVisc ) const; + + GEOSX_HOST_DEVICE + GEOSX_FORCE_INLINE + void computeViscosities( real64 pressure, + arraySlice1d< real64 > const & phaseVisc, + arraySlice1d< real64 > const & dPhaseVisc_dPressure, + arraySlice2d< real64 > const & dPhaseVisc_dGlobalCompFraction ) const; + + /// Phase ordering info + arrayView1d< integer const > m_phaseTypes; + arrayView1d< integer const > m_phaseOrder; + arrayView1d< integer const > m_hydrocarbonPhaseOrder; + + /// Surface mass density for each phase + arrayView1d< real64 const > m_surfacePhaseMassDensity; + + /// Table kernel wrappers to interpolate in the oil and gas (B vs p) tables + arrayView1d< TableFunction::KernelWrapper const > m_formationVolFactorTables; + + /// Table kernel wrappers to interpolate in the oil and gas (\mu vs p) tables + arrayView1d< TableFunction::KernelWrapper const > m_viscosityTables; + + /// Water reference pressure + real64 const m_waterRefPressure; + + /// Water formation volume factor + real64 const m_waterFormationVolFactor; + + /// Water compressibility + real64 const m_waterCompressibility; + + /// Water viscosity + real64 const m_waterViscosity; + + +}; + +class DeadOilFluid : public MultiFluidBase +{ +public: + + struct PhaseType + { + static constexpr integer OIL = 0; + static constexpr integer GAS = 1; + static constexpr integer WATER = 2; + static constexpr integer MAX_NUM_PHASES = 3; + }; + + DeadOilFluid( string const & name, Group * const parent ); + + virtual ~DeadOilFluid() override; + + virtual std::unique_ptr< ConstitutiveBase > + deliverClone( string const & name, + Group * const parent ) const override; + + static string catalogName() { return "DeadOilFluid"; } + + virtual string getCatalogName() const override { return catalogName(); } + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = DeadOilFluidUpdate; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() + { + return KernelWrapper( m_phaseTypes, + m_phaseOrder, + m_hydrocarbonPhaseOrder, + m_surfacePhaseMassDensity, + m_formationVolFactorTables, + m_viscosityTables, + m_waterRefPressure, + m_waterFormationVolFactor, + m_waterCompressibility, + m_waterViscosity, + m_componentMolarWeight, + m_useMass, + m_phaseFraction, + m_dPhaseFraction_dPressure, + m_dPhaseFraction_dTemperature, + m_dPhaseFraction_dGlobalCompFraction, + m_phaseDensity, + m_dPhaseDensity_dPressure, + m_dPhaseDensity_dTemperature, + m_dPhaseDensity_dGlobalCompFraction, + m_phaseMassDensity, + m_dPhaseMassDensity_dPressure, + m_dPhaseMassDensity_dTemperature, + m_dPhaseMassDensity_dGlobalCompFraction, + m_phaseViscosity, + m_dPhaseViscosity_dPressure, + m_dPhaseViscosity_dTemperature, + m_dPhaseViscosity_dGlobalCompFraction, + m_phaseCompFraction, + m_dPhaseCompFraction_dPressure, + m_dPhaseCompFraction_dTemperature, + m_dPhaseCompFraction_dGlobalCompFraction, + m_totalDensity, + m_dTotalDensity_dPressure, + m_dTotalDensity_dTemperature, + m_dTotalDensity_dGlobalCompFraction ); + } + + struct viewKeyStruct : MultiFluidBase::viewKeyStruct + { + static constexpr auto surfacePhaseMassDensitiesString = "surfacePhaseMassDensities"; + static constexpr auto formationVolumeFactorTableNamesString = "hydrocarbonFormationVolFactorTableNames"; + static constexpr auto viscosityTableNamesString = "hydrocarbonViscosityTableNames"; + static constexpr auto waterRefPressureString = "waterReferencePressure"; + static constexpr auto waterFormationVolumeFactorString = "waterFormationVolumeFactor"; + static constexpr auto waterCompressibilityString = "waterCompressibility"; + static constexpr auto waterViscosityString = "waterViscosity"; + + } viewKeysDeadOilFluid; + + +protected: + + virtual void postProcessInput() override; + + virtual void initializePostSubGroups( Group * const group ) override; + +private: + + /// create all the table kernel wrappers + void createAllKernelWrappers(); + + /// check that the table makes sense + void validateTable( TableFunction const & table ) const; + + // Input data + + /// Names of the formation volume factor tables + array1d< string > m_formationVolFactorTableNames; + + /// Names of the viscosity tables + array1d< string > m_viscosityTableNames; + + /// Surface densities + array1d< real64 > m_surfacePhaseMassDensity; + + /// Water reference pressure + real64 m_waterRefPressure; + + /// Water formation volume factor + real64 m_waterFormationVolFactor; + + /// Water compressibility + real64 m_waterCompressibility; + + /// Water viscosity + real64 m_waterViscosity; + + + /// Data after processing of input + + /// Phase ordering info + array1d< integer > m_phaseTypes; + array1d< integer > m_phaseOrder; + array1d< integer > m_hydrocarbonPhaseOrder; + + /// Table kernel wrappers to interpolate in the oil and gas (B vs p) tables + array1d< TableFunction::KernelWrapper > m_formationVolFactorTables; + + /// Table kernel wrappers to interpolate in the oil and gas (\mu vs p) tables + array1d< TableFunction::KernelWrapper > m_viscosityTables; + +}; + +GEOSX_HOST_DEVICE +GEOSX_FORCE_INLINE +void DeadOilFluidUpdate::computeDensities( real64 pressure, + arraySlice1d< real64, 0 > const & phaseMassDens ) const +{ + real64 fvf = 0.0; + real64 derivative = 0.0; + + // 1. Hydrocarbon phases: look up in the formation vol factor tables + + for( localIndex iph = 0; iph < m_hydrocarbonPhaseOrder.size(); ++iph ) + { + // get the phase index + localIndex const ip = m_hydrocarbonPhaseOrder[iph]; + // interpolate in the table to get the phase formation vol factor, and discard the derivative + m_formationVolFactorTables[iph].compute( &pressure, fvf, &derivative ); + + // we are ready to update the densities + phaseMassDens[ip] = m_surfacePhaseMassDensity[ip] / fvf; + } + + // 2. Water phase: use the constant formation volume factor and compressibility provided by the user + + using PT = DeadOilFluid::PhaseType; + integer const ipWater = m_phaseOrder[PT::WATER]; + + // if water is present + if( ipWater >= 0 ) + { + // note: double check, but I think std::exp is allowed in kernel + real64 const denom = m_waterFormationVolFactor * std::exp( -m_waterCompressibility * ( pressure - m_waterRefPressure ) ); + phaseMassDens[ipWater] = m_surfacePhaseMassDensity[ipWater] / denom; + } +} + +GEOSX_HOST_DEVICE +GEOSX_FORCE_INLINE +void DeadOilFluidUpdate::computeDensities( real64 pressure, + arraySlice1d< real64 > const & phaseMassDens, + arraySlice1d< real64 > const & dPhaseMassDens_dPres, + arraySlice2d< real64 > const & dPhaseMassDens_dGlobalCompFrac ) const +{ + real64 fvf = 0.0; + real64 derivative = 0.0; + + for( double & val : dPhaseMassDens_dGlobalCompFrac ) + { + val = 0.0; + } + + // 1. Hydrocarbon phases: look up in the formation vol factor tables + + for( localIndex iph = 0; iph < m_hydrocarbonPhaseOrder.size(); ++iph ) + { + // get the phase index + localIndex const ip = m_hydrocarbonPhaseOrder[iph]; + // interpolate in the table to get the phase formation vol factor and its derivative wrt pressure + m_formationVolFactorTables[iph].compute( &pressure, fvf, &derivative ); + + // we are ready to update the densities + real64 const fvfInv = 1.0 / fvf; + phaseMassDens[ip] = m_surfacePhaseMassDensity[ip] * fvfInv; + dPhaseMassDens_dPres[ip] = -derivative * phaseMassDens[ip] * fvfInv; + } + + // 2. Water phase: use the constant formation volume factor and compressibility provided by the user + + using PT = DeadOilFluid::PhaseType; + integer const ipWater = m_phaseOrder[PT::WATER]; + + // if water is present + if( ipWater >= 0 ) + { + // note: double check, but I think std::exp is allowed in kernel + real64 const expCompDeltaPres = std::exp( -m_waterCompressibility * ( pressure - m_waterRefPressure ) ); + real64 const dExpCompDeltaPres_dPres = -m_waterCompressibility * expCompDeltaPres; + real64 const denom = m_waterFormationVolFactor * expCompDeltaPres; + real64 const dDenom_dPres = m_waterFormationVolFactor * dExpCompDeltaPres_dPres; + real64 const denomInv = 1.0 / denom; + phaseMassDens[ipWater] = m_surfacePhaseMassDensity[ipWater] * denomInv; + dPhaseMassDens_dPres[ipWater] = -dDenom_dPres * phaseMassDens[ipWater] * denomInv; + } +} + +GEOSX_HOST_DEVICE +GEOSX_FORCE_INLINE +void DeadOilFluidUpdate::computeViscosities( real64 pressure, + arraySlice1d< real64, 0 > const & phaseVisc ) const +{ + // this derivative will be discarded + real64 derivative = 0.0; + + // 1. Hydrocarbon phases: look up in the viscosity tables + + for( localIndex iph = 0; iph < m_hydrocarbonPhaseOrder.size(); ++iph ) + { + // get the phase index + localIndex const ip = m_hydrocarbonPhaseOrder[iph]; + // interpolate in the table to get the phase viscosity, and discard the derivative + m_viscosityTables[iph].compute( &pressure, phaseVisc[ip], &derivative ); + } + + // 2. Water phase: use the constant viscosity provided by the user + + using PT = DeadOilFluid::PhaseType; + integer const ipWater = m_phaseOrder[PT::WATER]; + + // if water is present + if( ipWater >= 0 ) + { + // just assign the viscosity value + phaseVisc[ipWater] = m_waterViscosity; + } +} + +GEOSX_HOST_DEVICE +GEOSX_FORCE_INLINE +void DeadOilFluidUpdate::computeViscosities( real64 pressure, + arraySlice1d< real64 > const & phaseVisc, + arraySlice1d< real64 > const & dPhaseVisc_dPres, + arraySlice2d< real64 > const & dPhaseVisc_dGlobalCompFrac ) const +{ + for( double & val : dPhaseVisc_dGlobalCompFrac ) + { + val = 0.0; + } + + // 1. Hydrocarbon phases: look up in the viscosity tables + + for( localIndex iph = 0; iph < m_hydrocarbonPhaseOrder.size(); ++iph ) + { + // get the phase index + localIndex const ip = m_hydrocarbonPhaseOrder[iph]; + // interpolate in the table to get the phase viscosity and derivatives + m_viscosityTables[iph].compute( &pressure, phaseVisc[ip], &(dPhaseVisc_dPres)[ip] ); + } + + // 2. Water phase: use the constant viscosity provided by the user + + using PT = DeadOilFluid::PhaseType; + integer const ipWater = m_phaseOrder[PT::WATER]; + + // if water is present + if( ipWater >= 0 ) + { + // just assign the viscosity value + phaseVisc[ipWater] = m_waterViscosity; + dPhaseVisc_dPres[ipWater] = 0.0; + } +} + +GEOSX_HOST_DEVICE +GEOSX_FORCE_INLINE +void DeadOilFluidUpdate::compute( real64 pressure, + real64 temperature, + arraySlice1d< real64 const, 0 > const & composition, + arraySlice1d< real64, 0 > const & phaseFrac, + arraySlice1d< real64, 0 > const & phaseDens, + arraySlice1d< real64, 0 > const & phaseMassDens, + arraySlice1d< real64, 0 > const & phaseVisc, + arraySlice2d< real64, 1 > const & phaseCompFrac, + real64 & totalDens ) const +{ + GEOSX_UNUSED_VAR( temperature ); + localIndex const NC = m_componentMolarWeight.size(); + localIndex const NP = m_phaseTypes.size(); + + // 1. Read viscosities and formation volume factors from tables, update mass densities + computeViscosities( pressure, phaseVisc ); + computeDensities( pressure, phaseMassDens ); + + // 2. Update phaseDens (mass density if useMass == 1, molar density otherwise) + for( localIndex ip = 0; ip < NP; ++ip ) + { + real64 const mult = m_useMass ? 1.0 : 1.0 / m_componentMolarWeight[ip]; + phaseDens[ip] = phaseMassDens[ip] * mult; + } + + // 3. Update remaining variables: phaseFrac, phaseCompFrac using Dead-Oil assumptions + for( localIndex ip = 0; ip < NP; ++ip ) + { + phaseFrac[ip] = composition[ip]; + for( localIndex ic = 0; ic < NC; ++ic ) + { + phaseCompFrac[ip][ic] = (ip == ic) ? 1.0 : 0.0; + } + } + + // 4. Compute total fluid mass/molar density + totalDens = 0.0; + + // 4.1. Sum mass/molar fraction/density ratio over all phases to get the inverse of density + for( localIndex ip = 0; ip < NP; ++ip ) + { + totalDens += phaseFrac[ip] / phaseDens[ip]; + } + // 4.2. Invert the previous quantity to get actual density + totalDens = 1.0 / totalDens; +} + +GEOSX_HOST_DEVICE +GEOSX_FORCE_INLINE +void DeadOilFluidUpdate::compute( real64 pressure, + real64 temperature, + arraySlice1d< real64 const, 0 > const & composition, + arraySlice1d< real64 > const & phaseFraction, + arraySlice1d< real64 > const & dPhaseFraction_dPressure, + arraySlice1d< real64 > const & dPhaseFraction_dTemperature, + arraySlice2d< real64 > const & dPhaseFraction_dGlobalCompFraction, + arraySlice1d< real64 > const & phaseDensity, + arraySlice1d< real64 > const & dPhaseDensity_dPressure, + arraySlice1d< real64 > const & dPhaseDensity_dTemperature, + arraySlice2d< real64 > const & dPhaseDensity_dGlobalCompFraction, + arraySlice1d< real64 > const & phaseMassDensity, + arraySlice1d< real64 > const & dPhaseMassDensity_dPressure, + arraySlice1d< real64 > const & dPhaseMassDensity_dTemperature, + arraySlice2d< real64 > const & dPhaseMassDensity_dGlobalCompFraction, + arraySlice1d< real64 > const & phaseViscosity, + arraySlice1d< real64 > const & dPhaseViscosity_dPressure, + arraySlice1d< real64 > const & dPhaseViscosity_dTemperature, + arraySlice2d< real64 > const & dPhaseViscosity_dGlobalCompFraction, + arraySlice2d< real64 > const & phaseCompFraction, + arraySlice2d< real64 > const & dPhaseCompFraction_dPressure, + arraySlice2d< real64 > const & dPhaseCompFraction_dTemperature, + arraySlice3d< real64 > const & dPhaseCompFraction_dGlobalCompFraction, + real64 & totalDensity, + real64 & dTotalDensity_dPressure, + real64 & dTotalDensity_dTemperature, + arraySlice1d< real64, 0 > const & dTotalDensity_dGlobalCompFraction ) const +{ + GEOSX_UNUSED_VAR( temperature ); + GEOSX_UNUSED_VAR( dPhaseFraction_dTemperature ); + GEOSX_UNUSED_VAR( dPhaseDensity_dTemperature ); + GEOSX_UNUSED_VAR( dPhaseMassDensity_dTemperature ); + GEOSX_UNUSED_VAR( dPhaseViscosity_dTemperature ); + GEOSX_UNUSED_VAR( dPhaseCompFraction_dTemperature ); + GEOSX_UNUSED_VAR( dTotalDensity_dTemperature ); + + localIndex const NC = numComponents(); + localIndex const NP = numPhases(); + + // 1. Read viscosities and formation volume factors from tables, update mass densities + computeViscosities( pressure, + phaseViscosity, + dPhaseViscosity_dPressure, + dPhaseViscosity_dGlobalCompFraction ); + computeDensities( pressure, + phaseMassDensity, + dPhaseMassDensity_dPressure, + dPhaseMassDensity_dGlobalCompFraction ); + + // 2. Update phaseDens (mass density if useMass == 1, molar density otherwise) + for( localIndex ip = 0; ip < NP; ++ip ) + { + real64 const mult = m_useMass ? 1.0 : 1.0 / m_componentMolarWeight[ip]; + phaseDensity[ip] = phaseMassDensity[ip] * mult; + dPhaseDensity_dPressure[ip] = dPhaseMassDensity_dPressure[ip] * mult; + for( localIndex ic = 0; ic < NC; ++ic ) + { + dPhaseDensity_dGlobalCompFraction[ip][ic] = 0.0; + } + } + + // 3. Update remaining variables: phaseFrac, phaseCompFrac using Dead-Oil assumptions + for( localIndex ip = 0; ip < NP; ++ip ) + { + phaseFraction[ip] = composition[ip]; + dPhaseFraction_dPressure[ip] = 0.0; + for( localIndex ic = 0; ic < NC; ++ic ) + { + dPhaseFraction_dGlobalCompFraction[ip][ic] = (ip == ic) ? 1-composition[ip] : -composition[ip]; + + phaseCompFraction[ip][ic] = (ip == ic) ? 1.0 : 0.0; + dPhaseCompFraction_dPressure[ip][ic] = 0.0; + for( localIndex jc = 0; jc < NC; ++jc ) + { + dPhaseCompFraction_dGlobalCompFraction[ip][ic][jc] = 0.0; + } + } + } + + // 4. Compute total fluid mass/molar density and derivatives + totalDensity = 0.0; + dTotalDensity_dPressure = 0.0; + for( localIndex ic = 0; ic < NC; ++ic ) + { + dTotalDensity_dGlobalCompFraction[ic] = 0.0; + } + + // 4.1. Sum mass/molar fraction/density ratio over all phases to get the inverse of density + for( localIndex ip = 0; ip < NP; ++ip ) + { + real64 const densInv = 1.0 / phaseDensity[ip]; + real64 const value = phaseFraction[ip] * densInv; + + totalDensity += value; + dTotalDensity_dPressure += ( dPhaseFraction_dPressure[ip] - value * dPhaseDensity_dPressure[ip] ) * densInv; + for( localIndex ic = 0; ic < NC; ++ic ) + { + dTotalDensity_dGlobalCompFraction[ic] += ( dPhaseFraction_dGlobalCompFraction[ip][ic] + - value * dPhaseDensity_dGlobalCompFraction[ip][ic] ) * densInv; + } + } + + // 4.2. Invert the previous quantity to get actual density + totalDensity = 1.0 / totalDensity; + real64 const minusDens2 = -totalDensity * totalDensity; + dTotalDensity_dPressure *= minusDens2; + for( localIndex ic = 0; ic < NC; ++ic ) + { + dTotalDensity_dGlobalCompFraction[ic] *= minusDens2; + } +} + + +} //namespace constitutive + +} //namespace geosx + +#endif //GEOSX_CONSTITUTIVE_FLUID_DEADOILFLUID_HPP_ diff --git a/src/coreComponents/constitutive/fluid/multiFluidSelector.hpp b/src/coreComponents/constitutive/fluid/multiFluidSelector.hpp index 7b8ee2db83b..b96302b2571 100644 --- a/src/coreComponents/constitutive/fluid/multiFluidSelector.hpp +++ b/src/coreComponents/constitutive/fluid/multiFluidSelector.hpp @@ -19,6 +19,7 @@ #define GEOSX_CONSTITUTIVE_FLUID_MULTIFLUIDSELECTOR_HPP_ #include "constitutive/ConstitutivePassThruHandler.hpp" +#include "constitutive/fluid/DeadOilFluid.hpp" #include "constitutive/fluid/CompositionalMultiphaseFluid.hpp" #include "constitutive/fluid/BlackOilFluid.hpp" #include "constitutive/fluid/MultiPhaseMultiComponentFluid.hpp" @@ -33,7 +34,8 @@ template< typename LAMBDA > void constitutiveUpdatePassThru( MultiFluidBase const & fluid, LAMBDA && lambda ) { - ConstitutivePassThruHandler< BlackOilFluid, + ConstitutivePassThruHandler< DeadOilFluid, + BlackOilFluid, CompositionalMultiphaseFluid, MultiPhaseMultiComponentFluid >::execute( fluid, std::forward< LAMBDA >( lambda ) ); } @@ -42,7 +44,8 @@ template< typename LAMBDA > void constitutiveUpdatePassThru( MultiFluidBase & fluid, LAMBDA && lambda ) { - ConstitutivePassThruHandler< BlackOilFluid, + ConstitutivePassThruHandler< DeadOilFluid, + BlackOilFluid, CompositionalMultiphaseFluid, MultiPhaseMultiComponentFluid >::execute( fluid, std::forward< LAMBDA >( lambda ) ); } diff --git a/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp b/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp index 365be507df1..ae7472cf811 100644 --- a/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp +++ b/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp @@ -422,6 +422,117 @@ MultiFluidBase * makeDeadOilFluid( string const & name, Group * parent ) return fluid; } +MultiFluidBase * makeNativeDeadOilFluid( string const & name, Group * parent ) +{ + FunctionManager * functionManager = &FunctionManager::FunctionManager::instance(); + + // 1) First, define the tables (PVDO, PVDG) + + // 1D table with linear interpolation + localIndex const NaxisPVDO = 7; + localIndex 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] = 2000000; valuesPVDO_Bo[0] = 1.02; valuesPVDO_visc[0] = 0.000975; + coordinatesPVDO[0][1] = 5000000; valuesPVDO_Bo[1] = 1.03; valuesPVDO_visc[1] = 0.00091; + coordinatesPVDO[0][2] = 10000000; valuesPVDO_Bo[2] = 1.04; valuesPVDO_visc[2] = 0.00083; + coordinatesPVDO[0][3] = 20000000; valuesPVDO_Bo[3] = 1.05; valuesPVDO_visc[3] = 0.000695; + coordinatesPVDO[0][4] = 30000000; valuesPVDO_Bo[4] = 1.07; valuesPVDO_visc[4] = 0.000594; + coordinatesPVDO[0][5] = 40000000; valuesPVDO_Bo[5] = 1.08; valuesPVDO_visc[5] = 0.00051; + coordinatesPVDO[0][6] = 50000000; valuesPVDO_Bo[6] = 1.09; valuesPVDO_visc[6] = 0.000449; + + 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 = functionManager->createChild( "TableFunction", "PVDO_Bo" )->groupCast< TableFunction * >(); + tablePVDO_Bo->setTableCoordinates( coordinatesPVDO ); + tablePVDO_Bo->setTableValues( valuesPVDO_Bo ); + tablePVDO_Bo->reInitializeFunction(); + + tablePVDO_Bo->setInterpolationMethod( TableFunction::InterpolationType::Linear ); + + TableFunction * tablePVDO_visc = functionManager->createChild( "TableFunction", "PVDO_visc" )->groupCast< TableFunction * >(); + tablePVDO_visc->setTableCoordinates( coordinatesPVDO ); + tablePVDO_visc->setTableValues( valuesPVDO_visc ); + tablePVDO_visc->reInitializeFunction(); + + tablePVDO_visc->setInterpolationMethod( TableFunction::InterpolationType::Linear ); + + TableFunction * tablePVDG_Bg = functionManager->createChild( "TableFunction", "PVDG_Bg" )->groupCast< TableFunction * >(); + tablePVDG_Bg->setTableCoordinates( coordinatesPVDG ); + tablePVDG_Bg->setTableValues( valuesPVDG_Bg ); + tablePVDG_Bg->reInitializeFunction(); + + tablePVDG_Bg->setInterpolationMethod( TableFunction::InterpolationType::Linear ); + + TableFunction * tablePVDG_visc = functionManager->createChild( "TableFunction", "PVDG_visc" )->groupCast< TableFunction * >(); + tablePVDG_visc->setTableCoordinates( coordinatesPVDG ); + tablePVDG_visc->setTableValues( valuesPVDG_visc ); + tablePVDG_visc->reInitializeFunction(); + + tablePVDG_visc->setInterpolationMethod( TableFunction::InterpolationType::Linear ); + + // 2) Then, define the Dead-Oil constitutive model + + auto fluid = parent->registerGroup< DeadOilFluid >( name ); + + auto & compNames = fluid->getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString ); + compNames.resize( 3 ); + compNames[0] = "oil"; compNames[1] = "gas"; compNames[2] = "water"; + + auto & molarWgt = fluid->getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString ); + molarWgt.resize( 3 ); + molarWgt[0] = 114e-3; molarWgt[1] = 16e-3; molarWgt[2] = 18e-3; + + auto & phaseNames = fluid->getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString ); + phaseNames.resize( 3 ); + phaseNames[0] = "oil"; phaseNames[1] = "gas"; phaseNames[2] = "water"; + + auto & surfaceDens = fluid->getReference< array1d< real64 > >( DeadOilFluid::viewKeyStruct::surfacePhaseMassDensitiesString ); + surfaceDens.resize( 3 ); + surfaceDens[0] = 800.0; surfaceDens[1] = 0.9907; surfaceDens[2] = 1022.0; + + auto & FVFTableNames = fluid->getReference< string_array >( DeadOilFluid::viewKeyStruct::formationVolumeFactorTableNamesString ); + FVFTableNames.resize( 2 ); + FVFTableNames[0] = "PVDO_Bo"; FVFTableNames[1] = "PVDG_Bg"; + + auto & viscosityTableNames = fluid->getReference< string_array >( DeadOilFluid::viewKeyStruct::viscosityTableNamesString ); + viscosityTableNames.resize( 2 ); + viscosityTableNames[0] = "PVDO_visc"; viscosityTableNames[1] = "PVDG_visc"; + + auto & waterRefPressure = fluid->getReference< real64 >( DeadOilFluid::viewKeyStruct::waterRefPressureString ); + waterRefPressure = 30600000.1; + auto & waterFormationVolumeFactor = fluid->getReference< real64 >( DeadOilFluid::viewKeyStruct::waterFormationVolumeFactorString ); + waterFormationVolumeFactor = 1.03; + auto & waterCompressibility = fluid->getReference< real64 >( DeadOilFluid::viewKeyStruct::waterCompressibilityString ); + waterCompressibility = 0.00000000041; + auto & 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 ); @@ -560,6 +671,66 @@ TEST_F( DeadOilFluidTest, numericalDerivativesMass ) testNumericalDerivatives( *fluid, P, T, comp, eps, relTol, absTol ); } +class NativeDeadOilFluidTest : public ::testing::Test +{ +protected: + + virtual void SetUp() override + { + parent = std::make_unique< Group >( "parent", nullptr ); + parent->resize( 1 ); + fluid = makeNativeDeadOilFluid( "fluid", parent.get()); + + parent->initialize( parent.get() ); + parent->initializePostInitialConditions( parent.get() ); + } + + std::unique_ptr< Group > parent; + MultiFluidBase * fluid; +}; + +TEST_F( NativeDeadOilFluidTest, numericalDerivativesMolar ) +{ + fluid->setMassFlag( false ); + + // TODO test over a range of values + real64 const P1 = 5.4e6; + real64 const P2 = 1.24e7; + real64 const P3 = 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; + + testNumericalDerivatives( *fluid, P1, T, comp, eps, relTol ); + testNumericalDerivatives( *fluid, P2, T, comp, eps, relTol ); + testNumericalDerivatives( *fluid, P3, T, comp, eps, relTol ); +} + +TEST_F( NativeDeadOilFluidTest, numericalDerivativesMass ) +{ + fluid->setMassFlag( true ); + + // TODO test over a range of values + real64 const P1 = 5.4e6; + real64 const P2 = 1.24e7; + real64 const P3 = 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; + + testNumericalDerivatives( *fluid, P1, T, comp, eps, relTol, absTol ); + testNumericalDerivatives( *fluid, P2, T, comp, eps, relTol, absTol ); + testNumericalDerivatives( *fluid, P3, T, comp, eps, relTol, absTol ); +} + + int main( int argc, char * * argv ) { ::testing::InitGoogleTest( &argc, argv ); diff --git a/src/coreComponents/fileIO/schema/docs/Constitutive.rst b/src/coreComponents/fileIO/schema/docs/Constitutive.rst index 06ec905a795..2dc9672deef 100644 --- a/src/coreComponents/fileIO/schema/docs/Constitutive.rst +++ b/src/coreComponents/fileIO/schema/docs/Constitutive.rst @@ -14,6 +14,7 @@ Coulomb node :ref:`XML_Coulomb` DamageLinearElasticIsotropic node :ref:`XML_DamageLinearElasticIsotropic` DamageSpectralLinearElasticIsotropic node :ref:`XML_DamageSpectralLinearElasticIsotropic` DamageVolDevLinearElasticIsotropic node :ref:`XML_DamageVolDevLinearElasticIsotropic` +DeadOilFluid node :ref:`XML_DeadOilFluid` LinearElasticAnisotropic node :ref:`XML_LinearElasticAnisotropic` LinearElasticIsotropic node :ref:`XML_LinearElasticIsotropic` LinearElasticTransverseIsotropic node :ref:`XML_LinearElasticTransverseIsotropic` diff --git a/src/coreComponents/fileIO/schema/docs/Constitutive_other.rst b/src/coreComponents/fileIO/schema/docs/Constitutive_other.rst index 9735b878313..993b1bc54f1 100644 --- a/src/coreComponents/fileIO/schema/docs/Constitutive_other.rst +++ b/src/coreComponents/fileIO/schema/docs/Constitutive_other.rst @@ -14,6 +14,7 @@ Coulomb node :ref:`DATASTRUCTURE_Coulomb` DamageLinearElasticIsotropic node :ref:`DATASTRUCTURE_DamageLinearElasticIsotropic` DamageSpectralLinearElasticIsotropic node :ref:`DATASTRUCTURE_DamageSpectralLinearElasticIsotropic` DamageVolDevLinearElasticIsotropic node :ref:`DATASTRUCTURE_DamageVolDevLinearElasticIsotropic` +DeadOilFluid node :ref:`DATASTRUCTURE_DeadOilFluid` LinearElasticAnisotropic node :ref:`DATASTRUCTURE_LinearElasticAnisotropic` LinearElasticIsotropic node :ref:`DATASTRUCTURE_LinearElasticIsotropic` LinearElasticTransverseIsotropic node :ref:`DATASTRUCTURE_LinearElasticTransverseIsotropic` diff --git a/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst b/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst new file mode 100644 index 00000000000..1bbdce1fd8c --- /dev/null +++ b/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst @@ -0,0 +1,21 @@ + + +======================================= ============ ======== =============================================================================================================================================================================================================================== +Name Type Default Description +======================================= ============ ======== =============================================================================================================================================================================================================================== +componentMolarWeight real64_array required Component molar weights +componentNames string_array {} List of component names +hydrocarbonFormationVolFactorTableNames string_array required | Names of the formation volume factor tables (one per hydrocarbon phase, in the order provided in "phaseNames"). + | For instance, if "oil" is before "gas" in "phaseNames", the table order should be: oilTableName, gasTableName +hydrocarbonViscosityTableNames string_array required | Names of the viscosity tables (one per hydrocarbon phase, in the order provided in "phaseNames") + | For instance, if "oil" is before "gas" in "phaseNames", the table order should be: oilTableName, gasTableName +name string required A name is required for any non-unique nodes +phaseNames string_array required List of fluid phases +surfacePhaseMassDensities real64_array required List of surface mass densities for each phase +waterCompressibility real64 0 Water compressibility +waterFormationVolumeFactor real64 0 Water formation volume factor +waterReferencePressure real64 0 Water reference pressure +waterViscosity real64 0 Water viscosity +======================================= ============ ======== =============================================================================================================================================================================================================================== + + diff --git a/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst b/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst new file mode 100644 index 00000000000..5b3d3c27bb3 --- /dev/null +++ b/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst @@ -0,0 +1,33 @@ + + +====================================== ============================================================================================== ========================== +Name Type Description +====================================== ============================================================================================== ========================== +dPhaseCompFraction_dGlobalCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l, 3l, 4l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseCompFraction_dPressure LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseCompFraction_dTemperature LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseDensity_dPressure real64_array3d (no description available) +dPhaseDensity_dTemperature real64_array3d (no description available) +dPhaseFraction_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseFraction_dPressure real64_array3d (no description available) +dPhaseFraction_dTemperature real64_array3d (no description available) +dPhaseMassDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseMassDensity_dPressure real64_array3d (no description available) +dPhaseMassDensity_dTemperature real64_array3d (no description available) +dPhaseViscosity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseViscosity_dPressure real64_array3d (no description available) +dPhaseViscosity_dTemperature real64_array3d (no description available) +dTotalDensity_dGlobalCompFraction real64_array3d (no description available) +dTotalDensity_dPressure real64_array2d (no description available) +dTotalDensity_dTemperature real64_array2d (no description available) +phaseCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +phaseDensity real64_array3d (no description available) +phaseFraction real64_array3d (no description available) +phaseMassDensity real64_array3d (no description available) +phaseViscosity real64_array3d (no description available) +totalDensity real64_array2d (no description available) +useMass integer (no description available) +====================================== ============================================================================================== ========================== + + diff --git a/src/coreComponents/fileIO/schema/schema.xsd b/src/coreComponents/fileIO/schema/schema.xsd index 20fd04f5705..a4e973e0148 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd +++ b/src/coreComponents/fileIO/schema/schema.xsd @@ -1543,6 +1543,7 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> + @@ -1766,6 +1767,32 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/fileIO/schema/schema.xsd.other b/src/coreComponents/fileIO/schema/schema.xsd.other index 4bad1169808..24646392509 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd.other +++ b/src/coreComponents/fileIO/schema/schema.xsd.other @@ -673,6 +673,7 @@ + @@ -890,6 +891,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/managers/Functions/TableFunction.cpp b/src/coreComponents/managers/Functions/TableFunction.cpp index 0f623ca0f39..afc179583f7 100644 --- a/src/coreComponents/managers/Functions/TableFunction.cpp +++ b/src/coreComponents/managers/Functions/TableFunction.cpp @@ -114,6 +114,32 @@ void TableFunction::parseFile( array1d< T > & target, string const & filename, c inputStream.close(); } +void TableFunction::setInterpolationMethod( InterpolationType const method ) +{ + m_interpolationMethod = method; + reInitializeFunction(); +} + +void TableFunction::setTableCoordinates( array1d< real64_array > coordinates ) +{ + m_coordinates.resize( 0 ); + for( localIndex i = 0; i < coordinates.size(); ++i ) + { + for( localIndex j = 1; j < coordinates[i].size(); ++j ) + { + GEOSX_ERROR_IF( coordinates[i][j] - coordinates[i][j-1] <= 0, + "In the table, the coordinates must be strictly increasing, but axis " << i << "is not" ); + } + m_coordinates.appendArray( coordinates[i].begin(), coordinates[i].end() ); + } + reInitializeFunction(); +} + +void TableFunction::setTableValues( real64_array values ) +{ + m_values = std::move( values ); + reInitializeFunction(); +} void TableFunction::initializeFunction() { @@ -127,7 +153,7 @@ void TableFunction::initializeFunction() { // 1D Table m_dimensions = 1; - m_coordinates.emplace_back( m_tableCoordinates1D ); + m_coordinates.appendArray( m_tableCoordinates1D.begin(), m_tableCoordinates1D.end() ); m_size.emplace_back( m_tableCoordinates1D.size()); // Check to make sure that the table dimensions match @@ -137,12 +163,14 @@ void TableFunction::initializeFunction() { // ND Table m_dimensions = LvArray::integerConversion< localIndex >( m_coordinateFiles.size()); - m_coordinates.resize( m_dimensions ); parseFile( m_values, m_voxelFile, ',' ); + array1d< real64 > tmp; for( localIndex ii=0; ii 0 && m_values.size() > 0 ) // coordinates and values have been set + { + GEOSX_ERROR_IF( increment != m_values.size(), "Table dimensions do not match!" ); + } // Build a quick map to help with linear interpolation m_numCorners = static_cast< localIndex >(pow( 2, m_dimensions )); @@ -177,129 +207,90 @@ void TableFunction::reInitializeFunction() m_corners[jj][ii] = int(ii / pow( 2, jj )) % 2; } } + + // Create the kernel wrapper + m_kernelWrapper.create( m_interpolationMethod, + m_coordinates.toViewConst(), + m_values.toViewConst(), + m_dimensions, + m_size.toViewConst(), + m_indexIncrement.toViewConst(), + m_corners, + m_numCorners ); + } +TableFunction::KernelWrapper TableFunction::createKernelWrapper() const +{ + return TableFunction::KernelWrapper( m_interpolationMethod, + m_coordinates.toViewConst(), + m_values.toViewConst(), + m_dimensions, + m_size.toViewConst(), + m_indexIncrement.toViewConst(), + m_corners, + m_numCorners ); +} real64 TableFunction::evaluate( real64 const * const input ) const { - real64 result = 0.0; - - // Linear interpolation - if( m_interpolationMethod == InterpolationType::Linear ) - { - localIndex bounds[m_maxDimensions][2]; - real64 weights[m_maxDimensions][2]; - - // Determine position, weights - for( localIndex ii=0; ii= m_coordinates[ii][m_size[ii] - 1] ) - { - // Coordinate is to the right of this axis - bounds[ii][0] = m_size[ii] - 1; - bounds[ii][1] = bounds[ii][0]; - weights[ii][0] = 1; - weights[ii][1] = 0; - } - else - { - // Find the coordinate index - ///TODO make this fast - // Note: lower_bound uses a binary search... If we assume coordinates are - // evenly spaced, we can speed things up considerably - auto lower = std::lower_bound( m_coordinates[ii].begin(), m_coordinates[ii].end(), input[ii] ); - bounds[ii][1] = LvArray::integerConversion< localIndex >( std::distance( m_coordinates[ii].begin(), lower )); - bounds[ii][0] = bounds[ii][1] - 1; - - real64 dx = m_coordinates[ii][bounds[ii][1]] - m_coordinates[ii][bounds[ii][0]]; - weights[ii][0] = 1.0 - (input[ii] - m_coordinates[ii][bounds[ii][0]]) / dx; - weights[ii][1] = 1.0 - weights[ii][0]; - } - } + real64 scalarValue = 0; + stackArray1d< real64, maxDimensions > derivativesArray( m_dimensions ); - // Calculate the result - for( localIndex ii=0; ii= m_coordinates[ii][m_size[ii] - 1] ) - { - // Coordinate is to the right of the table axis - subIndex = m_size[ii] - 1; - } - else - { - // Coordinate is within the table axis - // Note: std::distance will return the index of the upper table vertex - auto lower = std::lower_bound( m_coordinates[ii].begin(), m_coordinates[ii].end(), input[ii] ); - subIndex = LvArray::integerConversion< localIndex >( std::distance( m_coordinates[ii].begin(), lower )); - - // Interpolation types: - // - Nearest returns the value of the closest table vertex - // - Upper returns the value of the next table vertex - // - Lower returns the value of the previous table vertex - if( m_interpolationMethod == InterpolationType::Nearest ) - { - if((input[ii] - m_coordinates[ii][subIndex - 1]) <= (m_coordinates[ii][subIndex] - input[ii])) - { - --subIndex; - } - } - else if( m_interpolationMethod == InterpolationType::Lower ) - { - if( subIndex > 0 ) - { - --subIndex; - } - } - } + // interpolate in table, return scalar value (derivatives are discarded) + m_kernelWrapper.compute( input, scalarValue, derivativesArray ); + return scalarValue; +} - // Increment the global table index - tableIndex += subIndex * m_indexIncrement[ii]; - } +TableFunction::KernelWrapper::KernelWrapper( TableFunction::InterpolationType interpolationMethod, + ArrayOfArraysView< real64 const > const & coordinates, + arrayView1d< real64 const > const & values, + localIndex dimensions, + arrayView1d< localIndex const > const & size, + arrayView1d< localIndex const > const & indexIncrement, + localIndex const (&corners)[TableFunction::maxDimensions][16], + localIndex const numCorners ) + : + m_interpolationMethod( interpolationMethod ), + m_coordinates( coordinates ), + m_values( values ), + m_dimensions( dimensions ), + m_size( size ), + m_indexIncrement( indexIncrement ), + m_numCorners( numCorners ) +{ + LvArray::tensorOps::copy< TableFunction::maxDimensions, 16 >( m_corners, corners ); +} - // Retrieve the nearest value - result = m_values[tableIndex]; - } +TableFunction::KernelWrapper::KernelWrapper() + : + m_interpolationMethod(), + m_coordinates(), + m_values(), + m_dimensions( 0 ), + m_size(), + m_indexIncrement(), + m_numCorners( 0 ) +{} - return result; +void TableFunction::KernelWrapper::create( TableFunction::InterpolationType interpolationMethod, + ArrayOfArraysView< real64 const > const & coordinates, + arrayView1d< real64 const > const & values, + localIndex dimensions, + arrayView1d< localIndex const > const & size, + arrayView1d< localIndex const > const & indexIncrement, + localIndex const (&corners)[TableFunction::maxDimensions][16], + localIndex const numCorners ) +{ + m_interpolationMethod = interpolationMethod; + m_coordinates = coordinates; + m_values = values; + m_dimensions = dimensions; + m_size = size; + m_indexIncrement = indexIncrement; + m_numCorners = numCorners; + + LvArray::tensorOps::copy< TableFunction::maxDimensions, 16 >( m_corners, corners ); } REGISTER_CATALOG_ENTRY( FunctionBase, TableFunction, string const &, Group * const ) diff --git a/src/coreComponents/managers/Functions/TableFunction.hpp b/src/coreComponents/managers/Functions/TableFunction.hpp index 62652c2d788..1ea89ad91b2 100644 --- a/src/coreComponents/managers/Functions/TableFunction.hpp +++ b/src/coreComponents/managers/Functions/TableFunction.hpp @@ -21,6 +21,7 @@ #include "common/EnumStrings.hpp" #include "managers/Functions/FunctionBase.hpp" +#include "LvArray/src/tensorOps.hpp" namespace geosx { @@ -33,6 +34,125 @@ namespace geosx class TableFunction : public FunctionBase { public: + + /// Enumerator of available interpolation types + enum class InterpolationType : integer + { + Linear, + Nearest, + Upper, + Lower + }; + + /// maximum dimensions for the coordinates in the table + static constexpr localIndex maxDimensions = 4; + + /** + * @class KernelWrapper + * + * A nested class encapsulating the kernel function doing the interpolation in the table + */ + class KernelWrapper + { +public: + + /** + * @brief The constructor of the kernel wrapper + * @param[in] interpolationMethod table interpolation method + * @param[in] coordinates array of table axes + * @param[in] values table values (in fortran order) + * @param[in] dimensions number of active table dimensions + * @param[in] size size of the table + * @param[in] indexIncrement array used to locate values within ND tables + * @param[in] corners corners of the box that surround the value in N dimensions + * @param[in] numCorners number of active table corners + */ + KernelWrapper( TableFunction::InterpolationType interpolationMethod, + ArrayOfArraysView< real64 const > const & coordinates, + arrayView1d< real64 const > const & values, + localIndex dimensions, + arrayView1d< localIndex const > const & size, + arrayView1d< localIndex const > const & indexIncrement, + localIndex const (&corners)[TableFunction::maxDimensions][16], + localIndex const numCorners ); + + /// Default constructor for the kernel wrapper + KernelWrapper(); + + /// Default copy constructor + KernelWrapper( KernelWrapper const & ) = default; + + /// Default move constructor + KernelWrapper( KernelWrapper && ) = default; + + /// Copy assignment operator + KernelWrapper & operator=( KernelWrapper const & ) = delete; + + /// Deleted move assignment operator + KernelWrapper & operator=( KernelWrapper && ) = delete; + + /** + * @brief The function that populates the wrapper + * @param[in] interpolationMethod table interpolation method + * @param[in] coordinates array of table axes + * @param[in] values table values (in fortran order) + * @param[in] dimensions number of active table dimensions + * @param[in] size size of the table + * @param[in] indexIncrement array used to locate values within ND tables + * @param[in] corners corners of the box that surround the value in N dimensions + * @param[in] numCorners number of active table corners + */ + void create( TableFunction::InterpolationType interpolationMethod, + ArrayOfArraysView< real64 const > const & coordinates, + arrayView1d< real64 const > const & values, + localIndex dimensions, + arrayView1d< localIndex const > const & size, + arrayView1d< localIndex const > const & indexIncrement, + localIndex const (&corners)[TableFunction::maxDimensions][16], + localIndex const numCorners ); + + /** + * @brief Main compute function to interpolate in the table and return the derivatives + * @param[in] input vector of input value + * @param[out] value interpolated value + * @param[out] derivatives vector of derivatives of interpolated value wrt the variables present in input + */ + template< typename IN_ARRAY, typename OUT_ARRAY > + GEOSX_HOST_DEVICE + void + compute( IN_ARRAY const & input, real64 & value, OUT_ARRAY && derivatives ) const; + +private: + + /// Table interpolation method + TableFunction::InterpolationType m_interpolationMethod; + + /// An array of table axes + ArrayOfArraysView< real64 const > m_coordinates; + + /// Table values (in fortran order) + arrayView1d< real64 const > m_values; + + /// Number of active table dimensions + localIndex m_dimensions; + + /// Size of the table + arrayView1d< localIndex const > m_size; + + /// Array used to locate values within ND tables + arrayView1d< localIndex const > m_indexIncrement; + + /** + * @brief The corners of the box that surround the value in N dimensions + * m_corners should be of size m_maxDimensions x (2^m_maxDimensions) + */ + localIndex m_corners[TableFunction::maxDimensions][16]; + + /// The number of active table corners + localIndex m_numCorners; + + }; + /** * @brief The constructor * @param[in] name the name of this object manager @@ -99,12 +219,12 @@ class TableFunction : public FunctionBase * @brief Get the table axes definitions * @return a reference to an array of arrays that define each table axis */ - array1d< real64_array > const & getCoordinates() const { return m_coordinates; } + ArrayOfArraysView< real64 const > getCoordinates() const { return m_coordinates.toViewConst(); } /** * @copydoc getCoordinates() const */ - array1d< real64_array > & getCoordinates() { return m_coordinates; } + ArrayOfArraysView< real64 > getCoordinates() { return m_coordinates.toView(); } /** * @brief Get the table values @@ -117,32 +237,35 @@ class TableFunction : public FunctionBase */ array1d< real64 > & getValues() { return m_values; } - /// Enumerator of available interpolation types - enum class InterpolationType : integer - { - Linear, - Nearest, - Upper, - Lower - }; + /** + * @brief Get the interpolation method + * @return The interpolation method + */ + InterpolationType getInterpolationMethod() const { return m_interpolationMethod; } /** * @brief Set the interpolation method * @param method The interpolation method */ - void setInterpolationMethod( InterpolationType const method ) { m_interpolationMethod = method; } + void setInterpolationMethod( InterpolationType const method ); /** * @brief Set the table coordinates * @param coordinates An array of arrays containing table coordinate definitions */ - void setTableCoordinates( array1d< real64_array > coordinates ) { m_coordinates = coordinates; } + void setTableCoordinates( array1d< real64_array > coordinates ); /** * @brief Set the table values * @param values An array of table values in fortran order */ - void setTableValues( real64_array values ) { m_values = values; } + void setTableValues( real64_array values ); + + /** + * @brief Create an instance of the kernel wrapper + * @return the kernel wrapper + */ + KernelWrapper createKernelWrapper() const; private: /// Coordinates for 1D table @@ -158,14 +281,11 @@ class TableFunction : public FunctionBase InterpolationType m_interpolationMethod; /// An array of table axes - array1d< real64_array > m_coordinates; + ArrayOfArrays< real64 > m_coordinates; /// Table values (in fortran order) real64_array m_values; - /// Maximum number of table dimensions - static localIndex constexpr m_maxDimensions = 4; - /// Number of active table dimensions localIndex m_dimensions; @@ -177,16 +297,170 @@ class TableFunction : public FunctionBase /** * @brief The corners of the box that surround the value in N dimensions - * m_corners should be of size m_maxDimensions x (2^m_maxDimensions) + * m_corners should be of size maxDimensions x (2^maxDimensions) */ - localIndex m_corners[m_maxDimensions][16]; + localIndex m_corners[maxDimensions][16]; /// The number of active table corners localIndex m_numCorners; + + /// Kernel wrapper to interpolate in table and return derivatives + KernelWrapper m_kernelWrapper; + }; -ENUM_STRINGS( TableFunction::InterpolationType, "linear", "nearest", "upper", "lower" ) +template< typename IN_ARRAY, typename OUT_ARRAY > +GEOSX_HOST_DEVICE +void +TableFunction::KernelWrapper::compute( IN_ARRAY const & input, real64 & value, OUT_ARRAY && derivatives ) const +{ + value = 0.0; + for( localIndex i = 0; i < m_dimensions; ++i ) + { + derivatives[i] = 0.0; + } + // Linear interpolation + if( m_interpolationMethod == TableFunction::InterpolationType::Linear ) + { + localIndex bounds[TableFunction::maxDimensions][2]{}; + real64 weights[TableFunction::maxDimensions][2]{}; + real64 dWeights_dInput[TableFunction::maxDimensions][2]{}; + real64 dCornerValue_dInput[TableFunction::maxDimensions]{}; + + // Determine position, weights + for( localIndex ii=0; ii= m_coordinates[ii][m_size[ii] - 1] ) + { + // Coordinate is to the right of this axis + bounds[ii][0] = m_size[ii] - 1; + bounds[ii][1] = bounds[ii][0]; + weights[ii][0] = 1; + weights[ii][1] = 0; + dWeights_dInput[ii][0] = 0; + dWeights_dInput[ii][1] = 0; + } + else + { + // Find the coordinate index + ///TODO make this fast + // Note: find uses a binary search... If we assume coordinates are + // evenly spaced, we can speed things up considerably + auto lower = LvArray::sortedArrayManipulation::find( m_coordinates[ii].begin(), m_coordinates.sizeOfArray( ii ), input[ii] ); + bounds[ii][1] = LvArray::integerConversion< localIndex >( lower ); + bounds[ii][0] = bounds[ii][1] - 1; + + real64 dx = m_coordinates[ii][bounds[ii][1]] - m_coordinates[ii][bounds[ii][0]]; + weights[ii][0] = 1.0 - (input[ii] - m_coordinates[ii][bounds[ii][0]]) / dx; + weights[ii][1] = 1.0 - weights[ii][0]; + dWeights_dInput[ii][0] = -1.0 / dx; + dWeights_dInput[ii][1] = -dWeights_dInput[ii][0]; + } + } + + // Calculate the result + for( localIndex ii=0; ii= m_coordinates[ii][m_size[ii] - 1] ) + { + // Coordinate is to the right of the table axis + subIndex = m_size[ii] - 1; + } + else + { + // Coordinate is within the table axis + // Note: std::distance will return the index of the upper table vertex + auto lower = LvArray::sortedArrayManipulation::find( m_coordinates[ii].begin(), m_coordinates.sizeOfArray( ii ), input[ii] ); + subIndex = LvArray::integerConversion< localIndex >( lower ); + + // Interpolation types: + // - Nearest returns the value of the closest table vertex + // - Upper returns the value of the next table vertex + // - Lower returns the value of the previous table vertex + if( m_interpolationMethod == TableFunction::InterpolationType::Nearest ) + { + if((input[ii] - m_coordinates[ii][subIndex - 1]) <= (m_coordinates[ii][subIndex] - input[ii])) + { + --subIndex; + } + } + else if( m_interpolationMethod == TableFunction::InterpolationType::Lower ) + { + if( subIndex > 0 ) + { + --subIndex; + } + } + } + + // Increment the global table index + tableIndex += subIndex * m_indexIncrement[ii]; + } + + // Retrieve the nearest value + value = m_values[tableIndex]; + } +} + +ENUM_STRINGS( TableFunction::InterpolationType, "linear", "nearest", "upper", "lower" ) } /* namespace geosx */ diff --git a/src/coreComponents/managers/unitTests/testFunctions.cpp b/src/coreComponents/managers/unitTests/testFunctions.cpp index 06ba96dada9..c0cb1f45e6d 100644 --- a/src/coreComponents/managers/unitTests/testFunctions.cpp +++ b/src/coreComponents/managers/unitTests/testFunctions.cpp @@ -12,6 +12,7 @@ * ------------------------------------------------------------------------------------------------------------ */ +#include "codingUtilities/UnitTestUtilities.hpp" #include "gtest/gtest.h" #include "managers/initialization.hpp" #include "managers/Functions/FunctionManager.hpp" @@ -26,7 +27,6 @@ using namespace geosx; - void evaluate1DFunction( FunctionBase * function, arrayView1d< real64 const > const & inputs, arrayView1d< real64 const > const & outputs ) @@ -41,15 +41,32 @@ void evaluate1DFunction( FunctionBase * function, } } +void checkDirectionalDerivative( real64 const (&input)[4], + real64 (& perturbedInput)[4], + real64 const & val, + real64 & perturbedVal, + real64 const (&derivatives)[4], + real64 (& perturbedDerivatives)[4], + real64 const & perturb, + real64 const & relTol, + localIndex const direction, + TableFunction::KernelWrapper kernelWrapper ) +{ + LvArray::tensorOps::copy< 4 >( perturbedInput, input ); + real64 const dInput = perturb * ( input[direction] + perturb ); + perturbedInput[direction] += dInput; + kernelWrapper.compute( perturbedInput, perturbedVal, perturbedDerivatives ); + geosx::testing::checkRelativeError( derivatives[direction], (perturbedVal-val)/dInput, relTol, geosx::testing::DEFAULT_ABS_TOL ); +} TEST( FunctionTests, 1DTable ) { FunctionManager * functionManager = &FunctionManager::FunctionManager::instance(); // 1D table, various interpolation methods - localIndex Naxis = 4; - localIndex Ntest = 6; + localIndex const Naxis = 4; + localIndex const Ntest = 6; // Setup table array1d< real64_array > coordinates; @@ -89,6 +106,7 @@ TEST( FunctionTests, 1DTable ) testExpected[4] = 3.0; testExpected[5] = 7.0; table_a->setInterpolationMethod( TableFunction::InterpolationType::Linear ); + table_a->reInitializeFunction(); evaluate1DFunction( table_a, testCoordinates, testExpected ); // Upper @@ -99,6 +117,7 @@ TEST( FunctionTests, 1DTable ) testExpected[4] = 7.0; testExpected[5] = 7.0; table_a->setInterpolationMethod( TableFunction::InterpolationType::Upper ); + table_a->reInitializeFunction(); evaluate1DFunction( table_a, testCoordinates, testExpected ); // Lower @@ -109,6 +128,7 @@ TEST( FunctionTests, 1DTable ) testExpected[4] = -5.0; testExpected[5] = 7.0; table_a->setInterpolationMethod( TableFunction::InterpolationType::Lower ); + table_a->reInitializeFunction(); evaluate1DFunction( table_a, testCoordinates, testExpected ); // Nearest @@ -119,6 +139,7 @@ TEST( FunctionTests, 1DTable ) testExpected[4] = 7.0; testExpected[5] = 7.0; table_a->setInterpolationMethod( TableFunction::InterpolationType::Nearest ); + table_a->reInitializeFunction(); evaluate1DFunction( table_a, testCoordinates, testExpected ); } @@ -131,11 +152,11 @@ TEST( FunctionTests, 2DTable ) // 2D table with linear interpolation // f(x, y) = 2*x - 3*y + 5 - localIndex Ndim = 2; - localIndex Nx = 3; - localIndex Ny = 4; - localIndex Ntest = 20; - string inputName = "coordinates"; + localIndex const Ndim = 2; + localIndex const Nx = 3; + localIndex const Ny = 4; + localIndex const Ntest = 20; + string const inputName = "coordinates"; // Setup table array1d< real64_array > coordinates; @@ -156,8 +177,8 @@ TEST( FunctionTests, 2DTable ) { for( localIndex ii=0; iireInitializeFunction(); // Setup a group for testing the batch mode function evaluation - string groupName = "testGroup"; + string const groupName = "testGroup"; dataRepository::Group testGroup( groupName, nullptr ); real64_array2d testCoordinates; @@ -208,8 +229,8 @@ TEST( FunctionTests, 2DTable ) testCoordinates[ii][jj] = distribution( generator ); } - real64 x = testCoordinates[ii][0]; - real64 y = testCoordinates[ii][1]; + real64 const x = testCoordinates[ii][0]; + real64 const y = testCoordinates[ii][1]; expected[ii] = (2.0*x) - (3.0*y) + 5.0; set.insert( ii ); } @@ -229,17 +250,17 @@ TEST( FunctionTests, 4DTable_multipleInputs ) { FunctionManager * functionManager = &FunctionManager::FunctionManager::instance(); - // 3D table with linear interpolation + // 4D table with linear interpolation // f(x, y, z, t) = 2.0 + 3*x - 5*y + 7*z + 11*t - localIndex Ndim = 4; - localIndex Nx = 3; - localIndex Ny = 4; - localIndex Nz = 5; - localIndex Nt = 2; - localIndex Ntest = 20; - localIndex Ntimes = 5; - string coordinatesName = "coordinates"; - string timeName = "time"; + localIndex const Ndim = 4; + localIndex const Nx = 3; + localIndex const Ny = 4; + localIndex const Nz = 5; + localIndex const Nt = 2; + localIndex const Ntest = 20; + localIndex const Ntimes = 5; + string const coordinatesName = "coordinates"; + string const timeName = "time"; // Setup table array1d< real64_array > coordinates; @@ -273,10 +294,10 @@ TEST( FunctionTests, 4DTable_multipleInputs ) { for( localIndex ii=0; iireInitializeFunction(); // Setup a group for testing the batch mode function evaluation - string groupName = "testGroup"; + string const groupName = "testGroup"; dataRepository::Group testGroup( groupName, nullptr ); real64_array2d testCoordinates; @@ -325,7 +346,7 @@ TEST( FunctionTests, 4DTable_multipleInputs ) // Build the inputs for( localIndex ii=0; ii coordinates; + coordinates.resize( Ndim ); + coordinates[0].resize( Nx ); + coordinates[0][0] = -1.0; + coordinates[0][1] = 0.0; + coordinates[0][2] = 1.0; + coordinates[1].resize( Ny ); + coordinates[1][0] = -1.0; + coordinates[1][1] = 0.0; + coordinates[1][2] = 0.5; + coordinates[1][3] = 1.0; + coordinates[2].resize( Nz ); + coordinates[2][0] = -1.0; + coordinates[2][1] = -0.4; + coordinates[2][2] = 0.3; + coordinates[2][3] = 0.5; + coordinates[2][4] = 1.0; + coordinates[3].resize( Nt ); + coordinates[3][0] = -1.0; + coordinates[3][1] = 0.34; + coordinates[3][2] = 0.5; + coordinates[3][3] = 1.0; + + real64_array values( Nx * Ny * Nz * Nt ); + for( localIndex mm=0, tablePosition=0; mmcreateChild( "TableFunction", "table_d" )->groupCast< TableFunction * >(); + table_d->setTableCoordinates( coordinates ); + table_d->setTableValues( values ); + table_d->setInterpolationMethod( TableFunction::InterpolationType::Linear ); + table_d->setInputVarNames( inputVarNames ); + table_d->reInitializeFunction(); + + real64 const relTol = 1e-4; + real64 const eps = std::numeric_limits< real64 >::epsilon(); + real64 const perturb = std::sqrt( eps ); + + real64 const start = coordinates[0][0]-0.1; // start outside the table, to check derivatives there too + real64 const end = coordinates[0][Nx-1]+0.1; // end outside the table + localIndex const nSamples = 7; // try not to fall on table coordinate, otherwise the finite-difference approximation won't work + real64 const delta = (end-start)/nSamples; + + real64 val = 0.0; + real64 perturbedVal = 0.0; + real64 input[4] = { start, start, start, start }; + real64 perturbedInput[4]{}; + real64 derivatives[4]{}; + real64 perturbedDerivatives[4]{}; + + TableFunction::KernelWrapper kernelWrapper = table_d->createKernelWrapper(); + for( localIndex mm=0; mmcreateChild( "SymbolicFunction", "table_d" )->groupCast< SymbolicFunction * >(); - table_d->setSymbolicExpression( expression ); - table_d->setInputVarNames( inputVarNames ); - table_d->setSymbolicVariableNames( inputVarNames ); - table_d->initializeFunction(); + SymbolicFunction * table_e = functionManager->createChild( "SymbolicFunction", "table_e" )->groupCast< SymbolicFunction * >(); + table_e->setSymbolicExpression( expression ); + table_e->setInputVarNames( inputVarNames ); + table_e->setSymbolicVariableNames( inputVarNames ); + table_e->initializeFunction(); // Setup a group for testing the batch mode function evaluation - string groupName = "testGroup"; + string const groupName = "testGroup"; dataRepository::Group testGroup( groupName, nullptr ); real64_array inputA; real64_array inputB; @@ -412,10 +548,10 @@ TEST( FunctionTests, 4DTable_symbolic ) // Build the inputs for( localIndex ii=0; iievaluate( &(testGroup), 0.0, set.toView(), output ); + table_e->evaluate( &(testGroup), 0.0, set.toView(), output ); // Compare results for( localIndex jj=0; jj Date: Thu, 25 Feb 2021 09:26:48 -0800 Subject: [PATCH 2/6] updated xml files for integratedTests and benchmarks --- .../deadoil_3ph_baker_1d.xml | 3 +- .../constitutive/fluid/BlackOilFluid.cpp | 21 +- .../constitutive/fluid/BlackOilFluid.hpp | 14 +- .../constitutive/fluid/DeadOilFluid.cpp | 188 +++++++++----- .../constitutive/fluid/DeadOilFluid.hpp | 8 +- .../constitutive/unitTests/testMultiFluid.cpp | 236 +++--------------- .../fileIO/schema/docs/BlackOilFluid.rst | 23 +- .../fileIO/schema/docs/DeadOilFluid.rst | 27 +- .../fileIO/schema/docs/DeadOilFluid_other.rst | 64 ++--- src/coreComponents/fileIO/schema/schema.xsd | 25 +- .../fileIO/schema/schema.xsd.other | 12 + .../fluidFlow/benchmarks/Egg/dead_oil_egg.xml | 68 +++-- .../fluidFlow/benchmarks/Egg/pvdg.txt | 15 -- .../fluidFlow/benchmarks/Egg/pvdo.txt | 12 +- .../fluidFlow/benchmarks/Egg/pvtw.txt | 5 +- .../SPE10/dead_oil_spe10_layers_83_84_85.xml | 44 ++-- .../fluidFlow/benchmarks/SPE10/pvdg.txt | 14 -- .../fluidFlow/benchmarks/SPE10/pvdo.txt | 8 +- .../fluidFlow/benchmarks/SPE10/pvtw.txt | 4 +- .../blackoil_2ph_1d.xml | 1 - .../deadoil_3ph_baker_1d.xml | 9 +- .../deadoil_3ph_corey_1d.xml | 9 +- .../deadoil_3ph_staircase_3d.xml | 9 +- ...l_3ph_staircase_gravity_segregation_3d.xml | 29 +-- .../compositionalMultiphaseFlow/pvdg.txt | 28 +-- .../compositionalMultiphaseFlow/pvdo.txt | 16 +- .../compositionalMultiphaseFlow/pvtw.txt | 4 +- .../dead_oil_wells_2d.xml | 3 +- .../compositionalMultiphaseWell/pvdo.txt | 16 +- 29 files changed, 350 insertions(+), 565 deletions(-) delete mode 100644 src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvdg.txt delete mode 100644 src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvdg.txt diff --git a/examples/AdaptiveTimeStepping/deadoil_3ph_baker_1d.xml b/examples/AdaptiveTimeStepping/deadoil_3ph_baker_1d.xml index 160731a9a88..606152b05f4 100644 --- a/examples/AdaptiveTimeStepping/deadoil_3ph_baker_1d.xml +++ b/examples/AdaptiveTimeStepping/deadoil_3ph_baker_1d.xml @@ -89,9 +89,8 @@ - ::concat( "\n* " ) ); } BlackOilFluid::~BlackOilFluid() @@ -97,23 +94,7 @@ void BlackOilFluid::createFluid() std::vector< double > densities( m_surfaceDensities.begin(), m_surfaceDensities.end() ); std::vector< double > molarWeights( m_componentMolarWeight.begin(), m_componentMolarWeight.end() ); - switch( m_fluidType ) - { - case FluidType::LiveOil: - { - m_fluid = pvt::MultiphaseSystemBuilder::buildLiveOil( phases, tableFiles, densities, molarWeights ); - break; - } - case FluidType::DeadOil: - { - m_fluid = pvt::MultiphaseSystemBuilder::buildDeadOil( phases, tableFiles, densities, molarWeights ); - break; - } - default: - { - GEOSX_ERROR( "Unknown fluid type" ); - } - } + m_fluid = pvt::MultiphaseSystemBuilder::buildLiveOil( phases, tableFiles, densities, molarWeights ); } REGISTER_CATALOG_ENTRY( ConstitutiveBase, BlackOilFluid, string const &, Group * const ) diff --git a/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp b/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp index bd51d5d1cac..f0687081200 100644 --- a/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp +++ b/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp @@ -1,5 +1,5 @@ /* - * ------------------------------------------------------------------------------------------------------------ + 1;5202;0c * ------------------------------------------------------------------------------------------------------------ * SPDX-License-Identifier: LGPL-2.1-only * * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC @@ -32,12 +32,6 @@ class BlackOilFluid : public MultiFluidPVTPackageWrapper { public: - enum class FluidType : integer - { - DeadOil, - LiveOil - }; - BlackOilFluid( string const & name, Group * const parent ); virtual ~BlackOilFluid() override; @@ -55,7 +49,6 @@ class BlackOilFluid : public MultiFluidPVTPackageWrapper { static constexpr char const * surfaceDensitiesString() { return "surfaceDensities"; } static constexpr char const * tableFilesString() { return "tableFiles"; } - static constexpr char const * fluidTypeString() { return "fluidType"; } }; protected: @@ -71,13 +64,8 @@ class BlackOilFluid : public MultiFluidPVTPackageWrapper // Black-oil table filenames path_array m_tableFiles; - // Type of black-oil fluid (live/dead) - FluidType m_fluidType; - }; -ENUM_STRINGS( BlackOilFluid::FluidType, "DeadOil", "LiveOil" ) - } /* namespace constitutive */ } /* namespace geosx */ diff --git a/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp b/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp index cc385f935e6..819fbe144e8 100644 --- a/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp @@ -18,10 +18,12 @@ #include "managers/GeosxState.hpp" #include "managers/Functions/FunctionManager.hpp" + namespace geosx { using namespace dataRepository; +using namespace stringutilities; namespace constitutive { @@ -54,36 +56,23 @@ DeadOilFluid::DeadOilFluid( string const & name, getWrapperBase( viewKeyStruct::componentMolarWeightString() ).setInputFlag( InputFlags::REQUIRED ); getWrapperBase( viewKeyStruct::phaseNamesString() ).setInputFlag( InputFlags::REQUIRED ); + registerWrapper( viewKeyStruct::tableFilesString(), &m_tableFiles ). + setInputFlag( InputFlags::REQUIRED ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( "List of filenames with input PVT tables" ); + registerWrapper( viewKeyStruct::surfacePhaseMassDensitiesString(), &m_surfacePhaseMassDensity ). setInputFlag( InputFlags::REQUIRED ). setDescription( "List of surface mass densities for each phase" ); registerWrapper( viewKeyStruct::formationVolumeFactorTableNamesString(), &m_formationVolFactorTableNames ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "Names of the formation volume factor tables (one per hydrocarbon phase, in the order provided in \"phaseNames\"). \n" - "For instance, if \"oil\" is before \"gas\" in \"phaseNames\", the table order should be: oilTableName, gasTableName" ); - + setSizedFromParent( 0 ); registerWrapper( viewKeyStruct::viscosityTableNamesString(), &m_viscosityTableNames ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "Names of the viscosity tables (one per hydrocarbon phase, in the order provided in \"phaseNames\") \n" - "For instance, if \"oil\" is before \"gas\" in \"phaseNames\", the table order should be: oilTableName, gasTableName" ); - - registerWrapper( viewKeyStruct::waterRefPressureString(), &m_waterRefPressure ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Water reference pressure" ); - - registerWrapper( viewKeyStruct::waterFormationVolumeFactorString(), &m_waterFormationVolFactor ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Water formation volume factor" ); - - registerWrapper( viewKeyStruct::waterCompressibilityString(), &m_waterCompressibility ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Water compressibility" ); - - registerWrapper( viewKeyStruct::waterViscosityString(), &m_waterViscosity ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Water viscosity" ); - + setSizedFromParent( 0 );; + registerWrapper( viewKeyStruct::waterRefPressureString(), &m_waterRefPressure ); + registerWrapper( viewKeyStruct::waterFormationVolumeFactorString(), &m_waterFormationVolFactor ); + registerWrapper( viewKeyStruct::waterCompressibilityString(), &m_waterCompressibility ); + registerWrapper( viewKeyStruct::waterViscosityString(), &m_waterViscosity ); } void DeadOilFluid::postProcessInput() @@ -105,44 +94,116 @@ void DeadOilFluid::postProcessInput() m_phaseTypes[ip] = phaseIndex; m_phaseOrder[phaseIndex] = LvArray::integerConversion< integer >( ip ); - if( phaseIndex == PhaseType::OIL || phaseIndex == PhaseType::GAS ) - { - m_hydrocarbonPhaseOrder.emplace_back( LvArray::integerConversion< integer >( ip ) ); - } } - // check the water phase parameters - integer const ipWater = m_phaseOrder[PhaseType::WATER]; - integer const ipGas = m_phaseOrder[PhaseType::GAS]; - if( ipWater >= 0 ) // if water is present - { - GEOSX_ERROR_IF( m_waterRefPressure <= 0.0, - "DeadOilFluid: a strictly positive value must be provided for: " << viewKeyStruct::waterRefPressureString() ); - GEOSX_ERROR_IF( m_waterFormationVolFactor <= 0.0, - "DeadOilFluid: a strictly positive value must be provided for: " << viewKeyStruct::waterFormationVolumeFactorString() ); - GEOSX_ERROR_IF( m_waterViscosity <= 0.0, - "DeadOilFluid: a strictly positive value must be provided for: " << viewKeyStruct::waterViscosityString() ); - GEOSX_ERROR_IF( m_waterCompressibility <= 0.0, - "DeadOilFluid: a strictly positive value must be provided for: " << viewKeyStruct::waterCompressibilityString() ); - } - else + GEOSX_ERROR_IF( m_tableFiles.size() != numFluidPhases(), + "The number of table files must be equal to the number of phases" ); + for( localIndex ip = 0; ip < numFluidPhases(); ++ip ) { - GEOSX_ERROR_IF( m_waterRefPressure > 0.0, - "DeadOilFluid: if water is absent, this keyword is not needed " << viewKeyStruct::waterRefPressureString() ); - GEOSX_ERROR_IF( m_waterFormationVolFactor > 0.0, - "DeadOilFluid: if water is absent, this keyword is not needed " << viewKeyStruct::waterFormationVolumeFactorString() ); - GEOSX_ERROR_IF( m_waterViscosity > 0.0, - "DeadOilFluid: if water is absent, this keyword is not needed " << viewKeyStruct::waterViscosityString() ); - GEOSX_ERROR_IF( m_waterCompressibility > 0.0, - "DeadOilFluid: if water is absent, this keyword is not needed " << viewKeyStruct::waterCompressibilityString() ); - } + string const filename = m_tableFiles[ip]; + + array1d< array1d< real64 > > pressureCoords; + pressureCoords.resize( 1 ); + array1d< real64 > formationVolFactor; + array1d< real64 > viscosity; + + std::ifstream is( filename ); + constexpr std::streamsize buf_size = 256; + char buf[buf_size]; + + localIndex counter = 0; + while( is.getline( buf, buf_size ) ) + { + string const str( buf ); + string_array const strs = Tokenize( str, " " ); + + if( strs.empty() ) + { + continue; + } + + if( strs[0].front() == '#' ) + { + continue; + } + if( streq( strs[0], "--" ) ) + { + continue; + } + + std::cout << "strs.size() = " << strs.size() << std::endl; + for( localIndex i = 0; i < strs.size(); ++i ) + { + std::cout << strs[i] << std::endl; + } + + GEOSX_ERROR_IF( (m_phaseTypes[ip] != PhaseType::WATER) + && strs.size() != 3, + "Three columns (pressure, formation volume factor, and viscosity) are expected for oil and gas" ); + GEOSX_ERROR_IF( (m_phaseTypes[ip] == PhaseType::WATER) + && strs.size() != 4, + "Four columns (pressure, formation volume factor, compressibility, and viscosity) are expected for water" ); + GEOSX_ERROR_IF( (m_phaseTypes[ip] == PhaseType::WATER) + && counter > 1, + "The water table must contain only one line" ); + + if( m_phaseTypes[ip] == PhaseType::WATER ) + { + m_waterRefPressure = stod( strs[0] ); + m_waterFormationVolFactor = stod( strs[1] ); + m_waterCompressibility = stod( strs[2] ); + m_waterViscosity = stod( strs[3] ); + + GEOSX_ERROR_IF( m_waterRefPressure <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterRefPressureString() ); + GEOSX_ERROR_IF( m_waterFormationVolFactor <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterFormationVolumeFactorString() ); + GEOSX_ERROR_IF( m_waterCompressibility <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterCompressibilityString() ); + GEOSX_ERROR_IF( m_waterViscosity <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterViscosityString() ); + } + else + { + pressureCoords[0].emplace_back( stod( strs[0] ) ); + formationVolFactor.emplace_back( stod( strs[1] ) ); + viscosity.emplace_back( stod( strs[2] ) ); + } + + counter++; + } - // check the table names - localIndex const numExpectedTables = (ipGas >= 0) ? 2 : 1; - GEOSX_ERROR_IF( m_formationVolFactorTableNames.size() != numExpectedTables, - "DeadOilFluid: one formation volume factor table must be provided for each hydrocarbon phase" ); - GEOSX_ERROR_IF( m_viscosityTableNames.size() != numExpectedTables, - "DeadOilFluid: one viscosity table must be provided for each hydrocarbon phase" ); + if( m_phaseTypes[ip] != PhaseType::WATER ) + { + m_hydrocarbonPhaseOrder.emplace_back( LvArray::integerConversion< integer >( ip ) ); + string const formationVolFactorTableName = (m_phaseTypes[ip] == PhaseType::OIL) ? "PVDO_Bo" : "PVDG_Bg"; + string const viscosityTableName = (m_phaseTypes[ip] == PhaseType::OIL) ? "PVDO_visco" : "PVDG_viscg"; + m_formationVolFactorTableNames.emplace_back( formationVolFactorTableName ); + m_viscosityTableNames.emplace_back( viscosityTableName ); + + GEOSX_ERROR_IF( pressureCoords[0].size() <= 1, + "DeadOilFluid: The oil and gas PVT tables must contain at least 2 values" ); + + FunctionManager & functionManager = getGlobalState().getFunctionManager(); + TableFunction & tablePVDX_B = + dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", formationVolFactorTableName ) ); + tablePVDX_B.setTableCoordinates( pressureCoords ); + tablePVDX_B.setTableValues( formationVolFactor ); + tablePVDX_B.reInitializeFunction(); + tablePVDX_B.setInterpolationMethod( TableFunction::InterpolationType::Linear ); + TableFunction & tablePVDX_visc = + dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", viscosityTableName ) ); + tablePVDX_visc.setTableCoordinates( pressureCoords ); + tablePVDX_visc.setTableValues( viscosity ); + tablePVDX_visc.reInitializeFunction(); + tablePVDX_visc.setInterpolationMethod( TableFunction::InterpolationType::Linear ); + } + is.close(); + } // check the size of the additional parameters GEOSX_ERROR_IF( m_surfacePhaseMassDensity.size() != m_phaseNames.size(), @@ -170,11 +231,9 @@ void DeadOilFluid::createAllKernelWrappers() // loop over the hydrocarbon phases for( localIndex iph = 0; iph < m_hydrocarbonPhaseOrder.size(); ++iph ) { - localIndex const ip = m_hydrocarbonPhaseOrder[iph]; - // grab the tables by name from the function manager - TableFunction const & FVFTable = functionManager.getGroup< TableFunction const >( m_formationVolFactorTableNames[ip] ); - TableFunction const & viscosityTable = functionManager.getGroup< TableFunction const >( m_viscosityTableNames[ip] ); + TableFunction const & FVFTable = functionManager.getGroup< TableFunction const >( m_formationVolFactorTableNames[iph] ); + TableFunction const & viscosityTable = functionManager.getGroup< TableFunction const >( m_viscosityTableNames[iph] ); validateTable( FVFTable ); validateTable( viscosityTable ); @@ -190,11 +249,6 @@ void DeadOilFluid::validateTable( TableFunction const & table ) const { array1d< real64 > const & property = table.getValues(); - GEOSX_ERROR_IF( table.getInterpolationMethod() != TableFunction::InterpolationType::Linear, - "DeadOilFluid: the interpolation method for the PVT tables must be linear" ); - GEOSX_ERROR_IF( property.size() <= 1, - "DeadOilFluid: each PVT table must contain at least 2 values" ); - for( localIndex i = 2; i < property.size(); ++i ) { GEOSX_ERROR_IF( (property[i] - property[i-1]) * (property[i-1] - property[i-2]) < 0, diff --git a/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp index 2a61c09c988..06c76860e2d 100644 --- a/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp @@ -121,7 +121,7 @@ class DeadOilFluidUpdate final : public MultiFluidBaseUpdate /// Deleted move assignment operator DeadOilFluidUpdate & operator=( DeadOilFluidUpdate && ) = delete; - + //GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE virtual void compute( real64 const pressure, @@ -330,6 +330,7 @@ class DeadOilFluid : public MultiFluidBase struct viewKeyStruct : MultiFluidBase::viewKeyStruct { static constexpr char const * surfacePhaseMassDensitiesString() { return "surfaceDensities"; } + static constexpr char const * tableFilesString() { return "tableFiles"; } static constexpr char const * formationVolumeFactorTableNamesString() { return "hydrocarbonFormationVolFactorTableNames"; } static constexpr char const * viscosityTableNamesString() { return "hydrocarbonViscosityTableNames"; } static constexpr char const * waterRefPressureString() { return "waterReferencePressure"; } @@ -356,6 +357,11 @@ class DeadOilFluid : public MultiFluidBase // Input data + // Black-oil table filenames + path_array m_tableFiles; + + // Fluid data + /// Names of the formation volume factor tables array1d< string > m_formationVolFactorTableNames; diff --git a/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp b/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp index 54aa4d369a3..4f71c7327ec 100644 --- a/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp +++ b/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp @@ -83,34 +83,32 @@ static const char * pvtw_str = "#\tPref[bar]\tBw[m3/sm3]\tCp[1/bar]\t Visc[cP /// Dead-oil tables written into temporary files during testing -static const char * pvdg_str = "# Pg(Pa)\tBg(m3/sm3)\tVisc(Pa.s)\n" - "\n" - "3000000\t\t0.04234\t\t0.00001344\n" - "6000000\t\t0.02046\t\t0.0000142\n" - "9000000\t\t0.01328\t\t0.00001526\n" - "12000000\t0.00977\t\t0.0000166\n" - "15000000\t0.00773\t\t0.00001818\n" - "18000000\t0.006426\t0.00001994\n" - "21000000\t0.005541\t0.00002181\n" - "24000000\t0.004919\t0.0000237\n" - "27000000\t0.004471\t0.00002559\n" - "29500000\t0.004194\t0.00002714\n" - "31000000\t0.004031\t0.00002806\n" - "33000000\t0.00391\t\t0.00002832\n" - "53000000\t0.003868\t0.00002935"; - -static const char * pvdo_str = "#P[Pa]\tBo[m3/sm3]\tVisc(Pa.s)\n" - "\n" - "2000000\t\t1.02\t0.000975\n" - "5000000\t\t1.03\t0.00091\n" - "10000000\t1.04\t0.00083\n" - "20000000\t1.05\t0.000695\n" - "30000000\t1.07\t0.000594\n" - "40000000\t1.08\t0.00051\n" - "50000000.7\t1.09\t0.000449"; - -static const char * pvdw_str = "#\tPref[bar]\tBw[m3/sm3]\tCp[1/bar]\t Visc[cP]\n" - "\t30600000.1\t1.03\t\t0.00000000041\t0.0003"; +static const char * pvdg_str = "# 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 * pvdo_str = "#P[Pa] Bo[m3/sm3] Visc(Pa.s)\n" + "2000000 1.02 0.000975\n" + "5000000 1.03 0.00091\n" + "10000000 1.04 0.00083\n" + "20000000 1.05 0.000695\n" + "30000000 1.07 0.000594\n" + "40000000 1.08 0.00051\n" + "50000000.7 1.09 0.000449"; + +static const char * pvdw_str = "# Pref[bar] Bw[m3/sm3] Cp[1/bar] Visc[cP]\n" + " 30600000.1 1.03 0.00000000041 0.0003"; void testNumericalDerivatives( MultiFluidBase & fluid, Group & parent, @@ -395,22 +393,17 @@ MultiFluidBase & makeLiveOilFluid( string const & name, Group * parent ) tableNames.resize( 3 ); tableNames[0] = "pvto.txt"; tableNames[1] = "pvtg.txt"; tableNames[2] = "pvtw.txt"; - BlackOilFluid::FluidType & fluidType = fluid.getReference< BlackOilFluid::FluidType >( BlackOilFluid::viewKeyStruct::fluidTypeString() ); - fluidType = BlackOilFluid::FluidType::LiveOil; - fluid.postProcessInputRecursive(); return fluid; } MultiFluidBase & makeDeadOilFluid( string const & name, Group * parent ) { - BlackOilFluid & fluid = parent->registerGroup< BlackOilFluid >( name ); - - // TODO we should actually create a fake XML node with data, but this seemed easier... + DeadOilFluid & fluid = parent->registerGroup< DeadOilFluid >( name ); string_array & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); compNames.resize( 3 ); - compNames[0] = "oil"; compNames[1] = "gas"; compNames[2] = "water"; + compNames[0] = "oil"; compNames[1] = "water"; compNames[2] = "gas"; array1d< real64 > & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); molarWgt.resize( 3 ); @@ -418,130 +411,15 @@ MultiFluidBase & makeDeadOilFluid( string const & name, Group * parent ) string_array & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); phaseNames.resize( 3 ); - phaseNames[0] = "oil"; phaseNames[1] = "gas"; phaseNames[2] = "water"; + phaseNames[0] = "oil"; phaseNames[1] = "water"; phaseNames[2] = "gas"; array1d< real64 > & surfaceDens = fluid.getReference< array1d< real64 > >( BlackOilFluid::viewKeyStruct::surfaceDensitiesString() ); surfaceDens.resize( 3 ); - surfaceDens[0] = 800.0; surfaceDens[1] = 0.9907; surfaceDens[2] = 1022.0; + surfaceDens[0] = 800.0; surfaceDens[1] = 1022.0; surfaceDens[2] = 0.9907; path_array & tableNames = fluid.getReference< path_array >( BlackOilFluid::viewKeyStruct::tableFilesString() ); tableNames.resize( 3 ); - tableNames[0] = "pvdo.txt"; tableNames[1] = "pvdg.txt"; tableNames[2] = "pvdw.txt"; - - BlackOilFluid::FluidType & fluidType = fluid.getReference< BlackOilFluid::FluidType >( BlackOilFluid::viewKeyStruct::fluidTypeString() ); - fluidType = BlackOilFluid::FluidType::DeadOil; - - fluid.postProcessInputRecursive(); - return fluid; -} - -MultiFluidBase & makeNativeDeadOilFluid( string const & name, Group * parent ) -{ - FunctionManager & functionManager = getGlobalState().getFunctionManager(); - - // 1) First, define the tables (PVDO, PVDG) - - // 1D table with linear interpolation - localIndex const NaxisPVDO = 7; - localIndex 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] = 2000000; valuesPVDO_Bo[0] = 1.02; valuesPVDO_visc[0] = 0.000975; - coordinatesPVDO[0][1] = 5000000; valuesPVDO_Bo[1] = 1.03; valuesPVDO_visc[1] = 0.00091; - coordinatesPVDO[0][2] = 10000000; valuesPVDO_Bo[2] = 1.04; valuesPVDO_visc[2] = 0.00083; - coordinatesPVDO[0][3] = 20000000; valuesPVDO_Bo[3] = 1.05; valuesPVDO_visc[3] = 0.000695; - coordinatesPVDO[0][4] = 30000000; valuesPVDO_Bo[4] = 1.07; valuesPVDO_visc[4] = 0.000594; - coordinatesPVDO[0][5] = 40000000; valuesPVDO_Bo[5] = 1.08; valuesPVDO_visc[5] = 0.00051; - coordinatesPVDO[0][6] = 50000000; valuesPVDO_Bo[6] = 1.09; valuesPVDO_visc[6] = 0.000449; - - 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 ); - tablePVDO_Bo.setTableValues( valuesPVDO_Bo ); - tablePVDO_Bo.reInitializeFunction(); - - tablePVDO_Bo.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - - TableFunction & tablePVDO_visc = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDO_visc" ) ); - tablePVDO_visc.setTableCoordinates( coordinatesPVDO ); - tablePVDO_visc.setTableValues( valuesPVDO_visc ); - tablePVDO_visc.reInitializeFunction(); - - tablePVDO_visc.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - - TableFunction & tablePVDG_Bg = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDG_Bg" ) ); - tablePVDG_Bg.setTableCoordinates( coordinatesPVDG ); - tablePVDG_Bg.setTableValues( valuesPVDG_Bg ); - tablePVDG_Bg.reInitializeFunction(); - - tablePVDG_Bg.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - - TableFunction & tablePVDG_visc = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDG_visc" ) ); - tablePVDG_visc.setTableCoordinates( coordinatesPVDG ); - tablePVDG_visc.setTableValues( valuesPVDG_visc ); - tablePVDG_visc.reInitializeFunction(); - - tablePVDG_visc.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - - // 2) Then, define the Dead-Oil constitutive model - - DeadOilFluid & fluid = parent->registerGroup< DeadOilFluid >( name ); - - auto & compNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::componentNamesString() ); - compNames.resize( 3 ); - compNames[0] = "oil"; compNames[1] = "gas"; compNames[2] = "water"; - - auto & molarWgt = fluid.getReference< array1d< real64 > >( MultiFluidBase::viewKeyStruct::componentMolarWeightString() ); - molarWgt.resize( 3 ); - molarWgt[0] = 114e-3; molarWgt[1] = 16e-3; molarWgt[2] = 18e-3; - - auto & phaseNames = fluid.getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); - phaseNames.resize( 3 ); - phaseNames[0] = "oil"; phaseNames[1] = "gas"; phaseNames[2] = "water"; - - auto & surfaceDens = fluid.getReference< array1d< real64 > >( DeadOilFluid::viewKeyStruct::surfacePhaseMassDensitiesString() ); - surfaceDens.resize( 3 ); - surfaceDens[0] = 800.0; surfaceDens[1] = 0.9907; surfaceDens[2] = 1022.0; - - auto & FVFTableNames = fluid.getReference< string_array >( DeadOilFluid::viewKeyStruct::formationVolumeFactorTableNamesString() ); - FVFTableNames.resize( 2 ); - FVFTableNames[0] = "PVDO_Bo"; FVFTableNames[1] = "PVDG_Bg"; - - auto & viscosityTableNames = fluid.getReference< string_array >( DeadOilFluid::viewKeyStruct::viscosityTableNamesString() ); - viscosityTableNames.resize( 2 ); - viscosityTableNames[0] = "PVDO_visc"; viscosityTableNames[1] = "PVDG_visc"; - - auto & waterRefPressure = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterRefPressureString() ); - waterRefPressure = 30600000.1; - auto & waterFormationVolumeFactor = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterFormationVolumeFactorString() ); - waterFormationVolumeFactor = 1.03; - auto & waterCompressibility = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterCompressibilityString() ); - waterCompressibility = 0.00000000041; - auto & waterViscosity = fluid.getReference< real64 >( DeadOilFluid::viewKeyStruct::waterViscosityString() ); - waterViscosity = 0.0003; + tableNames[0] = "pvdo.txt"; tableNames[1] = "pvdw.txt"; tableNames[2] = "pvdg.txt"; fluid.postProcessInputRecursive(); return fluid; @@ -647,56 +525,6 @@ TEST_F( DeadOilFluidTest, numericalDerivativesMolar ) { fluid->setMassFlag( false ); - // TODO test over a range of values - real64 const P = 5e6; - 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; - - testNumericalDerivatives( *fluid, parent, P, T, comp, eps, relTol ); -} - -TEST_F( DeadOilFluidTest, numericalDerivativesMass ) -{ - fluid->setMassFlag( true ); - - // TODO test over a range of values - real64 const P = 5e6; - 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-2; - real64 const absTol = 1e-14; - - testNumericalDerivatives( *fluid, parent, P, T, comp, eps, relTol, absTol ); -} - -class NativeDeadOilFluidTest : public CompositionalFluidTestBase -{ -public: - - NativeDeadOilFluidTest() - { - parent.resize( 1 ); - fluid = &makeNativeDeadOilFluid( "fluid", &parent ); - - parent.initialize(); - parent.initializePostInitialConditions(); - } - - ~NativeDeadOilFluidTest() - {} -}; - -TEST_F( NativeDeadOilFluidTest, numericalDerivativesMolar ) -{ - fluid->setMassFlag( false ); - real64 const P[3] = { 5.4e6, 1.24e7, 3.21e7 }; real64 const T = 297.15; array1d< real64 > comp( 3 ); @@ -711,7 +539,7 @@ TEST_F( NativeDeadOilFluidTest, numericalDerivativesMolar ) } } -TEST_F( NativeDeadOilFluidTest, numericalDerivativesMass ) +TEST_F( DeadOilFluidTest, numericalDerivativesMass ) { fluid->setMassFlag( true ); diff --git a/src/coreComponents/fileIO/schema/docs/BlackOilFluid.rst b/src/coreComponents/fileIO/schema/docs/BlackOilFluid.rst index d10188e874a..900e19353fd 100644 --- a/src/coreComponents/fileIO/schema/docs/BlackOilFluid.rst +++ b/src/coreComponents/fileIO/schema/docs/BlackOilFluid.rst @@ -1,17 +1,14 @@ -==================== ========================================== ======== ============================================================= -Name Type Default Description -==================== ========================================== ======== ============================================================= -componentMolarWeight real64_array required Component molar weights -componentNames string_array {} List of component names -fluidType geosx_constitutive_BlackOilFluid_FluidType required | Type of black-oil fluid. Valid options: - | * DeadOil - | * LiveOil -name string required A name is required for any non-unique nodes -phaseNames string_array required List of fluid phases -surfaceDensities real64_array required List of surface densities for each phase -tableFiles path_array required List of filenames with input PVT tables -==================== ========================================== ======== ============================================================= +==================== ============ ======== =========================================== +Name Type Default Description +==================== ============ ======== =========================================== +componentMolarWeight real64_array required Component molar weights +componentNames string_array {} List of component names +name string required A name is required for any non-unique nodes +phaseNames string_array required List of fluid phases +surfaceDensities real64_array required List of surface densities for each phase +tableFiles path_array required List of filenames with input PVT tables +==================== ============ ======== =========================================== diff --git a/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst b/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst index 10eb6506f18..93045746e8f 100644 --- a/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst +++ b/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst @@ -1,21 +1,14 @@ -======================================= ============ ======== =============================================================================================================================================================================================================================== -Name Type Default Description -======================================= ============ ======== =============================================================================================================================================================================================================================== -componentMolarWeight real64_array required Component molar weights -componentNames string_array {} List of component names -hydrocarbonFormationVolFactorTableNames string_array required | Names of the formation volume factor tables (one per hydrocarbon phase, in the order provided in "phaseNames"). - | For instance, if "oil" is before "gas" in "phaseNames", the table order should be: oilTableName, gasTableName -hydrocarbonViscosityTableNames string_array required | Names of the viscosity tables (one per hydrocarbon phase, in the order provided in "phaseNames") - | For instance, if "oil" is before "gas" in "phaseNames", the table order should be: oilTableName, gasTableName -name string required A name is required for any non-unique nodes -phaseNames string_array required List of fluid phases -surfaceDensities real64_array required List of surface mass densities for each phase -waterCompressibility real64 0 Water compressibility -waterFormationVolumeFactor real64 0 Water formation volume factor -waterReferencePressure real64 0 Water reference pressure -waterViscosity real64 0 Water viscosity -======================================= ============ ======== =============================================================================================================================================================================================================================== +==================== ============ ======== ============================================= +Name Type Default Description +==================== ============ ======== ============================================= +componentMolarWeight real64_array required Component molar weights +componentNames string_array {} List of component names +name string required A name is required for any non-unique nodes +phaseNames string_array required List of fluid phases +surfaceDensities real64_array required List of surface mass densities for each phase +tableFiles path_array required List of filenames with input PVT tables +==================== ============ ======== ============================================= diff --git a/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst b/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst index 5b3d3c27bb3..7d677303016 100644 --- a/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst +++ b/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst @@ -1,33 +1,39 @@ -====================================== ============================================================================================== ========================== -Name Type Description -====================================== ============================================================================================== ========================== -dPhaseCompFraction_dGlobalCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l, 3l, 4l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseCompFraction_dPressure LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseCompFraction_dTemperature LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseDensity_dPressure real64_array3d (no description available) -dPhaseDensity_dTemperature real64_array3d (no description available) -dPhaseFraction_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseFraction_dPressure real64_array3d (no description available) -dPhaseFraction_dTemperature real64_array3d (no description available) -dPhaseMassDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseMassDensity_dPressure real64_array3d (no description available) -dPhaseMassDensity_dTemperature real64_array3d (no description available) -dPhaseViscosity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseViscosity_dPressure real64_array3d (no description available) -dPhaseViscosity_dTemperature real64_array3d (no description available) -dTotalDensity_dGlobalCompFraction real64_array3d (no description available) -dTotalDensity_dPressure real64_array2d (no description available) -dTotalDensity_dTemperature real64_array2d (no description available) -phaseCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -phaseDensity real64_array3d (no description available) -phaseFraction real64_array3d (no description available) -phaseMassDensity real64_array3d (no description available) -phaseViscosity real64_array3d (no description available) -totalDensity real64_array2d (no description available) -useMass integer (no description available) -====================================== ============================================================================================== ========================== +======================================= ============================================================================================== ========================== +Name Type Description +======================================= ============================================================================================== ========================== +dPhaseCompFraction_dGlobalCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l, 3l, 4l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseCompFraction_dPressure LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseCompFraction_dTemperature LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseDensity_dPressure real64_array3d (no description available) +dPhaseDensity_dTemperature real64_array3d (no description available) +dPhaseFraction_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseFraction_dPressure real64_array3d (no description available) +dPhaseFraction_dTemperature real64_array3d (no description available) +dPhaseMassDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseMassDensity_dPressure real64_array3d (no description available) +dPhaseMassDensity_dTemperature real64_array3d (no description available) +dPhaseViscosity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseViscosity_dPressure real64_array3d (no description available) +dPhaseViscosity_dTemperature real64_array3d (no description available) +dTotalDensity_dGlobalCompFraction real64_array3d (no description available) +dTotalDensity_dPressure real64_array2d (no description available) +dTotalDensity_dTemperature real64_array2d (no description available) +hydrocarbonFormationVolFactorTableNames string_array (no description available) +hydrocarbonViscosityTableNames string_array (no description available) +phaseCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +phaseDensity real64_array3d (no description available) +phaseFraction real64_array3d (no description available) +phaseMassDensity real64_array3d (no description available) +phaseViscosity real64_array3d (no description available) +totalDensity real64_array2d (no description available) +useMass integer (no description available) +waterCompressibility real64 (no description available) +waterFormationVolumeFactor real64 (no description available) +waterReferencePressure real64 (no description available) +waterViscosity real64 (no description available) +======================================= ============================================================================================== ========================== diff --git a/src/coreComponents/fileIO/schema/schema.xsd b/src/coreComponents/fileIO/schema/schema.xsd index 3e8f2fbfa01..b71a81efc78 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd +++ b/src/coreComponents/fileIO/schema/schema.xsd @@ -1603,10 +1603,6 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> - - @@ -1616,11 +1612,6 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> - - - - - @@ -1815,24 +1806,12 @@ The expected format is "{ waterMax, oilMax }", in that order--> - - - - - - - - - - - - + + diff --git a/src/coreComponents/fileIO/schema/schema.xsd.other b/src/coreComponents/fileIO/schema/schema.xsd.other index 65ed62683a1..1d703090ebf 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd.other +++ b/src/coreComponents/fileIO/schema/schema.xsd.other @@ -950,6 +950,10 @@ + + + + @@ -964,6 +968,14 @@ + + + + + + + + diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/dead_oil_egg.xml b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/dead_oil_egg.xml index 6d054c35b08..e921fe40261 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/dead_oil_egg.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/dead_oil_egg.xml @@ -12,17 +12,19 @@ targetRegions="{ reservoir, wellRegion1, wellRegion2, wellRegion3, wellRegion4, wellRegion5, wellRegion6, wellRegion7, wellRegion8, wellRegion9, wellRegion10, wellRegion11, wellRegion12 }"> - - + solverType="fgmres" + preconditionerType="mgr" + krylovTol="1e-4" + krylovAdaptiveTol="1" + krylovWeakestTol="1e-2"/> + injectionStream="{ 0.0, 1.0 }"/> + injectionStream="{ 0.0, 1.0 }"/> + injectionStream="{ 0.0, 1.0 }"/> + injectionStream="{ 0.0, 1.0 }"/> + injectionStream="{ 0.0, 1.0 }"/> + injectionStream="{ 0.0, 1.0 }"/> + injectionStream="{ 0.0, 1.0 }"/> + injectionStream="{ 0.0, 1.0 }"/> @@ -235,7 +237,7 @@ distanceFromHead="2"/> + distanceFromHead="6"/> @@ -592,25 +594,25 @@ @@ -687,20 +689,19 @@ - + phaseNames="{ oil, water }" + surfaceDensities="{ 848.9, 1025.2 }" + componentMolarWeight="{ 114e-3, 18e-3 }" + tableFiles="{ pvdo.txt, pvtw.txt }"/> + phaseNames="{ oil, water }" + phaseMinVolumeFraction="{ 0.1, 0.2 }" + phaseRelPermExponent="{ 4.0, 3.0 }" + phaseRelPermMaxValue="{ 0.8, 0.75 }"/> - - + scale="0.9"/> + scale="0.1"/> diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvdg.txt b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvdg.txt deleted file mode 100644 index 791e9cb2bb0..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvdg.txt +++ /dev/null @@ -1,15 +0,0 @@ -# Pg(Pa) Bg(m3/sm3) Visc(Pa.s) -3000000 0.04234 0.00001344 -6000000 0.02046 0.0000142 -9000000 0.01328 0.00001526 -12000000 0.00977 0.0000166 -15000000 0.00773 0.00001818 -18000000 0.006426 0.00001994 -21000000 0.005541 0.00002181 -24000000 0.004919 0.0000237 -27000000 0.004471 0.00002559 -29500000 0.004194 0.00002714 -31000000 0.004031 0.00002806 -33000000 0.00391 0.00002832 -53000000 0.003868 0.00002935 -530000000 0.003768 0.00002945 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvdo.txt b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvdo.txt index 92a78f67ce5..ad7ab2e4c0c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvdo.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvdo.txt @@ -1,6 +1,6 @@ -# P[Pa] Bo[m3/sm3] Visc(Pa.s) -2068000 1.05 0.001 -5516000 1.02 0.001 -55158000 1.01 0.001 -95158000 1.00 0.001 -951580000 0.99 0.001 +# P[Pa] Bo[m3/sm3] Visc(Pa.s) +2068000 1.05 0.005 +5516000 1.02 0.005 +55158000 1.01 0.005 +95158000 1.00 0.005 +951580000 0.99 0.005 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvtw.txt b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvtw.txt index a2148c7926d..67aa2b7c19d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvtw.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/pvtw.txt @@ -1,3 +1,2 @@ -# Pref[Pa] Bw[m3/sm3] Cp[1/Pa] Visc[Pa.s] - 30600000.1 1.01 0.0000000001 0.001 - +#Pref[Pa] Bw[m3/sm3] Cp[1/Pa] Visc[Pa.s] +30600000.1 1.01 0.0000000001 0.001 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/dead_oil_spe10_layers_83_84_85.xml b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/dead_oil_spe10_layers_83_84_85.xml index 0ccea452ac8..69906c77049 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/dead_oil_spe10_layers_83_84_85.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/dead_oil_spe10_layers_83_84_85.xml @@ -16,10 +16,12 @@ maxTimeStepCuts="10" lineSearchAction="None" newtonMaxIter="20"/> - - + solverType="fgmres" + preconditionerType="mgr" + krylovTol="1e-4" + krylovAdaptiveTol="1" + krylovWeakestTol="1e-2"/> + injectionStream="{ 0.0, 1.0 }"/> @@ -203,7 +205,9 @@ + maxTime="2e6"> + + - + phaseNames="{ oil, water }" + surfaceDensities="{ 848.9, 1025.2 }" + componentMolarWeight="{ 114e-3, 18e-3 }" + tableFiles="{ pvdo.txt, pvtw.txt }"/> + phaseNames="{ oil, water }" + phaseMinVolumeFraction="{ 0.2, 0.2 }" + phaseRelPermExponent="{ 2.0, 2.0 }" + phaseRelPermMaxValue="{ 0.1, 1.0 }"/> - - diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvdg.txt b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvdg.txt deleted file mode 100644 index ab829cbbba1..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvdg.txt +++ /dev/null @@ -1,14 +0,0 @@ -# Pg(Pa) Bg(m3/sm3) Visc(Pa.s) -3000000 0.04234 0.00001344 -6000000 0.02046 0.0000142 -9000000 0.01328 0.00001526 -12000000 0.00977 0.0000166 -15000000 0.00773 0.00001818 -18000000 0.006426 0.00001994 -21000000 0.005541 0.00002181 -24000000 0.004919 0.0000237 -27000000 0.004471 0.00002559 -29500000 0.004194 0.00002714 -31000000 0.004031 0.00002806 -33000000 0.00391 0.00002832 -53000000 0.003868 0.00002935 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvdo.txt b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvdo.txt index 2dbc106fa01..04a70d8aac2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvdo.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvdo.txt @@ -1,4 +1,4 @@ -# P[Pa] Bo[m3/sm3] Visc(Pa.s) -2068000 1.05 0.001 -5516000 1.02 0.001 -55158000 1.01 0.001 +# P[Pa] Bo[m3/sm3] Visc(Pa.s) +2068000 1.05 0.00285 +5516000 1.02 0.00299 +55158000 1.01 0.00300 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvtw.txt b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvtw.txt index f19aa44f3f2..d4e05aa9dc9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvtw.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/SPE10/pvtw.txt @@ -1,2 +1,2 @@ -# Pref[Pa] Bw[m3/sm3] Cp[1/Pa] Visc[Pa.s] - 30600000.1 1.01 0.0000000004 0.001 +#Pref[Pa] Bw[m3/sm3] Cp[1/Pa] Visc[Pa.s] +30600000.1 1.01 0.0000000004 0.0003 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/blackoil_2ph_1d.xml b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/blackoil_2ph_1d.xml index 8b3ab95994d..2bc918d2dc1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/blackoil_2ph_1d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/blackoil_2ph_1d.xml @@ -103,7 +103,6 @@ + materialList="{ fluid, rock, relperm }"/> - + materialList="{ fluid, rock, relperm }"/> - + materialList="{ fluid, rock, relperm }"/> - @@ -80,13 +80,12 @@ - + phaseNames="{ oil, gas }" + surfaceDensities="{ 800.0, 0.9907 }" + componentMolarWeight="{ 114e-3, 16e-3 }" + tableFiles="{ pvdo.txt, pvdg.txt }"/> + phaseNames="{ oil, gas }" + phaseMinVolumeFraction="{ 0.05, 0.05 }" + phaseRelPermExponent="{ 1.5, 1.5 }" + phaseRelPermMaxValue="{ 0.9, 0.9 }"/> @@ -165,14 +164,6 @@ component="1" scale="0.4"/> - diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdg.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdg.txt index ab829cbbba1..b25fe86e44e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdg.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdg.txt @@ -1,14 +1,14 @@ -# Pg(Pa) Bg(m3/sm3) Visc(Pa.s) -3000000 0.04234 0.00001344 -6000000 0.02046 0.0000142 -9000000 0.01328 0.00001526 -12000000 0.00977 0.0000166 -15000000 0.00773 0.00001818 -18000000 0.006426 0.00001994 -21000000 0.005541 0.00002181 -24000000 0.004919 0.0000237 -27000000 0.004471 0.00002559 -29500000 0.004194 0.00002714 -31000000 0.004031 0.00002806 -33000000 0.00391 0.00002832 -53000000 0.003868 0.00002935 +# Pg(Pa) Bg(m3/sm3) Visc(Pa.s) +3000000 0.04234 0.00001344 +6000000 0.02046 0.0000142 +9000000 0.01328 0.00001526 +12000000 0.00977 0.0000166 +15000000 0.00773 0.00001818 +18000000 0.006426 0.00001994 +21000000 0.005541 0.00002181 +24000000 0.004919 0.0000237 +27000000 0.004471 0.00002559 +29500000 0.004194 0.00002714 +31000000 0.004031 0.00002806 +33000000 0.00391 0.00002832 +53000000 0.003868 0.00002935 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdo.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdo.txt index a05a6698b66..df0be305c70 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdo.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdo.txt @@ -1,8 +1,8 @@ -# P[Pa] Bo[m3/sm3] Visc(Pa.s) -2000000 1.02 0.000975 -5000000 1.03 0.00091 -10000000 1.04 0.00083 -20000000 1.05 0.000695 -30000000 1.07 0.000594 -40000000 1.08 0.00051 -50000000.7 1.09 0.000449 +# P[Pa] Bo[m3/sm3] Visc(Pa.s) +2000000 1.02 0.000975 +5000000 1.03 0.00091 +10000000 1.04 0.00083 +20000000 1.05 0.000695 +30000000 1.07 0.000594 +40000000 1.08 0.00051 +50000000.7 1.09 0.000449 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvtw.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvtw.txt index 610f671eced..4762e1c3636 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvtw.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvtw.txt @@ -1,2 +1,2 @@ -# Pref[bar] Bw[m3/sm3] Cp[1/bar] Visc[cP] - 30600000.1 1.03 0.00000000041 0.0003 \ No newline at end of file +#Pref[bar] Bw[m3/sm3] Cp[1/bar] Visc[cP] +30600000.1 1.03 0.00000000041 0.0003 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/dead_oil_wells_2d.xml b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/dead_oil_wells_2d.xml index a647abf2d59..8e18b4ca428 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/dead_oil_wells_2d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/dead_oil_wells_2d.xml @@ -200,9 +200,8 @@ - Date: Thu, 25 Feb 2021 12:35:12 -0800 Subject: [PATCH 3/6] small fixes --- .../constitutive/fluid/BlackOilFluid.hpp | 2 +- .../constitutive/fluid/DeadOilFluid.hpp | 4 ++-- .../managers/Functions/TableFunction.hpp | 11 +++++++---- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp b/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp index f0687081200..8ac8a2200df 100644 --- a/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp +++ b/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp @@ -1,5 +1,5 @@ /* - 1;5202;0c * ------------------------------------------------------------------------------------------------------------ + * ------------------------------------------------------------------------------------------------------------ * SPDX-License-Identifier: LGPL-2.1-only * * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC diff --git a/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp index 06c76860e2d..c3ef1de6374 100644 --- a/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp @@ -444,7 +444,7 @@ void DeadOilFluidUpdate::computeDensities( real64 pressure, real64 fvf = 0.0; real64 derivative = 0.0; - for( double & val : dPhaseMassDens_dGlobalCompFrac ) + for( real64 & val : dPhaseMassDens_dGlobalCompFrac ) { val = 0.0; } @@ -521,7 +521,7 @@ void DeadOilFluidUpdate::computeViscosities( real64 pressure, arraySlice1d< real64 > const & dPhaseVisc_dPres, arraySlice2d< real64 > const & dPhaseVisc_dGlobalCompFrac ) const { - for( double & val : dPhaseVisc_dGlobalCompFrac ) + for( real64 & val : dPhaseVisc_dGlobalCompFrac ) { val = 0.0; } diff --git a/src/coreComponents/managers/Functions/TableFunction.hpp b/src/coreComponents/managers/Functions/TableFunction.hpp index 248458e6415..a81a04f3593 100644 --- a/src/coreComponents/managers/Functions/TableFunction.hpp +++ b/src/coreComponents/managers/Functions/TableFunction.hpp @@ -47,6 +47,9 @@ class TableFunction : public FunctionBase /// maximum dimensions for the coordinates in the table static constexpr localIndex maxDimensions = 4; + /// maximum number of corners + static constexpr localIndex maxNumCorners = 16; + /** * @class KernelWrapper * @@ -73,7 +76,7 @@ class TableFunction : public FunctionBase localIndex dimensions, arrayView1d< localIndex const > const & size, arrayView1d< localIndex const > const & indexIncrement, - localIndex const (&corners)[TableFunction::maxDimensions][16], + localIndex const (&corners)[TableFunction::maxDimensions][maxNumCorners], localIndex const numCorners ); /// Default constructor for the kernel wrapper @@ -108,7 +111,7 @@ class TableFunction : public FunctionBase localIndex dimensions, arrayView1d< localIndex const > const & size, arrayView1d< localIndex const > const & indexIncrement, - localIndex const (&corners)[TableFunction::maxDimensions][16], + localIndex const (&corners)[TableFunction::maxDimensions][maxNumCorners], localIndex const numCorners ); /** @@ -146,7 +149,7 @@ class TableFunction : public FunctionBase * @brief The corners of the box that surround the value in N dimensions * m_corners should be of size m_maxDimensions x (2^m_maxDimensions) */ - localIndex m_corners[TableFunction::maxDimensions][16]; + localIndex m_corners[TableFunction::maxDimensions][maxNumCorners]; /// The number of active table corners localIndex m_numCorners; @@ -299,7 +302,7 @@ class TableFunction : public FunctionBase * @brief The corners of the box that surround the value in N dimensions * m_corners should be of size maxDimensions x (2^maxDimensions) */ - localIndex m_corners[maxDimensions][16]; + localIndex m_corners[maxDimensions][maxNumCorners]; /// The number of active table corners localIndex m_numCorners; From 10c3c5fae3a5e44987c4d6770180c27bfb12952a Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Thu, 25 Feb 2021 12:54:25 -0800 Subject: [PATCH 4/6] post-merge fixes --- .../deadoil_3ph_staircase_hybrid_3d.xml | 9 ++++----- .../dead_oil_wells_hybrid_2d.xml | 17 ++++++++--------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_staircase_hybrid_3d.xml b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_staircase_hybrid_3d.xml index 749fdef1ab5..a43237ff473 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_staircase_hybrid_3d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_staircase_hybrid_3d.xml @@ -7,7 +7,7 @@ logLevel="1" discretization="fluidHM" targetRegions="{ Channel }" - fluidNames="{ fluid1 }" + fluidNames="{ fluid }" solidNames="{ rock }" relPermNames="{ relperm }" temperature="300" @@ -93,7 +93,7 @@ + materialList="{ fluid, rock, relperm }"/> - + materialList="{ fluid, rock, relperm }"/> + materialList="{ fluid, relperm }"/> + materialList="{ fluid, relperm }"/> + materialList="{ fluid, relperm }"/> - Date: Tue, 2 Mar 2021 16:32:42 -0800 Subject: [PATCH 5/6] addressed review comments --- .../constitutive/fluid/DeadOilFluid.cpp | 410 ++++++++++++------ .../constitutive/fluid/DeadOilFluid.hpp | 53 ++- .../fluid/MultiFluidPVTPackageWrapper.cpp | 16 +- .../constitutive/unitTests/testMultiFluid.cpp | 139 ++++++ .../fileIO/schema/docs/DeadOilFluid.rst | 30 +- .../fileIO/schema/docs/DeadOilFluid_other.rst | 64 ++- src/coreComponents/fileIO/schema/schema.xsd | 20 +- .../fileIO/schema/schema.xsd.other | 12 - .../FieldSpecificationOps.hpp | 14 +- .../compositionalMultiphaseFlow/B_g_pvdg.txt | 13 + .../compositionalMultiphaseFlow/B_o_pvdo.txt | 7 + .../deadoil_3ph_baker_1d.xml | 44 +- .../deadoil_3ph_corey_1d.xml | 50 +-- .../deadoil_3ph_staircase_3d.xml | 58 ++- .../deadoil_3ph_staircase_hybrid_3d.xml | 14 +- .../compositionalMultiphaseFlow/pres_pvdg.txt | 13 + .../compositionalMultiphaseFlow/pres_pvdo.txt | 7 + .../compositionalMultiphaseFlow/pvdg.txt | 26 +- .../compositionalMultiphaseFlow/visc_pvdg.txt | 13 + .../compositionalMultiphaseFlow/visc_pvdo.txt | 7 + .../compositionalMultiphaseWell/pvdo.txt | 16 +- 21 files changed, 730 insertions(+), 296 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/B_g_pvdg.txt create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/B_o_pvdo.txt create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pres_pvdg.txt create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pres_pvdo.txt create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/visc_pvdg.txt create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/visc_pvdo.txt diff --git a/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp b/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp index 819fbe144e8..de902bae415 100644 --- a/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp @@ -56,23 +56,261 @@ DeadOilFluid::DeadOilFluid( string const & name, getWrapperBase( viewKeyStruct::componentMolarWeightString() ).setInputFlag( InputFlags::REQUIRED ); getWrapperBase( viewKeyStruct::phaseNamesString() ).setInputFlag( InputFlags::REQUIRED ); - registerWrapper( viewKeyStruct::tableFilesString(), &m_tableFiles ). - setInputFlag( InputFlags::REQUIRED ). - setRestartFlags( RestartFlags::NO_WRITE ). - setDescription( "List of filenames with input PVT tables" ); - registerWrapper( viewKeyStruct::surfacePhaseMassDensitiesString(), &m_surfacePhaseMassDensity ). setInputFlag( InputFlags::REQUIRED ). setDescription( "List of surface mass densities for each phase" ); + // 1) First option: specify PVT tables from one file per phase, read the files line by line, and populate the internal TableFunctions + + registerWrapper( viewKeyStruct::tableFilesString(), &m_tableFiles ). + setInputFlag( InputFlags::OPTIONAL ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( "List of filenames with input PVT tables (one per phase)" ); + + // 2) Second option: specify TableFunction names for oil and gas, + // and then specify water reference pressure, water Bw, water compressibility and water viscosity + registerWrapper( viewKeyStruct::formationVolumeFactorTableNamesString(), &m_formationVolFactorTableNames ). - setSizedFromParent( 0 ); + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "List of formation volume factor TableFunction names from the Functions block. \n" + "The user must provide one TableFunction per hydrocarbon phase, in the order provided in \"phaseNames\". \n" + "For instance, if \"oil\" is before \"gas\" in \"phaseNames\", the table order should be: oilTableName, gasTableName" ); registerWrapper( viewKeyStruct::viscosityTableNamesString(), &m_viscosityTableNames ). - setSizedFromParent( 0 );; - registerWrapper( viewKeyStruct::waterRefPressureString(), &m_waterRefPressure ); - registerWrapper( viewKeyStruct::waterFormationVolumeFactorString(), &m_waterFormationVolFactor ); - registerWrapper( viewKeyStruct::waterCompressibilityString(), &m_waterCompressibility ); - registerWrapper( viewKeyStruct::waterViscosityString(), &m_waterViscosity ); + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "List of viscosity TableFunction names from the Functions block. \n" + "The user must provide one TableFunction per hydrocarbon phase, in the order provided in \"phaseNames\". \n" + "For instance, if \"oil\" is before \"gas\" in \"phaseNames\", the table order should be: oilTableName, gasTableName" ); + registerWrapper( viewKeyStruct::waterRefPressureString(), &m_waterRefPressure ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Water reference pressure" ); + registerWrapper( viewKeyStruct::waterFormationVolumeFactorString(), &m_waterFormationVolFactor ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Water formation volume factor" ); + registerWrapper( viewKeyStruct::waterCompressibilityString(), &m_waterCompressibility ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Water compressibility" ); + registerWrapper( viewKeyStruct::waterViscosityString(), &m_waterViscosity ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Water viscosity" ); +} + +void DeadOilFluid::readTable( string const & filename, + array1d< array1d< real64 > > & values ) const +{ + std::ifstream is( filename ); + constexpr std::streamsize buf_size = 256; + char buf[buf_size]; + + while( is.getline( buf, buf_size ) ) + { + string const str( buf ); + string_array const strs = Tokenize( str, " " ); + + if( strs.empty() ) + { + continue; + } + if( strs[0].front() == '#' ) + { + continue; + } + if( streq( strs[0], "--" ) ) + { + continue; + } + + array1d< real64 > rowValues; + rowValues.resize( strs.size() ); + for( localIndex i = 0; i < strs.size(); ++i ) + { + rowValues[i] = stod( strs[i] ); + } + values.emplace_back( rowValues ); + } + is.close(); +} + +void DeadOilFluid::fillWaterData( array1d< array1d< real64 > > const & tableValues ) +{ + GEOSX_THROW_IF( tableValues.size() != 1, + "DeadOilFluid: the water table must contain one line and only one", + InputError ); + GEOSX_THROW_IF( tableValues[0].size() != 4, + "DeadOilFluid: four columns (pressure, formation volume factor, compressibility, and viscosity) are expected for water", + InputError ); + + GEOSX_THROW_IF( m_waterRefPressure > 0.0 || m_waterFormationVolFactor > 0.0 + || m_waterCompressibility > 0.0 || m_waterViscosity > 0.0, + "DeadOilFluid: input is redundant (user provided both water data and a water pvt file)", + InputError ); + + m_waterRefPressure = tableValues[0][0]; + m_waterFormationVolFactor = tableValues[0][1]; + m_waterCompressibility = tableValues[0][2]; + m_waterViscosity = tableValues[0][3]; + + GEOSX_THROW_IF( m_waterRefPressure <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterRefPressureString(), + InputError ); + GEOSX_THROW_IF( m_waterFormationVolFactor <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterFormationVolumeFactorString(), + InputError ); + GEOSX_THROW_IF( m_waterCompressibility <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterCompressibilityString(), + InputError ); + GEOSX_THROW_IF( m_waterViscosity <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterViscosityString(), + InputError ); +} + +void DeadOilFluid::fillHydrocarbonData( localIndex const ip, + array1d< array1d< real64 > > const & tableValues ) +{ + array1d< array1d< real64 > > pressureCoords; + pressureCoords.resize( 1 ); + pressureCoords[0].resize( tableValues.size() ); + array1d< real64 > formationVolFactor; + formationVolFactor.resize( tableValues.size() ); + array1d< real64 > viscosity; + viscosity.resize( tableValues.size() ); + + for( localIndex i = 0; i < tableValues.size(); ++i ) + { + GEOSX_THROW_IF( tableValues[i].size() != 3, + "DeadOilFluid: three columns (pressure, formation volume factor, and viscosity) are expected for oil and gas", + InputError ); + + pressureCoords[0][i] = tableValues[i][0]; + formationVolFactor[i] = tableValues[i][1]; + viscosity[i] = tableValues[i][2]; + } + + m_hydrocarbonPhaseOrder.emplace_back( LvArray::integerConversion< integer >( ip ) ); + string const formationVolFactorTableName = (m_phaseTypes[ip] == PhaseType::OIL) ? "PVDO_Bo" : "PVDG_Bg"; + string const viscosityTableName = (m_phaseTypes[ip] == PhaseType::OIL) ? "PVDO_visco" : "PVDG_viscg"; + m_formationVolFactorTableNames.emplace_back( formationVolFactorTableName ); + m_viscosityTableNames.emplace_back( viscosityTableName ); + + FunctionManager & functionManager = getGlobalState().getFunctionManager(); + TableFunction & tablePVDX_B = + dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", formationVolFactorTableName ) ); + tablePVDX_B.setTableCoordinates( pressureCoords ); + tablePVDX_B.setTableValues( formationVolFactor ); + tablePVDX_B.reInitializeFunction(); + tablePVDX_B.setInterpolationMethod( TableFunction::InterpolationType::Linear ); + TableFunction & tablePVDX_visc = + dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", viscosityTableName ) ); + tablePVDX_visc.setTableCoordinates( pressureCoords ); + tablePVDX_visc.setTableValues( viscosity ); + tablePVDX_visc.reInitializeFunction(); + tablePVDX_visc.setInterpolationMethod( TableFunction::InterpolationType::Linear ); +} + +void DeadOilFluid::readInputDataFromPVTFiles() +{ + GEOSX_THROW_IF( m_tableFiles.size() != numFluidPhases(), + "DeadOilFluid: the number of table files must be equal to the number of phases", + InputError ); + GEOSX_THROW_IF( m_formationVolFactorTableNames.size() > 0.0 || m_viscosityTableNames.size() > 0.0, + "DeadOilFluid: input is redundant (user provided both TableFunction names and pvt files)", + InputError ); + + array1d< array1d< real64 > > tableValues; + for( localIndex ip = 0; ip < numFluidPhases(); ++ip ) + { + tableValues.clear(); + readTable( m_tableFiles[ip], tableValues ); + + if( m_phaseTypes[ip] == PhaseType::WATER ) + { + fillWaterData( tableValues ); + } + else + { + fillHydrocarbonData( ip, tableValues ); + } + } +} + +void DeadOilFluid::useProvidedTableFunctions() +{ + GEOSX_THROW_IF( m_tableFiles.size() > 0, + "DeadOilFluid: input is redundant (user provided both TableFunction names and pvt files)", + InputError ); + + integer const ipWater = m_phaseOrder[PhaseType::WATER]; + integer const ipGas = m_phaseOrder[PhaseType::GAS]; + if( ipWater >= 0 ) // if water is present + { + GEOSX_THROW_IF( m_waterRefPressure <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterRefPressureString(), + InputError ); + GEOSX_THROW_IF( m_waterFormationVolFactor <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterFormationVolumeFactorString(), + InputError ); + GEOSX_THROW_IF( m_waterCompressibility <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterCompressibilityString(), + InputError ); + GEOSX_THROW_IF( m_waterViscosity <= 0.0, + "DeadOilFluid: a strictly positive value must be provided for: " + << viewKeyStruct::waterViscosityString(), + InputError ); + } + else + { + GEOSX_THROW_IF( m_waterRefPressure > 0.0, + "DeadOilFluid: if water is absent, this keyword is not needed " + << viewKeyStruct::waterRefPressureString(), + InputError ); + GEOSX_THROW_IF( m_waterFormationVolFactor > 0.0, + "DeadOilFluid: if water is absent, this keyword is not needed " + << viewKeyStruct::waterFormationVolumeFactorString(), + InputError ); + GEOSX_THROW_IF( m_waterViscosity > 0.0, + "DeadOilFluid: if water is absent, this keyword is not needed " + << viewKeyStruct::waterViscosityString(), + InputError ); + GEOSX_THROW_IF( m_waterCompressibility > 0.0, + "DeadOilFluid: if water is absent, this keyword is not needed " + << viewKeyStruct::waterCompressibilityString(), + InputError ); + } + + localIndex const numExpectedTables = (ipGas >= 0) ? 2 : 1; + GEOSX_THROW_IF( m_formationVolFactorTableNames.size() != numExpectedTables, + "DeadOilFluid: one formation volume factor table must be provided for each hydrocarbon phase", + InputError ); + GEOSX_THROW_IF( m_viscosityTableNames.size() != numExpectedTables, + "DeadOilFluid: one viscosity table must be provided for each hydrocarbon phase", + InputError ); + + for( localIndex ip = 0; ip < numFluidPhases(); ++ip ) + { + if( m_phaseTypes[ip] == PhaseType::OIL || m_phaseTypes[ip] == PhaseType::GAS ) + { + m_hydrocarbonPhaseOrder.emplace_back( LvArray::integerConversion< integer >( ip ) ); + } + } + + FunctionManager const & functionManager = getGlobalState().getFunctionManager(); + for( localIndex iph = 0; iph < m_hydrocarbonPhaseOrder.size(); ++iph ) + { + GEOSX_THROW_IF( !functionManager.hasGroup( m_formationVolFactorTableNames[iph] ), + "DeadOilFluid: the formation volume factor table " << m_formationVolFactorTableNames[iph] + << " could not be found", + InputError ); + GEOSX_THROW_IF( !functionManager.hasGroup( m_viscosityTableNames[iph] ), + "DeadOilFluid: the viscosity table " << m_viscosityTableNames[iph] + << " could not be found", + InputError ); + } } void DeadOilFluid::postProcessInput() @@ -88,128 +326,33 @@ void DeadOilFluid::postProcessInput() for( localIndex ip = 0; ip < numFluidPhases(); ++ip ) { auto it = phaseDict.find( m_phaseNames[ip] ); - GEOSX_ERROR_IF( it == phaseDict.end(), "DeadOilFluid: phase not supported: " << m_phaseNames[ip] ); + GEOSX_THROW_IF( it == phaseDict.end(), + "DeadOilFluid: phase not supported: " << m_phaseNames[ip], + InputError ); integer const phaseIndex = it->second; - GEOSX_ERROR_IF( phaseIndex >= PhaseType::MAX_NUM_PHASES, "DeadOilFluid: invalid phase index " << phaseIndex ); + GEOSX_THROW_IF( phaseIndex >= PhaseType::MAX_NUM_PHASES, "DeadOilFluid: invalid phase index " << phaseIndex, + InputError ); m_phaseTypes[ip] = phaseIndex; m_phaseOrder[phaseIndex] = LvArray::integerConversion< integer >( ip ); } - GEOSX_ERROR_IF( m_tableFiles.size() != numFluidPhases(), - "The number of table files must be equal to the number of phases" ); - for( localIndex ip = 0; ip < numFluidPhases(); ++ip ) - { - string const filename = m_tableFiles[ip]; - - array1d< array1d< real64 > > pressureCoords; - pressureCoords.resize( 1 ); - array1d< real64 > formationVolFactor; - array1d< real64 > viscosity; - - std::ifstream is( filename ); - constexpr std::streamsize buf_size = 256; - char buf[buf_size]; - - localIndex counter = 0; - while( is.getline( buf, buf_size ) ) - { - string const str( buf ); - string_array const strs = Tokenize( str, " " ); - - if( strs.empty() ) - { - continue; - } - - if( strs[0].front() == '#' ) - { - continue; - } - if( streq( strs[0], "--" ) ) - { - continue; - } - - std::cout << "strs.size() = " << strs.size() << std::endl; - for( localIndex i = 0; i < strs.size(); ++i ) - { - std::cout << strs[i] << std::endl; - } - - GEOSX_ERROR_IF( (m_phaseTypes[ip] != PhaseType::WATER) - && strs.size() != 3, - "Three columns (pressure, formation volume factor, and viscosity) are expected for oil and gas" ); - GEOSX_ERROR_IF( (m_phaseTypes[ip] == PhaseType::WATER) - && strs.size() != 4, - "Four columns (pressure, formation volume factor, compressibility, and viscosity) are expected for water" ); - GEOSX_ERROR_IF( (m_phaseTypes[ip] == PhaseType::WATER) - && counter > 1, - "The water table must contain only one line" ); - - if( m_phaseTypes[ip] == PhaseType::WATER ) - { - m_waterRefPressure = stod( strs[0] ); - m_waterFormationVolFactor = stod( strs[1] ); - m_waterCompressibility = stod( strs[2] ); - m_waterViscosity = stod( strs[3] ); - - GEOSX_ERROR_IF( m_waterRefPressure <= 0.0, - "DeadOilFluid: a strictly positive value must be provided for: " - << viewKeyStruct::waterRefPressureString() ); - GEOSX_ERROR_IF( m_waterFormationVolFactor <= 0.0, - "DeadOilFluid: a strictly positive value must be provided for: " - << viewKeyStruct::waterFormationVolumeFactorString() ); - GEOSX_ERROR_IF( m_waterCompressibility <= 0.0, - "DeadOilFluid: a strictly positive value must be provided for: " - << viewKeyStruct::waterCompressibilityString() ); - GEOSX_ERROR_IF( m_waterViscosity <= 0.0, - "DeadOilFluid: a strictly positive value must be provided for: " - << viewKeyStruct::waterViscosityString() ); - } - else - { - pressureCoords[0].emplace_back( stod( strs[0] ) ); - formationVolFactor.emplace_back( stod( strs[1] ) ); - viscosity.emplace_back( stod( strs[2] ) ); - } - - counter++; - } + GEOSX_THROW_IF( m_surfacePhaseMassDensity.size() != m_phaseNames.size(), + "DeadOilFluid: the number of surfacePhaseMassDensities is inconsistent with the phase names", + InputError ); + GEOSX_THROW_IF( m_componentMolarWeight.size() != m_phaseNames.size(), + "DeadOilFluid: the number of componentMolarWeights is inconsistent with the phase names", + InputError ); - if( m_phaseTypes[ip] != PhaseType::WATER ) - { - m_hydrocarbonPhaseOrder.emplace_back( LvArray::integerConversion< integer >( ip ) ); - string const formationVolFactorTableName = (m_phaseTypes[ip] == PhaseType::OIL) ? "PVDO_Bo" : "PVDG_Bg"; - string const viscosityTableName = (m_phaseTypes[ip] == PhaseType::OIL) ? "PVDO_visco" : "PVDG_viscg"; - m_formationVolFactorTableNames.emplace_back( formationVolFactorTableName ); - m_viscosityTableNames.emplace_back( viscosityTableName ); - - GEOSX_ERROR_IF( pressureCoords[0].size() <= 1, - "DeadOilFluid: The oil and gas PVT tables must contain at least 2 values" ); - - FunctionManager & functionManager = getGlobalState().getFunctionManager(); - TableFunction & tablePVDX_B = - dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", formationVolFactorTableName ) ); - tablePVDX_B.setTableCoordinates( pressureCoords ); - tablePVDX_B.setTableValues( formationVolFactor ); - tablePVDX_B.reInitializeFunction(); - tablePVDX_B.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - TableFunction & tablePVDX_visc = - dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", viscosityTableName ) ); - tablePVDX_visc.setTableCoordinates( pressureCoords ); - tablePVDX_visc.setTableValues( viscosity ); - tablePVDX_visc.reInitializeFunction(); - tablePVDX_visc.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - } - is.close(); + // at this point, we make the distinction between the two input options + if( !m_tableFiles.empty() ) + { + readInputDataFromPVTFiles(); + } + else + { + useProvidedTableFunctions(); } - - // check the size of the additional parameters - GEOSX_ERROR_IF( m_surfacePhaseMassDensity.size() != m_phaseNames.size(), - "DeadOilFluid: the number of surfacePhaseMassDensities is inconsistent with the phase names" ); - GEOSX_ERROR_IF( m_componentMolarWeight.size() != m_phaseNames.size(), - "DeadOilFluid: the number of componentMolarWeights is inconsistent with the phase names" ); } void DeadOilFluid::initializePostSubGroups() @@ -222,8 +365,9 @@ void DeadOilFluid::createAllKernelWrappers() { FunctionManager const & functionManager = getGlobalState().getFunctionManager(); - GEOSX_ERROR_IF( m_hydrocarbonPhaseOrder.size() != 1 && m_hydrocarbonPhaseOrder.size() != 2, - "DeadOilFluid: the number of hydrocarbon phases should be equal to 1 (oil) or 2 (oil+gas)" ); + GEOSX_THROW_IF( m_hydrocarbonPhaseOrder.size() != 1 && m_hydrocarbonPhaseOrder.size() != 2, + "DeadOilFluid: the number of hydrocarbon phases should be equal to 1 (oil) or 2 (oil+gas)", + InputError ); if( m_formationVolFactorTables.size() == 0 && m_viscosityTables.size() == 0 ) { @@ -249,10 +393,18 @@ void DeadOilFluid::validateTable( TableFunction const & table ) const { array1d< real64 > const & property = table.getValues(); + GEOSX_THROW_IF( table.getInterpolationMethod() != TableFunction::InterpolationType::Linear, + "DeadOilFluid: the interpolation method in the table must be linear", + InputError ); + GEOSX_THROW_IF( property.size() <= 1, + "DeadOilFluid: each PVT table must contain at least 2 values", + InputError ); + for( localIndex i = 2; i < property.size(); ++i ) { - GEOSX_ERROR_IF( (property[i] - property[i-1]) * (property[i-1] - property[i-2]) < 0, - "DeadOilFluid: the values in each PVT table must monotone" ); + GEOSX_THROW_IF( (property[i] - property[i-1]) * (property[i-1] - property[i-2]) < 0, + "DeadOilFluid: the values in each PVT table must monotone", + InputError ); } } diff --git a/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp index c3ef1de6374..6c57fa6d3db 100644 --- a/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp @@ -206,7 +206,7 @@ class DeadOilFluidUpdate final : public MultiFluidBaseUpdate //GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE void computeDensities( real64 pressure, - arraySlice1d< real64, 0 > const & phaseMassDens ) const; + arraySlice1d< real64 > const & phaseMassDens ) const; //GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE @@ -218,7 +218,7 @@ class DeadOilFluidUpdate final : public MultiFluidBaseUpdate //GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE void computeViscosities( real64 pressure, - arraySlice1d< real64, 0 > const & phaseVisc ) const; + arraySlice1d< real64 > const & phaseVisc ) const; //GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE @@ -349,6 +349,39 @@ class DeadOilFluid : public MultiFluidBase private: + /** + * @brief Read all the PVT table provided by the user in Eclipse format + */ + void readInputDataFromPVTFiles(); + + /** + * @brief Use the TableFunctions provided by the user to get the PVT data + */ + void useProvidedTableFunctions(); + + /** + * @brief Open the file "filename", read the table inside, and fill "values" with the table values + * @param[in] filename the name of the file containing the table + * @param[inout] values an array of arrays containing the table values + */ + void readTable( string const & filename, + array1d< array1d< real64 > > & values ) const; + + /** + * @brief Fill the water data (formation vol factor, compressibility, etc) + * @param[in] tableValues the values in the water table + */ + void fillWaterData( array1d< array1d< real64 > > const & tableValues ); + + /** + * @brief Fill the hydrocarbon data (pressure, formation vol factor, viscosity) + * @param[in] ip the index of the phase + * @param[in] tableValues the values in the oil or gas table + */ + void fillHydrocarbonData( localIndex const ip, + array1d< array1d< real64 > > const & tableValues ); + + /// create all the table kernel wrappers void createAllKernelWrappers(); @@ -402,7 +435,7 @@ class DeadOilFluid : public MultiFluidBase //GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE void DeadOilFluidUpdate::computeDensities( real64 pressure, - arraySlice1d< real64, 0 > const & phaseMassDens ) const + arraySlice1d< real64 > const & phaseMassDens ) const { real64 fvf = 0.0; real64 derivative = 0.0; @@ -554,12 +587,12 @@ void DeadOilFluidUpdate::computeViscosities( real64 pressure, GEOSX_FORCE_INLINE void DeadOilFluidUpdate::compute( real64 pressure, real64 temperature, - arraySlice1d< real64 const, 0 > const & composition, - arraySlice1d< real64, 0 > const & phaseFrac, - arraySlice1d< real64, 0 > const & phaseDens, - arraySlice1d< real64, 0 > const & phaseMassDens, - arraySlice1d< real64, 0 > const & phaseVisc, - arraySlice2d< real64, 1 > const & phaseCompFrac, + arraySlice1d< real64 const > const & composition, + arraySlice1d< real64 > const & phaseFrac, + arraySlice1d< real64 > const & phaseDens, + arraySlice1d< real64 > const & phaseMassDens, + arraySlice1d< real64 > const & phaseVisc, + arraySlice2d< real64 > const & phaseCompFrac, real64 & totalDens ) const { GEOSX_UNUSED_VAR( temperature ); @@ -603,7 +636,7 @@ void DeadOilFluidUpdate::compute( real64 pressure, GEOSX_FORCE_INLINE void DeadOilFluidUpdate::compute( real64 pressure, real64 temperature, - arraySlice1d< real64 const, 0 > const & composition, + arraySlice1d< real64 const > const & composition, arraySlice1d< real64 > const & phaseFraction, arraySlice1d< real64 > const & dPhaseFraction_dPressure, arraySlice1d< real64 > const & dPhaseFraction_dTemperature, diff --git a/src/coreComponents/constitutive/fluid/MultiFluidPVTPackageWrapper.cpp b/src/coreComponents/constitutive/fluid/MultiFluidPVTPackageWrapper.cpp index ab1f7404ad9..3f5ca015e79 100644 --- a/src/coreComponents/constitutive/fluid/MultiFluidPVTPackageWrapper.cpp +++ b/src/coreComponents/constitutive/fluid/MultiFluidPVTPackageWrapper.cpp @@ -83,12 +83,12 @@ MultiFluidPVTPackageWrapper::deliverClone( string const & name, void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure, real64 temperature, - arraySlice1d< real64 const, 0 > const & composition, - arraySlice1d< real64, 0 > const & phaseFrac, - arraySlice1d< real64, 0 > const & phaseDens, - arraySlice1d< real64, 0 > const & phaseMassDens, - arraySlice1d< real64, 0 > const & phaseVisc, - arraySlice2d< real64, 1 > const & phaseCompFrac, + arraySlice1d< real64 const > const & composition, + arraySlice1d< real64 > const & phaseFrac, + arraySlice1d< real64 > const & phaseDens, + arraySlice1d< real64 > const & phaseMassDens, + arraySlice1d< real64 > const & phaseVisc, + arraySlice2d< real64 > const & phaseCompFrac, real64 & totalDens ) const { localIndex const NC = m_componentMolarWeight.size(); @@ -199,7 +199,7 @@ void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure, void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure, real64 temperature, - arraySlice1d< real64 const, 0 > const & composition, + arraySlice1d< real64 const > const & composition, arraySlice1d< real64 > const & phaseFraction, arraySlice1d< real64 > const & dPhaseFraction_dPressure, arraySlice1d< real64 > const & dPhaseFraction_dTemperature, @@ -223,7 +223,7 @@ void MultiFluidPVTPackageWrapperUpdate::compute( real64 pressure, real64 & totalDensity, real64 & dTotalDensity_dPressure, real64 & dTotalDensity_dTemperature, - arraySlice1d< real64, 0 > const & dTotalDensity_dGlobalCompFraction ) const + arraySlice1d< real64 > const & dTotalDensity_dGlobalCompFraction ) const { // 0. make shortcut structs to avoid long names (TODO maybe remove) CompositionalVarContainer< 1 > phaseFrac{ diff --git a/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp b/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp index 4f71c7327ec..547b1301b83 100644 --- a/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp +++ b/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp @@ -210,6 +210,7 @@ void testNumericalDerivatives( MultiFluidBase & fluid, "Pres", phases, components ); + std::cout << phaseDens.value << std::endl; } // update temperature and check derivatives @@ -425,6 +426,113 @@ MultiFluidBase & makeDeadOilFluid( string const & name, Group * parent ) return fluid; } +MultiFluidBase & makeDeadOilFluidFromTable( string const & name, Group * parent ) +{ + FunctionManager & functionManager = getGlobalState().getFunctionManager(); + + // 1) First, define the tables (PVDO, PVDG) + + // 1D table with linear interpolation + localIndex const NaxisPVDO = 7; + localIndex 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] = 2000000; valuesPVDO_Bo[0] = 1.02; valuesPVDO_visc[0] = 0.000975; + coordinatesPVDO[0][1] = 5000000; valuesPVDO_Bo[1] = 1.03; valuesPVDO_visc[1] = 0.00091; + coordinatesPVDO[0][2] = 10000000; valuesPVDO_Bo[2] = 1.04; valuesPVDO_visc[2] = 0.00083; + coordinatesPVDO[0][3] = 20000000; valuesPVDO_Bo[3] = 1.05; valuesPVDO_visc[3] = 0.000695; + coordinatesPVDO[0][4] = 30000000; valuesPVDO_Bo[4] = 1.07; valuesPVDO_visc[4] = 0.000594; + coordinatesPVDO[0][5] = 40000000; valuesPVDO_Bo[5] = 1.08; valuesPVDO_visc[5] = 0.00051; + coordinatesPVDO[0][6] = 50000000; valuesPVDO_Bo[6] = 1.09; valuesPVDO_visc[6] = 0.000449; + + 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 ); + tablePVDO_Bo.setTableValues( valuesPVDO_Bo ); + tablePVDO_Bo.reInitializeFunction(); + tablePVDO_Bo.setInterpolationMethod( TableFunction::InterpolationType::Linear ); + + TableFunction & tablePVDO_visc = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDO_visc" ) ); + tablePVDO_visc.setTableCoordinates( coordinatesPVDO ); + tablePVDO_visc.setTableValues( valuesPVDO_visc ); + tablePVDO_visc.reInitializeFunction(); + tablePVDO_visc.setInterpolationMethod( TableFunction::InterpolationType::Linear ); + + TableFunction & tablePVDG_Bg = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDG_Bg" ) ); + tablePVDG_Bg.setTableCoordinates( coordinatesPVDG ); + tablePVDG_Bg.setTableValues( valuesPVDG_Bg ); + tablePVDG_Bg.reInitializeFunction(); + tablePVDG_Bg.setInterpolationMethod( TableFunction::InterpolationType::Linear ); + + TableFunction & tablePVDG_visc = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "PVDG_visc" ) ); + tablePVDG_visc.setTableCoordinates( coordinatesPVDG ); + tablePVDG_visc.setTableValues( valuesPVDG_visc ); + 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 ); @@ -558,6 +666,37 @@ TEST_F( DeadOilFluidTest, numericalDerivativesMass ) } } +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( localIndex i = 0; i < 3; ++i ) + { + testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, relTol ); + } +} int main( int argc, char * * argv ) { diff --git a/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst b/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst index 93045746e8f..c84a4c0613d 100644 --- a/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst +++ b/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst @@ -1,14 +1,24 @@ -==================== ============ ======== ============================================= -Name Type Default Description -==================== ============ ======== ============================================= -componentMolarWeight real64_array required Component molar weights -componentNames string_array {} List of component names -name string required A name is required for any non-unique nodes -phaseNames string_array required List of fluid phases -surfaceDensities real64_array required List of surface mass densities for each phase -tableFiles path_array required List of filenames with input PVT tables -==================== ============ ======== ============================================= +======================================= ============ ======== ===================================================================================================================================================================================================================================================================================================== +Name Type Default Description +======================================= ============ ======== ===================================================================================================================================================================================================================================================================================================== +componentMolarWeight real64_array required Component molar weights +componentNames string_array {} List of component names +hydrocarbonFormationVolFactorTableNames string_array {} | List of formation volume factor TableFunction names from the Functions block. + | The user must provide one TableFunction per hydrocarbon phase, in the order provided in "phaseNames". + | For instance, if "oil" is before "gas" in "phaseNames", the table order should be: oilTableName, gasTableName +hydrocarbonViscosityTableNames string_array {} | List of viscosity TableFunction names from the Functions block. + | The user must provide one TableFunction per hydrocarbon phase, in the order provided in "phaseNames". + | For instance, if "oil" is before "gas" in "phaseNames", the table order should be: oilTableName, gasTableName +name string required A name is required for any non-unique nodes +phaseNames string_array required List of fluid phases +surfaceDensities real64_array required List of surface mass densities for each phase +tableFiles path_array {} List of filenames with input PVT tables (one per phase) +waterCompressibility real64 0 Water compressibility +waterFormationVolumeFactor real64 0 Water formation volume factor +waterReferencePressure real64 0 Water reference pressure +waterViscosity real64 0 Water viscosity +======================================= ============ ======== ===================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst b/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst index 7d677303016..5b3d3c27bb3 100644 --- a/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst +++ b/src/coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst @@ -1,39 +1,33 @@ -======================================= ============================================================================================== ========================== -Name Type Description -======================================= ============================================================================================== ========================== -dPhaseCompFraction_dGlobalCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l, 3l, 4l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseCompFraction_dPressure LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseCompFraction_dTemperature LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseDensity_dPressure real64_array3d (no description available) -dPhaseDensity_dTemperature real64_array3d (no description available) -dPhaseFraction_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseFraction_dPressure real64_array3d (no description available) -dPhaseFraction_dTemperature real64_array3d (no description available) -dPhaseMassDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseMassDensity_dPressure real64_array3d (no description available) -dPhaseMassDensity_dTemperature real64_array3d (no description available) -dPhaseViscosity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -dPhaseViscosity_dPressure real64_array3d (no description available) -dPhaseViscosity_dTemperature real64_array3d (no description available) -dTotalDensity_dGlobalCompFraction real64_array3d (no description available) -dTotalDensity_dPressure real64_array2d (no description available) -dTotalDensity_dTemperature real64_array2d (no description available) -hydrocarbonFormationVolFactorTableNames string_array (no description available) -hydrocarbonViscosityTableNames string_array (no description available) -phaseCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) -phaseDensity real64_array3d (no description available) -phaseFraction real64_array3d (no description available) -phaseMassDensity real64_array3d (no description available) -phaseViscosity real64_array3d (no description available) -totalDensity real64_array2d (no description available) -useMass integer (no description available) -waterCompressibility real64 (no description available) -waterFormationVolumeFactor real64 (no description available) -waterReferencePressure real64 (no description available) -waterViscosity real64 (no description available) -======================================= ============================================================================================== ========================== +====================================== ============================================================================================== ========================== +Name Type Description +====================================== ============================================================================================== ========================== +dPhaseCompFraction_dGlobalCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l, 3l, 4l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseCompFraction_dPressure LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseCompFraction_dTemperature LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseDensity_dPressure real64_array3d (no description available) +dPhaseDensity_dTemperature real64_array3d (no description available) +dPhaseFraction_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseFraction_dPressure real64_array3d (no description available) +dPhaseFraction_dTemperature real64_array3d (no description available) +dPhaseMassDensity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseMassDensity_dPressure real64_array3d (no description available) +dPhaseMassDensity_dTemperature real64_array3d (no description available) +dPhaseViscosity_dGlobalCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +dPhaseViscosity_dPressure real64_array3d (no description available) +dPhaseViscosity_dTemperature real64_array3d (no description available) +dTotalDensity_dGlobalCompFraction real64_array3d (no description available) +dTotalDensity_dPressure real64_array2d (no description available) +dTotalDensity_dTemperature real64_array2d (no description available) +phaseCompFraction LvArray_Array< double, 4, camp_int_seq< long, 0l, 1l, 2l, 3l >, long, LvArray_ChaiBuffer > (no description available) +phaseDensity real64_array3d (no description available) +phaseFraction real64_array3d (no description available) +phaseMassDensity real64_array3d (no description available) +phaseViscosity real64_array3d (no description available) +totalDensity real64_array2d (no description available) +useMass integer (no description available) +====================================== ============================================================================================== ========================== diff --git a/src/coreComponents/fileIO/schema/schema.xsd b/src/coreComponents/fileIO/schema/schema.xsd index b065b8101a2..5b023be541d 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd +++ b/src/coreComponents/fileIO/schema/schema.xsd @@ -1847,12 +1847,28 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + + + - - + + + + + + + + + + diff --git a/src/coreComponents/fileIO/schema/schema.xsd.other b/src/coreComponents/fileIO/schema/schema.xsd.other index ad68c4950a2..a8e64a577ee 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd.other +++ b/src/coreComponents/fileIO/schema/schema.xsd.other @@ -959,10 +959,6 @@ - - - - @@ -977,14 +973,6 @@ - - - - - - - - diff --git a/src/coreComponents/managers/FieldSpecification/FieldSpecificationOps.hpp b/src/coreComponents/managers/FieldSpecification/FieldSpecificationOps.hpp index 9b8efc3a207..93ca7156d52 100644 --- a/src/coreComponents/managers/FieldSpecification/FieldSpecificationOps.hpp +++ b/src/coreComponents/managers/FieldSpecification/FieldSpecificationOps.hpp @@ -478,10 +478,22 @@ struct FieldSpecificationEqual : public FieldSpecificationOp< OpEqual > localIndex const numEntries = matrix.numNonZeros( localRow ); real64 diagonal = 0; + real64 const minDiagonal = 1e-15; for( localIndex j = 0; j < numEntries; ++j ) { if( columns[ j ] == dof ) - { diagonal = entries[ j ]; } + { + // check that the entry is large enough to enforce the boundary condition + if( entries[j] >= 0 && entries[j] < minDiagonal ) + { + entries[ j ] = minDiagonal; + } + else if( entries[j] < 0 && entries[j] > -minDiagonal ) + { + entries[ j ] = -minDiagonal; + } + diagonal = entries[j]; + } else { entries[ j ] = 0; } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/B_g_pvdg.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/B_g_pvdg.txt new file mode 100644 index 00000000000..c448c139d51 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/B_g_pvdg.txt @@ -0,0 +1,13 @@ +0.04234 +0.02046 +0.01328 +0.00977 +0.00773 +0.006426 +0.005541 +0.004919 +0.004471 +0.004194 +0.004031 +0.00391 +0.003868 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/B_o_pvdo.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/B_o_pvdo.txt new file mode 100644 index 00000000000..28d6dc2ca67 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/B_o_pvdo.txt @@ -0,0 +1,7 @@ +1.02 +1.03 +1.04 +1.05 +1.07 +1.08 +1.09 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_baker_1d.xml b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_baker_1d.xml index a89bb45e803..fa21ffc3238 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_baker_1d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_baker_1d.xml @@ -6,7 +6,7 @@ name="compflow" logLevel="1" discretization="fluidTPFA" - targetRegions="{ Region1 }" + targetRegions="{ region }" fluidNames="{ fluid }" solidNames="{ rock }" relPermNames="{ relperm }" @@ -27,7 +27,7 @@ @@ -123,7 +123,7 @@ component="0" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -132,7 +132,7 @@ component="1" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -141,7 +141,7 @@ component="2" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -149,7 +149,7 @@ name="referencePorosity" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="referencePorosity" scale="0.05"/> @@ -158,16 +158,16 @@ name="initialPressure" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="pressure" - scale="5e6"/> + scale="7.5e6"/> @@ -176,7 +176,7 @@ name="initialComposition_gas" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="globalCompFraction" component="1" scale="0.399"/> @@ -185,7 +185,7 @@ name="initialComposition_water" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="globalCompFraction" component="2" scale="0.001"/> @@ -193,7 +193,7 @@ @@ -202,7 +202,7 @@ @@ -210,7 +210,7 @@ @@ -218,24 +218,24 @@ - + @@ -243,7 +243,7 @@ @@ -251,7 +251,7 @@ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_corey_1d.xml b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_corey_1d.xml index 8629a057187..96491444315 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_corey_1d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_corey_1d.xml @@ -4,14 +4,14 @@ + targetRegions="{ region }"> @@ -27,7 +27,7 @@ @@ -121,7 +121,7 @@ component="0" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -130,7 +130,7 @@ component="1" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -139,7 +139,7 @@ component="2" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -147,25 +147,25 @@ name="referencePorosity" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="referencePorosity" scale="0.2"/> - + + scale="7.5e6"/> @@ -174,7 +174,7 @@ name="initialComposition_gas" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="globalCompFraction" component="1" scale="0.399"/> @@ -183,7 +183,7 @@ name="initialComposition_water" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="globalCompFraction" component="2" scale="0.001"/> @@ -191,7 +191,7 @@ @@ -200,7 +200,7 @@ @@ -208,7 +208,7 @@ @@ -216,24 +216,24 @@ - + @@ -241,7 +241,7 @@ @@ -249,7 +249,7 @@ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_staircase_3d.xml b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_staircase_3d.xml index 1401d7ad58e..773cd1551a6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_staircase_3d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_staircase_3d.xml @@ -30,7 +30,7 @@ + phaseNames="{ gas, water, oil }" + surfaceDensities="{ 0.9907, 1022.0, 800.0 }" + componentMolarWeight="{ 16e-3, 18e-3, 114e-3 }" + hydrocarbonFormationVolFactorTableNames="{ B_g_table, B_o_table }" + hydrocarbonViscosityTableNames="{ visc_g_table, visc_o_table }" + waterReferencePressure="30600000.1" + waterFormationVolumeFactor="1.03" + waterCompressibility="0.00000000041" + waterViscosity="0.0003"/> + scale="1.0e-16"/> + scale="1.0e-16"/> + scale="1.0e-16"/> - + - + - + @@ -288,8 +293,33 @@ component="2" scale="0.001"/> - + + + + + + + + + + + scale="1.0e-16"/> + scale="1.0e-16"/> + scale="1.0e-16"/> - + - + diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pres_pvdg.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pres_pvdg.txt new file mode 100644 index 00000000000..3749d67af82 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pres_pvdg.txt @@ -0,0 +1,13 @@ +3000000 +6000000 +9000000 +12000000 +15000000 +18000000 +21000000 +24000000 +27000000 +29500000 +31000000 +33000000 +53000000 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pres_pvdo.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pres_pvdo.txt new file mode 100644 index 00000000000..c962957b8f5 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pres_pvdo.txt @@ -0,0 +1,7 @@ +2000000 +5000000 +10000000 +20000000 +30000000 +40000000 +50000000.7 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdg.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdg.txt index b25fe86e44e..811a899bf10 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdg.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/pvdg.txt @@ -1,14 +1,14 @@ # Pg(Pa) Bg(m3/sm3) Visc(Pa.s) -3000000 0.04234 0.00001344 -6000000 0.02046 0.0000142 -9000000 0.01328 0.00001526 -12000000 0.00977 0.0000166 -15000000 0.00773 0.00001818 -18000000 0.006426 0.00001994 -21000000 0.005541 0.00002181 -24000000 0.004919 0.0000237 -27000000 0.004471 0.00002559 -29500000 0.004194 0.00002714 -31000000 0.004031 0.00002806 -33000000 0.00391 0.00002832 -53000000 0.003868 0.00002935 +3000000 0.04234 0.00005344 +6000000 0.02046 0.0000542 +9000000 0.01328 0.00005526 +12000000 0.00977 0.0000566 +15000000 0.00773 0.00005818 +18000000 0.006426 0.00005994 +21000000 0.005541 0.00006181 +24000000 0.004919 0.0000637 +27000000 0.004471 0.00006559 +29500000 0.004194 0.00006714 +31000000 0.004031 0.00006806 +33000000 0.00391 0.00006832 +53000000 0.003868 0.00006935 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/visc_pvdg.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/visc_pvdg.txt new file mode 100644 index 00000000000..ced86144c7c --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/visc_pvdg.txt @@ -0,0 +1,13 @@ +0.00005344 +0.0000542 +0.00005526 +0.0000566 +0.00005818 +0.00005994 +0.00006181 +0.0000637 +0.00006559 +0.00006714 +0.00006806 +0.00006832 +0.00006935 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/visc_pvdo.txt b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/visc_pvdo.txt new file mode 100644 index 00000000000..1a65ee3a6f7 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/visc_pvdo.txt @@ -0,0 +1,7 @@ +0.000975 +0.00091 +0.00083 +0.000695 +0.000594 +0.00051 +0.000449 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvdo.txt b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvdo.txt index b06ef2578f4..a40ff8fa545 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvdo.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvdo.txt @@ -1,9 +1,9 @@ # P[Pa] Bo[m3/sm3] Visc(Pa.s) -101325 1.00 0.001 -2000000 1.02 0.000975 -5000000 1.03 0.00091 -10000000 1.04 0.00083 -20000000 1.05 0.000695 -30000000 1.07 0.000594 -40000000 1.08 0.00051 -50000000.7 1.09 0.000449 +101325 1.2 0.001 +2000000 1.2 0.000975 +5000000 1.181 0.00091 +10000000 1.166 0.00083 +20000000 1.143 0.000695 +30000000 1.13 0.000594 +40000000 1.101 0.00051 +50000000.7 1.09 0.000449 From f7e6deece7c64a569c656519b23f59d4cea990b7 Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Tue, 9 Mar 2021 19:27:12 -0800 Subject: [PATCH 6/6] updated integratedTests submodule --- integratedTests | 2 +- .../constitutive/fluid/DeadOilFluid.cpp | 4 +- .../interfaces/hypre/HyprePreconditioner.cpp | 25 +++++--- .../meshUtilities/PerforationData.hpp | 4 -- .../fluidFlow/benchmarks/Egg/dead_oil_egg.xml | 3 +- .../deadoil_3ph_corey_1d.xml | 2 +- .../wells/CompositionalMultiphaseWell.cpp | 37 ++++++----- .../wells/CompositionalMultiphaseWell.hpp | 7 +- .../CompositionalMultiphaseWellKernels.hpp | 64 ++++++++++++++++--- .../compositionalMultiphaseWell/pvdg.txt | 28 ++++---- .../compositionalMultiphaseWell/pvtw.txt | 2 +- .../FieldCaseCo2InjTutorial.xml | 22 ++++--- 12 files changed, 128 insertions(+), 72 deletions(-) diff --git a/integratedTests b/integratedTests index 45858a1ba21..d977a52535e 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit 45858a1ba211ff195fcfaa48f51357c95c544997 +Subproject commit d977a52535ea595195cb6350694b77ad910a236e diff --git a/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp b/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp index de902bae415..276c8f5afba 100644 --- a/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp @@ -104,7 +104,7 @@ void DeadOilFluid::readTable( string const & filename, while( is.getline( buf, buf_size ) ) { string const str( buf ); - string_array const strs = Tokenize( str, " " ); + string_array const strs = stringutilities::tokenize( str, " " ); if( strs.empty() ) { @@ -114,7 +114,7 @@ void DeadOilFluid::readTable( string const & filename, { continue; } - if( streq( strs[0], "--" ) ) + if( strs[0] == "--" ) { continue; } diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp index e6e4e15e5e0..d2f9b191c2b 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp @@ -584,14 +584,19 @@ void HyprePreconditioner::createMGR( DofManager const * const dofManager ) mgr_nlevels = 3; /* options for solvers at each level */ - HYPRE_Int mgr_gsmooth_type = 16; // ILU(0) - HYPRE_Int mgr_num_gsmooth_sweeps = 1; + HYPRE_Int mgr_gsmooth_type = 16; + HYPRE_Int mgr_num_gsmooth_sweeps = 0; mgr_level_interp_type.resize( mgr_nlevels ); mgr_level_interp_type[0] = 2; mgr_level_interp_type[1] = 2; mgr_level_interp_type[2] = 2; + mgr_coarse_grid_method.resize( mgr_nlevels ); + mgr_coarse_grid_method[0] = 1; + mgr_coarse_grid_method[1] = 1; + mgr_coarse_grid_method[2] = 0; + mgr_level_frelax_method.resize( mgr_nlevels ); mgr_level_frelax_method[0] = 0; // Jacobi mgr_level_frelax_method[1] = 0; // Jacobi @@ -635,11 +640,8 @@ void HyprePreconditioner::createMGR( DofManager const * const dofManager ) mgr_cindexes[iLevel] = lv_cindexes[iLevel].data(); } - GEOSX_LAI_CHECK_ERROR( HYPRE_ILUCreate( &aux_precond ) ); - GEOSX_LAI_CHECK_ERROR( HYPRE_ILUSetType( aux_precond, 0 ) ); // Block Jacobi - ILU - GEOSX_LAI_CHECK_ERROR( HYPRE_ILUSetLevelOfFill( aux_precond, 0 ) ); - GEOSX_LAI_CHECK_ERROR( HYPRE_ILUSetMaxIter( aux_precond, 1 ) ); - GEOSX_LAI_CHECK_ERROR( HYPRE_ILUSetTol( aux_precond, 0.0 ) ); + + GEOSX_LAI_CHECK_ERROR( HYPRE_MGRDirectSolverCreate( &aux_precond ) ); GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetCpointsByPointMarkerArray( m_precond, mgr_bsize, mgr_nlevels, mgr_num_cindexes.data(), @@ -648,17 +650,20 @@ void HyprePreconditioner::createMGR( DofManager const * const dofManager ) GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetLevelFRelaxMethod( m_precond, mgr_level_frelax_method.data() ) ); GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetNonCpointsToFpoints( m_precond, 1 )); + GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetTruncateCoarseGridThreshold( m_precond, 1e-14 )); + GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetPMaxElmts( m_precond, 15 )); GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetLevelInterpType( m_precond, mgr_level_interp_type.data() ) ); + GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetCoarseGridMethod( m_precond, mgr_coarse_grid_method.data() ) ); GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetGlobalsmoothType( m_precond, mgr_gsmooth_type ) ); GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetMaxGlobalsmoothIters( m_precond, mgr_num_gsmooth_sweeps ) ); GEOSX_LAI_CHECK_ERROR( HYPRE_MGRSetCoarseSolver( m_precond, - (HYPRE_PtrToParSolverFcn)HYPRE_ILUSolve, - (HYPRE_PtrToParSolverFcn)HYPRE_ILUSetup, + (HYPRE_PtrToParSolverFcn)HYPRE_MGRDirectSolverSolve, + (HYPRE_PtrToParSolverFcn)HYPRE_MGRDirectSolverSetup, aux_precond ) ); - m_functions->aux_destroy = HYPRE_ILUDestroy; + m_functions->aux_destroy = HYPRE_MGRDirectSolverDestroy; } else if( m_parameters.mgr.strategy == "CompositionalMultiphaseHybridFVM" ) { diff --git a/src/coreComponents/meshUtilities/PerforationData.hpp b/src/coreComponents/meshUtilities/PerforationData.hpp index 732187db770..753baed9266 100644 --- a/src/coreComponents/meshUtilities/PerforationData.hpp +++ b/src/coreComponents/meshUtilities/PerforationData.hpp @@ -234,8 +234,6 @@ class PerforationData : public ObjectManagerBase static constexpr char const * locationString() { return "location"; } /// @return String key for the well transmissibility static constexpr char const * wellTransmissibilityString() { return "wellTransmissibility"; } - /// @return String key for the perforation direction - static constexpr char const * directionString() { return "direction"; } /// ViewKey for the global number of perforations dataRepository::ViewKey numPerforationsGlobal = { numPerforationsGlobalString() }; @@ -251,8 +249,6 @@ class PerforationData : public ObjectManagerBase dataRepository::ViewKey location = { locationString() }; /// ViewKey for the well transmissibility dataRepository::ViewKey wellTransmissibility = { wellTransmissibilityString() }; - /// ViewKey for the perf direction - dataRepository::ViewKey direction = { directionString() }; } /// ViewKey struct for the PerforationData class diff --git a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/dead_oil_egg.xml b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/dead_oil_egg.xml index a4b8c423e10..476e0533c96 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/dead_oil_egg.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/benchmarks/Egg/dead_oil_egg.xml @@ -12,7 +12,7 @@ targetRegions="{ reservoir, wellRegion1, wellRegion2, wellRegion3, wellRegion4, wellRegion5, wellRegion6, wellRegion7, wellRegion8, wellRegion9, wellRegion10, wellRegion11, wellRegion12 }"> ( keys::dGlobalCompFraction_dGlobalCompDensityString() ); m_dResCompFrac_dCompDens.setName( getName() + "/accessors/" + keys::dGlobalCompFraction_dGlobalCompDensityString() ); - m_resPhaseMob.clear(); - m_resPhaseMob = elemManager.constructArrayViewAccessor< real64, 2 >( keys::phaseMobilityString() ); - m_resPhaseMob.setName( getName() + "/accessors/" + keys::phaseMobilityString() ); - - m_dResPhaseMob_dPres.clear(); - m_dResPhaseMob_dPres = elemManager.constructArrayViewAccessor< real64, 2 >( keys::dPhaseMobility_dPressureString() ); - m_dResPhaseMob_dPres.setName( getName() + "/accessors/" + keys::dPhaseMobility_dPressureString() ); - - m_dResPhaseMob_dCompDens.clear(); - m_dResPhaseMob_dCompDens = elemManager.constructArrayViewAccessor< real64, 3 >( keys::dPhaseMobility_dGlobalCompDensityString() ); - m_dResPhaseMob_dCompDens.setName( getName() + "/accessors/" + keys::dPhaseMobility_dGlobalCompDensityString() ); - m_resPhaseVolFrac.clear(); m_resPhaseVolFrac = elemManager.constructArrayViewAccessor< real64, 2 >( keys::phaseVolumeFractionString() ); m_resPhaseVolFrac.setName( getName() + "/accessors/" + keys::phaseVolumeFractionString() ); @@ -1447,6 +1435,25 @@ void CompositionalMultiphaseWell::resetViews( DomainPartition & domain ) { using keys = MultiFluidBase::viewKeyStruct; + m_resPhaseDens.clear(); + m_resPhaseDens = elemManager.constructMaterialArrayViewAccessor< real64, 3 >( keys::phaseDensityString(), + flowSolver.targetRegionNames(), + flowSolver.fluidModelNames() ); + m_resPhaseDens.setName( getName() + "/accessors/" + keys::phaseDensityString() ); + + m_dResPhaseDens_dPres.clear(); + m_dResPhaseDens_dPres = elemManager.constructMaterialArrayViewAccessor< real64, 3 >( keys::dPhaseDensity_dPressureString(), + flowSolver.targetRegionNames(), + flowSolver.fluidModelNames() ); + m_dResPhaseDens_dPres.setName( getName() + "/accessors/" + keys::dPhaseDensity_dPressureString() ); + + m_dResPhaseDens_dComp.clear(); + m_dResPhaseDens_dComp = elemManager.constructMaterialArrayViewAccessor< real64, 4 >( keys::dPhaseDensity_dGlobalCompFractionString(), + flowSolver.targetRegionNames(), + flowSolver.fluidModelNames() ); + m_dResPhaseDens_dComp.setName( getName() + "/accessors/" + keys::dPhaseDensity_dGlobalCompFractionString() ); + + m_resPhaseMassDens.clear(); m_resPhaseMassDens = elemManager.constructMaterialArrayViewAccessor< real64, 3 >( keys::phaseMassDensityString(), flowSolver.targetRegionNames(), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index 6c9c8d817b8..4b8f33c43d7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -409,13 +409,12 @@ class CompositionalMultiphaseWell : public WellSolverBase ElementRegionManager::ElementViewAccessor< arrayView3d< real64 const > > m_dResCompFrac_dCompDens; - ElementRegionManager::ElementViewAccessor< arrayView2d< real64 const > > m_resPhaseMob; - ElementRegionManager::ElementViewAccessor< arrayView2d< real64 const > > m_dResPhaseMob_dPres; - ElementRegionManager::ElementViewAccessor< arrayView3d< real64 const > > m_dResPhaseMob_dCompDens; - /// views into reservoir material fields ElementRegionManager::ElementViewAccessor< arrayView3d< real64 const > > m_resPhaseDens; + ElementRegionManager::ElementViewAccessor< arrayView3d< real64 const > > m_dResPhaseDens_dPres; + ElementRegionManager::ElementViewAccessor< arrayView4d< real64 const > > m_dResPhaseDens_dComp; + ElementRegionManager::ElementViewAccessor< arrayView3d< real64 const > > m_resPhaseMassDens; ElementRegionManager::ElementViewAccessor< arrayView3d< real64 const > > m_resPhaseVisc; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index 3154b8e693c..16d0fc5cb6c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -676,12 +676,12 @@ struct PerforationKernel localIndex const numPhases, ElementViewConst< arrayView1d< real64 const > > const & resPres, ElementViewConst< arrayView1d< real64 const > > const & dResPres, - ElementViewConst< arrayView2d< real64 const > > const & resPhaseMob, - ElementViewConst< arrayView2d< real64 const > > const & dResPhaseMob_dPres, - ElementViewConst< arrayView3d< real64 const > > const & dResPhaseMob_dComp, ElementViewConst< arrayView2d< real64 const > > const & dResPhaseVolFrac_dPres, ElementViewConst< arrayView3d< real64 const > > const & dResPhaseVolFrac_dComp, ElementViewConst< arrayView3d< real64 const > > const & dResCompFrac_dCompDens, + ElementViewConst< arrayView3d< real64 const > > const & resPhaseDens, + ElementViewConst< arrayView3d< real64 const > > const & dResPhaseDens_dPres, + ElementViewConst< arrayView4d< real64 const > > const & dResPhaseDens_dComp, ElementViewConst< arrayView3d< real64 const > > const & resPhaseVisc, ElementViewConst< arrayView3d< real64 const > > const & dResPhaseVisc_dPres, ElementViewConst< arrayView4d< real64 const > > const & dResPhaseVisc_dComp, @@ -737,8 +737,10 @@ struct PerforationKernel stackArray2d< real64, 2 * maxNumComp > dPhaseCompFrac_dP( 2, NC ); stackArray3d< real64, 2 * maxNumComp * maxNumComp > dPhaseCompFrac_dC( 2, NC, NC ); + stackArray1d< real64, maxNumComp > dDens_dC( NC ); stackArray1d< real64, maxNumComp > dVisc_dC( NC ); stackArray1d< real64, maxNumComp > dRelPerm_dC( NC ); + stackArray1d< real64, maxNumComp > dMob_dC( NC ); real64 dPotDiff_dP[ 2 ] = { 0.0 }; stackArray2d< real64, 2 * maxNumComp > dPotDiff_dC( 2, NC ); @@ -824,24 +826,66 @@ struct PerforationKernel for( localIndex ip = 0; ip < NP; ++ip ) { + // density + real64 const resDens = resPhaseDens[er][esr][ei][0][ip]; + real64 const dResDens_dP = dResPhaseDens_dPres[er][esr][ei][0][ip]; + applyChainRule( NC, dResCompFrac_dCompDens[er][esr][ei], + dResPhaseDens_dComp[er][esr][ei][0][ip], + dDens_dC ); + + // viscosity + real64 const resVisc = resPhaseVisc[er][esr][ei][0][ip]; + real64 const dResVisc_dP = dResPhaseVisc_dPres[er][esr][ei][0][ip]; + applyChainRule( NC, dResCompFrac_dCompDens[er][esr][ei], + dResPhaseVisc_dComp[er][esr][ei][0][ip], + dVisc_dC ); + + // relative permeability + real64 const resRelPerm = resPhaseRelPerm[er][esr][ei][0][ip]; + real64 dResRelPerm_dP = 0.0; + for( localIndex jc = 0; jc < NC; ++jc ) + { + dRelPerm_dC[jc] = 0; + } + + for( localIndex jp = 0; jp < NP; ++jp ) + { + real64 const dResRelPerm_dS = dResPhaseRelPerm_dPhaseVolFrac[er][esr][ei][0][ip][jp]; + dResRelPerm_dP += dResRelPerm_dS * dResPhaseVolFrac_dPres[er][esr][ei][jp]; + + for( localIndex jc = 0; jc < NC; ++jc ) + { + dRelPerm_dC[jc] += dResRelPerm_dS * dResPhaseVolFrac_dComp[er][esr][ei][jp][jc]; + } + } + + // compute the reservoir phase mobility, including phase density + real64 const resPhaseMob = resDens * resRelPerm / resVisc; + real64 const dResPhaseMob_dPres = dResRelPerm_dP * resDens / resVisc + + resPhaseMob * (dResDens_dP / resDens - dResVisc_dP / resVisc); + for( localIndex jc = 0; jc < NC; ++jc ) + { + dMob_dC[jc] = dRelPerm_dC[jc] * resDens / resVisc + + resPhaseMob * (dDens_dC[jc] / resDens - dVisc_dC[jc] / resVisc); + } + // compute the phase flux and derivatives using upstream cell mobility - flux = resPhaseMob[er][esr][ei][ip] * potDiff; + flux = resPhaseMob * potDiff; dFlux_dP[CompositionalMultiphaseWell::SubRegionTag::RES] = - dResPhaseMob_dPres[er][esr][ei][ip] * potDiff - + resPhaseMob[er][esr][ei][ip] * dPotDiff_dP[CompositionalMultiphaseWell::SubRegionTag::RES]; + dResPhaseMob_dPres * potDiff + resPhaseMob * dPotDiff_dP[CompositionalMultiphaseWell::SubRegionTag::RES]; dFlux_dP[CompositionalMultiphaseWell::SubRegionTag::WELL] = - resPhaseMob[er][esr][ei][ip] * dPotDiff_dP[CompositionalMultiphaseWell::SubRegionTag::WELL]; + resPhaseMob * dPotDiff_dP[CompositionalMultiphaseWell::SubRegionTag::WELL]; for( localIndex ic = 0; ic < NC; ++ic ) { dFlux_dC[CompositionalMultiphaseWell::SubRegionTag::RES][ic] = - dResPhaseMob_dComp[er][esr][ei][ip][ic] * potDiff - + resPhaseMob[er][esr][ei][ip] * dPotDiff_dC[CompositionalMultiphaseWell::SubRegionTag::RES][ic]; + dMob_dC[ic] * potDiff + + resPhaseMob * dPotDiff_dC[CompositionalMultiphaseWell::SubRegionTag::RES][ic]; dFlux_dC[CompositionalMultiphaseWell::SubRegionTag::WELL][ic] = - resPhaseMob[er][esr][ei][ip] * dPotDiff_dC[CompositionalMultiphaseWell::SubRegionTag::WELL][ic]; + resPhaseMob * dPotDiff_dC[CompositionalMultiphaseWell::SubRegionTag::WELL][ic]; } // increment component fluxes diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvdg.txt b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvdg.txt index 513bdab9ab5..9480e70b5a2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvdg.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvdg.txt @@ -1,15 +1,15 @@ # Pg(Pa) Bg(m3/sm3) Visc(Pa.s) -101325 1.0 0.00001200 -3000000 0.04234 0.00001344 -6000000 0.02046 0.0000142 -9000000 0.01328 0.00001526 -12000000 0.00977 0.0000166 -15000000 0.00773 0.00001818 -18000000 0.006426 0.00001994 -21000000 0.005541 0.00002181 -24000000 0.004919 0.0000237 -27000000 0.004471 0.00002559 -29500000 0.004194 0.00002714 -31000000 0.004031 0.00002806 -33000000 0.00391 0.00002832 -53000000 0.003868 0.00002935 +101325 1.0 0.00001200 +3000000 0.04234 0.00001344 +6000000 0.02046 0.0000142 +9000000 0.01328 0.00001526 +12000000 0.00977 0.0000166 +15000000 0.00773 0.00001818 +18000000 0.006426 0.00001994 +21000000 0.005541 0.00002181 +24000000 0.004919 0.0000237 +27000000 0.004471 0.00002559 +29500000 0.004194 0.00002714 +31000000 0.004031 0.00002806 +33000000 0.00391 0.00002832 +53000000 0.003868 0.00002935 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtw.txt b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtw.txt index 610f671eced..ad340949293 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtw.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtw.txt @@ -1,2 +1,2 @@ # Pref[bar] Bw[m3/sm3] Cp[1/bar] Visc[cP] - 30600000.1 1.03 0.00000000041 0.0003 \ No newline at end of file +30600000.1 1.03 0.00000000041 0.0003 diff --git a/src/coreComponents/physicsSolvers/multiphysics/integratedTests/FieldCaseCo2InjTutorial.xml b/src/coreComponents/physicsSolvers/multiphysics/integratedTests/FieldCaseCo2InjTutorial.xml index c197eb4dfde..1760c396985 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/integratedTests/FieldCaseCo2InjTutorial.xml +++ b/src/coreComponents/physicsSolvers/multiphysics/integratedTests/FieldCaseCo2InjTutorial.xml @@ -12,12 +12,16 @@ initialDt="1e2" targetRegions="{ Reservoir, wellRegion }"> + solverType="fgmres" + preconditionerType="mgr" + krylovTol="1e-4" + krylovAdaptiveTol="1" + krylovWeakestTol="1e-2"/> @@ -95,7 +99,7 @@ + maxTime="5e7"> @@ -211,7 +215,7 @@ setNames="{ all }" objectPath="ElementRegions/Reservoir" fieldName="pressure" - scale="4e7"/> + scale="1.25e7"/> + scale="0.0"/> + scale="1.0"/>