remove tjb refs and add DerivOffset::dxxx index tags

src/coreComponents/constitutive/fluid/singlefluid/ProppantSlurryFluid.hpp

@@ -166,11 +166,11 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate
              m_dFluidVisc_dCompConc[k][q],
              isProppantBoundary,
              m_density[k][q],
-             m_dDensity[k][q][DerivOffset::dP],  // tjb add deriv:dp
+             m_dDensity[k][q][DerivOffset::dP],
              m_dDensity_dProppantConc[k][q],
              m_dDensity_dCompConc[k][q],
              m_viscosity[k][q],
-             m_dViscosity[k][q][DerivOffset::dP],// tjb add deriv:dp
+             m_dViscosity[k][q][DerivOffset::dP],
              m_dViscosity_dProppantConc[k][q],
              m_dViscosity_dCompConc[k][q] );
   }
@@ -262,11 +262,11 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate
                 arraySlice1d< real64 const > const & dFluidViscosity_dComponentConcentration,
                 integer const & isProppantBoundary,
                 real64 & density,
-                real64 & dDensity_dp,  // tjb
+                real64 & dDensity_dp,
                 real64 & dDensity_dProppantConcentration,
                 arraySlice1d< real64 > const & dDensity_dComponentConcentration,
                 real64 & viscosity,
-                real64 & dViscosity_dp,  // tjb
+                real64 & dViscosity_dp,
                 real64 & dViscosity_dProppantConcentration,
                 arraySlice1d< real64 > const & dViscosity_dComponentConcentration ) const;
 
@@ -479,11 +479,11 @@ ProppantSlurryFluidUpdate::
            arraySlice1d< real64 const > const & GEOS_UNUSED_PARAM( dFluidViscosity_dComponentConcentration ),
            integer const & isProppantBoundary,
            real64 & density,
-           real64 & dDensity_dp, //tjb
+           real64 & dDensity_dp,
            real64 & dDensity_dProppantConcentration,
            arraySlice1d< real64 > const & dDensity_dComponentConcentration,
            real64 & viscosity,
-           real64 & dViscosity_dp, //tjb
+           real64 & dViscosity_dp,
            real64 & dViscosity_dProppantConcentration,
            arraySlice1d< real64 > const & dViscosity_dComponentConcentration ) const
 {

src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.cpp

@@ -117,7 +117,8 @@ void SlurryFluidBase::allocateConstitutiveData( Group & parent,
   this->resize( parent.size() );
 
 
-  // tjb these are also sized in m_dDenisty in base class
+  // These are also sized in m_dDenisty in base class , only dP and dT are populated
+  // Future dev should incorporate concentration derivatives in dDensity
   m_dDensity_dProppantConc.resize( parent.size(), numConstitutivePointsPerParentIndex );
   m_dDensity_dCompConc.resize( parent.size(), numConstitutivePointsPerParentIndex, NC );
 

src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AquiferBCKernel.hpp

@@ -89,6 +89,7 @@ struct AquiferBCKernel
           arrayView1d< real64 > const & localRhs )
   {
     using Order = BoundaryStencil::Order;
+    using Deriv = constitutive::singlefluid::DerivativeOffset;
 
     BoundaryStencil::IndexContainerViewConstType const & seri = stencil.getElementRegionIndices();
     BoundaryStencil::IndexContainerViewConstType const & sesri = stencil.getElementSubRegionIndices();
@@ -121,7 +122,7 @@ struct AquiferBCKernel
                                 dAquiferVolFlux_dPres,
                                 aquiferDens,
                                 dens[er][esr][ei][0],
-                                dDens[er][esr][ei][0][0], // tjb
+                                dDens[er][esr][ei][0][Deriv::dP],
                                 dt,
                                 localFlux,
                                 localFluxJacobian );

src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/DirichletFluxComputeKernel.hpp

@@ -48,7 +48,7 @@ class DirichletFluxComputeKernel : public FluxComputeKernel< NUM_EQN, NUM_DOF,
                                                              BoundaryStencilWrapper >
 {
 public:
-
+  using Deriv = constitutive::singlefluid::DerivativeOffset;
   using AbstractBase = singlePhaseFVMKernels::FluxComputeKernelBase;
   using DofNumberAccessor = AbstractBase::DofNumberAccessor;
   using PermeabilityAccessors = AbstractBase::PermeabilityAccessors;
@@ -204,8 +204,7 @@ class DirichletFluxComputeKernel : public FluxComputeKernel< NUM_EQN, NUM_DOF,
 
     // Compute average density
     real64 const densMean = 0.5 * ( m_dens[er][esr][ei][0] + faceDens );
-    //real64 const dDens_dP = 0.5 * m_dDens_dPres[er][esr][ei][0];
-    real64 const dDens_dP = 0.5 * m_dDens[er][esr][ei][0][0]; // tjb
+    real64 const dDens_dP = 0.5 * m_dDens[er][esr][ei][0][Deriv::dP];
 
     // Evaluate potential difference
     real64 const potDif = (m_pres[er][esr][ei] - m_facePres[kf])

src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernelBase.hpp

@@ -63,7 +63,7 @@ class FluxComputeKernelBase
   using DofNumberAccessor = ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > >;
 
   using SingleFluidProp = constitutive::SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >;
-  using DerivOffset = constitutive::singlefluid::DerivativeOffset; // tjb why no tp
+  using DerivOffset = constitutive::singlefluid::DerivativeOffset;
   using SinglePhaseFlowAccessors =
     StencilAccessors< fields::ghostRank,
                       fields::flow::pressure,

src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp

@@ -39,7 +39,6 @@ namespace singlePhaseFluxKernelsHelper
 template< typename VIEWTYPE >
 using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
-
 GEOS_HOST_DEVICE
 inline
 void computeSinglePhaseFlux( localIndex const ( &seri )[2],
@@ -69,10 +68,7 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2],
   for( localIndex ke = 0; ke < 2; ++ke )
   {
     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][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,9 +108,6 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2],
     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]][DerivOffset::dP];
-    //tjb remove
-    //dMobility_dP[ke] = dMob_dPres[seri[ke]][sesri[ke]][sei[ke]];
-    //assert( fabs( dMob[seri[ke]][sesri[ke]][sei[ke]][DerivOffset::dP]-dMob_dPres[seri[ke]][sesri[ke]][sei[ke]] )<FLT_EPSILON );
   }
   else
   {
@@ -124,9 +117,6 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2],
     {
       mobility += mobWeights[ke] * mob[seri[ke]][sesri[ke]][sei[ke]];
       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]][DerivOffset::dP]-dMob_dPres[seri[ke]][sesri[ke]][sei[ke]] )<FLT_EPSILON );
     }
   }
 
@@ -168,12 +158,12 @@ void computeEnthalpyFlux( localIndex const ( &seri )[2],
                           ENERGYFLUX_DERIVATIVE_TYPE & dEnergyFlux_dT )
 {
   // Step 1: compute the derivatives of the mean density at the interface wrt temperature
-
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
   real64 dDensMean_dT[2]{0.0, 0.0};
 
   for( integer ke = 0; ke < 2; ++ke )
   {
-    real64 const dDens_dT = dDens[seri[ke]][sesri[ke]][sei[ke]][0][0]; // tjb tag
+    real64 const dDens_dT = dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT];
     dDensMean_dT[ke] = 0.5 * dDens_dT;
   }
 
@@ -245,17 +235,17 @@ void computeEnthalpyFlux( localIndex const ( &seri )[2],
     localIndex const k_up = 1 - localIndex( fmax( fmin( alpha, 1.0 ), 0.0 ) );
 
     enthalpyTimesMobWeight = enthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0];
-    dEnthalpy_dP[k_up] = dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][0]; // tjb tag
-    dEnthalpy_dT[k_up] = dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][0]; // tjb tag
+    dEnthalpy_dP[k_up] = dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][DerivOffset::dP];
+    dEnthalpy_dT[k_up] = dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][DerivOffset::dT];
   }
   else
   {
     real64 const mobWeights[2] = { alpha, 1.0 - alpha };
     for( integer ke = 0; ke < 2; ++ke )
     {
       enthalpyTimesMobWeight += mobWeights[ke] * enthalpy[seri[ke]][sesri[ke]][sei[ke]][0];
-      dEnthalpy_dP[ke] = mobWeights[ke] * dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][0]; // tjb
-      dEnthalpy_dT[ke] = mobWeights[ke] * dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][0]; // tjb
+      dEnthalpy_dP[ke] = mobWeights[ke] * dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dP];
+      dEnthalpy_dT[ke] = mobWeights[ke] * dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT];
     }
   }
 

src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernels.hpp

@@ -232,6 +232,8 @@ class ElementBasedAssemblyKernel
    * Can be converted from ElementRegionManager::ElementViewConstAccessor
    * by calling .toView() or .toViewConst() on an accessor instance
    */
+
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 0 >;
   template< typename VIEWTYPE >
   using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
@@ -379,7 +381,7 @@ class ElementBasedAssemblyKernel
         real64 const fGravCoef = m_faceGravCoef[m_elemToFaces[ei][jFaceLoc]];
 
         real64 const ccDens = m_elemDens[ei][0];
-        real64 const dCcDens_dPres = m_dElemDens[ei][0][0]; //tjb
+        real64 const dCcDens_dPres = m_dElemDens[ei][0][DerivOffset::dP];
         // no density evaluated at the face center
 
         // pressure difference

src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp

@@ -49,6 +49,8 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo
  * Can be converted from ElementRegionManager::ElementViewConstAccessor
  * by calling .toView() or .toViewConst() on an accessor instance
  */
+
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
   template< typename VIEWTYPE >
   using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
@@ -203,7 +205,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo
 
       // Compute the derivatives of the density wrt temperature
 
-      real64 const dDens_dT = 0.5 * m_dDens[er][esr][ei][0][1]; // tjb tag
+      real64 const dDens_dT = 0.5 * m_dDens[er][esr][ei][0][DerivOffset::dT];
       // Compute the derivatives of the phase potential difference wrt temperature
 
       real64 const dF_dT = -stack.transmissibility * dDens_dT * ( m_gravCoef[er][esr][ei] - m_faceGravCoef[kf] );
@@ -221,8 +223,8 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo
         real64 const dFlux_dP = mobility_up * dF_dP + dMobility_dP_up * f;
         real64 const dFlux_dT = mobility_up * dF_dT + m_dMob_dTemp[er][esr][ei] * f;
 
-        stack.dEnergyFlux_dP += dFlux_dP * enthalpy + flux * m_dEnthalpy[er][esr][ei][0][0];  // tjb add tags
-        stack.dEnergyFlux_dT += dFlux_dT * enthalpy + flux * m_dEnthalpy[er][esr][ei][0][1];  // tjb add tags
+        stack.dEnergyFlux_dP += dFlux_dP * enthalpy + flux * m_dEnthalpy[er][esr][ei][0][DerivOffset::dP];
+        stack.dEnergyFlux_dT += dFlux_dT * enthalpy + flux * m_dEnthalpy[er][esr][ei][0][DerivOffset::dT];
       }
       else
       {

src/coreComponents/physicsSolvers/fluidFlow/proppantTransport/ProppantTransportKernels.cpp

@@ -247,11 +247,11 @@ FluxKernel::
     localIndex const ei  = stencilElementIndices[i];
 
     edgeDensity += geometricWeight[i] * dens[ei][0];
-    dEdgeDens_dP[i] = geometricWeight[i] * dDens[ei][0][0];  // tjb
+    dEdgeDens_dP[i] = geometricWeight[i] * dDens[ei][0][DerivOffset::dP];
     dEdgeDens_dProppantC[i] = geometricWeight[i] * dDens_dProppantConc[ei][0];
 
     edgeViscosity += geometricWeight[i] * visc[ei][0];
-    dEdgeVisc_dP[i] = geometricWeight[i] * dVisc[ei][0][0];
+    dEdgeVisc_dP[i] = geometricWeight[i] * dVisc[ei][0][DerivOffset::dP];
     dEdgeVisc_dProppantC[i] = geometricWeight[i] * dVisc_dProppantConc[ei][0];
 
     proppantC[i] = proppantConc[ei];
@@ -822,7 +822,7 @@ void FluxKernel::
                        dComponentDens_dComponentConc[er][esr],
                        gravDepth[er][esr],
                        dens[er][esr],
-                       dDens[er][esr],  // tjb tag
+                       dDens[er][esr],
                        dDens_dProppantConc[er][esr],
                        dDens_dComponentConc[er][esr],
                        visc[er][esr],

src/coreComponents/physicsSolvers/fluidFlow/proppantTransport/ProppantTransportKernels.hpp

@@ -224,6 +224,7 @@ struct FluxKernel
                               fields::permeability::permeability,
                               fields::permeability::permeabilityMultiplier >;
 
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 0 >;
   /**
    * @brief The type for element-based non-constitutive data parameters.
    * Consists entirely of ArrayView's.

src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp

@@ -193,7 +193,7 @@ void SinglePhaseWell::updateBHPForConstraint( WellElementSubRegion & subRegion )
                               &refGravCoef] ( localIndex const )
   {
     currentBHP = pres[iwelemRef] + dens[iwelemRef][0] * ( refGravCoef - wellElemGravCoef[iwelemRef] );
-    dCurrentBHP_dPres = 1.0 + dDens[iwelemRef][0][0] * ( refGravCoef - wellElemGravCoef[iwelemRef] );  // tjb add tag
+    dCurrentBHP_dPres = 1.0 + dDens[iwelemRef][0][DerivOffset::dP] * ( refGravCoef - wellElemGravCoef[iwelemRef] );
   } );
 
   if( logLevel >= 2 )
@@ -289,7 +289,7 @@ void SinglePhaseWell::updateVolRateForConstraint( WellElementSubRegion & subRegi
 
       real64 const densInv = 1.0 / dens[iwelemRef][0];
       currentVolRate = connRate[iwelemRef] * densInv;
-      dCurrentVolRate_dPres = -( useSurfaceConditions ==  0 ) * dDens[iwelemRef][0][0] * currentVolRate * densInv;  // tjb add tag
+      dCurrentVolRate_dPres = -( useSurfaceConditions ==  0 ) * dDens[iwelemRef][0][DerivOffset::dP] * currentVolRate * densInv;
       dCurrentVolRate_dRate = densInv;
 
       if( logLevel >= 2 && useSurfaceConditions )

src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp

@@ -22,6 +22,8 @@
 
 #include "WellSolverBase.hpp"
 
+#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp"
+
 namespace geos
 {
 
@@ -45,7 +47,7 @@ class SinglePhaseWell : public WellSolverBase
 {
 public:
 
-
+  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
   /**
    * @brief main constructor for Group Objects
    * @param name the name of this instantiation of Group in the repository

src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp

@@ -1036,9 +1036,7 @@ class ElementBasedAssemblyKernel
     {
       real64 const phaseAmount = stack.volume * phaseVolFrac[ip] * phaseDens[ip];
       real64 const phaseAmount_n = stack.volume * phaseVolFrac_n[ip] * phaseDens_n[ip];
-      //remove tjb
-      real64 const dPhaseAmount_dP = stack.volume * ( dPhaseVolFrac[ip][Deriv::dP] * phaseDens[ip]
-                                                      + phaseVolFrac[ip] * dPhaseDens[ip][Deriv::dP] );
+
       dPhaseAmount[FLUID_PROP_COFFSET::dP]=stack.volume * ( dPhaseVolFrac[ip][Deriv::dP] * phaseDens[ip]
                                                             + phaseVolFrac[ip] * dPhaseDens[ip][Deriv::dP] );
 
@@ -1054,20 +1052,14 @@ class ElementBasedAssemblyKernel
                                                   + phaseDens[ip] * dPhaseVolFrac[ip][Deriv::dC+jc];
         dPhaseAmount[FLUID_PROP_COFFSET::dC+jc] *= stack.volume;
       }
-// tjb- remove when safe
-      for( integer ic = 0; ic < numComp; ic++ )
-      {
-        assert( fabs( dPhaseAmount[FLUID_PROP_COFFSET::dC+ic] -dPhaseAmount_dC[ic] ) < FLT_EPSILON );
-
-      }
       // ic - index of component whose conservation equation is assembled
       // (i.e. row number in local matrix)
       for( integer ic = 0; ic < numComp; ++ic )
       {
         real64 const phaseCompAmount = phaseAmount * phaseCompFrac[ip][ic];
         real64 const phaseCompAmount_n = phaseAmount_n * phaseCompFrac_n[ip][ic];
 
-        real64 const dPhaseCompAmount_dP = dPhaseAmount_dP * phaseCompFrac[ip][ic]
+        real64 const dPhaseCompAmount_dP = dPhaseAmount[FLUID_PROP_COFFSET::dP] * phaseCompFrac[ip][ic]
                                            + phaseAmount * dPhaseCompFrac[ip][ic][Deriv::dP];
 
         stack.localResidual[ic] += phaseCompAmount - phaseCompAmount_n;

src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp

@@ -253,6 +253,7 @@ PressureRelationKernel::
           CRSMatrixView< real64, globalIndex const > const & localMatrix,
           arrayView1d< real64 > const & localRhs )
 {
+  using Deriv = constitutive::singlefluid::DerivativeOffset;
   // static well control data
   bool const isProducer = wellControls.isProducer();
   WellControls::Control const currentControl = wellControls.getControl();
@@ -317,8 +318,8 @@ PressureRelationKernel::
 
       // compute avg density
       real64 const avgDensity = 0.5 * ( wellElemDensity[iwelem][0] + wellElemDensity[iwelemNext][0] );
-      real64 const dAvgDensity_dPresNext    = 0.5 * dWellElemDensity[iwelemNext][0][0]; // tjb add tag
-      real64 const dAvgDensity_dPresCurrent = 0.5 * dWellElemDensity[iwelem][0][0];// tjb add tag
+      real64 const dAvgDensity_dPresNext    = 0.5 * dWellElemDensity[iwelemNext][0][Deriv::dP];
+      real64 const dAvgDensity_dPresCurrent = 0.5 * dWellElemDensity[iwelem][0][Deriv::dP];
 
       // compute depth diff times acceleration
       real64 const gravD = wellElemGravCoef[iwelemNext] - wellElemGravCoef[iwelem];
@@ -481,7 +482,7 @@ PerforationKernel::
           arrayView1d< real64 > const & perfRate,
           arrayView2d< real64 > const & dPerfRate_dPres )
 {
-
+  using Deriv = constitutive::singlefluid::DerivativeOffset;
   forAll< parallelDevicePolicy<> >( size, [=] GEOS_HOST_DEVICE ( localIndex const iperf )
   {
 
@@ -495,15 +496,15 @@ PerforationKernel::
 
     compute( resPressure[er][esr][ei],
              resDensity[er][esr][ei][0],
-             dResDensity[er][esr][ei][0][0],  // tjb add tag
+             dResDensity[er][esr][ei][0][Deriv::dP],
              resViscosity[er][esr][ei][0],
-             dResViscosity[er][esr][ei][0][0], // tjb add tag
+             dResViscosity[er][esr][ei][0][Deriv::dP],
              wellElemGravCoef[iwelem],
              wellElemPressure[iwelem],
              wellElemDensity[iwelem][0],
-             dWellElemDensity[iwelem][0][0],  // tjb add tag
+             dWellElemDensity[iwelem][0][Deriv::dP],
              wellElemViscosity[iwelem][0],
-             dWellElemViscosity[iwelem][0][0], // tjb add tag
+             dWellElemViscosity[iwelem][0][Deriv::dP],
              perfGravCoef[iperf],
              perfTransmissibility[iperf],
              perfRate[iperf],
@@ -529,6 +530,7 @@ AccumulationKernel::
           CRSMatrixView< real64, globalIndex const > const & localMatrix,
           arrayView1d< real64 > const & localRhs )
 {
+  using Deriv = constitutive::singlefluid::DerivativeOffset;
   forAll< parallelDevicePolicy<> >( size, [=] GEOS_HOST_DEVICE ( localIndex const iwelem )
   {
 
@@ -541,7 +543,7 @@ AccumulationKernel::
     globalIndex const presDofColIndex = wellElemDofNumber[iwelem] + COFFSET::DPRES;
 
     real64 const localAccum = wellElemVolume[iwelem] * ( wellElemDensity[iwelem][0] - wellElemDensity_n[iwelem][0] );
-    real64 const localAccumJacobian = wellElemVolume[iwelem] * dWellElemDensity[iwelem][0][0]; // tjb add tag
+    real64 const localAccumJacobian = wellElemVolume[iwelem] * dWellElemDensity[iwelem][0][Deriv::dP];
 
     // add contribution to global residual and jacobian (no need for atomics here)
     localMatrix.addToRow< serialAtomic >( eqnRowIndex, &presDofColIndex, &localAccumJacobian, 1 );

src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp

@@ -22,6 +22,7 @@
 
 #include "constitutive/fluid/singlefluid/SingleFluidFields.hpp"
 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
+#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp"
 #include "common/DataTypes.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
 #include "mesh/ElementRegionManager.hpp"