diff --git a/.gitignore b/.gitignore index 99d92dc0e0b..381406c3b92 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ /build*/ +/localTests*/ /install*/ .cproject .project @@ -20,6 +21,7 @@ uberenv_libs spack-*.txt /src/docs/sphinx/datastructure +/scripts/myScript/ *~ __pycache__ diff --git a/host-configs/DTU.cmake b/host-configs/DTU.cmake new file mode 100644 index 00000000000..f85517f2597 --- /dev/null +++ b/host-configs/DTU.cmake @@ -0,0 +1,76 @@ +site_name(HOST_NAME) + +if(NOT DEFINED CONFIG_NAME) + set(CONFIG_NAME "${HOST_NAME}-DTU" CACHE STRING "Config name") +endif() + +message("CONFIG_NAME = ${CONFIG_NAME}") + +set(Python3_EXECUTABLE "/usr/bin/python3" CACHE PATH "") +set(CMAKE_C_COMPILER "gcc" CACHE PATH "") +set(CMAKE_C_FLAGS_RELEASE "-O3 -DNDEBUG -Wno-error=array-bounds -Wno-error=dangling-reference -Wno-error=maybe-uninitialized -lgomp -lpthread -lm -ldl" CACHE STRING "") +set(CMAKE_CXX_COMPILER "g++" CACHE PATH "") +set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG -Wno-error=array-bounds -Wno-error=dangling-reference -Wno-error=maybe-uninitialized -lgomp -lpthread -lm -ldl" CACHE STRING "") +set(CMAKE_Fortran_COMPILER "gfortran" CACHE PATH "") +set(CMAKE_Fortran_FLAGS_RELEASE "-O3 -DNDEBUG -Wno-error=array-bounds -Wno-error=dangling-reference -Wno-error=maybe-uninitialized -lgomp -lpthread -lm -ldl" CACHE STRING "") +set(CMAKE_VERBOSE_MAKEFILE ON) + +set(ENABLE_FORTRAN ON CACHE BOOL "") + +set(ENABLE_MPI ON CACHE BOOL "") +set(MPI_C_COMPILER "mpicc" CACHE PATH "") +set(MPI_CXX_COMPILER "mpic++" CACHE PATH "") +set(MPI_Fortran_COMPILER "mpifort" CACHE PATH "") +set(MPIEXEC "mpirun" CACHE PATH "") + +set(ENABLE_GTEST_DEATH_TESTS OFF CACHE BOOL "") + +set(ENABLE_PAMELA ON CACHE BOOL "") +set(ENABLE_PVTPackage ON CACHE BOOL "") + +set(GEOS_TPL_DIR "/home/behzadh/thirdPartyLibs/install-DTU-release" CACHE PATH "") +set(ENABLE_XML OFF CACHE BOOL "") +set(ENABLE_XML_UPDATES OFF CACHE BOOL "") +set(GEOS_ENABLE_TESTS OFF CACHE BOOL "") +set(ENABLE_EXAMPLES OFF CACHE BOOL "") +set(ENABLE_BENCHMARKS OFF CACHE BOOL "") +set(ENABLE_DOCS OFF CACHE BOOL "") +set(ENABLE_HYPRE ON CACHE BOOL "" FORCE) +set(ENABLE_VTK ON CACHE BOOL "" FORCE) +set(ENABLE_PVTPackage ON CACHE BOOL "" FORCE) + +# GTEST +set(ENABLE_GTEST_DEATH_TESTS OFF CACHE BOOL "") +set(gtest_disable_pthreads ON CACHE BOOL "") + +# disable most binaries and doc generation +set(ENABLE_TESTS OFF CACHE BOOL "" ) +set(DISABLE_UNIT_TESTS ON CACHE BOOL "" ) +set(ENABLE_EXAMPLES OFF CACHE BOOL "" ) +set(ENABLE_BENCHMARKS OFF CACHE BOOL "" ) +set(ENABLE_DOCS OFF CACHE BOOL "" ) + +set(ENABLE_MATHPRESSO OFF CACHE BOOL "") + +set(GEOSX_BUILD_OBJ_LIBS ON CACHE BOOL "") + +set(CUDA_ENABLED "OFF" CACHE BOOL "") +set(ENABLE_OPENMP "ON" CACHE BOOL "") + +option(ENABLE_CALIPER "Enables CALIPER" ON) + +include(${CMAKE_CURRENT_LIST_DIR}/tpls.cmake) + +set(ENABLE_MKL ON CACHE BOOL "") +set(MKL_ROOT "/opt/intel/oneapi/mkl/latest" ) +set(MKL_INCLUDE_DIRS "/opt/intel/oneapi/mkl/latest/include" CACHE STRING "") +set(MKL_LIBRARIES "/opt/intel/oneapi/mkl/latest/lib/intel64/libmkl_rt.so" CACHE STRING "") + +if(NOT MKL_LIBRARIES OR NOT MKL_INCLUDE_DIRS) + message(FATAL_ERROR "MKL not found") +endif() + +set_property(GLOBAL PROPERTY TARGET_SUPPORTS_SHARED_LIBS TRUE) + +set(CMAKE_SKIP_RULE_DEPENDENCY ON) +set(CMAKE_SKIP_INSTALL_RULES ON) \ No newline at end of file diff --git a/host-configs/DTU_d.cmake b/host-configs/DTU_d.cmake new file mode 100644 index 00000000000..8ba8c1687c3 --- /dev/null +++ b/host-configs/DTU_d.cmake @@ -0,0 +1,76 @@ +site_name(HOST_NAME) + +if(NOT DEFINED CONFIG_NAME) + set(CONFIG_NAME "${HOST_NAME}-DTU" CACHE STRING "Config name") +endif() + +message("CONFIG_NAME = ${CONFIG_NAME}") + +set(Python3_EXECUTABLE "/usr/bin/python3" CACHE PATH "") +set(CMAKE_C_COMPILER "gcc" CACHE PATH "") +set(CMAKE_C_FLAGS_RELEASE "-O3 -DNDEBUG -Wno-error=array-bounds -Wno-error=dangling-reference -Wno-error=maybe-uninitialized -lgomp -lpthread -lm -ldl" CACHE STRING "") +set(CMAKE_CXX_COMPILER "g++" CACHE PATH "") +set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG -Wno-error=array-bounds -Wno-error=dangling-reference -Wno-error=maybe-uninitialized -lgomp -lpthread -lm -ldl" CACHE STRING "") +set(CMAKE_Fortran_COMPILER "gfortran" CACHE PATH "") +set(CMAKE_Fortran_FLAGS_RELEASE "-O3 -DNDEBUG -Wno-error=array-bounds -Wno-error=dangling-reference -Wno-error=maybe-uninitialized -lgomp -lpthread -lm -ldl" CACHE STRING "") +set(CMAKE_VERBOSE_MAKEFILE ON) + +set(ENABLE_FORTRAN ON CACHE BOOL "") + +set(ENABLE_MPI ON CACHE BOOL "") +set(MPI_C_COMPILER "mpicc" CACHE PATH "") +set(MPI_CXX_COMPILER "mpic++" CACHE PATH "") +set(MPI_Fortran_COMPILER "mpifort" CACHE PATH "") +set(MPIEXEC "mpirun" CACHE PATH "") + +set(ENABLE_GTEST_DEATH_TESTS OFF CACHE BOOL "") + +set(ENABLE_PAMELA ON CACHE BOOL "") +set(ENABLE_PVTPackage ON CACHE BOOL "") + +set(GEOS_TPL_DIR "/home/behzadh/thirdPartyLibs/install" CACHE PATH "") +set(ENABLE_XML OFF CACHE BOOL "") +set(ENABLE_XML_UPDATES OFF CACHE BOOL "") +set(GEOS_ENABLE_TESTS OFF CACHE BOOL "") +set(ENABLE_EXAMPLES OFF CACHE BOOL "") +set(ENABLE_BENCHMARKS OFF CACHE BOOL "") +set(ENABLE_DOCS OFF CACHE BOOL "") +set(ENABLE_HYPRE ON CACHE BOOL "" FORCE) +set(ENABLE_VTK ON CACHE BOOL "" FORCE) +set(ENABLE_PVTPackage ON CACHE BOOL "" FORCE) + +# GTEST +set(ENABLE_GTEST_DEATH_TESTS OFF CACHE BOOL "") +set(gtest_disable_pthreads ON CACHE BOOL "") + +# disable most binaries and doc generation +set(ENABLE_TESTS OFF CACHE BOOL "" ) +set(DISABLE_UNIT_TESTS ON CACHE BOOL "" ) +set(ENABLE_EXAMPLES OFF CACHE BOOL "" ) +set(ENABLE_BENCHMARKS OFF CACHE BOOL "" ) +set(ENABLE_DOCS OFF CACHE BOOL "" ) + +set(ENABLE_MATHPRESSO OFF CACHE BOOL "") + +set(GEOSX_BUILD_OBJ_LIBS ON CACHE BOOL "") + +set(CUDA_ENABLED "OFF" CACHE BOOL "") +set(ENABLE_OPENMP "ON" CACHE BOOL "") + +option(ENABLE_CALIPER "Enables CALIPER" ON) + +include(${CMAKE_CURRENT_LIST_DIR}/tpls.cmake) + +set(ENABLE_MKL ON CACHE BOOL "") +set(MKL_ROOT "/opt/intel/oneapi/mkl/latest/") +find_library(MKL_LIBRARIES NAMES mkl_rt PATHS ${MKL_ROOT}/lib/intel64 NO_DEFAULT_PATH) +find_path(MKL_INCLUDE_DIRS NAMES mkl.h PATHS ${MKL_ROOT}/include NO_DEFAULT_PATH) + +if(NOT MKL_LIBRARIES OR NOT MKL_INCLUDE_DIRS) + message(FATAL_ERROR "MKL not found") +endif() + +set_property(GLOBAL PROPERTY TARGET_SUPPORTS_SHARED_LIBS TRUE) + +set(CMAKE_SKIP_RULE_DEPENDENCY ON) +set(CMAKE_SKIP_INSTALL_RULES ON) \ No newline at end of file diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp index f251d632117..0a17d25c29f 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp @@ -34,6 +34,111 @@ namespace constitutive namespace compositional { +namespace negative_two_phase_flash_internal +{ + +/** + * @brief Run a short successive-substitution flash using K-values initialized at @p wilsonPressure, + * while evaluating all fugacities at @p fugacityPressure (the true cell pressure). + * @return Equilibrium residual norm after the last (or converged) iteration. + */ +template< int USD2 > +GEOS_HOST_DEVICE +real64 probeEquilibriumResidual( integer const numComps, + real64 const wilsonPressure, + real64 const fugacityPressure, + real64 const temperature, + arraySlice1d< real64 const > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + FlashData const & flashData, + arraySlice1d< integer const > const & presentComponents, + real64 const flashTolerance, + integer const maxProbeIterations, + arraySlice1d< real64 > const & logLiquidFugacity, + arraySlice1d< real64 > const & logVapourFugacity, + arraySlice1d< real64 > const & fugacityRatios, + arraySlice1d< real64, USD2 > const & liquidComposition, + arraySlice1d< real64, USD2 > const & vapourComposition ) +{ + constexpr integer maxNumComps = MultiFluidConstants::MAX_NUM_COMPONENTS; + StackArray< real64, 1, maxNumComps > kProbe( numComps ); + if( flashData.liquidEos == EquationOfStateType::SoreideWhitson ) + { + KValueInitialization::computeSoreideWhitsonKvalue( numComps, + wilsonPressure, + temperature, + componentProperties, + composition, + kProbe.toSlice() ); + } + else + { + KValueInitialization::computeWilsonGasLiquidKvalue( numComps, + wilsonPressure, + temperature, + componentProperties, + kProbe.toSlice() ); + } + for( integer ic = 0; ic < numComps; ++ic ) + { + liquidComposition[ic] = composition[ic]; + vapourComposition[ic] = composition[ic]; + } + + real64 error = LvArray::NumericLimits< real64 >::max; + for( localIndex iterationCount = 0; iterationCount < maxProbeIterations; ++iterationCount ) + { + real64 const vapourFrac = RachfordRice::solve( kProbe.toSliceConst(), composition, presentComponents ); + + for( integer ic = 0; ic < numComps; ++ic ) + { + liquidComposition[ic] = composition[ic] / ( 1.0 + vapourFrac * ( kProbe[ic] - 1.0 ) ); + vapourComposition[ic] = kProbe[ic] * liquidComposition[ic]; + } + + normalizeComposition( numComps, liquidComposition ); + normalizeComposition( numComps, vapourComposition ); + + FugacityCalculator::computeLogFugacity( numComps, + fugacityPressure, + temperature, + liquidComposition.toSliceConst(), + componentProperties, + flashData.liquidEos, + flashData, + logLiquidFugacity ); + FugacityCalculator::computeLogFugacity( numComps, + fugacityPressure, + temperature, + vapourComposition.toSliceConst(), + componentProperties, + flashData.vapourEos, + flashData, + logVapourFugacity ); + + error = 0.0; + for( integer const ic : presentComponents ) + { + fugacityRatios[ic] = ( logLiquidFugacity[ic] - logVapourFugacity[ic] ) + log( liquidComposition[ic] ) - log( vapourComposition[ic] ); + error += (fugacityRatios[ic]*fugacityRatios[ic]); + } + error = LvArray::math::sqrt( error ); + + if( error < flashTolerance ) + { + break; + } + + for( integer ic = 0; ic < numComps; ++ic ) + { + kProbe[ic] *= exp( fugacityRatios[ic] ); + } + } + return error; +} + +} // namespace negative_two_phase_flash_internal + template< int USD1, int USD2 > GEOS_HOST_DEVICE bool NegativeTwoPhaseFlash::compute( integer const numComps, @@ -73,59 +178,152 @@ bool NegativeTwoPhaseFlash::compute( integer const numComps, } bool converged = false; - for( localIndex iterationCount = 0; iterationCount < maxIterations; ++iterationCount ) - { - // Solve Rachford-Rice Equation - vapourPhaseMoleFraction = RachfordRice::solve( kVapourLiquid.toSliceConst(), composition, presentComponents ); + real64 lastInnerError = 0.0; + real64 errAfterWilsonAtPressure = 0.0; + bool haveErrAfterWilsonAtPressure = false; - // Assign phase compositions - for( integer ic = 0; ic < numComps; ++ic ) + integer constexpr numRecoveryPasses = 7; + integer const maxProbeIterations = LvArray::math::min( maxIterations, integer{ 25 } ); + + for( integer recoveryPass = 0; recoveryPass < numRecoveryPasses && !converged; ++recoveryPass ) + { + if( recoveryPass > 0 ) { - liquidComposition[ic] = composition[ic] / ( 1.0 + vapourPhaseMoleFraction * ( kVapourLiquid[ic] - 1.0 ) ); - vapourComposition[ic] = kVapourLiquid[ic] * liquidComposition[ic]; - } + real64 wilsonReferencePressure = pressure; - normalizeComposition( numComps, liquidComposition ); - normalizeComposition( numComps, vapourComposition ); + if( recoveryPass == 6 ) + { + // Secant / Newton-style update of the Wilson reference pressure using finite differences. + // Fugacity evaluations remain at the true cell pressure. + real64 const relPerturbation = 5.0e-4; + real64 const wPlus = pressure * ( 1.0 + relPerturbation ); + real64 const wMinus = pressure * ( 1.0 - relPerturbation ); + + real64 const errPlus = negative_two_phase_flash_internal::probeEquilibriumResidual( + numComps, wPlus, pressure, temperature, composition, componentProperties, flashData, + presentComponents, flashTolerance, maxProbeIterations, + logLiquidFugacity, logVapourFugacity, fugacityRatios, + liquidComposition, vapourComposition ); + + real64 const errMinus = negative_two_phase_flash_internal::probeEquilibriumResidual( + numComps, wMinus, pressure, temperature, composition, componentProperties, flashData, + presentComponents, flashTolerance, maxProbeIterations, + logLiquidFugacity, logVapourFugacity, fugacityRatios, + liquidComposition, vapourComposition ); + + real64 const denom = errPlus - errMinus; + wilsonReferencePressure = pressure; + if( haveErrAfterWilsonAtPressure && + LvArray::math::abs( denom ) > flashTolerance * flashTolerance + pressure * LvArray::NumericLimits< real64 >::epsilon ) + { + wilsonReferencePressure = pressure - errAfterWilsonAtPressure * ( wPlus - wMinus ) / denom; + } + wilsonReferencePressure = LvArray::math::min( LvArray::math::max( wilsonReferencePressure, 0.85 * pressure ), 1.15 * pressure ); + } + else if( recoveryPass == 2 ) + { + wilsonReferencePressure = pressure * 0.999; + } + else if( recoveryPass == 3 ) + { + wilsonReferencePressure = pressure * 1.001; + } + else if( recoveryPass == 4 ) + { + wilsonReferencePressure = pressure * 0.99; + } + else if( recoveryPass == 5 ) + { + wilsonReferencePressure = pressure * 1.01; + } - FugacityCalculator::computeLogFugacity( numComps, - pressure, - temperature, - liquidComposition.toSliceConst(), - componentProperties, - flashData.liquidEos, - flashData, - logLiquidFugacity ); - FugacityCalculator::computeLogFugacity( numComps, - pressure, - temperature, - vapourComposition.toSliceConst(), - componentProperties, - flashData.vapourEos, - flashData, - logVapourFugacity ); + if( flashData.liquidEos == EquationOfStateType::SoreideWhitson ) + { + KValueInitialization::computeSoreideWhitsonKvalue( numComps, + wilsonReferencePressure, + temperature, + componentProperties, + composition, + kVapourLiquid ); + } + else + { + KValueInitialization::computeWilsonGasLiquidKvalue( numComps, + wilsonReferencePressure, + temperature, + componentProperties, + kVapourLiquid ); + } - // Compute fugacity ratios and calculate the error - real64 error = 0.0; - for( integer const ic : presentComponents ) - { - fugacityRatios[ic] = ( logLiquidFugacity[ic] - logVapourFugacity[ic] ) + log( liquidComposition[ic] ) - log( vapourComposition[ic] ); - error += (fugacityRatios[ic]*fugacityRatios[ic]); + for( integer ic = 0; ic < numComps; ++ic ) + { + liquidComposition[ic] = composition[ic]; + vapourComposition[ic] = composition[ic]; + } } - error = LvArray::math::sqrt( error ); - - // Compute fugacity ratios and check convergence - converged = (error < flashTolerance); - if( converged ) + converged = false; + for( localIndex iterationCount = 0; iterationCount < maxIterations; ++iterationCount ) { - break; + // Solve Rachford-Rice Equation + vapourPhaseMoleFraction = RachfordRice::solve( kVapourLiquid.toSliceConst(), composition, presentComponents ); + + // Assign phase compositions + for( integer ic = 0; ic < numComps; ++ic ) + { + liquidComposition[ic] = composition[ic] / ( 1.0 + vapourPhaseMoleFraction * ( kVapourLiquid[ic] - 1.0 ) ); + vapourComposition[ic] = kVapourLiquid[ic] * liquidComposition[ic]; + } + + normalizeComposition( numComps, liquidComposition ); + normalizeComposition( numComps, vapourComposition ); + + FugacityCalculator::computeLogFugacity( numComps, + pressure, + temperature, + liquidComposition.toSliceConst(), + componentProperties, + flashData.liquidEos, + flashData, + logLiquidFugacity ); + FugacityCalculator::computeLogFugacity( numComps, + pressure, + temperature, + vapourComposition.toSliceConst(), + componentProperties, + flashData.vapourEos, + flashData, + logVapourFugacity ); + + // Compute fugacity ratios and calculate the error + real64 error = 0.0; + for( integer const ic : presentComponents ) + { + fugacityRatios[ic] = ( logLiquidFugacity[ic] - logVapourFugacity[ic] ) + log( liquidComposition[ic] ) - log( vapourComposition[ic] ); + error += (fugacityRatios[ic]*fugacityRatios[ic]); + } + error = LvArray::math::sqrt( error ); + lastInnerError = error; + + // Compute fugacity ratios and check convergence + converged = (error < flashTolerance); + + if( converged ) + { + break; + } + + // Update K-values + for( integer ic = 0; ic < numComps; ++ic ) + { + kVapourLiquid[ic] *= exp( fugacityRatios[ic] ); + } } - // Update K-values - for( integer ic = 0; ic < numComps; ++ic ) + if( recoveryPass == 1 && !converged ) { - kVapourLiquid[ic] *= exp( fugacityRatios[ic] ); + errAfterWilsonAtPressure = lastInnerError; + haveErrAfterWilsonAtPressure = true; } } diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp index e52ab6a878e..533508bbb92 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp @@ -38,6 +38,11 @@ ThermalCompressibleSinglePhaseFluid::ThermalCompressibleSinglePhaseFluid( string setInputFlag( InputFlags::OPTIONAL ). setDescription( "Fluid thermal expansion coefficient. Unit: 1/K" ); + registerWrapper( viewKeyStruct::viscosityExpansivityString(), &m_viscosityExpansivity ). + setApplyDefaultValue( 0.0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Fluid viscosity thermal expansion coefficient at the reference temperature. Unit: 1/K" ); + registerWrapper( viewKeyStruct::specificHeatCapacityString(), &m_specificHeatCapacity ). setApplyDefaultValue( 0.0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -80,6 +85,7 @@ void ThermalCompressibleSinglePhaseFluid::postInputInitialization() }; checkNonnegative( m_thermalExpansionCoeff, viewKeyStruct::thermalExpansionCoeffString() ); + checkNonnegative( m_viscosityExpansivity, viewKeyStruct::viscosityExpansivityString() ); checkNonnegative( m_specificHeatCapacity, viewKeyStruct::specificHeatCapacityString() ); checkNonnegative( m_referenceInternalEnergy, viewKeyStruct::referenceInternalEnergyString() ); @@ -98,7 +104,7 @@ ThermalCompressibleSinglePhaseFluid::KernelWrapper ThermalCompressibleSinglePhaseFluid::createKernelWrapper() { return KernelWrapper( KernelWrapper::DensRelationType( m_referencePressure, m_referenceTemperature, m_referenceDensity, m_compressibility, -m_thermalExpansionCoeff ), - KernelWrapper::ViscRelationType( m_referencePressure, m_referenceViscosity, m_viscosibility ), + KernelWrapper::ViscRelationType( m_referencePressure, m_referenceTemperature, m_referenceViscosity, m_viscosibility, -m_viscosityExpansivity ), KernelWrapper::IntEnergyRelationType( m_referenceTemperature, m_referenceInternalEnergy, m_specificHeatCapacity/m_referenceInternalEnergy ), m_density.value, m_density.derivs, diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp index 2f110c295ef..ee9e4610ac1 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp @@ -43,7 +43,7 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate public: using DensRelationType = ExponentialRelation< real64, DENS_EAT, 3 >; - using ViscRelationType = ExponentialRelation< real64, VISC_EAT >; + using ViscRelationType = ExponentialRelation< real64, VISC_EAT, 3 >; using IntEnergyRelationType = ExponentialRelation< real64, INTENERGY_EAT >; using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; @@ -114,8 +114,7 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate real64 & dEnthalpy_dPressure, real64 & dEnthalpy_dTemperature ) const override { - m_viscRelation.compute( pressure, viscosity, dViscosity_dPressure ); - dViscosity_dTemperature = 0.0; + m_viscRelation.compute( pressure, temperature, viscosity, dViscosity_dPressure, dViscosity_dTemperature ); m_densRelation.compute( pressure, temperature, density, dDensity_dPressure, dDensity_dTemperature ); @@ -240,6 +239,7 @@ class ThermalCompressibleSinglePhaseFluid : public CompressibleSinglePhaseFluid struct viewKeyStruct : public CompressibleSinglePhaseFluid::viewKeyStruct { static constexpr char const * thermalExpansionCoeffString() { return "thermalExpansionCoeff"; } + static constexpr char const * viscosityExpansivityString() { return "viscosityExpansivity"; } static constexpr char const * specificHeatCapacityString() { return "specificHeatCapacity"; } static constexpr char const * referenceTemperatureString() { return "referenceTemperature"; } static constexpr char const * referenceInternalEnergyString() { return "referenceInternalEnergy"; } @@ -255,6 +255,9 @@ class ThermalCompressibleSinglePhaseFluid : public CompressibleSinglePhaseFluid /// scalar fluid thermal expansion coefficient real64 m_thermalExpansionCoeff; + /// scalar fluid viscosity thermal expansion coefficient + real64 m_viscosityExpansivity; + /// scalar fluid volumetric heat capacity coefficient real64 m_specificHeatCapacity; diff --git a/src/coreComponents/constitutive/relativePermeability/TableRelativePermeabilityHelpers.cpp b/src/coreComponents/constitutive/relativePermeability/TableRelativePermeabilityHelpers.cpp index a93ab87c8a8..6148bd37bf0 100644 --- a/src/coreComponents/constitutive/relativePermeability/TableRelativePermeabilityHelpers.cpp +++ b/src/coreComponents/constitutive/relativePermeability/TableRelativePermeabilityHelpers.cpp @@ -68,10 +68,10 @@ TableRelativePermeabilityHelpers::validateRelativePermeabilityTable( TableFuncti // note that the TableFunction class has already checked that the coordinates are monotone // check phase relative permeability - GEOS_THROW_IF( !isZero( relPerm[i] ) && (relPerm[i] - relPerm[i-1]) < 1e-15, - GEOS_FMT( "TableFunction '{}' values must be strictly increasing (|Delta kr| > 1e-15 between two non-zero values)", - fullConstitutiveName ), - InputError, relPermTable.getDataContext() ); + // GEOS_THROW_IF( !isZero( relPerm[i] ) && (relPerm[i] - relPerm[i-1]) < 1e-15, + // GEOS_FMT( "TableFunction '{}' values must be strictly increasing (|Delta kr| > 1e-15 between two non-zero values)", + // fullConstitutiveName ), + // InputError, relPermTable.getDataContext() ); if( isZero( relPerm[i-1] ) && !isZero( relPerm[i] ) ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/GlobalComponentFractionKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/GlobalComponentFractionKernel.hpp index c5317e3a659..b00b0f0de50 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/GlobalComponentFractionKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/GlobalComponentFractionKernel.hpp @@ -28,6 +28,9 @@ namespace geos namespace isothermalCompositionalMultiphaseBaseKernels { +/// Minimum |sum(global component density)| for forming global component fractions; matches minDensForDivision in AccumulationKernel.hpp. +static constexpr real64 minTotalDensityForGlobalCompFraction = 1e-10; + /******************************** GlobalComponentFractionKernel ********************************/ /** @@ -77,16 +80,32 @@ class GlobalComponentFractionKernel : public PropertyKernelBase< NUM_COMP > totalDensity += compDens[ic]; } - real64 const totalDensityInv = 1.0 / totalDensity; - - for( integer ic = 0; ic < numComp; ++ic ) + real64 const absTotalDensity = totalDensity >= 0.0 ? totalDensity : -totalDensity; + if( absTotalDensity < minTotalDensityForGlobalCompFraction ) + { + real64 const invN = 1.0 / numComp; + for( integer ic = 0; ic < numComp; ++ic ) + { + compFrac[ic] = invN; + for( integer jc = 0; jc < numComp; ++jc ) + { + dCompFrac_dCompDens[ic][jc] = 0.0; + } + } + } + else { - compFrac[ic] = compDens[ic] * totalDensityInv; - for( integer jc = 0; jc < numComp; ++jc ) + real64 const totalDensityInv = 1.0 / totalDensity; + + for( integer ic = 0; ic < numComp; ++ic ) { - dCompFrac_dCompDens[ic][jc] = -compFrac[ic] * totalDensityInv; + compFrac[ic] = compDens[ic] * totalDensityInv; + for( integer jc = 0; jc < numComp; ++jc ) + { + dCompFrac_dCompDens[ic][jc] = -compFrac[ic] * totalDensityInv; + } + dCompFrac_dCompDens[ic][ic] += totalDensityInv; } - dCompFrac_dCompDens[ic][ic] += totalDensityInv; } compFractionKernelOp( compFrac, dCompFrac_dCompDens ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp index acf38c45ce1..7ee033b97ca 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp @@ -64,7 +64,7 @@ class MultiphasePoromechanics : /// maxNumTrialSupportPointPerElem by definition. When the FE_TYPE is not a Virtual Element, this /// will be the actual number of nodes per element. static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; - static constexpr int maxNumComponents = 3; + static constexpr int maxNumComponents = 4; using Base::numDofPerTestSupportPoint; using Base::numDofPerTrialSupportPoint; using Base::m_dofNumber; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp index 8aa3b74379a..99f5ea53dbb 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp @@ -62,7 +62,7 @@ class ThermalMultiphasePoromechanics : /// will be the actual number of nodes per element. static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; /// Maximum number of fluid components - static constexpr int maxNumComponents = 3; + static constexpr int maxNumComponents = 4; using Base::numDofPerTestSupportPoint; using Base::numDofPerTrialSupportPoint; using Base::m_dofNumber;