use deriv containers in mobility / accum kernels

GEOS-DEV · Dec 12, 2024 · 40511fc · 40511fc
1 parent 3ebe93b
commit 40511fc
Show file tree

Hide file tree

Showing 10 changed files with 65 additions and 15 deletions.
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.cpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.cpp
@@ -54,10 +54,12 @@ SingleFluidBase::SingleFluidBase( string const & name, Group * const parent )
   registerField( fields::singlefluid::dInternalEnergy_dPressure{}, &m_dInternalEnergy_dPressure );
   registerField( fields::singlefluid::dInternalEnergy_dTemperature{}, &m_dInternalEnergy_dTemperature );
 
-  registerField( fields::singlefluid::enthalpy{}, &m_enthalpy.value );
+  registerField( fields::singlefluid::enthalpy{}, &m_enthalpy.value );  
   registerField( fields::singlefluid::dEnthalpy{}, &m_enthalpy.derivs );
-  //registerField( fields::singlefluid::dEnthalpy_dPressure{}, &m_dEnthalpy_dPressure );
-  //registerField( fields::singlefluid::dEnthalpy_dTemperature{}, &m_dEnthalpy_dTemperature );
+#if 1 
+  registerField( fields::singlefluid::dEnthalpy_dPressure{}, &m_dEnthalpy_dPressure );
+  registerField( fields::singlefluid::dEnthalpy_dTemperature{}, &m_dEnthalpy_dTemperature );
+#endif
 
 }
 

diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidFields.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidFields.hpp
@@ -163,7 +163,7 @@ DECLARE_FIELD( dEnthalpy,
                LEVEL_0,
                WRITE_AND_READ,
                "dEnthalpy" );
-#if 0
+#if 1
 DECLARE_FIELD( dEnthalpy_dPressure,
                "dEnthalpy_dPressure",
                array2d< real64 >,

diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
@@ -156,6 +156,9 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate
              m_dDens_dPres[k][q],
              m_viscosity[k][q],
              m_dVisc_dPres[k][q] );
+    // tjb
+    m_dDensity[k][q][0]   = m_dDens_dPres[k][q];
+    m_dViscosity[k][q][0] = m_dVisc_dPres[k][q];
   }
 
   GEOS_HOST_DEVICE
@@ -179,6 +182,16 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate
              m_enthalpy[k][q],
              m_dEnthalpy_dPres[k][q],
              m_dEnthalpy_dTemp[k][q] );
+    // tjb
+    m_dDensity[k][q][0] = m_dDens_dPres[k][q];
+    m_dDensity[k][q][1] = m_dDens_dTemp[k][q];
+    m_dViscosity[k][q][0] = m_dVisc_dPres[k][q];
+    m_dViscosity[k][q][1] = m_dVisc_dTemp[k][q];
+    m_dInternalEnergy[k][q][0] = m_dIntEnergy_dPres[k][q];
+    m_dInternalEnergy[k][q][1] = m_dIntEnergy_dTemp[k][q];
+    m_dEnthalpy[k][q][0] = m_dEnthalpy_dPres[k][q];
+    m_dEnthalpy[k][q][1] = m_dEnthalpy_dTemp[k][q];
+
   }
 
 private:

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp
@@ -167,13 +167,13 @@ class AccumulationKernel
     stack.localResidual[0] = stack.poreVolume * m_density[ei][0] - m_mass_n[ei];
     // Derivative of residual wrt to pressure in the cell
     //std::cout << m_dDensity_dPres[ei][0]<< " " <<  m_dDensity[ei][0][DerivOffset::dP] << std::endl;
-    // tjb
+    //tjb  use DerivOffset::dP
     //    std::cout.flush();
-    //assert(fabs(m_dDensity_dPres[ei][0]-m_dDensity[ei][0][DerivOffset::dP])<FLT_EPSILON);
+    assert(fabs(m_dDensity_dPres[ei][0]-m_dDensity[ei][0][0])<FLT_EPSILON);
 
 
-    stack.localJacobian[0][0] = stack.dPoreVolume_dPres * m_density[ei][0] + m_dDensity_dPres[ei][0] * stack.poreVolume;
-    //stack.localJacobian[0][0] = stack.dPoreVolume_dPres * m_density[ei][0] + m_dDensity[ei][0][DerivOffset::dP] * stack.poreVolume;
+    //stack.localJacobian[0][0] = stack.dPoreVolume_dPres * m_density[ei][0] + m_dDensity_dPres[ei][0] * stack.poreVolume;
+    stack.localJacobian[0][0] = stack.dPoreVolume_dPres * m_density[ei][0] + m_dDensity[ei][0][0] * stack.poreVolume;
     // Customize the kernel with this lambda
     kernelOp();
   }

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp
@@ -70,7 +70,7 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2],
   {
     densMean        += 0.5 * dens[seri[ke]][sesri[ke]][sei[ke]][0];
     dDensMean_dP[ke] = 0.5 * dDens_dPres[seri[ke]][sesri[ke]][sei[ke]][0];
-    //dDensMean_dP[ke] = 0.5 * dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dP];
+    assert(fabs(dDens_dPres[seri[ke]][sesri[ke]][sei[ke]][0]-dDens[seri[ke]][sesri[ke]][sei[ke]][0][0])<FLT_EPSILON);
   }
 
   // compute potential difference

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/MobilityKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/MobilityKernel.hpp
@@ -59,10 +59,12 @@ struct MobilityKernel
   static void
   compute( real64 const & dens,
            real64 const & dDens_dP,  // tjb
+           real64 const & dDens_dT,  // tjb
            real64 const & dDens_dPres,
            real64 const & dDens_dTemp,
            real64 const & visc,
            real64 const & dVisc_dP,  // tjb
+           real64 const & dVisc_dT,  // tjb
            real64 const & dVisc_dPres,
            real64 const & dVisc_dTemp,
            real64 & mob,
@@ -73,7 +75,11 @@ struct MobilityKernel
     dMob_dPres = dDens_dPres / visc - mob / visc * dVisc_dPres;
     dMob_dPres = dDens_dP / visc - mob / visc * dVisc_dP;  // tjb keep
     assert( fabs( dDens_dP -dDens_dPres ) < FLT_EPSILON );
+    assert( fabs( dVisc_dP -dVisc_dPres ) < FLT_EPSILON );
     dMob_dTemp = dDens_dTemp / visc - mob / visc * dVisc_dTemp;
+    dMob_dTemp = dDens_dT / visc - mob / visc * dVisc_dT;  // tjb keep
+    assert( fabs( dDens_dT - dDens_dTemp ) < FLT_EPSILON );
+    assert( fabs( dVisc_dT - dVisc_dTemp ) < FLT_EPSILON );
   }
 
 // Value-only (no derivatives) version
@@ -131,10 +137,12 @@ struct MobilityKernel
     {
       compute( dens[a][0],
                dDens[a][0][0], // tjb use deriv::dp
+               dDens[a][0][1], // tjb use deriv::dt
                dDens_dPres[a][0],
                dDens_dTemp[a][0],
                visc[a][0],
                dVisc[a][0][0], // tjb use deriv::dp
+               dVisc[a][0][1], // tjb use deriv::dt
                dVisc_dPres[a][0],
                dVisc_dTemp[a][0],
                mob[a],

diff --git a/...oreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp b/...oreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp
@@ -51,6 +51,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
   using Base::m_porosity;
   using Base::m_dPoro_dPres;
   using Base::m_density;
+  using Base::m_dDensity;
   using Base::m_dDensity_dPres;
   using Base::m_localMatrix;
   using Base::m_localRhs;
@@ -76,6 +77,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
     m_dDensity_dTemp( fluid.dDensity_dTemperature() ),
     m_dPoro_dTemp( solid.getDporosity_dTemperature() ),
     m_internalEnergy( fluid.internalEnergy() ),
+    m_dInternalEnergy( fluid.dInternalEnergy() ),
     m_dInternalEnergy_dPres( fluid.dInternalEnergy_dPressure() ),
     m_dInternalEnergy_dTemp( fluid.dInternalEnergy_dTemperature() ),
     m_rockInternalEnergy( solid.getInternalEnergy() ),
@@ -163,15 +165,28 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
 
       // 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];
-
+      // tjb
+       assert(fabs(m_dDensity_dPres[ei][0]-m_dDensity[ei][0][0])<FLT_EPSILON);
+      assert(fabs(m_dInternalEnergy_dPres[ei][0]-m_dInternalEnergy[ei][0][0])<FLT_EPSILON);
       real64 const dFluidEnergy_dP = stack.dPoreVolume_dPres * m_density[ei][0] * m_internalEnergy[ei][0]
+                                     + stack.poreVolume * m_dDensity[ei][0][0] * m_internalEnergy[ei][0]
+                                     + stack.poreVolume * m_density[ei][0] * m_dInternalEnergy[ei][0][0];
+      // tjb delete
+      real64 xx = 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.poreVolume * m_density[ei][0] * m_dInternalEnergy_dPres[ei][0];
+      assert(fabs(dFluidEnergy_dP-xx)<FLT_EPSILON);
+      // tjb
+       assert(fabs(m_dDensity_dTemp[ei][0]-m_dDensity[ei][0][1])<FLT_EPSILON);
+      assert(fabs(m_dInternalEnergy_dTemp[ei][0]-m_dInternalEnergy[ei][0][1])<FLT_EPSILON);
+      real64 const dFluidEnergy_dT = stack.poreVolume * m_dDensity[ei][0][1] * m_internalEnergy[ei][0]
+                                     + stack.poreVolume * m_density[ei][0] * m_dInternalEnergy[ei][0][1]
                                      + stack.dPoreVolume_dTemp * m_density[ei][0] * m_internalEnergy[ei][0];
-
+      // tjb - delete
+      real64 yy   = 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];
+      assert(fabs(dFluidEnergy_dT-yy)<FLT_EPSILON);
       // local accumulation
       stack.localResidual[numEqn-1] += fluidEnergy;
 
@@ -216,6 +231,9 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
 
   /// Views on fluid internal energy
   arrayView2d< real64 const > const m_internalEnergy;
+  arrayView3d< real64 const > const m_dInternalEnergy;
+
+  //tjb remove
   arrayView2d< real64 const > const m_dInternalEnergy_dPres;
   arrayView2d< real64 const > const m_dInternalEnergy_dTemp;
 

diff --git a/...onents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp b/...onents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp
@@ -89,6 +89,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo
     StencilMaterialAccessors< constitutive::SingleFluidBase,
                               fields::singlefluid::dDensity_dTemperature,
                               fields::singlefluid::enthalpy,
+                              fields::singlefluid::dEnthalpy,
                               fields::singlefluid::dEnthalpy_dPressure,
                               fields::singlefluid::dEnthalpy_dTemperature >;
 
@@ -145,6 +146,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo
     m_dMob_dTemp( thermalSinglePhaseFlowAccessors.get( fields::flow::dMobility_dTemperature {} ) ),
     m_dDens_dTemp( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dDensity_dTemperature {} ) ),
     m_enthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::enthalpy {} ) ),
+    m_dEnthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) ),
     m_dEnthalpy_dPres( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy_dPressure {} ) ),
     m_dEnthalpy_dTemp( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy_dTemperature {} ) ),
     m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) ),
@@ -298,6 +300,8 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo
 
   /// Views on enthalpies
   ElementViewConst< arrayView2d< real64 const > > const m_enthalpy;
+  ElementViewConst< arrayView3d< real64 const > > const m_dEnthalpy;
+  // tjb remove
   ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dPres;
   ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dTemp;
 

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp
@@ -84,6 +84,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E
     StencilMaterialAccessors< constitutive::SingleFluidBase,
                               fields::singlefluid::dDensity_dTemperature,
                               fields::singlefluid::enthalpy,
+                              fields::singlefluid::dEnthalpy,
                               fields::singlefluid::dEnthalpy_dPressure,
                               fields::singlefluid::dEnthalpy_dTemperature >;
 
@@ -133,6 +134,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E
     m_dMob_dTemp( thermalSinglePhaseFlowAccessors.get( fields::flow::dMobility_dTemperature {} ) ),
     m_dDens_dTemp( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dDensity_dTemperature {} ) ),
     m_enthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::enthalpy {} ) ),
+    m_dEnthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) ),
     m_dEnthalpy_dPres( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy_dPressure {} ) ),
     m_dEnthalpy_dTemp( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy_dTemperature {} ) ),
     m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) ),
@@ -420,6 +422,8 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E
 
   /// Views on enthalpies
   ElementViewConst< arrayView2d< real64 const > > const m_enthalpy;
+  ElementViewConst< arrayView3d< real64 const > > const m_dEnthalpy;
+  // tjb remove
   ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dPres;
   ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dTemp;
 

diff --git a/...rs/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp b/...rs/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp
@@ -332,6 +332,7 @@ class ConnectorBasedAssemblyKernel : public singlePhasePoromechanicsEmbeddedFrac
 
   /// Views on enthalpies
   ElementViewConst< arrayView2d< real64 const > > const m_enthalpy;
+  ElementViewConst< arrayView3d< real64 const > > const m_dEnthalpy;
   ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dPres;
   ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dTemp;