diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 4e065cceae2..8678f64f335 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3359-14024-3c9ebb4 + baseline: integratedTests/baseline_integratedTests-pr3193-14118-8ee1c34 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 2e3a8b70338..34b169b92d9 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3193 (2025-10-10) +===================== +Enable geothermal gradient in HydrostaticEquilibrium for single-phase flow. + PR #3359 (2025-10-07) ===================== Add functions to connect well perforation to surface elements. diff --git a/inputFiles/poromechanics/PoroElastic_hybridHexPrism_co2_base.xml b/inputFiles/poromechanics/PoroElastic_hybridHexPrism_co2_base.xml index 85b5592eb52..86462b23ca7 100644 --- a/inputFiles/poromechanics/PoroElastic_hybridHexPrism_co2_base.xml +++ b/inputFiles/poromechanics/PoroElastic_hybridHexPrism_co2_base.xml @@ -44,7 +44,7 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Isothermal_smoke.xml b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Isothermal_smoke.xml new file mode 100644 index 00000000000..cfba8984843 --- /dev/null +++ b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Isothermal_smoke.xml @@ -0,0 +1,47 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Thermal_base.xml b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Thermal_base.xml new file mode 100755 index 00000000000..bbd526c564a --- /dev/null +++ b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Thermal_base.xml @@ -0,0 +1,243 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Thermal_smoke.xml b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Thermal_smoke.xml new file mode 100755 index 00000000000..be25c992d05 --- /dev/null +++ b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_Thermal_smoke.xml @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_base.xml b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_base.xml index 6de436890ed..8b38ef13757 100644 --- a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_base.xml +++ b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_base.xml @@ -61,15 +61,19 @@ + targetRegions="{ Reservoir }" + logLevel="1"> + newtonMaxIter="8" + lineSearchAction="None" + logLevel="1"/> + krylovTol="1.0e-10" + logLevel="1"/> @@ -156,6 +160,8 @@ + diff --git a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_smoke.xml b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_smoke.xml index df0089a6540..f8a726797e3 100644 --- a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_smoke.xml +++ b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_smoke.xml @@ -21,6 +21,11 @@ + + + diff --git a/inputFiles/singlePhaseFlow/singlePhaseFlow.ats b/inputFiles/singlePhaseFlow/singlePhaseFlow.ats index e6795bbe42f..1392414bd66 100644 --- a/inputFiles/singlePhaseFlow/singlePhaseFlow.ats +++ b/inputFiles/singlePhaseFlow/singlePhaseFlow.ats @@ -89,7 +89,15 @@ decks = [ partitions=((1, 1, 1), (2, 2, 1), (3, 3, 1)), restart_step=10, check_step=20, - restartcheck_params=RestartcheckParameters(**restartcheck_params)) + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="FieldCaseTutorial3_Isothermal_smoke", + description= + '3D Single phase flow test with burdens', + partitions=((1, 1, 1), (2, 2, 1)), + restart_step=0, + check_step=1, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), ] generate_geos_tests(decks) diff --git a/inputFiles/singlePhaseFlow/thermalSinglePhaseFlow.ats b/inputFiles/singlePhaseFlow/thermalSinglePhaseFlow.ats index 886d8add60a..b47f93c6527 100644 --- a/inputFiles/singlePhaseFlow/thermalSinglePhaseFlow.ats +++ b/inputFiles/singlePhaseFlow/thermalSinglePhaseFlow.ats @@ -36,7 +36,15 @@ decks = [ partitions=((1, 1, 1), (2, 2, 1)), restart_step=1, check_step=2, - restartcheck_params=RestartcheckParameters(**restartcheck_params)) + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="FieldCaseTutorial3_Thermal_smoke", + description= + '3D thermal single phase flow with geothermal gradient', + partitions=((1, 1, 1), (2, 2, 1)), + restart_step=0, + check_step=1, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), ] generate_geos_tests(decks) diff --git a/inputFiles/thermoPoromechanics/ThermoPoroElastic_SinglePhase_ThermalGradient_base.xml b/inputFiles/thermoPoromechanics/ThermoPoroElastic_SinglePhase_ThermalGradient_base.xml index 3f2e79211b2..b3891e72d94 100755 --- a/inputFiles/thermoPoromechanics/ThermoPoroElastic_SinglePhase_ThermalGradient_base.xml +++ b/inputFiles/thermoPoromechanics/ThermoPoroElastic_SinglePhase_ThermalGradient_base.xml @@ -160,7 +160,8 @@ name="INITIALISATION" datumElevation="-3250" datumPressure="3.763746e7" - objectPath="ElementRegions"/> + objectPath="ElementRegions" + temperatureVsElevationTableName="initTempTable"/> ::concat( "\n* " ) ); @@ -108,14 +108,15 @@ void CompressibleSinglePhaseFluid::postInputInitialization() checkPositive( m_referenceViscosity, viewKeyStruct::referenceViscosityString() ); // Due to the way update wrapper is currently implemented, we can only support one model type - auto const checkModelType = [&]( ExponentApproximationType const value, auto const & attribute ) + auto const checkModelType = [&]( ExponentApproximationType const value, ExponentApproximationType const expectedValue, auto const & attribute ) { - GEOS_THROW_IF_NE_MSG( value, ExponentApproximationType::Linear, - GEOS_FMT( "{}: invalid model type in attribute '{}' (only linear currently supported)", getFullName(), attribute ), - InputError ); + GEOS_THROW_IF( value != expectedValue, + GEOS_FMT( "{}: invalid model type in attribute '{}' (only {} currently supported)", + getFullName(), attribute, EnumStrings< ExponentApproximationType >::toString( expectedValue ) ), + InputError ); }; - checkModelType( m_densityModelType, viewKeyStruct::densityModelTypeString() ); - checkModelType( m_viscosityModelType, viewKeyStruct::viscosityModelTypeString() ); + checkModelType( m_densityModelType, ExponentApproximationType::Full, viewKeyStruct::densityModelTypeString() ); + checkModelType( m_viscosityModelType, ExponentApproximationType::Linear, viewKeyStruct::viscosityModelTypeString() ); // Set default values for derivatives (cannot be done in base class) // TODO: reconsider the necessity of this diff --git a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp index 4446616bb9c..09faa09ab49 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp @@ -153,8 +153,8 @@ class CompressibleSinglePhaseFluid : public SingleFluidBase virtual void allocateConstitutiveData( dataRepository::Group & parent, localIndex const numPts ) override; - /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT, not just linear) - using KernelWrapper = CompressibleSinglePhaseUpdate< ExponentApproximationType::Linear, ExponentApproximationType::Linear >; + /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT combinations, not just this combination) + using KernelWrapper = CompressibleSinglePhaseUpdate< ExponentApproximationType::Full, ExponentApproximationType::Linear >; /** * @brief Create an update kernel wrapper. diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp index 351d9743ace..57baebcd702 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp @@ -131,6 +131,41 @@ class SingleFluidBaseUpdate dViscosity_dPressure ); } + template< typename FLUIDWRAPPER > + GEOS_HOST_DEVICE + static void computeThermalValues( FLUIDWRAPPER const fluidWrapper, + real64 const pressure, + real64 const temperature, + real64 & density, + real64 & viscosity ) + { + real64 dDensity_dPressure = 0.0; + real64 dViscosity_dPressure = 0.0; + real64 dDensity_dTemperature = 0.0; + real64 dViscosity_dTemperature = 0.0; + real64 internalEnergy = 0.0; + real64 dInternalEnergy_dPressure = 0.0; + real64 dInternalEnergy_dTemperature = 0.0; + real64 enthalpy = 0.0; + real64 dEnthalpy_dPressure = 0.0; + real64 dEnthalpy_dTemperature = 0.0; + + fluidWrapper.compute( pressure, + temperature, + density, + dDensity_dPressure, + dDensity_dTemperature, + viscosity, + dViscosity_dPressure, + dViscosity_dTemperature, + internalEnergy, + dInternalEnergy_dPressure, + dInternalEnergy_dTemperature, + enthalpy, + dEnthalpy_dPressure, + dEnthalpy_dTemperature ); + } + //START_SPHINX_INCLUDE_02 private: /** diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp index 0e32f717116..2e51f625ed5 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp @@ -32,7 +32,6 @@ namespace constitutive ThermalCompressibleSinglePhaseFluid::ThermalCompressibleSinglePhaseFluid( string const & name, Group * const parent ): CompressibleSinglePhaseFluid( name, parent ) { - m_densityModelType = ExponentApproximationType::Full; m_numDOF=2; registerWrapper( viewKeyStruct::thermalExpansionCoeffString(), &m_thermalExpansionCoeff ). setApplyDefaultValue( 0.0 ). @@ -85,13 +84,14 @@ void ThermalCompressibleSinglePhaseFluid::postInputInitialization() checkNonnegative( m_referenceInternalEnergy, viewKeyStruct::referenceInternalEnergyString() ); // Due to the way update wrapper is currently implemented, we can only support one model type - auto const checkModelType = [&]( ExponentApproximationType const value, auto const & attribute ) + auto const checkModelType = [&]( ExponentApproximationType const value, ExponentApproximationType const expectedValue, auto const & attribute ) { - GEOS_THROW_IF( value != ExponentApproximationType::Linear && value != ExponentApproximationType::Full, - GEOS_FMT( "{}: invalid model type in attribute '{}' (only linear or fully exponential currently supported)", getFullName(), attribute ), + GEOS_THROW_IF( value != expectedValue, + GEOS_FMT( "{}: invalid model type in attribute '{}' (only {} currently supported)", + getFullName(), attribute, EnumStrings< ExponentApproximationType >::toString( expectedValue ) ), InputError ); }; - checkModelType( m_internalEnergyModelType, viewKeyStruct::internalEnergyModelTypeString() ); + checkModelType( m_internalEnergyModelType, ExponentApproximationType::Linear, viewKeyStruct::internalEnergyModelTypeString() ); } ThermalCompressibleSinglePhaseFluid::KernelWrapper diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp index d16dc3b37a7..574164785b7 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp @@ -204,7 +204,7 @@ class ThermalCompressibleSinglePhaseFluid : public CompressibleSinglePhaseFluid using CompressibleSinglePhaseFluid::m_densityModelType; - /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT, not just linear) + /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT combinations, not just this combination) using KernelWrapper = ThermalCompressibleSinglePhaseUpdate< ExponentApproximationType::Full, ExponentApproximationType::Linear, ExponentApproximationType::Linear >; /** diff --git a/src/coreComponents/fieldSpecification/EquilibriumInitialCondition.cpp b/src/coreComponents/fieldSpecification/EquilibriumInitialCondition.cpp index 296dd53bb67..30a41770e9d 100644 --- a/src/coreComponents/fieldSpecification/EquilibriumInitialCondition.cpp +++ b/src/coreComponents/fieldSpecification/EquilibriumInitialCondition.cpp @@ -93,12 +93,6 @@ EquilibriumInitialCondition::EquilibriumInitialCondition( string const & name, G void EquilibriumInitialCondition::postInputInitialization() { - GEOS_THROW_IF( ( m_temperatureVsElevationTableName.empty() != m_componentFractionVsElevationTableNames.empty() ), - getCatalogName() << " " << getDataContext() << ": both " << - viewKeyStruct::componentFractionVsElevationTableNamesString() << " and " << - viewKeyStruct::temperatureVsElevationTableNameString() << " must be provided for a multiphase simulation", - InputError ); - FunctionManager const & functionManager = FunctionManager::getInstance(); if( !m_componentFractionVsElevationTableNames.empty() ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index bb4d66c6122..0df05edeea1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -44,6 +44,7 @@ #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/HydrostaticPressureKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalHydrostaticPressureKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/FluidUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolidInternalEnergyUpdateKernel.hpp" @@ -421,6 +422,26 @@ void SinglePhaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) " - Use a gravityVector aligned with the z-axis, such as (0.0,0.0,-9.81)\n" << " - Remove the hydrostatic equilibrium initial condition from the XML file", InputError ); + + // ensure that the temperature tables are defined for thermal simulations + GEOS_THROW_IF( m_isThermal && bc.getTemperatureVsElevationTableName().empty(), + getCatalogName() << " " << bc.getDataContext() + << ": " << EquilibriumInitialCondition::viewKeyStruct::temperatureVsElevationTableNameString() + << " must be provided for a thermal simulation", + InputError ); + + //ensure that compositions are empty + GEOS_THROW_IF( !bc.getComponentFractionVsElevationTableNames().empty(), + getCatalogName() << " " << bc.getDataContext() + << ": " << EquilibriumInitialCondition::viewKeyStruct::componentFractionVsElevationTableNamesString() + << " must not be provided for a single phase simulation.", + InputError ); + + GEOS_THROW_IF( !bc.getComponentNames().empty(), + getCatalogName() << " " << bc.getDataContext() + << ": " << EquilibriumInitialCondition::viewKeyStruct::componentNamesString() + << " must not be provided for a single phase simulation.", + InputError ); } ); if( equilCounter == 0 ) @@ -492,6 +513,16 @@ void SinglePhaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) // Step 3.2: retrieve the fluid model to compute densities // we end up with the same issue as in applyDirichletBC: there is not a clean way to retrieve the fluid info + FunctionManager & functionManager = FunctionManager::getInstance(); + TableFunction::KernelWrapper tempTableWrapper; + // Creation of Wrapper in case of TemperatureVsElevationTableName + if( m_isThermal ) + { + string const tempTableName = fs.getTemperatureVsElevationTableName(); + TableFunction const & tempTable = functionManager.getGroup< TableFunction >( tempTableName ); + tempTableWrapper = tempTable.createKernelWrapper(); + } + // filter out region not in target Group const & region = subRegion.getParent().getParent(); auto it = regionFilter.find( region.getName() ); @@ -518,18 +549,32 @@ void SinglePhaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) typename FluidType::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); // note: inside this kernel, serialPolicy is used, and elevation/pressure values don't go to the GPU - bool const equilHasConverged = - HydrostaticPressureKernel::launch( numPointsInTable, - maxNumEquilIterations, - equilTolerance, - gravVector, - minElevation, - elevationIncrement, - datumElevation, - datumPressure, - fluidWrapper, - elevationValues.toNestedView(), - pressureValues.toView() ); + bool const equilHasConverged = m_isThermal ? + thermalSinglePhaseBaseKernels:: + HydrostaticPressureKernel::launch( numPointsInTable, + maxNumEquilIterations, + equilTolerance, + gravVector, + minElevation, + elevationIncrement, + datumElevation, + datumPressure, + fluidWrapper, + tempTableWrapper, + elevationValues.toNestedView(), + pressureValues.toView() ) : + singlePhaseBaseKernels:: + HydrostaticPressureKernel::launch( numPointsInTable, + maxNumEquilIterations, + equilTolerance, + gravVector, + minElevation, + elevationIncrement, + datumElevation, + datumPressure, + fluidWrapper, + elevationValues.toNestedView(), + pressureValues.toView() ); GEOS_THROW_IF( !equilHasConverged, getCatalogName() << " " << getDataContext() << @@ -539,8 +584,6 @@ void SinglePhaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) // Step 3.4: create hydrostatic pressure table - FunctionManager & functionManager = FunctionManager::getInstance(); - string const tableName = fs.getName() + "_" + subRegion.getName() + "_table"; TableFunction * const presTable = dynamicCast< TableFunction * >( functionManager.createChild( TableFunction::catalogName(), tableName ) ); presTable->setTableCoordinates( elevationValues, { units::Distance } ); @@ -552,14 +595,20 @@ void SinglePhaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) // TODO: this last step should probably be delayed to wait for the creation of FaceElements arrayView2d< real64 const > const elemCenter = subRegion.getElementCenter(); arrayView1d< real64 > const pres = subRegion.getField< flow::pressure >(); + arrayView1d< real64 > const temp = subRegion.getField< flow::temperature >(); - RAJA::ReduceMin< parallelDeviceReduce, real64 > minPressure( LvArray::NumericLimits< real64 >::max ); - forAll< parallelDevicePolicy< > >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + RAJA::ReduceMin< parallelHostReduce, real64 > minPressure( LvArray::NumericLimits< real64 >::max ); + + forAll< parallelHostPolicy >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { localIndex const k = targetSet[i]; real64 const elevation = elemCenter[k][2]; pres[k] = presTableWrapper.compute( &elevation ); + if( m_isThermal ) + { + temp[k] = tempTableWrapper.compute( &elevation ); + } minPressure.min( pres[k] ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalHydrostaticPressureKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalHydrostaticPressureKernel.hpp new file mode 100644 index 00000000000..d6615b3b424 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalHydrostaticPressureKernel.hpp @@ -0,0 +1,223 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ThermalHydrostaticPressureKernel.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALHYDROSTATICPRESSUREKERNEL_HPP +#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALHYDROSTATICPRESSUREKERNEL_HPP + +#include "common/DataTypes.hpp" +#include "common/GEOS_RAJA_Interface.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" + +namespace geos +{ + +namespace thermalSinglePhaseBaseKernels +{ + +/******************************** HydrostaticPressureKernel ********************************/ + +struct HydrostaticPressureKernel +{ + + template< typename FLUID_WRAPPER > + static bool + computeHydrostaticPressure( integer const maxNumEquilIterations, + real64 const & equilTolerance, + real64 const (&gravVector)[ 3 ], + FLUID_WRAPPER fluidWrapper, + TableFunction::KernelWrapper tempTableWrapper, + real64 const & refElevation, + real64 const & refPres, + real64 const & refDens, + real64 const & newElevation, + real64 & newPres, + real64 & newDens ) + { + // Step 1: guess the pressure with the refDensity + + real64 const gravCoef = gravVector[2] * ( refElevation - newElevation ); + real64 const temp = tempTableWrapper.compute( &newElevation ); + real64 pres0 = refPres - refDens * gravCoef; + real64 pres1 = 0.0; + + // Step 2: compute the mass density at this elevation using the guess, and update pressure + + real64 dens = 0.0; + real64 visc = 0.0; + constitutive::SingleFluidBaseUpdate::computeThermalValues( fluidWrapper, + pres0, + temp, + dens, + visc ); + pres1 = refPres - 0.5 * ( refDens + dens ) * gravCoef; + + // Step 3: fixed-point iteration until convergence + + bool equilHasConverged = false; + for( localIndex eqIter = 0; eqIter < maxNumEquilIterations; ++eqIter ) + { + + // check convergence + equilHasConverged = ( LvArray::math::abs( pres0 - pres1 ) < equilTolerance ); + pres0 = pres1; + + // if converged, move on + if( equilHasConverged ) + { + break; + } + + // compute the density at this elevation using the previous pressure, and compute the new pressure + constitutive::SingleFluidBaseUpdate::computeThermalValues( fluidWrapper, + pres0, + temp, + dens, + visc ); + pres1 = refPres - 0.5 * ( refDens + dens ) * gravCoef; + } + + // Step 4: save the hydrostatic pressure and the corresponding density + + newPres = pres1; + newDens = dens; + + return equilHasConverged; + } + + + template< typename FLUID_WRAPPER > + static bool + launch( localIndex const size, + integer const maxNumEquilIterations, + real64 const equilTolerance, + real64 const (&gravVector)[ 3 ], + real64 const & minElevation, + real64 const & elevationIncrement, + real64 const & datumElevation, + real64 const & datumPres, + FLUID_WRAPPER fluidWrapper, + TableFunction::KernelWrapper tempTableWrapper, + arrayView1d< arrayView1d< real64 > const > elevationValues, + arrayView1d< real64 > pressureValues ) + { + bool hasConverged = true; + + // Step 1: compute the mass density at the datum elevation + + real64 datumDens = 0.0; + real64 datumVisc = 0.0; + real64 const datumTemp = tempTableWrapper.compute( &datumElevation ); + + constitutive::SingleFluidBaseUpdate::computeThermalValues( fluidWrapper, + datumPres, + datumTemp, + datumDens, + datumVisc ); + + // Step 2: find the closest elevation to datumElevation + + forAll< parallelHostPolicy >( size, [=] ( localIndex const i ) + { + real64 const elevation = minElevation + i * elevationIncrement; + elevationValues[0][i] = elevation; + } ); + integer const iRef = LvArray::sortedArrayManipulation::find( elevationValues[0].begin(), + elevationValues[0].size(), + datumElevation ); + + + // Step 3: compute the mass density and pressure at the reference elevation + + array1d< real64 > dens( pressureValues.size() ); + + bool const hasConvergedRef = + computeHydrostaticPressure( maxNumEquilIterations, + equilTolerance, + gravVector, + fluidWrapper, + tempTableWrapper, + datumElevation, + datumPres, + datumDens, + elevationValues[0][iRef], + pressureValues[iRef], + dens[iRef] ); + if( !hasConvergedRef ) + { + return false; + } + + + // Step 4: for each elevation above the reference elevation, compute the pressure + + localIndex const numEntriesAboveRef = size - iRef - 1; + forAll< serialPolicy >( numEntriesAboveRef, [=, &hasConverged] ( localIndex const i ) + { + bool const hasConvergedAboveRef = + computeHydrostaticPressure( maxNumEquilIterations, + equilTolerance, + gravVector, + fluidWrapper, + tempTableWrapper, + elevationValues[0][iRef+i], + pressureValues[iRef+i], + dens[iRef+i], + elevationValues[0][iRef+i+1], + pressureValues[iRef+i+1], + dens[iRef+i+1] ); + if( !hasConvergedAboveRef ) + { + hasConverged = false; + } + + + } ); + + // Step 5: for each elevation below the reference elevation, compute the pressure + + localIndex const numEntriesBelowRef = iRef; + forAll< serialPolicy >( numEntriesBelowRef, [=, &hasConverged] ( localIndex const i ) + { + bool const hasConvergedBelowRef = + computeHydrostaticPressure( maxNumEquilIterations, + equilTolerance, + gravVector, + fluidWrapper, + tempTableWrapper, + elevationValues[0][iRef-i], + pressureValues[iRef-i], + dens[iRef-i], + elevationValues[0][iRef-i-1], + pressureValues[iRef-i-1], + dens[iRef-i-1] ); + if( !hasConvergedBelowRef ) + { + hasConverged = false; + } + } ); + + return hasConverged; + } +}; + +} // namespace thermalSinglePhaseBaseKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALHYDROSTATICPRESSUREKERNEL_HPP