diff --git a/src/coreComponents/codingUtilities/Utilities.hpp b/src/coreComponents/codingUtilities/Utilities.hpp index 868d096c28b..c50a2bfc05d 100644 --- a/src/coreComponents/codingUtilities/Utilities.hpp +++ b/src/coreComponents/codingUtilities/Utilities.hpp @@ -59,6 +59,15 @@ bool isZero( T const val, T const tol = LvArray::NumericLimits< T >::epsilon ) return -tol <= val && val <= tol; } +template< typename ARRAY_TYPE > +GEOS_FORCE_INLINE GEOS_HOST_DEVICE +bool hasNonZero( ARRAY_TYPE const & array ) +{ + return std::any_of( array.begin(), array.end(), []( real64 value ) { + return !isZero( value ); // Check if the value is non-zero + } ); +} + template< typename T > GEOS_FORCE_INLINE GEOS_HOST_DEVICE constexpr bool isOdd( T x ) diff --git a/src/coreComponents/functions/TableFunction.cpp b/src/coreComponents/functions/TableFunction.cpp index 3fc62833d41..a26baaa6654 100644 --- a/src/coreComponents/functions/TableFunction.cpp +++ b/src/coreComponents/functions/TableFunction.cpp @@ -196,6 +196,11 @@ real64 TableFunction::evaluate( real64 const * const input ) const return m_kernelWrapper.compute( input ); } +real64 TableFunction::getCoord( real64 const * const input, localIndex const dim, InterpolationType interpolationMethod ) const +{ + return m_kernelWrapper.getCoord( input, dim, interpolationMethod ); +} + TableFunction::KernelWrapper::KernelWrapper( InterpolationType const interpolationMethod, ArrayOfArraysView< real64 const > const & coordinates, arrayView1d< real64 const > const & values ) diff --git a/src/coreComponents/functions/TableFunction.hpp b/src/coreComponents/functions/TableFunction.hpp index b7be867d123..e759abb3fc1 100644 --- a/src/coreComponents/functions/TableFunction.hpp +++ b/src/coreComponents/functions/TableFunction.hpp @@ -173,6 +173,18 @@ class TableFunction : public FunctionBase real64 interpolateRound( IN_ARRAY const & input ) const; + /** + * @brief ... + * @param[in] input vector of input value + * @param[in] dim table dimension + * @param[in] interpolationMethod interpolation method + * @return coordinate value + */ + template< typename IN_ARRAY > + GEOS_HOST_DEVICE + real64 + getCoord( IN_ARRAY const & input, localIndex dim, InterpolationType interpolationMethod ) const; + /** * @brief Interpolate in the table with derivatives using linear method. * @param[in] input vector of input value @@ -240,6 +252,15 @@ class TableFunction : public FunctionBase */ virtual real64 evaluate( real64 const * const input ) const override final; + /** + * @brief Method to get coordinates + * @param input a scalar input + * @param dim the table dimension + * @param interpolationMethod the interpolation method + * @return the coordinate + */ + real64 getCoord( real64 const * const input, localIndex dim, InterpolationType interpolationMethod ) const; + /** * @brief Check if the given coordinate is in the bounds of the table coordinates in the * specified dimension, throw an exception otherwise. @@ -554,6 +575,57 @@ TableFunction::KernelWrapper::interpolateRound( IN_ARRAY const & input ) const return m_values[tableIndex]; } +template< typename IN_ARRAY > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 +TableFunction::KernelWrapper::getCoord( IN_ARRAY const & input, localIndex const dim, InterpolationType interpolationMethod ) const +{ + // Determine the index to the nearest table entry + localIndex subIndex; + arraySlice1d< real64 const > const coords = m_coordinates[dim]; + // Determine the index along each table axis + if( input[dim] <= coords[0] ) + { + // Coordinate is to the left of the table axis + subIndex = 0; + } + else if( input[dim] >= coords[coords.size() - 1] ) + { + // Coordinate is to the right of the table axis + subIndex = coords.size() - 1; + } + else + { + // Coordinate is within the table axis + // Note: find() will return the index of the upper table vertex + auto const lower = LvArray::sortedArrayManipulation::find( coords.begin(), coords.size(), input[dim] ); + subIndex = LvArray::integerConversion< localIndex >( lower ); + + // Interpolation types: + // - Nearest returns the value of the closest table vertex + // - Upper returns the value of the next table vertex + // - Lower returns the value of the previous table vertex + if( interpolationMethod == TableFunction::InterpolationType::Nearest ) + { + if( ( input[dim] - coords[subIndex - 1]) <= ( coords[subIndex] - input[dim]) ) + { + --subIndex; + } + } + else if( interpolationMethod == TableFunction::InterpolationType::Lower ) + { + if( subIndex > 0 ) + { + --subIndex; + } + } + } + + // Retrieve the nearest coordinate + return coords[subIndex]; +} + template< typename IN_ARRAY, typename OUT_ARRAY > GEOS_HOST_DEVICE GEOS_FORCE_INLINE diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index e0eefbde8c2..7a3bb36484f 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -303,7 +303,7 @@ bool PhysicsSolverBase::execute( real64 const time_n, if( dtRemaining > 0.0 ) { - nextDt = setNextDt( dtAccepted, domain ); + nextDt = setNextDt( time_n + (dt - dtRemaining), dtAccepted, domain ); if( nextDt < dtRemaining ) { // better to do two equal steps than one big and one small (even potentially tiny) @@ -333,7 +333,7 @@ bool PhysicsSolverBase::execute( real64 const time_n, " reached. Consider increasing maxSubSteps." ); // Decide what to do with the next Dt for the event running the solver. - m_nextDt = setNextDt( nextDt, domain ); + m_nextDt = setNextDt( time_n + dt, nextDt, domain ); // Increase counter to indicate how many cycles since the last timestep cut if( m_numTimestepsSinceLastDtCut >= 0 ) @@ -364,15 +364,16 @@ void PhysicsSolverBase::logEndOfCycleInformation( integer const cycleNumber, GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, "------------------------------------------------------------------\n" ); } -real64 PhysicsSolverBase::setNextDt( real64 const & currentDt, +real64 PhysicsSolverBase::setNextDt( real64 const & GEOS_UNUSED_PARAM( currentTime ), + real64 const & currentDt, DomainPartition & domain ) { integer const minTimeStepIncreaseInterval = m_nonlinearSolverParameters.minTimeStepIncreaseInterval(); - real64 const nextDtNewton = setNextDtBasedOnNewtonIter( currentDt ); - if( m_nonlinearSolverParameters.getLogLevel() > 0 ) - GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on Newton iterations = {}", getName(), nextDtNewton )); + real64 const nextDtNewton = setNextDtBasedOnIterNumber( currentDt ); + if( m_nonlinearSolverParameters.getLogLevel() > 0 && nextDtNewton < LvArray::NumericLimits< real64 >::max ) + GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on number of iterations = {}", getName(), nextDtNewton )); real64 const nextDtStateChange = setNextDtBasedOnStateChange( currentDt, domain ); - if( m_nonlinearSolverParameters.getLogLevel() > 0 ) + if( m_nonlinearSolverParameters.getLogLevel() > 0 && nextDtStateChange < LvArray::NumericLimits< real64 >::max ) GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on state change = {}", getName(), nextDtStateChange )); if( ( m_numTimestepsSinceLastDtCut >= 0 ) && ( m_numTimestepsSinceLastDtCut < minTimeStepIncreaseInterval ) ) @@ -429,7 +430,7 @@ real64 PhysicsSolverBase::setNextDtBasedOnStateChange( real64 const & currentDt, return LvArray::NumericLimits< real64 >::max; // i.e., not implemented } -real64 PhysicsSolverBase::setNextDtBasedOnNewtonIter( real64 const & currentDt ) +real64 PhysicsSolverBase::setNextDtBasedOnIterNumber( real64 const & currentDt ) { integer & newtonIter = m_nonlinearSolverParameters.m_numNewtonIterations; integer const iterDecreaseLimit = m_nonlinearSolverParameters.timeStepDecreaseIterLimit(); diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index 95a87f1a87f..a7f24c23314 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -216,19 +216,21 @@ class PhysicsSolverBase : public ExecutableGroup /** * @brief function to set the next time step size + * @param[in] currentTime the current time * @param[in] currentDt the current time step size * @param[in] domain the domain object * @return the prescribed time step size */ - virtual real64 setNextDt( real64 const & currentDt, + virtual real64 setNextDt( real64 const & currentTime, + real64 const & currentDt, DomainPartition & domain ); /** - * @brief function to set the next time step size based on Newton convergence + * @brief function to set the next time step size based on convergence * @param[in] currentDt the current time step size * @return the prescribed time step size */ - virtual real64 setNextDtBasedOnNewtonIter( real64 const & currentDt ); + virtual real64 setNextDtBasedOnIterNumber( real64 const & currentDt ); /** * @brief function to set the next dt based on state change diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 67ee8bf2f9a..61a88e68230 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -1340,10 +1340,12 @@ void CompositionalMultiphaseFVM::assembleHydrofracFluxTerms( real64 const GEOS_U } ); } -real64 CompositionalMultiphaseFVM::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain ) +real64 CompositionalMultiphaseFVM::setNextDt( real64 const & currentTime, + real64 const & currentDt, + DomainPartition & domain ) { if( m_targetFlowCFL < 0 ) - return CompositionalMultiphaseBase::setNextDt( currentDt, domain ); + return CompositionalMultiphaseBase::setNextDt( currentTime, currentDt, domain ); else return setNextDtBasedOnCFL( currentDt, domain ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp index 14d66667035..955abb0e74e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp @@ -190,11 +190,13 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase /** * @brief function to set the next time step size + * @param[in] currentTime the current time * @param[in] currentDt the current time step size * @param[in] domain the domain object * @return the prescribed time step size */ - real64 setNextDt( real64 const & currentDt, + real64 setNextDt( real64 const & currentTime, + real64 const & currentDt, DomainPartition & domain ) override; struct viewKeyStruct : CompositionalMultiphaseBase::viewKeyStruct diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index e33e4ec0b17..241898b888a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -274,7 +274,7 @@ void CompositionalMultiphaseWell::registerDataOnMesh( Group & meshBodies ) integer const numPhase = m_numPhases; // format: time,bhp,total_rate,total_vol_rate,phase0_vol_rate,phase1_vol_rate,... std::ofstream outputFile( m_ratesOutputDir + "/" + wellControlsName + ".csv" ); - outputFile << "Time [s],BHP [Pa],Total rate [" << massUnit << "/s],Total " << conditionKey << " Volumetric rate [" << unitKey << "m3/s]"; + outputFile << "Time [s],Time step [s],BHP [Pa],Total rate [" << massUnit << "/s],Total " << conditionKey << " Volumetric rate [" << unitKey << "m3/s]"; for( integer ip = 0; ip < numPhase; ++ip ) outputFile << ",Phase" << ip << " " << conditionKey << " volumetric rate [" << unitKey << "m3/s]"; outputFile << std::endl; @@ -443,7 +443,7 @@ void CompositionalMultiphaseWell::validateInjectionStreams( WellElementSubRegion } void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n, - real64 const & dt, + real64 const & GEOS_UNUSED_PARAM( dt ), WellElementSubRegion const & subRegion ) { string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString()); @@ -452,9 +452,9 @@ void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n // now that we know we are single-phase, we can check a few things in the constraints WellControls const & wellControls = getWellControls( subRegion ); WellControls::Control const currentControl = wellControls.getControl(); - real64 const & targetTotalRate = wellControls.getTargetTotalRate( time_n + dt ); - real64 const & targetPhaseRate = wellControls.getTargetPhaseRate( time_n + dt ); - real64 const & targetMassRate = wellControls.getTargetMassRate( time_n + dt ); + real64 const & targetTotalRate = wellControls.getTargetTotalRate( time_n ); + real64 const & targetPhaseRate = wellControls.getTargetPhaseRate( time_n ); + real64 const & targetMassRate = wellControls.getTargetMassRate( time_n ); GEOS_THROW_IF( wellControls.isInjector() && currentControl == WellControls::Control::PHASEVOLRATE, "WellControls " << wellControls.getDataContext() << @@ -978,7 +978,7 @@ real64 CompositionalMultiphaseWell::updateSubRegionState( WellElementSubRegion & return maxPhaseVolChange; } -void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, real64 const & time_n, real64 const & dt ) +void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, real64 const & time_n ) { GEOS_MARK_FUNCTION; @@ -1005,13 +1005,11 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea WellElementSubRegion & subRegion ) { WellControls const & wellControls = getWellControls( subRegion ); + PerforationData const & perforationData = *subRegion.getPerforationData(); + arrayView2d< real64 const > const compPerfRate = perforationData.getField< fields::well::compPerforationRate >(); - if( time_n <= 0.0 || - ( !wellControls.isWellOpen( time_n ) && wellControls.isWellOpen( time_n + dt ) ) ) + if( time_n <= 0.0 || ( wellControls.isWellOpen( time_n ) && !hasNonZero( compPerfRate ) ) ) { - - PerforationData const & perforationData = *subRegion.getPerforationData(); - // get well primary variables on well elements arrayView1d< real64 > const & wellElemPressure = subRegion.getField< fields::well::pressure >(); arrayView1d< real64 > const & wellElemTemp = subRegion.getField< fields::well::temperature >(); @@ -1092,7 +1090,7 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea launch( subRegion.size(), m_targetPhaseIndex, wellControls, - 0.0, // initialization done at t = 0 + time_n, // initialization done at time_n wellElemPhaseDens, wellElemTotalDens, connRate ); @@ -1125,7 +1123,7 @@ void CompositionalMultiphaseWell::assembleFluxTerms( real64 const & time, WellElementSubRegion & subRegion ) { WellControls const & well_controls = getWellControls( subRegion ); - if( well_controls.isWellOpen( time + dt ) && !m_keepVariablesConstantDuringInitStep ) + if( well_controls.isWellOpen( time ) && !m_keepVariablesConstantDuringInitStep ) { string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString()); MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); @@ -1195,7 +1193,7 @@ void CompositionalMultiphaseWell::assembleAccumulationTerms( real64 const & time int numPhases = fluid.numFluidPhases(); int numComponents = fluid.numFluidComponents(); WellControls const & wellControls = getWellControls( subRegion ); - if( wellControls.isWellOpen( time+ dt ) && !m_keepVariablesConstantDuringInitStep ) + if( wellControls.isWellOpen( time ) && !m_keepVariablesConstantDuringInitStep ) { if( isThermal() ) { @@ -1323,7 +1321,7 @@ CompositionalMultiphaseWell::calculateResidualNorm( real64 const & time_n, subRegion, fluid, wellControls, - time_n + dt, + time_n, dt, m_nonlinearSolverParameters.m_minNormalizer, subRegionResidualNorm ); @@ -1351,7 +1349,7 @@ CompositionalMultiphaseWell::calculateResidualNorm( real64 const & time_n, subRegion, fluid, wellControls, - time_n + dt, + time_n, dt, m_nonlinearSolverParameters.m_minNormalizer, subRegionResidualNorm ); @@ -1657,7 +1655,9 @@ CompositionalMultiphaseWell::checkSystemSolution( DomainPartition & domain, return MpiWrapper::min( localCheck ); } -void CompositionalMultiphaseWell::computePerforationRates( real64 const & time_n, real64 const & dt, DomainPartition & domain ) +void CompositionalMultiphaseWell::computePerforationRates( real64 const & time_n, + real64 const & GEOS_UNUSED_PARAM( dt ), + DomainPartition & domain ) { GEOS_MARK_FUNCTION; @@ -1675,7 +1675,7 @@ void CompositionalMultiphaseWell::computePerforationRates( real64 const & time_n { PerforationData * const perforationData = subRegion.getPerforationData(); WellControls const & wellControls = getWellControls( subRegion ); - if( wellControls.isWellOpen( time_n + dt ) && !m_keepVariablesConstantDuringInitStep ) + if( wellControls.isWellOpen( time_n ) && !m_keepVariablesConstantDuringInitStep ) { bool const disableReservoirToWellFlow = wellControls.isInjector() and !wellControls.isCrossflowEnabled(); @@ -1897,7 +1897,7 @@ void CompositionalMultiphaseWell::resetStateToBeginningOfStep( DomainPartition & } void CompositionalMultiphaseWell::assemblePressureRelations( real64 const & time_n, - real64 const & dt, + real64 const & GEOS_UNUSED_PARAM( dt ), DomainPartition const & domain, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, @@ -1919,7 +1919,7 @@ void CompositionalMultiphaseWell::assemblePressureRelations( real64 const & time WellControls & wellControls = getWellControls( subRegion ); - if( wellControls.isWellOpen( time_n + dt ) && !m_keepVariablesConstantDuringInitStep ) + if( wellControls.isWellOpen( time_n ) && !m_keepVariablesConstantDuringInitStep ) { string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); @@ -1954,7 +1954,7 @@ void CompositionalMultiphaseWell::assemblePressureRelations( real64 const & time subRegion.getTopWellElementIndex(), m_targetPhaseIndex, wellControls, - time_n + dt, // controls evaluated with BHP/rate of the end of step + time_n, // controls evaluated with BHP/rate of the beginning of step wellElemDofNumber, wellElemGravCoef, nextWellElemIndex, @@ -1970,32 +1970,30 @@ void CompositionalMultiphaseWell::assemblePressureRelations( real64 const & time // TODO: move the switch logic into wellControls // TODO: implement a more general switch when more then two constraints per well type are allowed - real64 const timeAtEndOfStep = time_n + dt; - if( wellControls.getControl() == WellControls::Control::BHP ) { if( wellControls.isProducer() ) { - wellControls.switchToPhaseRateControl( wellControls.getTargetPhaseRate( timeAtEndOfStep ) ); + wellControls.switchToPhaseRateControl( wellControls.getTargetPhaseRate( time_n ) ); GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, GEOS_FMT( "Control switch for well {} from BHP constraint to phase volumetric rate constraint", subRegion.getName() ) ); } else if( wellControls.getInputControl() == WellControls::Control::MASSRATE ) { - wellControls.switchToMassRateControl( wellControls.getTargetMassRate( timeAtEndOfStep ) ); + wellControls.switchToMassRateControl( wellControls.getTargetMassRate( time_n ) ); GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, GEOS_FMT( "Control switch for well {} from BHP constraint to mass rate constraint", subRegion.getName()) ); } else { - wellControls.switchToTotalRateControl( wellControls.getTargetTotalRate( timeAtEndOfStep ) ); + wellControls.switchToTotalRateControl( wellControls.getTargetTotalRate( time_n ) ); GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, GEOS_FMT( "Control switch for well {} from BHP constraint to total volumetric rate constraint", subRegion.getName()) ); } } else { - wellControls.switchToBHPControl( wellControls.getTargetBHP( timeAtEndOfStep ) ); + wellControls.switchToBHPControl( wellControls.getTargetBHP( time_n ) ); GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, GEOS_FMT( "Control switch for well {} from rate constraint to BHP constraint", subRegion.getName() ) ); } @@ -2117,10 +2115,10 @@ void CompositionalMultiphaseWell::printRates( real64 const & time_n, if( m_writeCSV > 0 ) { outputFile.open( m_ratesOutputDir + "/" + wellControlsName + ".csv", std::ios_base::app ); - outputFile << time_n + dt; + outputFile << time_n << "," << dt; } - if( !wellControls.isWellOpen( time_n + dt ) ) + if( !wellControls.isWellOpen( time_n ) ) { GEOS_LOG( GEOS_FMT( "{}: well is shut", wellControlsName ) ); if( outputFile.is_open()) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index 72614596376..bead7cb12f7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -357,7 +357,7 @@ class CompositionalMultiphaseWell : public WellSolverBase * @brief Initialize all the primary and secondary variables in all the wells * @param domain the domain containing the well manager to access individual wells */ - void initializeWells( DomainPartition & domain, real64 const & time_n, real64 const & dt ) override; + void initializeWells( DomainPartition & domain, real64 const & time_n ) override; virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 75797a3490c..87b81ec8c24 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -121,13 +121,13 @@ string SinglePhaseWell::resElementDofName() const } void SinglePhaseWell::validateWellConstraints( real64 const & time_n, - real64 const & dt, + real64 const & GEOS_UNUSED_PARAM( dt ), WellElementSubRegion const & subRegion ) { WellControls const & wellControls = getWellControls( subRegion ); WellControls::Control const currentControl = wellControls.getControl(); - real64 const targetTotalRate = wellControls.getTargetTotalRate( time_n + dt ); - real64 const targetPhaseRate = wellControls.getTargetPhaseRate( time_n + dt ); + real64 const targetTotalRate = wellControls.getTargetTotalRate( time_n ); + real64 const targetPhaseRate = wellControls.getTargetPhaseRate( time_n ); GEOS_THROW_IF( currentControl == WellControls::Control::PHASEVOLRATE, "WellControls " << wellControls.getDataContext() << ": Phase rate control is not available for SinglePhaseWell", @@ -334,18 +334,11 @@ real64 SinglePhaseWell::updateSubRegionState( WellElementSubRegion & subRegion ) return 0.0; // change in phasevolume fraction doesnt apply } -void SinglePhaseWell::initializeWells( DomainPartition & domain, real64 const & time_n, real64 const & dt ) +void SinglePhaseWell::initializeWells( DomainPartition & domain, real64 const & time_n ) { GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( time_n ); - GEOS_UNUSED_VAR( dt ); - // different functionality than for compositional - // compositional has better treatment of well initialization in cases where well starts later in sim - // logic will be incorporated into single phase in different push - if( time_n > 0 ) - { - return; - } + // loop over the wells forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & meshLevel, @@ -381,49 +374,153 @@ void SinglePhaseWell::initializeWells( DomainPartition & domain, real64 const & arrayView1d< real64 const > const & perfGravCoef = perforationData.getField< fields::well::gravityCoefficient >(); - // TODO: change the way we access the flowSolver here - SinglePhaseBase const & flowSolver = getParent().getGroup< SinglePhaseBase >( getFlowSolverName() ); - PresInitializationKernel::SinglePhaseFlowAccessors resSinglePhaseFlowAccessors( meshLevel.getElemManager(), flowSolver.getName() ); - PresInitializationKernel::SingleFluidAccessors resSingleFluidAccessors( meshLevel.getElemManager(), flowSolver.getName() ); - - // 1) Loop over all perforations to compute an average density - // 2) Initialize the reference pressure - // 3) Estimate the pressures in the well elements using the average density - PresInitializationKernel:: - launch( perforationData.size(), - subRegion.size(), - perforationData.getNumPerforationsGlobal(), - wellControls, - 0.0, // initialization done at t = 0 - resSinglePhaseFlowAccessors.get( fields::flow::pressure{} ), - resSingleFluidAccessors.get( fields::singlefluid::density{} ), - resElementRegion, - resElementSubRegion, - resElementIndex, - perfGravCoef, - wellElemGravCoef, - wellElemPressure ); - - // 4) Recompute the pressure-dependent properties - // Note: I am leaving that here because I would like to use the perforationRates (computed in UpdateState) - // to better initialize the rates - updateSubRegionState( subRegion ); - - string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); - SingleFluidBase & fluid = subRegion.getConstitutiveModel< SingleFluidBase >( fluidName ); - arrayView2d< real64 const > const & wellElemDens = fluid.density(); - // 5) Estimate the well rates - RateInitializationKernel::launch( subRegion.size(), - wellControls, - 0.0, // initialization done at t = 0 - wellElemDens, - connRate ); + if( time_n <= 0.0 || (wellControls.isWellOpen( time_n ) && !hasNonZero( connRate ) ) ) + { + // TODO: change the way we access the flowSolver here + SinglePhaseBase const & flowSolver = getParent().getGroup< SinglePhaseBase >( getFlowSolverName() ); + PresInitializationKernel::SinglePhaseFlowAccessors resSinglePhaseFlowAccessors( meshLevel.getElemManager(), flowSolver.getName() ); + PresInitializationKernel::SingleFluidAccessors resSingleFluidAccessors( meshLevel.getElemManager(), flowSolver.getName() ); + + // 1) Loop over all perforations to compute an average density + // 2) Initialize the reference pressure + // 3) Estimate the pressures in the well elements using the average density + PresInitializationKernel:: + launch( perforationData.size(), + subRegion.size(), + perforationData.getNumPerforationsGlobal(), + wellControls, + 0.0, // initialization done at t = 0 + resSinglePhaseFlowAccessors.get( fields::flow::pressure{} ), + resSingleFluidAccessors.get( fields::singlefluid::density{} ), + resElementRegion, + resElementSubRegion, + resElementIndex, + perfGravCoef, + wellElemGravCoef, + wellElemPressure ); + + // 4) Recompute the pressure-dependent properties + // Note: I am leaving that here because I would like to use the perforationRates (computed in UpdateState) + // to better initialize the rates + updateSubRegionState( subRegion ); + + string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); + SingleFluidBase & fluid = subRegion.getConstitutiveModel< SingleFluidBase >( fluidName ); + arrayView2d< real64 const > const & wellElemDens = fluid.density(); + + // 5) Estimate the well rates + RateInitializationKernel::launch( subRegion.size(), + wellControls, + 0.0, // initialization done at t = 0 + wellElemDens, + connRate ); + } } ); } ); } +void SinglePhaseWell::shutDownWell( real64 const time_n, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + string const wellDofKey = dofManager.getKey( wellElementDofName() ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & regionNames ) + { + + ElementRegionManager const & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, + [&]( localIndex const, + WellElementSubRegion const & subRegion ) + { + + // if the well is open, we don't have to do anything, so we just return + WellControls const & wellControls = getWellControls( subRegion ); + if( wellControls.isWellOpen( time_n ) ) + { + return; + } + + globalIndex const rankOffset = dofManager.rankOffset(); + + arrayView1d< integer const > const ghostRank = + subRegion.getReference< array1d< integer > >( ObjectManagerBase::viewKeyStruct::ghostRankString() ); + arrayView1d< globalIndex const > const dofNumber = + subRegion.getReference< array1d< globalIndex > >( wellDofKey ); + + arrayView1d< real64 const > const pres = + subRegion.getField< fields::well::pressure >(); + arrayView1d< real64 const > const connRate = + subRegion.getField< fields::well::connectionRate >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + if( ghostRank[ei] >= 0 ) + { + return; + } + + globalIndex const dofIndex = dofNumber[ei]; + localIndex const localRow = dofIndex - rankOffset; + real64 rhsValue; + + // 4.1. Apply pressure value to the matrix/rhs + FieldSpecificationEqual::SpecifyFieldValue( dofIndex, + rankOffset, + localMatrix, + rhsValue, + pres[ei], // freeze the current pressure value + pres[ei] ); + localRhs[localRow] = rhsValue; + + // 4.2. Apply rate value to the matrix/rhs + FieldSpecificationEqual::SpecifyFieldValue( dofIndex + 1, + rankOffset, + localMatrix, + rhsValue, + connRate[ei], // freeze the current pressure value + connRate[ei] ); + localRhs[localRow + 1] = rhsValue; + + } ); + } ); + } ); +} + +void SinglePhaseWell::assembleSystem( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + string const wellDofKey = dofManager.getKey( wellElementDofName()); + + // assemble the accumulation term in the mass balance equations + assembleAccumulationTerms( time, dt, domain, dofManager, localMatrix, localRhs ); + + // then assemble the pressure relations between well elements + assemblePressureRelations( time, dt, domain, dofManager, localMatrix, localRhs ); + // then compute the perforation rates (later assembled by the coupled solver) + computePerforationRates( time, dt, domain ); + + // then assemble the flux terms in the mass balance equations + // get a reference to the degree-of-freedom numbers + // then assemble the flux terms in the mass balance equations + assembleFluxTerms( time, dt, domain, dofManager, localMatrix, localRhs ); + + // then apply a special treatment to the wells that are shut + shutDownWell( time, domain, dofManager, localMatrix, localRhs ); +} void SinglePhaseWell::assembleFluxTerms( real64 const & time_n, real64 const & dt, @@ -473,7 +570,7 @@ void SinglePhaseWell::assembleFluxTerms( real64 const & time_n, } void SinglePhaseWell::assemblePressureRelations( real64 const & time_n, - real64 const & dt, + real64 const & GEOS_UNUSED_PARAM( dt ), DomainPartition const & domain, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, @@ -521,7 +618,7 @@ void SinglePhaseWell::assemblePressureRelations( real64 const & time_n, subRegion.isLocallyOwned(), subRegion.getTopWellElementIndex(), wellControls, - time_n + dt, // controls evaluated with BHP/rate of the end of the time interval + time_n, // controls evaluated with BHP/rate of the beginning of the time interval wellElemDofNumber, wellElemGravCoef, nextWellElemIndex, @@ -536,16 +633,14 @@ void SinglePhaseWell::assemblePressureRelations( real64 const & time_n, // Note: if BHP control is not viable, we switch to TOTALVOLRATE // if TOTALVOLRATE is not viable, we switch to BHP - real64 const timeAtEndOfStep = time_n + dt; - if( wellControls.getControl() == WellControls::Control::BHP ) { - wellControls.switchToTotalRateControl( wellControls.getTargetTotalRate( timeAtEndOfStep ) ); + wellControls.switchToTotalRateControl( wellControls.getTargetTotalRate( time_n ) ); GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, GEOS_FMT( "Control switch for well {} from BHP constraint to rate constraint", subRegion.getName()) ); } else { - wellControls.switchToBHPControl( wellControls.getTargetBHP( timeAtEndOfStep ) ); + wellControls.switchToBHPControl( wellControls.getTargetBHP( time_n ) ); GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, GEOS_FMT( "Control switch for well {} from rate constraint to BHP constraint", subRegion.getName()) ); } } @@ -739,7 +834,7 @@ SinglePhaseWell::calculateResidualNorm( real64 const & time_n, subRegion, fluid, wellControls, - time_n + dt, + time_n, dt, m_nonlinearSolverParameters.m_minNormalizer, subRegionResidualNorm ); @@ -936,7 +1031,7 @@ void SinglePhaseWell::implicitStepComplete( real64 const & time_n, } void SinglePhaseWell::printRates( real64 const & time_n, - real64 const & dt, + real64 const & GEOS_UNUSED_PARAM( dt ), DomainPartition & domain ) { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -974,10 +1069,10 @@ void SinglePhaseWell::printRates( real64 const & time_n, if( m_writeCSV > 0 ) { outputFile.open( m_ratesOutputDir + "/" + wellControlsName + ".csv", std::ios_base::app ); - outputFile << time_n + dt; + outputFile << time_n; } - if( !wellControls.isWellOpen( time_n + dt ) ) + if( !wellControls.isWellOpen( time_n ) ) { GEOS_LOG( GEOS_FMT( "{}: well is shut", wellControlsName ) ); if( outputFile.is_open()) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp index a27be685524..014dddcf8ad 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp @@ -167,6 +167,22 @@ class SinglePhaseWell : public WellSolverBase */ virtual real64 updateSubRegionState( WellElementSubRegion & subRegion ) override; + /** + * @brief function to assemble the linear system matrix and rhs + * @param time the time at the beginning of the step + * @param dt the desired timestep + * @param domain the domain partition + * @param dofManager degree-of-freedom manager associated with the linear system + * @param matrix the system matrix + * @param rhs the system right-hand side vector + */ + virtual void assembleSystem( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override; + /** * @brief assembles the flux terms for all connections between well elements * @param time_n previous time value @@ -224,6 +240,19 @@ class SinglePhaseWell : public WellSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) override; + /* + * @brief apply a special treatment to the wells that are shut + * @param time_n the time at the previous converged time step + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param matrix the system matrix + * @param rhs the system right-hand side vector + */ + void shutDownWell( real64 const time_n, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); struct viewKeyStruct : WellSolverBase::viewKeyStruct { static constexpr char const * dofFieldString() { return "singlePhaseWellVars"; } @@ -253,7 +282,7 @@ class SinglePhaseWell : public WellSolverBase * @brief Initialize all the primary and secondary variables in all the wells * @param domain the domain containing the well manager to access individual wells */ - void initializeWells( DomainPartition & domain, real64 const & time_n, real64 const & dt ) override; + void initializeWells( DomainPartition & domain, real64 const & time_n ) override; /** * @brief Make sure that the well constraints are compatible diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index 47877476c2f..bced9c29913 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -483,4 +483,27 @@ bool WellControls::isWellOpen( real64 const & currentTime ) const return isOpen; } +void WellControls::setNextDtFromTables( real64 const currentTime, real64 & nextDt ) +{ + setNextDtFromTable( m_targetBHPTable, currentTime, nextDt ); + setNextDtFromTable( m_targetMassRateTable, currentTime, nextDt ); + setNextDtFromTable( m_targetPhaseRateTable, currentTime, nextDt ); + setNextDtFromTable( m_targetTotalRateTable, currentTime, nextDt ); + setNextDtFromTable( m_statusTable, currentTime, nextDt ); +} + +void WellControls::setNextDtFromTable( TableFunction const * table, real64 const currentTime, real64 & nextDt ) +{ + if( table ) + { + // small epsilon to make sure we land on the other side of table interval and pick up the right rate + real64 const eps = 1e-6; + real64 const dtLimit = (table->getCoord( ¤tTime, 0, TableFunction::InterpolationType::Upper ) - currentTime) * ( 1.0 + eps ); + if( dtLimit > eps && dtLimit < nextDt ) + { + nextDt = dtLimit; + } + } +} + } //namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index 9e202096385..9084155853e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -199,11 +199,6 @@ class WellControls : public dataRepository::Group { return m_rateSign * m_targetPhaseRateTable->evaluate( ¤tTime ); } - /** - * @brief Get the target phase name - * @return the target phase name - */ - const string & getTargetPhaseName() const { return m_targetPhaseName; } /** * @brief Get the target mass rate @@ -214,6 +209,11 @@ class WellControls : public dataRepository::Group return m_rateSign * m_targetMassRateTable->evaluate( ¤tTime ); } + /** + * @brief Get the target phase name + * @return the target phase name + */ + const string & getTargetPhaseName() const { return m_targetPhaseName; } /** * @brief Const accessor for the composition of the injection stream @@ -276,6 +276,13 @@ class WellControls : public dataRepository::Group */ real64 getInitialPressureCoefficient() const { return m_initialPressureCoefficient; } + /** + * @brief set next time step based on tables intervals + * @param[in] currentTime the current time + * @param[inout] nextDt the time step + */ + void setNextDtFromTables( real64 const currentTime, real64 & nextDt ); + ///@} /** @@ -335,6 +342,8 @@ class WellControls : public dataRepository::Group virtual void postInputInitialization() override; + void setNextDtFromTable( TableFunction const * table, real64 const currentTime, real64 & nextDt ); + private: /// Well type (as Type enum) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp index 1b7ed01bb4e..6a15b7a3e30 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp @@ -61,6 +61,11 @@ WellSolverBase::WellSolverBase( string const & name, setInputFlag( dataRepository::InputFlags::OPTIONAL ). setDescription( "Write rates into a CSV file" ); + this->registerWrapper( viewKeyStruct::timeStepFromTablesFlagString(), &m_timeStepFromTables ). + setApplyDefaultValue( 0 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "Choose time step to honor rates/bhp tables time intervals" ); + addLogLevel< logInfo::WellControl >(); addLogLevel< logInfo::Crossflow >(); } @@ -195,12 +200,12 @@ void WellSolverBase::setupDofs( DomainPartition const & domain, } void WellSolverBase::implicitStepSetup( real64 const & time_n, - real64 const & dt, + real64 const & GEOS_UNUSED_PARAM( dt ), DomainPartition & domain ) { // Initialize the primary and secondary variables for the first time step - initializeWells( domain, time_n, dt ); + initializeWells( domain, time_n ); } void WellSolverBase::updateState( DomainPartition & domain ) @@ -308,4 +313,29 @@ WellControls const & WellSolverBase::getWellControls( WellElementSubRegion const return this->getGroup< WellControls >( subRegion.getWellControlsName()); } +real64 WellSolverBase::setNextDt( real64 const & currentTime, const real64 & currentDt, geos::DomainPartition & domain ) +{ + real64 nextDt = PhysicsSolverBase::setNextDt( currentTime, currentDt, domain ); + + if( m_timeStepFromTables ) + { + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, + WellElementSubRegion & subRegion ) + { + WellControls & wellControls = getWellControls( subRegion ); + real64 const nextDt_orig = nextDt; + wellControls.setNextDtFromTables( currentTime, nextDt ); + if( m_nonlinearSolverParameters.getLogLevel() > 0 && nextDt < nextDt_orig ) + GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on tables coordinates = {}", getName(), nextDt )); + } ); + } ); + } + + return nextDt; +} + } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp index a02af416759..dd113892a28 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp @@ -254,6 +254,17 @@ class WellSolverBase : public PhysicsSolverBase real64 const & dt, DomainPartition & domain ) = 0; + /** + * @brief function to set the next time step size + * @param[in] currentTime the current time + * @param[in] currentDt the current time step size + * @param[in] domain the domain object + * @return the prescribed time step size + */ + virtual real64 setNextDt( real64 const & currentTime, + real64 const & currentDt, + DomainPartition & domain ) override; + /** * @brief Utility function to keep the well variables during a time step (used in poromechanics simulations) * @param[in] keepVariablesConstantDuringInitStep flag to tell the solver to freeze its primary variables during a time step @@ -267,6 +278,7 @@ class WellSolverBase : public PhysicsSolverBase static constexpr char const * fluidNamesString() { return "fluidNames"; } static constexpr char const * isThermalString() { return "isThermal"; } static constexpr char const * writeCSVFlagString() { return "writeCSV"; } + static constexpr char const * timeStepFromTablesFlagString() { return "timeStepFromTables"; } }; private: @@ -292,7 +304,7 @@ class WellSolverBase : public PhysicsSolverBase * @brief Initialize all the primary and secondary variables in all the wells * @param domain the domain containing the well manager to access individual wells */ - virtual void initializeWells( DomainPartition & domain, real64 const & time_n, real64 const & dt ) = 0; + virtual void initializeWells( DomainPartition & domain, real64 const & time_n ) = 0; /** * @brief Make sure that the well constraints are compatible @@ -330,6 +342,9 @@ class WellSolverBase : public PhysicsSolverBase integer m_writeCSV; string const m_ratesOutputDir; + // flag to enable time step selection base on rates/bhp tables coordinates + integer m_timeStepFromTables; + /// flag to freeze the initial state during initialization in coupled problems integer m_keepVariablesConstantDuringInitStep; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp index 01fb94d2e72..f362c2b9959 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp @@ -303,7 +303,7 @@ PressureRelationKernel:: localIndex const iwelemControl, integer const targetPhaseIndex, WellControls const & wellControls, - real64 const & timeAtEndOfStep, + real64 const & time, arrayView1d< globalIndex const > const & wellElemDofNumber, arrayView1d< real64 const > const & wellElemGravCoef, arrayView1d< localIndex const > const & nextWellElemIndex, @@ -319,10 +319,10 @@ PressureRelationKernel:: bool const isProducer = wellControls.isProducer(); WellControls::Control const currentControl = wellControls.getControl(); WellControls::Control const inputControl = wellControls.getInputControl(); - real64 const targetBHP = wellControls.getTargetBHP( timeAtEndOfStep ); - real64 const targetTotalRate = wellControls.getTargetTotalRate( timeAtEndOfStep ); - real64 const targetPhaseRate = wellControls.getTargetPhaseRate( timeAtEndOfStep ); - real64 const targetMassRate = wellControls.getTargetMassRate( timeAtEndOfStep ); + real64 const targetBHP = wellControls.getTargetBHP( time ); + real64 const targetTotalRate = wellControls.getTargetTotalRate( time ); + real64 const targetPhaseRate = wellControls.getTargetPhaseRate( time ); + real64 const targetMassRate = wellControls.getTargetMassRate( time ); // dynamic well control data real64 const & currentBHP = @@ -451,7 +451,7 @@ PressureRelationKernel:: localIndex const iwelemControl, \ integer const targetPhaseIndex, \ WellControls const & wellControls, \ - real64 const & timeAtEndOfStep, \ + real64 const & time, \ arrayView1d< globalIndex const > const & wellElemDofNumber, \ arrayView1d< real64 const > const & wellElemGravCoef, \ arrayView1d< localIndex const > const & nextWellElemIndex, \ @@ -700,16 +700,16 @@ RateInitializationKernel:: launch( localIndex const subRegionSize, integer const targetPhaseIndex, WellControls const & wellControls, - real64 const & currentTime, + real64 const & time, arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, arrayView2d< real64 const, multifluid::USD_FLUID > const & totalDens, arrayView1d< real64 > const & connRate ) { WellControls::Control const control = wellControls.getControl(); bool const isProducer = wellControls.isProducer(); - real64 const targetTotalRate = wellControls.getTargetTotalRate( currentTime ); - real64 const targetPhaseRate = wellControls.getTargetPhaseRate( currentTime ); - real64 const targetMassRate = wellControls.getTargetMassRate( currentTime ); + real64 const targetTotalRate = wellControls.getTargetTotalRate( time ); + real64 const targetPhaseRate = wellControls.getTargetPhaseRate( time ); + real64 const targetMassRate = wellControls.getTargetMassRate( time ); // Estimate the connection rates forAll< parallelDevicePolicy<> >( subRegionSize, [=] GEOS_HOST_DEVICE ( localIndex const iwelem ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp index 5722613372d..cb5962042e6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp @@ -203,7 +203,7 @@ struct PressureRelationKernel localIndex const iwelemControl, integer const targetPhaseIndex, WellControls const & wellControls, - real64 const & timeAtEndOfStep, + real64 const & time, arrayView1d< globalIndex const > const & wellElemDofNumber, arrayView1d< real64 const > const & wellElemGravCoef, arrayView1d< localIndex const > const & nextWellElemIndex, @@ -507,7 +507,7 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa WellElementSubRegion const & subRegion, constitutive::MultiFluidBase const & fluid, WellControls const & wellControls, - real64 const timeAtEndOfStep, + real64 const time, real64 const dt, real64 const minNormalizer ) : Base( rankOffset, @@ -523,10 +523,10 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa m_iwelemControl( subRegion.getTopWellElementIndex() ), m_isProducer( wellControls.isProducer() ), m_currentControl( wellControls.getControl() ), - m_targetBHP( wellControls.getTargetBHP( timeAtEndOfStep ) ), - m_targetTotalRate( wellControls.getTargetTotalRate( timeAtEndOfStep ) ), - m_targetPhaseRate( wellControls.getTargetPhaseRate( timeAtEndOfStep ) ), - m_targetMassRate( wellControls.getTargetMassRate( timeAtEndOfStep ) ), + m_targetBHP( wellControls.getTargetBHP( time ) ), + m_targetTotalRate( wellControls.getTargetTotalRate( time ) ), + m_targetPhaseRate( wellControls.getTargetPhaseRate( time ) ), + m_targetMassRate( wellControls.getTargetMassRate( time ) ), m_volume( subRegion.getElementVolume() ), m_phaseDens_n( fluid.phaseDensity_n() ), m_totalDens_n( fluid.totalDensity_n() ) @@ -705,7 +705,7 @@ class ResidualNormKernelFactory * @param[in] subRegion the well element subregion * @param[in] fluid the fluid model * @param[in] wellControls the controls - * @param[in] timeAtEndOfStep the time at the end of the step (time_n + dt) + * @param[in] time the time * @param[in] dt the time step size * @param[out] residualNorm the residual norm on the subRegion */ @@ -720,7 +720,7 @@ class ResidualNormKernelFactory WellElementSubRegion const & subRegion, constitutive::MultiFluidBase const & fluid, WellControls const & wellControls, - real64 const timeAtEndOfStep, + real64 const time, real64 const dt, real64 const minNormalizer, real64 (& residualNorm)[1] ) @@ -729,7 +729,7 @@ class ResidualNormKernelFactory arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); ResidualNormKernel kernel( rankOffset, localResidual, dofNumber, ghostRank, - numComp, numDof, targetPhaseIndex, subRegion, fluid, wellControls, timeAtEndOfStep, dt, minNormalizer ); + numComp, numDof, targetPhaseIndex, subRegion, fluid, wellControls, time, dt, minNormalizer ); ResidualNormKernel::launchLinf< POLICY >( subRegion.size(), kernel, residualNorm ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp index e3dae7a4616..116be9fa414 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp @@ -243,7 +243,7 @@ PressureRelationKernel:: bool const isLocallyOwned, localIndex const iwelemControl, WellControls const & wellControls, - real64 const & timeAtEndOfStep, + real64 const & time, arrayView1d< globalIndex const > const & wellElemDofNumber, arrayView1d< real64 const > const & wellElemGravCoef, arrayView1d< localIndex const > const & nextWellElemIndex, @@ -256,8 +256,8 @@ PressureRelationKernel:: // static well control data bool const isProducer = wellControls.isProducer(); WellControls::Control const currentControl = wellControls.getControl(); - real64 const targetBHP = wellControls.getTargetBHP( timeAtEndOfStep ); - real64 const targetRate = wellControls.getTargetTotalRate( timeAtEndOfStep ); + real64 const targetBHP = wellControls.getTargetBHP( time ); + real64 const targetRate = wellControls.getTargetTotalRate( time ); // dynamic well control data real64 const & currentBHP = diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp index 8269b9f3cce..7b8c0be89e2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp @@ -145,7 +145,7 @@ struct PressureRelationKernel bool const isLocallyOwned, localIndex const iwelemControl, WellControls const & wellControls, - real64 const & timeAtEndOfStep, + real64 const & time, arrayView1d< globalIndex const > const & wellElemDofNumber, arrayView1d< real64 const > const & wellElemGravCoef, arrayView1d< localIndex const > const & nextWellElemIndex, @@ -327,7 +327,7 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa WellElementSubRegion const & subRegion, constitutive::SingleFluidBase const & fluid, WellControls const & wellControls, - real64 const timeAtEndOfStep, + real64 const time, real64 const dt, real64 const minNormalizer ) : Base( rankOffset, @@ -339,8 +339,8 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa m_isLocallyOwned( subRegion.isLocallyOwned() ), m_iwelemControl( subRegion.getTopWellElementIndex() ), m_currentControl( wellControls.getControl() ), - m_targetBHP( wellControls.getTargetBHP( timeAtEndOfStep ) ), - m_targetRate( wellControls.getTargetTotalRate( timeAtEndOfStep ) ), + m_targetBHP( wellControls.getTargetBHP( time ) ), + m_targetRate( wellControls.getTargetTotalRate( time ) ), m_volume( subRegion.getElementVolume() ), m_density_n( fluid.density_n() ) {} @@ -443,7 +443,7 @@ class ResidualNormKernelFactory * @param[in] subRegion the well element subregion * @param[in] fluid the fluid model * @param[in] wellControls the controls - * @param[in] timeAtEndOfStep the time at the end of the step (time_n + dt) + * @param[in] time the time * @param[in] dt the time step size * @param[out] residualNorm the residual norm on the subRegion */ @@ -455,7 +455,7 @@ class ResidualNormKernelFactory WellElementSubRegion const & subRegion, constitutive::SingleFluidBase const & fluid, WellControls const & wellControls, - real64 const timeAtEndOfStep, + real64 const time, real64 const dt, real64 const minNormalizer, real64 (& residualNorm)[1] ) @@ -463,7 +463,7 @@ class ResidualNormKernelFactory arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); - ResidualNormKernel kernel( rankOffset, localResidual, dofNumber, ghostRank, subRegion, fluid, wellControls, timeAtEndOfStep, dt, minNormalizer ); + ResidualNormKernel kernel( rankOffset, localResidual, dofNumber, ghostRank, subRegion, fluid, wellControls, time, dt, minNormalizer ); ResidualNormKernel::launchLinf< POLICY >( subRegion.size(), kernel, residualNorm ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp index 6e225367a36..29a35e388dd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp @@ -165,7 +165,7 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa WellElementSubRegion const & subRegion, MultiFluidBase const & fluid, WellControls const & wellControls, - real64 const timeAtEndOfStep, + real64 const time, real64 const dt, real64 const minNormalizer ) : Base( rankOffset, @@ -180,10 +180,10 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa m_iwelemControl( subRegion.getTopWellElementIndex() ), m_isProducer( wellControls.isProducer() ), m_currentControl( wellControls.getControl() ), - m_targetBHP( wellControls.getTargetBHP( timeAtEndOfStep ) ), - m_targetTotalRate( wellControls.getTargetTotalRate( timeAtEndOfStep ) ), - m_targetPhaseRate( wellControls.getTargetPhaseRate( timeAtEndOfStep ) ), - m_targetMassRate( wellControls.getTargetMassRate( timeAtEndOfStep ) ), + m_targetBHP( wellControls.getTargetBHP( time ) ), + m_targetTotalRate( wellControls.getTargetTotalRate( time ) ), + m_targetPhaseRate( wellControls.getTargetPhaseRate( time ) ), + m_targetMassRate( wellControls.getTargetMassRate( time ) ), m_volume( subRegion.getElementVolume() ), m_phaseDens_n( fluid.phaseDensity_n() ), m_totalDens_n( fluid.totalDensity_n() ), @@ -390,7 +390,7 @@ class ResidualNormKernelFactory * @param[in] subRegion the well element subregion * @param[in] fluid the fluid model * @param[in] wellControls the controls - * @param[in] timeAtEndOfStep the time at the end of the step (time_n + dt) + * @param[in] time the time * @param[in] dt the time step size * @param[out] residualNorm the residual norm on the subRegion */ @@ -404,7 +404,7 @@ class ResidualNormKernelFactory WellElementSubRegion const & subRegion, MultiFluidBase const & fluid, WellControls const & wellControls, - real64 const timeAtEndOfStep, + real64 const time, real64 const dt, real64 const minNormalizer, real64 (& residualNorm)[2] ) @@ -418,7 +418,7 @@ class ResidualNormKernelFactory arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); kernelType kernel( rankOffset, localResidual, dofNumber, ghostRank, - targetPhaseIndex, subRegion, fluid, wellControls, timeAtEndOfStep, dt, minNormalizer ); + targetPhaseIndex, subRegion, fluid, wellControls, time, dt, minNormalizer ); kernelType::template launchLinf< POLICY >( subRegion.size(), kernel, residualNorm ); } ); } diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp index 3c21d1b6ec8..12ce1f92812 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp @@ -118,7 +118,7 @@ real64 ExplicitQDRateAndState::solverStep( real64 const & time_n, else { // Retry with updated time step - dtAdaptive = setNextDt( dtAdaptive, domain ); + dtAdaptive = setNextDt( time_n, dtAdaptive, domain ); } } // return last successful adaptive time step (passed along to setNextDt) @@ -290,9 +290,11 @@ void ExplicitQDRateAndState::evalTimestep( DomainPartition & domain ) } } -real64 ExplicitQDRateAndState::setNextDt( real64 const & currentDt, DomainPartition & domain ) +real64 ExplicitQDRateAndState::setNextDt( real64 const & currentTime, + real64 const & currentDt, + DomainPartition & domain ) { - GEOS_UNUSED_VAR( domain ); + GEOS_UNUSED_VAR( currentTime, domain ); real64 const nextDt = m_stepUpdateFactor*currentDt; if( m_successfulStep ) { diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp index 69769af9381..41948de83ec 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp @@ -52,7 +52,8 @@ class ExplicitQDRateAndState : public QDRateAndStateBase DomainPartition & domain ) override final; - virtual real64 setNextDt( real64 const & currentDt, + virtual real64 setNextDt( real64 const & currentTime, + real64 const & currentDt, DomainPartition & domain ) override final; /** diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp index 3d8ae4ec756..1d7eeac3c63 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp @@ -104,9 +104,11 @@ void ImplicitQDRateAndState::updateSlip( ElementSubRegionBase & subRegion, real6 } ); } -real64 ImplicitQDRateAndState::setNextDt( real64 const & currentDt, DomainPartition & domain ) +real64 ImplicitQDRateAndState::setNextDt( real64 const & currentTime, + real64 const & currentDt, + DomainPartition & domain ) { - GEOS_UNUSED_VAR( currentDt ); + GEOS_UNUSED_VAR( currentTime, currentDt ); real64 maxSlipRate = 0.0; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp index 907612d0aef..02c112c8331 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp @@ -42,7 +42,8 @@ class ImplicitQDRateAndState : public QDRateAndStateBase constexpr static char const * targetSlipIncrementString() { return "targetSlipIncrement"; } }; - virtual real64 setNextDt( real64 const & currentDt, + virtual real64 setNextDt( real64 const & currentTime, + real64 const & currentDt, DomainPartition & domain ) override final; /** diff --git a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp index afa810d945a..b6fffd08c59 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp @@ -315,7 +315,7 @@ assembleCouplingTerms( real64 const time_n, ( wellControls.isInjector() ) && wellControls.isCrossflowEnabled() && getLogLevel() >= 1; // since detect crossflow requires communication, we detect it only if the logLevel is sufficiently high - if( !wellControls.isWellOpen( time_n + dt ) ) + if( !wellControls.isWellOpen( time_n ) ) { return; } diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp index 8517edb3dfe..12801fda2e3 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp @@ -309,26 +309,15 @@ class CoupledSolver : public PhysicsSolverBase } virtual real64 - setNextDtBasedOnStateChange( real64 const & currentDt, - DomainPartition & domain ) override + setNextDt( real64 const & currentTime, + real64 const & currentDt, + DomainPartition & domain ) override { - real64 nextDt = LvArray::NumericLimits< real64 >::max; + real64 nextDt = PhysicsSolverBase::setNextDt( currentTime, currentDt, domain ); forEachArgInTuple( m_solvers, [&]( auto & solver, auto ) { real64 const singlePhysicsNextDt = - solver->setNextDtBasedOnStateChange( currentDt, domain ); - nextDt = LvArray::math::min( singlePhysicsNextDt, nextDt ); - } ); - return nextDt; - } - - virtual real64 setNextDtBasedOnNewtonIter( real64 const & currentDt ) override - { - real64 nextDt = PhysicsSolverBase::setNextDtBasedOnNewtonIter( currentDt ); - forEachArgInTuple( m_solvers, [&]( auto & solver, auto ) - { - real64 const singlePhysicsNextDt = - solver->setNextDtBasedOnNewtonIter( currentDt ); + solver->setNextDt( currentTime, currentDt, domain ); nextDt = LvArray::math::min( singlePhysicsNextDt, nextDt ); } ); return nextDt; diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp index 3175f518666..c93c9c813fb 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp @@ -916,15 +916,16 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::resetStateToBeginningOfStep( D } template< typename POROMECHANICS_SOLVER > -real64 HydrofractureSolver< POROMECHANICS_SOLVER >::setNextDt( real64 const & currentDt, +real64 HydrofractureSolver< POROMECHANICS_SOLVER >::setNextDt( real64 const & currentTime, + real64 const & currentDt, DomainPartition & domain ) { - GEOS_UNUSED_VAR( domain ); + GEOS_UNUSED_VAR( currentTime, domain ); real64 nextDt = 0.0; if( m_numResolves[0] == 0 && m_numResolves[1] == 0 ) { - nextDt = this->setNextDtBasedOnNewtonIter( currentDt ); + nextDt = this->setNextDtBasedOnIterNumber( currentDt ); } else { diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp index 05a96830440..3faf5e11c36 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp @@ -122,7 +122,8 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) override; - virtual real64 setNextDt( real64 const & currentDt, + virtual real64 setNextDt( real64 const & currentTime, + real64 const & currentDt, DomainPartition & domain ) override; virtual void updateState( DomainPartition & domain ) override final; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.cpp index f83652d7455..c361c350159 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.cpp @@ -262,7 +262,7 @@ assembleCouplingTerms( real64 const time_n, PerforationData const * const perforationData = subRegion.getPerforationData(); WellControls const & wellControls = Base::wellSolver()->getWellControls( subRegion ); - if( !wellControls.isWellOpen( time_n + dt ) ) + if( !wellControls.isWellOpen( time_n ) ) { return; }