diff --git a/src/coreComponents/finiteVolume/FluxApproximationBase.hpp b/src/coreComponents/finiteVolume/FluxApproximationBase.hpp index fafce71f6f..449b0ee3db 100644 --- a/src/coreComponents/finiteVolume/FluxApproximationBase.hpp +++ b/src/coreComponents/finiteVolume/FluxApproximationBase.hpp @@ -39,7 +39,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 + HU2PH ///< HU simplified 2-phase version }; /** @@ -48,7 +49,8 @@ enum class UpwindingScheme : integer ENUM_STRINGS( UpwindingScheme, "PPU", "C1PPU", - "IHU" ); + "IHU", + "HU2PH" ); /** * @struct UpwindingParameters @@ -57,7 +59,7 @@ ENUM_STRINGS( UpwindingScheme, */ struct UpwindingParameters { - /// PPU or C1-PPU or IHU + /// PPU or C1-PPU or IHU or HU2PH UpwindingScheme upwindingScheme; /// C1-PPU smoothing tolerance diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt index 651d9135e0..a98227700f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt @@ -62,6 +62,7 @@ set( physicsSolvers_headers fluidFlow/kernels/compositional/GlobalComponentFractionKernel.hpp fluidFlow/kernels/compositional/HydrostaticPressureKernel.hpp fluidFlow/kernels/compositional/IHUPhaseFlux.hpp + fluidFlow/kernels/compositional/HU2PhaseFlux.hpp fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp fluidFlow/kernels/compositional/PhaseComponentFlux.hpp fluidFlow/kernels/compositional/PhaseMobilityKernel.hpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index aeb88ef68e..8261b44a64 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -172,6 +172,13 @@ void CompositionalMultiphaseFVM::initializePreSubGroups() { GEOS_ERROR( "A discretization deriving from FluxApproximationBase must be selected with CompositionalMultiphaseFlow" ); } + else + { + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName ); + GEOS_ERROR_IF( fluxApprox.upwindingParams().upwindingScheme == UpwindingScheme::HU2PH && m_numPhases != 2, + GEOS_FMT( "{}: upwinding scheme {} only supports 2-phase flow", + getName(), EnumStrings< UpwindingScheme >::toString( UpwindingScheme::HU2PH ))); + } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp index 3d79301623..c781f02f1e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp @@ -26,7 +26,6 @@ #include "constitutive/capillaryPressure/layouts.hpp" #include "mesh/ElementRegionManager.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp" namespace geos @@ -67,7 +66,6 @@ struct C1PPUPhaseFlux * @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 @@ -91,21 +89,15 @@ struct C1PPUPhaseFlux 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< 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] ) + real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] ) { real64 dPotGrad_dTrans = 0.0; real64 dPresGrad_dP[numFluxSupportPoints]{}; @@ -214,13 +206,6 @@ struct C1PPUPhaseFlux } potGrad = potGrad * Ttrans; - - // choose upstream cell for composition upwinding - k_up = (phaseFlux >= 0) ? 0 : 1; - - //distribute on phaseComponentFlux here - PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei, phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, phaseFlux - , dPhaseFlux_dP, dPhaseFlux_dC, compFlux, dCompFlux_dP, dCompFlux_dC ); } }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp index 4588295330..4d36066b97 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp @@ -35,6 +35,8 @@ #include "physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/HU2PhaseFlux.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp" namespace geos { @@ -228,6 +230,7 @@ class FluxComputeKernel : public FluxComputeKernelBase StackVariables & stack, FUNC && compFluxKernelOp = NoOpFunc{} ) const { + using namespace isothermalCompositionalMultiphaseFVMKernelUtilities; // first, compute the transmissibilities at this face m_stencilWrapper.computeWeights( iconn, @@ -269,12 +272,11 @@ class FluxComputeKernel : public FluxComputeKernelBase real64 phaseFlux = 0.0; real64 dPhaseFlux_dP[numFluxSupportPoints]{}; real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{}; - - localIndex k_up = -1; + real64 dPhaseFlux_dTrans = 0.0; // not really used if( m_kernelFlags.isSet( KernelFlags::C1PPU ) ) { - isothermalCompositionalMultiphaseFVMKernelUtilities::C1PPUPhaseFlux::compute< numComp, numFluxSupportPoints > + C1PPUPhaseFlux::compute< numComp, numFluxSupportPoints > ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), @@ -286,22 +288,17 @@ class FluxComputeKernel : public FluxComputeKernelBase m_gravCoef, m_phaseMob, m_dPhaseMob, m_phaseVolFrac, 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 ); + dPhaseFlux_dC ); } else if( m_kernelFlags.isSet( KernelFlags::IHU ) ) { - isothermalCompositionalMultiphaseFVMKernelUtilities::IHUPhaseFlux::compute< numComp, numFluxSupportPoints > + IHUPhaseFlux::compute< numComp, numFluxSupportPoints > ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), @@ -313,22 +310,39 @@ class FluxComputeKernel : public FluxComputeKernelBase m_gravCoef, m_phaseMob, m_dPhaseMob, m_phaseVolFrac, 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 ); + dPhaseFlux_dC ); + } + else if( m_kernelFlags.isSet( KernelFlags::HU2PH ) ) + { + HU2PhaseFlux::compute< numComp, numFluxSupportPoints > + ( m_numPhases, + ip, + m_kernelFlags.isSet( KernelFlags::CapPressure ), + m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ), + seri, sesri, sei, + trans, + dTrans_dPres, + m_pres, + m_gravCoef, + m_phaseMob, m_dPhaseMob, + m_phaseVolFrac, m_dPhaseVolFrac, + m_dCompFrac_dCompDens, + m_phaseMassDens, m_dPhaseMassDens, + m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac, + potGrad, + phaseFlux, + dPhaseFlux_dP, + dPhaseFlux_dC ); } else { - isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux::compute< numComp, numFluxSupportPoints > + PPUPhaseFlux::compute< numComp, numFluxSupportPoints > ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), @@ -340,21 +354,25 @@ class FluxComputeKernel : public FluxComputeKernelBase m_gravCoef, m_phaseMob, m_dPhaseMob, m_phaseVolFrac, 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, - dCompFlux_dTrans ); + dPhaseFlux_dTrans ); } + // choose upstream cell for composition upwinding + localIndex const k_up = (phaseFlux >= 0) ? 0 : 1; + + // distribute on phaseComponentFlux here + PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei, + m_phaseCompFrac, m_dPhaseCompFrac, m_dCompFrac_dCompDens, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans, + compFlux, dCompFlux_dP, dCompFlux_dC, dCompFlux_dTrans ); + // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives // possible use: assemble the derivatives wrt temperature, and the flux term of the energy equation for this phase compFluxKernelOp( ip, m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp index c7b85da8d9..bc46059cec 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp @@ -56,9 +56,9 @@ enum class KernelFlags /// Flag indicating whether C1-PPU is used or not C1PPU = 1 << 5, // 32 /// Flag indicating whether IHU is used or not - IHU = 1 << 6 // 64 - /// Add more flags like that if needed: - // Flag8 = 1 << 7 //128 + IHU = 1 << 6, // 64 + /// Flag indicating whether HU 2-phase simplified version is used or not + HU2PH = 1 << 7 // 128 }; /******************************** FluxComputeKernelBase ********************************/ @@ -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 new file mode 100644 index 0000000000..391629ec5a --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/HU2PhaseFlux.hpp @@ -0,0 +1,666 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file IHU2PhaseFlux.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHU2PHASEFLUX_HPP +#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHU2PHASEFLUX_HPP + +#include "common/DataLayouts.hpp" +#include "common/DataTypes.hpp" +#include "constitutive/fluid/multifluid/Layouts.hpp" +#include "constitutive/capillaryPressure/layouts.hpp" +#include "mesh/ElementRegionManager.hpp" + + +namespace geos +{ + +namespace isothermalCompositionalMultiphaseFVMKernelUtilities +{ + +template< typename VIEWTYPE > +using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + +using Deriv = constitutive::multifluid::DerivativeOffset; + +/*** HU 2 phase simplified version ***/ + +struct HU2PhaseFlux +{ + + static constexpr double minTotMob = 1e-12; + + /** + * @brief Simplified 2-phase version of hybrid upwinding + * @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 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, + 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 & GEOS_UNUSED_PARAM( potGrad ), + real64 ( &phaseFlux ), + real64 ( & dPhaseFlux_dP )[numFluxSupportPoints], + real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] ) + { + // viscous part + computeViscousFlux< numComp, numFluxSupportPoints >( ip, numPhase, + hasCapPressure, + checkPhasePresenceInGravity, + seri, sesri, sei, + trans, dTrans_dPres, + pres, gravCoef, + dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, + phaseMob, dPhaseMob, + phaseVolFrac, dPhaseVolFrac, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + 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, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + + // capillary part + if( hasCapPressure ) + { + computeCapillaryFlux< numComp, numFluxSupportPoints >( ip, numPhase, + seri, sesri, sei, + trans, dTrans_dPres, + phaseMob, dPhaseMob, + dPhaseVolFrac, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + } + } + +protected: + + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + computeViscousFlux( integer const & ip, + integer const & numPhase, + 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< arrayView3d< real64 const > > const & dCompFrac_dCompDens, + ElementViewConst< arrayView3d< real64 const > > const & phaseMassDens, + ElementViewConst< arrayView4d< real64 const > > const & dPhaseMassDens, + ElementViewConst< arrayView2d< real64 const, 1 > > const & phaseMob, + ElementViewConst< arrayView3d< real64 const > > const & dPhaseMob, + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, + ElementViewConst< arrayView3d< real64 const > > const & dPhaseVolFrac, + ElementViewConst< arrayView3d< real64 const > > const & phaseCapPressure, + ElementViewConst< arrayView4d< real64 const > > const & dPhaseCapPressure_dPhaseVolFrac, + real64 ( &phaseFlux ), + real64 ( & dPhaseFlux_dP )[numFluxSupportPoints], + real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] ) + { + // form total velocity and derivatives (TODO move it OUT!) + real64 totFlux{}; + real64 dTotFlux_dP[numFluxSupportPoints]{}; + real64 dTotFlux_dC[numFluxSupportPoints][numComp]{}; + + computeTotalFlux( numPhase, + hasCapPressure, + checkPhasePresenceInGravity, + seri, sesri, sei, + trans, dTrans_dPres, + pres, gravCoef, + phaseMob, dPhaseMob, + phaseVolFrac, dPhaseVolFrac, + dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + totFlux, dTotFlux_dP, dTotFlux_dC ); + + localIndex k_up = -1; + real64 mob{}; + real64 dMob_dP{}; + real64 dMob_dC[numComp]{}; + + real64 totMob{}; + real64 dTotMob_dP{}; + real64 dTotMob_dC[numComp]{}; + + for( localIndex jp = 0; jp < numPhase; ++jp ) + { + if( jp == ip ) // ip will be computed later + continue; + + // upwind based on totFlux sign + upwindMobility< numComp, numFluxSupportPoints >( jp, + seri, + sesri, + sei, + totFlux, + phaseMob, + dPhaseMob, + k_up, + mob, + dMob_dP, + dMob_dC ); + // accumulate total mobility + UpwindHelpers::addToValueAndDerivatives( mob, dMob_dP, dMob_dC, totMob, dTotMob_dP, dTotMob_dC ); + } + + // upwind based on totFlux sign + upwindMobility< numComp, numFluxSupportPoints >( ip, + seri, + sesri, + sei, + totFlux, + phaseMob, + dPhaseMob, + k_up, + mob, + dMob_dP, + dMob_dC ); + // accumulate total mobility + UpwindHelpers::addToValueAndDerivatives( mob, dMob_dP, dMob_dC, totMob, dTotMob_dP, dTotMob_dC ); + + // safeguard + totMob = LvArray::math::max( totMob, minTotMob ); + real64 const invTotMob = 1 / totMob; + + // fractional flow for viscous part as \lambda_i^{up}/\sum_{NP}(\lambda_j^{up}) + + real64 const fractionalFlow = mob * invTotMob; + real64 dFractionalFlow_dP{}; + real64 dFractionalFlow_dC[numComp]{}; + UpwindHelpers::addToDerivativesScaled( dMob_dP, dMob_dC, invTotMob, dFractionalFlow_dP, dFractionalFlow_dC ); + UpwindHelpers::addToDerivativesScaled( dTotMob_dP, dTotMob_dC, -fractionalFlow * invTotMob, dFractionalFlow_dP, dFractionalFlow_dC ); + + /// Assembling the viscous flux (and derivatives) from fractional flow and total velocity as \phi_{\mu} = f_i^{up,\mu} uT + + real64 const viscousPhaseFlux = fractionalFlow * totFlux; + real64 dViscousPhaseFlux_dP[numFluxSupportPoints]{}; + real64 dViscousPhaseFlux_dC[numFluxSupportPoints][numComp]{}; + + // fractionalFlow derivatives + UpwindHelpers::addToDerivativesScaled( dFractionalFlow_dP, dFractionalFlow_dC, totFlux, dViscousPhaseFlux_dP[k_up], dViscousPhaseFlux_dC[k_up] ); + + // Ut derivatives + UpwindHelpers::addToDerivativesScaled( dTotFlux_dP, dTotFlux_dC, fractionalFlow, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC ); + + // accumulate in the flux and its derivatives (need to be very careful doing that) + UpwindHelpers::addToValueAndDerivatives( viscousPhaseFlux, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + } + + template< localIndex numComp, localIndex numFluxSupportPoints > + 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, + ElementViewConst< arrayView3d< real64 const > > const & phaseMassDens, + ElementViewConst< arrayView4d< real64 const > > const & dPhaseMassDens, + real64 & phaseFlux, + real64 (& dPhaseFlux_dP)[numFluxSupportPoints], + real64 (& dPhaseFlux_dC)[numFluxSupportPoints][numComp] ) + { + /// Assembling the gravitational flux (and derivatives) + real64 gravPhaseFlux{}; + real64 dGravPhaseFlux_dP[numFluxSupportPoints]{}; + real64 dGravPhaseFlux_dC[numFluxSupportPoints][numComp]{}; + + real64 pot_i{}; + real64 dPot_i_dP[numFluxSupportPoints]{}; + real64 dPot_i_dC[numFluxSupportPoints][numComp]{}; + computeGravityPotential< numComp, numFluxSupportPoints >( ip, + seri, + sesri, + sei, + checkPhasePresenceInGravity, + trans, + dTrans_dPres, + gravCoef, + dCompFrac_dCompDens, + phaseVolFrac, + phaseMassDens, + dPhaseMassDens, + pot_i, + dPot_i_dP, + dPot_i_dC ); + + for( localIndex jp = 0; jp < numPhase; ++jp ) + { + if( ip != jp ) + { + real64 pot_j{}; + real64 dPot_j_dP[numFluxSupportPoints]{}; + real64 dPot_j_dC[numFluxSupportPoints][numComp]{}; + computeGravityPotential< numComp, numFluxSupportPoints >( jp, + seri, + sesri, + sei, + checkPhasePresenceInGravity, + trans, + dTrans_dPres, + gravCoef, + dCompFrac_dCompDens, + phaseVolFrac, + phaseMassDens, + dPhaseMassDens, + pot_j, + dPot_j_dP, + dPot_j_dC ); + + computePotDiffFlux( ip, jp, + seri, sesri, sei, + pot_i, dPot_i_dP, dPot_i_dC, + pot_j, dPot_j_dP, dPot_j_dC, + phaseMob, dPhaseMob, + gravPhaseFlux, dGravPhaseFlux_dP, dGravPhaseFlux_dC ); + } + } + + // update phaseFlux from gravitational + UpwindHelpers::addToValueAndDerivatives( gravPhaseFlux, dGravPhaseFlux_dP, dGravPhaseFlux_dC, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + } + + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + computeCapillaryFlux( integer const & ip, integer const & numPhase, + localIndex const (&seri)[numFluxSupportPoints], + localIndex const (&sesri)[numFluxSupportPoints], + localIndex const (&sei)[numFluxSupportPoints], + real64 const (&trans)[2], real64 const (&dTrans_dPres)[2], + ElementViewConst< arrayView2d< real64 const, 1 > > const & phaseMob, + ElementViewConst< arrayView3d< real64 const > > const & dPhaseMob, + ElementViewConst< arrayView3d< real64 const > > const & dPhaseVolFrac, + ElementViewConst< arrayView3d< real64 const > > const & phaseCapPressure, + ElementViewConst< arrayView4d< real64 const > > const & dPhaseCapPressure_dPhaseVolFrac, + real64 & phaseFlux, + real64 (& dPhaseFlux_dP)[numFluxSupportPoints], + real64 (& dPhaseFlux_dC)[numFluxSupportPoints][numComp] ) + { + /// Assembling the capillary flux (and derivatives) + real64 capPhaseFlux{}; + real64 dCapPhaseFlux_dP[numFluxSupportPoints]{}; + real64 dCapPhaseFlux_dC[numFluxSupportPoints][numComp]{}; + + real64 pot_i{}; + real64 dPot_i_dP[numFluxSupportPoints]{}; + real64 dPot_i_dC[numFluxSupportPoints][numComp]{}; + + computeCapillaryPotential< numComp, numFluxSupportPoints >( ip, + numPhase, + seri, + sesri, + sei, + trans, + dTrans_dPres, + dPhaseVolFrac, + phaseCapPressure, + dPhaseCapPressure_dPhaseVolFrac, + pot_i, + dPot_i_dP, + dPot_i_dC ); + + for( localIndex jp = 0; jp < numPhase; ++jp ) + { + if( ip != jp ) + { + real64 pot_j{}; + real64 dPot_j_dP[numFluxSupportPoints]{}; + real64 dPot_j_dC[numFluxSupportPoints][numComp]{}; + computeCapillaryPotential< numComp, numFluxSupportPoints >( jp, + numPhase, + seri, + sesri, + sei, + trans, + dTrans_dPres, + dPhaseVolFrac, + phaseCapPressure, + dPhaseCapPressure_dPhaseVolFrac, + pot_j, + dPot_j_dP, + dPot_j_dC ); + + computePotDiffFlux( ip, jp, + seri, sesri, sei, + pot_i, dPot_i_dP, dPot_i_dC, + pot_j, dPot_j_dP, dPot_j_dC, + phaseMob, dPhaseMob, + capPhaseFlux, dCapPhaseFlux_dP, dCapPhaseFlux_dC ); + } + } + + // update phaseFlux from capillary flux + UpwindHelpers::addToValueAndDerivatives( capPhaseFlux, dCapPhaseFlux_dP, dCapPhaseFlux_dC, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + } + + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + computeTotalFlux( integer const & numPhase, + const integer & hasCapPressure, + const integer & 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, 1 > > const & phaseMob, + ElementViewConst< arrayView3d< real64 const > > const & dPhaseMob, + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, + ElementViewConst< arrayView3d< real64 const > > const & dPhaseVolFrac, + ElementViewConst< arrayView3d< real64 const > > const & dCompFrac_dCompDens, + ElementViewConst< arrayView3d< real64 const > > const & phaseMassDens, + ElementViewConst< arrayView4d< real64 const > > const & dPhaseMassDens, + ElementViewConst< arrayView3d< real64 const > > const & phaseCapPressure, + ElementViewConst< arrayView4d< real64 const > > const & dPhaseCapPressure_dPhaseVolFrac, + real64 & totFlux, real64 (& dTotFlux_dP)[numFluxSupportPoints], real64 (& dTotFlux_dC)[numFluxSupportPoints][numComp] ) + { + for( integer jp = 0; jp < numPhase; ++jp ) + { + // working arrays for phase flux + real64 potGrad{}; + real64 phaseFlux{}; + real64 dPhaseFlux_dP[numFluxSupportPoints]{}; + real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{}; + real64 dPhaseFlux_dTrans{}; + PPUPhaseFlux::compute( numPhase, jp, + hasCapPressure, + checkPhasePresenceInGravity, + seri, sesri, sei, + trans, dTrans_dPres, + pres, gravCoef, + phaseMob, dPhaseMob, + phaseVolFrac, dPhaseVolFrac, + dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans ); + + UpwindHelpers::addToValueAndDerivatives( phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, + totFlux, dTotFlux_dP, dTotFlux_dC ); + } + } + + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + upwindMobility( localIndex const ip, + localIndex const (&seri)[numFluxSupportPoints], + localIndex const (&sesri)[numFluxSupportPoints], + localIndex const (&sei)[numFluxSupportPoints], + real64 const pot, + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, + ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, + localIndex & upwindDir, + real64 & mobility, + real64 & dMobility_dP, + real64 ( & dMobility_dC)[numComp] ) + { + upwindDir = (pot > 0) ? 0 : 1; + UpwindHelpers::assignToZero( mobility, dMobility_dP, dMobility_dC ); + UpwindHelpers::assignMobilityAndDerivatives( ip, upwindDir, seri, sesri, sei, phaseMob, dPhaseMob, mobility, dMobility_dP, dMobility_dC ); + } + + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void computeGravityPotential( localIndex const ip, + 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, + real64 ( & dGravPot_dP )[numFluxSupportPoints], + real64 ( & dGravPot_dC )[numFluxSupportPoints][numComp] ) + { + // init + UpwindHelpers::assignToZero( gravPot, dGravPot_dP, dGravPot_dC ); + + // get average density + real64 densMean{}; + 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 ); + + // compute potential difference MPFA-style + for( localIndex i = 0; i < numFluxSupportPoints; ++i ) + { + localIndex const er = seri[i]; + localIndex const esr = sesri[i]; + localIndex const ei = sei[i]; + + real64 const gravD = trans[i] * gravCoef[er][esr][ei]; + real64 const dGravD_dP = dTrans_dP[i] * gravCoef[er][esr][ei]; + gravPot += densMean * gravD; + dGravPot_dP[i] += densMean * dGravD_dP; + + // need to add contributions from both cells the mean density depends on + UpwindHelpers::addToDerivativesScaled( dDensMean_dP, dDensMean_dC, gravD, dGravPot_dP, dGravPot_dC ); + } + + } + + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void computeCapillaryPotential( localIndex const ip, + localIndex const numPhase, + localIndex const (&seri)[numFluxSupportPoints], + localIndex const (&sesri)[numFluxSupportPoints], + localIndex const (&sei)[numFluxSupportPoints], + real64 const (&trans)[2], + real64 const (&dTrans_dPres)[2], + ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, + ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, + ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, + real64 & capPot, + real64 ( & dCapPot_dP )[numFluxSupportPoints], + real64 ( & dCapPot_dC )[numFluxSupportPoints][numComp] ) + { + // init + UpwindHelpers::assignToZero( capPot, dCapPot_dP, dCapPot_dC ); + + // need to add contributions from both cells + for( localIndex i = 0; i < numFluxSupportPoints; ++i ) + { + localIndex const er = seri[i]; + localIndex const esr = sesri[i]; + localIndex const ei = sei[i]; + + real64 const capPressure = phaseCapPressure[er][esr][ei][0][ip]; + real64 dCapPressure_dP{}; + real64 dCapPressure_dC[numComp]{}; + for( integer jp = 0; jp < numPhase; ++jp ) + { + real64 const dCapPressure_dS = dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp]; + dCapPressure_dP += dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dP]; + + for( integer jc = 0; jc < numComp; ++jc ) + { + dCapPressure_dC[jc] += dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dC+jc]; + } + } + + capPot += trans[i] * capPressure; + dCapPot_dP[i] += trans[i] * dCapPressure_dP + dTrans_dPres[i] * capPressure; + for( integer jc = 0; jc < numComp; ++jc ) + { + dCapPot_dC[i][jc] += trans[i] * dCapPressure_dC[jc]; + } + } + } + + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + computePotDiffFlux( integer const & ip, integer const & jp, + localIndex const (&seri)[numFluxSupportPoints], + localIndex const (&sesri)[numFluxSupportPoints], + localIndex const (&sei)[numFluxSupportPoints], + real64 const & pot_i, real64 const ( &dPot_i_dP )[numFluxSupportPoints], real64 const (&dPot_i_dC )[numFluxSupportPoints][numComp], + real64 const & pot_j, real64 const ( &dPot_j_dP )[numFluxSupportPoints], real64 const ( &dPot_j_dC)[numFluxSupportPoints][numComp], + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, + ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, + real64 & phaseFlux, real64 ( & dPhaseFlux_dP )[numFluxSupportPoints], real64 ( & dPhaseFlux_dC)[numFluxSupportPoints][numComp] ) + { + // upwind based on pot diff sign + real64 const potDiff = pot_j - pot_i; + real64 dPotDiff_dP[numFluxSupportPoints]{}; + real64 dPotDiff_dC[numFluxSupportPoints][numComp]{}; + UpwindHelpers::addToDerivativesScaled( dPot_j_dP, dPot_j_dC, 1.0, dPotDiff_dP, dPotDiff_dC ); + UpwindHelpers::addToDerivativesScaled( dPot_i_dP, dPot_i_dC, -1.0, dPotDiff_dP, dPotDiff_dC ); + + localIndex k_up_i = -1; + real64 mob_i{}; + real64 dMob_i_dP{}; + real64 dMob_i_dC[numComp]{}; + upwindMobility< numComp, numFluxSupportPoints >( ip, + seri, + sesri, + sei, + potDiff, + phaseMob, + dPhaseMob, + k_up_i, + mob_i, + dMob_i_dP, + dMob_i_dC ); + localIndex k_up_j = -1; + real64 mob_j{}; + real64 dMob_j_dP{}; + real64 dMob_j_dC[numComp]{}; + upwindMobility< numComp, numFluxSupportPoints >( jp, + seri, + sesri, + sei, + -potDiff, + phaseMob, + dPhaseMob, + k_up_j, + mob_j, + dMob_j_dP, + dMob_j_dC ); + + // safeguard + real64 const mobTot = LvArray::math::max( mob_i + mob_j, minTotMob ); + real64 const mobTotInv = 1 / mobTot; + real64 dMobTot_dP[numFluxSupportPoints]{}; + real64 dMobTot_dC[numFluxSupportPoints][numComp]{}; + UpwindHelpers::addToDerivatives( dMob_i_dP, dMob_i_dC, dMobTot_dP[k_up_i], dMobTot_dC[k_up_i] ); + UpwindHelpers::addToDerivatives( dMob_j_dP, dMob_j_dC, dMobTot_dP[k_up_j], dMobTot_dC[k_up_j] ); + + // Assembling flux phase-wise as \phi_{i,g} = \sum_{k\nei} \lambda_k^{up,g} f_k^{up,g} (Pot_i - Pot_k) + phaseFlux += mob_i * mob_j * mobTotInv * potDiff; + + // mob_i derivatives + UpwindHelpers::addToDerivativesScaled( dMob_i_dP, dMob_i_dC, mob_j * mobTotInv * potDiff, dPhaseFlux_dP[k_up_i], dPhaseFlux_dC[k_up_i] ); + + // mob_j derivatives + UpwindHelpers::addToDerivativesScaled( dMob_j_dP, dMob_j_dC, mob_i * mobTotInv * potDiff, dPhaseFlux_dP[k_up_j], dPhaseFlux_dC[k_up_j] ); + + // mobTot derivatives + real64 const mobTotInv2 = mobTotInv * mobTotInv; + UpwindHelpers::addToDerivativesScaled( dMobTot_dP, dMobTot_dC, -mob_i * mob_j * mobTotInv2 * potDiff, dPhaseFlux_dP, dPhaseFlux_dC ); + + // potDiff derivatives + UpwindHelpers::addToDerivativesScaled( dPotDiff_dP, dPotDiff_dC, mob_i * mob_j * mobTotInv, dPhaseFlux_dP, dPhaseFlux_dC ); + } + +}; + +} // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities + +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHU2PHASEFLUX_HPP diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp index 97f5bdd5a2..607fea61b3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp @@ -42,7 +42,143 @@ using Deriv = constitutive::multifluid::DerivativeOffset; namespace UpwindHelpers { -static constexpr double minTotMob = 1e-12; +template< localIndex numComp > +GEOS_HOST_DEVICE +static void assignToZero( real64 & deriv_dP, real64 ( & deriv_dC )[numComp] ) +{ + deriv_dP = 0.0; + for( localIndex ic = 0; ic < numComp; ++ic ) + { + deriv_dC[ic] = 0.0; + } +} + +template< localIndex numComp > +GEOS_HOST_DEVICE +static void assignToZero( real64 & value, real64 & deriv_dP, real64 ( & deriv_dC )[numComp] ) +{ + value = 0.0; + assignToZero( deriv_dP, deriv_dC ); +} + +template< localIndex numComp, localIndex numFluxSupportPoints > +GEOS_HOST_DEVICE +static void assignToZero( real64 & value, real64 ( & deriv_dP )[numFluxSupportPoints], real64 ( & deriv_dC )[numFluxSupportPoints][numComp] ) +{ + value = 0; + for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) + { + assignToZero( deriv_dP[ke], deriv_dC[ke] ); + } +} + +template< localIndex numComp > +GEOS_HOST_DEVICE +static void addToDerivativesScaled( real64 const ( &dDeriv_dP ), real64 const ( &dDeriv_dC )[numComp], + real64 const & factor, + real64 ( &deriv_dP ), real64 ( & deriv_dC )[numComp] ) +{ + deriv_dP += dDeriv_dP * factor; + for( localIndex ic = 0; ic < numComp; ++ic ) + deriv_dC[ic] += dDeriv_dC[ic] * factor; +} + +template< localIndex numComp, localIndex numFluxSupportPoints > +GEOS_HOST_DEVICE +static void addToDerivativesScaled( real64 const ( &dDeriv_dP )[numFluxSupportPoints], real64 const ( &dDeriv_dC )[numFluxSupportPoints][numComp], + real64 const & factor, + real64 ( & deriv_dP )[numFluxSupportPoints], real64 ( & deriv_dC )[numFluxSupportPoints][numComp] ) +{ + for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) + { + addToDerivativesScaled( dDeriv_dP[ke], dDeriv_dC[ke], factor, deriv_dP[ke], deriv_dC[ke] ); + } +} + +template< localIndex numComp > +GEOS_HOST_DEVICE +static void addToDerivatives( real64 const ( &dDeriv_dP ), real64 const ( &dDeriv_dC )[numComp], + real64 ( &deriv_dP ), real64 ( & deriv_dC )[numComp] ) +{ + addToDerivativesScaled( dDeriv_dP, dDeriv_dC, 1.0, deriv_dP, deriv_dC ); +} + +template< localIndex numComp, localIndex numFluxSupportPoints > +GEOS_HOST_DEVICE +static void addToDerivatives( real64 const ( &dDeriv_dP )[numFluxSupportPoints], real64 const ( &dDeriv_dC )[numFluxSupportPoints][numComp], + real64 ( & deriv_dP )[numFluxSupportPoints], real64 ( & deriv_dC )[numFluxSupportPoints][numComp] ) +{ + for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) + { + addToDerivativesScaled( dDeriv_dP, dDeriv_dC, 1.0, deriv_dP, deriv_dC ); + } +} + +template< localIndex numComp > +GEOS_HOST_DEVICE +static void addToValueAndDerivatives( real64 const & dValue, real64 const ( &dDeriv_dP ), real64 const ( &dDeriv_dC )[numComp], + real64 & value, real64 ( &deriv_dP ), real64 ( & deriv_dC )[numComp] ) +{ + value += dValue; + addToDerivatives( dDeriv_dP, dDeriv_dC, deriv_dP, deriv_dC ); +} + +template< localIndex numComp, localIndex numFluxSupportPoints > +GEOS_HOST_DEVICE +static void addToValueAndDerivatives( real64 const & dValue, real64 const ( &dDeriv_dP )[numFluxSupportPoints], real64 const ( &dDeriv_dC )[numFluxSupportPoints][numComp], + real64 & value, real64 ( & deriv_dP )[numFluxSupportPoints], real64 ( & deriv_dC )[numFluxSupportPoints][numComp] ) +{ + value += dValue; + for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) + { + addToDerivatives( dDeriv_dP[ke], dDeriv_dC[ke], deriv_dP[ke], deriv_dC[ke] ); + } +} + +template< localIndex numComp, localIndex numFluxSupportPoints > +GEOS_HOST_DEVICE +static void assignMobilityAndDerivatives( localIndex const & ip, localIndex const & upwindDir, + localIndex const (&seri)[numFluxSupportPoints], + localIndex const (&sesri)[numFluxSupportPoints], + localIndex const (&sei)[numFluxSupportPoints], + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, + ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, + real64 & mobility, real64 & dMobility_dP, real64 ( & dMobility_dC )[numComp] ) +{ + localIndex const er_up = seri[upwindDir]; + localIndex const esr_up = sesri[upwindDir]; + localIndex const ei_up = sei[upwindDir]; + + if( std::fabs( phaseMob[er_up][esr_up][ei_up][ip] ) > 1e-20 ) + { + mobility += phaseMob[er_up][esr_up][ei_up][ip]; + dMobility_dP += dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP]; + for( localIndex ic = 0; ic < numComp; ++ic ) + { + dMobility_dC[ic] += dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dC + ic]; + } + } +} + +template< localIndex numComp, localIndex numFluxSupportPoints > +GEOS_HOST_DEVICE +static void computeFractionalFlowAndDerivatives( localIndex const & k_up, real64 const & mob, real64 const & dMob_dP, real64 const ( &dMob_dC )[numComp], + real64 const & totMob, real64 const ( &dTotMob_dP )[numFluxSupportPoints], real64 const ( &dTotMob_dC )[numFluxSupportPoints][numComp], + real64 & fractionalFlow, real64 ( & dFractionalFlow_dP )[numFluxSupportPoints], real64 ( & dFractionalFlow_dC )[numFluxSupportPoints][numComp] ) +{ + // guard against no flow region + // fractional flow too low to let the upstream phase flow + if( std::fabs( mob ) > 1e-20 ) + { + real64 const invTotMob = 1 / totMob; + + fractionalFlow = mob * invTotMob; + + addToDerivativesScaled( dMob_dP, dMob_dC, invTotMob, dFractionalFlow_dP[k_up], dFractionalFlow_dC[k_up] ); + + addToDerivativesScaled( dTotMob_dP, dTotMob_dC, -fractionalFlow * invTotMob, dFractionalFlow_dP, dFractionalFlow_dC ); + } +} template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND > GEOS_HOST_DEVICE @@ -54,7 +190,7 @@ upwindMobilityViscous( localIndex const numPhase, localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], - real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place + real64 const totFlux, // in fine should be a ElemnetViewConst once seq form are in place ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, @@ -68,19 +204,9 @@ upwindMobilityViscous( localIndex const numPhase, integer const hasCapPressure, localIndex & upwindDir, real64 & mobility, - real64( &dMobility_dP), - real64 ( & dMobility_dC)[numComp] - ) + real64 (&dMobility_dP), + real64 (& dMobility_dC)[numComp] ) { - - //reinit - mobility = 0.0; - dMobility_dP = 0.0; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMobility_dC[ic] = 0.0; - } - UPWIND scheme; scheme.template getUpwindDirectionViscous< numComp, numFluxSupportPoints, UPWIND >( numPhase, ip, @@ -102,19 +228,9 @@ upwindMobilityViscous( localIndex const numPhase, hasCapPressure, upwindDir ); - localIndex const er_up = seri[upwindDir]; - localIndex const esr_up = sesri[upwindDir]; - localIndex const ei_up = sei[upwindDir]; - - if( std::fabs( phaseMob[er_up][esr_up][ei_up][ip] ) > 1e-20 ) - { - mobility = phaseMob[er_up][esr_up][ei_up][ip]; - dMobility_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP]; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMobility_dC[ic] = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dC + ic]; - } - } + //reinit + assignToZero( mobility, dMobility_dP, dMobility_dC ); + assignMobilityAndDerivatives( ip, upwindDir, seri, sesri, sei, phaseMob, dPhaseMob, mobility, dMobility_dP, dMobility_dC ); } template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND > @@ -143,19 +259,9 @@ upwindMobilityGravity( localIndex const numPhase, integer const checkPhasePresenceInGravity, localIndex & upwindDir, real64 & mobility, - real64( &dMobility_dP), - real64 ( & dMobility_dC)[numComp] - ) + real64 ( &dMobility_dP ), + real64 ( & dMobility_dC )[numComp] ) { - - //reinit - mobility = 0.0; - dMobility_dP = 0.0; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMobility_dC[ic] = 0.0; - } - UPWIND scheme; scheme.template getUpwindDirectionGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase, ip, @@ -179,19 +285,9 @@ upwindMobilityGravity( localIndex const numPhase, checkPhasePresenceInGravity, upwindDir ); - localIndex const er_up = seri[upwindDir]; - localIndex const esr_up = sesri[upwindDir]; - localIndex const ei_up = sei[upwindDir]; - - if( std::fabs( phaseMob[er_up][esr_up][ei_up][ip] ) > 1e-20 ) - { - mobility = phaseMob[er_up][esr_up][ei_up][ip]; - dMobility_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP]; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMobility_dC[ic] = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dC + ic]; - } - } + //reinit + assignToZero( mobility, dMobility_dP, dMobility_dC ); + assignMobilityAndDerivatives( ip, upwindDir, seri, sesri, sei, phaseMob, dPhaseMob, mobility, dMobility_dP, dMobility_dC ); } template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND > @@ -218,19 +314,9 @@ upwindMobilityCapillary( localIndex const numPhase, integer const hasCapPressure, localIndex & upwindDir, real64 & mobility, - real64( &dMobility_dP), - real64 ( & dMobility_dC)[numComp] - ) + real64 ( &dMobility_dP ), + real64 ( & dMobility_dC )[numComp] ) { - - //reinit - mobility = 0.0; - dMobility_dP = 0.0; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMobility_dC[ic] = 0.0; - } - UPWIND scheme; scheme.template getUpwindDirectionCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase, ip, @@ -252,19 +338,9 @@ upwindMobilityCapillary( localIndex const numPhase, hasCapPressure, upwindDir ); - localIndex const er_up = seri[upwindDir]; - localIndex const esr_up = sesri[upwindDir]; - localIndex const ei_up = sei[upwindDir]; - - if( std::fabs( phaseMob[er_up][esr_up][ei_up][ip] ) > 1e-20 ) - { - mobility = phaseMob[er_up][esr_up][ei_up][ip]; - dMobility_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP]; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMobility_dC[ic] = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dC + ic]; - } - } + //reinit + assignToZero( mobility, dMobility_dP, dMobility_dC ); + assignMobilityAndDerivatives( ip, upwindDir, seri, sesri, sei, phaseMob, dPhaseMob, mobility, dMobility_dP, dMobility_dC ); } template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND > @@ -277,7 +353,6 @@ computeFractionalFlowViscous( localIndex const numPhase, localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], - localIndex const & k_up_ppu, real64 const totFlux, real64 const totMob, real64 const (&dTotMob_dP)[numFluxSupportPoints], @@ -293,33 +368,10 @@ computeFractionalFlowViscous( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, integer const hasCapPressure, - localIndex & k_up_main, real64 & fractionalFlow, - real64 ( & dFractionalFlow_dP)[numFluxSupportPoints], - real64 ( & dFractionalFlow_dC)[numFluxSupportPoints][numComp] ) + real64 (& dFractionalFlow_dP)[numFluxSupportPoints], + real64 (& dFractionalFlow_dC)[numFluxSupportPoints][numComp] ) { - // get var to memorized the numerator mobility properly upwinded - real64 mainMob{}; - real64 dMMob_dP{}; - real64 dMMob_dC[numComp]{}; - -// real64 totMob{}; -// real64 dTotMob_dP[numFluxSupportPoints]{}; -// real64 dTotMob_dC[numFluxSupportPoints][numComp]{}; - - //reinit - //fractional flow too low to let the upstream phase flow - k_up_main = -1; //to throw error if unmodified - fractionalFlow = 0; - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dFractionalFlow_dP[ke] = 0; - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[ke][jc] = 0; - } - } - localIndex k_up; real64 mob{}; real64 dMob_dP{}; @@ -349,38 +401,13 @@ computeFractionalFlowViscous( localIndex const numPhase, dMob_dP, dMob_dC ); - k_up_main = k_up; - mainMob = mob; - dMMob_dP = dMob_dP; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMMob_dC[ic] = dMob_dC[ic]; - } - - //guard against no flow region - if( std::fabs( mainMob ) > 1e-20 ) - { - fractionalFlow = mainMob / LvArray::math::max( totMob, minTotMob ); - dFractionalFlow_dP[k_up_main] = dMMob_dP / LvArray::math::max( totMob, minTotMob ); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[k_up_main][jc] = dMMob_dC[jc] / totMob; - - } - - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dFractionalFlow_dP[ke] -= fractionalFlow * dTotMob_dP[k_up_ppu] / LvArray::math::max( totMob, minTotMob ); - - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[ke][jc] -= fractionalFlow * dTotMob_dC[k_up_ppu][jc] / LvArray::math::max( totMob, minTotMob ); - } - } - } + // reinit + assignToZero( fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC ); + computeFractionalFlowAndDerivatives( k_up, mob, dMob_dP, dMob_dC, + totMob, dTotMob_dP, dTotMob_dC, + fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC ); } - template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND > GEOS_HOST_DEVICE static void @@ -391,7 +418,6 @@ computeFractionalFlowGravity( localIndex const numPhase, localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], - localIndex const & k_up_ppu, real64 const totFlux, real64 const totMob, real64 const (&dTotMob_dP)[numFluxSupportPoints], @@ -409,31 +435,10 @@ computeFractionalFlowGravity( localIndex const numPhase, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, integer const hasCapPressure, integer const checkPhasePresenceInGravity, - localIndex & k_up_main, real64 & fractionalFlow, - real64 ( & dFractionalFlow_dP)[numFluxSupportPoints], - real64 ( & dFractionalFlow_dC)[numFluxSupportPoints][numComp] - ) + real64 (& dFractionalFlow_dP)[numFluxSupportPoints], + real64 (& dFractionalFlow_dC)[numFluxSupportPoints][numComp] ) { - // get var to memorized the numerator mobility properly upwinded - real64 mainMob{}; - real64 dMMob_dP{}; - real64 dMMob_dC[numComp]{}; - - //reinit - //fractional flow too low to let the upstream phase flow - k_up_main = -1; //to throw error if unmodified - fractionalFlow = 0; - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dFractionalFlow_dP[ke] = 0; - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[ke][jc] = 0; - } - } - - localIndex k_up; real64 mob{}; real64 dMob_dP{}; @@ -465,35 +470,11 @@ computeFractionalFlowGravity( localIndex const numPhase, dMob_dP, dMob_dC ); - k_up_main = k_up; - mainMob = mob; - dMMob_dP = dMob_dP; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMMob_dC[ic] = dMob_dC[ic]; - } - - //guard against no flow region - if( std::fabs( mainMob ) > 1e-20 ) - { - fractionalFlow = mainMob / LvArray::math::max( totMob, minTotMob ); - dFractionalFlow_dP[k_up_main] = dMMob_dP / LvArray::math::max( totMob, minTotMob ); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[k_up_main][jc] = dMMob_dC[jc] / totMob; - - } - - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dFractionalFlow_dP[ke] -= fractionalFlow * dTotMob_dP[k_up_ppu] / LvArray::math::max( totMob, minTotMob ); - - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[ke][jc] -= fractionalFlow * dTotMob_dC[k_up_ppu][jc] / LvArray::math::max( totMob, minTotMob ); - } - } - } + // reinit + assignToZero( fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC ); + computeFractionalFlowAndDerivatives( k_up, mob, dMob_dP, dMob_dC, + totMob, dTotMob_dP, dTotMob_dC, + fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC ); } template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND > @@ -506,7 +487,6 @@ computeFractionalFlowCapillary( localIndex const numPhase, localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], - localIndex const & k_up_ppu, real64 const totFlux, real64 const totMob, real64 const (&dTotMob_dP)[numFluxSupportPoints], @@ -522,30 +502,11 @@ computeFractionalFlowCapillary( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, integer const hasCapPressure, - localIndex & k_up_main, real64 & fractionalFlow, real64 ( & dFractionalFlow_dP)[numFluxSupportPoints], real64 ( & dFractionalFlow_dC)[numFluxSupportPoints][numComp] ) { - // get var to memorized the numerator mobility properly upwinded - real64 mainMob{}; - real64 dMMob_dP{}; - real64 dMMob_dC[numComp]{}; - - //reinit - //fractional flow too low to let the upstream phase flow - k_up_main = -1; //to throw error if unmodified - fractionalFlow = 0; - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dFractionalFlow_dP[ke] = 0; - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[ke][jc] = 0; - } - } - localIndex k_up; real64 mob{}; real64 dMob_dP{}; @@ -575,35 +536,11 @@ computeFractionalFlowCapillary( localIndex const numPhase, dMob_dP, dMob_dC ); - k_up_main = k_up; - mainMob = mob; - dMMob_dP = dMob_dP; - for( localIndex ic = 0; ic < numComp; ++ic ) - { - dMMob_dC[ic] = dMob_dC[ic]; - } - - //guard against no flow region - if( std::fabs( mainMob ) > 1e-20 ) - { - fractionalFlow = mainMob / LvArray::math::max( totMob, minTotMob ); - dFractionalFlow_dP[k_up_main] = dMMob_dP / LvArray::math::max( totMob, minTotMob ); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[k_up_main][jc] = dMMob_dC[jc] / LvArray::math::max( totMob, minTotMob ); - - } - - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dFractionalFlow_dP[ke] -= fractionalFlow * dTotMob_dP[k_up_ppu] / LvArray::math::max( totMob, minTotMob ); - - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dFractionalFlow_dC[ke][jc] -= fractionalFlow * dTotMob_dC[k_up_ppu][jc] / LvArray::math::max( totMob, minTotMob ); - } - } - } + // reinit + assignToZero( fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC ); + computeFractionalFlowAndDerivatives( k_up, mob, dMob_dP, dMob_dC, + totMob, dTotMob_dP, dTotMob_dC, + fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC ); } /** @@ -615,7 +552,6 @@ computeFractionalFlowCapillary( localIndex const numPhase, */ struct computePotentialViscous { - template< localIndex numComp, localIndex numFluxSupportPoints > GEOS_HOST_DEVICE static void compute( localIndex const GEOS_UNUSED_PARAM( numPhase ), @@ -641,12 +577,10 @@ 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 - } }; @@ -682,30 +616,21 @@ 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 - pot = 0.0; - for( localIndex i = 0; i < numFluxSupportPoints; ++i ) - { - dPot_dPres[i] = 0.0; - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPot_dComp[i][jc] = 0.0; - dProp_dComp[jc] = 0.0; - } - } - - calculateMeanDensity( checkPhasePresenceInGravity, ip, seri, sesri, sei, - phaseVolFrac, dCompFrac_dCompDens, - phaseMassDens, dPhaseMassDens, dProp_dComp, - densMean, dDensMean_dPres, dDensMean_dComp ); + isothermalCompositionalMultiphaseFVMKernels::helpers:: + calculateMeanDensity( ip, seri, sesri, sei, + checkPhasePresenceInGravity, + phaseVolFrac, dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, + densMean, dDensMean_dPres, dDensMean_dComp ); // compute potential difference MPFA-style for( localIndex i = 0; i < numFluxSupportPoints; ++i ) @@ -717,78 +642,13 @@ struct computePotentialGravity real64 const gravD = transmissibility[i] * gravCoef[er][esr][ei]; real64 const dGravD_dP = dTrans_dPres[i] * gravCoef[er][esr][ei]; pot += densMean * gravD; + dPot_dPres[i] += densMean * dGravD_dP; // need to add contributions from both cells the mean density depends on - for( localIndex j = 0; j < numFluxSupportPoints; ++j ) - { - dPot_dPres[j] += dDensMean_dPres[j] * gravD + densMean * dGravD_dP; - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPot_dComp[j][jc] += dDensMean_dComp[j][jc] * gravD; - } - } + addToDerivativesScaled( dDensMean_dPres, dDensMean_dComp, gravD, dPot_dPres, dPot_dComp ); } } - 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 @@ -821,10 +681,10 @@ 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 ); for( localIndex i = 0; i < numFluxSupportPoints; ++i ) { @@ -832,20 +692,22 @@ struct computePotentialCapillary localIndex const esr = sesri[i]; localIndex const ei = sei[i]; - pot += transmissibility[i] * phaseCapPressure[er][esr][ei][0][ip]; + real64 const capPres = phaseCapPressure[er][esr][ei][0][ip]; + + pot += transmissibility[i] * capPres; + + dPot_dPres[i] += dTrans_dPres[i] * capPres; + // need to add contributions from both cells for( localIndex jp = 0; jp < numPhase; ++jp ) { - real64 const dCapPressure_dS = dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp]; - dPot_dPres[i] += - transmissibility[i] * dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dP] - + dTrans_dPres[i] * phaseCapPressure[er][esr][ei][0][jp]; + + dPot_dPres[i] += transmissibility[i] * dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dP]; for( localIndex jc = 0; jc < numComp; ++jc ) { - dPot_dComp[i][jc] += transmissibility[i] * dCapPressure_dS * - dPhaseVolFrac[er][esr][ei][jp][Deriv::dC + jc]; + dPot_dComp[i][jc] += transmissibility[i] * dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dC + jc]; } } } @@ -863,7 +725,6 @@ static void computePotentialFluxesGravity( localIndex const numPhase, localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], - localIndex const & k_up_ppu, real64 const totFlux, real64 const totMob, real64 const (&dTotMob_dP)[numFluxSupportPoints], @@ -881,8 +742,6 @@ static void computePotentialFluxesGravity( localIndex const numPhase, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, localIndex const hasCapPressure, integer const checkPhasePresenceInGravity, - localIndex( &k_up), - localIndex (&k_up_o), real64 & phaseFlux, real64 (& dPhaseFlux_dP)[numFluxSupportPoints], real64 ( & dPhaseFlux_dC)[numFluxSupportPoints][numComp] ) @@ -895,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, @@ -917,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 @@ -929,7 +786,6 @@ static void computePotentialFluxesGravity( localIndex const numPhase, sei, transmissibility, dTrans_dPres, - k_up_ppu, totFlux, totMob, dTotMob_dP, @@ -947,7 +803,6 @@ static void computePotentialFluxesGravity( localIndex const numPhase, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, checkPhasePresenceInGravity, - k_up, fflow, dFflow_dP, dFflow_dC ); @@ -961,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, @@ -983,10 +837,10 @@ 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; real64 mobOther{}; real64 dMobOther_dP{}; real64 dMobOther_dC[numComp]{}; @@ -1023,31 +877,16 @@ static void computePotentialFluxesGravity( localIndex const numPhase, // Assembling gravitational flux phase-wise as \phi_{i,g} = \sum_{k\nei} \lambda_k^{up,g} f_k^{up,g} (G_i - G_k) phaseFlux -= fflow * mobOther * (pot - potOther); - dPhaseFlux_dP[k_up_o] -= fflow * dMobOther_dP * (pot - potOther); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPhaseFlux_dC[k_up_o][jc] -= fflow * dMobOther_dC[jc] * (pot - potOther); - } - //mob related part of dFflow_dP is only upstream defined but totMob related is defined everywhere - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dPhaseFlux_dP[ke] -= dFflow_dP[ke] * mobOther * (pot - potOther); + // mobOther part + UpwindHelpers::addToDerivativesScaled( dMobOther_dP, dMobOther_dC, -fflow * (pot - potOther), dPhaseFlux_dP[k_up_o], dPhaseFlux_dC[k_up_o] ); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPhaseFlux_dC[ke][jc] -= dFflow_dC[ke][jc] * mobOther * (pot - potOther); - } - } + // mob related part of dFflow_dP is only upstream defined but totMob related is defined everywhere + UpwindHelpers::addToDerivativesScaled( dFflow_dP, dFflow_dC, -mobOther * (pot - potOther), dPhaseFlux_dP, dPhaseFlux_dC ); - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dPhaseFlux_dP[ke] -= fflow * mobOther * (dPot_dP[ke] - dPotOther_dP[ke]); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPhaseFlux_dC[ke][jc] -= fflow * mobOther * (dPot_dC[ke][jc] - dPotOther_dC[ke][jc]); - } - } + // potential difference part + UpwindHelpers::addToDerivativesScaled( dPot_dP, dPot_dC, -fflow * mobOther, dPhaseFlux_dP, dPhaseFlux_dC ); + UpwindHelpers::addToDerivativesScaled( dPotOther_dP, dPotOther_dC, fflow * mobOther, dPhaseFlux_dP, dPhaseFlux_dC ); } } @@ -1062,7 +901,6 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], - localIndex const & k_up_ppu, real64 const totFlux, real64 const totMob, real64 const (&dTotMob_dP)[numFluxSupportPoints], @@ -1078,13 +916,10 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, localIndex const hasCapPressure, - localIndex( &k_up), - localIndex (&k_up_o), real64 & phaseFlux, real64 (& dPhaseFlux_dP)[numFluxSupportPoints], real64 ( & dPhaseFlux_dC)[numFluxSupportPoints][numComp] ) { - real64 fflow{}; real64 dFflow_dP[numFluxSupportPoints]{}; real64 dFflow_dC[numFluxSupportPoints][numComp]{}; @@ -1092,27 +927,25 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, real64 pot{}; real64 dPot_dP[numFluxSupportPoints]{}; real64 dPot_dC[numFluxSupportPoints][numComp]{}; - real64 dProp_dC[numComp]{}; - UpwindHelpers::computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase, - ip, - seri, - sesri, - sei, - transmissibility, - dTrans_dPres, - totFlux, - gravCoef, - dCompFrac_dCompDens, - phaseMassDens, - dPhaseMassDens, - dPhaseVolFrac, - phaseCapPressure, - dPhaseCapPressure_dPhaseVolFrac, - pot, - dPot_dP, - dPot_dC, - dProp_dC ); + computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase, + ip, + seri, + sesri, + sei, + transmissibility, + dTrans_dPres, + totFlux, + gravCoef, + dCompFrac_dCompDens, + phaseMassDens, + dPhaseMassDens, + dPhaseVolFrac, + phaseCapPressure, + dPhaseCapPressure_dPhaseVolFrac, + pot, + dPot_dP, + 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 @@ -1123,7 +956,6 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, sei, transmissibility, dTrans_dPres, - k_up_ppu, totFlux, totMob, dTotMob_dP, @@ -1139,7 +971,6 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - k_up, fflow, dFflow_dP, dFflow_dC ); @@ -1149,34 +980,32 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, { if( ip != jp ) { - 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::computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase, - jp, - seri, - sesri, - sei, - transmissibility, - dTrans_dPres, - totFlux, - gravCoef, - dCompFrac_dCompDens, - phaseMassDens, - dPhaseMassDens, - dPhaseVolFrac, - phaseCapPressure, - dPhaseCapPressure_dPhaseVolFrac, - potOther, - dPotOther_dP, - dPotOther_dC, - dPropOther_dC ); + computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase, + jp, + seri, + sesri, + sei, + transmissibility, + dTrans_dPres, + totFlux, + gravCoef, + dCompFrac_dCompDens, + phaseMassDens, + dPhaseMassDens, + dPhaseVolFrac, + phaseCapPressure, + dPhaseCapPressure_dPhaseVolFrac, + potOther, + dPotOther_dP, + dPotOther_dC ); //Eventually get the mobility of the second phase + localIndex k_up_o = -1; real64 mobOther{}; real64 dMobOther_dP{}; real64 dMobOther_dC[numComp]{}; @@ -1211,31 +1040,16 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, // Assembling gravitational flux phase-wise as \phi_{i,g} = \sum_{k\nei} \lambda_k^{up,g} f_k^{up,g} (G_i - G_k) phaseFlux -= fflow * mobOther * (pot - potOther); - dPhaseFlux_dP[k_up_o] -= fflow * dMobOther_dP * (pot - potOther); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPhaseFlux_dC[k_up_o][jc] -= fflow * dMobOther_dC[jc] * (pot - potOther); - } - //mob related part of dFflow_dP is only upstream defined but totMob related is defined everywhere - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dPhaseFlux_dP[ke] -= dFflow_dP[ke] * mobOther * (pot - potOther); + // mobOther part + UpwindHelpers::addToDerivativesScaled( dMobOther_dP, dMobOther_dC, -fflow * (pot - potOther), dPhaseFlux_dP[k_up_o], dPhaseFlux_dC[k_up_o] ); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPhaseFlux_dC[ke][jc] -= dFflow_dC[ke][jc] * mobOther * (pot - potOther); - } - } + // mob related part of dFflow_dP is only upstream defined but totMob related is defined everywhere + UpwindHelpers::addToDerivativesScaled( dFflow_dP, dFflow_dC, -mobOther * (pot - potOther), dPhaseFlux_dP, dPhaseFlux_dC ); - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dPhaseFlux_dP[ke] -= fflow * mobOther * (dPot_dP[ke] - dPotOther_dP[ke]); - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPhaseFlux_dC[ke][jc] -= fflow * mobOther * (dPot_dC[ke][jc] - dPotOther_dC[ke][jc]); - } - } + // potential difference part + UpwindHelpers::addToDerivativesScaled( dPot_dP, dPot_dC, -fflow * mobOther, dPhaseFlux_dP, dPhaseFlux_dC ); + UpwindHelpers::addToDerivativesScaled( dPotOther_dP, dPotOther_dC, fflow * mobOther, dPhaseFlux_dP, dPhaseFlux_dC ); } } } @@ -1294,7 +1108,6 @@ class UpwindScheme ) { real64 pot{}; - /// each derived concrete class has to define a computePotential method that is calling UpwindScheme::potential method with a specific /// lamda defining how to get these potentials UPWIND::template computePotentialViscous< numComp, numFluxSupportPoints >( numPhase, @@ -1316,7 +1129,6 @@ class UpwindScheme dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, pot ); - //all definition has been changed to fit pot>0 => first cell is upstream upwindDir = (pot > 0) ? 0 : 1; } @@ -1387,8 +1199,7 @@ class UpwindScheme localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], - real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in - // place + real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, @@ -1449,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; @@ -1472,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]; @@ -1517,14 +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, @@ -1543,8 +1348,7 @@ class HybridUpwind : public UpwindScheme dPhaseCapPressure_dPhaseVolFrac, potential, dPot_dP, - dPot_dC, - dProp_dC ); + dPot_dC ); } template< localIndex numComp, localIndex numFluxSupportPoints > @@ -1582,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, @@ -1605,8 +1408,7 @@ class HybridUpwind : public UpwindScheme dPhaseCapPressure_dPhaseVolFrac, potential_, dPotential_dP_, - dPotential_dC_, - dProp_dC ); + dPotential_dC_ ); } ); } @@ -1647,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, @@ -1668,8 +1470,7 @@ class HybridUpwind : public UpwindScheme dPhaseCapPressure_dPhaseVolFrac, potential_, dPotential_dP_, - dPotential_dC_, - dProp_dC ); + dPotential_dC_ ); } ); } @@ -1682,6 +1483,8 @@ struct IHUPhaseFlux using UPWIND_SCHEME = HybridUpwind; + static constexpr double minTotMob = 1e-12; + /** * @brief Form the Implicit Hybrid Upwind from pressure gradient and gravitational head * @tparam numComp number of components @@ -1704,7 +1507,6 @@ struct IHUPhaseFlux * @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 @@ -1728,21 +1530,15 @@ struct IHUPhaseFlux 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< 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] ) + real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] ) { //loop over all phases to form total velocity @@ -1754,13 +1550,9 @@ struct IHUPhaseFlux 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]; - real64 dDummy_dTrans[numComp]; + real64 dPhaseFlux_dTrans; for( integer jp = 0; jp < numPhase; ++jp ) { @@ -1771,35 +1563,44 @@ struct IHUPhaseFlux pres, gravCoef, phaseMob, dPhaseMob, phaseVolFrac, dPhaseVolFrac, - phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - k_up_ppu, potGrad, - phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, - dummy, dDummy_dP, dDummy_dC, dDummy_dTrans ); + potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans ); - totFlux += phaseFlux; - - phaseFlux = 0.; + // accumulate into total flux + UpwindHelpers::addToValueAndDerivatives( phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, + totFlux, dTotFlux_dP, dTotFlux_dC ); +/* + // old wrong way: 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.; } } + */ + + // new correct way: + // choose upstream cell + localIndex const k_up = (phaseFlux >= 0) ? 0 : 1; + // accumulate into total mobility + UpwindHelpers::assignMobilityAndDerivatives( jp, k_up, seri, sesri, sei, + phaseMob, dPhaseMob, + totMob, dTotMob_dP[k_up], dTotMob_dC[k_up] ); } - //fractional flow loop with IHU - //maybe needed to have density out for upwinding + // safeguard + totMob = LvArray::math::max( totMob, minTotMob ); + + // assign to zero + UpwindHelpers::assignToZero( phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + + // fractional flow loop with IHU + // maybe needed to have density out for upwinding // choose upstream cell // create local work arrays @@ -1821,7 +1622,6 @@ struct IHUPhaseFlux sei, trans, dTrans_dPres, - k_up_ppu, totFlux, totMob, dTotMob_dP, @@ -1837,57 +1637,27 @@ struct IHUPhaseFlux phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - k_up, fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC ); - /// Assembling the viscous flux (and derivatives) from fractional flow and total velocity as \phi_{\mu} = f_i^{up,\mu} uT viscousPhaseFlux = fractionalFlow * totFlux; - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dViscousPhaseFlux_dP[ke] += dFractionalFlow_dP[ke] * totFlux; - - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dViscousPhaseFlux_dC[ke][jc] += dFractionalFlow_dC[ke][jc] * totFlux; - } - } + UpwindHelpers::addToDerivativesScaled( dFractionalFlow_dP, dFractionalFlow_dC, + totFlux, + dViscousPhaseFlux_dP, dViscousPhaseFlux_dC ); //NON-FIXED UT -- to be canceled out if considered fixed - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dViscousPhaseFlux_dP[ke] += fractionalFlow * dTotFlux_dP[ke]; - - - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dViscousPhaseFlux_dC[ke][jc] += fractionalFlow * dTotFlux_dC[ke][jc]; - } - } - //distribute on phaseComponentFlux here - PhaseComponentFlux::compute( ip, k_up, - seri, sesri, sei, - phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, - viscousPhaseFlux, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC, - compFlux, dCompFlux_dP, dCompFlux_dC ); + UpwindHelpers::addToDerivativesScaled( dTotFlux_dP, dTotFlux_dC, + fractionalFlow, + dViscousPhaseFlux_dP, dViscousPhaseFlux_dC ); // accumulate in the flux and its derivatives - phaseFlux += viscousPhaseFlux; - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dPhaseFlux_dP[ke] += dViscousPhaseFlux_dP[ke]; - - - for( localIndex ic = 0; ic < numComp; ++ic ) - dPhaseFlux_dC[ke][ic] += dViscousPhaseFlux_dC[ke][ic]; - } + UpwindHelpers::addToValueAndDerivatives( viscousPhaseFlux, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); /// Assembling the gravitational flux (and derivatives) from fractional flow and total velocity as \phi_{g} = f_i^{up,g} uT - localIndex k_up_g = -1; - localIndex k_up_og = -1; real64 gravitationalPhaseFlux{}; real64 gravitationalPhaseFlux_dP[numFluxSupportPoints]{}; @@ -1902,7 +1672,6 @@ struct IHUPhaseFlux sei, trans, dTrans_dPres, - k_up_ppu, totFlux, totMob, dTotMob_dP, @@ -1920,37 +1689,17 @@ struct IHUPhaseFlux dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, checkPhasePresenceInGravity, - k_up_g, - k_up_og, gravitationalPhaseFlux, gravitationalPhaseFlux_dP, gravitationalPhaseFlux_dC ); - - - //distribute on phaseComponentFlux here - PhaseComponentFlux::compute( ip, k_up_g, - seri, sesri, sei, - phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, - gravitationalPhaseFlux, gravitationalPhaseFlux_dP, gravitationalPhaseFlux_dC, - compFlux, dCompFlux_dP, dCompFlux_dC ); - - //update phaseFlux from gravitational - phaseFlux += gravitationalPhaseFlux; - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dPhaseFlux_dP[ke] += gravitationalPhaseFlux_dP[ke]; - for( localIndex ic = 0; ic < numComp; ++ic ) - dPhaseFlux_dC[ke][ic] += gravitationalPhaseFlux_dC[ke][ic]; - } - + UpwindHelpers::addToValueAndDerivatives( gravitationalPhaseFlux, gravitationalPhaseFlux_dP, gravitationalPhaseFlux_dC, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); if( hasCapPressure ) { /// Assembling the capillary flux (and derivatives) from fractional flow and total velocity as \phi_{g} = f_i^{up,g} uT - localIndex k_up_pc = -1; - localIndex k_up_opc = -1; real64 capillaryPhaseFlux{}; real64 capillaryPhaseFlux_dP[numFluxSupportPoints]{}; @@ -1965,7 +1714,6 @@ struct IHUPhaseFlux sei, trans, dTrans_dPres, - k_up_ppu, totFlux, totMob, dTotMob_dP, @@ -1981,32 +1729,15 @@ struct IHUPhaseFlux phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - k_up_pc, - k_up_opc, capillaryPhaseFlux, capillaryPhaseFlux_dP, capillaryPhaseFlux_dC ); - //distribute on phaseComponentFlux here - PhaseComponentFlux::compute( ip, k_up_pc, - seri, sesri, sei, - phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, - capillaryPhaseFlux, capillaryPhaseFlux_dP, capillaryPhaseFlux_dC, - compFlux, dCompFlux_dP, dCompFlux_dC ); - - //update phaseFlux from capillary - phaseFlux += capillaryPhaseFlux; - for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) - { - dPhaseFlux_dP[ke] += capillaryPhaseFlux_dP[ke]; - for( localIndex ic = 0; ic < numComp; ++ic ) - dPhaseFlux_dC[ke][ic] += capillaryPhaseFlux_dC[ke][ic]; - - } + UpwindHelpers::addToValueAndDerivatives( capillaryPhaseFlux, capillaryPhaseFlux_dP, capillaryPhaseFlux_dC, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); }//end if cappres - } }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp index 14fc73f029..0815489aed 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp @@ -63,15 +63,11 @@ struct PPUPhaseFlux * @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 - * @param compFlux component flux - * @param dCompFlux_dP derivative of phase flux wrt pressure - * @param dCompFlux_dC derivative of phase flux wrt comp density - * @param dCompFlux_dTrans derivative of phase flux wrt transmissibility + * @param dPhaseFlux_dTrans derivative of phase flux wrt transmissibility */ template< integer numComp, integer numFluxSupportPoints > GEOS_HOST_DEVICE @@ -91,23 +87,27 @@ struct PPUPhaseFlux 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< 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], - real64 ( & dCompFlux_dTrans )[numComp] ) + real64 & dPhaseFlux_dTrans ) { + // assign to zero + for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) + { + dPhaseFlux_dP[ke] = 0; + for( integer jc = 0; jc < numComp; ++jc ) + { + dPhaseFlux_dC[ke][jc] = 0; + } + } + real64 dPotGrad_dTrans = 0; real64 dPresGrad_dP[numFluxSupportPoints]{}; real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; @@ -124,8 +124,7 @@ struct PPUPhaseFlux // *** upwinding *** // choose upstream cell - k_up = (potGrad >= 0) ? 0 : 1; - + localIndex const 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]; @@ -135,35 +134,27 @@ struct PPUPhaseFlux // compute phase flux using upwind mobility phaseFlux = mobility * potGrad; - real64 const dPhaseFlux_dTrans = mobility * dPotGrad_dTrans; + dPhaseFlux_dTrans = mobility * dPotGrad_dTrans; // pressure gradient depends on all points in the stencil for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { - dPhaseFlux_dP[ke] += dPresGrad_dP[ke] - dGravHead_dP[ke]; - dPhaseFlux_dP[ke] *= mobility; + dPhaseFlux_dP[ke] += mobility * (dPresGrad_dP[ke] - dGravHead_dP[ke]); for( integer jc = 0; jc < numComp; ++jc ) { - dPhaseFlux_dC[ke][jc] += dPresGrad_dC[ke][jc] - dGravHead_dC[ke][jc]; - dPhaseFlux_dC[ke][jc] *= mobility; + dPhaseFlux_dC[ke][jc] += mobility * (dPresGrad_dC[ke][jc] - dGravHead_dC[ke][jc]); } } real64 const dMob_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP]; - arraySlice1d< real64 const, compflow::USD_PHASE_DC - 2 > dPhaseMobSub = - dPhaseMob[er_up][esr_up][ei_up][ip]; + arraySlice1d< real64 const, compflow::USD_PHASE_DC - 2 > dMob_dC = dPhaseMob[er_up][esr_up][ei_up][ip]; // add contribution from upstream cell mobility derivatives dPhaseFlux_dP[k_up] += dMob_dP * potGrad; for( integer jc = 0; jc < numComp; ++jc ) { - dPhaseFlux_dC[k_up][jc] += dPhaseMobSub[Deriv::dC+jc] * potGrad; + dPhaseFlux_dC[k_up][jc] += dMob_dC[Deriv::dC+jc] * potGrad; } - - //distribute on phaseComponentFlux here - PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei, phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, - phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans, compFlux, dCompFlux_dP, dCompFlux_dC, dCompFlux_dTrans ); - } }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp index f9ca248b4d..a3e31acb83 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp @@ -80,22 +80,22 @@ struct PotGrad } } - // create local work arrays - real64 densMean = 0.0; - real64 dDensMean_dP[numFluxSupportPoints]{}; - real64 dDensMean_dC[numFluxSupportPoints][numComp]{}; - real64 presGrad = 0.0; real64 dPresGrad_dTrans = 0.0; real64 gravHead = 0.0; real64 dGravHead_dTrans = 0.0; real64 dCapPressure_dC[numComp]{}; - calculateMeanDensity( checkPhasePresenceInGravity, ip, - seri, sesri, sei, - phaseVolFrac, dCompFrac_dCompDens, - phaseMassDens, dPhaseMassDens, - densMean, dDensMean_dP, dDensMean_dC ); + // 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 ); /// compute the TPFA potential difference for( integer i = 0; i < numFluxSupportPoints; i++ ) @@ -164,68 +164,6 @@ struct PotGrad dPotGrad_dTrans = dPresGrad_dTrans - dGravHead_dTrans; } - 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 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp index 74332aaf86..a065a3e132 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp @@ -24,6 +24,7 @@ #include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" #include "constitutive/relativePermeability/RelativePermeabilityFields.hpp" + namespace geos { diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp index 33b07aaaaa..8be7d96cb2 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp @@ -21,7 +21,7 @@ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP #include "physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp" -//#include "physicsSolvers/fluidFlow/FluxKernelsHelper.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp" namespace geos { @@ -221,9 +221,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl real64 phaseFlux = 0.0; real64 dPhaseFlux_dP[numFluxSupportPoints]{}; real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{}; - - localIndex k_up = -1; - + real64 dPhaseFlux_dTrans = 0.0; isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux::compute< numComp, numFluxSupportPoints > ( m_numPhases, @@ -236,27 +234,31 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl m_pres, m_gravCoef, m_phaseMob, m_dPhaseMob, - m_phaseVolFrac, - m_dPhaseVolFrac, - m_phaseCompFrac, m_dPhaseCompFrac, + m_phaseVolFrac, m_dPhaseVolFrac, 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, - dCompFlux_dTrans ); + dPhaseFlux_dTrans ); + + // choose upstream cell for composition upwinding + localIndex const k_up = (phaseFlux >= 0) ? 0 : 1; + + // distribute on phaseComponentFlux here + isothermalCompositionalMultiphaseFVMKernelUtilities::PhaseComponentFlux:: + compute( ip, k_up, seri, sesri, sei, + m_phaseCompFrac, m_dPhaseCompFrac, m_dCompFrac_dCompDens, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans, + compFlux, dCompFlux_dP, dCompFlux_dC, dCompFlux_dTrans ); // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives // possible use: assemble the derivatives wrt temperature, and the flux term of the energy equation for this phase compFluxKernelOp( ip, k, seri, sesri, sei, connectionIndex, k_up, seri[k_up], sesri[k_up], sei[k_up], potGrad, - phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans ); } // loop over phases @@ -295,7 +297,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl stack.dFlux_dAperture[ic][k[1]][k[0]] -= dFlux_dAper[0]; stack.dFlux_dAperture[ic][k[1]][k[1]] -= dFlux_dAper[1]; } -// + connectionIndex++; } }