use DerivOffset:dP dT

GEOS-DEV · Dec 20, 2024 · d36cae9 · d36cae9
1 parent 77979c0
commit d36cae9
Show file tree

Hide file tree

Showing 6 changed files with 51 additions and 37 deletions.
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ProppantSlurryFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ProppantSlurryFluid.hpp
@@ -35,6 +35,8 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate
 {
 public:
   using SingleFluidProp = SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >;
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC<0>;
+
   /**
    * @brief
    * @param compressibility
@@ -168,12 +170,12 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate
              m_dFluidVisc_dCompConc[k][q],
              isProppantBoundary,
              m_density[k][q],
-             m_dDensity[k][q][0],  // tjb add deriv:dp
+             m_dDensity[k][q][DerivOffset::dP],  // tjb add deriv:dp
              m_dDensity_dPressure[k][q],
              m_dDensity_dProppantConc[k][q],
              m_dDensity_dCompConc[k][q],
              m_viscosity[k][q],
-             m_dViscosity[k][q][0],// tjb add deriv:dp
+             m_dViscosity[k][q][DerivOffset::dP],// tjb add deriv:dp
              m_dViscosity_dPressure[k][q],
              m_dViscosity_dProppantConc[k][q],
              m_dViscosity_dCompConc[k][q] );

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp
@@ -47,13 +47,18 @@ class AccumulationKernel
 
 public:
   using SingleFluidProp = constitutive::SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >;
-  using DerivOffset = constitutive::singlefluid::DerivativeOffset;
+
   /// Compute time value for the number of degrees of freedom
   static constexpr integer numDof = NUM_DOF;
 
   /// Compute time value for the number of equations
   static constexpr integer numEqn = NUM_DOF;
 
+
+  /// Note: Derivative lineup only supports dP & dT, not component terms 
+  static constexpr integer isThermal = NUM_DOF-1;
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC<isThermal>;
+
   /**
    * @brief Constructor
    * @param[in] rankOffset the offset of my MPI rank
@@ -169,11 +174,11 @@ class AccumulationKernel
     //std::cout << m_dDensity_dPres[ei][0]<< " " <<  m_dDensity[ei][0][DerivOffset::dP] << std::endl;
     //tjb  use DerivOffset::dP
     //    std::cout.flush();
-    assert( fabs( m_dDensity_dPres[ei][0]-m_dDensity[ei][0][0] )<FLT_EPSILON );
+    assert( fabs( m_dDensity_dPres[ei][0]-m_dDensity[ei][0][DerivOffset::dP] )<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][0] * stack.poreVolume;
+    stack.localJacobian[0][0] = stack.dPoreVolume_dPres * m_density[ei][0] + m_dDensity[ei][0][DerivOffset::dP] * 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
@@ -62,7 +62,7 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2],
                              real64 ( & dFlux_dP )[2],
                              real64 & dFlux_dTrans )
 {
-  //using DerivOffset = constitutive::singlefluid::DerivativeOffset;
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC<0>;
   // average density
   real64 densMean = 0.0;
   real64 dDensMean_dP[2];
@@ -72,8 +72,8 @@ 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];
     // tjb - add derivid
-    dDensMean_dP[ke] = 0.5 * dDens[seri[ke]][sesri[ke]][sei[ke]][0][0];
-    assert( fabs( dDens_dPres[seri[ke]][sesri[ke]][sei[ke]][0]-dDens[seri[ke]][sesri[ke]][sei[ke]][0][0] )<FLT_EPSILON );
+    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][DerivOffset::dP] )<FLT_EPSILON );
   }
 
   // compute potential difference
@@ -112,10 +112,10 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2],
     // happy path: single upwind direction
     localIndex const ke = 1 - localIndex( fmax( fmin( alpha, 1.0 ), 0.0 ) );
     mobility = mob[seri[ke]][sesri[ke]][sei[ke]];
-    dMobility_dP[ke] = dMob[seri[ke]][sesri[ke]][sei[ke]][0];
+    dMobility_dP[ke] = dMob[seri[ke]][sesri[ke]][sei[ke]][DerivOffset::dP];
     //tjb remove
     dMobility_dP[ke] = dMob_dPres[seri[ke]][sesri[ke]][sei[ke]];
-    assert( fabs( dMob[seri[ke]][sesri[ke]][sei[ke]][0]-dMob_dPres[seri[ke]][sesri[ke]][sei[ke]] )<FLT_EPSILON );
+    assert( fabs( dMob[seri[ke]][sesri[ke]][sei[ke]][DerivOffset::dP]-dMob_dPres[seri[ke]][sesri[ke]][sei[ke]] )<FLT_EPSILON );
   }
   else
   {
@@ -124,10 +124,10 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2],
     for( localIndex ke = 0; ke < 2; ++ke )
     {
       mobility += mobWeights[ke] * mob[seri[ke]][sesri[ke]][sei[ke]];
-      dMobility_dP[ke] = mobWeights[ke] * dMob[seri[ke]][sesri[ke]][sei[ke]][0];
+      dMobility_dP[ke] = mobWeights[ke] * dMob[seri[ke]][sesri[ke]][sei[ke]][DerivOffset::dP];
       //tjb remove
       dMobility_dP[ke] = mobWeights[ke] * dMob_dPres[seri[ke]][sesri[ke]][sei[ke]];
-      assert( fabs( dMob[seri[ke]][sesri[ke]][sei[ke]][0]-dMob_dPres[seri[ke]][sesri[ke]][sei[ke]] )<FLT_EPSILON );
+      assert( fabs( dMob[seri[ke]][sesri[ke]][sei[ke]][DerivOffset::dP]-dMob_dPres[seri[ke]][sesri[ke]][sei[ke]] )<FLT_EPSILON );
     }
   }
 

diff --git a/...oreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp b/...oreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp
@@ -56,6 +56,9 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
   using Base::m_localMatrix;
   using Base::m_localRhs;
 
+  /// Note: Derivative lineup only supports dP & dT, not component terms 
+  static constexpr integer isThermal = NUM_DOF-1;
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC<isThermal>;
   /**
    * @brief Constructor
    * @param[in] rankOffset the offset of my MPI rank
@@ -166,28 +169,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 );
+      assert( fabs( m_dDensity_dPres[ei][0]-m_dDensity[ei][0][DerivOffset::dP] )<FLT_EPSILON );
+      assert( fabs( m_dInternalEnergy_dPres[ei][0]-m_dInternalEnergy[ei][0][DerivOffset::dP] )<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];
+                                     + stack.poreVolume * m_dDensity[ei][0][DerivOffset::dP] * m_internalEnergy[ei][0]
+                                     + stack.poreVolume * m_density[ei][0] * m_dInternalEnergy[ei][0][DerivOffset::dP];
       // 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];
-      GEOS_UNUSED_VAR(xx);
+      GEOS_UNUSED_VAR( xx );
       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]
+      assert( fabs( m_dDensity_dTemp[ei][0]-m_dDensity[ei][0][DerivOffset::dT] )<FLT_EPSILON );
+      assert( fabs( m_dInternalEnergy_dTemp[ei][0]-m_dInternalEnergy[ei][0][DerivOffset::dT] )<FLT_EPSILON );
+      real64 const dFluidEnergy_dT = stack.poreVolume * m_dDensity[ei][0][DerivOffset::dT] * m_internalEnergy[ei][0]
+                                     + stack.poreVolume * m_density[ei][0] * m_dInternalEnergy[ei][0][DerivOffset::dT]
                                      + 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];
-      GEOS_UNUSED_VAR(yy);
+      GEOS_UNUSED_VAR( yy );
       assert( fabs( dFluidEnergy_dT-yy )<FLT_EPSILON );
       // local accumulation
       stack.localResidual[numEqn-1] += fluidEnergy;

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp
@@ -189,6 +189,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E
   void computeFlux( localIndex const iconn,
                     StackVariables & stack ) const
   {
+    using DerivOffset = constitutive::singlefluid::DerivativeOffsetC<1>;
     // ***********************************************
     // First, we call the base computeFlux to compute:
     //  1) compFlux and its derivatives (including derivatives wrt temperature),
@@ -218,8 +219,8 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E
       {
         //real64 const dDens_dT = m_dDens_dTemp[seri[ke]][sesri[ke]][sei[ke]][0];
         //tjb derivid
-        real64 const dDens_dT = m_dDens[seri[ke]][sesri[ke]][sei[ke]][0][1];
-        assert( fabs( m_dDens_dTemp[seri[ke]][sesri[ke]][sei[ke]][0]-m_dDens[seri[ke]][sesri[ke]][sei[ke]][0][1] )<FLT_EPSILON );
+        real64 const dDens_dT = m_dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT];
+        assert( fabs( m_dDens_dTemp[seri[ke]][sesri[ke]][sei[ke]][0]-m_dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT] )<FLT_EPSILON );
         dDensMean_dT[ke] = 0.5 * dDens_dT;
       }
 
@@ -265,18 +266,19 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E
       if( alpha <= 0.0 || alpha >= 1.0 )
       {
         localIndex const k_up = 1 - localIndex( fmax( fmin( alpha, 1.0 ), 0.0 ) );
-        dMob_dT[k_up] = m_dMob[seri[k_up]][sesri[k_up]][sei[k_up]][1];
+        dMob_dT[k_up] = m_dMob[seri[k_up]][sesri[k_up]][sei[k_up]][DerivOffset::dT];
         // tjb remove
-         assert( fabs( m_dMob_dTemp[seri[k_up]][sesri[k_up]][sei[k_up]]-m_dMob[seri[k_up]][sesri[k_up]][sei[k_up]][1] )<FLT_EPSILON );
+        assert( fabs( m_dMob_dTemp[seri[k_up]][sesri[k_up]][sei[k_up]]-m_dMob[seri[k_up]][sesri[k_up]][sei[k_up]][DerivOffset::dT] )<FLT_EPSILON );
         dMob_dT[k_up] = m_dMob_dTemp[seri[k_up]][sesri[k_up]][sei[k_up]];
       }
       else
       {
         real64 const mobWeights[2] = { alpha, 1.0 - alpha };
         for( integer ke = 0; ke < 2; ++ke )
         {
-          dMob_dT[ke] = mobWeights[ke] * m_dMob[seri[ke]][sesri[ke]][sei[ke]][1];
-          assert( fabs( m_dMob_dTemp[seri[ke]][sesri[ke]][sei[ke]]-m_dMob[seri[ke]][sesri[ke]][sei[ke]][1] )<FLT_EPSILON );
+          dMob_dT[ke] = mobWeights[ke] * m_dMob[seri[ke]][sesri[ke]][sei[ke]][DerivOffset::dT];
+          // tjb - remove
+          assert( fabs( m_dMob_dTemp[seri[ke]][sesri[ke]][sei[ke]]-m_dMob[seri[ke]][sesri[ke]][sei[ke]][DerivOffset::dT] )<FLT_EPSILON );
           dMob_dT[ke] = mobWeights[ke] * m_dMob_dTemp[seri[ke]][sesri[ke]][sei[ke]];
         }
       }
@@ -306,10 +308,10 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E
 
         enthalpy = m_enthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0];
         //tjb
-        assert( fabs( m_dEnthalpy_dPres[seri[k_up]][sesri[k_up]][sei[k_up]][0]-m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][0] )<FLT_EPSILON );
-        assert( fabs( m_dEnthalpy_dTemp[seri[k_up]][sesri[k_up]][sei[k_up]][0]-m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][1] )<FLT_EPSILON );
-        dEnthalpy_dP[k_up] = m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][0];
-        dEnthalpy_dT[k_up] = m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][1];
+        assert( fabs( m_dEnthalpy_dPres[seri[k_up]][sesri[k_up]][sei[k_up]][0]-m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][DerivOffset::dP] )<FLT_EPSILON );
+        assert( fabs( m_dEnthalpy_dTemp[seri[k_up]][sesri[k_up]][sei[k_up]][0]-m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][DerivOffset::dT] )<FLT_EPSILON );
+        dEnthalpy_dP[k_up] = m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][DerivOffset::dP];
+        dEnthalpy_dT[k_up] = m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][DerivOffset::dT];
       }
       else
       {
@@ -318,10 +320,10 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E
         {
           enthalpy += mobWeights[ke] * m_enthalpy[seri[ke]][sesri[ke]][sei[ke]][0];
           //tjb
-          assert( fabs( m_dEnthalpy_dPres[seri[ke]][sesri[ke]][sei[ke]][0]-m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][0] )<FLT_EPSILON );
-          assert( fabs( m_dEnthalpy_dTemp[seri[ke]][sesri[ke]][sei[ke]][0]-m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][1] )<FLT_EPSILON );
-          dEnthalpy_dP[ke] = mobWeights[ke] * m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][0];
-          dEnthalpy_dT[ke] = mobWeights[ke] * m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][1];
+          assert( fabs( m_dEnthalpy_dPres[seri[ke]][sesri[ke]][sei[ke]][0]-m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dP] )<FLT_EPSILON );
+          assert( fabs( m_dEnthalpy_dTemp[seri[ke]][sesri[ke]][sei[ke]][0]-m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT] )<FLT_EPSILON );
+          dEnthalpy_dP[ke] = mobWeights[ke] * m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dP];
+          dEnthalpy_dT[ke] = mobWeights[ke] * m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT];
         }
       }
 

diff --git a/src/coreComponents/unitTests/fluidFlowTests/testSinglePhaseMobilityKernel.cpp b/src/coreComponents/unitTests/fluidFlowTests/testSinglePhaseMobilityKernel.cpp
@@ -30,8 +30,10 @@ TEST( SinglePhaseBaseKernels, mobility )
   int constexpr NTEST = 3;
 
   real64 const dens[NTEST]        = { 800.0, 1000.0, 1500.0 };
+
   real64 const dDens_dPres[NTEST] = { 1e-5, 1e-10, 0.0    };
   real64 const visc[NTEST]        = { 5.0, 2.0, 1.0    };
+
   real64 const dVisc_dPres[NTEST] = { 1e-7, 0.0, 0.0    };
 
   for( int i = 0; i < NTEST; ++i )
@@ -41,7 +43,7 @@ TEST( SinglePhaseBaseKernels, mobility )
     real64 mob;
     real64 dMob_dPres;
 
-    MobilityKernel::compute( dens[i], dDens_dPres[i], visc[i], dVisc_dPres[i], mob, dMob_dPres );
+    MobilityKernel::compute( dens[i], dDens_dPres[i], dDens_dPres[i], visc[i], dVisc_dPres[i], dVisc_dPres[i], mob, dMob_dPres );
 
     // compute etalon
     real64 const mob_et = dens[i] / visc[i];