20 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_C1PPUPHASEFLUX_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_C1PPUPHASEFLUX_HPP
25 #include "constitutive/fluid/multifluid/Layouts.hpp"
26 #include "constitutive/capillaryPressure/layouts.hpp"
34 namespace isothermalCompositionalMultiphaseFVMKernelUtilities
38 static constexpr
real64 epsC1PPU = 5000;
40 template<
typename VIEWTYPE >
41 using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
43 using Deriv = constitutive::multifluid::DerivativeOffset;
74 template<
integer numComp,
integer numFluxSupportPo
ints >
80 integer const checkPhasePresenceInGravity,
81 localIndex const ( &seri )[numFluxSupportPoints],
82 localIndex const ( &sesri )[numFluxSupportPoints],
83 localIndex const ( &sei )[numFluxSupportPoints],
84 real64 const ( &trans )[2],
85 real64 const ( &dTrans_dPres )[2],
99 real64 ( & dPhaseFlux_dP )[numFluxSupportPoints],
100 real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] )
102 real64 dPotGrad_dTrans = 0.0;
103 real64 dPresGrad_dP[numFluxSupportPoints]{};
104 real64 dPresGrad_dC[numFluxSupportPoints][numComp]{};
105 real64 dGravHead_dP[numFluxSupportPoints]{};
106 real64 dGravHead_dC[numFluxSupportPoints][numComp]{};
107 PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, checkPhasePresenceInGravity,
108 seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef,
109 phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens,
110 phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac,
111 potGrad, dPotGrad_dTrans, dPresGrad_dP,
112 dPresGrad_dC, dGravHead_dP, dGravHead_dC );
116 for(
integer i = 0; i < numFluxSupportPoints; i++ )
122 real64 const gravD = trans[i] * gravCoef[er][esr][ei];
133 real64 Ttrans = fabs( trans[0] );
134 potGrad = potGrad / Ttrans;
136 real64 const mobility_i = phaseMob[seri[0]][sesri[0]][sei[0]][ip];
137 real64 const mobility_j = phaseMob[seri[1]][sesri[1]][sei[1]][ip];
141 if( fabs( gravHead ) <= 1e-20 )
143 real64 const tmpSqrt = sqrt( potGrad * potGrad + smoEps * smoEps );
144 real64 const smoMax = 0.5 * (-potGrad + tmpSqrt);
146 phaseFlux = Ttrans * ( potGrad * mobility_i - smoMax * (mobility_j - mobility_i) );
154 real64 const dMob_dP = dPhaseMob[seri[0]][sesri[0]][sei[0]][ip][Deriv::dP];
155 dPhaseFlux_dP[0] += Ttrans * potGrad * dMob_dP;
161 dPhaseMobSub = dPhaseMob[seri[0]][sesri[0]][sei[0]][ip];
162 for(
integer jc = 0; jc < numComp; ++jc )
164 dPhaseFlux_dC[0][jc] += Ttrans * potGrad * dPhaseMobSub[Deriv::dC + jc];
168 real64 const tmpInv = 1.0 / tmpSqrt;
169 real64 const dSmoMax_x = 0.5 * (1.0 - potGrad * tmpInv);
172 real64 const dMobDiff_sign[numFluxSupportPoints] = {-1.0, 1.0};
173 for(
integer ke = 0; ke < numFluxSupportPoints; ++ke )
177 real64 const dPotGrad_dP = dPresGrad_dP[ke] - dGravHead_dP[ke];
180 dPhaseFlux_dP[ke] += dPotGrad_dP * mobility_i;
183 real64 const dSmoMax_dP = -dPotGrad_dP * dSmoMax_x;
184 dPhaseFlux_dP[ke] += -dSmoMax_dP * (mobility_j - mobility_i);
186 real64 const dMob_dP = dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][ip][Deriv::dP];
187 dPhaseFlux_dP[ke] += -Ttrans * smoMax * dMobDiff_sign[ke] * dMob_dP;
192 dPhaseMobSub = dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][ip];
194 for(
integer jc = 0; jc < numComp; ++jc )
196 real64 const dPotGrad_dC = dPresGrad_dC[ke][jc] - dGravHead_dC[ke][jc];
199 dPhaseFlux_dC[ke][jc] += dPotGrad_dC * mobility_i;
202 real64 const dSmoMax_dC = -dPotGrad_dC * dSmoMax_x;
203 dPhaseFlux_dC[ke][jc] += -dSmoMax_dC * (mobility_j - mobility_i);
204 dPhaseFlux_dC[ke][jc] += -Ttrans * smoMax * dMobDiff_sign[ke] * dPhaseMobSub[Deriv::dC + jc];
208 potGrad = potGrad * Ttrans;
#define GEOS_HOST_DEVICE
Marks a host-device function.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
ArraySlice< T, 1, USD > arraySlice1d
Alias for 1D array slice.
std::int32_t integer
Signed integer type.
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
static GEOS_HOST_DEVICE void compute(integer const numPhase, integer const ip, integer const hasCapPressure, integer const checkPhasePresenceInGravity, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&trans)[2], real64 const (&dTrans_dPres)[2], ElementViewConst< arrayView1d< real64 const > > const &pres, ElementViewConst< arrayView1d< real64 const > > const &gravCoef, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseMob, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseVolFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &dCompFrac_dCompDens, ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &phaseMassDens, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &dPhaseMassDens, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &dPhaseCapPressure_dPhaseVolFrac, real64 &potGrad, real64(&phaseFlux), real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp])
Form the PhasePotentialUpwind from pressure gradient and gravitational head.