Skip to content

Moved fluid variables at time n to multifluid constitutive model #1827

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

Merged
francoishamon merged 6 commits into from
Mar 25, 2022
Merged

Moved fluid variables at time n to multifluid constitutive model #1827

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
835eadb
moved old variables to multifluid constitutive model
francoishamon Mar 17, 2022
d7492b3
Merge branch 'develop' into refactor/hamon/moveOldToConstitutive
francoishamon Mar 17, 2022
0bfcafb
Merge branch 'develop' into refactor/hamon/moveOldToConstitutive
francoishamon Mar 23, 2022
6d472b4
fixed compilation
francoishamon Mar 23, 2022
4a043cf
Merge branch 'develop' into refactor/hamon/moveOldToConstitutive
francoishamon Mar 25, 2022
b547abe
updated integratedTests submodule
francoishamon Mar 25, 2022
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
2 changes: 1 addition & 1 deletion integratedTests
Open in desktop
Submodule integratedTests updated from d79950 to 4bbb2f
54 changes: 54 additions & 0 deletions src/coreComponents/constitutive/fluid/MultiFluidBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ MultiFluidBase::MultiFluidBase( string const & name, Group * const parent )
registerExtrinsicData( extrinsicMeshData::multifluid::dPhaseFraction{}, &m_phaseFraction.derivs );

registerExtrinsicData( extrinsicMeshData::multifluid::phaseDensity{}, &m_phaseDensity.value );
registerExtrinsicData( extrinsicMeshData::multifluid::phaseDensityOld{}, &m_phaseDensityOld );
registerExtrinsicData( extrinsicMeshData::multifluid::dPhaseDensity{}, &m_phaseDensity.derivs );

registerExtrinsicData( extrinsicMeshData::multifluid::phaseMassDensity{}, &m_phaseMassDensity.value );
Expand All @@ -63,15 +64,19 @@ MultiFluidBase::MultiFluidBase( string const & name, Group * const parent )
registerExtrinsicData( extrinsicMeshData::multifluid::dPhaseViscosity{}, &m_phaseViscosity.derivs );

registerExtrinsicData( extrinsicMeshData::multifluid::phaseEnthalpy{}, &m_phaseEnthalpy.value );
registerExtrinsicData( extrinsicMeshData::multifluid::phaseEnthalpyOld{}, &m_phaseEnthalpyOld );
registerExtrinsicData( extrinsicMeshData::multifluid::dPhaseEnthalpy{}, &m_phaseEnthalpy.derivs );

registerExtrinsicData( extrinsicMeshData::multifluid::phaseInternalEnergy{}, &m_phaseInternalEnergy.value );
registerExtrinsicData( extrinsicMeshData::multifluid::phaseInternalEnergyOld{}, &m_phaseInternalEnergyOld );
registerExtrinsicData( extrinsicMeshData::multifluid::dPhaseInternalEnergy{}, &m_phaseInternalEnergy.derivs );

registerExtrinsicData( extrinsicMeshData::multifluid::phaseCompFraction{}, &m_phaseCompFraction.value );
registerExtrinsicData( extrinsicMeshData::multifluid::phaseCompFractionOld{}, &m_phaseCompFractionOld );
registerExtrinsicData( extrinsicMeshData::multifluid::dPhaseCompFraction{}, &m_phaseCompFraction.derivs );

registerExtrinsicData( extrinsicMeshData::multifluid::totalDensity{}, &m_totalDensity.value );
registerExtrinsicData( extrinsicMeshData::multifluid::totalDensityOld{}, &m_totalDensityOld );
registerExtrinsicData( extrinsicMeshData::multifluid::dTotalDensity{}, &m_totalDensity.derivs );

registerExtrinsicData( extrinsicMeshData::multifluid::initialTotalMassDensity{}, &m_initialTotalMassDensity );
Expand All @@ -88,6 +93,7 @@ void MultiFluidBase::resizeFields( localIndex const size, localIndex const numPt
m_phaseFraction.derivs.resize( size, numPts, numPhase, numDof );

m_phaseDensity.value.resize( size, numPts, numPhase );
m_phaseDensityOld.resize( size, numPts, numPhase );
m_phaseDensity.derivs.resize( size, numPts, numPhase, numDof );

m_phaseMassDensity.value.resize( size, numPts, numPhase );
Expand All @@ -97,15 +103,19 @@ void MultiFluidBase::resizeFields( localIndex const size, localIndex const numPt
m_phaseViscosity.derivs.resize( size, numPts, numPhase, numDof );

m_phaseEnthalpy.value.resize( size, numPts, numPhase );
m_phaseEnthalpyOld.resize( size, numPts, numPhase );
m_phaseEnthalpy.derivs.resize( size, numPts, numPhase, numDof );

m_phaseInternalEnergy.value.resize( size, numPts, numPhase );
m_phaseInternalEnergyOld.resize( size, numPts, numPhase );
m_phaseInternalEnergy.derivs.resize( size, numPts, numPhase, numDof );

m_phaseCompFraction.value.resize( size, numPts, numPhase, numComp );
m_phaseCompFractionOld.resize( size, numPts, numPhase, numComp );
m_phaseCompFraction.derivs.resize( size, numPts, numPhase, numComp, numDof );

m_totalDensity.value.resize( size, numPts );
m_totalDensityOld.resize( size, numPts );
m_totalDensity.derivs.resize( size, numPts, numDof );

m_initialTotalMassDensity.resize( size, numPts );
Expand Down Expand Up @@ -193,8 +203,52 @@ void MultiFluidBase::initializeState( arrayView2d< real64 const, compflow::USD_P
}
}
} );

// initialize the "old" variables
saveConvergedState();
}

void MultiFluidBase::saveConvergedState() const
{
localIndex const numElem = m_initialTotalMassDensity.size( 0 );
localIndex const numGauss = m_initialTotalMassDensity.size( 1 );
integer const numPhase = m_phaseMassDensity.value.size( 2 );
integer const numComp = m_phaseCompFraction.value.size( 3 );
Comment on lines +213 to +216
Copy link
Contributor

@TotoGaz TotoGaz Mar 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some day we can get rid of these hard coded indices. 😉


FluidProp::ViewTypeConst const totalDensity = m_totalDensity.toViewConst();
PhaseProp::ViewTypeConst const phaseDensity = m_phaseDensity.toViewConst();
PhaseProp::ViewTypeConst const phaseEnthalpy = m_phaseEnthalpy.toViewConst();
PhaseProp::ViewTypeConst const phaseInternalEnergy = m_phaseInternalEnergy.toViewConst();
PhaseComp::ViewTypeConst const phaseCompFraction = m_phaseCompFraction.toViewConst();

arrayView2d< real64, multifluid::USD_FLUID > const totalDensityOld = m_totalDensityOld.toView();
arrayView3d< real64, multifluid::USD_PHASE > const phaseDensityOld = m_phaseDensityOld.toView();
arrayView3d< real64, multifluid::USD_PHASE > const phaseEnthalpyOld = m_phaseEnthalpyOld.toView();
arrayView3d< real64, multifluid::USD_PHASE > const phaseInternalEnergyOld = m_phaseInternalEnergyOld.toView();
arrayView4d< real64, multifluid::USD_PHASE_COMP > const phaseCompFractionOld = m_phaseCompFractionOld.toView();

forAll< parallelDevicePolicy<> >( numElem, [=] GEOSX_HOST_DEVICE ( localIndex const k )
{
for( localIndex q = 0; q < numGauss; ++q )
{
totalDensityOld[k][q] = totalDensity.value[k][q];
for( integer ip = 0; ip < numPhase; ++ip )
{
phaseDensityOld[k][q][ip] = phaseDensity.value[k][q][ip];
phaseEnthalpyOld[k][q][ip] = phaseEnthalpy.value[k][q][ip];
phaseInternalEnergyOld[k][q][ip] = phaseInternalEnergy.value[k][q][ip];
for( integer ic = 0; ic < numComp; ++ic )
{
phaseCompFractionOld[k][q][ip][ic] = phaseCompFraction.value[k][q][ip][ic];
}
Comment on lines +237 to +243
Copy link
Contributor

@TotoGaz TotoGaz Mar 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So we do not store the derivatives, right?

Copy link
Collaborator

@CusiniM CusiniM Mar 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

well, for the quantities at the previous time-step we only need the value non its derivatives.

}
}
} );


}


} // namespace constitutive

} // namespace geosx
28 changes: 28 additions & 0 deletions src/coreComponents/constitutive/fluid/MultiFluidBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,9 @@ class MultiFluidBase : public ConstitutiveBase
arrayView3d< real64 const, multifluid::USD_PHASE > phaseDensity() const
{ return m_phaseDensity.value; }

arrayView3d< real64 const, multifluid::USD_PHASE > phaseDensityOld() const
{ return m_phaseDensityOld; }

arrayView4d< real64 const, multifluid::USD_PHASE_DC > dPhaseDensity() const
{ return m_phaseDensity.derivs; }

Expand All @@ -139,12 +142,18 @@ class MultiFluidBase : public ConstitutiveBase
arrayView4d< real64 const, multifluid::USD_PHASE_COMP > phaseCompFraction() const
{ return m_phaseCompFraction.value; }

arrayView4d< real64 const, multifluid::USD_PHASE_COMP > phaseCompFractionOld() const
{ return m_phaseCompFractionOld; }

arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > dPhaseCompFraction() const
{ return m_phaseCompFraction.derivs; }

arrayView2d< real64 const, multifluid::USD_FLUID > totalDensity() const
{ return m_totalDensity.value; }

arrayView2d< real64 const, multifluid::USD_FLUID > totalDensityOld() const
{ return m_totalDensityOld; }

arrayView3d< real64 const, multifluid::USD_FLUID_DC > dTotalDensity() const
{ return m_totalDensity.derivs; }

Expand All @@ -154,12 +163,18 @@ class MultiFluidBase : public ConstitutiveBase
arrayView3d< real64 const, multifluid::USD_PHASE > phaseEnthalpy() const
{ return m_phaseEnthalpy.value; }

arrayView3d< real64 const, multifluid::USD_PHASE > phaseEnthalpyOld() const
{ return m_phaseEnthalpyOld; }

arrayView4d< real64 const, multifluid::USD_PHASE_DC > dPhaseEnthalpy() const
{ return m_phaseEnthalpy.derivs; }

arrayView3d< real64 const, multifluid::USD_PHASE > phaseInternalEnergy() const
{ return m_phaseInternalEnergy.value; }

arrayView3d< real64 const, multifluid::USD_PHASE > phaseInternalEnergyOld() const
{ return m_phaseInternalEnergyOld; }

arrayView4d< real64 const, multifluid::USD_PHASE_DC > dPhaseInternalEnergy() const
{ return m_phaseInternalEnergy.derivs; }

Expand All @@ -169,6 +184,11 @@ class MultiFluidBase : public ConstitutiveBase
*/
virtual void initializeState( arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFraction ) const;

/**
* @brief Save the phase densities, component fractions, enthalpies and internal energies (for accumulation)
*/
virtual void saveConvergedState() const override;

struct viewKeyStruct : ConstitutiveBase::viewKeyStruct
{
static constexpr char const * componentNamesString() { return "componentNames"; }
Expand Down Expand Up @@ -537,6 +557,14 @@ class MultiFluidBase : public ConstitutiveBase
PhaseComp m_phaseCompFraction;
FluidProp m_totalDensity;

// backup data

array3d< real64, multifluid::LAYOUT_PHASE > m_phaseDensityOld;
array3d< real64, multifluid::LAYOUT_PHASE > m_phaseEnthalpyOld;
array3d< real64, multifluid::LAYOUT_PHASE > m_phaseInternalEnergyOld;
array4d< real64, multifluid::LAYOUT_PHASE_COMP > m_phaseCompFractionOld;
array2d< real64, multifluid::LAYOUT_FLUID > m_totalDensityOld;

// initial data (used to compute the body force in the poromechanics solver)

array2d< real64, multifluid::LAYOUT_FLUID > m_initialTotalMassDensity;
Expand Down
41 changes: 40 additions & 1 deletion src/coreComponents/constitutive/fluid/MultiFluidExtrinsicData.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,14 @@ EXTRINSIC_MESH_DATA_TRAIT( phaseDensity,
WRITE_AND_READ,
"Phase density" );

EXTRINSIC_MESH_DATA_TRAIT( phaseDensityOld,
"phaseDensityOld",
array3dLayoutPhase,
0,
NOPLOT,
WRITE_AND_READ,
"Phase density at the previous converged time step" );

EXTRINSIC_MESH_DATA_TRAIT( dPhaseDensity,
"dPhaseDensity",
array4dLayoutPhase_dC,
Expand Down Expand Up @@ -110,6 +118,14 @@ EXTRINSIC_MESH_DATA_TRAIT( phaseEnthalpy,
NO_WRITE,
"Phase enthalpy" );

EXTRINSIC_MESH_DATA_TRAIT( phaseEnthalpyOld,
"phaseEnthalpyOld",
array3dLayoutPhase,
0,
NOPLOT,
WRITE_AND_READ,
"Phase enthalpy at the previous converged time step" );

EXTRINSIC_MESH_DATA_TRAIT( dPhaseEnthalpy,
"dPhaseEnthalpy",
array4dLayoutPhase_dC,
Expand All @@ -126,6 +142,14 @@ EXTRINSIC_MESH_DATA_TRAIT( phaseInternalEnergy,
NO_WRITE,
"Phase internal energy" );

EXTRINSIC_MESH_DATA_TRAIT( phaseInternalEnergyOld,
"phaseInternalEnergyOld",
array3dLayoutPhase,
0,
NOPLOT,
WRITE_AND_READ,
"Phase internal energy at the previous converged time step" );

EXTRINSIC_MESH_DATA_TRAIT( dPhaseInternalEnergy,
"dPhaseInternalEnergy",
array4dLayoutPhase_dC,
Expand All @@ -142,6 +166,14 @@ EXTRINSIC_MESH_DATA_TRAIT( phaseCompFraction,
WRITE_AND_READ,
"Phase component fraction" );

EXTRINSIC_MESH_DATA_TRAIT( phaseCompFractionOld,
"phaseCompFractionOld",
array4dLayoutPhaseComp,
0,
NOPLOT,
WRITE_AND_READ,
"Phase component fraction at the previous converged time step" );

EXTRINSIC_MESH_DATA_TRAIT( dPhaseCompFraction,
"dPhaseCompFraction",
array5dLayoutPhaseComp_dC,
Expand All @@ -158,6 +190,14 @@ EXTRINSIC_MESH_DATA_TRAIT( totalDensity,
WRITE_AND_READ,
"Total density" );

EXTRINSIC_MESH_DATA_TRAIT( totalDensityOld,
"totalDensityOld",
array2dLayoutFluid,
0,
NOPLOT,
WRITE_AND_READ,
"Total density at the previous converged time step" );

EXTRINSIC_MESH_DATA_TRAIT( initialTotalMassDensity,
"initialTotalMassDensity",
array2dLayoutFluid,
Expand All @@ -174,7 +214,6 @@ EXTRINSIC_MESH_DATA_TRAIT( dTotalDensity,
NO_WRITE,
"Derivative of total density with respect to pressure, temperature, and global component fractions" );


}

}
Expand Down
4 changes: 2 additions & 2 deletions src/coreComponents/constitutive/solid/PorousSolid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,10 +148,10 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity,
real64 const & solidDensity,
real64 const & initialFluidTotalMassDensity,
arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & fluidPhaseDensity,
arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & fluidPhaseDensityOld,
arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & fluidPhaseDensityOld,
arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const & dFluidPhaseDensity,
arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_COMP - 2 > const & fluidPhaseCompFrac,
arraySlice2d< real64 const, compflow::USD_PHASE_COMP - 1 > const & fluidPhaseCompFracOld,
arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_COMP - 2 > const & fluidPhaseCompFracOld,
arraySlice3d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC -2 > const & dFluidPhaseCompFrac,
arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & fluidPhaseMassDensity,
arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const & dFluidPhaseMassDensity,
Expand Down
43 changes: 11 additions & 32 deletions src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -226,13 +226,8 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies )

subRegion.registerExtrinsicData< phaseVolumeFractionOld >( getName() ).
reference().resizeDimension< 1 >( m_numPhases );
subRegion.registerExtrinsicData< totalDensityOld >( getName() );
subRegion.registerExtrinsicData< phaseDensityOld >( getName() ).
reference().resizeDimension< 1 >( m_numPhases );
subRegion.registerExtrinsicData< phaseMobilityOld >( getName() ).
reference().resizeDimension< 1 >( m_numPhases );
subRegion.registerExtrinsicData< phaseComponentFractionOld >( getName() ).
reference().resizeDimension< 1, 2 >( m_numPhases, m_numComponents );

} );

Expand Down Expand Up @@ -1021,7 +1016,6 @@ void CompositionalMultiphaseBase::backupFields( MeshLevel & mesh,
{
GEOSX_MARK_FUNCTION;

integer const numComp = m_numComponents;
integer const numPhase = m_numPhases;

// backup some fields used in time derivative approximation
Expand All @@ -1033,26 +1027,13 @@ void CompositionalMultiphaseBase::backupFields( MeshLevel & mesh,

arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac =
subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFraction >();
arrayView2d< real64 const, compflow::USD_PHASE > const & phaseMob =
arrayView2d< real64 const, compflow::USD_PHASE > const phaseMob =
subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseMobility >();

string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() );
MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName );
arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity();
arrayView3d< real64 const, multifluid::USD_PHASE > const phaseDens = fluid.phaseDensity();
arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const phaseCompFrac = fluid.phaseCompFraction();

arrayView1d< real64 > const totalDensOld =
subRegion.getExtrinsicData< extrinsicMeshData::flow::totalDensityOld >();

arrayView2d< real64, compflow::USD_PHASE > const phaseDensOld =
subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseDensityOld >();
arrayView2d< real64, compflow::USD_PHASE > const phaseVolFracOld =
subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFractionOld >();
arrayView2d< real64, compflow::USD_PHASE > const phaseMobOld =
subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseMobilityOld >();
arrayView3d< real64, compflow::USD_PHASE_COMP > const phaseCompFracOld =
subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseComponentFractionOld >();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei )
{
Expand All @@ -1061,16 +1042,9 @@ void CompositionalMultiphaseBase::backupFields( MeshLevel & mesh,

for( integer ip = 0; ip < numPhase; ++ip )
{
phaseDensOld[ei][ip] = phaseDens[ei][0][ip];
phaseVolFracOld[ei][ip] = phaseVolFrac[ei][ip];
phaseMobOld[ei][ip] = phaseMob[ei][ip];

for( integer ic = 0; ic < numComp; ++ic )
{
phaseCompFracOld[ei][ip][ic] = phaseCompFrac[ei][0][ip][ic];
}
}
totalDensOld[ei] = totalDens[ei][0];
} );
} );
}
Expand Down Expand Up @@ -1647,20 +1621,25 @@ void CompositionalMultiphaseBase::implicitStepComplete( real64 const & time,
}
} );

// Step 3: save the converged solid state
// Step 3: save the converged fluid state
string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() );
MultiFluidBase const & fluidMaterial = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName );
fluidMaterial.saveConvergedState();

// Step 4: save the converged solid state
string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() );
CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName );
porousMaterial.saveConvergedState();

// Step 4: save converged state for the relperm model to handle hysteresis
// Step 5: save converged state for the relperm model to handle hysteresis
arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac =
subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFraction >();
string const & relPermName = subRegion.getReference< string >( viewKeyStruct::relPermNamesString() );
RelativePermeabilityBase & relPermMaterial =
RelativePermeabilityBase const & relPermMaterial =
getConstitutiveModel< RelativePermeabilityBase >( subRegion, relPermName );
relPermMaterial.saveConvergedPhaseVolFractionState( phaseVolFrac );

// Step 5: if capillary pressure is supported, send the converged porosity and permeability to the capillary pressure model
// Step 6: if capillary pressure is supported, send the converged porosity and permeability to the capillary pressure model
// note: this is needed when the capillary pressure depends on porosity and permeability (Leverett J-function for instance)
if( m_hasCapPressure )
{
Expand All @@ -1677,7 +1656,7 @@ void CompositionalMultiphaseBase::implicitStepComplete( real64 const & time,
capPressureMaterial.saveConvergedRockState( porosity, permeability );
}

// Step 6: if the thermal option is on, send the converged porosity and phase volume fraction to the thermal conductivity model
// Step 7: if the thermal option is on, send the converged porosity and phase volume fraction to the thermal conductivity model
// note: this is needed because the phaseVolFrac-weighted thermal conductivity treats phaseVolumeFraction explicitly for now
if( m_isThermal )
{
Expand Down
Loading