diff --git a/examples/AdaptiveTimeStepping/deadoil_3ph_baker_1d.xml b/examples/AdaptiveTimeStepping/deadoil_3ph_baker_1d.xml index 4c7a2197426..fa18ecde681 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..8ac8a2200df 100644 --- a/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp +++ b/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp @@ -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 new file mode 100644 index 00000000000..276c8f5afba --- /dev/null +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.cpp @@ -0,0 +1,432 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "DeadOilFluid.hpp" + +#include "constitutive/fluid/MultiFluidUtils.hpp" +#include "managers/GeosxState.hpp" +#include "managers/Functions/FunctionManager.hpp" + + +namespace geosx +{ + +using namespace dataRepository; +using namespace stringutilities; + +namespace constitutive +{ + +constexpr integer DeadOilFluid::PhaseType::GAS; +constexpr integer DeadOilFluid::PhaseType::OIL; +constexpr integer DeadOilFluid::PhaseType::WATER; + +namespace +{ + +std::unordered_map< string, integer > 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" ); + + // 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 ). + 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 ). + 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 = stringutilities::tokenize( str, " " ); + + if( strs.empty() ) + { + continue; + } + if( strs[0].front() == '#' ) + { + continue; + } + if( 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() +{ + 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_THROW_IF( it == phaseDict.end(), + "DeadOilFluid: phase not supported: " << m_phaseNames[ip], + InputError ); + integer const phaseIndex = it->second; + 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_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 ); + + // at this point, we make the distinction between the two input options + if( !m_tableFiles.empty() ) + { + readInputDataFromPVTFiles(); + } + else + { + useProvidedTableFunctions(); + } +} + +void DeadOilFluid::initializePostSubGroups() +{ + MultiFluidBase::initializePostSubGroups(); + createAllKernelWrappers(); +} + +void DeadOilFluid::createAllKernelWrappers() +{ + FunctionManager const & functionManager = getGlobalState().getFunctionManager(); + + 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 ) + { + + // loop over the hydrocarbon phases + for( localIndex iph = 0; iph < m_hydrocarbonPhaseOrder.size(); ++iph ) + { + // grab the tables by name from the function manager + TableFunction const & FVFTable = functionManager.getGroup< TableFunction const >( m_formationVolFactorTableNames[iph] ); + TableFunction const & viscosityTable = functionManager.getGroup< TableFunction const >( m_viscosityTableNames[iph] ); + 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_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_THROW_IF( (property[i] - property[i-1]) * (property[i-1] - property[i-2]) < 0, + "DeadOilFluid: the values in each PVT table must monotone", + InputError ); + } +} + +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..6c57fa6d3db --- /dev/null +++ b/src/coreComponents/constitutive/fluid/DeadOilFluid.hpp @@ -0,0 +1,753 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 + * This kernel can be called on the GPU (after adding GEOSX_HOST_DEVICE flags in this class and in MultiFluidBase) + */ +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 > 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 > 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 = default; + + 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 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"; } + static constexpr char const * waterFormationVolumeFactorString() { return "waterFormationVolumeFactor"; } + static constexpr char const * waterCompressibilityString() { return "waterCompressibility"; } + static constexpr char const * waterViscosityString() { return "waterViscosity"; } + + }; + + +protected: + + virtual void postProcessInput() override; + + virtual void initializePostSubGroups() override; + +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(); + + /// check that the table makes sense + void validateTable( TableFunction const & table ) const; + + // Input data + + // Black-oil table filenames + path_array m_tableFiles; + + // Fluid 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 > 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( real64 & 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( real64 & 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 > 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 ); + localIndex const nComps = m_componentMolarWeight.size(); + localIndex const nPhases = 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 < nPhases; ++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 < nPhases; ++ip ) + { + phaseFrac[ip] = composition[ip]; + for( localIndex ic = 0; ic < nComps; ++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 < nPhases; ++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 > 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 nComps = numComponents(); + localIndex const nPhases = 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 < nPhases; ++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 < nComps; ++ic ) + { + dPhaseDensity_dGlobalCompFraction[ip][ic] = 0.0; + } + } + + // 3. Update remaining variables: phaseFrac, phaseCompFrac using Dead-Oil assumptions + for( localIndex ip = 0; ip < nPhases; ++ip ) + { + phaseFraction[ip] = composition[ip]; + dPhaseFraction_dPressure[ip] = 0.0; + for( localIndex ic = 0; ic < nComps; ++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 < nComps; ++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 < nComps; ++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 < nPhases; ++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 < nComps; ++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 < nComps; ++ic ) + { + dTotalDensity_dGlobalCompFraction[ic] *= minusDens2; + } +} + +} //namespace constitutive + +} //namespace geosx + +#endif //GEOSX_CONSTITUTIVE_FLUID_DEADOILFLUID_HPP_ 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/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 2df414a3dbd..547b1301b83 100644 --- a/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp +++ b/src/coreComponents/constitutive/unitTests/testMultiFluid.cpp @@ -13,12 +13,15 @@ */ // Source includes -#include "managers/initialization.hpp" #include "common/DataTypes.hpp" #include "common/TimingMacros.hpp" #include "constitutive/fluid/multiFluidSelector.hpp" #include "constitutive/fluid/MultiFluidUtils.hpp" #include "physicsSolvers/fluidFlow/unitTests/testCompFlowUtils.hpp" +#include "managers/initialization.hpp" +#include "managers/Functions/FunctionManager.hpp" +#include "managers/GeosxState.hpp" + // TPL includes #include @@ -80,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, @@ -209,6 +210,7 @@ void testNumericalDerivatives( MultiFluidBase & fluid, "Pres", phases, components ); + std::cout << phaseDens.value << std::endl; } // update temperature and check derivatives @@ -392,22 +394,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 ); @@ -415,18 +412,122 @@ 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"; + tableNames[0] = "pvdo.txt"; tableNames[1] = "pvdw.txt"; tableNames[2] = "pvdg.txt"; - BlackOilFluid::FluidType & fluidType = fluid.getReference< BlackOilFluid::FluidType >( BlackOilFluid::viewKeyStruct::fluidTypeString() ); - fluidType = BlackOilFluid::FluidType::DeadOil; + fluid.postProcessInputRecursive(); + 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; @@ -532,8 +633,7 @@ TEST_F( DeadOilFluidTest, numericalDerivativesMolar ) { fluid->setMassFlag( false ); - // TODO test over a range of values - real64 const P = 5e6; + 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; @@ -541,31 +641,68 @@ TEST_F( DeadOilFluidTest, numericalDerivativesMolar ) real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); real64 const relTol = 1e-4; - testNumericalDerivatives( *fluid, parent, P, T, comp, eps, relTol ); + for( localIndex i = 0; i < 3; ++i ) + { + testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, relTol ); + } } TEST_F( DeadOilFluidTest, numericalDerivativesMass ) { fluid->setMassFlag( true ); - // TODO test over a range of values - real64 const P = 5e6; + 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-2; + real64 const relTol = 1e-4; real64 const absTol = 1e-14; - testNumericalDerivatives( *fluid, parent, P, T, comp, eps, relTol, absTol ); + for( localIndex i = 0; i < 3; ++i ) + { + testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, relTol, absTol ); + } +} + +class DeadOilFluidFromTableTest : public CompositionalFluidTestBase +{ +public: + + DeadOilFluidFromTableTest() + { + parent.resize( 1 ); + fluid = &makeDeadOilFluidFromTable( "fluid", &parent ); + + parent.initialize(); + parent.initializePostInitialConditions(); + } +}; + +TEST_F( DeadOilFluidFromTableTest, numericalDerivativesMolar ) +{ + fluid->setMassFlag( false ); + + real64 const P[3] = { 5.4e6, 1.24e7, 3.21e7 }; + real64 const T = 297.15; + array1d< real64 > comp( 3 ); + comp[0] = 0.1; comp[1] = 0.3; comp[2] = 0.6; + + real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); + real64 const relTol = 1e-4; + + for( localIndex i = 0; i < 3; ++i ) + { + testNumericalDerivatives( *fluid, parent, P[i], T, comp, eps, relTol ); + } } int main( int argc, char * * argv ) { ::testing::InitGoogleTest( &argc, argv ); - geosx::basicSetup( argc, argv ); + geosx::GeosxState state( geosx::basicSetup( argc, argv ) ); int const result = RUN_ALL_TESTS(); diff --git a/src/coreComponents/constitutive/unitTests/testRelPerm.cpp b/src/coreComponents/constitutive/unitTests/testRelPerm.cpp index 86a1fcfa0a4..cfb1a2aa3dd 100644 --- a/src/coreComponents/constitutive/unitTests/testRelPerm.cpp +++ b/src/coreComponents/constitutive/unitTests/testRelPerm.cpp @@ -167,7 +167,7 @@ RelativePermeabilityBase & makeVanGenuchtenBakerRelPermThreePhase( string const RelativePermeabilityBase & makeTableRelPermTwoPhase( string const & name, Group & parent ) { - FunctionManager * functionManager = &getGlobalState().getFunctionManager(); + FunctionManager & functionManager = getGlobalState().getFunctionManager(); // 1) First, define the tables @@ -191,14 +191,14 @@ RelativePermeabilityBase & makeTableRelPermTwoPhase( string const & name, Group values[i] = coordinates[0][i]*coordinates[0][i]; } - TableFunction & table_w = dynamicCast< TableFunction & >( *functionManager->createChild( "TableFunction", "water_swof" ) ); + TableFunction & table_w = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "water_swof" ) ); table_w.setTableCoordinates( coordinates ); table_w.setTableValues( values ); table_w.reInitializeFunction(); table_w.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - TableFunction & table_o = dynamicCast< TableFunction & >( *functionManager->createChild( "TableFunction", "oil_swof" ) ); + TableFunction & table_o = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "oil_swof" ) ); table_o.setTableCoordinates( coordinates ); table_o.setTableValues( values ); table_o.reInitializeFunction(); @@ -223,7 +223,7 @@ RelativePermeabilityBase & makeTableRelPermTwoPhase( string const & name, Group RelativePermeabilityBase & makeTableRelPermThreePhase( string const & name, Group & parent ) { - FunctionManager * functionManager = &getGlobalState().getFunctionManager(); + FunctionManager & functionManager = getGlobalState().getFunctionManager(); // 1) First, define the tables @@ -249,14 +249,14 @@ RelativePermeabilityBase & makeTableRelPermThreePhase( string const & name, Grou values[i] = coordinates[0][i]*coordinates[0][i]; } - TableFunction & table_ow_w = dynamicCast< TableFunction & >( *functionManager->createChild( "TableFunction", "water_swof" ) ); + TableFunction & table_ow_w = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "water_swof" ) ); table_ow_w.setTableCoordinates( coordinates ); table_ow_w.setTableValues( values ); table_ow_w.reInitializeFunction(); table_ow_w.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - TableFunction & table_ow_o = dynamicCast< TableFunction & >( *functionManager->createChild( "TableFunction", "oil_swof" ) ); + TableFunction & table_ow_o = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "oil_swof" ) ); table_ow_o.setTableCoordinates( coordinates ); table_ow_o.setTableValues( values ); table_ow_o.reInitializeFunction(); @@ -278,14 +278,14 @@ RelativePermeabilityBase & makeTableRelPermThreePhase( string const & name, Grou values[i] = coordinates[0][i]*coordinates[0][i]*coordinates[0][i]; } - TableFunction & table_og_g = dynamicCast< TableFunction & >( *functionManager->createChild( "TableFunction", "gas_sgof" ) ); + TableFunction & table_og_g = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "gas_sgof" ) ); table_og_g.setTableCoordinates( coordinates ); table_og_g.setTableValues( values ); table_og_g.reInitializeFunction(); table_og_g.setInterpolationMethod( TableFunction::InterpolationType::Linear ); - TableFunction & table_og_o = dynamicCast< TableFunction & >( *functionManager->createChild( "TableFunction", "oil_sgof" ) ); + TableFunction & table_og_o = dynamicCast< TableFunction & >( *functionManager.createChild( "TableFunction", "oil_sgof" ) ); table_og_o.setTableCoordinates( coordinates ); table_og_o.setTableValues( values ); table_og_o.reInitializeFunction(); 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/Constitutive.rst b/src/coreComponents/fileIO/schema/docs/Constitutive.rst index 911789f350c..d2c6312e746 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` DamageElasticIsotropic node :ref:`XML_DamageElasticIsotropic` DamageSpectralElasticIsotropic node :ref:`XML_DamageSpectralElasticIsotropic` DamageVolDevElasticIsotropic node :ref:`XML_DamageVolDevElasticIsotropic` +DeadOilFluid node :ref:`XML_DeadOilFluid` DruckerPrager node :ref:`XML_DruckerPrager` ElasticIsotropic node :ref:`XML_ElasticIsotropic` ElasticTransverseIsotropic node :ref:`XML_ElasticTransverseIsotropic` diff --git a/src/coreComponents/fileIO/schema/docs/Constitutive_other.rst b/src/coreComponents/fileIO/schema/docs/Constitutive_other.rst index c92dca56db8..9bc5975aa00 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` DamageElasticIsotropic node :ref:`DATASTRUCTURE_DamageElasticIsotropic` DamageSpectralElasticIsotropic node :ref:`DATASTRUCTURE_DamageSpectralElasticIsotropic` DamageVolDevElasticIsotropic node :ref:`DATASTRUCTURE_DamageVolDevElasticIsotropic` +DeadOilFluid node :ref:`DATASTRUCTURE_DeadOilFluid` DruckerPrager node :ref:`DATASTRUCTURE_DruckerPrager` ElasticIsotropic node :ref:`DATASTRUCTURE_ElasticIsotropic` ElasticTransverseIsotropic node :ref:`DATASTRUCTURE_ElasticTransverseIsotropic` diff --git a/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst b/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst new file mode 100644 index 00000000000..c84a4c0613d --- /dev/null +++ b/src/coreComponents/fileIO/schema/docs/DeadOilFluid.rst @@ -0,0 +1,24 @@ + + +======================================= ============ ======== ===================================================================================================================================================================================================================================================================================================== +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 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 052ca6f22b4..8d698b94889 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd +++ b/src/coreComponents/fileIO/schema/schema.xsd @@ -1632,6 +1632,7 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> + @@ -1655,10 +1656,6 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> - - @@ -1668,11 +1665,6 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> - - - - - @@ -1862,6 +1854,36 @@ 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 ed16d722796..a8e64a577ee 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd.other +++ b/src/coreComponents/fileIO/schema/schema.xsd.other @@ -695,6 +695,7 @@ + @@ -921,6 +922,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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/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/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; 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 52ea7578e83..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,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"/> @@ -588,29 +590,30 @@ @@ -687,20 +690,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 7196c77db8b..123259f7a49 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/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/blackoil_2ph_1d.xml b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/blackoil_2ph_1d.xml index 795e38bb695..f0a1b4bbb35 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/blackoil_2ph_1d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/blackoil_2ph_1d.xml @@ -104,7 +104,6 @@ + materialList="{ fluid, rock, relperm }"/> - @@ -133,7 +132,7 @@ component="1" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -142,7 +141,7 @@ component="2" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -150,7 +149,7 @@ name="referencePorosity" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="referencePorosity" scale="0.05"/> @@ -159,16 +158,16 @@ name="initialPressure" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="pressure" - scale="5e6"/> + scale="7.5e6"/> @@ -177,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"/> @@ -186,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"/> @@ -194,7 +193,7 @@ @@ -203,7 +202,7 @@ @@ -211,7 +210,7 @@ @@ -219,24 +218,24 @@ - + @@ -244,7 +243,7 @@ @@ -252,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 acd64252ee0..20627e153d8 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 @@ -6,12 +6,12 @@ name="compflow" logLevel="1" discretization="fluidTPFA" - fluidNames="{ fluid1 }" + fluidNames="{ fluid }" solidNames="{ rock }" relPermNames="{ relperm }" temperature="300" useMass="0" - targetRegions="{ Region1 }"> + targetRegions="{ region }"> @@ -27,7 +27,7 @@ + materialList="{ fluid, rock, relperm }"/> - @@ -131,7 +130,7 @@ component="1" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -140,7 +139,7 @@ component="2" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/Region1/block1" + objectPath="ElementRegions/region/block1" fieldName="permeability" scale="1.0e-16"/> @@ -148,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"/> @@ -175,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"/> @@ -184,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"/> @@ -192,7 +191,7 @@ @@ -201,7 +200,7 @@ @@ -209,7 +208,7 @@ @@ -217,24 +216,24 @@ - + @@ -242,7 +241,7 @@ @@ -250,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 04f4e6c642b..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 @@ -8,7 +8,7 @@ logLevel="1" discretization="fluidTPFA" targetRegions="{ Channel }" - fluidNames="{ fluid1 }" + fluidNames="{ fluid }" solidNames="{ rock }" relPermNames="{ relperm }" temperature="300" @@ -30,7 +30,7 @@ + materialList="{ fluid, rock, relperm }"/> - + + scale="1.0e-16"/> + scale="1.0e-16"/> + scale="1.0e-16"/> - + - + - + @@ -289,8 +293,33 @@ component="2" scale="0.001"/> - + + + + + + + + + + @@ -81,13 +81,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 }"/> @@ -166,14 +165,6 @@ component="1" scale="0.4"/> - 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..aa126cf3c87 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" @@ -28,7 +28,7 @@ + materialList="{ fluid, rock, relperm }"/> - + 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 ab829cbbba1..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 +# Pg(Pa) Bg(m3/sm3) Visc(Pa.s) +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/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/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/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 26c230ba649..46d517fbfe0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -1237,12 +1237,12 @@ void CompositionalMultiphaseWell::computePerforationRates( WellElementSubRegion numFluidPhases(), m_resPres.toNestedViewConst(), m_deltaResPres.toNestedViewConst(), - m_resPhaseMob.toNestedViewConst(), - m_dResPhaseMob_dPres.toNestedViewConst(), - m_dResPhaseMob_dCompDens.toNestedViewConst(), m_dResPhaseVolFrac_dPres.toNestedViewConst(), m_dResPhaseVolFrac_dCompDens.toNestedViewConst(), m_dResCompFrac_dCompDens.toNestedViewConst(), + m_resPhaseDens.toNestedViewConst(), + m_dResPhaseDens_dPres.toNestedViewConst(), + m_dResPhaseDens_dComp.toNestedViewConst(), m_resPhaseVisc.toNestedViewConst(), m_dResPhaseVisc_dPres.toNestedViewConst(), m_dResPhaseVisc_dComp.toNestedViewConst(), @@ -1419,18 +1419,6 @@ void CompositionalMultiphaseWell::resetViews( DomainPartition & domain ) m_dResCompFrac_dCompDens = elemManager.constructArrayViewAccessor< real64, 3 >( 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/dead_oil_wells_2d.xml b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/dead_oil_wells_2d.xml index 4106d64884a..38f16756cf5 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 @@ - + materialList="{ fluid, rock, relperm }"/> + materialList="{ fluid, relperm }"/> + materialList="{ fluid, relperm }"/> + materialList="{ fluid, relperm }"/> - + 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"/> diff --git a/src/docs/sphinx/CompleteXMLSchema.rst b/src/docs/sphinx/CompleteXMLSchema.rst index e061b304d34..9906bd9cecc 100644 --- a/src/docs/sphinx/CompleteXMLSchema.rst +++ b/src/docs/sphinx/CompleteXMLSchema.rst @@ -176,6 +176,13 @@ Element: DamageVolDevElasticIsotropic .. include:: ../../coreComponents/fileIO/schema/docs/DamageVolDevElasticIsotropic.rst +.. _XML_DeadOilFluid: + +Element: DeadOilFluid +===================== +.. include:: ../../coreComponents/fileIO/schema/docs/DeadOilFluid.rst + + .. _XML_Dirichlet: Element: Dirichlet @@ -929,6 +936,13 @@ Datastructure: DamageVolDevElasticIsotropic .. include:: ../../coreComponents/fileIO/schema/docs/DamageVolDevElasticIsotropic_other.rst +.. _DATASTRUCTURE_DeadOilFluid: + +Datastructure: DeadOilFluid +=========================== +.. include:: ../../coreComponents/fileIO/schema/docs/DeadOilFluid_other.rst + + .. _DATASTRUCTURE_Dirichlet: Datastructure: Dirichlet