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

fix: use mass and energy consistently for single phase solvers #3485

Open
wants to merge 8 commits into
base: develop
Choose a base branch
from
Open

fix: use mass and energy consistently for single phase solvers #3485

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
8b64891
wip
Dec 8, 2024
acf8b26
Update SinglePhasePoromechanicsEFEM.hpp
paveltomin Dec 9, 2024
1451422
Update SinglePhasePoromechanicsEFEM_impl.hpp
paveltomin Dec 9, 2024
3012d47
no need for solid, this is a fracture
Dec 9, 2024
105335e
bug fix
Dec 9, 2024
6fc9ab2
bug fix
Dec 9, 2024
7d30086
revert
Dec 9, 2024
6e49781
separate fixed stress
Dec 11, 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
2 changes: 2 additions & 0 deletions src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,6 +168,8 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies )
if( m_isThermal )
{
subRegion.registerField< fields::flow::energy >( getName() );
subRegion.registerField< fields::flow::dEnergy_dPressure >( getName() );
subRegion.registerField< fields::flow::dEnergy_dTemperature >( getName() );
subRegion.registerField< fields::flow::energy_n >( getName() );
}
} );
Expand Down
34 changes: 13 additions & 21 deletions src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -241,46 +241,38 @@ DECLARE_FIELD( temperatureScalingFactor,
NO_WRITE,
"Scaling factors for temperature" );

DECLARE_FIELD( mass,
Copy link
Contributor Author

Choose a reason for hiding this comment

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

mass went to single phase file since it is used there, not for compositional

"mass",
DECLARE_FIELD( energy,
"energy",
array1d< real64 >,
0,
LEVEL_0,
WRITE_AND_READ,
"Mass" );
"Energy" );

DECLARE_FIELD( mass_n,
"mass_n",
DECLARE_FIELD( dEnergy_dPressure,
"dEnergy_dPressure",
array1d< real64 >,
0,
NOPLOT,
WRITE_AND_READ,
"Mass at the previous converged time step" );
NO_WRITE,
"Derivative of energy with respect to pressure" );

DECLARE_FIELD( energy,
"energy",
DECLARE_FIELD( dEnergy_dTemperature,
"dEnergy_dTemperature",
array1d< real64 >,
0,
LEVEL_0,
WRITE_AND_READ,
"Energy" );
NOPLOT,
NO_WRITE,
"Derivative of energy with respect to temperature" );

DECLARE_FIELD( energy_n,
"energy_n",
array1d< real64 >,
0,
NOPLOT,
WRITE_AND_READ,
NO_WRITE,
"Energy at the previous converged time step" );

DECLARE_FIELD( massCreated,
"massCreated",
array1d< real64 >,
0,
LEVEL_1,
WRITE_AND_READ,
"The amount of remaining mass that was introduced when the SurfaceElement was created." );

}

}
Expand Down
72 changes: 55 additions & 17 deletions src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,8 +76,6 @@ SinglePhaseBase::SinglePhaseBase( const string & name,

void SinglePhaseBase::registerDataOnMesh( Group & meshBodies )
{
using namespace fields::flow;

FlowSolverBase::registerDataOnMesh( meshBodies );

m_numDofPerCell = m_isThermal ? 2 : 1;
Expand All @@ -93,18 +91,19 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies )
[&]( localIndex const,
ElementSubRegionBase & subRegion )
{
subRegion.registerField< deltaVolume >( getName() );
subRegion.registerField< fields::flow::deltaVolume >( getName() );

subRegion.registerField< mobility >( getName() );
subRegion.registerField< dMobility_dPressure >( getName() );
subRegion.registerField< fields::flow::mobility >( getName() );
subRegion.registerField< fields::flow::dMobility_dPressure >( getName() );

subRegion.registerField< fields::flow::mass >( getName() );
subRegion.registerField< fields::flow::dMass_dPressure >( getName() );
subRegion.registerField< fields::flow::mass_n >( getName() );


if( m_isThermal )
{
subRegion.registerField< dMobility_dTemperature >( getName() );
subRegion.registerField< fields::flow::dMobility_dTemperature >( getName() );
subRegion.registerField< fields::flow::dMass_dTemperature >( getName() );
}
} );

Expand All @@ -117,11 +116,11 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies )

FaceManager & faceManager = mesh.getFaceManager();
{
faceManager.registerField< facePressure >( getName() );
faceManager.registerField< fields::flow::facePressure >( getName() );

if( m_isThermal )
{
faceManager.registerField< faceTemperature >( getName() );
faceManager.registerField< fields::flow::faceTemperature >( getName() );
}
}
} );
Expand Down Expand Up @@ -272,11 +271,13 @@ void SinglePhaseBase::updateMass( ElementSubRegionBase & subRegion ) const
GEOS_MARK_FUNCTION;

arrayView1d< real64 > const mass = subRegion.getField< fields::flow::mass >();
arrayView1d< real64 > const dMass_dP = subRegion.getField< fields::flow::dMass_dPressure >();
arrayView1d< real64 > const mass_n = subRegion.getField< fields::flow::mass_n >();

CoupledSolidBase const & porousSolid =
getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
arrayView2d< real64 const > const porosity = porousSolid.getPorosity();
arrayView2d< real64 const > const dPorosity_dP = porousSolid.getDporosity_dPressure();
arrayView2d< real64 const > const porosity_n = porousSolid.getPorosity_n();

arrayView1d< real64 const > const volume = subRegion.getElementVolume();
Expand All @@ -285,39 +286,78 @@ void SinglePhaseBase::updateMass( ElementSubRegionBase & subRegion ) const
SingleFluidBase & fluid =
getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
arrayView2d< real64 const > const density = fluid.density();
arrayView2d< real64 const > const dDensity_dP = fluid.dDensity_dPressure();
arrayView2d< real64 const > const density_n = fluid.density_n();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
mass[ei] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * density[ei][0];
real64 const vol = volume[ei] + deltaVolume[ei];
mass[ei] = porosity[ei][0] * density[ei][0] * vol;
dMass_dP[ei] = ( dPorosity_dP[ei][0] * density[ei][0] + porosity[ei][0] * dDensity_dP[ei][0] ) * vol;
if( isZero( mass_n[ei] ) ) // this is a hack for hydrofrac cases
{
mass_n[ei] = porosity_n[ei][0] * volume[ei] * density_n[ei][0]; // initialize newly created element mass
}
} );

if( m_isThermal )
{
arrayView1d< real64 > const dMass_dT = subRegion.getField< fields::flow::dMass_dTemperature >();
arrayView2d< real64 const > const dPorosity_dT = porousSolid.getDporosity_dTemperature();
arrayView2d< real64 const > const dDensity_dT = fluid.dDensity_dTemperature();
forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
real64 const vol = volume[ei] + deltaVolume[ei];
dMass_dT[ei] = ( dPorosity_dT[ei][0] * density[ei][0] + dDensity_dT[ei][0] * dDensity_dP[ei][0] ) * vol;
} );
}
}

void SinglePhaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const
{
GEOS_MARK_FUNCTION;

arrayView1d< real64 > const energy = subRegion.getField< fields::flow::energy >();
arrayView1d< real64 > const dEnergy_dP = subRegion.getField< fields::flow::dEnergy_dPressure >();
arrayView1d< real64 > const dEnergy_dT = subRegion.getField< fields::flow::dEnergy_dTemperature >();

arrayView1d< real64 const > const volume = subRegion.getElementVolume();
arrayView1d< real64 > const deltaVolume = subRegion.getField< fields::flow::deltaVolume >();

CoupledSolidBase const & porousSolid =
getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
arrayView2d< real64 const > const porosity = porousSolid.getPorosity();
arrayView2d< real64 const > const dPorosity_dP = porousSolid.getDporosity_dPressure();
arrayView2d< real64 const > const dPorosity_dT = porousSolid.getDporosity_dTemperature();
arrayView2d< real64 const > const rockInternalEnergy = porousSolid.getInternalEnergy();

arrayView1d< real64 const > const volume = subRegion.getElementVolume();
arrayView1d< real64 > const deltaVolume = subRegion.getField< fields::flow::deltaVolume >();
arrayView2d< real64 const > const dRockInternalEnergy_dT = porousSolid.getDinternalEnergy_dTemperature();

SingleFluidBase & fluid =
getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
arrayView2d< real64 const > const density = fluid.density();
arrayView2d< real64 const > const dDensity_dP = fluid.dDensity_dPressure();
arrayView2d< real64 const > const dDensity_dT = fluid.dDensity_dTemperature();
arrayView2d< real64 const > const fluidInternalEnergy = fluid.internalEnergy();
arrayView2d< real64 const > const dFluidInternalEnergy_dP = fluid.dInternalEnergy_dPressure();
arrayView2d< real64 const > const dFluidInternalEnergy_dT = fluid.dInternalEnergy_dTemperature();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
energy[ei] = ( volume[ei] + deltaVolume[ei] ) *
( porosity[ei][0] * density[ei][0] * fluidInternalEnergy[ei][0] + ( 1.0 - porosity[ei][0] ) * rockInternalEnergy[ei][0] );
real64 const vol = volume[ei] + deltaVolume[ei];
energy[ei] = vol *
( porosity[ei][0] * density[ei][0] * fluidInternalEnergy[ei][0] +
( 1.0 - porosity[ei][0] ) * rockInternalEnergy[ei][0] );
dEnergy_dP[ei] = vol *
( dPorosity_dP[ei][0] * density[ei][0] * fluidInternalEnergy[ei][0] +
porosity[ei][0] * dDensity_dP[ei][0] * fluidInternalEnergy[ei][0] +
porosity[ei][0] * density[ei][0] * dFluidInternalEnergy_dP[ei][0] -
dPorosity_dP[ei][0] * rockInternalEnergy[ei][0] );
dEnergy_dT[ei] = vol *
( dPorosity_dT[ei][0] * density[ei][0] * fluidInternalEnergy[ei][0] +
porosity[ei][0] * dDensity_dT[ei][0] * fluidInternalEnergy[ei][0] +
porosity[ei][0] * density[ei][0] * dFluidInternalEnergy_dT[ei][0] -
dPorosity_dT[ei][0] * rockInternalEnergy[ei][0] +
( 1.0 - porosity[ei][0] ) * dRockInternalEnergy_dT[ei][0] );
} );
}

Expand Down Expand Up @@ -594,8 +634,6 @@ void SinglePhaseBase::initializeFluidState( MeshLevel & mesh, arrayView1d< strin

// 2. save the initial density (for use in the single-phase poromechanics solver to compute the deltaBodyForce)
fluid.initializeState();

updateMass( subRegion );
Copy link
Contributor Author

Choose a reason for hiding this comment

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

this call is redundant, updateMass already called in updateFluidState few lines above

} );
}

Expand Down
11 changes: 0 additions & 11 deletions src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -438,13 +438,6 @@ void SinglePhaseBase::accumulationAssemblyLaunch( DofManager const & dofManager,
CRSMatrixView< real64, globalIndex const > const & localMatrix,
arrayView1d< real64 > const & localRhs )
{
geos::constitutive::SingleFluidBase const & fluid =
getConstitutiveModel< geos::constitutive::SingleFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) );
//START_SPHINX_INCLUDE_COUPLEDSOLID
geos::constitutive::CoupledSolidBase const & solid =
getConstitutiveModel< geos::constitutive::CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
//END_SPHINX_INCLUDE_COUPLEDSOLID

string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );

if( m_isThermal )
Expand All @@ -454,8 +447,6 @@ void SinglePhaseBase::accumulationAssemblyLaunch( DofManager const & dofManager,
createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(),
dofKey,
subRegion,
fluid,
solid,
localMatrix,
localRhs );
}
Expand All @@ -466,8 +457,6 @@ void SinglePhaseBase::accumulationAssemblyLaunch( DofManager const & dofManager,
createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(),
dofKey,
subRegion,
fluid,
solid,
localMatrix,
localRhs );
}
Expand Down
40 changes: 40 additions & 0 deletions src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,46 @@ namespace fields
namespace flow
{

DECLARE_FIELD( mass,
"mass",
array1d< real64 >,
0,
LEVEL_0,
WRITE_AND_READ,
"Fluid mass" );

DECLARE_FIELD( dMass_dPressure,
"dMass_dPressure",
array1d< real64 >,
0,
NOPLOT,
NO_WRITE,
"Derivative of mass with respect to pressure" );

DECLARE_FIELD( dMass_dTemperature,
"dMass_dTemperature",
array1d< real64 >,
0,
NOPLOT,
NO_WRITE,
"Derivative of mass with respect to temperature" );

DECLARE_FIELD( mass_n,
"mass_n",
array1d< real64 >,
0,
NOPLOT,
NO_WRITE,
"Fluid mass at the previous converged time step" );

DECLARE_FIELD( massCreated,
"massCreated",
array1d< real64 >,
0,
LEVEL_1,
WRITE_AND_READ,
"The amount of remaining mass that was introduced when the SurfaceElement was created." );

DECLARE_FIELD( mobility,
"mobility",
array1d< real64 >,
Expand Down
Loading
Loading