From 8b648917d2967d3f83a7c3668fcb428d84e231ed Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Sun, 8 Dec 2024 17:13:35 -0600 Subject: [PATCH] wip --- .../fluidFlow/FlowSolverBase.cpp | 2 + .../fluidFlow/FlowSolverBaseFields.hpp | 34 ++-- .../fluidFlow/SinglePhaseBase.cpp | 72 +++++++-- .../fluidFlow/SinglePhaseBase.hpp | 11 -- .../fluidFlow/SinglePhaseBaseFields.hpp | 40 +++++ .../singlePhase/AccumulationKernels.hpp | 65 ++------ .../ThermalAccumulationKernels.hpp | 146 +++--------------- .../SinglePhasePoromechanics.hpp | 13 +- .../SinglePhasePoromechanicsEFEM.hpp | 7 +- .../SinglePhasePoromechanicsEFEM_impl.hpp | 12 +- .../SinglePhasePoromechanics_impl.hpp | 18 +-- .../ThermalSinglePhasePoromechanics.hpp | 22 ++- .../ThermalSinglePhasePoromechanicsEFEM.hpp | 17 +- ...ermalSinglePhasePoromechanicsEFEM_impl.hpp | 35 ++--- .../ThermalSinglePhasePoromechanics_impl.hpp | 58 ++----- 15 files changed, 217 insertions(+), 335 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index 8cd3ff1b942..3a1096f8fb4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -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() ); } } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp index d9f01e908b9..0709f4d2be8 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 812c342de0a..e7032a8dd22 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] + dDensity_dT[ei][0] * dDensity_dP[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 9e40ecda15f..d277955b997 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 8b3f80b95c9..10b80fcda09 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 d7bcb6068f1..2ee219432c1 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 792558b5b2e..5160784b114 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/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp index 711adb82808..451dce8508a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp @@ -183,11 +183,7 @@ class SinglePhasePoromechanics : GEOS_HOST_DEVICE void computeFluidIncrement( localIndex const k, localIndex const q, - real64 const & porosity, - real64 const & porosity_n, real64 const & dPorosity_dVolStrain, - real64 const & dPorosity_dPressure, - real64 const & dPorosity_dTemperature, StackVariables & stack ) const; /** @@ -252,12 +248,13 @@ class SinglePhasePoromechanics : protected: - /// Fluid density + /// Fluid density and derivatives arrayView2d< real64 const > const m_fluidDensity; - /// Fluid density at the previous converged time step - arrayView2d< real64 const > const m_fluidDensity_n; - /// Derivative of fluid density wrt pressure arrayView2d< real64 const > const m_dFluidDensity_dPressure; + /// Fluid mass and derivatives + arrayView1d< real64 const > const m_fluidMass; + arrayView1d< real64 const > const m_dFluidMass_dPressure; + arrayView1d< real64 const > const m_fluidMass_n; integer const m_performStressInitialization; }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp index e88ce14509f..b2b8a5d1bff 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_fluidMass_n; + arrayView1d< real64 const > const m_dFluidMass_dPressure; + /// 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 9229f1ff774..3373ba4aaef 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_fluidMass_n( embeddedSurfSubRegion.template getField< fields::flow::mass_n >() ), + m_dFluidMass_dPressure( embeddedSurfSubRegion.template getField< fields::flow::dMass_dPressure >() ), 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/SinglePhasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp index bb79e9e66b4..57939028a0a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp @@ -69,8 +69,10 @@ SinglePhasePoromechanics( NodeManager const & nodeManager, inputFlowDofKey, fluidModelKey ), m_fluidDensity( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ), - m_fluidDensity_n( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density_n() ), m_dFluidDensity_dPressure( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dPressure() ), + m_fluidMass( elementSubRegion.template getField< fields::flow::mass >() ), + m_dFluidMass_dPressure( elementSubRegion.template getField< fields::flow::dMass_dPressure >() ), + m_fluidMass_n( elementSubRegion.template getField< fields::flow::mass_n >() ), m_performStressInitialization( performStressInitialization ) {} @@ -122,11 +124,7 @@ smallStrainUpdate( localIndex const k, // Step 3: compute fluid mass increment computeFluidIncrement( k, q, - porosity, - porosity_n, dPorosity_dVolStrain, - dPorosity_dPressure, - dPorosity_dTemperature, stack ); } @@ -166,18 +164,12 @@ GEOS_FORCE_INLINE void SinglePhasePoromechanics< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >:: computeFluidIncrement( localIndex const k, localIndex const q, - real64 const & porosity, - real64 const & porosity_n, real64 const & dPorosity_dVolStrain, - real64 const & dPorosity_dPressure, - real64 const & dPorosity_dTemperature, StackVariables & stack ) const { - GEOS_UNUSED_VAR( dPorosity_dTemperature ); - - stack.fluidMassIncrement = porosity * m_fluidDensity( k, q ) - porosity_n * m_fluidDensity_n( k, q ); + stack.fluidMassIncrement = m_fluidMass[k] - m_fluidMass_n[k]; stack.dFluidMassIncrement_dVolStrainIncrement = dPorosity_dVolStrain * m_fluidDensity( k, q ); - stack.dFluidMassIncrement_dPressure = dPorosity_dPressure * m_fluidDensity( k, q ) + porosity * m_dFluidDensity_dPressure( k, q ); + stack.dFluidMassIncrement_dPressure = m_dFluidMass_dPressure[k]; } template< typename SUBREGION_TYPE, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp index 378285326bd..1e5252bafea 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp @@ -66,8 +66,6 @@ class ThermalSinglePhasePoromechanics : using Base::m_pressure; using Base::m_pressure_n; using Base::m_fluidDensity; - using Base::m_fluidDensity_n; - using Base::m_dFluidDensity_dPressure; using Base::m_solidDensity; using Base::m_flowDofNumber; using Base::m_dt; @@ -212,11 +210,7 @@ class ThermalSinglePhasePoromechanics : GEOS_HOST_DEVICE void computeFluidIncrement( localIndex const k, localIndex const q, - real64 const & porosity, - real64 const & porosity_n, real64 const & dPorosity_dVolStrain, - real64 const & dPorosity_dPressure, - real64 const & dPorosity_dTemperature, StackVariables & stack ) const; /** @@ -279,23 +273,27 @@ class ThermalSinglePhasePoromechanics : protected: + /// Views on fluid mass 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 rock internal energy - arrayView2d< real64 const > const m_rockInternalEnergy_n; arrayView2d< real64 const > const m_rockInternalEnergy; - arrayView2d< real64 const > const m_dRockInternalEnergy_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; }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp index b90d9006683..13d37dd564d 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 b838e72f7b2..82249efd27e 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,15 @@ 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; - 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; + // Energy balance accumulation + stack.energyIncrement = m_energy[embSurfIndex] - m_energy_n[embSurfIndex]; + stack.dEnergyIncrement_dJump = // TODO solid part + m_fluidDensity( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) * m_surfaceArea[ embSurfIndex ]; + 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 +188,6 @@ complete( localIndex const k, if( massBalanceEquationIndex >= 0 && massBalanceEquationIndex < m_matrix.numRows() ) { - m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( massBalanceEquationIndex, &fractureTemperatureDof, &stack.dFluidMassIncrement_dTemperature, @@ -204,7 +196,6 @@ complete( localIndex const k, if( energyBalanceEquationIndex >= 0 && energyBalanceEquationIndex < m_matrix.numRows() ) { - m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( energyBalanceEquationIndex, &stack.jumpColIndices[0], &stack.dEnergyIncrement_dJump, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp index 3f20f8f7618..f40c01bf440 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp @@ -65,19 +65,17 @@ ThermalSinglePhasePoromechanics( NodeManager const & nodeManager, inputFlowDofKey, performStressInitialization, fluidModelKey ), + m_dFluidMass_dTemperature( elementSubRegion.template getField< fields::flow::dMass_dTemperature >() ), m_dFluidDensity_dTemperature( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dTemperature() ), - m_fluidInternalEnergy_n( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy_n() ), m_fluidInternalEnergy( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy() ), - m_dFluidInternalEnergy_dPressure( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( - fluidModelKey ) ).dInternalEnergy_dPressure() ), - m_dFluidInternalEnergy_dTemperature( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( - fluidModelKey ) ).dInternalEnergy_dTemperature() ), - m_rockInternalEnergy_n( inputConstitutiveType.getInternalEnergy_n() ), m_rockInternalEnergy( inputConstitutiveType.getInternalEnergy() ), - m_dRockInternalEnergy_dTemperature( inputConstitutiveType.getDinternalEnergy_dTemperature() ), - m_temperature_n( elementSubRegion.template getField< fields::flow::temperature_n >() ), - m_temperature( elementSubRegion.template getField< fields::flow::temperature >() ) + 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( elementSubRegion.template getField< fields::flow::temperature >() ), + m_temperature_n( elementSubRegion.template getField< fields::flow::temperature_n >() ) {} template< typename SUBREGION_TYPE, @@ -143,11 +141,7 @@ smallStrainUpdate( localIndex const k, // Step 3: compute fluid mass increment computeFluidIncrement( k, q, - porosity, - porosity_n, dPorosity_dVolStrain, - dPorosity_dPressure, - dPorosity_dTemperature, stack ); } @@ -189,45 +183,23 @@ GEOS_FORCE_INLINE void ThermalSinglePhasePoromechanics< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >:: computeFluidIncrement( localIndex const k, localIndex const q, - real64 const & porosity, - real64 const & porosity_n, real64 const & dPorosity_dVolStrain, - real64 const & dPorosity_dPressure, - real64 const & dPorosity_dTemperature, StackVariables & stack ) const { // Step 1: compute fluid mass increment and its derivatives wrt vol strain and pressure Base::computeFluidIncrement( k, q, - porosity, - porosity_n, dPorosity_dVolStrain, - dPorosity_dPressure, - dPorosity_dTemperature, stack ); // Step 2: compute derivative of fluid mass increment wrt temperature - stack.dFluidMassIncrement_dTemperature = dPorosity_dTemperature * m_fluidDensity( k, q ) + porosity * m_dFluidDensity_dTemperature( k, q ); - - // Step 3: compute fluid energy increment and its derivatives wrt vol strain, pressure, and temperature - real64 const fluidMass = porosity * m_fluidDensity( k, q ); - real64 const fluidEnergy = fluidMass * m_fluidInternalEnergy( k, q ); - real64 const fluidEnergy_n = porosity_n * m_fluidDensity_n( k, q ) * m_fluidInternalEnergy_n( k, q ); - stack.energyIncrement = fluidEnergy - fluidEnergy_n; - - stack.dEnergyIncrement_dVolStrainIncrement = stack.dFluidMassIncrement_dVolStrainIncrement * m_fluidInternalEnergy( k, q ); - stack.dEnergyIncrement_dPressure = stack.dFluidMassIncrement_dPressure * m_fluidInternalEnergy( k, q ) - + fluidMass * m_dFluidInternalEnergy_dPressure( k, q ); - stack.dEnergyIncrement_dTemperature = stack.dFluidMassIncrement_dTemperature * m_fluidInternalEnergy( k, q ) - + fluidMass * m_dFluidInternalEnergy_dTemperature( k, q ); - - - // Step 4: assemble the solid part of the accumulation term - real64 const oneMinusPoro = 1 - porosity; - - stack.energyIncrement += oneMinusPoro * m_rockInternalEnergy( k, 0 ) - ( 1 - porosity_n ) * m_rockInternalEnergy_n( k, 0 ); - stack.dEnergyIncrement_dVolStrainIncrement += -dPorosity_dVolStrain * m_rockInternalEnergy( k, 0 ); - stack.dEnergyIncrement_dPressure += -dPorosity_dPressure * m_rockInternalEnergy( k, 0 ); - stack.dEnergyIncrement_dTemperature += -dPorosity_dTemperature * m_rockInternalEnergy( k, 0 ) + oneMinusPoro * m_dRockInternalEnergy_dTemperature( k, 0 ); + stack.dFluidMassIncrement_dTemperature = m_dFluidMass_dTemperature[k]; + + // Step 3: compute energy increment and its derivatives wrt pressure, and temperature, vol strain + stack.energyIncrement = m_energy[k] - m_energy_n[k]; + stack.dEnergyIncrement_dPressure = m_dEnergy_dPressure[k]; + stack.dEnergyIncrement_dTemperature = m_dEnergy_dTemperature[k]; + stack.dEnergyIncrement_dVolStrainIncrement = stack.dFluidMassIncrement_dVolStrainIncrement * m_fluidInternalEnergy( k, q ) - + dPorosity_dVolStrain * m_rockInternalEnergy( k, 0 ); } template< typename SUBREGION_TYPE,