Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/byer3/well dt #3473

Merged
merged 3 commits into from
Dec 3, 2024
Merged

Feature/byer3/well dt #3473

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fa91c83
Always initialize well if previously closed
tjb-ltk Dec 3, 2024
4a7868e
re-introduce Jacobian shutin-logic for singlephase wells, inadverntly…
tjb-ltk Dec 3, 2024
d74337f
add init logic to singlephase well, fix compo version
tjb-ltk Dec 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 9 additions & 5 deletions src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -969,6 +969,12 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea
// TODO: change the way we access the flowSolver here
CompositionalMultiphaseBase const & flowSolver = getParent().getGroup< CompositionalMultiphaseBase >( getFlowSolverName() );

auto hasNonZeroCompRate = []( arrayView2d< real64 > const & arr ) {
return std::any_of( arr.begin(), arr.end(), []( real64 value ) {
return value > 0 || value < 0; // Check if the value is non-zero
} );
};

// loop over the wells
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
Expand All @@ -984,13 +990,11 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea
WellElementSubRegion & subRegion )
{
WellControls const & wellControls = getWellControls( subRegion );
PerforationData & perforationData = *subRegion.getPerforationData();
arrayView2d< real64 > 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 ) && !hasNonZeroCompRate( 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 >();
Expand Down
193 changes: 149 additions & 44 deletions src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -339,13 +339,13 @@ void SinglePhaseWell::initializeWells( DomainPartition & domain, real64 const &
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;
}

auto hasNonZeroCompRate = []( arrayView1d< real64 > const & arr ) {
return std::any_of( arr.begin(), arr.end(), []( real64 value ) {
return value > 0 || value < 0; // Check if the value is non-zero
} );
};

// loop over the wells
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & meshLevel,
Expand Down Expand Up @@ -381,49 +381,154 @@ 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 ) && !hasNonZeroCompRate( 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,
real64 const dt,
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, dt, domain, dofManager, localMatrix, localRhs );
}

void SinglePhaseWell::assembleFluxTerms( real64 const & time_n,
real64 const & dt,
Expand Down
31 changes: 31 additions & 0 deletions src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -224,6 +240,21 @@ 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 dt the time step size
* @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,
real64 const dt,
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"; }
Expand Down
Loading