From 9a5804e2caf3085d9270e79766784cda67a536cc Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Mon, 20 Jul 2020 18:38:45 -0700 Subject: [PATCH 1/7] implemented global damping based on compDens relative change --- .../fluid/MultiPhaseMultiComponentFluid.hpp | 2 +- .../docs/CompositionalMultiphaseFlow.rst | 41 +++--- .../docs/CompositionalMultiphaseWell.rst | 35 ++--- src/coreComponents/fileIO/schema/schema.xsd | 4 + .../fluidFlow/CompositionalMultiphaseFlow.cpp | 129 +++++++++++++++++- .../fluidFlow/CompositionalMultiphaseFlow.hpp | 30 +++- .../fluidFlow/CompositionalMultiphaseWell.cpp | 113 ++++++++++++++- .../fluidFlow/CompositionalMultiphaseWell.hpp | 27 +++- .../CompositionalMultiphaseWellKernels.hpp | 55 +++++++- ...case_compositional_multiphase_wells_3d.xml | 2 + .../multiphysics/ReservoirSolverBase.cpp | 13 ++ .../multiphysics/ReservoirSolverBase.hpp | 5 + 12 files changed, 402 insertions(+), 54 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/MultiPhaseMultiComponentFluid.hpp b/src/coreComponents/constitutive/fluid/MultiPhaseMultiComponentFluid.hpp index 42e32f80e31..82bf8097a95 100644 --- a/src/coreComponents/constitutive/fluid/MultiPhaseMultiComponentFluid.hpp +++ b/src/coreComponents/constitutive/fluid/MultiPhaseMultiComponentFluid.hpp @@ -166,7 +166,7 @@ class MultiPhaseMultiComponentFluidUpdate final : public MultiFluidBaseUpdate m_phaseDensity[k][q], m_dPhaseDensity_dPressure[k][q], m_dPhaseDensity_dTemperature[k][q], - m_dPhaseFraction_dGlobalCompFraction[k][q], + m_dPhaseDensity_dGlobalCompFraction[k][q], m_phaseViscosity[k][q], m_dPhaseViscosity_dPressure[k][q], m_dPhaseViscosity_dTemperature[k][q], diff --git a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst index 320d28be1f7..225966dc1c3 100644 --- a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst +++ b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst @@ -1,24 +1,25 @@ -========================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== -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 +============================ ============ ======== ====================================================================================================================================================================================================================================================================================================================== +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 +maxRelComponentDensityChange real64 1 Maximum (relative) change in a component density 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..f0fac2ba354 100644 --- a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseWell.rst +++ b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseWell.rst @@ -1,21 +1,22 @@ -========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== -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 +============================ ============ ======== ======================================================================================================================================================================================================================================================================================================================== +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 +maxRelComponentDensityChange real64 1 Maximum (relative) change in a component density 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 2b8bb4bd3b1..3c56fbc776d 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd +++ b/src/coreComponents/fileIO/schema/schema.xsd @@ -830,6 +830,8 @@ + + @@ -883,6 +885,8 @@ + + diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp index a886fa3ba4f..35403f4e72a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp @@ -47,7 +47,10 @@ CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, FlowSolverBase( name, parent ), m_numPhases( 0 ), m_numComponents( 0 ), - m_capPressureFlag( 0 ) + m_capPressureFlag( 0 ), + m_maxRelCompDensChange( 1.0 ), + m_compDensCheckTol( 1e-10 ), + m_minScalingFactor( 1e-2 ) { //START_SPHINX_INCLUDE_00 this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> @@ -69,6 +72,12 @@ CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, setInputFlag( InputFlags::OPTIONAL )-> setDescription( "Name of the capillary pressure constitutive model to use" ); + this->registerWrapper( viewKeyStruct::maxRelCompDensChangeString, &m_maxRelCompDensChange )-> + setSizedFromParent( 0 )-> + setInputFlag( InputFlags::OPTIONAL )-> + setApplyDefaultValue( 1.0 )-> + setDescription( "Maximum (relative) change in a component density between two Newton iterations" ); + m_linearSolverParameters.get().mgr.strategy = "CompositionalMultiphaseFlow"; } @@ -78,6 +87,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_maxRelCompDensChange, 1.0, + "The maximum relative change in component density must smaller or equal to 1.0" ); + GEOSX_ERROR_IF_LT_MSG( m_maxRelCompDensChange, 0.0, + "The maximum relative change in component density must larger or equal to 0.0" ); } void CompositionalMultiphaseFlow::RegisterDataOnMesh( Group * const MeshBodies ) @@ -1199,13 +1213,81 @@ 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_maxRelCompDensChange >= 1.0 ) + { + // no rescaling wanted, we just return 1.0; + return 1.0; + } + + real64 constexpr eps = 1e-10; + real64 const maxRelCompDensChange = m_maxRelCompDensChange; + + 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 ) + { + + // 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 prevCompDens = compDens[ei][ic] + dCompDens[ei][ic]; + real64 const absCompDensChange = fabs( localSolution[lid] ); + real64 const maxAbsCompDensChange = maxRelCompDensChange * prevCompDens; + // TODO: see if this is robust for the transition prevCompDens \approx 0 => compDens > 0 + if( absCompDensChange > maxAbsCompDensChange && maxAbsCompDensChange > eps ) + { + minVal.min( maxAbsCompDensChange / absCompDensChange ); + } + } + } + } ); + + if( minVal.get() < scalingFactor ) + { + scalingFactor = minVal.get(); + } + } ); + + return LvArray::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 ) { localIndex const NC = m_numComponents; - + real64 const checkTol = m_compDensCheckTol; + MeshLevel const & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); globalIndex const rankOffset = dofManager.rankOffset(); @@ -1222,9 +1304,9 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d 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, integer > check( 1 ); + RAJA::ReduceMin< serialReduce, integer > check( 1 ); - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) + forAll< serialPolicy >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { if( elemGhostRank[ei] < 0 ) { @@ -1236,7 +1318,7 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d 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 ); + check.min( newDens >= - checkTol ); } } } ); @@ -1252,8 +1334,7 @@ void CompositionalMultiphaseFlow::ApplySystemSolution( DofManager const & dofMan real64 const scalingFactor, DomainPartition & domain ) { - MeshLevel & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); - + MeshLevel & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); dofManager.addVectorToField( localSolution, viewKeyStruct::dofFieldString, viewKeyStruct::deltaPressureString, @@ -1266,6 +1347,8 @@ void CompositionalMultiphaseFlow::ApplySystemSolution( DofManager const & dofMan scalingFactor, 1, m_numDofPerCell ); + ChopNegativeDensities( domain ); + std::map< string, string_array > fieldNames; fieldNames["elems"].emplace_back( string( viewKeyStruct::deltaPressureString ) ); fieldNames["elems"].emplace_back( string( viewKeyStruct::deltaGlobalCompDensityString ) ); @@ -1277,6 +1360,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..80e805ad1c8 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, @@ -310,6 +315,10 @@ class CompositionalMultiphaseFlow : public FlowSolverBase static constexpr auto phaseViscosityString = "phaseViscosity"; static constexpr auto phaseRelativePermeabilityString = "phaseRelativePermeability"; static constexpr auto phaseCapillaryPressureString = "phaseCapillaryPressure"; + + // accepted relative change in component density between two Newton iterations + static constexpr auto maxRelCompDensChangeString = "maxRelComponentDensityChange"; + } viewKeysCompMultiphaseFlow; struct groupKeyStruct : SolverBase::groupKeyStruct @@ -362,7 +371,7 @@ class CompositionalMultiphaseFlow : public FlowSolverBase DomainPartition & domain, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const; - + protected: virtual void PostProcessInput() override; @@ -393,7 +402,12 @@ class CompositionalMultiphaseFlow : public FlowSolverBase */ void ResetViews( MeshLevel & mesh ) override; - + /** + * @brief sets all the negative component densities (if any) to zero + * @param domain the physical domain object + */ + void ChopNegativeDensities( DomainPartition & domain ); + /// the max number of fluid phases localIndex m_numPhases; @@ -415,6 +429,18 @@ 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_maxAbsCompFracChange; + + /// maximum relative change in a component density between two Newton iterations + real64 m_maxRelCompDensChange; + + /// tolerance used in CheckSystemSolution to check acceptable values of component densities + real64 const m_compDensCheckTol; + + /// minimum value of the scaling factor obtained by enforcing maxRelCompDensChange and maxAbsCompFracChange + real64 const m_minScalingFactor; + ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_pressure; ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_deltaPressure; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWell.cpp index 6b527c0f9d2..0c00f3bc9f1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/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_maxRelCompDensChange( 1.0 ), + m_compDensCheckTol( 1e-10 ), + m_minScalingFactor( 1e-2 ) { this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> setInputFlag( InputFlags::REQUIRED )-> @@ -65,6 +68,13 @@ 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::maxRelCompDensChangeString, &m_maxRelCompDensChange )-> + setSizedFromParent( 0 )-> + setInputFlag( InputFlags::OPTIONAL )-> + setApplyDefaultValue( 1.0 )-> + setDescription( "Maximum (relative) change in a component density between two Newton iterations" ); + } void CompositionalMultiphaseWell::PostProcessInput() @@ -76,6 +86,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_maxRelCompDensChange, 1.0, + "The maximum relative change in component density must smaller or equal to 1.0" ); + GEOSX_ERROR_IF_LT_MSG( m_maxRelCompDensChange, 0.0, + "The maximum relative change in component density must larger or equal to 0.0" ); + } void CompositionalMultiphaseWell::RegisterDataOnMesh( Group * const meshBodies ) @@ -627,6 +643,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_maxRelCompDensChange >= 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_maxRelCompDensChange ); + + + if (subRegionScalingFactor < scalingFactor) + { + scalingFactor = subRegionScalingFactor; + } + } ); + + return LvArray::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); +} + bool CompositionalMultiphaseWell::CheckSystemSolution( DomainPartition const & domain, DofManager const & dofManager, @@ -670,6 +740,7 @@ CompositionalMultiphaseWell::CheckSystemSolution( DomainPartition const & domain dWellElemPressure, wellElemCompDens, dWellElemCompDens, + m_compDensCheckTol, scalingFactor ); if( subRegionSolutionCheck == 0 ) @@ -776,6 +847,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 +866,9 @@ CompositionalMultiphaseWell::ApplySystemSolution( DofManager const & dofManager, scalingFactor, m_numDofPerWellElement - 1, m_numDofPerWellElement ); + 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 +883,42 @@ 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]; + // 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/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWell.hpp index fe6a38afc51..a30231c0aba 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/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, @@ -209,7 +214,7 @@ class CompositionalMultiphaseWell : public WellSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) override; - + arrayView1d< string const > const & relPermModelNames() const { return m_relPermModelNames; } struct viewKeyStruct : WellSolverBase::viewKeyStruct @@ -250,6 +255,10 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr auto compPerforationRateString = "compPerforationRate"; static constexpr auto dCompPerforationRate_dPresString = "dCompPerforationRate_dPres"; static constexpr auto dCompPerforationRate_dCompString = "dCompPerforationRate_dComp"; + + // accepted relative and absolute changes between two Newton iterations + static constexpr auto maxRelCompDensChangeString = "maxRelComponentDensityChange"; + } viewKeysCompMultiphaseWell; struct groupKeyStruct : SolverBase::groupKeyStruct @@ -306,11 +315,11 @@ 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 + * @brief sets all the negative component densities (if any) to zero + * @param domain the physical domain object */ - void RecordWellData( WellElementSubRegion const & subRegion ); - + void ChopNegativeDensities( DomainPartition & domain ); + /// the max number of fluid phases localIndex m_numPhases; @@ -326,7 +335,15 @@ class CompositionalMultiphaseWell : public WellSolverBase /// list of relative permeability model names per target region array1d< string > m_relPermModelNames; + /// maximum relative change in a component density between two Newton iterations + real64 m_maxRelCompDensChange; + + /// tolerance used in CheckSystemSolution to check acceptable values of component densities + real64 const m_compDensCheckTol; + /// minimum value of the scaling factor obtained by enforcing maxRelCompDensChange + real64 const m_minScalingFactor; + /// views into reservoir primary variable fields ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_resPressure; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWellKernels.hpp index fb278d5be37..4363b0f2d11 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseWellKernels.hpp @@ -1209,6 +1209,53 @@ 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 maxRelCompDensChange ) + { + real64 constexpr eps = 1e-10; + + RAJA::ReduceMin< REDUCE_POLICY, real64 > minVal( 1.0 ); + + forAll< POLICY >( wellElemDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) + { + if( wellElemGhostRank[iwelem] < 0 ) + { + + 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 prevCompDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic]; + real64 const absCompDensChange = fabs( localSolution[lid] ); + real64 const maxAbsCompDensChange = maxRelCompDensChange * prevCompDens; + + // if the relative change is too large, chop back + if( absCompDensChange > maxAbsCompDensChange && maxAbsCompDensChange > eps ) + { + minVal.min( maxAbsCompDensChange / absCompDensChange ); + } + } + } + } ); + return minVal.get(); + } + +}; + + /******************************** SolutionCheckKernel ********************************/ struct SolutionCheckKernel @@ -1224,6 +1271,7 @@ struct SolutionCheckKernel arrayView1d< real64 const > const & dWellElemPressure, arrayView2d< real64 const > const & wellElemCompDens, arrayView2d< real64 const > const & dWellElemCompDens, + real64 const densCheckTol, real64 const scalingFactor ) { RAJA::ReduceMin< REDUCE_POLICY, localIndex > minVal( 1 ); @@ -1236,6 +1284,8 @@ 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 ); @@ -1247,7 +1297,10 @@ struct SolutionCheckKernel lid = wellElemDofNumber[iwelem] + ic + 1 - rankOffset; real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic] + scalingFactor * localSolution[lid]; - if( newDens < 0.0 ) + + // the component density must be positive, up to a tolerance + // (slightly) negative densities are chopped back in CompositionalMultiphaseWell::ApplySystemSolution + if( newDens < -densCheckTol ) { minVal.min( 0 ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseWell/staircase_compositional_multiphase_wells_3d.xml b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseWell/staircase_compositional_multiphase_wells_3d.xml index bcc2afc820a..289724688f4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseWell/staircase_compositional_multiphase_wells_3d.xml +++ b/src/coreComponents/physicsSolvers/fluidFlow/integratedTests/compositionalMultiphaseWell/staircase_compositional_multiphase_wells_3d.xml @@ -26,6 +26,7 @@ solidNames="{ rock }" relPermNames="{ relperm }" temperature="297.15" + maxRelComponentDensityChange="0.9" useMass="0"/> 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::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, From 0c603297650db257255bce75af25a1bffb54ce36 Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Mon, 20 Jul 2020 20:10:37 -0700 Subject: [PATCH 2/7] uncrustified --- .../fluidFlow/CompositionalMultiphaseFlow.cpp | 49 ++++++++++--------- .../fluidFlow/CompositionalMultiphaseFlow.hpp | 25 +++++----- .../wells/CompositionalMultiphaseWell.cpp | 27 +++++----- .../wells/CompositionalMultiphaseWell.hpp | 26 +++++----- .../CompositionalMultiphaseWellKernels.hpp | 17 ++++--- .../multiphysics/ReservoirSolverBase.cpp | 2 +- 6 files changed, 75 insertions(+), 71 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp index 35403f4e72a..d0b5e94185d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp @@ -50,7 +50,7 @@ CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, m_capPressureFlag( 0 ), m_maxRelCompDensChange( 1.0 ), m_compDensCheckTol( 1e-10 ), - m_minScalingFactor( 1e-2 ) + m_minScalingFactor( 5e-2 ) { //START_SPHINX_INCLUDE_00 this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> @@ -77,7 +77,7 @@ CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, setInputFlag( InputFlags::OPTIONAL )-> setApplyDefaultValue( 1.0 )-> setDescription( "Maximum (relative) change in a component density between two Newton iterations" ); - + m_linearSolverParameters.get().mgr.strategy = "CompositionalMultiphaseFlow"; } @@ -89,9 +89,9 @@ void CompositionalMultiphaseFlow::PostProcessInput() m_capPressureFlag = CheckModelNames( m_capPressureModelNames, viewKeyStruct::capPressureNamesString, true ); GEOSX_ERROR_IF_GT_MSG( m_maxRelCompDensChange, 1.0, - "The maximum relative change in component density must smaller or equal to 1.0" ); + "The maximum relative change in component density must smaller or equal to 1.0" ); GEOSX_ERROR_IF_LT_MSG( m_maxRelCompDensChange, 0.0, - "The maximum relative change in component density must larger or equal to 0.0" ); + "The maximum relative change in component density must larger or equal to 0.0" ); } void CompositionalMultiphaseFlow::RegisterDataOnMesh( Group * const MeshBodies ) @@ -1218,19 +1218,19 @@ real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition co arrayView1d< real64 const > const & localSolution ) { GEOSX_MARK_FUNCTION; - + // check if we want to rescale the Newton update if( m_maxRelCompDensChange >= 1.0 ) { // no rescaling wanted, we just return 1.0; return 1.0; - } + } + + real64 constexpr eps = 1e-10; + real64 const maxRelCompDensChange = m_maxRelCompDensChange; - real64 constexpr eps = 1e-10; - real64 const maxRelCompDensChange = m_maxRelCompDensChange; - localIndex const NC = m_numComponents; - + MeshLevel const & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); globalIndex const rankOffset = dofManager.rankOffset(); @@ -1246,22 +1246,23 @@ real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition co 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 ) { - // compute the change in component densities and component fractions + // 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 prevCompDens = compDens[ei][ic] + dCompDens[ei][ic]; + real64 const prevCompDens = compDens[ei][ic] + dCompDens[ei][ic]; real64 const absCompDensChange = fabs( localSolution[lid] ); - real64 const maxAbsCompDensChange = maxRelCompDensChange * prevCompDens; - // TODO: see if this is robust for the transition prevCompDens \approx 0 => compDens > 0 + real64 const maxAbsCompDensChange = maxRelCompDensChange * prevCompDens; + // TODO: see if this is robust for the transition prevCompDens \approx 0 => compDens > 0 + // TODO: this may produce a very conservative chopping strategy, revisit later if too restrictive if( absCompDensChange > maxAbsCompDensChange && maxAbsCompDensChange > eps ) { minVal.min( maxAbsCompDensChange / absCompDensChange ); @@ -1276,7 +1277,7 @@ real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition co } } ); - return LvArray::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); + return LvArray::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); } @@ -1287,7 +1288,7 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d { localIndex const NC = m_numComponents; real64 const checkTol = m_compDensCheckTol; - + MeshLevel const & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); globalIndex const rankOffset = dofManager.rankOffset(); @@ -1304,9 +1305,9 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d arrayView2d< real64 const > const & compDens = subRegion.getReference< array2d< real64 > >( viewKeyStruct::globalCompDensityString ); arrayView2d< real64 const > const & dCompDens = subRegion.getReference< array2d< real64 > >( viewKeyStruct::deltaGlobalCompDensityString ); - RAJA::ReduceMin< serialReduce, integer > check( 1 ); + RAJA::ReduceMin< parallelDeviceReduce, integer > check( 1 ); - forAll< serialPolicy >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { if( elemGhostRank[ei] < 0 ) { @@ -1318,7 +1319,7 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d for( localIndex ic = 0; ic < NC; ++ic ) { real64 const newDens = compDens[ei][ic] + dCompDens[ei][ic] + scalingFactor * localSolution[localRow + ic + 1]; - check.min( newDens >= - checkTol ); + check.min( newDens >= -checkTol ); } } } ); @@ -1334,7 +1335,7 @@ void CompositionalMultiphaseFlow::ApplySystemSolution( DofManager const & dofMan real64 const scalingFactor, DomainPartition & domain ) { - MeshLevel & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); + MeshLevel & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); dofManager.addVectorToField( localSolution, viewKeyStruct::dofFieldString, viewKeyStruct::deltaPressureString, @@ -1347,6 +1348,8 @@ void CompositionalMultiphaseFlow::ApplySystemSolution( DofManager const & dofMan scalingFactor, 1, m_numDofPerCell ); + // some densities may be slightly negative, they are chopped here + // TODO: find a way to do that in CheckSystemSolution or ScalingForSystemSolution ChopNegativeDensities( domain ); std::map< string, string_array > fieldNames; @@ -1363,7 +1366,7 @@ 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 ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp index 80e805ad1c8..604a997aa81 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp @@ -317,8 +317,8 @@ class CompositionalMultiphaseFlow : public FlowSolverBase static constexpr auto phaseCapillaryPressureString = "phaseCapillaryPressure"; // accepted relative change in component density between two Newton iterations - static constexpr auto maxRelCompDensChangeString = "maxRelComponentDensityChange"; - + static constexpr auto maxRelCompDensChangeString = "maxRelComponentDensityChange"; + } viewKeysCompMultiphaseFlow; struct groupKeyStruct : SolverBase::groupKeyStruct @@ -371,7 +371,13 @@ class CompositionalMultiphaseFlow : public FlowSolverBase DomainPartition & domain, 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; @@ -402,12 +408,6 @@ class CompositionalMultiphaseFlow : public FlowSolverBase */ void ResetViews( MeshLevel & mesh ) override; - /** - * @brief sets all the negative component densities (if any) to zero - * @param domain the physical domain object - */ - void ChopNegativeDensities( DomainPartition & domain ); - /// the max number of fluid phases localIndex m_numPhases; @@ -429,18 +429,15 @@ 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_maxAbsCompFracChange; - /// maximum relative change in a component density between two Newton iterations real64 m_maxRelCompDensChange; - /// tolerance used in CheckSystemSolution to check acceptable values of component densities + /// tolerance used in CheckSystemSolution to check acceptable values of component densities real64 const m_compDensCheckTol; /// minimum value of the scaling factor obtained by enforcing maxRelCompDensChange and maxAbsCompFracChange real64 const m_minScalingFactor; - + ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_pressure; ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_deltaPressure; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 9d0cda93cb0..8b40849d2f2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -54,7 +54,7 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, m_useMass( false ), m_maxRelCompDensChange( 1.0 ), m_compDensCheckTol( 1e-10 ), - m_minScalingFactor( 1e-2 ) + m_minScalingFactor( 5e-2 ) { this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> setInputFlag( InputFlags::REQUIRED )-> @@ -74,7 +74,7 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, setInputFlag( InputFlags::OPTIONAL )-> setApplyDefaultValue( 1.0 )-> setDescription( "Maximum (relative) change in a component density between two Newton iterations" ); - + } void CompositionalMultiphaseWell::PostProcessInput() @@ -88,10 +88,10 @@ void CompositionalMultiphaseWell::PostProcessInput() "(referenced from well solver " << getName() << ")" ); GEOSX_ERROR_IF_GT_MSG( m_maxRelCompDensChange, 1.0, - "The maximum relative change in component density must smaller or equal to 1.0" ); + "The maximum relative change in component density must smaller or equal to 1.0" ); GEOSX_ERROR_IF_LT_MSG( m_maxRelCompDensChange, 0.0, - "The maximum relative change in component density must larger or equal to 0.0" ); - + "The maximum relative change in component density must larger or equal to 0.0" ); + } void CompositionalMultiphaseWell::RegisterDataOnMesh( Group * const meshBodies ) @@ -655,8 +655,8 @@ CompositionalMultiphaseWell::ScalingForSystemSolution( DomainPartition const & d { // no rescaling wanted, we just return 1.0; return 1.0; - } - + } + MeshLevel const & meshLevel = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); real64 scalingFactor = 1.0; @@ -686,15 +686,15 @@ CompositionalMultiphaseWell::ScalingForSystemSolution( DomainPartition const & d wellElemCompDens, dWellElemCompDens, m_maxRelCompDensChange ); - - if (subRegionScalingFactor < scalingFactor) + + if( subRegionScalingFactor < scalingFactor ) { scalingFactor = subRegionScalingFactor; } } ); - return LvArray::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); + return LvArray::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); } bool @@ -866,6 +866,8 @@ CompositionalMultiphaseWell::ApplySystemSolution( DofManager const & dofManager, scalingFactor, m_numDofPerWellElement - 1, m_numDofPerWellElement ); + // some densities may be slightly negative, they are chopped here + // TODO: find a way to do that in CheckSystemSolution or ScalingForSystemSolution ChopNegativeDensities( domain ); // synchronize @@ -905,9 +907,10 @@ void CompositionalMultiphaseWell::ChopNegativeDensities( DomainPartition & domai { for( localIndex ic = 0; ic < NC; ++ic ) { - // get the latest component density (i.e., after the update of dWellElemCompDens) + // get the latest component density (i.e., after the update of dWellElemCompDens) real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic]; - // if the new density is negative, chop back to zero + // 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]; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index f152b4d2a0e..0e0e8e147d6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -214,7 +214,13 @@ class CompositionalMultiphaseWell : public WellSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, 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 @@ -257,8 +263,8 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr auto dCompPerforationRate_dCompString = "dCompPerforationRate_dComp"; // accepted relative and absolute changes between two Newton iterations - static constexpr auto maxRelCompDensChangeString = "maxRelComponentDensityChange"; - + static constexpr auto maxRelCompDensChangeString = "maxRelComponentDensityChange"; + } viewKeysCompMultiphaseWell; struct groupKeyStruct : SolverBase::groupKeyStruct @@ -314,12 +320,6 @@ class CompositionalMultiphaseWell : public WellSolverBase */ void ResizeFields( WellElementSubRegion & subRegion ); - /** - * @brief sets all the negative component densities (if any) to zero - * @param domain the physical domain object - */ - void ChopNegativeDensities( DomainPartition & domain ); - /// the max number of fluid phases localIndex m_numPhases; @@ -337,13 +337,13 @@ class CompositionalMultiphaseWell : public WellSolverBase /// maximum relative change in a component density between two Newton iterations real64 m_maxRelCompDensChange; - - /// tolerance used in CheckSystemSolution to check acceptable values of component densities - real64 const m_compDensCheckTol; + + /// tolerance used in CheckSystemSolution to check acceptable values of component densities + real64 const m_compDensCheckTol; /// minimum value of the scaling factor obtained by enforcing maxRelCompDensChange real64 const m_minScalingFactor; - + /// views into reservoir primary variable fields ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_resPressure; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index 686f32fabd3..fcd19ce63a9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -1224,8 +1224,8 @@ struct SolutionScalingKernel arrayView2d< real64 const > const & dWellElemCompDens, real64 const maxRelCompDensChange ) { - real64 constexpr eps = 1e-10; - + real64 constexpr eps = 1e-10; + RAJA::ReduceMin< REDUCE_POLICY, real64 > minVal( 1.0 ); forAll< POLICY >( wellElemDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) @@ -1236,13 +1236,14 @@ struct SolutionScalingKernel 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 prevCompDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic]; real64 const absCompDensChange = fabs( localSolution[lid] ); real64 const maxAbsCompDensChange = maxRelCompDensChange * prevCompDens; - + // if the relative change is too large, chop back + // TODO: this may produce a very conservative chopping strategy, revisit later if too restrictive if( absCompDensChange > maxAbsCompDensChange && maxAbsCompDensChange > eps ) { minVal.min( maxAbsCompDensChange / absCompDensChange ); @@ -1284,8 +1285,8 @@ 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 + + // the pressure must be positive if( newPres < 0.0 ) { minVal.min( 0 ); @@ -1298,8 +1299,8 @@ struct SolutionCheckKernel real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic] + scalingFactor * localSolution[lid]; - // the component density must be positive, up to a tolerance - // (slightly) negative densities are chopped back in CompositionalMultiphaseWell::ApplySystemSolution + // the component density must be positive, up to a tolerance + // (slightly) negative densities are chopped back in CompositionalMultiphaseWell::ApplySystemSolution if( newDens < -densCheckTol ) { minVal.min( 0 ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp b/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp index 7a8a3c11bdb..5ddfdb04576 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp @@ -401,7 +401,7 @@ real64 ReservoirSolverBase::ScalingForSystemSolution( DomainPartition const & do GEOSX_LOG_LEVEL_RANK_0( 2, "Scaling factor for the reservoir: " << flowScalingFactor << "; for the well(s): " << wellScalingFactor ); - return LvArray::min( flowScalingFactor, wellScalingFactor ); + return LvArray::min( flowScalingFactor, wellScalingFactor ); } From c48208608779d1712984dce857f3e47ccacd0fba Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Thu, 23 Jul 2020 23:04:24 -0700 Subject: [PATCH 3/7] addressed Sergey and Matteo's comments --- .../docs/CompositionalMultiphaseFlow.rst | 43 ++++++------ .../docs/CompositionalMultiphaseWell.rst | 37 +++++----- src/coreComponents/fileIO/schema/schema.xsd | 12 ++-- .../fluidFlow/CompositionalMultiphaseFlow.cpp | 70 +++++++++++++------ .../fluidFlow/CompositionalMultiphaseFlow.hpp | 19 ++--- .../wells/CompositionalMultiphaseWell.cpp | 39 +++++++---- .../wells/CompositionalMultiphaseWell.hpp | 27 +++---- .../CompositionalMultiphaseWellKernels.hpp | 47 ++++++++----- 8 files changed, 170 insertions(+), 124 deletions(-) diff --git a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst index 225966dc1c3..375fc56129a 100644 --- a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst +++ b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseFlow.rst @@ -1,25 +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 -maxRelComponentDensityChange real64 1 Maximum (relative) change in a component density 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` -============================ ============ ======== ====================================================================================================================================================================================================================================================================================================================== +============================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== +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 f0fac2ba354..a0dca92ba59 100644 --- a/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseWell.rst +++ b/src/coreComponents/fileIO/schema/docs/CompositionalMultiphaseWell.rst @@ -1,22 +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 -maxRelComponentDensityChange real64 1 Maximum (relative) change in a component density 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` -============================ ============ ======== ======================================================================================================================================================================================================================================================================================================================== +============================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== +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 3c56fbc776d..39a67e3004a 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd +++ b/src/coreComponents/fileIO/schema/schema.xsd @@ -816,6 +816,8 @@ + + @@ -830,8 +832,8 @@ - - + + @@ -875,6 +877,8 @@ + + @@ -885,8 +889,8 @@ - - + + diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp index d0b5e94185d..2a41707f299 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp @@ -48,9 +48,9 @@ CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, m_numPhases( 0 ), m_numComponents( 0 ), m_capPressureFlag( 0 ), - m_maxRelCompDensChange( 1.0 ), - m_compDensCheckTol( 1e-10 ), - m_minScalingFactor( 5e-2 ) + m_maxCompFracChange( 1.0 ), + m_minScalingFactor( 0.1 ), + m_allowCompDensChopping( 1 ) { //START_SPHINX_INCLUDE_00 this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> @@ -72,11 +72,17 @@ CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, setInputFlag( InputFlags::OPTIONAL )-> setDescription( "Name of the capillary pressure constitutive model to use" ); - this->registerWrapper( viewKeyStruct::maxRelCompDensChangeString, &m_maxRelCompDensChange )-> + this->registerWrapper( viewKeyStruct::maxCompFracChangeString, &m_maxCompFracChange )-> setSizedFromParent( 0 )-> setInputFlag( InputFlags::OPTIONAL )-> setApplyDefaultValue( 1.0 )-> - setDescription( "Maximum (relative) change in a component density between two Newton iterations" ); + 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"; @@ -88,10 +94,10 @@ void CompositionalMultiphaseFlow::PostProcessInput() CheckModelNames( m_relPermModelNames, viewKeyStruct::relPermNamesString ); m_capPressureFlag = CheckModelNames( m_capPressureModelNames, viewKeyStruct::capPressureNamesString, true ); - GEOSX_ERROR_IF_GT_MSG( m_maxRelCompDensChange, 1.0, - "The maximum relative change in component density must smaller or equal to 1.0" ); - GEOSX_ERROR_IF_LT_MSG( m_maxRelCompDensChange, 0.0, - "The maximum relative change in component density must larger or equal to 0.0" ); + 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 ) @@ -1220,14 +1226,14 @@ real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition co GEOSX_MARK_FUNCTION; // check if we want to rescale the Newton update - if( m_maxRelCompDensChange >= 1.0 ) + if( m_maxCompFracChange >= 1.0 ) { // no rescaling wanted, we just return 1.0; return 1.0; } real64 constexpr eps = 1e-10; - real64 const maxRelCompDensChange = m_maxRelCompDensChange; + real64 const maxCompFracChange = m_maxCompFracChange; localIndex const NC = m_numComponents; @@ -1251,6 +1257,11 @@ real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition co { 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 ) @@ -1258,12 +1269,16 @@ real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition co localIndex const lid = dofNumber[ei] + ic + 1 - rankOffset; // compute scaling factor based on relative change in component densities - real64 const prevCompDens = compDens[ei][ic] + dCompDens[ei][ic]; real64 const absCompDensChange = fabs( localSolution[lid] ); - real64 const maxAbsCompDensChange = maxRelCompDensChange * prevCompDens; - // TODO: see if this is robust for the transition prevCompDens \approx 0 => compDens > 0 - // TODO: this may produce a very conservative chopping strategy, revisit later if too restrictive - if( absCompDensChange > maxAbsCompDensChange && maxAbsCompDensChange > eps ) + 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 ); } @@ -1287,7 +1302,7 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d real64 const scalingFactor ) { localIndex const NC = m_numComponents; - real64 const checkTol = m_compDensCheckTol; + integer const allowCompDensChopping = m_allowCompDensChopping; MeshLevel const & mesh = *domain.getMeshBody( 0 )->getMeshLevel( 0 ); @@ -1316,10 +1331,16 @@ 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 is not allowed, the time step fails if a component density is negative + // otherwise, 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 >= -checkTol ); + 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 ); + } } } } ); @@ -1348,9 +1369,12 @@ void CompositionalMultiphaseFlow::ApplySystemSolution( DofManager const & dofMan scalingFactor, 1, m_numDofPerCell ); - // some densities may be slightly negative, they are chopped here - // TODO: find a way to do that in CheckSystemSolution or ScalingForSystemSolution - ChopNegativeDensities( domain ); + // 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 ) ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp index 604a997aa81..9dcfbce873a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.hpp @@ -285,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"; @@ -316,9 +319,6 @@ class CompositionalMultiphaseFlow : public FlowSolverBase static constexpr auto phaseRelativePermeabilityString = "phaseRelativePermeability"; static constexpr auto phaseCapillaryPressureString = "phaseCapillaryPressure"; - // accepted relative change in component density between two Newton iterations - static constexpr auto maxRelCompDensChangeString = "maxRelComponentDensityChange"; - } viewKeysCompMultiphaseFlow; struct groupKeyStruct : SolverBase::groupKeyStruct @@ -429,14 +429,15 @@ class CompositionalMultiphaseFlow : public FlowSolverBase /// name of the cap pressure constitutive model array1d< string > m_capPressureModelNames; - /// maximum relative change in a component density between two Newton iterations - real64 m_maxRelCompDensChange; + /// 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; - /// tolerance used in CheckSystemSolution to check acceptable values of component densities - real64 const m_compDensCheckTol; + /// flag indicating whether local (cell-wise) chopping of negative compositions is allowed + integer m_allowCompDensChopping; - /// minimum value of the scaling factor obtained by enforcing maxRelCompDensChange and maxAbsCompFracChange - real64 const m_minScalingFactor; ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_pressure; ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > m_deltaPressure; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 8b40849d2f2..f4dc22c17f0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -52,9 +52,9 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, m_numComponents( 0 ), m_temperature( 0.0 ), m_useMass( false ), - m_maxRelCompDensChange( 1.0 ), - m_compDensCheckTol( 1e-10 ), - m_minScalingFactor( 5e-2 ) + m_maxCompFracChange( 1.0 ), + m_minScalingFactor( 0.1 ), + m_allowCompDensChopping( 1 ) { this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> setInputFlag( InputFlags::REQUIRED )-> @@ -69,11 +69,17 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, setInputFlag( InputFlags::REQUIRED )-> setDescription( "Names of relative permeability constitutive models to use" ); - this->registerWrapper( viewKeyStruct::maxRelCompDensChangeString, &m_maxRelCompDensChange )-> + this->registerWrapper( viewKeyStruct::maxCompFracChangeString, &m_maxCompFracChange )-> setSizedFromParent( 0 )-> setInputFlag( InputFlags::OPTIONAL )-> setApplyDefaultValue( 1.0 )-> - setDescription( "Maximum (relative) change in a component density between two Newton iterations" ); + 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" ); } @@ -87,10 +93,10 @@ void CompositionalMultiphaseWell::PostProcessInput() "Flow solver " << GetFlowSolverName() << " not found or incompatible type " "(referenced from well solver " << getName() << ")" ); - GEOSX_ERROR_IF_GT_MSG( m_maxRelCompDensChange, 1.0, - "The maximum relative change in component density must smaller or equal to 1.0" ); - GEOSX_ERROR_IF_LT_MSG( m_maxRelCompDensChange, 0.0, - "The maximum relative change in component density must larger or equal to 0.0" ); + 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" ); } @@ -651,7 +657,7 @@ CompositionalMultiphaseWell::ScalingForSystemSolution( DomainPartition const & d GEOSX_MARK_FUNCTION; // check if we want to rescale the Newton update - if( m_maxRelCompDensChange >= 1.0 ) + if( m_maxCompFracChange >= 1.0 ) { // no rescaling wanted, we just return 1.0; return 1.0; @@ -685,7 +691,7 @@ CompositionalMultiphaseWell::ScalingForSystemSolution( DomainPartition const & d wellElemGhostRank, wellElemCompDens, dWellElemCompDens, - m_maxRelCompDensChange ); + m_maxCompFracChange ); if( subRegionScalingFactor < scalingFactor ) @@ -740,7 +746,7 @@ CompositionalMultiphaseWell::CheckSystemSolution( DomainPartition const & domain dWellElemPressure, wellElemCompDens, dWellElemCompDens, - m_compDensCheckTol, + m_allowCompDensChopping, scalingFactor ); if( subRegionSolutionCheck == 0 ) @@ -866,9 +872,12 @@ CompositionalMultiphaseWell::ApplySystemSolution( DofManager const & dofManager, scalingFactor, m_numDofPerWellElement - 1, m_numDofPerWellElement ); - // some densities may be slightly negative, they are chopped here - // TODO: find a way to do that in CheckSystemSolution or ScalingForSystemSolution - ChopNegativeDensities( domain ); + // 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; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index 0e0e8e147d6..c8f5cc675cb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -229,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; @@ -247,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 = @@ -262,9 +260,6 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr auto dCompPerforationRate_dPresString = "dCompPerforationRate_dPres"; static constexpr auto dCompPerforationRate_dCompString = "dCompPerforationRate_dComp"; - // accepted relative and absolute changes between two Newton iterations - static constexpr auto maxRelCompDensChangeString = "maxRelComponentDensityChange"; - } viewKeysCompMultiphaseWell; struct groupKeyStruct : SolverBase::groupKeyStruct @@ -335,14 +330,14 @@ class CompositionalMultiphaseWell : public WellSolverBase /// list of relative permeability model names per target region array1d< string > m_relPermModelNames; - /// maximum relative change in a component density between two Newton iterations - real64 m_maxRelCompDensChange; + /// maximum (absolute) change in a component fraction between two Newton iterations + real64 m_maxCompFracChange; - /// tolerance used in CheckSystemSolution to check acceptable values of component densities - real64 const m_compDensCheckTol; + /// minimum value of the scaling factor obtained by enforcing maxCompFracChange + real64 m_minScalingFactor; - /// minimum value of the scaling factor obtained by enforcing maxRelCompDensChange - real64 const 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 fcd19ce63a9..86baef48bad 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -1222,7 +1222,7 @@ struct SolutionScalingKernel arrayView1d< integer const > const & wellElemGhostRank, arrayView2d< real64 const > const & wellElemCompDens, arrayView2d< real64 const > const & dWellElemCompDens, - real64 const maxRelCompDensChange ) + real64 const maxCompFracChange ) { real64 constexpr eps = 1e-10; @@ -1233,18 +1233,27 @@ struct SolutionScalingKernel 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 prevCompDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic]; real64 const absCompDensChange = fabs( localSolution[lid] ); - real64 const maxAbsCompDensChange = maxRelCompDensChange * prevCompDens; - - // if the relative change is too large, chop back - // TODO: this may produce a very conservative chopping strategy, revisit later if too restrictive - if( absCompDensChange > maxAbsCompDensChange && maxAbsCompDensChange > eps ) + 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 ); } @@ -1272,7 +1281,7 @@ struct SolutionCheckKernel arrayView1d< real64 const > const & dWellElemPressure, arrayView2d< real64 const > const & wellElemCompDens, arrayView2d< real64 const > const & dWellElemCompDens, - real64 const densCheckTol, + integer const allowCompDensChopping, real64 const scalingFactor ) { RAJA::ReduceMin< REDUCE_POLICY, localIndex > minVal( 1 ); @@ -1292,18 +1301,20 @@ struct SolutionCheckKernel 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, negative component densities will be chopped (i.e., set to zero in ApplySystemSolution) + if( !allowCompDensChopping ) { - lid = wellElemDofNumber[iwelem] + ic + 1 - rankOffset; - real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic] - + scalingFactor * localSolution[lid]; - - // the component density must be positive, up to a tolerance - // (slightly) negative densities are chopped back in CompositionalMultiphaseWell::ApplySystemSolution - if( newDens < -densCheckTol ) + for( localIndex ic = 0; ic < numComponents; ++ic ) { - minVal.min( 0 ); + lid = wellElemDofNumber[iwelem] + ic + 1 - rankOffset; + real64 const newDens = wellElemCompDens[iwelem][ic] + dWellElemCompDens[iwelem][ic] + + scalingFactor * localSolution[lid]; + + if( newDens < 0 ) + { + minVal.min( 0 ); + } } } } From a6731e134e721ca1643e584c02185cee6e5e4675 Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Thu, 23 Jul 2020 23:25:10 -0700 Subject: [PATCH 4/7] changed a well integrated test to use CO2; fixed PVT bug --- .../SPE10/dead_oil_spe10_layers_83_84_85.xml | 34 ++++----- .../deadoil_3ph_baker_1d.xml | 14 ++-- .../deadoil_3ph_corey_1d.xml | 14 ++-- .../deadoil_3ph_staircase_3d.xml | 20 +++--- .../compositionalMultiphaseWell/co2flash.txt | 1 + .../dead_oil_wells_2d.xml | 12 ++-- .../compositionalMultiphaseWell/pvtgas.txt | 2 + .../compositionalMultiphaseWell/pvtliquid.txt | 2 + ...ells_3d.xml => staircase_co2_wells_3d.xml} | 69 ++++++------------- src/externalComponents/PVTPackage | 2 +- 10 files changed, 72 insertions(+), 98 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/co2flash.txt create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtgas.txt create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/pvtliquid.txt rename src/coreComponents/physicsSolvers/fluidFlow/wells/integratedTests/compositionalMultiphaseWell/{staircase_compositional_multiphase_wells_3d.xml => staircase_co2_wells_3d.xml} (77%) 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 45fa7ac919d..7dac476cd02 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"> @@ -262,7 +263,7 @@ @@ -275,13 +276,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/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 77% 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 289724688f4..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,9 +25,9 @@ fluidNames="{ fluid1 }" solidNames="{ rock }" relPermNames="{ relperm }" - temperature="297.15" - maxRelComponentDensityChange="0.9" - useMass="0"/> + temperature="368.15" + maxCompFractionChange="0.2" + useMass="1"/> + wellTemperature="368.15" + maxCompFractionChange="0.2" + useMass="1"> + injectionStream="{ 0.995, 0.005 }"/> @@ -147,20 +147,13 @@ - + phaseNames="{ gas, water }" + componentNames="{ co2, water }" + componentMolarWeight="{ 44e-3, 18e-3 }" + phasePVTParaFiles="{ pvtgas.txt, pvtliquid.txt }" + flashModelParaFile="co2flash.txt"/> @@ -217,43 +210,25 @@ setNames="{ all }" objectPath="ElementRegions/Channel" fieldName="pressure" - scale="1e6"/> + scale="9e6"/> + scale="0.005"/> - - - - + scale="0.995"/> diff --git a/src/externalComponents/PVTPackage b/src/externalComponents/PVTPackage index 0c4bb093e81..8d9ded9ee1b 160000 --- a/src/externalComponents/PVTPackage +++ b/src/externalComponents/PVTPackage @@ -1 +1 @@ -Subproject commit 0c4bb093e81e3e64163392d9f08b1c81361ba5f7 +Subproject commit 8d9ded9ee1b1e112ff5533e4bd1b392740c9b2bd From 1f3450c2eb8dbba63ba9f9542a4b52edb79eae54 Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Mon, 3 Aug 2020 10:10:23 -0700 Subject: [PATCH 5/7] fixed bug and division by zero in properties --- .../PVTFunctions/CO2SolubilityFunction.cpp | 12 ++++++++++-- .../BrooksCoreyRelativePermeability.hpp | 2 +- .../fluidFlow/CompositionalMultiphaseFlow.cpp | 17 ++++++++++++++--- .../wells/CompositionalMultiphaseWell.cpp | 2 +- .../CompositionalMultiphaseWellKernels.hpp | 18 +++++++++++++++++- 5 files changed, 43 insertions(+), 8 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp b/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp index 566909fb065..358178a063a 100644 --- a/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp +++ b/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp @@ -335,7 +335,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 +345,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 - 1e-10 ) + { + Y = compFraction[m_CO2Index] * 1e10; + } + else + { + Y = compFraction[m_CO2Index] / (1.0 - compFraction[m_CO2Index]); + } if( Y < X ) { @@ -356,7 +362,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/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp index 2a41707f299..e06361f3cb5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp @@ -49,7 +49,7 @@ CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, m_numComponents( 0 ), m_capPressureFlag( 0 ), m_maxCompFracChange( 1.0 ), - m_minScalingFactor( 0.1 ), + m_minScalingFactor( 0.01 ), m_allowCompDensChopping( 1 ) { //START_SPHINX_INCLUDE_00 @@ -1332,8 +1332,9 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d check.min( newPres >= 0.0 ); } - // if component density is not allowed, the time step fails if a component density is negative - // otherwise, negative component densities will be chopped (i.e., set to zero in ApplySystemSolution) + // 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 ) { for( localIndex ic = 0; ic < NC; ++ic ) @@ -1342,6 +1343,16 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d 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 >= 1e-6 ); + } } } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index f4dc22c17f0..c77e681acc1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -53,7 +53,7 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, m_temperature( 0.0 ), m_useMass( false ), m_maxCompFracChange( 1.0 ), - m_minScalingFactor( 0.1 ), + m_minScalingFactor( 0.01 ), m_allowCompDensChopping( 1 ) { this->registerWrapper( viewKeyStruct::temperatureString, &m_temperature )-> diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index 86baef48bad..2a8bce62dee 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -1302,7 +1302,8 @@ struct SolutionCheckKernel } // if component density is not allowed, the time step fails if a component density is negative - // otherwise, negative component densities will be chopped (i.e., set to zero in ApplySystemSolution) + // 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 ) @@ -1317,6 +1318,21 @@ struct SolutionCheckKernel } } } + else + { + 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 < 1e-6 ) + { + minVal.min( 0 ); + } + } } } ); return minVal.get(); From 0606cf9184dcee26e216cebebebd36434380ff93 Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Mon, 3 Aug 2020 10:57:54 -0700 Subject: [PATCH 6/7] fixed post-merge issues --- integratedTests | 2 +- .../physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp | 2 +- .../fluidFlow/wells/CompositionalMultiphaseWell.cpp | 2 +- .../physicsSolvers/multiphysics/ReservoirSolverBase.cpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/integratedTests b/integratedTests index 574d8a7cc7e..c9cb5ddfdad 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit 574d8a7cc7ef0e310df61e4337392f5ab659dc0b +Subproject commit c9cb5ddfdad3520ff56c5b59e16ba63303ab530e diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp index 56231eaa89b..d51ea168e37 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp @@ -1292,7 +1292,7 @@ real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition co } } ); - return LvArray::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); + return LvArray::math::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index c77e681acc1..c336c8ea4b1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -700,7 +700,7 @@ CompositionalMultiphaseWell::ScalingForSystemSolution( DomainPartition const & d } } ); - return LvArray::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); + return LvArray::math::max( MpiWrapper::Min( scalingFactor, MPI_COMM_GEOSX ), m_minScalingFactor ); } bool diff --git a/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp b/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp index 5ddfdb04576..4ca4213b9e4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/ReservoirSolverBase.cpp @@ -401,7 +401,7 @@ real64 ReservoirSolverBase::ScalingForSystemSolution( DomainPartition const & do GEOSX_LOG_LEVEL_RANK_0( 2, "Scaling factor for the reservoir: " << flowScalingFactor << "; for the well(s): " << wellScalingFactor ); - return LvArray::min( flowScalingFactor, wellScalingFactor ); + return LvArray::math::min( flowScalingFactor, wellScalingFactor ); } From ab49ee41ef68274df462df321958fff705918a42 Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Thu, 6 Aug 2020 15:57:11 -0700 Subject: [PATCH 7/7] consolidated the constants --- .../fluid/PVTFunctions/CO2SolubilityFunction.cpp | 12 ++++-------- .../fluidFlow/CompositionalMultiphaseFlow.cpp | 8 ++++++-- .../wells/CompositionalMultiphaseWellKernels.hpp | 8 ++++++-- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp b/src/coreComponents/constitutive/fluid/PVTFunctions/CO2SolubilityFunction.cpp index 358178a063a..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[] = @@ -345,9 +341,9 @@ void CO2SolubilityFunction::Partition( EvalVarArgs const & pressure, EvalVarArgs //Y = C/W = z/(1-z) - if( compFraction[m_CO2Index].m_var > 1.0 - 1e-10 ) + if( compFraction[m_CO2Index].m_var > 1.0 - minForDivision ) { - Y = compFraction[m_CO2Index] * 1e10; + Y = compFraction[m_CO2Index] / minForDivision; } else { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp index d51ea168e37..5dd2c703dab 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFlow.cpp @@ -41,6 +41,8 @@ using namespace dataRepository; using namespace constitutive; using namespace CompositionalMultiphaseFlowKernels; +static constexpr real64 minDensForDivision = 1e-10; + CompositionalMultiphaseFlow::CompositionalMultiphaseFlow( const string & name, Group * const parent ) : @@ -1232,7 +1234,7 @@ real64 CompositionalMultiphaseFlow::ScalingForSystemSolution( DomainPartition co return 1.0; } - real64 constexpr eps = 1e-10; + real64 constexpr eps = minDensForDivision; real64 const maxCompFracChange = m_maxCompFracChange; localIndex const NC = m_numComponents; @@ -1301,6 +1303,8 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d arrayView1d< real64 const > const & localSolution, real64 const scalingFactor ) { + real64 constexpr eps = minDensForDivision; + localIndex const NC = m_numComponents; integer const allowCompDensChopping = m_allowCompDensChopping; @@ -1351,7 +1355,7 @@ bool CompositionalMultiphaseFlow::CheckSystemSolution( DomainPartition const & d real64 const newDens = compDens[ei][ic] + dCompDens[ei][ic] + scalingFactor * localSolution[localRow + ic + 1]; totalDens += (newDens > 0.0) ? newDens : 0.0; } - check.min( totalDens >= 1e-6 ); + check.min( totalDens >= eps ); } } } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index 2a8bce62dee..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 @@ -1224,7 +1226,7 @@ struct SolutionScalingKernel arrayView2d< real64 const > const & dWellElemCompDens, real64 const maxCompFracChange ) { - real64 constexpr eps = 1e-10; + real64 constexpr eps = minDensForDivision; RAJA::ReduceMin< REDUCE_POLICY, real64 > minVal( 1.0 ); @@ -1284,6 +1286,8 @@ struct SolutionCheckKernel 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 ) @@ -1328,7 +1332,7 @@ struct SolutionCheckKernel + scalingFactor * localSolution[lid]; totalDens += (newDens > 0.0) ? newDens : 0.0; } - if( totalDens < 1e-6 ) + if( totalDens < eps ) { minVal.min( 0 ); }