Skip to content

Commit

Permalink
Some formulation options (#2546)
Browse files Browse the repository at this point in the history 
Add formulation options for multiphase flow solver:

useTotalMassEquation - enabled by default, when disabled - no total mass equation will be formed

useSimpleAccumulation - disabled by default, when enabled - the accumulation term will be assembled using simple form: PV * rho_c, mathematically equivalent to the current version, numerically might be slightly different. Simplifies accumulation linearization.

minCompDens - default value 1e-10 (same as before), allows experiments with minimum component ammount.
  • Loading branch information
@paveltomin
paveltomin authored Nov 6, 2023
1 parent 95aea4c commit 23abbb8
Show file tree
Hide file tree
Showing 26 changed files with 563 additions and 242 deletions.
2 changes: 2 additions & 0 deletions inputFiles/compositionalMultiphaseFlow/co2_flux_3d.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
discretization="fluidTPFA"
temperature="368.15"
useMass="1"
useTotalMassEquation="0"
useSimpleAccumulation="1"
targetRegions="{ region }">
<NonlinearSolverParameters
newtonTol="1.0e-5"
Expand Down
2 changes: 1 addition & 1 deletion integratedTests
Submodule integratedTests updated 244 files
21 changes: 21 additions & 0 deletions src/coreComponents/codingUtilities/Utilities.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,6 +392,27 @@ struct NoOpFunc
operator()( Ts && ... ) const {}
};

template< typename FLAGS_ENUM >
struct BitFlags
{
GEOS_HOST_DEVICE
void set( FLAGS_ENUM flag )
{
m_FlagValue |= (integer)flag;
}

GEOS_HOST_DEVICE
bool isSet( FLAGS_ENUM flag ) const
{
return (m_FlagValue & (integer)flag) == (integer)flag;
}

private:

using FlagsType = std::underlying_type_t< FLAGS_ENUM >;
FlagsType m_FlagValue{};

};

} // namespace geos

Expand Down
64 changes: 52 additions & 12 deletions src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,10 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name,
m_hasDispersion( 0 ),
m_keepFlowVariablesConstantDuringInitStep( 0 ),
m_minScalingFactor( 0.01 ),
m_allowCompDensChopping( 1 )
m_allowCompDensChopping( 1 ),
m_useTotalMassEquation( 1 ),
m_useSimpleAccumulation( 0 ),
m_minCompDens( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision )
{
//START_SPHINX_INCLUDE_00
this->registerWrapper( viewKeyStruct::inputTemperatureString(), &m_inputTemperature ).
Expand Down Expand Up @@ -122,6 +125,24 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name,
setApplyDefaultValue( 1 ).
setDescription( "Flag indicating whether local (cell-wise) chopping of negative compositions is allowed" );

this->registerWrapper( viewKeyStruct::useTotalMassEquationString(), &m_useTotalMassEquation ).
setSizedFromParent( 0 ).
setInputFlag( InputFlags::OPTIONAL ).
setApplyDefaultValue( 1 ).
setDescription( "Flag indicating whether total mass equation is used" );

this->registerWrapper( viewKeyStruct::useSimpleAccumulationString(), &m_useSimpleAccumulation ).
setSizedFromParent( 0 ).
setInputFlag( InputFlags::OPTIONAL ).
setApplyDefaultValue( 0 ).
setDescription( "Flag indicating whether simple accumulation form is used" );

this->registerWrapper( viewKeyStruct::minCompDensString(), &m_minCompDens ).
setSizedFromParent( 0 ).
setInputFlag( InputFlags::OPTIONAL ).
setApplyDefaultValue( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision ).
setDescription( "Minimum allowed global component density" );

}

void CompositionalMultiphaseBase::postProcessInput()
Expand Down Expand Up @@ -222,6 +243,7 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies )
m_numPhases = referenceFluid.numFluidPhases();
m_numComponents = referenceFluid.numFluidComponents();
}

// n_c components + one pressure ( + one temperature if needed )
m_numDofPerCell = m_isThermal ? m_numComponents + 2 : m_numComponents + 1;

Expand Down Expand Up @@ -1300,6 +1322,8 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom
createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
m_numPhases,
dofManager.rankOffset(),
m_useTotalMassEquation,
m_useSimpleAccumulation,
dofKey,
subRegion,
fluid,
Expand All @@ -1314,6 +1338,8 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom
createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
m_numPhases,
dofManager.rankOffset(),
m_useTotalMassEquation,
m_useSimpleAccumulation,
dofKey,
subRegion,
fluid,
Expand Down Expand Up @@ -1469,12 +1495,14 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time,

integer const fluidComponentId = fs.getComponent();
integer const numFluidComponents = m_numComponents;
integer const useTotalMassEquation = m_useTotalMassEquation;
forAll< parallelDevicePolicy<> >( targetSet.size(), [sizeScalingFactor,
targetSet,
rankOffset,
ghostRank,
fluidComponentId,
numFluidComponents,
useTotalMassEquation,
dofNumber,
rhsContributionArrayView,
localRhs] GEOS_HOST_DEVICE ( localIndex const a )
Expand All @@ -1486,15 +1514,23 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time,
return;
}

// for all "fluid components", we add the value to the total mass balance equation
globalIndex const totalMassBalanceRow = dofNumber[ei] - rankOffset;
localRhs[totalMassBalanceRow] += rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor
// here

// for all "fluid components" except the last one, we add the value to the component mass balance equation (shifted appropriately)
if( fluidComponentId < numFluidComponents - 1 )
if( useTotalMassEquation > 0 )
{
globalIndex const compMassBalanceRow = totalMassBalanceRow + fluidComponentId + 1; // component mass bal equations are shifted
// for all "fluid components", we add the value to the total mass balance equation
globalIndex const totalMassBalanceRow = dofNumber[ei] - rankOffset;
localRhs[totalMassBalanceRow] += rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the
// sizeScalingFactor
// here
if( fluidComponentId < numFluidComponents - 1 )
{
globalIndex const compMassBalanceRow = totalMassBalanceRow + fluidComponentId + 1; // component mass bal equations are shifted
localRhs[compMassBalanceRow] += rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the
// sizeScalingFactor here
}
}
else
{
globalIndex const compMassBalanceRow = dofNumber[ei] - rankOffset + fluidComponentId;
localRhs[compMassBalanceRow] += rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the
// sizeScalingFactor here
}
Expand Down Expand Up @@ -1919,9 +1955,9 @@ void CompositionalMultiphaseBase::chopNegativeDensities( DomainPartition & domai
{
for( integer ic = 0; ic < numComp; ++ic )
{
if( compDens[ei][ic] < minDensForDivision )
if( compDens[ei][ic] < m_minCompDens )
{
compDens[ei][ic] = minDensForDivision;
compDens[ei][ic] = m_minCompDens;
}
}
}
Expand Down Expand Up @@ -2191,6 +2227,7 @@ void CompositionalMultiphaseBase::saveIterationState( ElementSubRegionBase & sub

void CompositionalMultiphaseBase::updateState( DomainPartition & domain )
{
real64 maxDeltaPhaseVolFrac = 0.0;
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
MeshLevel & mesh,
arrayView1d< string const > const & regionNames )
Expand All @@ -2202,14 +2239,17 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain )
// update porosity, permeability, and solid internal energy
updatePorosityAndPermeability( subRegion );
// update all fluid properties
updateFluidState( subRegion );
real64 const deltaPhaseVolFrac = updateFluidState( subRegion );
maxDeltaPhaseVolFrac = LvArray::math::max( maxDeltaPhaseVolFrac, deltaPhaseVolFrac );
// for thermal, update solid internal energy
if( m_isThermal )
{
updateSolidInternalEnergyModel( subRegion );
}
} );
} );

GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "{}: Max deltaPhaseVolFrac = {}", getName(), maxDeltaPhaseVolFrac ) );
}

} // namespace geos
16 changes: 14 additions & 2 deletions src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -221,8 +221,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase
DofManager const & dofManager,
CRSMatrixView< real64, globalIndex const > const & localMatrix,
arrayView1d< real64 > const & localRhs ) const = 0;


/**@}*/

struct viewKeyStruct : FlowSolverBase::viewKeyStruct
Expand Down Expand Up @@ -253,6 +251,9 @@ class CompositionalMultiphaseBase : public FlowSolverBase
static constexpr char const * maxRelativePresChangeString() { return "maxRelativePressureChange"; }
static constexpr char const * maxRelativeTempChangeString() { return "maxRelativeTemperatureChange"; }
static constexpr char const * allowLocalCompDensChoppingString() { return "allowLocalCompDensityChopping"; }
static constexpr char const * useTotalMassEquationString() { return "useTotalMassEquation"; }
static constexpr char const * useSimpleAccumulationString() { return "useSimpleAccumulation"; }
static constexpr char const * minCompDensString() { return "minCompDens"; }

};

Expand Down Expand Up @@ -357,6 +358,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase

virtual void initializePostInitialConditionsPreSubGroups() override;

integer useTotalMassEquation() const { return m_useTotalMassEquation; }

protected:

virtual void postProcessInput() override;
Expand Down Expand Up @@ -445,6 +448,15 @@ class CompositionalMultiphaseBase : public FlowSolverBase
/// flag indicating whether local (cell-wise) chopping of negative compositions is allowed
integer m_allowCompDensChopping;

/// flag indicating whether total mass equation is used
integer m_useTotalMassEquation;

/// flag indicating whether simple accumulation form is used
integer m_useSimpleAccumulation;

/// minimum allowed global component density
real64 m_minCompDens;

/// name of the fluid constitutive model used as a reference for component/phase description
string m_referenceFluidModelName;

Expand Down
12 changes: 10 additions & 2 deletions src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,8 @@ void CompositionalMultiphaseFVM::initializePreSubGroups()
CompositionalMultiphaseBase::initializePreSubGroups();

m_linearSolverParameters.get().mgr.strategy = m_isThermal
? LinearSolverParameters::MGR::StrategyType::thermalCompositionalMultiphaseFVM
: LinearSolverParameters::MGR::StrategyType::compositionalMultiphaseFVM;
? LinearSolverParameters::MGR::StrategyType::thermalCompositionalMultiphaseFVM
: LinearSolverParameters::MGR::StrategyType::compositionalMultiphaseFVM;

DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
Expand Down Expand Up @@ -145,6 +145,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt,
dofManager.rankOffset(),
elemDofKey,
m_hasCapPressure,
m_useTotalMassEquation,
getName(),
mesh.getElemManager(),
stencilWrapper,
Expand All @@ -161,6 +162,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt,
dofManager.rankOffset(),
elemDofKey,
m_hasCapPressure,
m_useTotalMassEquation,
fluxApprox.upwindingParams(),
getName(),
mesh.getElemManager(),
Expand All @@ -184,6 +186,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt,
elemDofKey,
m_hasDiffusion,
m_hasDispersion,
m_useTotalMassEquation,
getName(),
mesh.getElemManager(),
stencilWrapper,
Expand All @@ -201,6 +204,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt,
elemDofKey,
m_hasDiffusion,
m_hasDispersion,
m_useTotalMassEquation,
getName(),
mesh.getElemManager(),
stencilWrapper,
Expand Down Expand Up @@ -248,6 +252,7 @@ void CompositionalMultiphaseFVM::assembleStabilizedFluxTerms( real64 const dt,
dofManager.rankOffset(),
elemDofKey,
m_hasCapPressure,
m_useTotalMassEquation,
getName(),
mesh.getElemManager(),
stencilWrapper,
Expand Down Expand Up @@ -895,6 +900,7 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n,
createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
m_numPhases,
dofManager.rankOffset(),
m_useTotalMassEquation,
elemDofKey,
getName(),
faceManager,
Expand All @@ -912,6 +918,7 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n,
createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
m_numPhases,
dofManager.rankOffset(),
m_useTotalMassEquation,
elemDofKey,
getName(),
faceManager,
Expand Down Expand Up @@ -992,6 +999,7 @@ void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time,
m_numPhases,
waterPhaseIndex,
allowAllPhasesIntoAquifer,
m_useTotalMassEquation,
stencil,
dofManager.rankOffset(),
elemDofNumber.toNestedViewConst(),
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -407,6 +407,7 @@ void CompositionalMultiphaseHybridFVM::assembleFluxTerms( real64 const dt,
dofManager.rankOffset(),
lengthTolerance,
dt,
m_useTotalMassEquation,
localMatrix,
localRhs );

Expand Down
Loading

0 comments on commit 23abbb8

Please sign in to comment.