diff --git a/src/coreComponents/constitutive/solid/CompressibleSolid.hpp b/src/coreComponents/constitutive/solid/CompressibleSolid.hpp index 83e0690e9a..be4ec85d6c 100644 --- a/src/coreComponents/constitutive/solid/CompressibleSolid.hpp +++ b/src/coreComponents/constitutive/solid/CompressibleSolid.hpp @@ -55,11 +55,7 @@ class CompressibleSolidUpdates : public CoupledSolidUpdates< NullModel, PORO_TYP virtual void updateStateFromPressureAndTemperature( localIndex const k, localIndex const q, real64 const & pressure, - real64 const & GEOS_UNUSED_PARAM( pressure_k ), - real64 const & GEOS_UNUSED_PARAM( pressure_n ), - real64 const & temperature, - real64 const & GEOS_UNUSED_PARAM( temperature_k ), - real64 const & GEOS_UNUSED_PARAM( temperature_n ) ) const override final + real64 const & temperature ) const override final { m_porosityUpdate.updateFromPressureAndTemperature( k, q, pressure, temperature ); real64 const porosity = m_porosityUpdate.getPorosity( k, q ); diff --git a/src/coreComponents/constitutive/solid/CoupledSolid.hpp b/src/coreComponents/constitutive/solid/CoupledSolid.hpp index 2c2e508e12..cb7b7da5b7 100644 --- a/src/coreComponents/constitutive/solid/CoupledSolid.hpp +++ b/src/coreComponents/constitutive/solid/CoupledSolid.hpp @@ -85,11 +85,20 @@ class CoupledSolidUpdates virtual void updateStateFromPressureAndTemperature( localIndex const k, localIndex const q, real64 const & pressure, - real64 const & pressure_k, - real64 const & pressure_n, - real64 const & temperature, - real64 const & temperature_k, - real64 const & temperature_n ) const + real64 const & temperature ) const + { + GEOS_UNUSED_VAR( k, q, pressure, temperature ); + } + + GEOS_HOST_DEVICE + virtual void updateStateFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n ) const { GEOS_UNUSED_VAR( k, q, pressure, pressure_k, pressure_n, diff --git a/src/coreComponents/constitutive/solid/PorousSolid.hpp b/src/coreComponents/constitutive/solid/PorousSolid.hpp index b4a9429eb4..7d301222be 100644 --- a/src/coreComponents/constitutive/solid/PorousSolid.hpp +++ b/src/coreComponents/constitutive/solid/PorousSolid.hpp @@ -54,14 +54,14 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, {} GEOS_HOST_DEVICE - virtual void updateStateFromPressureAndTemperature( localIndex const k, - localIndex const q, - real64 const & pressure, - real64 const & pressure_k, - real64 const & pressure_n, - real64 const & temperature, - real64 const & temperature_k, - real64 const & temperature_n ) const override final + virtual void updateStateFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n ) const override final { updateBiotCoefficientAndAssignModuli( k ); diff --git a/src/coreComponents/constitutive/solid/porosity/BiotPorosity.hpp b/src/coreComponents/constitutive/solid/porosity/BiotPorosity.hpp index 8ea2df08b2..54730c0b73 100644 --- a/src/coreComponents/constitutive/solid/porosity/BiotPorosity.hpp +++ b/src/coreComponents/constitutive/solid/porosity/BiotPorosity.hpp @@ -105,7 +105,9 @@ class BiotPorosityUpdates : public PorosityBaseUpdates dPorosity_dPressure = biotSkeletonModulusInverse; dPorosity_dTemperature = -porosityThermalExpansion; - savePorosity( k, q, porosity, dPorosity_dPressure, dPorosity_dTemperature ); + m_newPorosity[k][q] = porosity; + m_dPorosity_dPressure[k][q] = dPorosity_dPressure; + m_dPorosity_dTemperature[k][q] = dPorosity_dTemperature; } GEOS_HOST_DEVICE diff --git a/src/coreComponents/constitutive/solid/porosity/PorosityBase.hpp b/src/coreComponents/constitutive/solid/porosity/PorosityBase.hpp index 72cffcf623..07f97ac5b9 100644 --- a/src/coreComponents/constitutive/solid/porosity/PorosityBase.hpp +++ b/src/coreComponents/constitutive/solid/porosity/PorosityBase.hpp @@ -59,30 +59,6 @@ class PorosityBaseUpdates m_referencePorosity ( referencePorosity ) {} - /** - * @brief Helper to save porosity back to m_newPorosity array - * - * This is mostly defined for improving code readability. - * - * @param[in] k Element index. - * @param[in] q Quadrature point index. - * @param[in] porosity porosity to be saved to m_newPorosity[k][q] - * @param[in] dPorosity_dPressure porosity derivative w.r.t pressure to be saved to m_dPorosity_dPressure[k][q] - * @param[in] dPorosity_dTemperature porosity derivative w.r.t temperature to be saved to m_dPorosity_dTemperature[k][q] - */ - GEOS_HOST_DEVICE - GEOS_FORCE_INLINE - void savePorosity( localIndex const k, - localIndex const q, - real64 const & porosity, - real64 const & dPorosity_dPressure, - real64 const & dPorosity_dTemperature ) const - { - m_newPorosity[k][q] = porosity; - m_dPorosity_dPressure[k][q] = dPorosity_dPressure; - m_dPorosity_dTemperature[k][q] = dPorosity_dTemperature; - } - GEOS_HOST_DEVICE inline real64 getPorosity( localIndex const k, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index c081f20d4c..7a33376faa 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -45,11 +45,7 @@ template< typename POROUSWRAPPER_TYPE > void updatePorosityAndPermeabilityFromPressureAndTemperature( POROUSWRAPPER_TYPE porousWrapper, CellElementSubRegion & subRegion, arrayView1d< real64 const > const & pressure, - arrayView1d< real64 const > const & pressure_k, - arrayView1d< real64 const > const & pressure_n, - arrayView1d< real64 const > const & temperature, - arrayView1d< real64 const > const & temperature_k, - arrayView1d< real64 const > const & temperature_n ) + arrayView1d< real64 const > const & temperature ) { forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const k ) { @@ -57,15 +53,35 @@ void updatePorosityAndPermeabilityFromPressureAndTemperature( POROUSWRAPPER_TYPE { porousWrapper.updateStateFromPressureAndTemperature( k, q, pressure[k], - pressure_k[k], - pressure_n[k], - temperature[k], - temperature_k[k], - temperature_n[k] ); + temperature[k] ); } } ); } +template< typename POROUSWRAPPER_TYPE > +void updatePorosityAndPermeabilityFixedStress( POROUSWRAPPER_TYPE porousWrapper, + CellElementSubRegion & subRegion, + arrayView1d< real64 const > const & pressure, + arrayView1d< real64 const > const & pressure_k, + arrayView1d< real64 const > const & pressure_n, + arrayView1d< real64 const > const & temperature, + arrayView1d< real64 const > const & temperature_k, + arrayView1d< real64 const > const & temperature_n ) +{ + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const k ) + { + for( localIndex q = 0; q < porousWrapper.numGauss(); ++q ) + { + porousWrapper.updateStateFixedStress( k, q, + pressure[k], + pressure_k[k], + pressure_n[k], + temperature[k], + temperature_k[k], + temperature_n[k] ); + } + } ); +} template< typename POROUSWRAPPER_TYPE > void updatePorosityAndPermeabilityFromPressureAndAperture( POROUSWRAPPER_TYPE porousWrapper, @@ -168,6 +184,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() ); } } ); @@ -583,10 +601,7 @@ void FlowSolverBase::updatePorosityAndPermeability( CellElementSubRegion & subRe GEOS_MARK_FUNCTION; arrayView1d< real64 const > const & pressure = subRegion.getField< fields::flow::pressure >(); - arrayView1d< real64 const > const & pressure_n = subRegion.getField< fields::flow::pressure_n >(); - arrayView1d< real64 const > const & temperature = subRegion.getField< fields::flow::temperature >(); - arrayView1d< real64 const > const & temperature_n = subRegion.getField< fields::flow::temperature_n >(); string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() ); CoupledSolidBase & porousSolid = subRegion.template getConstitutiveModel< CoupledSolidBase >( solidName ); @@ -596,13 +611,15 @@ void FlowSolverBase::updatePorosityAndPermeability( CellElementSubRegion & subRe typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousWrapper = castedPorousSolid.createKernelUpdates(); if( m_isFixedStressPoromechanicsUpdate ) { + arrayView1d< real64 const > const & pressure_n = subRegion.getField< fields::flow::pressure_n >(); arrayView1d< real64 const > const & pressure_k = subRegion.getField< fields::flow::pressure_k >(); + arrayView1d< real64 const > const & temperature_n = subRegion.getField< fields::flow::temperature_n >(); arrayView1d< real64 const > const & temperature_k = subRegion.getField< fields::flow::temperature_k >(); - updatePorosityAndPermeabilityFromPressureAndTemperature( porousWrapper, subRegion, pressure, pressure_k, pressure_n, temperature, temperature_k, temperature_n ); + updatePorosityAndPermeabilityFixedStress( porousWrapper, subRegion, pressure, pressure_k, pressure_n, temperature, temperature_k, temperature_n ); } else { - updatePorosityAndPermeabilityFromPressureAndTemperature( porousWrapper, subRegion, pressure, pressure_n, pressure_n, temperature, temperature_n, temperature_n ); + updatePorosityAndPermeabilityFromPressureAndTemperature( porousWrapper, subRegion, pressure, temperature ); } } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp index ddb20e7bf5..7c3306a8b3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp @@ -241,46 +241,38 @@ DECLARE_FIELD( temperatureScalingFactor, NO_WRITE, "Scaling factors for temperature" ); -DECLARE_FIELD( mass, - "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." ); - } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 98934988d8..1f95b1c247 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -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; @@ -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() ); } } ); @@ -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() ); } } } ); @@ -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(); @@ -285,14 +286,31 @@ 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] + porosity[ei][0] * dDensity_dT[ei][0] ) * vol; + } ); + } } void SinglePhaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const @@ -300,24 +318,46 @@ 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] ); } ); } @@ -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 ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index db72beaa75..9e8725a381 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -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 ) @@ -454,8 +447,6 @@ void SinglePhaseBase::accumulationAssemblyLaunch( DofManager const & dofManager, createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), dofKey, subRegion, - fluid, - solid, localMatrix, localRhs ); } @@ -466,8 +457,6 @@ void SinglePhaseBase::accumulationAssemblyLaunch( DofManager const & dofManager, createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), dofKey, subRegion, - fluid, - solid, localMatrix, localRhs ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp index b68848e634..81bcc6b997 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp @@ -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 >, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp index c2e2845ff7..8df69de36f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp @@ -25,6 +25,7 @@ #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" #include "codingUtilities/Utilities.hpp" namespace geos @@ -56,28 +57,20 @@ class AccumulationKernel * @param[in] rankOffset the offset of my MPI rank * @param[in] dofKey the string key to retrieve the degress of freedom numbers * @param[in] subRegion the element subregion - * @param[in] fluid the fluid model - * @param[in] solid the solid model * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector */ AccumulationKernel( globalIndex const rankOffset, string const dofKey, SUBREGION_TYPE const & subRegion, - constitutive::SingleFluidBase const & fluid, - constitutive::CoupledSolidBase const & solid, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) : m_rankOffset( rankOffset ), m_dofNumber( subRegion.template getReference< array1d< globalIndex > >( dofKey ) ), m_elemGhostRank( subRegion.ghostRank() ), - m_volume( subRegion.getElementVolume() ), - m_deltaVolume( subRegion.template getField< fields::flow::deltaVolume >() ), - m_porosity( solid.getPorosity() ), - m_dPoro_dPres( solid.getDporosity_dPressure() ), - m_density( fluid.density() ), - m_dDensity_dPres( fluid.dDensity_dPressure() ), + m_mass( subRegion.template getField< fields::flow::mass >() ), + m_dMass_dPres( subRegion.template getField< fields::flow::dMass_dPressure >() ), m_mass_n( subRegion.template getField< fields::flow::mass_n >() ), m_localMatrix( localMatrix ), m_localRhs( localRhs ) @@ -91,14 +84,6 @@ class AccumulationKernel { public: - // Pore volume information - - /// Pore volume at time n+1 - real64 poreVolume = 0.0; - - /// Derivative of pore volume with respect to pressure - real64 dPoreVolume_dPres = 0.0; - // Residual information /// Index of the local row corresponding to this element @@ -134,10 +119,6 @@ class AccumulationKernel void setup( localIndex const ei, StackVariables & stack ) const { - // initialize the pore volume - stack.poreVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * m_porosity[ei][0]; - stack.dPoreVolume_dPres = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dPres[ei][0]; - // set row index and degrees of freedom indices for this element stack.localRow = m_dofNumber[ei] - m_rankOffset; for( integer idof = 0; idof < numDof; ++idof ) @@ -160,10 +141,10 @@ class AccumulationKernel FUNC && kernelOp = NoOpFunc{} ) const { // Residual contribution is mass conservation in the cell - stack.localResidual[0] = stack.poreVolume * m_density[ei][0] - m_mass_n[ei]; + stack.localResidual[0] = m_mass[ei] - m_mass_n[ei]; // Derivative of residual wrt to pressure in the cell - stack.localJacobian[0][0] = stack.dPoreVolume_dPres * m_density[ei][0] + m_dDensity_dPres[ei][0] * stack.poreVolume; + stack.localJacobian[0][0] = m_dMass_dPres[ei]; // Customize the kernel with this lambda kernelOp(); @@ -227,19 +208,9 @@ class AccumulationKernel /// View on the ghost ranks arrayView1d< integer const > const m_elemGhostRank; - /// View on the element volumes - arrayView1d< real64 const > const m_volume; - arrayView1d< real64 const > const m_deltaVolume; - - /// Views on the porosity - arrayView2d< real64 const > const m_porosity; - arrayView2d< real64 const > const m_dPoro_dPres; - - /// Views on density - arrayView2d< real64 const > const m_density; - arrayView2d< real64 const > const m_dDensity_dPres; - /// View on mass + arrayView1d< real64 const > const m_mass; + arrayView1d< real64 const > const m_dMass_dPres; arrayView1d< real64 const > const m_mass_n; /// View on the local CRS matrix @@ -273,11 +244,9 @@ class SurfaceElementAccumulationKernel : public AccumulationKernel< SurfaceEleme SurfaceElementAccumulationKernel( globalIndex const rankOffset, string const dofKey, SurfaceElementSubRegion const & subRegion, - constitutive::SingleFluidBase const & fluid, - constitutive::CoupledSolidBase const & solid, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) - : Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ) + : Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ) , m_creationMass( subRegion.getField< fields::flow::massCreated >() ) {} @@ -291,13 +260,11 @@ class SurfaceElementAccumulationKernel : public AccumulationKernel< SurfaceEleme void computeAccumulation( localIndex const ei, Base::StackVariables & stack ) const { - Base::computeAccumulation( ei, stack, [&] () + Base::computeAccumulation( ei, stack ); + if( Base::m_mass_n[ei] > 1.1 * m_creationMass[ei] ) { - if( Base::m_mass_n[ei] > 1.1 * m_creationMass[ei] ) - { - stack.localResidual[0] += m_creationMass[ei] * 0.25; - } - } ); + stack.localResidual[0] += m_creationMass[ei] * 0.25; + } } protected: @@ -329,25 +296,23 @@ class AccumulationKernelFactory createAndLaunch( globalIndex const rankOffset, string const dofKey, SUBREGION_TYPE const & subRegion, - constitutive::SingleFluidBase const & fluid, - constitutive::CoupledSolidBase const & solid, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { if constexpr ( std::is_base_of_v< CellElementSubRegion, SUBREGION_TYPE > ) { integer constexpr NUM_DOF = 1; - AccumulationKernel< CellElementSubRegion, NUM_DOF > kernel( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + AccumulationKernel< CellElementSubRegion, NUM_DOF > kernel( rankOffset, dofKey, subRegion, localMatrix, localRhs ); AccumulationKernel< CellElementSubRegion, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel ); } else if constexpr ( std::is_base_of_v< SurfaceElementSubRegion, SUBREGION_TYPE > ) { - SurfaceElementAccumulationKernel kernel( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + SurfaceElementAccumulationKernel kernel( rankOffset, dofKey, subRegion, localMatrix, localRhs ); SurfaceElementAccumulationKernel::launch< POLICY >( subRegion.size(), kernel ); } else { - GEOS_UNUSED_VAR( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + GEOS_UNUSED_VAR( rankOffset, dofKey, subRegion, localMatrix, localRhs ); GEOS_ERROR( "Unsupported subregion type: " << typeid(SUBREGION_TYPE).name() ); } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp index 2cba35dac4..3e20187ad4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp @@ -46,12 +46,6 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU using Base::m_rankOffset; using Base::m_dofNumber; using Base::m_elemGhostRank; - using Base::m_volume; - using Base::m_deltaVolume; - using Base::m_porosity; - using Base::m_dPoro_dPres; - using Base::m_density; - using Base::m_dDensity_dPres; using Base::m_localMatrix; using Base::m_localRhs; @@ -68,19 +62,14 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU AccumulationKernel( globalIndex const rankOffset, string const dofKey, SUBREGION_TYPE const & subRegion, - constitutive::SingleFluidBase const & fluid, - constitutive::CoupledSolidBase const & solid, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) - : Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ), - m_dDensity_dTemp( fluid.dDensity_dTemperature() ), - m_dPoro_dTemp( solid.getDporosity_dTemperature() ), - m_internalEnergy( fluid.internalEnergy() ), - m_dInternalEnergy_dPres( fluid.dInternalEnergy_dPressure() ), - m_dInternalEnergy_dTemp( fluid.dInternalEnergy_dTemperature() ), - m_rockInternalEnergy( solid.getInternalEnergy() ), - m_dRockInternalEnergy_dTemp( solid.getDinternalEnergy_dTemperature() ), - m_energy_n( subRegion.template getField< fields::flow::energy_n >() ) + : Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ), + m_energy( subRegion.template getField< fields::flow::energy >() ), + m_dEnergy_dPres( subRegion.template getField< fields::flow::dEnergy_dPressure >() ), + m_dEnergy_dTemp( subRegion.template getField< fields::flow::dEnergy_dTemperature >() ), + m_energy_n( subRegion.template getField< fields::flow::energy_n >() ), + m_dMass_dTemp( subRegion.template getField< fields::flow::dMass_dTemperature >() ) {} /** @@ -88,60 +77,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack */ struct StackVariables : public Base::StackVariables - { -public: - - GEOS_HOST_DEVICE - StackVariables() - : Base::StackVariables() - {} - - using Base::StackVariables::poreVolume; - using Base::StackVariables::dPoreVolume_dPres; - using Base::StackVariables::localRow; - using Base::StackVariables::dofIndices; - using Base::StackVariables::localResidual; - using Base::StackVariables::localJacobian; - - /// Derivative of pore volume with respect to temperature - real64 dPoreVolume_dTemp = 0.0; - - // Solid energy - - /// Solid energy at time n+1 - real64 solidEnergy = 0.0; - - /// Derivative of solid internal energy with respect to pressure - real64 dSolidEnergy_dPres = 0.0; - - /// Derivative of solid internal energy with respect to temperature - real64 dSolidEnergy_dTemp = 0.0; - }; - - - /** - * @brief Performs the setup phase for the kernel. - * @param[in] ei the element index - * @param[in] stack the stack variables - */ - GEOS_HOST_DEVICE - void setup( localIndex const ei, - StackVariables & stack ) const - { - Base::setup( ei, stack ); - - stack.dPoreVolume_dTemp = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dTemp[ei][0]; - - // initialize the solid volume - real64 const solidVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * ( 1.0 - m_porosity[ei][0] ); - real64 const dSolidVolume_dPres = -( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dPres[ei][0]; - real64 const dSolidVolume_dTemp = -( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dTemp[ei][0]; - - // initialize the solid internal energy - stack.solidEnergy = solidVolume * m_rockInternalEnergy[ei][0]; - stack.dSolidEnergy_dPres = dSolidVolume_dPres * m_rockInternalEnergy[ei][0]; - stack.dSolidEnergy_dTemp = solidVolume * m_dRockInternalEnergy_dTemp[ei][0] + dSolidVolume_dTemp * m_rockInternalEnergy[ei][0]; - } + {}; /** * @brief Compute the local accumulation contributions to the residual and Jacobian @@ -154,36 +90,15 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU void computeAccumulation( localIndex const ei, StackVariables & stack ) const { - stack.localResidual[numEqn-1] = -m_energy_n[ei]; - - Base::computeAccumulation( ei, stack, [&] () - { - // Step 1: assemble the derivatives of the mass balance equation w.r.t temperature - stack.localJacobian[0][numDof-1] = stack.poreVolume * m_dDensity_dTemp[ei][0] + stack.dPoreVolume_dTemp * m_density[ei][0]; - - // Step 2: assemble the fluid part of the accumulation term of the energy equation - real64 const fluidEnergy = stack.poreVolume * m_density[ei][0] * m_internalEnergy[ei][0]; - - real64 const dFluidEnergy_dP = stack.dPoreVolume_dPres * m_density[ei][0] * m_internalEnergy[ei][0] - + stack.poreVolume * m_dDensity_dPres[ei][0] * m_internalEnergy[ei][0] - + stack.poreVolume * m_density[ei][0] * m_dInternalEnergy_dPres[ei][0]; - - real64 const dFluidEnergy_dT = stack.poreVolume * m_dDensity_dTemp[ei][0] * m_internalEnergy[ei][0] - + stack.poreVolume * m_density[ei][0] * m_dInternalEnergy_dTemp[ei][0] - + stack.dPoreVolume_dTemp * m_density[ei][0] * m_internalEnergy[ei][0]; - - // local accumulation - stack.localResidual[numEqn-1] += fluidEnergy; + Base::computeAccumulation( ei, stack ); - // derivatives w.r.t. pressure and temperature - stack.localJacobian[numEqn-1][0] = dFluidEnergy_dP; - stack.localJacobian[numEqn-1][numDof-1] = dFluidEnergy_dT; - } ); + // assemble the derivatives of the mass balance equation w.r.t temperature + stack.localJacobian[0][numDof-1] = m_dMass_dTemp[ei]; - // Step 3: assemble the solid part of the accumulation term of the energy equation - stack.localResidual[numEqn-1] += stack.solidEnergy; - stack.localJacobian[numEqn-1][0] += stack.dSolidEnergy_dPres; - stack.localJacobian[numEqn-1][numDof-1] += stack.dSolidEnergy_dTemp; + // assemble the accumulation term of the energy equation + stack.localResidual[numEqn-1] = m_energy[ei] - m_energy_n[ei]; + stack.localJacobian[numEqn-1][0] += m_dEnergy_dPres[ei]; + stack.localJacobian[numEqn-1][numDof-1] += m_dEnergy_dTemp[ei]; } /** @@ -208,24 +123,15 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU protected: - /// View on derivative of fluid density w.r.t temperature - arrayView2d< real64 const > const m_dDensity_dTemp; - - /// View on derivative of porosity w.r.t temperature - arrayView2d< real64 const > const m_dPoro_dTemp; - - /// Views on fluid internal energy - arrayView2d< real64 const > const m_internalEnergy; - arrayView2d< real64 const > const m_dInternalEnergy_dPres; - arrayView2d< real64 const > const m_dInternalEnergy_dTemp; - - /// Views on rock internal energy - arrayView2d< real64 const > const m_rockInternalEnergy; - arrayView2d< real64 const > const m_dRockInternalEnergy_dTemp; - /// View on energy + arrayView1d< real64 const > const m_energy; + arrayView1d< real64 const > const m_dEnergy_dPres; + arrayView1d< real64 const > const m_dEnergy_dTemp; arrayView1d< real64 const > const m_energy_n; + /// View on mass derivative + arrayView1d< real64 const > const m_dMass_dTemp; + }; /** @@ -252,11 +158,9 @@ class SurfaceElementAccumulationKernel : public AccumulationKernel< SurfaceEleme SurfaceElementAccumulationKernel( globalIndex const rankOffset, string const dofKey, SurfaceElementSubRegion const & subRegion, - constitutive::SingleFluidBase const & fluid, - constitutive::CoupledSolidBase const & solid, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) - : Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ), + : Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ), m_creationMass( subRegion.getField< fields::flow::massCreated >() ) {} @@ -306,25 +210,23 @@ class AccumulationKernelFactory createAndLaunch( globalIndex const rankOffset, string const dofKey, SUBREGION_TYPE const & subRegion, - constitutive::SingleFluidBase const & fluid, - constitutive::CoupledSolidBase const & solid, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { if constexpr ( std::is_base_of_v< CellElementSubRegion, SUBREGION_TYPE > ) { integer constexpr NUM_DOF = 2; - AccumulationKernel< CellElementSubRegion, NUM_DOF > kernel( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + AccumulationKernel< CellElementSubRegion, NUM_DOF > kernel( rankOffset, dofKey, subRegion, localMatrix, localRhs ); AccumulationKernel< CellElementSubRegion, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel ); } else if constexpr ( std::is_base_of_v< SurfaceElementSubRegion, SUBREGION_TYPE > ) { - SurfaceElementAccumulationKernel kernel( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + SurfaceElementAccumulationKernel kernel( rankOffset, dofKey, subRegion, localMatrix, localRhs ); SurfaceElementAccumulationKernel::launch< POLICY >( subRegion.size(), kernel ); } else { - GEOS_UNUSED_VAR( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + GEOS_UNUSED_VAR( rankOffset, dofKey, subRegion, localMatrix, localRhs ); GEOS_ERROR( "Unsupported subregion type: " << typeid(SUBREGION_TYPE).name() ); } } diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp index ac266ec21d..d26c58f6c6 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp @@ -247,10 +247,13 @@ class SinglePhasePoromechanicsEFEM : arrayView1d< globalIndex const > const m_wDofNumber; + /// 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; + /// The rank global densities - arrayView2d< real64 const > const m_solidDensity; arrayView2d< real64 const > const m_fluidDensity; - arrayView2d< real64 const > const m_fluidDensity_n; arrayView2d< real64 const > const m_dFluidDensity_dPressure; /// The rank-global fluid pressure array. diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp index efc8ae09ea..578a3d898f 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp @@ -75,9 +75,10 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, m_matrixPresDofNumber( elementSubRegion.template getReference< array1d< globalIndex > >( inputFlowDofKey ) ), m_fracturePresDofNumber( embeddedSurfSubRegion.template getReference< array1d< globalIndex > >( inputFlowDofKey ) ), m_wDofNumber( jumpDofNumber ), - m_solidDensity( inputConstitutiveType.getDensity() ), + m_fluidMass( embeddedSurfSubRegion.template getField< fields::flow::mass >() ), + m_dFluidMass_dPressure( embeddedSurfSubRegion.template getField< fields::flow::dMass_dPressure >() ), + m_fluidMass_n( embeddedSurfSubRegion.template getField< fields::flow::mass_n >() ), 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_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dPressure() ), m_matrixPressure( elementSubRegion.template getField< fields::flow::pressure >() ), @@ -321,12 +322,9 @@ complete( localIndex const k, real64 const localJumpFracPressureJacobian = m_surfaceArea[embSurfIndex]; // Mass balance accumulation - real64 const newVolume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex ); - real64 const newMass = m_fluidDensity( embSurfIndex, 0 ) * newVolume; - real64 const oldMass = m_fluidDensity_n( embSurfIndex, 0 ) * m_elementVolumeFrac( embSurfIndex ); - real64 const localFlowResidual = ( newMass - oldMass ); + real64 const localFlowResidual = m_fluidMass[embSurfIndex] - m_fluidMass_n[embSurfIndex]; real64 const localFlowJumpJacobian = m_fluidDensity( embSurfIndex, 0 ) * m_surfaceArea[ embSurfIndex ]; - real64 const localFlowFlowJacobian = m_dFluidDensity_dPressure( embSurfIndex, 0 ) * newVolume; + real64 const localFlowFlowJacobian = m_dFluidMass_dPressure[ embSurfIndex ]; for( localIndex i = 0; i < nUdof; ++i ) { diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp index f996e34f0e..d0ee383c02 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp @@ -57,7 +57,6 @@ class ThermalSinglePhasePoromechanicsEFEM : using Base::m_matrixPresDofNumber; using Base::m_wDofNumber; using Base::m_fluidDensity; - using Base::m_fluidDensity_n; using Base::m_dFluidDensity_dPressure; using Base::m_porosity_n; using Base::m_surfaceArea; @@ -66,8 +65,6 @@ class ThermalSinglePhasePoromechanicsEFEM : using Base::m_cellsToEmbeddedSurfaces; using Base::m_dt; - - ThermalSinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, EdgeManager const & edgeManager, FaceManager const & faceManager, @@ -160,18 +157,24 @@ class ThermalSinglePhasePoromechanicsEFEM : private: + /// Views on fluid density derivative wrt temperature + arrayView1d< real64 const > const m_dFluidMass_dTemperature; + /// Views on fluid density derivative wrt temperature arrayView2d< real64 const > const m_dFluidDensity_dTemperature; /// Views on fluid internal energy - arrayView2d< real64 const > const m_fluidInternalEnergy_n; arrayView2d< real64 const > const m_fluidInternalEnergy; - arrayView2d< real64 const > const m_dFluidInternalEnergy_dPressure; - arrayView2d< real64 const > const m_dFluidInternalEnergy_dTemperature; + + /// Views on energy + arrayView1d< real64 const > const m_energy; + arrayView1d< real64 const > const m_dEnergy_dPressure; + arrayView1d< real64 const > const m_dEnergy_dTemperature; + arrayView1d< real64 const > const m_energy_n; /// Views on temperature - arrayView1d< real64 const > const m_temperature_n; arrayView1d< real64 const > const m_temperature; + arrayView1d< real64 const > const m_temperature_n; /// The rank-global fluid pressure array. arrayView1d< real64 const > const m_matrixTemperature; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp index c538093411..f08bc3be48 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp @@ -66,18 +66,15 @@ ThermalSinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, inputDt, inputGravityVector, fluidModelKey ), - m_dFluidDensity_dTemperature( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( - fluidModelKey ) ).dDensity_dTemperature() ), - m_fluidInternalEnergy_n( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy_n() ), + m_dFluidMass_dTemperature( embeddedSurfSubRegion.template getField< fields::flow::dMass_dTemperature >() ), m_fluidInternalEnergy( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy() ), - m_dFluidInternalEnergy_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( - fluidModelKey ) ).dInternalEnergy_dPressure() ), - m_dFluidInternalEnergy_dTemperature( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( - fluidModelKey ) ).dInternalEnergy_dTemperature() ), - m_temperature_n( embeddedSurfSubRegion.template getField< fields::flow::temperature_n >() ), + m_energy( elementSubRegion.template getField< fields::flow::energy >() ), + m_dEnergy_dPressure( elementSubRegion.template getField< fields::flow::dEnergy_dPressure >() ), + m_dEnergy_dTemperature( elementSubRegion.template getField< fields::flow::dEnergy_dTemperature >() ), + m_energy_n( elementSubRegion.template getField< fields::flow::energy_n >() ), m_temperature( embeddedSurfSubRegion.template getField< fields::flow::temperature >() ), + m_temperature_n( embeddedSurfSubRegion.template getField< fields::flow::temperature_n >() ), m_matrixTemperature( elementSubRegion.template getField< fields::flow::temperature >() ) - {} @@ -174,19 +171,14 @@ complete( localIndex const k, } localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0]; - // Energy balance accumulation - real64 const volume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex ); - real64 const volume_n = m_elementVolumeFrac( embSurfIndex ); - real64 const fluidEnergy = m_fluidDensity( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) * volume; - real64 const fluidEnergy_n = m_fluidDensity_n( embSurfIndex, 0 ) * m_fluidInternalEnergy_n( embSurfIndex, 0 ) * volume_n; - stack.dFluidMassIncrement_dTemperature = m_dFluidDensity_dTemperature( embSurfIndex, 0 ) * volume; + stack.dFluidMassIncrement_dTemperature = m_dFluidMass_dTemperature[ embSurfIndex ]; - stack.energyIncrement = fluidEnergy - fluidEnergy_n; + // 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_dFluidDensity_dPressure( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) * volume; - stack.dEnergyIncrement_dTemperature = ( m_dFluidDensity_dTemperature( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) + - m_fluidDensity( embSurfIndex, 0 ) * m_dFluidInternalEnergy_dTemperature( embSurfIndex, 0 ) ) * volume; + stack.dEnergyIncrement_dPressure = m_dEnergy_dPressure[ embSurfIndex ]; + stack.dEnergyIncrement_dTemperature = m_dEnergy_dTemperature[ embSurfIndex ]; globalIndex const fracturePressureDof = m_fracturePresDofNumber[ embSurfIndex ]; globalIndex const fractureTemperatureDof = m_fracturePresDofNumber[ embSurfIndex ] + 1; @@ -195,7 +187,6 @@ complete( localIndex const k, if( massBalanceEquationIndex >= 0 && massBalanceEquationIndex < m_matrix.numRows() ) { - m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( massBalanceEquationIndex, &fractureTemperatureDof, &stack.dFluidMassIncrement_dTemperature, @@ -204,7 +195,6 @@ complete( localIndex const k, if( energyBalanceEquationIndex >= 0 && energyBalanceEquationIndex < m_matrix.numRows() ) { - m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( energyBalanceEquationIndex, &stack.jumpColIndices[0], &stack.dEnergyIncrement_dJump,