Skip to content

Commit

Permalink
merge fix
Browse files Browse the repository at this point in the history
  • Loading branch information
Pavel Tomin authored and Pavel Tomin committed Jan 23, 2025
1 parent ba262ba commit 6be1912
Show file tree
Hide file tree
Showing 10 changed files with 58 additions and 92 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -185,8 +185,6 @@ 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
Original file line number Diff line number Diff line change
Expand Up @@ -257,14 +257,6 @@ DECLARE_FIELD( energy_n,
NO_WRITE,
"Energy at the previous converged time step" );

DECLARE_FIELD( dEnergy,
"dEnergy",
array1d< real64 >, // TODO: change to array2d< real64 >
0,
NOPLOT,
NO_WRITE,
"Derivatives of energy" );

}

}
Expand Down
56 changes: 30 additions & 26 deletions src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,12 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies )

subRegion.registerField< fields::flow::mass >( getName() );
subRegion.registerField< fields::flow::mass_n >( getName() );
subRegion.registerField< fields::flow::dMass >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell );;
subRegion.registerField< fields::flow::dMass >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell );

if( m_isThermal )
{
subRegion.registerField< fields::flow::dEnergy >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell );
}
} );

elemManager.forElementSubRegions< SurfaceElementSubRegion >( regionNames,
Expand Down Expand Up @@ -256,9 +261,11 @@ void SinglePhaseBase::updateMass( ElementSubRegionBase & subRegion ) const
{
GEOS_MARK_FUNCTION;

using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; // TODO check if this is correct

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 >();
arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const dMass = subRegion.getField< fields::flow::dMass >();

CoupledSolidBase const & porousSolid =
getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
Expand All @@ -267,19 +274,19 @@ void SinglePhaseBase::updateMass( ElementSubRegionBase & subRegion ) const
arrayView2d< real64 const > const porosity_n = porousSolid.getPorosity_n();

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

SingleFluidBase & fluid =
getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
arrayView2d< real64 const, singlefluid::USD_FLUID > const density = fluid.density();
arrayView2d< real64 const, singlefluid::USD_FLUID > const density_n = fluid.density_n();
arrayView2d< real64 const > const dDensity_dP = fluid.dDensity_dPressure();
arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const dDensity = fluid.dDensity();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
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;
dMass[ei][DerivOffset::dP] = ( dPorosity_dP[ei][0] * density[ei][0] + porosity[ei][0] * dDensity[ei][0][DerivOffset::dP] ) * 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
Expand All @@ -288,13 +295,11 @@ void SinglePhaseBase::updateMass( ElementSubRegionBase & subRegion ) const

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] + porosity[ei][0] * dDensity_dT[ei][0] ) * vol;
dMass[ei][DerivOffset::dT] = ( dPorosity_dT[ei][0] * density[ei][0] + porosity[ei][0] * dDensity[ei][0][DerivOffset::dT] ) * vol;
} );
}
}
Expand All @@ -303,12 +308,13 @@ void SinglePhaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const
{
GEOS_MARK_FUNCTION;

using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;

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 >();
arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const dEnergy = subRegion.getField< fields::flow::dEnergy >();

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

CoupledSolidBase const & porousSolid =
getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
Expand All @@ -322,28 +328,26 @@ void SinglePhaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const
getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const density = fluid.density();
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const fluidInternalEnergy = fluid.internalEnergy();
arrayView2d< real64 const > const dDensity_dP = fluid.dDensity_dPressure();
arrayView2d< real64 const > const dDensity_dT = fluid.dDensity_dTemperature();
arrayView2d< real64 const > const dFluidInternalEnergy_dP = fluid.dInternalEnergy_dPressure();
arrayView2d< real64 const > const dFluidInternalEnergy_dT = fluid.dInternalEnergy_dTemperature();
arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const dDensity = fluid.dDensity();
arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const dFluidInternalEnergy = fluid.dInternalEnergy();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
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] );
dEnergy[ei][DerivOffset::dP] = vol *
( dPorosity_dP[ei][0] * density[ei][0] * fluidInternalEnergy[ei][0] +
porosity[ei][0] * dDensity[ei][0][DerivOffset::dP] * fluidInternalEnergy[ei][0] +
porosity[ei][0] * density[ei][0] * dFluidInternalEnergy[ei][0][DerivOffset::dP] -
dPorosity_dP[ei][0] * rockInternalEnergy[ei][0] );
dEnergy[ei][DerivOffset::dT] = vol *
( dPorosity_dT[ei][0] * density[ei][0] * fluidInternalEnergy[ei][0] +
porosity[ei][0] * dDensity[ei][0][DerivOffset::dT] * fluidInternalEnergy[ei][0] +
porosity[ei][0] * density[ei][0] * dFluidInternalEnergy[ei][0][DerivOffset::dT] -
dPorosity_dT[ei][0] * rockInternalEnergy[ei][0] +
( 1.0 - porosity[ei][0] ) * dRockInternalEnergy_dT[ei][0] );
} );
}

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,14 @@ DECLARE_FIELD( dMobility,
NO_WRITE,
"dMobility" );

DECLARE_FIELD( dEnergy,
"dEnergy",
array2dLayoutFluid,
0,
NOPLOT,
NO_WRITE,
"Derivatives of energy" );

}

}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -151,7 +151,7 @@ class AccumulationKernel
stack.localResidual[0] = m_mass[ei] - m_mass_n[ei];

// Derivative of residual wrt to pressure in the cell
stack.localJacobian[0][0] = m_dMass[ei][0][DerivOffset::dP];
stack.localJacobian[0][0] = m_dMass[ei][DerivOffset::dP];

// Customize the kernel with this lambda
kernelOp();
Expand Down Expand Up @@ -218,7 +218,7 @@ class AccumulationKernel
/// View on mass
arrayView1d< real64 const > const m_mass;
arrayView1d< real64 const > const m_mass_n;
arrayView2d< real64 const > const m_dMass; // TODO check...
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_dMass;

/// View on the local CRS matrix
CRSMatrixView< real64, globalIndex const > const m_localMatrix;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
using Base::m_elemGhostRank;
using Base::m_localMatrix;
using Base::m_localRhs;
using Base::m_dMass;

/// Note: Derivative lineup only supports dP & dT, not component terms
static constexpr integer isThermal = NUM_DOF-1;
Expand All @@ -70,8 +71,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
: Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ),
m_energy( subRegion.template getField< fields::flow::energy >() ),
m_energy_n( subRegion.template getField< fields::flow::energy_n >() ),
m_dEnergy( subRegion.template getField< fields::flow::dEnergy >() ),
m_dMass( subRegion.template getField< fields::flow::dMass >() )
m_dEnergy( subRegion.template getField< fields::flow::dEnergy >() )
{}

/**
Expand Down Expand Up @@ -128,10 +128,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
/// View on energy
arrayView1d< real64 const > const m_energy;
arrayView1d< real64 const > const m_energy_n;
arrayView2d< real64 const > const m_dEnergy;

/// View on mass derivative
arrayView2d< real64 const > const m_dMass;
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_dEnergy;

};

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,6 @@ class SinglePhasePoromechanicsEFEM :
using Base::m_elemsToNodes;
using Base::m_constitutiveUpdate;
using Base::m_finiteElementSpace;
using Base::m_dt;


SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
Expand Down Expand Up @@ -250,24 +249,18 @@ class SinglePhasePoromechanicsEFEM :

/// The rank global fluid mass
arrayView1d< real64 const > const m_fluidMass;
arrayView1d< real64 const > const m_dFluidMass_dPressure;
arrayView1d< real64 const > const m_fluidMass_n;
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_dFluidMass;

/// The rank global densities
arrayView2d< real64 const > const m_solidDensity;
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_fluidDensity;
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_fluidDensity_n;
arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const m_dFluidDensity;

/// The rank-global fluid pressure array.
arrayView1d< real64 const > const m_matrixPressure;

/// The rank-global fluid pressure array.
arrayView1d< real64 const > const m_fracturePressure;

/// The rank-global delta-fluid pressure array.
arrayView2d< real64 const > const m_porosity_n;

arrayView2d< real64 const > const m_tractionVec;

arrayView3d< real64 const > const m_dTraction_dJump;
Expand All @@ -286,10 +279,6 @@ class SinglePhasePoromechanicsEFEM :

arrayView1d< real64 const > const m_elementVolumeCell;

arrayView1d< real64 const > const m_elementVolumeFrac;

arrayView1d< real64 const > const m_deltaVolume;

SortedArrayView< localIndex const > const m_fracturedElems;

ArrayOfArraysView< localIndex const > const m_cellsToEmbeddedSurfaces;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -77,14 +77,10 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
m_wDofNumber( jumpDofNumber ),
m_fluidMass( embeddedSurfSubRegion.template getField< fields::flow::mass >() ),
m_fluidMass_n( embeddedSurfSubRegion.template getField< fields::flow::mass_n >() ),
m_dFluidMass( embeddedSurfSubRegion.template getField< fields::flow::dMass_dPressure >() ),
m_dFluidMass( embeddedSurfSubRegion.template getField< fields::flow::dMass >() ),
m_fluidDensity( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ),
m_fluidDensity_n( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density_n() ),
m_dFluidDensity( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
fluidModelKey ) ).dDensity() ),
m_matrixPressure( elementSubRegion.template getField< fields::flow::pressure >() ),
m_fracturePressure( embeddedSurfSubRegion.template getField< fields::flow::pressure >() ),
m_porosity_n( inputConstitutiveType.getPorosity_n() ),
m_tractionVec( embeddedSurfSubRegion.getField< fields::contact::traction >() ),
m_dTraction_dJump( embeddedSurfSubRegion.getField< fields::contact::dTraction_dJump >() ),
m_dTraction_dPressure( embeddedSurfSubRegion.getField< fields::contact::dTraction_dPressure >() ),
Expand All @@ -94,8 +90,6 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
m_surfaceCenter( embeddedSurfSubRegion.getElementCenter() ),
m_surfaceArea( embeddedSurfSubRegion.getElementArea() ),
m_elementVolumeCell( elementSubRegion.getElementVolume() ),
m_elementVolumeFrac( embeddedSurfSubRegion.getElementVolume() ),
m_deltaVolume( embeddedSurfSubRegion.template getField< fields::flow::deltaVolume >() ),
m_fracturedElems( elementSubRegion.fracturedElementsList() ),
m_cellsToEmbeddedSurfaces( elementSubRegion.embeddedSurfacesList().toViewConst() ),
m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] },
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -58,13 +58,10 @@ class ThermalSinglePhasePoromechanicsEFEM :
using Base::m_fracturePresDofNumber;
using Base::m_matrixPresDofNumber;
using Base::m_wDofNumber;
using Base::m_dFluidDensity;
using Base::m_porosity_n;
using Base::m_dFluidMass;
using Base::m_fluidDensity;
using Base::m_surfaceArea;
using Base::m_elementVolumeFrac;
using Base::m_deltaVolume;
using Base::m_cellsToEmbeddedSurfaces;
using Base::m_dt;

ThermalSinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
EdgeManager const & edgeManager,
Expand Down Expand Up @@ -158,25 +155,15 @@ class ThermalSinglePhasePoromechanicsEFEM :

private:

/// Views on fluid density derivatives
arrayView2d< real64 const > const m_dFluidMass;

/// Views on fluid density derivatives
arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const m_dFluidDensity;

/// Views on fluid internal energy
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_fluidInternalEnergy;

/// Views on energy
arrayView1d< real64 const > const m_energy;
arrayView1d< real64 const > const m_energy_n;
arrayView2d< real64 const > const m_dEnergy;
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_dEnergy;

/// Views on temperature
arrayView1d< real64 const > const m_temperature;
arrayView1d< real64 const > const m_temperature_n;
/// Views on fluid internal energy
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_fluidInternalEnergy;

/// The rank-global fluid pressure array.
/// The rank-global fluid temperature array.
arrayView1d< real64 const > const m_matrixTemperature;
};

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -66,13 +66,10 @@ ThermalSinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
inputDt,
inputGravityVector,
fluidModelKey ),
m_dFluidMass_dTemperature( embeddedSurfSubRegion.template getField< fields::flow::dMass_dTemperature >() ),
m_fluidInternalEnergy( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy() ),
m_energy( elementSubRegion.template getField< fields::flow::energy >() ),
m_energy_n( elementSubRegion.template getField< fields::flow::energy_n >() ),
m_dEnergy( elementSubRegion.template getField< fields::flow::dEnergy >() ),
m_temperature( embeddedSurfSubRegion.template getField< fields::flow::temperature >() ),
m_temperature_n( embeddedSurfSubRegion.template getField< fields::flow::temperature_n >() ),
m_fluidInternalEnergy( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy() ),
m_matrixTemperature( elementSubRegion.template getField< fields::flow::temperature >() )
{}

Expand Down Expand Up @@ -171,16 +168,16 @@ complete( localIndex const k,

localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0];

stack.dFluidMassIncrement_dTemperature = m_dFluidMass[ embSurfIndex ][ DerivOffset::dT ];
stack.dFluidMassIncrement_dTemperature = m_dFluidMass[embSurfIndex][DerivOffset::dT];

// Energy balance accumulation
stack.energyIncrement = m_energy[embSurfIndex] - m_energy_n[embSurfIndex];
stack.dEnergyIncrement_dJump = m_fluidDensity( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) * m_surfaceArea[ embSurfIndex ];
stack.dEnergyIncrement_dPressure = m_dEnergy[ embSurfIndex ][ DerivOffset::dP ];
stack.dEnergyIncrement_dTemperature = m_dEnergy[ embSurfIndex ][ DerivOffset::dT ];
stack.dEnergyIncrement_dJump = m_fluidDensity[embSurfIndex][0] * m_fluidInternalEnergy[embSurfIndex][0] * m_surfaceArea[embSurfIndex];
stack.dEnergyIncrement_dPressure = m_dEnergy[embSurfIndex][DerivOffset::dP];
stack.dEnergyIncrement_dTemperature = m_dEnergy[embSurfIndex][DerivOffset::dT];

globalIndex const fracturePressureDof = m_fracturePresDofNumber[ embSurfIndex ];
globalIndex const fractureTemperatureDof = m_fracturePresDofNumber[ embSurfIndex ] + 1;
globalIndex const fracturePressureDof = m_fracturePresDofNumber[embSurfIndex];
globalIndex const fractureTemperatureDof = m_fracturePresDofNumber[embSurfIndex] + 1;
localIndex const massBalanceEquationIndex = fracturePressureDof - m_dofRankOffset;
localIndex const energyBalanceEquationIndex = massBalanceEquationIndex + 1;

Expand Down

0 comments on commit 6be1912

Please sign in to comment.