diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp index ba4268e5eaa..dd5d1e61b1b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp @@ -176,6 +176,81 @@ class FluxComputeKernelBase BitFlags< KernelFlags > const m_kernelFlags; }; +namespace helpers +{ + template< typename VIEWTYPE > + using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + +template< localIndex numComp, localIndex numFluxSupportPoints > +GEOS_HOST_DEVICE +static void calculateMeanDensity( localIndex const ip, + localIndex const (&seri)[numFluxSupportPoints], + localIndex const (&sesri)[numFluxSupportPoints], + localIndex const (&sei)[numFluxSupportPoints], + integer const checkPhasePresenceInGravity, + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, + 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, + real64 & densMean, real64 (& dDensMean_dPres)[numFluxSupportPoints], real64 (& dDensMean_dComp)[numFluxSupportPoints][numComp] ) +{ + using Deriv = constitutive::multifluid::DerivativeOffset; + + densMean = 0; + integer denom = 0; + real64 dDens_dC[numComp]{}; + for( localIndex i = 0; i < numFluxSupportPoints; ++i ) + { + localIndex const er = seri[i]; + localIndex const esr = sesri[i]; + localIndex const ei = sei[i]; + + bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); + if( checkPhasePresenceInGravity && !phaseExists ) + { + dDensMean_dPres[i] = 0.0; + for( localIndex jc = 0; jc < numComp; ++jc ) + { + dDensMean_dComp[i][jc] = 0.0; + } + continue; + } + + // density + real64 const density = phaseMassDens[er][esr][ei][0][ip]; + real64 const dDens_dPres = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP]; + + applyChainRule( numComp, + dCompFrac_dCompDens[er][esr][ei], + dPhaseMassDens[er][esr][ei][0][ip], + dDens_dC, + Deriv::dC ); + + // average density and derivatives + densMean += density; + dDensMean_dPres[i] = dDens_dPres; + for( localIndex jc = 0; jc < numComp; ++jc ) + { + dDensMean_dComp[i][jc] = dDens_dC[jc]; + } + denom++; + } + if( denom > 1 ) + { + densMean /= denom; + for( localIndex i = 0; i < numFluxSupportPoints; ++i ) + { + dDensMean_dPres[i] /= denom; + for( integer jc = 0; jc < numComp; ++jc ) + { + dDensMean_dComp[i][jc] /= denom; + } + } + } +} + +} + } // namespace isothermalCompositionalMultiphaseFVMKernels } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/HU2PhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/HU2PhaseFlux.hpp index 04cd797c29f..cdc43b58f86 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/HU2PhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/HU2PhaseFlux.hpp @@ -78,7 +78,7 @@ struct HU2PhaseFlux compute( integer const numPhase, integer const ip, integer const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -103,7 +103,7 @@ struct HU2PhaseFlux // viscous part computeViscousFlux< numComp, numFluxSupportPoints >( ip, numPhase, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, @@ -115,8 +115,10 @@ struct HU2PhaseFlux phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); // gravity part computeGravityFlux< numComp, numFluxSupportPoints >( ip, numPhase, + checkPhasePresenceInGravity, seri, sesri, sei, trans, dTrans_dPres, gravCoef, + phaseVolFrac, phaseMob, dPhaseMob, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, @@ -143,7 +145,7 @@ struct HU2PhaseFlux computeViscousFlux( integer const & ip, integer const & numPhase, integer const & hasCapPressure, - integer const & useNewGravity, + integer const & checkPhasePresenceInGravity, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], @@ -170,7 +172,7 @@ struct HU2PhaseFlux computeTotalFlux( numPhase, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, @@ -259,11 +261,13 @@ struct HU2PhaseFlux GEOS_HOST_DEVICE static void computeGravityFlux( integer const & ip, integer const & numPhase, + 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 & gravCoef, + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, ElementViewConst< arrayView2d< real64 const, 1 > > const & phaseMob, ElementViewConst< arrayView3d< real64 const > > const & dPhaseMob, ElementViewConst< arrayView3d< real64 const > > const & dCompFrac_dCompDens, @@ -285,10 +289,12 @@ struct HU2PhaseFlux seri, sesri, sei, + checkPhasePresenceInGravity, trans, dTrans_dPres, gravCoef, dCompFrac_dCompDens, + phaseVolFrac, phaseMassDens, dPhaseMassDens, pot_i, @@ -306,10 +312,12 @@ struct HU2PhaseFlux seri, sesri, sei, + checkPhasePresenceInGravity, trans, dTrans_dPres, gravCoef, dCompFrac_dCompDens, + phaseVolFrac, phaseMassDens, dPhaseMassDens, pot_j, @@ -410,7 +418,7 @@ struct HU2PhaseFlux static void computeTotalFlux( integer const & numPhase, const integer & hasCapPressure, - const integer & useNewGravity, + const integer & checkPhasePresenceInGravity, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], @@ -437,7 +445,7 @@ struct HU2PhaseFlux real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{}; PPUPhaseFlux::compute( numPhase, jp, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, @@ -479,10 +487,12 @@ struct HU2PhaseFlux localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], + integer const checkPhasePresenceInGravity, real64 const (&trans)[2], real64 const (&dTrans_dP)[2], ElementViewConst< arrayView1d< real64 const > > const & gravCoef, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, real64 & gravPot, @@ -492,50 +502,16 @@ struct HU2PhaseFlux // init UpwindHelpers::assignToZero( gravPot, dGravPot_dP, dGravPot_dC ); - // get average density TODO change after #3337 is merged - + // get average density real64 densMean{}; real64 dDensMean_dP[numFluxSupportPoints]{}; real64 dDensMean_dC[numFluxSupportPoints][numComp]{}; - - integer denom = 0; - for( localIndex i = 0; i < numFluxSupportPoints; ++i ) - { - localIndex const er = seri[i]; - localIndex const esr = sesri[i]; - localIndex const ei = sei[i]; - - // density - real64 const density = phaseMassDens[er][esr][ei][0][ip]; - real64 const dDens_dPres = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP]; - real64 dDens_dC[numComp]{}; - applyChainRule( numComp, - dCompFrac_dCompDens[er][esr][ei], - dPhaseMassDens[er][esr][ei][0][ip], - dDens_dC, - Deriv::dC ); - - // average density and derivatives - densMean += density; - dDensMean_dP[i] = dDens_dPres; - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dDensMean_dC[i][jc] = dDens_dC[jc]; - } - denom++; - } - if( denom > 1 ) - { - densMean /= denom; - for( integer i = 0; i < numFluxSupportPoints; ++i ) - { - dDensMean_dP[i] /= denom; - for( integer jc = 0; jc < numComp; ++jc ) - { - dDensMean_dC[i][jc] /= denom; - } - } - } + isothermalCompositionalMultiphaseFVMKernels::helpers:: + calculateMeanDensity( ip, seri, sesri, sei, + checkPhasePresenceInGravity, + phaseVolFrac, dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, + densMean, dDensMean_dP, dDensMean_dC ); // compute potential difference MPFA-style for( localIndex i = 0; i < numFluxSupportPoints; ++i ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp index 6127cd66050..559457abec2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp @@ -577,8 +577,7 @@ struct computePotentialViscous GEOS_UNUSED_PARAM( dPhaseCapPressure_dPhaseVolFrac ), real64 & pot, real64( &GEOS_UNUSED_PARAM( dPot_dPres ))[numFluxSupportPoints], - real64( &GEOS_UNUSED_PARAM( dPot_dComp ))[numFluxSupportPoints][numComp], - real64( &GEOS_UNUSED_PARAM( dProp_dComp ))[numComp] ) + real64( &GEOS_UNUSED_PARAM( dPot_dComp ))[numFluxSupportPoints][numComp] ) { pot = totFlux; //could be relevant for symmetry to include derivative @@ -617,27 +616,20 @@ struct computePotentialGravity GEOS_UNUSED_PARAM( dPhaseCapPressure_dPhaseVolFrac ), real64 & pot, real64 ( & dPot_dPres )[numFluxSupportPoints], - real64 (& dPot_dComp )[numFluxSupportPoints][numComp], - real64 ( & dProp_dComp )[numComp] ) + real64 ( & dPot_dComp )[numFluxSupportPoints][numComp] ) { + //init + assignToZero( pot, dPot_dPres, dPot_dComp ); + //working arrays real64 densMean{}; real64 dDensMean_dPres[numFluxSupportPoints]{}; real64 dDensMean_dComp[numFluxSupportPoints][numComp]{}; - - //init - assignToZero( pot, dPot_dPres, dPot_dComp ); - for( localIndex i = 0; i < numFluxSupportPoints; ++i ) - { - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dProp_dComp[jc] = 0.0; - } - } - - calculateMeanDensity( checkPhasePresenceInGravity, ip, seri, sesri, sei, + isothermalCompositionalMultiphaseFVMKernels::helpers:: + calculateMeanDensity( ip, seri, sesri, sei, + checkPhasePresenceInGravity, phaseVolFrac, dCompFrac_dCompDens, - phaseMassDens, dPhaseMassDens, dProp_dComp, + phaseMassDens, dPhaseMassDens, densMean, dDensMean_dPres, dDensMean_dComp ); // compute potential difference MPFA-style @@ -657,65 +649,6 @@ struct computePotentialGravity } } - template< localIndex numComp, localIndex numFluxSupportPoints > - GEOS_HOST_DEVICE - static void calculateMeanDensity( integer const checkPhasePresenceInGravity, - localIndex const ip, - localIndex const (&seri)[numFluxSupportPoints], - localIndex const (&sesri)[numFluxSupportPoints], - localIndex const (&sei)[numFluxSupportPoints], - ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, - 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, - real64 (& dProp_dComp)[numComp], - real64 & densMean, real64 (& dDensMean_dPres)[numFluxSupportPoints], real64 (& dDensMean_dComp)[numFluxSupportPoints][numComp] ) - { - integer denom = 0; - for( localIndex i = 0; i < numFluxSupportPoints; ++i ) - { - localIndex const er = seri[i]; - localIndex const esr = sesri[i]; - localIndex const ei = sei[i]; - - bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( checkPhasePresenceInGravity && !phaseExists ) - { - continue; - } - - // density - real64 const density = phaseMassDens[er][esr][ei][0][ip]; - real64 const dDens_dPres = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP]; - - applyChainRule( numComp, - dCompFrac_dCompDens[er][esr][ei], - dPhaseMassDens[er][esr][ei][0][ip], - dProp_dComp, - Deriv::dC ); - - // average density and derivatives - densMean += density; - dDensMean_dPres[i] = dDens_dPres; - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dDensMean_dComp[i][jc] = dProp_dComp[jc]; - } - denom++; - } - if( denom > 1 ) - { - densMean /= denom; - for( localIndex i = 0; i < numFluxSupportPoints; ++i ) - { - dDensMean_dPres[i] /= denom; - for( integer jc = 0; jc < numComp; ++jc ) - { - dDensMean_dComp[i][jc] /= denom; - } - } - } - } }; /*! @copydoc computePotential @@ -748,8 +681,7 @@ struct computePotentialCapillary ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, real64 & pot, real64 ( & dPot_dPres)[numFluxSupportPoints], - real64 (& dPot_dComp)[numFluxSupportPoints][numComp], - real64( &GEOS_UNUSED_PARAM( dProp_dComp ))[numComp] ) + real64 (& dPot_dComp)[numFluxSupportPoints][numComp] ) { //init assignToZero( pot, dPot_dPres, dPot_dComp ); @@ -822,7 +754,6 @@ static void computePotentialFluxesGravity( localIndex const numPhase, real64 pot{}; real64 dPot_dP[numFluxSupportPoints]{}; real64 dPot_dC[numFluxSupportPoints][numComp]{}; - real64 dProp_dC[numComp]{}; // UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, @@ -844,8 +775,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, dPhaseCapPressure_dPhaseVolFrac, pot, dPot_dP, - dPot_dC, - dProp_dC ); + dPot_dC ); // and the fractional flow for gravitational part as \lambda_i^{up}/\sum_{numPhase}(\lambda_k^{up}) with up decided upon // the Upwind strategy @@ -886,7 +816,6 @@ static void computePotentialFluxesGravity( localIndex const numPhase, real64 potOther{}; real64 dPotOther_dP[numFluxSupportPoints]{}; real64 dPotOther_dC[numFluxSupportPoints][numComp]{}; - real64 dPropOther_dC[numComp]{}; //Fetch pot for phase j!=i defined as \rho_j g dz/dx UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, @@ -908,8 +837,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, dPhaseCapPressure_dPhaseVolFrac, potOther, dPotOther_dP, - dPotOther_dC, - dPropOther_dC ); + dPotOther_dC ); //Eventually get the mobility of the second phase localIndex k_up_o = -1; @@ -999,7 +927,6 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, real64 pot{}; real64 dPot_dP[numFluxSupportPoints]{}; real64 dPot_dC[numFluxSupportPoints][numComp]{}; - real64 dProp_dC[numComp]{}; computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase, ip, @@ -1018,8 +945,7 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, dPhaseCapPressure_dPhaseVolFrac, pot, dPot_dP, - dPot_dC, - dProp_dC ); + dPot_dC ); // and the fractional flow for gravitational part as \lambda_i^{up}/\sum_{numPhase}(\lambda_k^{up}) with up decided upon // the Upwind strategy @@ -1057,7 +983,6 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, real64 potOther{}; real64 dPotOther_dP[numFluxSupportPoints]{}; real64 dPotOther_dC[numFluxSupportPoints][numComp]{}; - real64 dPropOther_dC[numComp]{}; //Fetch pot for phase j!=i defined as \rho_j g dz/dx computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase, @@ -1077,8 +1002,7 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, dPhaseCapPressure_dPhaseVolFrac, potOther, dPotOther_dP, - dPotOther_dC, - dPropOther_dC ); + dPotOther_dC ); //Eventually get the mobility of the second phase localIndex k_up_o = -1; @@ -1336,9 +1260,8 @@ class UpwindScheme real64 pot{}; real64 pot_dP[numFluxSupportPoints]{}; real64 pot_dC[numFluxSupportPoints][numComp]{}; - real64 dProp_dC[numComp]{}; - fn( ip, pot, pot_dP, pot_dC, dProp_dC ); + fn( ip, pot, pot_dP, pot_dC ); localIndex const k_up = 0; localIndex const k_dw = 1; @@ -1359,9 +1282,8 @@ class UpwindScheme real64 potOther{}; real64 potOther_dP[numFluxSupportPoints]{}; real64 potOther_dC[numFluxSupportPoints][numComp]{}; - real64 dPropOther_dC[numComp]{}; - fn( jp, potOther, potOther_dP, potOther_dC, dPropOther_dC ); + fn( jp, potOther, potOther_dP, potOther_dC ); real64 const mob_up = phaseMob[er_up][esr_up][ei_up][jp]; real64 const mob_dw = phaseMob[er_dw][esr_dw][ei_dw][jp]; @@ -1404,12 +1326,10 @@ class HybridUpwind : public UpwindScheme ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, integer const GEOS_UNUSED_PARAM( hasCapPressure ), - real64 & potential - ) + real64 & potential ) { real64 dPot_dP[numFluxSupportPoints]{}; real64 dPot_dC[numFluxSupportPoints][numComp]{}; - real64 dProp_dC[numComp]{}; UpwindHelpers::computePotentialViscous::compute< numComp, numFluxSupportPoints >( numPhase, ip, @@ -1428,8 +1348,7 @@ class HybridUpwind : public UpwindScheme dPhaseCapPressure_dPhaseVolFrac, potential, dPot_dP, - dPot_dC, - dProp_dC ); + dPot_dC ); } template< localIndex numComp, localIndex numFluxSupportPoints > @@ -1467,8 +1386,7 @@ class HybridUpwind : public UpwindScheme [&]( localIndex ipp, real64 & potential_, real64 (& dPotential_dP_)[numFluxSupportPoints], - real64 (& dPotential_dC_)[numFluxSupportPoints][numComp], - real64 (& dProp_dC)[numComp] ) { + real64 (& dPotential_dC_)[numFluxSupportPoints][numComp] ) { UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, @@ -1490,8 +1408,7 @@ class HybridUpwind : public UpwindScheme dPhaseCapPressure_dPhaseVolFrac, potential_, dPotential_dP_, - dPotential_dC_, - dProp_dC ); + dPotential_dC_ ); } ); } @@ -1532,8 +1449,8 @@ class HybridUpwind : public UpwindScheme [&]( localIndex ipp, real64 & potential_, real64 (& dPotential_dP_)[numFluxSupportPoints], - real64 (& dPotential_dC_)[numFluxSupportPoints][numComp], - real64 (& dProp_dC)[numComp] ) { + real64 (& dPotential_dC_)[numFluxSupportPoints][numComp] ) + { UpwindHelpers::computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase, @@ -1553,8 +1470,7 @@ class HybridUpwind : public UpwindScheme dPhaseCapPressure_dPhaseVolFrac, potential_, dPotential_dP_, - dPotential_dC_, - dProp_dC ); + dPotential_dC_ ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp index b783b061712..067090f1ed6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp @@ -79,17 +79,17 @@ struct PotGrad } } - // create local work arrays - real64 densMean = 0.0; - real64 dDensMean_dP[numFluxSupportPoints]{}; - real64 dDensMean_dC[numFluxSupportPoints][numComp]{}; - real64 presGrad = 0.0; real64 gravHead = 0.0; real64 dCapPressure_dC[numComp]{}; - calculateMeanDensity( checkPhasePresenceInGravity, ip, - seri, sesri, sei, + // create local work arrays + real64 densMean = 0.0; + real64 dDensMean_dP[numFluxSupportPoints]{}; + real64 dDensMean_dC[numFluxSupportPoints][numComp]{}; + isothermalCompositionalMultiphaseFVMKernels::helpers:: + calculateMeanDensity( ip, seri, sesri, sei, + checkPhasePresenceInGravity, phaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, densMean, dDensMean_dP, dDensMean_dC ); @@ -158,68 +158,6 @@ struct PotGrad } - template< integer numComp, integer numFluxSupportPoints > - GEOS_HOST_DEVICE - static void - calculateMeanDensity( integer const checkPhasePresenceInGravity, - integer const ip, - localIndex const ( &seri )[numFluxSupportPoints], - localIndex const ( &sesri )[numFluxSupportPoints], - localIndex const ( &sei )[numFluxSupportPoints], - ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, - 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, - real64 & densMean, real64 ( & dDensMean_dP)[numFluxSupportPoints], real64 ( & dDensMean_dC )[numFluxSupportPoints][numComp] ) - { - real64 dDens_dC[numComp]{}; - - integer denom = 0; - for( integer i = 0; i < numFluxSupportPoints; ++i ) - { - localIndex const er = seri[i]; - localIndex const esr = sesri[i]; - localIndex const ei = sei[i]; - - bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( checkPhasePresenceInGravity && !phaseExists ) - { - continue; - } - - // density - real64 const density = phaseMassDens[er][esr][ei][0][ip]; - real64 const dDens_dP = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP]; - - applyChainRule( numComp, - dCompFrac_dCompDens[er][esr][ei], - dPhaseMassDens[er][esr][ei][0][ip], - dDens_dC, - Deriv::dC ); - - // average density and derivatives - densMean += density; - dDensMean_dP[i] = dDens_dP; - for( integer jc = 0; jc < numComp; ++jc ) - { - dDensMean_dC[i][jc] = dDens_dC[jc]; - } - denom++; - } - if( denom > 1 ) - { - densMean /= denom; - for( integer i = 0; i < numFluxSupportPoints; ++i ) - { - dDensMean_dP[i] /= denom; - for( integer jc = 0; jc < numComp; ++jc ) - { - dDensMean_dC[i][jc] /= denom; - } - } - } - } - }; } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities