diff --git a/integratedTests b/integratedTests index c03ad98d8ca..d2f67012021 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit c03ad98d8cad213d468110e64849f79033794cfa +Subproject commit d2f67012021172f3b9d9c3cdd831c5151c37faaf diff --git a/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp b/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp index 566909fb065..fedd82e6ff7 100644 --- a/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp +++ b/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp @@ -28,17 +28,13 @@ using namespace stringutilities; namespace PVTProps { +constexpr real64 minForDivision = 1e-10; + constexpr real64 T_K_f = 273.15; constexpr real64 P_Pa_f = 1e+5; - constexpr real64 P_c = 73.773 * P_Pa_f; constexpr real64 T_c = 304.1282; -//constexpr real64 rho_c = 467.6; - -//constexpr real64 R = 188.9241; - constexpr real64 Rgas = 8.314467; - constexpr real64 V_c = Rgas*T_c/P_c; constexpr real64 acoef[] = @@ -335,7 +331,6 @@ void CO2SolubilityFunction::Partition( EvalVarArgs const & pressure, EvalVarArgs solubility = m_CO2SolubilityTable->Value( P, T ); real64 const waterMW = m_componentMolarWeight[m_waterIndex]; - // real64 CO2MW = m_componentMolarWeight[m_CO2Index]; solubility *= waterMW; @@ -346,7 +341,14 @@ void CO2SolubilityFunction::Partition( EvalVarArgs const & pressure, EvalVarArgs //Y = C/W = z/(1-z) - Y = compFraction[m_CO2Index] / (1.0 - compFraction[m_CO2Index]); + if( compFraction[m_CO2Index].m_var > 1.0 - minForDivision ) + { + Y = compFraction[m_CO2Index] / minForDivision; + } + else + { + Y = compFraction[m_CO2Index] / (1.0 - compFraction[m_CO2Index]); + } if( Y < X ) { @@ -356,7 +358,9 @@ void CO2SolubilityFunction::Partition( EvalVarArgs const & pressure, EvalVarArgs phaseFraction[m_phaseGasIndex] = 0.0; for( localIndex c = 0; c < m_componentNames.size(); ++c ) + { phaseCompFraction[m_phaseLiquidIndex][c] = compFraction[c]; + } } else diff --git a/src/coreComponents/constitutive/relativePermeability/BrooksCoreyRelativePermeability.hpp b/src/coreComponents/constitutive/relativePermeability/BrooksCoreyRelativePermeability.hpp index e7c5457ba82..16c44197ce9 100644 --- a/src/coreComponents/constitutive/relativePermeability/BrooksCoreyRelativePermeability.hpp +++ b/src/coreComponents/constitutive/relativePermeability/BrooksCoreyRelativePermeability.hpp @@ -173,7 +173,7 @@ BrooksCoreyRelativePermeabilityUpdate:: } else { - phaseRelPerm[ip] = (satScaled < 0.0) ? 0.0 : scale; + phaseRelPerm[ip] = (satScaled <= 0.0) ? 0.0 : scale; } } } diff --git a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst index 320d28be1f7..375fc56129a 100644 --- a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst +++ b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst @@ -1,24 +1,26 @@ -========================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== -Name Type Default Description -========================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== -capPressureNames string_array {} Name of the capillary pressure constitutive model to use -cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] -discretization string required Name of discretization object to use for this solver. -fluidNames string_array required Names of fluid constitutive models for each region. -initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. -inputFluxEstimate real64 1 Initial estimate of the input flux used only for residual scaling. This should be essentially equivalent to the input flux * dt. -logLevel integer 0 Log level -meanPermCoeff real64 1 Coefficient to move between harmonic mean (1.0) and arithmetic mean (0.0) for the calculation of permeability between elements. -name string required A name is required for any non-unique nodes -relPermNames string_array required Name of the relative permeability constitutive model to use -solidNames string_array required Names of solid constitutive models for each region. -targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. -temperature real64 required Temperature -useMass integer 0 Use mass formulation instead of molar -LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` -NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` -========================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== +============================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +============================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== +allowLocalCompDensityChopping integer 1 Flag indicating whether local (cell-wise) chopping of negative compositions is allowed +capPressureNames string_array {} Name of the capillary pressure constitutive model to use +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +discretization string required Name of discretization object to use for this solver. +fluidNames string_array required Names of fluid constitutive models for each region. +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +inputFluxEstimate real64 1 Initial estimate of the input flux used only for residual scaling. This should be essentially equivalent to the input flux * dt. +logLevel integer 0 Log level +maxCompFractionChange real64 1 Maximum (absolute) change in a component fraction between two Newton iterations +meanPermCoeff real64 1 Coefficient to move between harmonic mean (1.0) and arithmetic mean (0.0) for the calculation of permeability between elements. +name string required A name is required for any non-unique nodes +relPermNames string_array required Name of the relative permeability constitutive model to use +solidNames string_array required Names of solid constitutive models for each region. +targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +temperature real64 required Temperature +useMass integer 0 Use mass formulation instead of molar +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +============================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseWell.rst b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseWell.rst index 0c7bda74961..a0dca92ba59 100644 --- a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseWell.rst +++ b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseWell.rst @@ -1,21 +1,23 @@ -========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== -Name Type Default Description -========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== -cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] -discretization string none Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -fluidNames string_array required Name of fluid constitutive object to use for this solver. -initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. -logLevel integer 0 Log level -name string required A name is required for any non-unique nodes -relPermNames string_array required Names of relative permeability constitutive models to use -targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. -useMass integer 0 Use mass formulation instead of molar -wellTemperature real64 required Temperature -LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` -NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` -WellControls node :ref:`XML_WellControls` -========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== +============================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +============================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== +allowLocalCompDensityChopping integer 1 Flag indicating whether local (cell-wise) chopping of negative compositions is allowed +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +discretization string none Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +fluidNames string_array required Name of fluid constitutive object to use for this solver. +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +logLevel integer 0 Log level +maxCompFractionChange real64 1 Maximum (absolute) change in a component fraction between two Newton iterations +name string required A name is required for any non-unique nodes +relPermNames string_array required Names of relative permeability constitutive models to use +targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +useMass integer 0 Use mass formulation instead of molar +wellTemperature real64 required Temperature +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +WellControls node :ref:`XML_WellControls` +============================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/fileIO/schema/schema.xsd b/src/coreComponents/fileIO/schema/schema.xsd index eaca798fe43..1b09986f4e2 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd +++ b/src/coreComponents/fileIO/schema/schema.xsd @@ -761,6 +761,8 @@ + + @@ -775,6 +777,8 @@ + + @@ -818,6 +822,8 @@ + + @@ -828,6 +834,8 @@ + + diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp index a9124673b32..5dd2c703dab 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp @@ -41,13 +41,18 @@ using namespace dataRepository; using namespace constitutive; using namespace CompositionalMultiphaseFlowKernels; +static constexpr real64 minDensForDivision = 1e-10; + CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, Group * const parent ) : FlowSolverBase( name, parent ), m_numPhases( 0 ), m_numComponents( 0 ), - m_capPressureFlag( 0 ) + m_capPressureFlag( 0 ), + m_maxCompFracChange( 1.0 ), + m_minScalingFactor( 0.01 ), + m_allowCompDensChopping( 1 ) { //START_SPHINX_INCLUDE_00 this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> @@ -69,6 +74,18 @@ CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, setInputFlag( InputFlags::OPTIONAL )-> setDescription( "Name of the capillary pressure constitutive model to use" ); + this->registerWrapper( viewKeyStruct::maxCompFracChangeString, &m_maxCompFracChange )-> + setSizedFromParent( 0 )-> + setInputFlag( InputFlags::OPTIONAL )-> + setApplyDefaultValue( 1.0 )-> + setDescription( "Maximum (absolute) change in a component fraction between two Newton iterations" ); + + this->registerWrapper( viewKeyStruct::allowLocalCompDensChoppingString, &m_allowCompDensChopping )-> + setSizedFromParent( 0 )-> + setInputFlag( InputFlags::OPTIONAL )-> + setApplyDefaultValue( 1 )-> + setDescription( "Flag indicating whether local (cell-wise) chopping of negative compositions is allowed" ); + m_linearSolverParameters.get().mgr.strategy = "CompositionalMultiphaseFlow"; } @@ -78,6 +95,11 @@ void CompositionalMultiphaseFlow::PostProcessInput() FlowSolverBase::PostProcessInput(); CheckModelNames( m_relPermModelNames, viewKeyStruct::relPermNamesString ); m_capPressureFlag = CheckModelNames( m_capPressureModelNames, viewKeyStruct::capPressureNamesString, true ); + + GEOSX_ERROR_IF_GT_MSG( m_maxCompFracChange, 1.0, + "The maximum absolute change in component fraction must smaller or equal to 1.0" ); + GEOSX_ERROR_IF_LT_MSG( m_maxCompFracChange, 0.0, + "The maximum absolute change in component fraction must larger or equal to 0.0" ); } void CompositionalMultiphaseFlow::RegisterDataOnMesh( Group * const MeshBodies ) @@ -1199,12 +1221,92 @@ void CompositionalMultiphaseFlow::SolveSystem( DofManager const & dofManager, SolverBase::SolveSystem( dofManager, matrix, rhs, solution ); } +real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) +{ + GEOSX_MARK_FUNCTION; + + // check if we want to rescale the Newton update + if( m_maxCompFracChange >= 1.0 ) + { + // no rescaling wanted, we just return 1.0; + return 1.0; + } + + real64 constexpr eps = minDensForDivision; + real64 const maxCompFracChange = m_maxCompFracChange; + + localIndex const NC = m_numComponents; + + MeshLevel const & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); + + globalIndex const rankOffset = dofManager.rankOffset(); + string const dofKey = dofManager.getKey( viewKeyStruct::dofFieldString ); + real64 scalingFactor = 1.0; + + forTargetSubRegions( mesh, [&]( localIndex const, ElementSubRegionBase const & subRegion ) + { + arrayView1d< globalIndex const > const & dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const & elemGhostRank = subRegion.ghostRank(); + + arrayView2d< real64 const > const & compDens = subRegion.getReference< array2d< real64 > >( viewKeyStruct::globalCompDensityString ); + arrayView2d< real64 const > const & dCompDens = subRegion.getReference< array2d< real64 > >( viewKeyStruct::deltaGlobalCompDensityString ); + + RAJA::ReduceMin< parallelDeviceReduce, real64 > minVal( 1.0 ); + + forAll< parallelDevicePolicy<> >( dofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) + { + if( elemGhostRank[ei] < 0 ) + { + real64 prevTotalDens = 0; + for( localIndex ic = 0; ic < NC; ++ic ) + { + prevTotalDens += compDens[ei][ic] + dCompDens[ei][ic]; + } + + // compute the change in component densities and component fractions + for( localIndex ic = 0; ic < NC; ++ic ) + { + localIndex const lid = dofNumber[ei] + ic + 1 - rankOffset; + + // compute scaling factor based on relative change in component densities + real64 const absCompDensChange = fabs( localSolution[lid] ); + real64 const maxAbsCompDensChange = maxCompFracChange * prevTotalDens; + + // This actually checks the change in component fraction, using a lagged total density + // Indeed we can rewrite the following check as: + // | prevCompDens / prevTotalDens - newCompDens / prevTotalDens | > maxCompFracChange + // Note that the total density in the second term is lagged (i.e, we use prevTotalDens) + // because I found it more robust than using directly newTotalDens (which can vary also + // wildly when the compDens change is large) + if( absCompDensChange > maxAbsCompDensChange && absCompDensChange > eps ) + { + minVal.min( maxAbsCompDensChange / absCompDensChange ); + } + } + } + } ); + + if( minVal.get() < scalingFactor ) + { + scalingFactor = minVal.get(); + } + } ); + + return LvArray::math::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); +} + + bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & domain, DofManager const & dofManager, arrayView1d< real64 const > const & localSolution, real64 const scalingFactor ) { + real64 constexpr eps = minDensForDivision; + localIndex const NC = m_numComponents; + integer const allowCompDensChopping = m_allowCompDensChopping; MeshLevel const & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); @@ -1233,10 +1335,27 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d real64 const newPres = pres[ei] + dPres[ei] + scalingFactor * localSolution[localRow]; check.min( newPres >= 0.0 ); } - for( localIndex ic = 0; ic < NC; ++ic ) + + // if component density chopping is not allowed, the time step fails if a component density is negative + // otherwise, we just check that the total density is positive, and negative component densities + // will be chopped (i.e., set to zero) in ApplySystemSolution) + if( !allowCompDensChopping ) { - real64 const newDens = compDens[ei][ic] + dCompDens[ei][ic] + scalingFactor * localSolution[localRow + ic + 1]; - check.min( newDens >= 0.0 ); + for( localIndex ic = 0; ic < NC; ++ic ) + { + real64 const newDens = compDens[ei][ic] + dCompDens[ei][ic] + scalingFactor * localSolution[localRow + ic + 1]; + check.min( newDens >= 0.0 ); + } + } + else + { + real64 totalDens = 0.0; + for( localIndex ic = 0; ic < NC; ++ic ) + { + real64 const newDens = compDens[ei][ic] + dCompDens[ei][ic] + scalingFactor * localSolution[localRow + ic + 1]; + totalDens += (newDens > 0.0) ? newDens : 0.0; + } + check.min( totalDens >= eps ); } } } ); @@ -1253,7 +1372,6 @@ void CompositionalMultiphaseFlow::ApplySystemSolution( DofManager const & dofMan DomainPartition & domain ) { MeshLevel & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); - dofManager.addVectorToField( localSolution, viewKeyStruct::dofFieldString, viewKeyStruct::deltaPressureString, @@ -1266,6 +1384,13 @@ void CompositionalMultiphaseFlow::ApplySystemSolution( DofManager const & dofMan scalingFactor, 1, m_numDofPerCell ); + // if component density chopping is allowed, some component densities may be negative after the update + // these negative component densities are set to zero in this function + if( m_allowCompDensChopping ) + { + ChopNegativeDensities( domain ); + } + std::map< string, string_array > fieldNames; fieldNames["elems"].emplace_back( string( viewKeyStruct::deltaPressureString ) ); fieldNames["elems"].emplace_back( string( viewKeyStruct::deltaGlobalCompDensityString ) ); @@ -1277,6 +1402,38 @@ void CompositionalMultiphaseFlow::ApplySystemSolution( DofManager const & dofMan } ); } +void CompositionalMultiphaseFlow::ChopNegativeDensities( DomainPartition & domain ) +{ + MeshLevel & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); + + localIndex const NC = m_numComponents; + forTargetSubRegions( mesh, [&]( localIndex const, ElementSubRegionBase & subRegion ) + { + arrayView1d< integer const > const & ghostRank = + subRegion.getReference< array1d< integer > >( ObjectManagerBase::viewKeyStruct::ghostRankString ); + + arrayView2d< real64 const > const & compDens = + subRegion.getReference< array2d< real64 > >( viewKeyStruct::globalCompDensityString ); + arrayView2d< real64 > const & dCompDens = + subRegion.getReference< array2d< real64 > >( viewKeyStruct::deltaGlobalCompDensityString ); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) + { + if( ghostRank[ei] < 0 ) + { + for( localIndex ic = 0; ic < NC; ++ic ) + { + real64 const newDens = compDens[ei][ic] + dCompDens[ei][ic]; + if( newDens < 0 ) + { + dCompDens[ei][ic] = -compDens[ei][ic]; + } + } + } + } ); + } ); +} + void CompositionalMultiphaseFlow::ResetStateToBeginningOfStep( DomainPartition & domain ) { MeshLevel & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp index bd0738f120c..9dcfbce873a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp @@ -140,6 +140,11 @@ class CompositionalMultiphaseFlow : public FlowSolverBase ParallelVector & rhs, ParallelVector & solution ) override; + virtual real64 + ScalingForSystemSolution( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) override; + virtual bool CheckSystemSolution( DomainPartition const & domain, DofManager const & dofManager, @@ -280,6 +285,9 @@ class CompositionalMultiphaseFlow : public FlowSolverBase static constexpr auto relPermNamesString = "relPermNames"; static constexpr auto capPressureNamesString = "capPressureNames"; + static constexpr auto maxCompFracChangeString = "maxCompFractionChange"; + static constexpr auto allowLocalCompDensChoppingString = "allowLocalCompDensityChopping"; + static constexpr auto facePressureString = "facePressure"; static constexpr auto bcPressureString = "bcPressure"; @@ -310,6 +318,7 @@ class CompositionalMultiphaseFlow : public FlowSolverBase static constexpr auto phaseViscosityString = "phaseViscosity"; static constexpr auto phaseRelativePermeabilityString = "phaseRelativePermeability"; static constexpr auto phaseCapillaryPressureString = "phaseCapillaryPressure"; + } viewKeysCompMultiphaseFlow; struct groupKeyStruct : SolverBase::groupKeyStruct @@ -363,6 +372,12 @@ class CompositionalMultiphaseFlow : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const; + /** + * @brief Sets all the negative component densities (if any) to zero. + * @param domain the physical domain object + */ + void ChopNegativeDensities( DomainPartition & domain ); + protected: virtual void PostProcessInput() override; @@ -393,7 +408,6 @@ class CompositionalMultiphaseFlow : public FlowSolverBase */ void ResetViews( MeshLevel & mesh ) override; - /// the max number of fluid phases localIndex m_numPhases; @@ -415,6 +429,16 @@ class CompositionalMultiphaseFlow : public FlowSolverBase /// name of the cap pressure constitutive model array1d< string > m_capPressureModelNames; + /// maximum (absolute) change in a component fraction between two Newton iterations + real64 m_maxCompFracChange; + + /// minimum value of the scaling factor obtained by enforcing maxCompFracChange + real64 m_minScalingFactor; + + /// flag indicating whether local (cell-wise) chopping of negative compositions is allowed + integer m_allowCompDensChopping; + + ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_pressure; ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_deltaPressure; 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 022aed06463..f5580054e06 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 @@ -8,7 +8,7 @@ flowSolverName="compositionalMultiphaseFlow" wellSolverName="compositionalMultiphaseWell" logLevel="1" - initialDt="1e1" + initialDt="2e3" targetRegions="{ reservoir, wellRegion1, wellRegion2, wellRegion3, wellRegion4, wellRegion5 }"> @@ -65,9 +66,9 @@ name="wellControls5" type="injector" control="liquidRate" - targetBHP="1e7" - targetRate="1e-2" - injectionStream="{ 0.1, 0.1, 0.8 }"/> + targetBHP="1e8" + targetRate="4e1" + injectionStream="{ 0.0, 0.0, 1.0 }"/> @@ -77,8 +78,8 @@ + fieldsToImport="{ PERM, PORO }" + fieldNamesInGEOSX="{ permeability, referencePorosity }"/> + maxTime="17280000"> @@ -261,7 +262,7 @@ @@ -274,13 +275,6 @@ - + scale="0.9"/> + scale="0.05"/> + scale="0.05"/> diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_baker_1d.xml b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_baker_1d.xml index 114fa34a74f..e65bff43efb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_baker_1d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseFlow/deadoil_3ph_baker_1d.xml @@ -48,29 +48,29 @@ + maxTime="2e7"> 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 b6620fe0149..32588f60d4e 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 @@ -48,29 +48,29 @@ + maxTime="2e7"> 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 ec8a2c6b5c2..e924ad5bc72 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 @@ -55,36 +55,36 @@ + maxTime="2e8"> diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index fe65b3f6e93..c336c8ea4b1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -51,7 +51,10 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, m_numPhases( 0 ), m_numComponents( 0 ), m_temperature( 0.0 ), - m_useMass( false ) + m_useMass( false ), + m_maxCompFracChange( 1.0 ), + m_minScalingFactor( 0.01 ), + m_allowCompDensChopping( 1 ) { this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> setInputFlag( InputFlags::REQUIRED )-> @@ -65,6 +68,19 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, this->registerWrapper( viewKeyStruct::relPermNamesString, &m_relPermModelNames )-> setInputFlag( InputFlags::REQUIRED )-> setDescription( "Names of relative permeability constitutive models to use" ); + + this->registerWrapper( viewKeyStruct::maxCompFracChangeString, &m_maxCompFracChange )-> + setSizedFromParent( 0 )-> + setInputFlag( InputFlags::OPTIONAL )-> + setApplyDefaultValue( 1.0 )-> + setDescription( "Maximum (absolute) change in a component fraction between two Newton iterations" ); + + this->registerWrapper( viewKeyStruct::allowLocalCompDensChoppingString, &m_allowCompDensChopping )-> + setSizedFromParent( 0 )-> + setInputFlag( InputFlags::OPTIONAL )-> + setApplyDefaultValue( 1 )-> + setDescription( "Flag indicating whether local (cell-wise) chopping of negative compositions is allowed" ); + } void CompositionalMultiphaseWell::PostProcessInput() @@ -76,6 +92,12 @@ void CompositionalMultiphaseWell::PostProcessInput() GEOSX_ERROR_IF( flowSolver == nullptr, "Flow solver " << GetFlowSolverName() << " not found or incompatible type " "(referenced from well solver " << getName() << ")" ); + + GEOSX_ERROR_IF_GT_MSG( m_maxCompFracChange, 1.0, + "The maximum absolute change in component fraction must smaller or equal to 1.0" ); + GEOSX_ERROR_IF_LT_MSG( m_maxCompFracChange, 0.0, + "The maximum absolute change in component fraction must larger or equal to 0.0" ); + } void CompositionalMultiphaseWell::RegisterDataOnMesh( Group * const meshBodies ) @@ -627,6 +649,60 @@ CompositionalMultiphaseWell::CalculateResidualNorm( DomainPartition const & doma return sqrt( MpiWrapper::Sum( localResidualNorm, MPI_COMM_GEOSX ) ); } +real64 +CompositionalMultiphaseWell::ScalingForSystemSolution( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) +{ + GEOSX_MARK_FUNCTION; + + // check if we want to rescale the Newton update + if( m_maxCompFracChange >= 1.0 ) + { + // no rescaling wanted, we just return 1.0; + return 1.0; + } + + MeshLevel const & meshLevel = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); + + real64 scalingFactor = 1.0; + forTargetSubRegions< WellElementSubRegion >( meshLevel, [&]( localIndex const, + WellElementSubRegion const & subRegion ) + { + // get the degree of freedom numbers on well elements and ghosting info + string const wellDofKey = dofManager.getKey( WellElementDofName() ); + arrayView1d< globalIndex const > const & wellElemDofNumber = + subRegion.getReference< array1d< globalIndex > >( wellDofKey ); + arrayView1d< integer const > const & wellElemGhostRank = + subRegion.getReference< array1d< integer > >( ObjectManagerBase::viewKeyStruct::ghostRankString ); + + // get a reference to the primary variables on well elements + arrayView2d< real64 const > const & wellElemCompDens = + subRegion.getReference< array2d< real64 > >( viewKeyStruct::globalCompDensityString ); + arrayView2d< real64 const > const & dWellElemCompDens = + subRegion.getReference< array2d< real64 > >( viewKeyStruct::deltaGlobalCompDensityString ); + + real64 const subRegionScalingFactor = + SolutionScalingKernel::Launch< parallelDevicePolicy<>, + parallelDeviceReduce >( localSolution, + dofManager.rankOffset(), + NumFluidComponents(), + wellElemDofNumber, + wellElemGhostRank, + wellElemCompDens, + dWellElemCompDens, + m_maxCompFracChange ); + + + if( subRegionScalingFactor < scalingFactor ) + { + scalingFactor = subRegionScalingFactor; + } + } ); + + return LvArray::math::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); +} + bool CompositionalMultiphaseWell::CheckSystemSolution( DomainPartition const & domain, DofManager const & dofManager, @@ -670,6 +746,7 @@ CompositionalMultiphaseWell::CheckSystemSolution( DomainPartition const & domain dWellElemPressure, wellElemCompDens, dWellElemCompDens, + m_allowCompDensChopping, scalingFactor ); if( subRegionSolutionCheck == 0 ) @@ -776,6 +853,7 @@ CompositionalMultiphaseWell::ApplySystemSolution( DofManager const & dofManager, real64 const scalingFactor, DomainPartition & domain ) { + // update all the fields using the global damping coefficients dofManager.addVectorToField( localSolution, WellElementDofName(), viewKeyStruct::deltaPressureString, @@ -794,6 +872,14 @@ CompositionalMultiphaseWell::ApplySystemSolution( DofManager const & dofManager, scalingFactor, m_numDofPerWellElement - 1, m_numDofPerWellElement ); + // if component density chopping is allowed, some component densities may be negative after the update + // these negative component densities are set to zero in this function + if( m_allowCompDensChopping ) + { + ChopNegativeDensities( domain ); + } + + // synchronize std::map< string, string_array > fieldNames; fieldNames["elems"].emplace_back( string( viewKeyStruct::deltaPressureString ) ); fieldNames["elems"].emplace_back( string( viewKeyStruct::deltaGlobalCompDensityString ) ); @@ -808,6 +894,43 @@ CompositionalMultiphaseWell::ApplySystemSolution( DofManager const & dofManager, } +void CompositionalMultiphaseWell::ChopNegativeDensities( DomainPartition & domain ) +{ + MeshLevel & meshLevel = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); + localIndex const NC = m_numComponents; + + forTargetSubRegions< WellElementSubRegion >( meshLevel, [&]( localIndex const, + WellElementSubRegion & subRegion ) + { + arrayView1d< integer const > const & wellElemGhostRank = + subRegion.getReference< array1d< integer > >( ObjectManagerBase::viewKeyStruct::ghostRankString ); + + arrayView2d< real64 const > const & wellElemCompDens = + subRegion.getReference< array2d< real64 > >( viewKeyStruct::globalCompDensityString ); + arrayView2d< real64 > const & dWellElemCompDens = + subRegion.getReference< array2d< real64 > >( viewKeyStruct::deltaGlobalCompDensityString ); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) + { + if( wellElemGhostRank[iwelem] < 0 ) + { + for( localIndex ic = 0; ic < NC; ++ic ) + { + // get the latest component density (i.e., after the update of dWellElemCompDens) + real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic]; + // we allowed for some densities to be slightly negative in CheckSystemSolution + // if the new density is negative, chop back to zero + if( newDens < 0 ) + { + dWellElemCompDens[iwelem][ic] = -wellElemCompDens[iwelem][ic]; + } + } + } + } ); + } ); +} + + void CompositionalMultiphaseWell::ResetStateToBeginningOfStep( DomainPartition & domain ) { MeshLevel & meshLevel = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index 00b43e1a535..c8f5cc675cb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -111,6 +111,11 @@ class CompositionalMultiphaseWell : public WellSolverBase DofManager const & dofManager, arrayView1d< real64 const > const & localRhs ) override; + virtual real64 + ScalingForSystemSolution( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) override; + virtual bool CheckSystemSolution( DomainPartition const & domain, DofManager const & dofManager, @@ -210,6 +215,12 @@ class CompositionalMultiphaseWell : public WellSolverBase arrayView1d< real64 > const & localRhs ) override; + /** + * @brief Sets all the negative component densities (if any) to zero. + * @param domain the physical domain object + */ + void ChopNegativeDensities( DomainPartition & domain ); + arrayView1d< string const > const & relPermModelNames() const { return m_relPermModelNames; } struct viewKeyStruct : WellSolverBase::viewKeyStruct @@ -218,9 +229,12 @@ class CompositionalMultiphaseWell : public WellSolverBase // inputs static constexpr auto temperatureString = "wellTemperature"; - static constexpr auto useMassFlagString = "useMass"; + static constexpr auto useMassFlagString = CompositionalMultiphaseFlow::viewKeyStruct::useMassFlagString; - static constexpr auto relPermNamesString = "relPermNames"; + static constexpr auto relPermNamesString = CompositionalMultiphaseFlow::viewKeyStruct::relPermNamesString; + + static constexpr auto maxCompFracChangeString = CompositionalMultiphaseFlow::viewKeyStruct::maxCompFracChangeString; + static constexpr auto allowLocalCompDensChoppingString = CompositionalMultiphaseFlow::viewKeyStruct::allowLocalCompDensChoppingString; // primary solution field static constexpr auto pressureString = CompositionalMultiphaseFlow::viewKeyStruct::pressureString; @@ -236,11 +250,6 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr auto dPhaseVolumeFraction_dGlobalCompDensityString = CompositionalMultiphaseFlow::viewKeyStruct::dPhaseVolumeFraction_dGlobalCompDensityString; - // mixture density - static constexpr auto mixtureDensityString = "wellElementMixtureDensity"; - static constexpr auto dMixtureDensity_dPressureString = "dWellElementMixtureDensity_dPres"; - static constexpr auto dMixtureDensity_dGlobalCompDensityString = "dWellElementMixtureDensity_dComp"; - // global component fractions static constexpr auto globalCompFractionString = CompositionalMultiphaseFlow::viewKeyStruct::globalCompFractionString; static constexpr auto dGlobalCompFraction_dGlobalCompDensityString = @@ -250,6 +259,7 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr auto compPerforationRateString = "compPerforationRate"; static constexpr auto dCompPerforationRate_dPresString = "dCompPerforationRate_dPres"; static constexpr auto dCompPerforationRate_dCompString = "dCompPerforationRate_dComp"; + } viewKeysCompMultiphaseWell; struct groupKeyStruct : SolverBase::groupKeyStruct @@ -305,12 +315,6 @@ class CompositionalMultiphaseWell : public WellSolverBase */ void ResizeFields( WellElementSubRegion & subRegion ); - /** - * @brief Save all the rates and pressures in the well for reporting purposes - * @param well the well with its perforations - */ - void RecordWellData( WellElementSubRegion const & subRegion ); - /// the max number of fluid phases localIndex m_numPhases; @@ -326,6 +330,14 @@ class CompositionalMultiphaseWell : public WellSolverBase /// list of relative permeability model names per target region array1d< string > m_relPermModelNames; + /// maximum (absolute) change in a component fraction between two Newton iterations + real64 m_maxCompFracChange; + + /// minimum value of the scaling factor obtained by enforcing maxCompFracChange + real64 m_minScalingFactor; + + /// flag indicating whether local (cell-wise) chopping of negative compositions is allowed + integer m_allowCompDensChopping; /// views into reservoir primary variable fields diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index ca623864525..c4314e31419 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -31,6 +31,8 @@ namespace geosx namespace CompositionalMultiphaseWellKernels { +static constexpr real64 minDensForDivision = 1e-10; + /******************************** ControlEquationHelper ********************************/ struct ControlEquationHelper @@ -1209,6 +1211,63 @@ struct ResidualNormKernel }; +/******************************** SolutionScalingKernel ********************************/ + +struct SolutionScalingKernel +{ + template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > + static real64 + Launch( LOCAL_VECTOR const localSolution, + globalIndex const rankOffset, + localIndex const numComponents, + arrayView1d< globalIndex const > const & wellElemDofNumber, + arrayView1d< integer const > const & wellElemGhostRank, + arrayView2d< real64 const > const & wellElemCompDens, + arrayView2d< real64 const > const & dWellElemCompDens, + real64 const maxCompFracChange ) + { + real64 constexpr eps = minDensForDivision; + + RAJA::ReduceMin< REDUCE_POLICY, real64 > minVal( 1.0 ); + + forAll< POLICY >( wellElemDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) + { + if( wellElemGhostRank[iwelem] < 0 ) + { + + real64 prevTotalDens = 0; + for( localIndex ic = 0; ic < numComponents; ++ic ) + { + prevTotalDens += wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic]; + } + + for( localIndex ic = 0; ic < numComponents; ++ic ) + { + localIndex const lid = wellElemDofNumber[iwelem] + ic + 1 - rankOffset; + + // compute scaling factor based on relative change in component densities + real64 const absCompDensChange = fabs( localSolution[lid] ); + real64 const maxAbsCompDensChange = maxCompFracChange * prevTotalDens; + + // This actually checks the change in component fraction, using a lagged total density + // Indeed we can rewrite the following check as: + // | prevCompDens / prevTotalDens - newCompDens / prevTotalDens | > maxCompFracChange + // Note that the total density in the second term is lagged (i.e, we use prevTotalDens) + // because I found it more robust than using directly newTotalDens (which can vary also + // wildly when the compDens change is large) + if( absCompDensChange > maxAbsCompDensChange && absCompDensChange > eps ) + { + minVal.min( maxAbsCompDensChange / absCompDensChange ); + } + } + } + } ); + return minVal.get(); + } + +}; + + /******************************** SolutionCheckKernel ********************************/ struct SolutionCheckKernel @@ -1224,8 +1283,11 @@ struct SolutionCheckKernel arrayView1d< real64 const > const & dWellElemPressure, arrayView2d< real64 const > const & wellElemCompDens, arrayView2d< real64 const > const & dWellElemCompDens, + integer const allowCompDensChopping, real64 const scalingFactor ) { + real64 constexpr eps = minDensForDivision; + RAJA::ReduceMin< REDUCE_POLICY, localIndex > minVal( 1 ); forAll< POLICY >( wellElemDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) @@ -1236,18 +1298,41 @@ struct SolutionCheckKernel localIndex lid = wellElemDofNumber[iwelem] + CompositionalMultiphaseWell::ColOffset::DPRES - rankOffset; real64 const newPres = wellElemPressure[iwelem] + dWellElemPressure[iwelem] + scalingFactor * localSolution[lid]; + + // the pressure must be positive if( newPres < 0.0 ) { minVal.min( 0 ); } - // comp densities - for( localIndex ic = 0; ic < numComponents; ++ic ) + // if component density is not allowed, the time step fails if a component density is negative + // otherwise, we just check that the total density is positive, and negative component densities + // will be chopped (i.e., set to zero) in ApplySystemSolution + if( !allowCompDensChopping ) + { + for( localIndex ic = 0; ic < numComponents; ++ic ) + { + lid = wellElemDofNumber[iwelem] + ic + 1 - rankOffset; + real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic] + + scalingFactor * localSolution[lid]; + + if( newDens < 0 ) + { + minVal.min( 0 ); + } + } + } + else { - lid = wellElemDofNumber[iwelem] + ic + 1 - rankOffset; - real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic] - + scalingFactor * localSolution[lid]; - if( newDens < 0.0 ) + real64 totalDens = 0.0; + for( localIndex ic = 0; ic < numComponents; ++ic ) + { + lid = wellElemDofNumber[iwelem] + ic + 1 - rankOffset; + real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic] + + scalingFactor * localSolution[lid]; + totalDens += (newDens > 0.0) ? newDens : 0.0; + } + if( totalDens < eps ) { minVal.min( 0 ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/co2flash.txt b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/co2flash.txt new file mode 100644 index 00000000000..712dfc99728 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/co2flash.txt @@ -0,0 +1 @@ +FlashModel CO2Solubility 1e6 1.5e7 5e4 94 96 1 0 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 0ecc5360c69..d249164da3a 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 @@ -51,8 +51,8 @@ type="injector" control="liquidRate" targetBHP="1e8" - targetRate="1e-4" - injectionStream="{ 0.1, 0.1, 0.8 }"/> + targetRate="1e-1" + injectionStream="{ 0.0, 0.0, 1.0 }"/> @@ -134,21 +134,21 @@ + maxTime="1e6"> diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtgas.txt b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtgas.txt new file mode 100644 index 00000000000..fd3d781edae --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtgas.txt @@ -0,0 +1,2 @@ +DensityFun SpanWagnerCO2Density 1e6 1.5e7 5e4 94 96 1 +ViscosityFun FenghourCO2Viscosity 1e6 1.5e7 5e4 94 96 1 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtliquid.txt b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtliquid.txt new file mode 100644 index 00000000000..f5852c95fd9 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtliquid.txt @@ -0,0 +1,2 @@ +DensityFun BrineCO2Density 1e6 1.5e7 5e4 94 96 1 0 +ViscosityFun BrineViscosity 0 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/staircase_compositional_multiphase_wells_3d.xml b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/staircase_co2_wells_3d.xml similarity index 78% rename from src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/staircase_compositional_multiphase_wells_3d.xml rename to src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/staircase_co2_wells_3d.xml index bcc2afc820a..a2014343fdb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/staircase_compositional_multiphase_wells_3d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/staircase_co2_wells_3d.xml @@ -25,8 +25,9 @@ fluidNames="{ fluid1 }" solidNames="{ rock }" relPermNames="{ relperm }" - temperature="297.15" - useMass="0"/> + temperature="368.15" + maxCompFractionChange="0.2" + useMass="1"/> + wellTemperature="368.15" + maxCompFractionChange="0.2" + useMass="1"> + injectionStream="{ 0.995, 0.005 }"/> @@ -145,20 +147,13 @@ - + phaseNames="{ gas, water }" + componentNames="{ co2, water }" + componentMolarWeight="{ 44e-3, 18e-3 }" + phasePVTParaFiles="{ pvtgas.txt, pvtliquid.txt }" + flashModelParaFile="co2flash.txt"/> @@ -215,43 +210,25 @@ setNames="{ all }" objectPath="ElementRegions/Channel" fieldName="pressure" - scale="1e6"/> + scale="9e6"/> + scale="0.005"/> - - - - + scale="0.995"/> diff --git a/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp b/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp index a0e0410b271..4ca4213b9e4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp @@ -391,5 +391,18 @@ void ReservoirSolverBase::ImplicitStepComplete( real64 const & time_n, void ReservoirSolverBase::ResetViews( DomainPartition * const GEOSX_UNUSED_PARAM( domain ) ) {} +real64 ReservoirSolverBase::ScalingForSystemSolution( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) +{ + real64 const flowScalingFactor = m_flowSolver->ScalingForSystemSolution( domain, dofManager, localSolution ); + real64 const wellScalingFactor = m_wellSolver->ScalingForSystemSolution( domain, dofManager, localSolution ); + + GEOSX_LOG_LEVEL_RANK_0( 2, "Scaling factor for the reservoir: " << flowScalingFactor + << "; for the well(s): " << wellScalingFactor ); + + return LvArray::math::min( flowScalingFactor, wellScalingFactor ); +} + } /* namespace geosx */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.hpp index 7b6ff6604ba..ff962948f11 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.hpp @@ -118,6 +118,11 @@ class ReservoirSolverBase : public SolverBase ParallelVector & rhs, ParallelVector & solution ) override; + virtual real64 + ScalingForSystemSolution( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) override; + virtual bool CheckSystemSolution( DomainPartition const & domain, DofManager const & dofManager, diff --git a/src/externalComponents/PVTPackage b/src/externalComponents/PVTPackage index 0c4bb093e81..ac3e7d18495 160000 --- a/src/externalComponents/PVTPackage +++ b/src/externalComponents/PVTPackage @@ -1 +1 @@ -Subproject commit 0c4bb093e81e3e64163392d9f08b1c81361ba5f7 +Subproject commit ac3e7d184955ac8a2a55c6183f85e072d1950b8f