From 55761ee3dfe916a4f7fff13a9802d9cd8186f99b Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 17 Jul 2024 16:31:13 -0700 Subject: [PATCH 1/7] initial setup for 2 phase PPU + interface condition flux --- .../finiteVolume/FluxApproximationBase.hpp | 8 +- ...positionalMultiphaseFVMKernelUtilities.hpp | 75 +++++++++++++++++++ ...ermalCompositionalMultiphaseFVMKernels.hpp | 33 +++++++- 3 files changed, 111 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/finiteVolume/FluxApproximationBase.hpp b/src/coreComponents/finiteVolume/FluxApproximationBase.hpp index 3830cf52cc1..85ef747e57b 100644 --- a/src/coreComponents/finiteVolume/FluxApproximationBase.hpp +++ b/src/coreComponents/finiteVolume/FluxApproximationBase.hpp @@ -38,7 +38,8 @@ enum class UpwindingScheme : integer { PPU, ///< PPU upwinding C1PPU, ///< C1-PPU upwinding from https://doi.org/10.1016/j.advwatres.2017.07.028 - IHU ///< IHU as in https://link.springer.com/content/pdf/10.1007/s10596-019-09835-6.pdf + IHU, ///< IHU as in https://link.springer.com/content/pdf/10.1007/s10596-019-09835-6.pdf + PPUIC }; /** @@ -47,7 +48,8 @@ enum class UpwindingScheme : integer ENUM_STRINGS( UpwindingScheme, "PPU", "C1PPU", - "IHU" ); + "IHU", + "PPUIC" ); /** * @struct UpwindingParameters @@ -56,7 +58,7 @@ ENUM_STRINGS( UpwindingScheme, */ struct UpwindingParameters { - /// PPU or C1-PPU or IHU + /// PPU or C1-PPU or IHU or PPU+interface conditions UpwindingScheme upwindingScheme; /// C1-PPU smoothing tolerance diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp index 9824d8a17e4..cab56544360 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp @@ -562,6 +562,81 @@ struct C1PPUPhaseFlux +struct PPUICPhaseFlux +{ + /** + * @brief Form the PhasePotentialUpwind from pressure gradient and gravitational head + * @tparam numComp number of components + * @tparam numFluxSupportPoints number of flux support points + * @param numPhase number of phases + * @param ip phase index + * @param hasCapPressure flag indicating if there is capillary pressure + * @param seri arraySlice of the stencil-implied element region index + * @param sesri arraySlice of the stencil-implied element subregion index + * @param sei arraySlice of the stencil-implied element index + * @param trans transmissibility at the connection + * @param dTrans_dPres derivative of transmissibility wrt pressure + * @param pres pressure + * @param gravCoef gravitational coefficient + * @param phaseMob phase mobility + * @param dPhaseMob derivative of phase mobility wrt pressure, temperature, comp density + * @param dPhaseVolFrac derivative of phase volume fraction wrt pressure, temperature, comp density + * @param dCompFrac_dCompDens derivative of component fraction wrt component density + * @param phaseMassDens phase mass density + * @param dPhaseMassDens derivative of phase mass density wrt pressure, temperature, comp fraction + * @param phaseCapPressure phase capillary pressure + * @param dPhaseCapPressure_dPhaseVolFrac derivative of phase capillary pressure wrt phase volume fraction + * @param k_up uptream index for this phase + * @param potGrad potential gradient for this phase + * @param phaseFlux phase flux + * @param dPhaseFlux_dP derivative of phase flux wrt pressure + * @param dPhaseFlux_dC derivative of phase flux wrt comp density + */ + template< integer numComp, integer numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + compute( integer const numPhase, + integer const ip, + integer const hasCapPressure, + 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< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, + ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, + 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, + localIndex & k_up, + real64 & potGrad, + real64 ( &phaseFlux ), + real64 ( & dPhaseFlux_dP )[numFluxSupportPoints], + real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp], + real64 ( & compFlux )[numComp], + real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp], + real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] ) + { + // Interface condition solver goes here + // Can use PPU phaseflux for reference + // TODO Ryan- put in checks to ensure using correct constitutive law (maybe on kernel launch?) + // Some of the helpers below might also be useful, looks like whoever implemented IHU used them + // + // A good step 1 would be to make sure all the constitutive info we are gonna need is here, if not we can pass more by using accessors in FaceBasedAssemblyKernel (Ryan will help) + } +}; + + + + + /************************* HELPERS ******************/ namespace UpwindHelpers { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp index 7efdb9aa8d3..8d91879bee8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp @@ -66,9 +66,10 @@ enum class FaceBasedAssemblyKernelFlags /// Flag indicating whether C1-PPU is used or not C1PPU = 1 << 2, // 4 /// Flag indicating whether IHU is used or not - IHU = 1 << 3 // 8 + IHU = 1 << 3, // 8 + /// Flag indicating whether interface conditions are used (PPU only for now) + PPUIC = 1 << 4 // 16 /// Add more flags like that if needed: - // Flag5 = 1 << 4, // 16 // Flag6 = 1 << 5, // 32 // Flag7 = 1 << 6, // 64 // Flag8 = 1 << 7 //128 @@ -673,6 +674,32 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase dCompFlux_dP, dCompFlux_dC ); } + else if( m_kernelFlags.isSet( FaceBasedAssemblyKernelFlags::PPUIC ) ) + { + isothermalCompositionalMultiphaseFVMKernelUtilities::PPUICPhaseFlux::compute< numComp, numFluxSupportPoints > + ( m_numPhases, + ip, + m_kernelFlags.isSet( FaceBasedAssemblyKernelFlags::CapPressure ), + seri, sesri, sei, + trans, + dTrans_dPres, + m_pres, + m_gravCoef, + m_phaseMob, m_dPhaseMob, + m_dPhaseVolFrac, + m_phaseCompFrac, m_dPhaseCompFrac, + m_dCompFrac_dCompDens, + m_phaseMassDens, m_dPhaseMassDens, + m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac, + k_up, + potGrad, + phaseFlux, + dPhaseFlux_dP, + dPhaseFlux_dC, + compFlux, + dCompFlux_dP, + dCompFlux_dC ); + } else { isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux::compute< numComp, numFluxSupportPoints > @@ -902,6 +929,8 @@ class FaceBasedAssemblyKernelFactory kernelFlags.set( FaceBasedAssemblyKernelFlags::C1PPU ); else if( upwindingParams.upwindingScheme == UpwindingScheme::IHU ) kernelFlags.set( FaceBasedAssemblyKernelFlags::IHU ); + else if( upwindingParams.upwindingScheme == UpwindingScheme::PPUIC) + kernelFlags.set( FaceBasedAssemblyKernelFlags::PPUIC ); using kernelType = FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; From 01267d068e543ab9fff6fbe3e9aee72098dc5ce1 Mon Sep 17 00:00:00 2001 From: Ammara-14 Date: Wed, 31 Jul 2024 15:01:21 -0700 Subject: [PATCH 2/7] testing git --- ...hermalCompositionalMultiphaseFVMKernelUtilities.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp index 6f699009df6..d8008828328 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp @@ -631,6 +631,16 @@ struct PPUICPhaseFlux // Some of the helpers below might also be useful, looks like whoever implemented IHU used them // // A good step 1 would be to make sure all the constitutive info we are gonna need is here, if not we can pass more by using accessors in FaceBasedAssemblyKernel (Ryan will help) + // + + + + // Step 1: convert from density to component fraction + // + // step 2: loval solver + // + // step 3: global solver fluix and jacobian + } }; From e2b4389ef438f842a0406e53fa049d455b2716bb Mon Sep 17 00:00:00 2001 From: Ammara-14 Date: Wed, 31 Jul 2024 15:50:30 -0700 Subject: [PATCH 3/7] another commit --- .../IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp index d8008828328..49cf803577f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp @@ -640,6 +640,9 @@ struct PPUICPhaseFlux // step 2: loval solver // // step 3: global solver fluix and jacobian + // + // Another change for tutorial + // } }; From 8ecf871048e3c3bf8640414d69ab7e33a3bec31f Mon Sep 17 00:00:00 2001 From: Ammara-14 Date: Thu, 1 Aug 2024 10:52:22 -0700 Subject: [PATCH 4/7] Added Total Flux Computation --- ...positionalMultiphaseFVMKernelUtilities.hpp | 89 +++++++++++++++++-- 1 file changed, 83 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp index 49cf803577f..9ff051c2b47 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp @@ -627,22 +627,99 @@ struct PPUICPhaseFlux { // Interface condition solver goes here // Can use PPU phaseflux for reference + + // TODO Ryan- put in checks to ensure using correct constitutive law (maybe on kernel launch?) // Some of the helpers below might also be useful, looks like whoever implemented IHU used them // // A good step 1 would be to make sure all the constitutive info we are gonna need is here, if not we can pass more by using accessors in FaceBasedAssemblyKernel (Ryan will help) // + // Step 1: convert from density to component fraction + // Ammar suggest to use the provided derivative 1/dCompFrac_dCompDens since the algebraic expansions show that this is equal to dCompDens_dCompFrac + // step 2: loval solver + + // First we calculate the total flux using PPU flux: + + real64 dPresGrad_dP[numFluxSupportPoints]{}; + real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; + real64 dGravHead_dP[numFluxSupportPoints]{}; + real64 dGravHead_dC[numFluxSupportPoints][numComp]{}; + PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, seri, sesri, sei, trans, dTrans_dPres, pres, + gravCoef, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP, + dPresGrad_dC, dGravHead_dP, dGravHead_dC ); + + // *** upwinding *** + + // choose upstream cell + k_up = (potGrad >= 0) ? 0 : 1; + + localIndex const er_up = seri[k_up]; + localIndex const esr_up = sesri[k_up]; + localIndex const ei_up = sei[k_up]; + + real64 const mobility = phaseMob[er_up][esr_up][ei_up][ip]; + + //loop over all phases to form total velocity + real64 totFlux{}; + real64 dTotFlux_dP[numFluxSupportPoints]{}; + real64 dTotFlux_dC[numFluxSupportPoints][numComp]{}; + + //store totMob upwinded by PPU for later schemes + real64 totMob{}; + real64 dTotMob_dP[numFluxSupportPoints]{}; + real64 dTotMob_dC[numFluxSupportPoints][numComp]{}; + localIndex k_up_ppu = -1; + + //unelegant but need dummy when forming PPU total velocity + real64 dummy[numComp]; + real64 dDummy_dP[numFluxSupportPoints][numComp]; + real64 dDummy_dC[numFluxSupportPoints][numComp][numComp]; + + + for( integer jp = 0; jp < numPhase; ++jp ) + { + PPUPhaseFlux::compute( numPhase, jp, hasCapPressure, + seri, sesri, sei, + trans, dTrans_dPres, + pres, gravCoef, + phaseMob, dPhaseMob, + dPhaseVolFrac, + phaseCompFrac, dPhaseCompFrac, + dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + k_up_ppu, potGrad, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, + dummy, dDummy_dP, dDummy_dC ); + + totFlux += phaseFlux; + + phaseFlux = 0.; + for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) + { + dTotFlux_dP[ke] += dPhaseFlux_dP[ke]; + totMob += phaseMob[seri[ke]][sesri[ke]][sei[ke]][jp]; + dTotMob_dP[ke] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dP]; + dPhaseFlux_dP[ke] = 0.; + + for( localIndex jc = 0; jc < numComp; ++jc ) + { + dTotFlux_dC[ke][jc] += dPhaseFlux_dC[ke][jc]; + dTotMob_dC[ke][jc] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dC + jc]; + dPhaseFlux_dC[ke][jc] = 0.; + } + } + + + } - // Step 1: convert from density to component fraction - // - // step 2: loval solver - // - // step 3: global solver fluix and jacobian // - // Another change for tutorial + // step 3: global solver flux and jacobian // + } }; From 3433c2497eba31a6df8b91c849ab29c77573067e Mon Sep 17 00:00:00 2001 From: Ammara-14 Date: Fri, 2 Aug 2024 11:25:02 -0700 Subject: [PATCH 5/7] Added initialization of local solver - it shows linking error --- ...positionalMultiphaseFVMKernelUtilities.hpp | 215 +++++++++++++++++- 1 file changed, 203 insertions(+), 12 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp index 9ff051c2b47..7fb35e6d600 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp @@ -19,6 +19,9 @@ #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_ISOTHERMALCOMPOSITIONALMULTIPHASEFVMKERNELUTILITIES_HPP_ #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_ISOTHERMALCOMPOSITIONALMULTIPHASEFVMKERNELUTILITIES_HPP_ +#define Pe_max 5000000.0 +#define Pe_min 4000000.0 +#define eps2 0.05 #include "common/DataLayouts.hpp" #include "common/DataTypes.hpp" @@ -33,6 +36,136 @@ namespace geos namespace isothermalCompositionalMultiphaseFVMKernelUtilities { +real64 computePCalpha ( real64 S ) { + real64 Pe = Pe_max; + return (Pe * (1.0 - S)); +} + +real64 computePCalphaInv ( real64 Pc ) { + real64 Pe = Pe_max; + return (1.0 - Pc / Pe); +} + +real64 computePCalpha_tilde ( real64 S ) { + real64 Pe = Pe_max; + real64 Pc_tilde = Pe * (1.0 - S); + if (S < (1.0 - Pe_min / Pe_max) ) { + Pc_tilde = Pe_min; + } + return Pc_tilde; +} + +real64 computePCbeta ( real64 S ) { + real64 Pe = Pe_min; + return (Pe * (1.0 - S)); +} + +real64 computePCbetaInv ( real64 Pc ) { + real64 Pe = Pe_min; + return (1.0 - Pc / Pe); +} + +real64 computeMobilityL ( real64 S ) { + real64 mu = 0.0004; + return (S*S / mu); +} + +real64 computeMobilityV ( real64 S ) { + real64 mu = 0.00002; + return (std::pow((1.0 - S),2) / mu); +} + +real64 computedPCdSalpha ( real64 S ) { + GEOS_UNUSED_VAR(S); + real64 Pe = Pe_max; + return (Pe * (- 1.0)); +} + +real64 computedPCdSalpha_tilde ( real64 S ) { + real64 Pe = Pe_max; + real64 dPcdS_tilde = Pe * (- 1.0); + if (S < (1.0 - Pe_min / Pe_max)) { + dPcdS_tilde = 0.0; + } + return dPcdS_tilde; +} + +real64 computedPCdSbeta ( real64 S ) { + GEOS_UNUSED_VAR(S); + real64 Pe = Pe_min; + return (Pe * (- 1.0)); +} + +real64 computedPCbetaInv ( real64 Pc ) { + GEOS_UNUSED_VAR(Pc); + real64 Pe = Pe_min; + return (- 1.0 / Pe); +} + +real64 computedMobilitydSL ( real64 S ) { + // GEOS_UNUSED_VAR(S); + real64 mu = 0.0004; + return (2.0 * S / mu); +} + +real64 computedMobilitydSV ( real64 S ) { + // GEOS_UNUSED_VAR(S); + real64 mu = 0.00002; + return (- 2.0 * (1.0 - S) / mu); +} + +real64 computeVli ( real64 Z ) { + real64 Vli = -0.4291 * std::pow(Z,3) + 1.4209 * Z * Z - 1.9911 * Z + 1; + return Vli; +} + +real64 computedVlidZ ( real64 Z ) { + real64 dVli = -0.4291 * 3.0 * std::pow(Z,2) + 1.4209 * 2.0 * Z - 1.9911; + return dVli; +} + +real64 computeVlj ( real64 Z ) { + real64 Vlj = -0.4541 * std::pow(Z,3) + 1.4637 * Z * Z - 2.0092 * Z + 1; + return Vlj; +} + +real64 computedVljdZ ( real64 Z ) { + real64 dVlj = -0.4541 * 3.0 * std::pow(Z,2) + 1.4637 * 2.0 * Z - 2.0092; + return dVlj; +} + +real64 computerhoTi ( real64 Z ) { + real64 rhoTi = -359.23 * std::pow(Z,3) + 921.51 * Z * Z - 911.94 * Z + 1017; + return rhoTi; +} + +real64 computedrhoTidZ ( real64 Z ) { + real64 drhoTi = -359.23 * 3.0 * std::pow(Z,2) + 921.51 * 2.0 * Z - 911.94; + return drhoTi; +} + +real64 computerhoTj ( real64 Z ) { + real64 rhoTj = -711.67 * std::pow(Z,3) + 1657.1 * Z * Z - 1419.5 * Z + 1009.6; + return rhoTj; +} + +real64 computedrhoTjdZ ( real64 Z ) { + real64 drhoTj = -711.67 * 3.0 * std::pow(Z,2) + 1657.1 * 2.0 * Z - 1419.5; + return drhoTj; +} + +real64 safeUpdate ( real64 S1, real64 S0) { +real64 S_kink1 = 1.0; real64 S_kink2 = 0.0; +if ( S0 < S_kink1 && S1 > (S_kink1 + eps2)) { + S1 = S_kink1 - eps2; +} +if ( S0 > (S_kink2 + eps2) && S1 < (S_kink2 + eps2)) { + S1 = S_kink2 + eps2; +} +return S1; +} + + // TODO make input parameter static constexpr real64 epsC1PPU = 5000; @@ -635,9 +768,9 @@ struct PPUICPhaseFlux // A good step 1 would be to make sure all the constitutive info we are gonna need is here, if not we can pass more by using accessors in FaceBasedAssemblyKernel (Ryan will help) // // Step 1: convert from density to component fraction - // Ammar suggest to use the provided derivative 1/dCompFrac_dCompDens since the algebraic expansions show that this is equal to dCompDens_dCompFrac + // Ammar- suggest to use the provided derivative to multiply by 1/dCompFrac_dCompDens since the algebraic expansions show that this is equal to dCompDens_dCompFrac - // step 2: loval solver + // step 2: local solver // First we calculate the total flux using PPU flux: @@ -666,12 +799,6 @@ struct PPUICPhaseFlux real64 dTotFlux_dP[numFluxSupportPoints]{}; real64 dTotFlux_dC[numFluxSupportPoints][numComp]{}; - //store totMob upwinded by PPU for later schemes - real64 totMob{}; - real64 dTotMob_dP[numFluxSupportPoints]{}; - real64 dTotMob_dC[numFluxSupportPoints][numComp]{}; - localIndex k_up_ppu = -1; - //unelegant but need dummy when forming PPU total velocity real64 dummy[numComp]; real64 dDummy_dP[numFluxSupportPoints][numComp]; @@ -690,7 +817,7 @@ struct PPUICPhaseFlux dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - k_up_ppu, potGrad, + k_up, potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dummy, dDummy_dP, dDummy_dC ); @@ -700,14 +827,11 @@ struct PPUICPhaseFlux for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) { dTotFlux_dP[ke] += dPhaseFlux_dP[ke]; - totMob += phaseMob[seri[ke]][sesri[ke]][sei[ke]][jp]; - dTotMob_dP[ke] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dP]; dPhaseFlux_dP[ke] = 0.; for( localIndex jc = 0; jc < numComp; ++jc ) { dTotFlux_dC[ke][jc] += dPhaseFlux_dC[ke][jc]; - dTotMob_dC[ke][jc] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dC + jc]; dPhaseFlux_dC[ke][jc] = 0.; } } @@ -715,6 +839,73 @@ struct PPUICPhaseFlux } + //////////// + + // calculate quantities on primary connected cells + int i = 0; + localIndex const er_i = seri[i]; + localIndex const esr_i = sesri[i]; + localIndex const ei_i = sei[i]; + + localIndex const er_j = seri[i+1]; + localIndex const esr_j = sesri[i+1]; + localIndex const ei_j = sei[i+1]; + + // density + real64 const rho_l_i = phaseMassDens[er_i][esr_i][ei_i][0][ip]; + real64 const rho_l_j = phaseMassDens[er_j][esr_j][ei_j][0][ip]; + + real64 const pi = pres[er_i][esr_i][ei_i]; + real64 const pj = pres[er_j][esr_j][ei_j]; + + real64 const Ti = trans[i]; + real64 const Tj = trans[i+1]; + + // Newton Loop + // missing Zi and Zj, now hard-coded: + real64 const Zi = 0.1; real64 const Zj = 0.5; + // first guess for Z^alpha + real64 Z_alpha = (Zi + Zj) / 2.0; + + if (Zi > 1.0 - eps2) { + Z_alpha = 1.0 - eps2; + } + if (Zj > 1.0 - eps2) { + Z_alpha = 1.0 - eps2; + } + + real64 rho_t_i = computerhoTi ( Zi ); + real64 v_l_i = computeVli ( Zi ); + + real64 Si = v_l_i * rho_t_i / rho_l_i; + if (Si < 0.0) { + Si = 0.0; + } + if (Si > 1.0) { + Si = 1.0; + } + real64 pc_i = computePCalpha (Si); + real64 dpcdSi = computedPCdSalpha (Si); + + real64 rho_t_j = computerhoTj ( Zj ); + real64 v_l_j = computeVlj ( Zj ); + real64 Sj = v_l_j * rho_t_j / rho_l_j; + if (Sj < 0.0) { + Sj = 0.0; + } + if (Sj > 1.0) { + Sj = 1.0; + } + real64 pc_j = computePCbeta (Sj); + real64 dpcdSj = computedPCdSbeta (Sj); + + real64 Ut = totFlux; + + real64 tol = 1e-8; int max_iter = 50; + +// Just a test: + + phaseFlux = Z_alpha * 2.0; // // step 3: global solver flux and jacobian From bab3938f23cc5765575b2bb852d2944c6051aa8c Mon Sep 17 00:00:00 2001 From: Ammara-14 Date: Tue, 6 Aug 2024 15:35:46 -0700 Subject: [PATCH 6/7] Added initialization of local solver and fixed linking error --- ...positionalMultiphaseFVMKernelUtilities.hpp | 242 +++++++++--------- 1 file changed, 122 insertions(+), 120 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp index 7fb35e6d600..446fa51675e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp @@ -35,135 +35,137 @@ namespace geos namespace isothermalCompositionalMultiphaseFVMKernelUtilities { - +template< typename VIEWTYPE > real64 computePCalpha ( real64 S ) { real64 Pe = Pe_max; return (Pe * (1.0 - S)); } -real64 computePCalphaInv ( real64 Pc ) { - real64 Pe = Pe_max; - return (1.0 - Pc / Pe); -} - -real64 computePCalpha_tilde ( real64 S ) { - real64 Pe = Pe_max; - real64 Pc_tilde = Pe * (1.0 - S); - if (S < (1.0 - Pe_min / Pe_max) ) { - Pc_tilde = Pe_min; - } - return Pc_tilde; -} - -real64 computePCbeta ( real64 S ) { - real64 Pe = Pe_min; - return (Pe * (1.0 - S)); -} - -real64 computePCbetaInv ( real64 Pc ) { - real64 Pe = Pe_min; - return (1.0 - Pc / Pe); -} - -real64 computeMobilityL ( real64 S ) { - real64 mu = 0.0004; - return (S*S / mu); -} - -real64 computeMobilityV ( real64 S ) { - real64 mu = 0.00002; - return (std::pow((1.0 - S),2) / mu); -} - +// real64 computePCalphaInv ( real64 Pc ) { +// real64 Pe = Pe_max; +// return (1.0 - Pc / Pe); +// } + +// real64 computePCalpha_tilde ( real64 S ) { +// real64 Pe = Pe_max; +// real64 Pc_tilde = Pe * (1.0 - S); +// if (S < (1.0 - Pe_min / Pe_max) ) { +// Pc_tilde = Pe_min; +// } +// return Pc_tilde; +// } + +// real64 computePCbeta ( real64 S ) { +// real64 Pe = Pe_min; +// return (Pe * (1.0 - S)); +// } + +// real64 computePCbetaInv ( real64 Pc ) { +// real64 Pe = Pe_min; +// return (1.0 - Pc / Pe); +// } + +// real64 computeMobilityL ( real64 S ) { +// real64 mu = 0.0004; +// return (S*S / mu); +// } + +// real64 computeMobilityV ( real64 S ) { +// real64 mu = 0.00002; +// return (std::pow((1.0 - S),2) / mu); +// } +template< typename VIEWTYPE > real64 computedPCdSalpha ( real64 S ) { GEOS_UNUSED_VAR(S); real64 Pe = Pe_max; return (Pe * (- 1.0)); } -real64 computedPCdSalpha_tilde ( real64 S ) { - real64 Pe = Pe_max; - real64 dPcdS_tilde = Pe * (- 1.0); - if (S < (1.0 - Pe_min / Pe_max)) { - dPcdS_tilde = 0.0; - } - return dPcdS_tilde; -} - -real64 computedPCdSbeta ( real64 S ) { - GEOS_UNUSED_VAR(S); - real64 Pe = Pe_min; - return (Pe * (- 1.0)); -} - -real64 computedPCbetaInv ( real64 Pc ) { - GEOS_UNUSED_VAR(Pc); - real64 Pe = Pe_min; - return (- 1.0 / Pe); -} - -real64 computedMobilitydSL ( real64 S ) { - // GEOS_UNUSED_VAR(S); - real64 mu = 0.0004; - return (2.0 * S / mu); -} - -real64 computedMobilitydSV ( real64 S ) { - // GEOS_UNUSED_VAR(S); - real64 mu = 0.00002; - return (- 2.0 * (1.0 - S) / mu); -} +// real64 computedPCdSalpha_tilde ( real64 S ) { +// real64 Pe = Pe_max; +// real64 dPcdS_tilde = Pe * (- 1.0); +// if (S < (1.0 - Pe_min / Pe_max)) { +// dPcdS_tilde = 0.0; +// } +// return dPcdS_tilde; +// } + +// real64 computedPCdSbeta ( real64 S ) { +// GEOS_UNUSED_VAR(S); +// real64 Pe = Pe_min; +// return (Pe * (- 1.0)); +// } + +// real64 computedPCbetaInv ( real64 Pc ) { +// GEOS_UNUSED_VAR(Pc); +// real64 Pe = Pe_min; +// return (- 1.0 / Pe); +// } + +// real64 computedMobilitydSL ( real64 S ) { +// // GEOS_UNUSED_VAR(S); +// real64 mu = 0.0004; +// return (2.0 * S / mu); +// } + +// real64 computedMobilitydSV ( real64 S ) { +// // GEOS_UNUSED_VAR(S); +// real64 mu = 0.00002; +// return (- 2.0 * (1.0 - S) / mu); +// } +template< typename VIEWTYPE > real64 computeVli ( real64 Z ) { real64 Vli = -0.4291 * std::pow(Z,3) + 1.4209 * Z * Z - 1.9911 * Z + 1; return Vli; } -real64 computedVlidZ ( real64 Z ) { - real64 dVli = -0.4291 * 3.0 * std::pow(Z,2) + 1.4209 * 2.0 * Z - 1.9911; - return dVli; -} +// real64 computedVlidZ ( real64 Z ) { +// real64 dVli = -0.4291 * 3.0 * std::pow(Z,2) + 1.4209 * 2.0 * Z - 1.9911; +// return dVli; +// } -real64 computeVlj ( real64 Z ) { - real64 Vlj = -0.4541 * std::pow(Z,3) + 1.4637 * Z * Z - 2.0092 * Z + 1; - return Vlj; -} +// real64 computeVlj ( real64 Z ) { +// real64 Vlj = -0.4541 * std::pow(Z,3) + 1.4637 * Z * Z - 2.0092 * Z + 1; +// return Vlj; +// } -real64 computedVljdZ ( real64 Z ) { - real64 dVlj = -0.4541 * 3.0 * std::pow(Z,2) + 1.4637 * 2.0 * Z - 2.0092; - return dVlj; -} +// real64 computedVljdZ ( real64 Z ) { +// real64 dVlj = -0.4541 * 3.0 * std::pow(Z,2) + 1.4637 * 2.0 * Z - 2.0092; +// return dVlj; +// } +template< typename VIEWTYPE > real64 computerhoTi ( real64 Z ) { real64 rhoTi = -359.23 * std::pow(Z,3) + 921.51 * Z * Z - 911.94 * Z + 1017; return rhoTi; } -real64 computedrhoTidZ ( real64 Z ) { - real64 drhoTi = -359.23 * 3.0 * std::pow(Z,2) + 921.51 * 2.0 * Z - 911.94; - return drhoTi; -} - -real64 computerhoTj ( real64 Z ) { - real64 rhoTj = -711.67 * std::pow(Z,3) + 1657.1 * Z * Z - 1419.5 * Z + 1009.6; - return rhoTj; -} - -real64 computedrhoTjdZ ( real64 Z ) { - real64 drhoTj = -711.67 * 3.0 * std::pow(Z,2) + 1657.1 * 2.0 * Z - 1419.5; - return drhoTj; -} - -real64 safeUpdate ( real64 S1, real64 S0) { -real64 S_kink1 = 1.0; real64 S_kink2 = 0.0; -if ( S0 < S_kink1 && S1 > (S_kink1 + eps2)) { - S1 = S_kink1 - eps2; -} -if ( S0 > (S_kink2 + eps2) && S1 < (S_kink2 + eps2)) { - S1 = S_kink2 + eps2; -} -return S1; -} +// real64 computedrhoTidZ ( real64 Z ) { +// real64 drhoTi = -359.23 * 3.0 * std::pow(Z,2) + 921.51 * 2.0 * Z - 911.94; +// return drhoTi; +// } + +// real64 computerhoTj ( real64 Z ) { +// real64 rhoTj = -711.67 * std::pow(Z,3) + 1657.1 * Z * Z - 1419.5 * Z + 1009.6; +// return rhoTj; +// } + +// real64 computedrhoTjdZ ( real64 Z ) { +// real64 drhoTj = -711.67 * 3.0 * std::pow(Z,2) + 1657.1 * 2.0 * Z - 1419.5; +// return drhoTj; +// } + +// real64 safeUpdate ( real64 S1, real64 S0) { +// real64 S_kink1 = 1.0; real64 S_kink2 = 0.0; +// if ( S0 < S_kink1 && S1 > (S_kink1 + eps2)) { +// S1 = S_kink1 - eps2; +// } +// if ( S0 > (S_kink2 + eps2) && S1 < (S_kink2 + eps2)) { +// S1 = S_kink2 + eps2; +// } +// return S1; +// } // TODO make input parameter @@ -874,8 +876,8 @@ struct PPUICPhaseFlux Z_alpha = 1.0 - eps2; } - real64 rho_t_i = computerhoTi ( Zi ); - real64 v_l_i = computeVli ( Zi ); + real64 rho_t_i = computerhoTi >( Zi ); + real64 v_l_i = computeVli >( Zi ); real64 Si = v_l_i * rho_t_i / rho_l_i; if (Si < 0.0) { @@ -884,20 +886,20 @@ struct PPUICPhaseFlux if (Si > 1.0) { Si = 1.0; } - real64 pc_i = computePCalpha (Si); - real64 dpcdSi = computedPCdSalpha (Si); - - real64 rho_t_j = computerhoTj ( Zj ); - real64 v_l_j = computeVlj ( Zj ); - real64 Sj = v_l_j * rho_t_j / rho_l_j; - if (Sj < 0.0) { - Sj = 0.0; - } - if (Sj > 1.0) { - Sj = 1.0; - } - real64 pc_j = computePCbeta (Sj); - real64 dpcdSj = computedPCdSbeta (Sj); + real64 pc_i = computePCalpha >(Si); + real64 dpcdSi = computedPCdSalpha >(Si); + + // real64 rho_t_j = computerhoTj > ( Zj ); + // real64 v_l_j = computeVlj > ( Zj ); + // real64 Sj = v_l_j * rho_t_j / rho_l_j; + // if (Sj < 0.0) { + // Sj = 0.0; + // } + // if (Sj > 1.0) { + // Sj = 1.0; + // } + // real64 pc_j = computePCbeta > (Sj); + // real64 dpcdSj = computedPCdSbeta > (Sj); real64 Ut = totFlux; From ea3c18e9b2dac2cf48cd6afb7faac444115a4f93 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Tue, 6 Aug 2024 17:42:23 -0700 Subject: [PATCH 7/7] Added accessor for globalCompDensity --- ...othermalCompositionalMultiphaseFVMKernelUtilities.hpp | 3 +++ .../IsothermalCompositionalMultiphaseFVMKernels.hpp | 9 +++++++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp index 446fa51675e..59c672eee9d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp @@ -751,6 +751,7 @@ struct PPUICPhaseFlux 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, + ElementViewConst< arrayView2d< real64 const, compflow::USD_COMP > > const & globalCompDensity, localIndex & k_up, real64 & potGrad, real64 ( &phaseFlux ), @@ -914,6 +915,8 @@ struct PPUICPhaseFlux // + GEOS_UNUSED_VAR(mobility, rho_l_j, pi, pj, Ti, Tj, pc_i, dpcdSi, Ut, tol, max_iter, globalCompDensity, compFlux, dCompFlux_dP, dCompFlux_dC); + } }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp index f40f4eb2369..431a25570ff 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp @@ -304,9 +304,10 @@ class FaceBasedAssemblyKernelBase using DofNumberAccessor = ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > >; using CompFlowAccessors = - StencilAccessors< fields::ghostRank, + StencilAccessors< fields::ghostRank, fields::flow::gravityCoefficient, fields::flow::pressure, + fields::flow::globalCompDensity, fields::flow::dGlobalCompFraction_dGlobalCompDensity, fields::flow::phaseVolumeFraction, fields::flow::dPhaseVolumeFraction, @@ -474,7 +475,8 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase m_stencilWrapper( stencilWrapper ), m_seri( stencilWrapper.getElementRegionIndices() ), m_sesri( stencilWrapper.getElementSubRegionIndices() ), - m_sei( stencilWrapper.getElementIndices() ) + m_sei( stencilWrapper.getElementIndices() ), + m_globalCompDensity( compFlowAccessors.get( fields::flow::globalCompDensity {} ) ) // TODO: Maybe only store this for IC solver where its needed? { } /** @@ -692,6 +694,7 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase m_dCompFrac_dCompDens, m_phaseMassDens, m_dPhaseMassDens, m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac, + m_globalCompDensity, k_up, potGrad, phaseFlux, @@ -870,6 +873,8 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase typename STENCILWRAPPER::IndexContainerViewConstType const m_sesri; typename STENCILWRAPPER::IndexContainerViewConstType const m_sei; + ElementViewConst< arrayView2d < real64 const, compflow::USD_COMP > > const m_globalCompDensity; + }; /**