diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt index 14d4bde4ec7..0f172acafbe 100644 --- a/src/coreComponents/physicsSolvers/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/CMakeLists.txt @@ -84,7 +84,6 @@ set( physicsSolvers_sources PhysicsSolverManager.cpp SolverBase.cpp fluidFlow/CompositionalMultiphaseBase.cpp - fluidFlow/CompositionalMultiphaseBaseKernels.cpp fluidFlow/CompositionalMultiphaseFVM.cpp fluidFlow/CompositionalMultiphaseFVMKernels.cpp fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -93,7 +92,6 @@ set( physicsSolvers_sources fluidFlow/proppantTransport/ProppantTransport.cpp fluidFlow/proppantTransport/ProppantTransportKernels.cpp fluidFlow/SinglePhaseBase.cpp - fluidFlow/SinglePhaseBaseKernels.cpp fluidFlow/SinglePhaseFVM.cpp fluidFlow/SinglePhaseHybridFVM.cpp fluidFlow/SinglePhaseProppantBase.cpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 82df06e9f3a..80548de9d38 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -353,26 +353,10 @@ void CompositionalMultiphaseBase::updateComponentFraction( ObjectManagerBase & d { GEOSX_MARK_FUNCTION; - // outputs + ComponentFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + dataGroup ); - arrayView2d< real64, compflow::USD_COMP > const & compFrac = dataGroup.getExtrinsicData< extrinsicMeshData::flow::globalCompFraction >(); - - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity >(); - - // inputs - arrayView2d< real64 const, compflow::USD_COMP > const compDens = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::globalCompDensity >(); - - arrayView2d< real64 const, compflow::USD_COMP > const dCompDens = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::deltaGlobalCompDensity >(); - - KernelLaunchSelector1< ComponentFractionKernel >( m_numComponents, - dataGroup.size(), - compDens, - dCompDens, - compFrac, - dCompFrac_dCompDens ); } void CompositionalMultiphaseBase::updatePhaseVolumeFraction( ObjectManagerBase & dataGroup, @@ -380,52 +364,14 @@ void CompositionalMultiphaseBase::updatePhaseVolumeFraction( ObjectManagerBase & { GEOSX_MARK_FUNCTION; - // outputs - - arrayView2d< real64, compflow::USD_PHASE > const phaseVolFrac = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFraction >(); - - arrayView2d< real64, compflow::USD_PHASE > const dPhaseVolFrac_dPres = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dPressure >(); - - arrayView3d< real64, compflow::USD_PHASE_DC > const dPhaseVolFrac_dComp = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dGlobalCompDensity >(); - - // inputs - - arrayView3d< real64 const, compflow::USD_COMP_DC > const dCompFrac_dCompDens = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity >(); - - arrayView2d< real64 const, compflow::USD_COMP > const compDens = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::globalCompDensity >(); - - arrayView2d< real64 const, compflow::USD_COMP > const dCompDens = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::deltaGlobalCompDensity >(); - MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, m_fluidModelNames[targetIndex] ); - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseFrac = fluid.phaseFraction(); - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseFrac_dPres = fluid.dPhaseFraction_dPressure(); - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseFrac_dComp = fluid.dPhaseFraction_dGlobalCompFraction(); - - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens = fluid.phaseDensity(); - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres = fluid.dPhaseDensity_dPressure(); - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp = fluid.dPhaseDensity_dGlobalCompFraction(); - - KernelLaunchSelector2< PhaseVolumeFractionKernel >( m_numComponents, m_numPhases, - dataGroup.size(), - compDens, - dCompDens, - dCompFrac_dCompDens, - phaseDens, - dPhaseDens_dPres, - dPhaseDens_dComp, - phaseFrac, - dPhaseFrac_dPres, - dPhaseFrac_dComp, - phaseVolFrac, - dPhaseVolFrac_dPres, - dPhaseVolFrac_dComp ); + PhaseVolumeFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid ); + } void CompositionalMultiphaseBase::updateFluidModel( ObjectManagerBase & dataGroup, localIndex const targetIndex ) const @@ -998,19 +944,16 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidModelNames()[targetIndex] ); CoupledSolidBase const & solid = getConstitutiveModel< CoupledSolidBase >( subRegion, m_solidModelNames[targetIndex] ); - bool const isIsothermal = true; ElementBasedAssemblyKernelFactory:: - createAndLaunch< parallelDevicePolicy<> > - ( isIsothermal, - m_numComponents, - m_numPhases, - dofManager.rankOffset(), - dofKey, - subRegion, - fluid, - solid, - localMatrix, - localRhs ); + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.cpp deleted file mode 100644 index 26a0215019e..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.cpp +++ /dev/null @@ -1,318 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2018-2020 TotalEnergies - * Copyright (c) 2019- GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file CompositionalMultiphaseBaseKernels.cpp - */ - -#include "CompositionalMultiphaseBaseKernels.hpp" -#include "CompositionalMultiphaseUtilities.hpp" - -namespace geosx -{ - -namespace CompositionalMultiphaseBaseKernels -{ - -/******************************** ComponentFractionKernel ********************************/ - -template< localIndex NC > -GEOSX_HOST_DEVICE -void -ComponentFractionKernel:: - compute( arraySlice1d< real64 const, compflow::USD_COMP - 1 > const compDens, - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const dCompDens, - arraySlice1d< real64, compflow::USD_COMP - 1 > const compFrac, - arraySlice2d< real64, compflow::USD_COMP_DC - 1 > const dCompFrac_dCompDens ) -{ - real64 totalDensity = 0.0; - - for( localIndex ic = 0; ic < NC; ++ic ) - { - totalDensity += compDens[ic] + dCompDens[ic]; - } - - real64 const totalDensityInv = 1.0 / totalDensity; - - for( localIndex ic = 0; ic < NC; ++ic ) - { - compFrac[ic] = (compDens[ic] + dCompDens[ic]) * totalDensityInv; - for( localIndex jc = 0; jc < NC; ++jc ) - { - dCompFrac_dCompDens[ic][jc] = -compFrac[ic] * totalDensityInv; - } - dCompFrac_dCompDens[ic][ic] += totalDensityInv; - } -} - -template< localIndex NC > -void -ComponentFractionKernel:: - launch( localIndex const size, - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, - arrayView2d< real64, compflow::USD_COMP > const & compFrac, - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens ) -{ - forAll< parallelDevicePolicy<> >( size, [=] GEOSX_HOST_DEVICE ( localIndex const a ) - { - compute< NC >( compDens[a], - dCompDens[a], - compFrac[a], - dCompFrac_dCompDens[a] ); - } ); -} - -template< localIndex NC > -void -ComponentFractionKernel:: - launch( SortedArrayView< localIndex const > const & targetSet, - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, - arrayView2d< real64, compflow::USD_COMP > const & compFrac, - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens ) -{ - forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const i ) - { - localIndex const a = targetSet[ i ]; - compute< NC >( compDens[a], - dCompDens[a], - compFrac[a], - dCompFrac_dCompDens[a] ); - } ); -} - -#define INST_ComponentFractionKernel( NC ) \ - template \ - void ComponentFractionKernel:: \ - launch< NC >( localIndex const size, \ - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, \ - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, \ - arrayView2d< real64, compflow::USD_COMP > const & compFrac, \ - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens ); \ - template \ - void ComponentFractionKernel:: \ - launch< NC >( SortedArrayView< localIndex const > const & targetSet, \ - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, \ - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, \ - arrayView2d< real64, compflow::USD_COMP > const & compFrac, \ - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens ) - -INST_ComponentFractionKernel( 1 ); -INST_ComponentFractionKernel( 2 ); -INST_ComponentFractionKernel( 3 ); -INST_ComponentFractionKernel( 4 ); -INST_ComponentFractionKernel( 5 ); - -#undef INST_ComponentFractionKernel - -/******************************** PhaseVolumeFractionKernel ********************************/ - -template< localIndex NC, localIndex NP > -GEOSX_HOST_DEVICE -void -PhaseVolumeFractionKernel:: - compute( arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & compDens, - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & dCompDens, - arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dCompFrac_dCompDens, - arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & phaseDens, - arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & dPhaseDens_dPres, - arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const & dPhaseDens_dComp, - arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & phaseFrac, - arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & dPhaseFrac_dPres, - arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const & dPhaseFrac_dComp, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & phaseVolFrac, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & dPhaseVolFrac_dPres, - arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const & dPhaseVolFrac_dComp ) -{ - real64 work[NC]; - - // compute total density from component partial densities - real64 totalDensity = 0.0; - real64 const dTotalDens_dCompDens = 1.0; - for( localIndex ic = 0; ic < NC; ++ic ) - { - totalDensity += compDens[ic] + dCompDens[ic]; - } - - for( localIndex ip = 0; ip < NP; ++ip ) - { - - // set the saturation to zero if the phase is absent - bool const phaseExists = (phaseFrac[ip] > 0); - if( !phaseExists ) - { - phaseVolFrac[ip] = 0.; - dPhaseVolFrac_dPres[ip] = 0.; - for( localIndex jc = 0; jc < NC; ++jc ) - { - dPhaseVolFrac_dComp[ip][jc] = 0.; - } - continue; - } - - // Expression for volume fractions: S_p = (nu_p / rho_p) * rho_t - real64 const phaseDensInv = 1.0 / phaseDens[ip]; - - // compute saturation and derivatives except multiplying by the total density - phaseVolFrac[ip] = phaseFrac[ip] * phaseDensInv; - - dPhaseVolFrac_dPres[ip] = - (dPhaseFrac_dPres[ip] - phaseVolFrac[ip] * dPhaseDens_dPres[ip]) * phaseDensInv; - - for( localIndex jc = 0; jc < NC; ++jc ) - { - dPhaseVolFrac_dComp[ip][jc] = - (dPhaseFrac_dComp[ip][jc] - phaseVolFrac[ip] * dPhaseDens_dComp[ip][jc]) * phaseDensInv; - } - - // apply chain rule to convert derivatives from global component fractions to densities - applyChainRuleInPlace( NC, dCompFrac_dCompDens, dPhaseVolFrac_dComp[ip], work ); - - // now finalize the computation by multiplying by total density - for( localIndex jc = 0; jc < NC; ++jc ) - { - dPhaseVolFrac_dComp[ip][jc] *= totalDensity; - dPhaseVolFrac_dComp[ip][jc] += phaseVolFrac[ip] * dTotalDens_dCompDens; - } - - phaseVolFrac[ip] *= totalDensity; - dPhaseVolFrac_dPres[ip] *= totalDensity; - } -} - -template< localIndex NC, localIndex NP > -void PhaseVolumeFractionKernel:: - launch( localIndex const size, - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & phaseDens, - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & dPhaseDens_dPres, - arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & phaseFrac, - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & dPhaseFrac_dPres, - arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > const & dPhaseFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp ) -{ - forAll< parallelDevicePolicy<> >( size, [=] GEOSX_HOST_DEVICE ( localIndex const a ) - { - compute< NC, NP >( compDens[a], - dCompDens[a], - dCompFrac_dCompDens[a], - phaseDens[a][0], - dPhaseDens_dPres[a][0], - dPhaseDens_dComp[a][0], - phaseFrac[a][0], - dPhaseFrac_dPres[a][0], - dPhaseFrac_dComp[a][0], - phaseVolFrac[a], - dPhaseVolFrac_dPres[a], - dPhaseVolFrac_dComp[a] ); - } ); -} - -template< localIndex NC, localIndex NP > -void PhaseVolumeFractionKernel:: - launch( SortedArrayView< localIndex const > const & targetSet, - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & phaseDens, - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & dPhaseDens_dPres, - arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & phaseFrac, - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & dPhaseFrac_dPres, - arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > const & dPhaseFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp ) -{ - forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const i ) - { - localIndex const a = targetSet[ i ]; - compute< NC, NP >( compDens[a], - dCompDens[a], - dCompFrac_dCompDens[a], - phaseDens[a][0], - dPhaseDens_dPres[a][0], - dPhaseDens_dComp[a][0], - phaseFrac[a][0], - dPhaseFrac_dPres[a][0], - dPhaseFrac_dComp[a][0], - phaseVolFrac[a], - dPhaseVolFrac_dPres[a], - dPhaseVolFrac_dComp[a] ); - } ); -} - -#define INST_PhaseVolumeFractionKernel( NC, NP ) \ - template \ - void \ - PhaseVolumeFractionKernel:: \ - launch< NC, NP >( localIndex const size, \ - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, \ - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, \ - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, \ - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & phaseDens, \ - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & dPhaseDens_dPres, \ - arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, \ - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & phaseFrac, \ - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & dPhaseFrac_dPres, \ - arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > const & dPhaseFrac_dComp, \ - arrayView2d< real64, compflow::USD_PHASE > const & phaseVolFrac, \ - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, \ - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp ); \ - template \ - void \ - PhaseVolumeFractionKernel:: \ - launch< NC, NP >( SortedArrayView< localIndex const > const & targetSet, \ - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, \ - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, \ - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, \ - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & phaseDens, \ - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & dPhaseDens_dPres, \ - arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, \ - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & phaseFrac, \ - arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const & dPhaseFrac_dPres, \ - arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > const & dPhaseFrac_dComp, \ - arrayView2d< real64, compflow::USD_PHASE > const & phaseVolFrac, \ - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, \ - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp ) - -INST_PhaseVolumeFractionKernel( 1, 1 ); -INST_PhaseVolumeFractionKernel( 2, 1 ); -INST_PhaseVolumeFractionKernel( 3, 1 ); -INST_PhaseVolumeFractionKernel( 4, 1 ); -INST_PhaseVolumeFractionKernel( 5, 1 ); - -INST_PhaseVolumeFractionKernel( 1, 2 ); -INST_PhaseVolumeFractionKernel( 2, 2 ); -INST_PhaseVolumeFractionKernel( 3, 2 ); -INST_PhaseVolumeFractionKernel( 4, 2 ); -INST_PhaseVolumeFractionKernel( 5, 2 ); - -INST_PhaseVolumeFractionKernel( 1, 3 ); -INST_PhaseVolumeFractionKernel( 2, 3 ); -INST_PhaseVolumeFractionKernel( 3, 3 ); -INST_PhaseVolumeFractionKernel( 4, 3 ); -INST_PhaseVolumeFractionKernel( 5, 3 ); - -#undef INST_PhaseVolumeFractionKernel - -} // namespace CompositionalMultiphaseBaseKernels - -} // namespace geosx diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp index 535e5e86dbf..84105c61bb1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp @@ -24,12 +24,11 @@ #include "common/GEOS_RAJA_Interface.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" #include "constitutive/fluid/MultiFluidBase.hpp" -#include "constitutive/fluid/layouts.hpp" #include "functions/TableFunction.hpp" +#include "mesh/ElementSubRegionBase.hpp" +#include "mesh/ObjectManagerBase.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp" -#include "mesh/ElementRegionManager.hpp" - namespace geosx { @@ -41,94 +40,407 @@ using namespace constitutive; static constexpr real64 minDensForDivision = 1e-10; -/******************************** ComponentFractionKernel ********************************/ +/******************************** PropertyKernelBase ********************************/ /** - * @brief Functions to compute component fractions from global component densities (mass or molar) + * @brief Internal struct to provide no-op defaults used in the inclusion + * of lambda functions into kernel component functions. + * @struct NoOpFuncs */ -struct ComponentFractionKernel +struct NoOpFunc { - template< localIndex NC > + template< typename ... Ts > GEOSX_HOST_DEVICE - static void - compute( arraySlice1d< real64 const, compflow::USD_COMP - 1 > compDens, - arraySlice1d< real64 const, compflow::USD_COMP - 1 > dCompDens, - arraySlice1d< real64, compflow::USD_COMP - 1 > compFrac, - arraySlice2d< real64, compflow::USD_COMP_DC - 1 > dCompFrac_dCompDens ); + constexpr void + operator()( Ts && ... ) const {} +}; - template< localIndex NC > +/** + * @class PropertyKernelBase + * @tparam NUM_COMP number of fluid components + * @brief Define the base interface for the property update kernels + */ +template< integer NUM_COMP > +class PropertyKernelBase +{ +public: + + /// Compile time value for the number of components + static constexpr integer numComp = NUM_COMP; + + /** + * @brief Performs the kernel launch + * @tparam POLICY the policy used in the RAJA kernels + * @tparam KERNEL_TYPE the kernel type + * @param[in] numElems the number of elements + * @param[inout] kernelComponent the kernel component providing access to the compute function + */ + template< typename POLICY, typename KERNEL_TYPE > static void - launch( localIndex const size, - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, - arrayView2d< real64, compflow::USD_COMP > const & compFrac, - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens ); + launch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + forAll< POLICY >( numElems, [=] GEOSX_HOST_DEVICE ( localIndex const ei ) + { + kernelComponent.compute( ei ); + } ); + } - template< localIndex NC > + /** + * @brief Performs the kernel launch on a sorted array + * @tparam POLICY the policy used in the RAJA kernels + * @tparam KERNEL_TYPE the kernel type + * @param[in] targetSet the indices of the elements in which we compute the property + * @param[inout] kernelComponent the kernel component providing access to the compute function + */ + template< typename POLICY, typename KERNEL_TYPE > static void launch( SortedArrayView< localIndex const > const & targetSet, - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, - arrayView2d< real64, compflow::USD_COMP > const & compFrac, - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens ); + KERNEL_TYPE const & kernelComponent ) + { + forAll< POLICY >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const i ) + { + localIndex const ei = targetSet[ i ]; + kernelComponent.compute( ei ); + } ); + } + }; -/******************************** PhaseVolumeFractionKernel ********************************/ +namespace internal +{ + +template< typename T, typename LAMBDA > +void kernelLaunchSelectorCompSwitch( T value, LAMBDA && lambda ) +{ + static_assert( std::is_integral< T >::value, "kernelLaunchSelectorCompSwitch: type should be integral" ); + + switch( value ) + { + case 1: + { lambda( std::integral_constant< T, 1 >() ); return; } + case 2: + { lambda( std::integral_constant< T, 2 >() ); return; } + case 3: + { lambda( std::integral_constant< T, 3 >() ); return; } + case 4: + { lambda( std::integral_constant< T, 4 >() ); return; } + case 5: + { lambda( std::integral_constant< T, 5 >() ); return; } + default: + { GEOSX_ERROR( "Unsupported number of components: " << value ); } + } +} + +} // namespace internal + + +/******************************** ComponentFractionKernel ********************************/ /** - * @brief Functions to compute phase volume fractions (saturations) and derivatives + * @class ComponentFractionKernel + * @tparam NUM_COMP number of fluid components + * @brief Define the interface for the update kernel in charge of computing the phase volume fractions */ -struct PhaseVolumeFractionKernel +template< integer NUM_COMP > +class ComponentFractionKernel : public PropertyKernelBase< NUM_COMP > { - template< localIndex NC, localIndex NP > +public: + + using Base = PropertyKernelBase< NUM_COMP >; + using Base::numComp; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + ComponentFractionKernel( ObjectManagerBase & subRegion ) + : Base(), + m_compDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::globalCompDensity >() ), + m_dCompDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::deltaGlobalCompDensity >() ), + m_compFrac( subRegion.getExtrinsicData< extrinsicMeshData::flow::globalCompFraction >() ), + m_dCompFrac_dCompDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity >() ) + {} + + /** + * @brief Compute the phase volume fractions in an element + * @tparam FUNC the type of the function that can be used to customize the kernel + * @param[in] ei the element index + * @param[in] phaseVolFractionKernelOp the function used to customize the kernel + */ + template< typename FUNC = NoOpFunc > GEOSX_HOST_DEVICE + void compute( localIndex const ei, + FUNC && compFractionKernelOp = NoOpFunc{} ) const + { + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const compDens = m_compDens[ei]; + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const dCompDens = m_dCompDens[ei]; + arraySlice1d< real64, compflow::USD_COMP - 1 > const compFrac = m_compFrac[ei]; + arraySlice2d< real64, compflow::USD_COMP_DC - 1 > const dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei]; + + real64 totalDensity = 0.0; + + for( integer ic = 0; ic < numComp; ++ic ) + { + totalDensity += compDens[ic] + dCompDens[ic]; + } + + real64 const totalDensityInv = 1.0 / totalDensity; + + for( integer ic = 0; ic < numComp; ++ic ) + { + compFrac[ic] = (compDens[ic] + dCompDens[ic]) * totalDensityInv; + for( integer jc = 0; jc < numComp; ++jc ) + { + dCompFrac_dCompDens[ic][jc] = -compFrac[ic] * totalDensityInv; + } + dCompFrac_dCompDens[ic][ic] += totalDensityInv; + } + + compFractionKernelOp( compFrac, dCompFrac_dCompDens ); + } + +protected: + + // inputs + + // Views on component densities + arrayView2d< real64 const, compflow::USD_COMP > m_compDens; + arrayView2d< real64 const, compflow::USD_COMP > m_dCompDens; + + // outputs + + // Views on component fraction + arrayView2d< real64, compflow::USD_COMP > m_compFrac; + arrayView3d< real64, compflow::USD_COMP_DC > m_dCompFrac_dCompDens; + +}; + +/** + * @class ComponentFractionKernelFactory + */ +class ComponentFractionKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numComp the number of fluid components + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + template< typename POLICY > static void - compute( arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & compDens, - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & dCompDens, - arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dCompFrac_dCompDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseDens_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseDens_dComp, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseFrac, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseFrac_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseFrac_dComp, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & phaseVolFrac, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & dPhaseVolFrac_dPres, - arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const & dPhaseVolFrac_dComp ); - - template< localIndex NC, localIndex NP > - static void - launch( localIndex const size, - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseFrac, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseFrac_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp ); - - template< localIndex NC, localIndex NP > - static void - launch( SortedArrayView< localIndex const > const & targetSet, - arrayView2d< real64 const, compflow::USD_COMP > const & compDens, - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseFrac, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseFrac_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp ); + createAndLaunch( integer const numComp, + ObjectManagerBase & subRegion ) + { + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + ComponentFractionKernel< NUM_COMP > kernel( subRegion ); + ComponentFractionKernel< NUM_COMP >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + }; +/******************************** PhaseVolumeFractionKernel ********************************/ + +/** + * @class PhaseVolumeFractionKernel + * @tparam NUM_COMP number of fluid components + * @tparam NUM_PHASE number of fluid phases + * @brief Define the interface for the property kernel in charge of computing the phase volume fractions + */ +template< integer NUM_COMP, integer NUM_PHASE > +class PhaseVolumeFractionKernel : public PropertyKernelBase< NUM_COMP > +{ +public: + + using Base = PropertyKernelBase< NUM_COMP >; + using Base::numComp; + + /// Compile time value for the number of phases + static constexpr integer numPhase = NUM_PHASE; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + PhaseVolumeFractionKernel( ObjectManagerBase & subRegion, + MultiFluidBase const & fluid ) + : Base(), + m_phaseVolFrac( subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFraction >() ), + m_dPhaseVolFrac_dPres( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dPressure >() ), + m_dPhaseVolFrac_dComp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dGlobalCompDensity >() ), + m_compDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::globalCompDensity >() ), + m_dCompDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::deltaGlobalCompDensity >() ), + m_dCompFrac_dCompDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity >() ), + m_phaseFrac( fluid.phaseFraction() ), + m_dPhaseFrac_dPres( fluid.dPhaseFraction_dPressure() ), + m_dPhaseFrac_dComp( fluid.dPhaseFraction_dGlobalCompFraction() ), + m_phaseDens( fluid.phaseDensity() ), + m_dPhaseDens_dPres( fluid.dPhaseDensity_dPressure() ), + m_dPhaseDens_dComp( fluid.dPhaseDensity_dGlobalCompFraction() ) + {} + + /** + * @brief Compute the phase volume fractions in an element + * @tparam FUNC the type of the function that can be used to customize the kernel + * @param[in] ei the element index + * @param[in] phaseVolFractionKernelOp the function used to customize the kernel + */ + template< typename FUNC = NoOpFunc > + GEOSX_HOST_DEVICE + void compute( localIndex const ei, + FUNC && phaseVolFractionKernelOp = NoOpFunc{} ) const + { + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const compDens = m_compDens[ei]; + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const dCompDens = m_dCompDens[ei]; + arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const phaseDens = m_phaseDens[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const dPhaseDens_dPres = m_dPhaseDens_dPres[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseDens_dComp = m_dPhaseDens_dComp[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const phaseFrac = m_phaseFrac[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const dPhaseFrac_dPres = m_dPhaseFrac_dPres[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseFrac_dComp = m_dPhaseFrac_dComp[ei][0]; + arraySlice1d< real64, compflow::USD_PHASE - 1 > const phaseVolFrac = m_phaseVolFrac[ei]; + arraySlice1d< real64, compflow::USD_PHASE - 1 > const dPhaseVolFrac_dPres = m_dPhaseVolFrac_dPres[ei]; + arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const dPhaseVolFrac_dComp = m_dPhaseVolFrac_dComp[ei]; + + real64 work[numComp]{}; + + // compute total density from component partial densities + real64 totalDensity = 0.0; + real64 const dTotalDens_dCompDens = 1.0; + for( integer ic = 0; ic < numComp; ++ic ) + { + totalDensity += compDens[ic] + dCompDens[ic]; + } + + for( integer ip = 0; ip < numPhase; ++ip ) + { + + // set the saturation to zero if the phase is absent + bool const phaseExists = (phaseFrac[ip] > 0); + if( !phaseExists ) + { + phaseVolFrac[ip] = 0.; + dPhaseVolFrac_dPres[ip] = 0.; + for( integer jc = 0; jc < numComp; ++jc ) + { + dPhaseVolFrac_dComp[ip][jc] = 0.; + } + continue; + } + + // Expression for volume fractions: S_p = (nu_p / rho_p) * rho_t + real64 const phaseDensInv = 1.0 / phaseDens[ip]; + + // compute saturation and derivatives except multiplying by the total density + phaseVolFrac[ip] = phaseFrac[ip] * phaseDensInv; + + dPhaseVolFrac_dPres[ip] = + (dPhaseFrac_dPres[ip] - phaseVolFrac[ip] * dPhaseDens_dPres[ip]) * phaseDensInv; + + for( integer jc = 0; jc < numComp; ++jc ) + { + dPhaseVolFrac_dComp[ip][jc] = + (dPhaseFrac_dComp[ip][jc] - phaseVolFrac[ip] * dPhaseDens_dComp[ip][jc]) * phaseDensInv; + } + + // apply chain rule to convert derivatives from global component fractions to densities + applyChainRuleInPlace( numComp, dCompFrac_dCompDens, dPhaseVolFrac_dComp[ip], work ); + + // call the lambda in the phase loop to allow the reuse of the phaseVolFrac and totalDensity + // possible use: assemble the derivatives wrt temperature + phaseVolFractionKernelOp( ip, phaseVolFrac[ip], totalDensity ); + + // now finalize the computation by multiplying by total density + for( integer jc = 0; jc < numComp; ++jc ) + { + dPhaseVolFrac_dComp[ip][jc] *= totalDensity; + dPhaseVolFrac_dComp[ip][jc] += phaseVolFrac[ip] * dTotalDens_dCompDens; + } + + phaseVolFrac[ip] *= totalDensity; + dPhaseVolFrac_dPres[ip] *= totalDensity; + } + } + +protected: + + // outputs + + /// Views on phase volume fractions + arrayView2d< real64, compflow::USD_PHASE > m_phaseVolFrac; + arrayView2d< real64, compflow::USD_PHASE > m_dPhaseVolFrac_dPres; + arrayView3d< real64, compflow::USD_PHASE_DC > m_dPhaseVolFrac_dComp; + + // inputs + + /// Views on component densities + arrayView2d< real64 const, compflow::USD_COMP > m_compDens; + arrayView2d< real64 const, compflow::USD_COMP > m_dCompDens; + arrayView3d< real64 const, compflow::USD_COMP_DC > m_dCompFrac_dCompDens; + + /// Views on phase fractions + arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseFrac; + arrayView3d< real64 const, multifluid::USD_PHASE > m_dPhaseFrac_dPres; + arrayView4d< real64 const, multifluid::USD_PHASE_DC > m_dPhaseFrac_dComp; + + /// Views on phase densities + arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseDens; + arrayView3d< real64 const, multifluid::USD_PHASE > m_dPhaseDens_dPres; + arrayView4d< real64 const, multifluid::USD_PHASE_DC > m_dPhaseDens_dComp; + +}; + +/** + * @class PhaseVolumeFractionKernelFactory + */ +class PhaseVolumeFractionKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numComp the number of fluid components + * @param[in] numPhase the number of fluid phases + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + template< typename POLICY > + static void + createAndLaunch( integer const numComp, + integer const numPhase, + ObjectManagerBase & subRegion, + MultiFluidBase const & fluid ) + { + if( numPhase == 2 ) + { + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseVolumeFractionKernel< NUM_COMP, 2 > kernel( subRegion, fluid ); + PhaseVolumeFractionKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseVolumeFractionKernel< NUM_COMP, 3 > kernel( subRegion, fluid ); + PhaseVolumeFractionKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + } +}; /******************************** FluidUpdateKernel ********************************/ @@ -281,40 +593,25 @@ struct CapillaryPressureUpdateKernel /******************************** ElementBasedAssemblyKernel ********************************/ -/** - * @brief Internal struct to provide no-op defaults used in the inclusion - * of lambda functions into kernel component functions. - * @struct NoOpFuncs - */ -struct NoOpFunc -{ - template< typename ... Ts > - GEOSX_HOST_DEVICE - constexpr void - operator()( Ts && ... ) const {} -}; - - - /** * @class ElementBasedAssemblyKernel * @tparam NUM_COMP number of fluid components * @tparam NUM_DOF number of degrees of freedom * @brief Define the interface for the assembly kernel in charge of accumulation and volume balance */ -template< localIndex NUM_COMP, localIndex NUM_DOF > +template< integer NUM_COMP, integer NUM_DOF > class ElementBasedAssemblyKernel { public: /// Compile time value for the number of components - static constexpr localIndex numComp = NUM_COMP; + static constexpr integer numComp = NUM_COMP; /// Compute time value for the number of degrees of freedom - static constexpr localIndex numDof = NUM_DOF; + static constexpr integer numDof = NUM_DOF; /// Compute time value for the number of equations - static constexpr localIndex numEqn = NUM_DOF; + static constexpr integer numEqn = NUM_DOF; /** * @brief Constructor @@ -620,7 +917,7 @@ class ElementBasedAssemblyKernel protected: /// Number of fluid phases - localIndex const m_numPhases; + integer const m_numPhases; /// Offset for my MPI rank globalIndex const m_rankOffset; @@ -667,34 +964,6 @@ class ElementBasedAssemblyKernel }; -namespace internal -{ - -template< typename T, typename LAMBDA > -void kernelLaunchSelectorCompSwitch( T value, LAMBDA && lambda ) -{ - static_assert( std::is_integral< T >::value, "kernelLaunchSelectorCompSwitch: type should be integral" ); - - switch( value ) - { - case 1: - { lambda( std::integral_constant< T, 1 >() ); return; } - case 2: - { lambda( std::integral_constant< T, 2 >() ); return; } - case 3: - { lambda( std::integral_constant< T, 3 >() ); return; } - case 4: - { lambda( std::integral_constant< T, 4 >() ); return; } - case 5: - { lambda( std::integral_constant< T, 5 >() ); return; } - default: - { GEOSX_ERROR( "Unsupported number of components: " << value ); } - } -} - -} // namespace internal - - /** * @class ElementBasedAssemblyKernelFactory */ @@ -705,7 +974,6 @@ class ElementBasedAssemblyKernelFactory /** * @brief Create a new kernel and launch * @tparam POLICY the policy used in the RAJA kernel - * @param[in] isIsothermal flag specifying whether the assembly is isothermal or non-isothermal * @param[in] numComps the number of fluid components * @param[in] numPhases the number of fluid phases * @param[in] rankOffset the offset of my MPI rank @@ -718,9 +986,8 @@ class ElementBasedAssemblyKernelFactory */ template< typename POLICY > static void - createAndLaunch( bool const isIsothermal, - localIndex const numComps, - localIndex const numPhases, + createAndLaunch( integer const numComps, + integer const numPhases, globalIndex const rankOffset, string const dofKey, ElementSubRegionBase const & subRegion, @@ -731,28 +998,11 @@ class ElementBasedAssemblyKernelFactory { internal::kernelLaunchSelectorCompSwitch( numComps, [&] ( auto NC ) { - if( isIsothermal ) - { - localIndex constexpr NUM_COMP = NC(); - localIndex constexpr NUM_DOF = NC()+1; - ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > - kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); - ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF >::template - launch< POLICY, ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > >( subRegion.size(), kernel ); - } - else - { - GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); - /* - TODO: uncomment and move when the thermal kernel is ready - localIndex constexpr NUM_COMP = NC(); - localIndex constexpr NUM_DOF = NC()+2; - ThermalElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > - kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); - ThermalElementBasedAssemblyKernel< NUM_COMP, NUM_DOF >::template - launch< POLICY, ThermalElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > >( subRegion.size(), kernel ); - */ - } + integer constexpr NUM_COMP = NC(); + integer constexpr NUM_DOF = NC()+1; + ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > + kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel ); } ); } @@ -766,7 +1016,7 @@ struct ResidualNormKernel template< typename POLICY, typename REDUCE_POLICY > static void launch( arrayView1d< real64 const > const & localResidual, globalIndex const rankOffset, - localIndex const numComponents, + integer const numComponents, arrayView1d< globalIndex const > const & dofNumber, arrayView1d< integer const > const & ghostRank, arrayView1d< real64 const > const & refPoro, @@ -783,7 +1033,7 @@ struct ResidualNormKernel localIndex const localRow = dofNumber[ei] - rankOffset; real64 const normalizer = totalDensOld[ei] * refPoro[ei] * volume[ei]; - for( localIndex idof = 0; idof < numComponents + 1; ++idof ) + for( integer idof = 0; idof < numComponents + 1; ++idof ) { real64 const val = localResidual[localRow + idof] / normalizer; localSum += val * val; @@ -833,7 +1083,7 @@ struct SolutionCheckKernel // will be chopped (i.e., set to zero) in ApplySystemSolution) if( !allowCompDensChopping ) { - for( localIndex ic = 0; ic < numComponents; ++ic ) + for( integer ic = 0; ic < numComponents; ++ic ) { real64 const newDens = compDens[ei][ic] + dCompDens[ei][ic] + scalingFactor * localSolution[localRow + ic + 1]; check.min( newDens >= 0.0 ); @@ -842,7 +1092,7 @@ struct SolutionCheckKernel else { real64 totalDens = 0.0; - for( localIndex ic = 0; ic < numComponents; ++ic ) + for( integer ic = 0; ic < numComponents; ++ic ) { real64 const newDens = compDens[ei][ic] + dCompDens[ei][ic] + scalingFactor * localSolution[localRow + ic + 1]; totalDens += (newDens > 0.0) ? newDens : 0.0; @@ -873,8 +1123,8 @@ struct HydrostaticPressureKernel template< typename FLUID_WRAPPER > static ReturnType - computeHydrostaticPressure( localIndex const numComps, - localIndex const numPhases, + computeHydrostaticPressure( integer const numComps, + integer const numPhases, integer const ipInit, integer const maxNumEquilIterations, real64 const & equilTolerance, @@ -905,7 +1155,7 @@ struct HydrostaticPressureKernel real64 const gravCoef = gravVector[2] * ( refElevation - newElevation ); real64 const temp = tempTableWrapper.compute( &newElevation ); - for( localIndex ic = 0; ic < numComps; ++ic ) + for( integer ic = 0; ic < numComps; ++ic ) { compFrac[0][ic] = compFracTableWrappers[ic].compute( &newElevation ); } @@ -931,7 +1181,7 @@ struct HydrostaticPressureKernel // Step 4: fixed-point iteration until convergence bool equilHasConverged = false; - for( localIndex eqIter = 0; eqIter < maxNumEquilIterations; ++eqIter ) + for( integer eqIter = 0; eqIter < maxNumEquilIterations; ++eqIter ) { // check convergence @@ -944,7 +1194,7 @@ struct HydrostaticPressureKernel // make sure that the fluid is single-phase, other we have to issue a warning (for now) // if only one phase is mobile, we are in good shape (unfortunately it is hard to access relperm from here) localIndex numberOfPhases = 0; - for( localIndex ip = 0; ip < numPhases; ++ip ) + for( integer ip = 0; ip < numPhases; ++ip ) { if( phaseFrac[0][0][ip] > MIN_FOR_PHASE_PRESENCE ) { @@ -975,7 +1225,7 @@ struct HydrostaticPressureKernel // Step 5: save the hydrostatic pressure and the corresponding density newPres = pres1; - for( localIndex ip = 0; ip < numPhases; ++ip ) + for( integer ip = 0; ip < numPhases; ++ip ) { newPhaseMassDens[ip] = phaseMassDens[0][0][ip]; } @@ -997,8 +1247,8 @@ struct HydrostaticPressureKernel template< typename FLUID_WRAPPER > static ReturnType launch( localIndex const size, - localIndex const numComps, - localIndex const numPhases, + integer const numComps, + integer const numPhases, integer const ipInit, integer const maxNumEquilIterations, real64 const equilTolerance, @@ -1028,7 +1278,7 @@ struct HydrostaticPressureKernel real64 datumTotalDens = 0.0; real64 const datumTemp = tempTableWrapper.compute( &datumElevation ); - for( localIndex ic = 0; ic < numComps; ++ic ) + for( integer ic = 0; ic < numComps; ++ic ) { datumCompFrac[0][ic] = compFracTableWrappers[ic].compute( &datumElevation ); } @@ -1058,7 +1308,7 @@ struct HydrostaticPressureKernel array2d< real64 > phaseMassDens( pressureValues.size(), numPhases ); // temporary array without permutation to compile on Lassen array1d< real64 > datumPhaseMassDensTmp( numPhases ); - for( localIndex ip = 0; ip < numPhases; ++ip ) + for( integer ip = 0; ip < numPhases; ++ip ) { datumPhaseMassDensTmp[ip] = datumPhaseMassDens[0][0][ip]; } @@ -1163,7 +1413,7 @@ struct HydrostaticPressureKernel /******************************** Kernel launch machinery ********************************/ template< typename KERNELWRAPPER, typename ... ARGS > -void KernelLaunchSelector1( localIndex const numComp, ARGS && ... args ) +void KernelLaunchSelector1( integer const numComp, ARGS && ... args ) { internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { @@ -1172,7 +1422,7 @@ void KernelLaunchSelector1( localIndex const numComp, ARGS && ... args ) } template< typename KERNELWRAPPER, typename ... ARGS > -void KernelLaunchSelector2( localIndex const numComp, localIndex const numPhase, ARGS && ... args ) +void KernelLaunchSelector2( integer const numComp, integer const numPhase, ARGS && ... args ) { // Ideally this would be inside the dispatch, but it breaks on Summit with GCC 9.1.0 and CUDA 11.0.3. if( numPhase == 2 ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 87c4582a124..3c607a9767f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -519,65 +519,17 @@ void CompositionalMultiphaseFVM::updatePhaseMobility( ObjectManagerBase & dataGr { GEOSX_MARK_FUNCTION; - // note that for convenience, the phase mobility computed here also includes phase density - - // outputs - - arrayView2d< real64, compflow::USD_PHASE > const phaseMob = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::phaseMobility >(); - - arrayView2d< real64, compflow::USD_PHASE > const dPhaseMob_dPres = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dPressure >(); - - arrayView3d< real64, compflow::USD_PHASE_DC > const dPhaseMob_dComp = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dGlobalCompDensity >(); - - // inputs - - arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFraction >(); - - arrayView2d< real64 const, compflow::USD_PHASE > const dPhaseVolFrac_dPres = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dPressure >(); - - arrayView3d< real64 const, compflow::USD_PHASE_DC > const dPhaseVolFrac_dComp = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dGlobalCompDensity >(); - - arrayView3d< real64 const, compflow::USD_COMP_DC > const dCompFrac_dCompDens = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity >(); + // note that the phase mobility computed here also includes phase density MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, m_fluidModelNames[targetIndex] ); - - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens = fluid.phaseDensity(); - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres = fluid.dPhaseDensity_dPressure(); - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp = fluid.dPhaseDensity_dGlobalCompFraction(); - - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc = fluid.phaseViscosity(); - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres = fluid.dPhaseViscosity_dPressure(); - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp = fluid.dPhaseViscosity_dGlobalCompFraction(); - RelativePermeabilityBase const & relperm = getConstitutiveModel< RelativePermeabilityBase >( dataGroup, m_relPermModelNames[targetIndex] ); - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm = relperm.phaseRelPerm(); - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac = relperm.dPhaseRelPerm_dPhaseVolFraction(); - - KernelLaunchSelector2< PhaseMobilityKernel >( m_numComponents, m_numPhases, - dataGroup.size(), - dCompFrac_dCompDens, - phaseDens, - dPhaseDens_dPres, - dPhaseDens_dComp, - phaseVisc, - dPhaseVisc_dPres, - dPhaseVisc_dComp, - phaseRelPerm, - dPhaseRelPerm_dPhaseVolFrac, - phaseVolFrac, - dPhaseVolFrac_dPres, - dPhaseVolFrac_dComp, - phaseMob, - dPhaseMob_dPres, - dPhaseMob_dComp ); + PhaseMobilityKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid, + relperm ); } void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp index b377f156ae4..46339fed0c4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp @@ -33,236 +33,13 @@ namespace geosx namespace CompositionalMultiphaseFVMKernels { -/******************************** PhaseMobilityKernel ********************************/ - -template< localIndex NC, localIndex NP > -GEOSX_HOST_DEVICE -void -PhaseMobilityKernel:: - compute( arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dCompFrac_dCompDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseDens_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseDens_dComp, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseVisc, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseVisc_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseVisc_dComp, - arraySlice1d< real64 const, relperm::USD_RELPERM - 2 > const & phaseRelPerm, - arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > const & dPhaseRelPerm_dPhaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & phaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & dPhaseVolFrac_dPres, - arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const & dPhaseVolFrac_dComp, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & phaseMob, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & dPhaseMob_dPres, - arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const & dPhaseMob_dComp ) -{ - real64 dRelPerm_dC[NC]; - real64 dDens_dC[NC]; - real64 dVisc_dC[NC]; - - for( localIndex ip = 0; ip < NP; ++ip ) - { - - // compute the phase mobility only if the phase is present - bool const phaseExists = (phaseVolFrac[ip] > 0); - if( !phaseExists ) - { - phaseMob[ip] = 0.; - dPhaseMob_dPres[ip] = 0.; - for( localIndex jc = 0; jc < NC; ++jc ) - { - dPhaseMob_dComp[ip][jc] = 0.; - } - continue; - } - - real64 const density = phaseDens[ip]; - real64 const dDens_dP = dPhaseDens_dPres[ip]; - applyChainRule( NC, dCompFrac_dCompDens, dPhaseDens_dComp[ip], dDens_dC ); - - real64 const viscosity = phaseVisc[ip]; - real64 const dVisc_dP = dPhaseVisc_dPres[ip]; - applyChainRule( NC, dCompFrac_dCompDens, dPhaseVisc_dComp[ip], dVisc_dC ); - - real64 const relPerm = phaseRelPerm[ip]; - real64 dRelPerm_dP = 0.0; - for( localIndex ic = 0; ic < NC; ++ic ) - { - dRelPerm_dC[ic] = 0.0; - } - - for( localIndex jp = 0; jp < NP; ++jp ) - { - real64 const dRelPerm_dS = dPhaseRelPerm_dPhaseVolFrac[ip][jp]; - dRelPerm_dP += dRelPerm_dS * dPhaseVolFrac_dPres[jp]; - - for( localIndex jc = 0; jc < NC; ++jc ) - { - dRelPerm_dC[jc] += dRelPerm_dS * dPhaseVolFrac_dComp[jp][jc]; - } - } - - real64 const mobility = relPerm * density / viscosity; - - phaseMob[ip] = mobility; - dPhaseMob_dPres[ip] = dRelPerm_dP * density / viscosity - + mobility * (dDens_dP / density - dVisc_dP / viscosity); - - // compositional derivatives - for( localIndex jc = 0; jc < NC; ++jc ) - { - dPhaseMob_dComp[ip][jc] = dRelPerm_dC[jc] * density / viscosity - + mobility * (dDens_dC[jc] / density - dVisc_dC[jc] / viscosity); - } - } -} - -template< localIndex NC, localIndex NP > -void PhaseMobilityKernel:: - launch( localIndex const size, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ) -{ - forAll< parallelDevicePolicy<> >( size, [=] GEOSX_HOST_DEVICE ( localIndex const a ) - { - compute< NC, NP >( dCompFrac_dCompDens[a], - phaseDens[a][0], - dPhaseDens_dPres[a][0], - dPhaseDens_dComp[a][0], - phaseVisc[a][0], - dPhaseVisc_dPres[a][0], - dPhaseVisc_dComp[a][0], - phaseRelPerm[a][0], - dPhaseRelPerm_dPhaseVolFrac[a][0], - phaseVolFrac[a], - dPhaseVolFrac_dPres[a], - dPhaseVolFrac_dComp[a], - phaseMob[a], - dPhaseMob_dPres[a], - dPhaseMob_dComp[a] ); - } ); -} - -template< localIndex NC, localIndex NP > -void PhaseMobilityKernel:: - launch( SortedArrayView< localIndex const > const & targetSet, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ) -{ - forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const i ) - { - localIndex const a = targetSet[ i ]; - compute< NC, NP >( dCompFrac_dCompDens[a], - phaseDens[a][0], - dPhaseDens_dPres[a][0], - dPhaseDens_dComp[a][0], - phaseVisc[a][0], - dPhaseVisc_dPres[a][0], - dPhaseVisc_dComp[a][0], - phaseRelPerm[a][0], - dPhaseRelPerm_dPhaseVolFrac[a][0], - phaseVolFrac[a], - dPhaseVolFrac_dPres[a], - dPhaseVolFrac_dComp[a], - phaseMob[a], - dPhaseMob_dPres[a], - dPhaseMob_dComp[a] ); - } ); -} - -#define INST_PhaseMobilityKernel( NC, NP ) \ - template \ - void \ - PhaseMobilityKernel:: \ - launch< NC, NP >( localIndex const size, \ - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres, \ - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, \ - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, \ - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, \ - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, \ - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, \ - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, \ - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, \ - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ); \ - template \ - void \ - PhaseMobilityKernel:: \ - launch< NC, NP >( SortedArrayView< localIndex const > const & targetSet, \ - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres, \ - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, \ - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, \ - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, \ - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, \ - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, \ - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, \ - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, \ - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ) - -INST_PhaseMobilityKernel( 1, 1 ); -INST_PhaseMobilityKernel( 2, 1 ); -INST_PhaseMobilityKernel( 3, 1 ); -INST_PhaseMobilityKernel( 4, 1 ); -INST_PhaseMobilityKernel( 5, 1 ); - -INST_PhaseMobilityKernel( 1, 2 ); -INST_PhaseMobilityKernel( 2, 2 ); -INST_PhaseMobilityKernel( 3, 2 ); -INST_PhaseMobilityKernel( 4, 2 ); -INST_PhaseMobilityKernel( 5, 2 ); - -INST_PhaseMobilityKernel( 1, 3 ); -INST_PhaseMobilityKernel( 2, 3 ); -INST_PhaseMobilityKernel( 3, 3 ); -INST_PhaseMobilityKernel( 4, 3 ); -INST_PhaseMobilityKernel( 5, 3 ); - -#undef INST_PhaseMobilityKernel - - /******************************** FluxKernel ********************************/ -template< localIndex NC, localIndex MAX_NUM_ELEMS, localIndex MAX_STENCIL_SIZE > +template< integer NC, localIndex MAX_NUM_ELEMS, localIndex MAX_STENCIL_SIZE > GEOSX_HOST_DEVICE void FluxKernel:: - compute( localIndex const numPhases, + compute( integer const numPhases, localIndex const stencilSize, localIndex const numFluxElems, arraySlice1d< localIndex const > const seri, @@ -292,14 +69,14 @@ FluxKernel:: arraySlice1d< real64 > const localFlux, arraySlice2d< real64 > const localFluxJacobian ) { - localIndex constexpr NDOF = NC + 1; + integer constexpr NDOF = NC + 1; stackArray1d< real64, NC > compFlux( NC ); stackArray2d< real64, MAX_STENCIL_SIZE * NC > dCompFlux_dP( stencilSize, NC ); stackArray3d< real64, MAX_STENCIL_SIZE * NC * NC > dCompFlux_dC( stencilSize, NC, NC ); // loop over phases, compute and upwind phase flux and sum contributions to each component's flux - for( localIndex ip = 0; ip < numPhases; ++ip ) + for( integer ip = 0; ip < numPhases; ++ip ) { // clear working arrays real64 densMean{}; @@ -343,7 +120,7 @@ FluxKernel:: // average density and derivatives densMean += 0.5 * density; dDensMean_dP[i] = 0.5 * dDens_dP; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dDensMean_dC[i][jc] = 0.5 * dProp_dC[jc]; } @@ -362,7 +139,7 @@ FluxKernel:: real64 capPressure = 0.0; real64 dCapPressure_dP = 0.0; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dCapPressure_dC[ic] = 0.0; } @@ -371,12 +148,12 @@ FluxKernel:: { capPressure = phaseCapPressure[er][esr][ei][0][ip]; - for( localIndex jp = 0; jp < numPhases; ++jp ) + for( integer jp = 0; jp < numPhases; ++jp ) { real64 const dCapPressure_dS = dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp]; dCapPressure_dP += dCapPressure_dS * dPhaseVolFrac_dPres[er][esr][ei][jp]; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dCapPressure_dC[jc] += dCapPressure_dS * dPhaseVolFrac_dComp[er][esr][ei][jp][jc]; } @@ -385,7 +162,7 @@ FluxKernel:: presGrad += transmissibility[i] * (pres[er][esr][ei] + dPres[er][esr][ei] - capPressure); dPresGrad_dP[i] += transmissibility[i] * (1 - dCapPressure_dP) + dTrans_dPres[i] * (pres[er][esr][ei] + dPres[er][esr][ei] - capPressure); - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dPresGrad_dC[i][jc] += -transmissibility[i] * dCapPressure_dC[jc]; } @@ -402,7 +179,7 @@ FluxKernel:: for( localIndex j = 0; j < numFluxElems; ++j ) { dGravHead_dP[j] += dDensMean_dP[j] * gravD + dGravD_dP * densMean; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dGravHead_dC[j][jc] += dDensMean_dC[j][jc] * gravD; } @@ -420,14 +197,14 @@ FluxKernel:: // choose upstream cell localIndex const k_up = (potGrad >= 0) ? 0 : 1; - localIndex er_up = seri[k_up]; - localIndex esr_up = sesri[k_up]; - localIndex ei_up = sei[k_up]; + localIndex const er_up = seri[k_up]; + localIndex const esr_up = sesri[k_up]; + localIndex const ei_up = sei[k_up]; real64 const mobility = phaseMob[er_up][esr_up][ei_up][ip]; // skip the phase flux if phase not present or immobile upstream - if( std::fabs( mobility ) < 1e-20 ) // TODO better constant + if( LvArray::math::abs( mobility ) < 1e-20 ) // TODO better constant { continue; } @@ -436,7 +213,7 @@ FluxKernel:: for( localIndex ke = 0; ke < stencilSize; ++ke ) { dPhaseFlux_dP[ke] += dPresGrad_dP[ke]; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dPhaseFlux_dC[ke][jc] += dPresGrad_dC[ke][jc]; } @@ -458,7 +235,7 @@ FluxKernel:: for( localIndex ke = 0; ke < stencilSize; ++ke ) { dPhaseFlux_dP[ke] *= mobility; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dPhaseFlux_dC[ke][jc] *= mobility; } @@ -469,7 +246,7 @@ FluxKernel:: // add contribution from upstream cell mobility derivatives dPhaseFlux_dP[k_up] += dMob_dP * potGrad; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dPhaseFlux_dC[k_up][jc] += dPhaseMob_dCompSub[jc] * potGrad; } @@ -480,7 +257,7 @@ FluxKernel:: arraySlice2d< real64 const, multifluid::USD_PHASE_COMP_DC-3 > dPhaseCompFrac_dCompSub = dPhaseCompFrac_dComp[er_up][esr_up][ei_up][0][ip]; // compute component fluxes and derivatives using upstream cell composition - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { real64 const ycp = phaseCompFracSub[ic]; compFlux[ic] += phaseFlux * ycp; @@ -489,7 +266,7 @@ FluxKernel:: for( localIndex ke = 0; ke < stencilSize; ++ke ) { dCompFlux_dP[ke][ic] += dPhaseFlux_dP[ke] * ycp; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dCompFlux_dC[ke][ic][jc] += dPhaseFlux_dC[ke][jc] * ycp; } @@ -501,7 +278,7 @@ FluxKernel:: // convert derivatives of component fraction w.r.t. component fractions to derivatives w.r.t. component // densities applyChainRule( NC, dCompFrac_dCompDens[er_up][esr_up][ei_up], dPhaseCompFrac_dCompSub[ic], dProp_dC ); - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dCompFlux_dC[k_up][ic][jc] += phaseFlux * dProp_dC[jc]; } @@ -511,7 +288,7 @@ FluxKernel:: // *** end of upwinding // populate local flux vector and derivatives - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { localFlux[ic] = dt * compFlux[ic]; localFlux[NC + ic] = -dt * compFlux[ic]; @@ -522,7 +299,7 @@ FluxKernel:: localFluxJacobian[ic][localDofIndexPres] = dt * dCompFlux_dP[ke][ic]; localFluxJacobian[NC + ic][localDofIndexPres] = -dt * dCompFlux_dP[ke][ic]; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { localIndex const localDofIndexComp = localDofIndexPres + jc + 1; localFluxJacobian[ic][localDofIndexComp] = dt * dCompFlux_dC[ke][ic][jc]; @@ -532,10 +309,10 @@ FluxKernel:: } } -template< localIndex NC, typename STENCILWRAPPER_TYPE > +template< integer NC, typename STENCILWRAPPER_TYPE > void FluxKernel:: - launch( localIndex const numPhases, + launch( integer const numPhases, STENCILWRAPPER_TYPE const & stencilWrapper, globalIndex const rankOffset, ElementViewConst< arrayView1d< globalIndex const > > const & dofNumber, @@ -632,7 +409,7 @@ FluxKernel:: { globalIndex const offset = dofNumber[seri( iconn, i )][sesri( iconn, i )][sei( iconn, i )]; - for( localIndex jdof = 0; jdof < NDOF; ++jdof ) + for( integer jdof = 0; jdof < NDOF; ++jdof ) { dofColIndices[i * NDOF + jdof] = offset + jdof; } @@ -653,7 +430,7 @@ FluxKernel:: GEOSX_ASSERT_GE( localRow, 0 ); GEOSX_ASSERT_GT( localMatrix.numRows(), localRow + NC ); - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { RAJA::atomicAdd( parallelDeviceAtomic{}, &localRhs[localRow + ic], localFlux[i * NC + ic] ); localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow + ic, @@ -669,7 +446,7 @@ FluxKernel:: #define INST_FluxKernel( NC, STENCILWRAPPER_TYPE ) \ template \ void FluxKernel:: \ - launch< NC, STENCILWRAPPER_TYPE >( localIndex const numPhases, \ + launch< NC, STENCILWRAPPER_TYPE >( integer const numPhases, \ STENCILWRAPPER_TYPE const & stencilWrapper, \ globalIndex const rankOffset, \ ElementViewConst< arrayView1d< globalIndex const > > const & dofNumber, \ @@ -727,11 +504,11 @@ INST_FluxKernel( 5, FaceElementToCellStencilWrapper ); /******************************** CFLFluxKernel ********************************/ -template< localIndex NC, localIndex NUM_ELEMS, localIndex MAX_STENCIL_SIZE > +template< integer NC, localIndex NUM_ELEMS, localIndex MAX_STENCIL_SIZE > GEOSX_HOST_DEVICE void CFLFluxKernel:: - compute( localIndex const numPhases, + compute( integer const numPhases, localIndex const stencilSize, real64 const & dt, arraySlice1d< localIndex const > const seri, @@ -750,7 +527,7 @@ CFLFluxKernel:: ElementView< arrayView2d< real64, compflow::USD_COMP > > const & compOutflux ) { // loop over phases, compute and upwind phase flux and sum contributions to each component's flux - for( localIndex ip = 0; ip < numPhases; ++ip ) + for( integer ip = 0; ip < numPhases; ++ip ) { // clear working arrays @@ -810,7 +587,7 @@ CFLFluxKernel:: RAJA::atomicAdd( parallelDeviceAtomic{}, &phaseOutflux[er_up][esr_up][ei_up][ip], absPhaseFlux ); // increment the component (mass/molar) outflux of the upstream cell - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { real64 const absCompFlux = phaseCompFrac[er_up][esr_up][ei_up][0][ip][ic] * phaseDens[er_up][esr_up][ei_up][0][ip] @@ -820,10 +597,10 @@ CFLFluxKernel:: } } -template< localIndex NC, typename STENCILWRAPPER_TYPE > +template< integer NC, typename STENCILWRAPPER_TYPE > void CFLFluxKernel:: - launch( localIndex const numPhases, + launch( integer const numPhases, real64 const & dt, STENCILWRAPPER_TYPE const & stencilWrapper, ElementViewConst< arrayView1d< real64 const > > const & pres, @@ -883,7 +660,7 @@ CFLFluxKernel:: #define INST_CFLFluxKernel( NC, STENCILWRAPPER_TYPE ) \ template \ void CFLFluxKernel:: \ - launch< NC, STENCILWRAPPER_TYPE >( localIndex const numPhases, \ + launch< NC, STENCILWRAPPER_TYPE >( integer const numPhases, \ real64 const & dt, \ STENCILWRAPPER_TYPE const & stencil, \ ElementViewConst< arrayView1d< real64 const > > const & pres, \ @@ -927,7 +704,7 @@ INST_CFLFluxKernel( 5, FaceElementToCellStencilWrapper ); /******************************** CFLKernel ********************************/ -template< localIndex NP > +template< integer NP > GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE void @@ -984,20 +761,20 @@ CFLKernel:: { // from Keith Coats, IMPES stability: Selection of stable timesteps (2003) real64 totalMob = 0.0; - for( localIndex ip = 0; ip < numMobilePhases; ++ip ) + for( integer ip = 0; ip < numMobilePhases; ++ip ) { totalMob += mob[ip]; } real64 f[2][2]{}; - for( localIndex i = 0; i < 2; ++i ) + for( integer i = 0; i < 2; ++i ) { - for( localIndex j = 0; j < 2; ++j ) + for( integer j = 0; j < 2; ++j ) { f[i][j] = ( i == j )*totalMob - mob[i]; f[i][j] /= (totalMob * mob[j]); real64 sum = 0; - for( localIndex k = 0; k < 3; ++k ) + for( integer k = 0; k < 3; ++k ) { sum += dPhaseRelPerm_dPhaseVolFrac[k][j] / phaseVisc[k] * phaseOutflux[j]; @@ -1012,7 +789,7 @@ CFLKernel:: } -template< localIndex NC > +template< integer NC > GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE void @@ -1026,7 +803,7 @@ CFLKernel:: compCFLNumber = 0.0; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { if( compFrac[ic] > minComponentFraction ) { @@ -1040,7 +817,7 @@ CFLKernel:: } } -template< localIndex NC, localIndex NP > +template< integer NC, integer NP > void CFLKernel:: launch( localIndex const size, @@ -1127,12 +904,12 @@ INST_CFLKernel( 5, 3 ); /******************************** AquiferBCKernel ********************************/ -template< localIndex NC > +template< integer NC > GEOSX_HOST_DEVICE void AquiferBCKernel:: - compute( localIndex const numPhases, - localIndex const ipWater, + compute( integer const numPhases, + integer const ipWater, bool const allowAllPhasesIntoAquifer, real64 const & aquiferVolFlux, real64 const & dAquiferVolFlux_dPres, @@ -1208,11 +985,11 @@ AquiferBCKernel:: } } -template< localIndex NC > +template< integer NC > void AquiferBCKernel:: - launch( localIndex const numPhases, - localIndex const ipWater, + launch( integer const numPhases, + integer const ipWater, bool const allowAllPhasesIntoAquifer, BoundaryStencil const & stencil, globalIndex const rankOffset, @@ -1250,7 +1027,7 @@ AquiferBCKernel:: forAll< parallelDevicePolicy<> >( stencil.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iconn ) { - constexpr localIndex NDOF = NC + 1; + constexpr integer NDOF = NC + 1; // working arrays globalIndex dofColIndices[NDOF]{}; @@ -1296,7 +1073,7 @@ AquiferBCKernel:: // populate dof indices globalIndex const offset = dofNumber[er][esr][ei]; - for( localIndex jdof = 0; jdof < NDOF; ++jdof ) + for( integer jdof = 0; jdof < NDOF; ++jdof ) { dofColIndices[jdof] = offset + jdof; } @@ -1315,7 +1092,7 @@ AquiferBCKernel:: GEOSX_ASSERT_GE( localRow, 0 ); GEOSX_ASSERT_GT( localMatrix.numRows(), localRow + NC ); - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { RAJA::atomicAdd( parallelDeviceAtomic{}, &localRhs[localRow + ic], localFlux[ic] ); localMatrix.addToRow< parallelDeviceAtomic >( localRow + ic, @@ -1330,8 +1107,8 @@ AquiferBCKernel:: #define INST_AquiferBCKernel( NC ) \ template \ void AquiferBCKernel:: \ - launch< NC >( localIndex const numPhases, \ - localIndex const ipWater, \ + launch< NC >( integer const numPhases, \ + integer const ipWater, \ bool const allowAllPhasesIntoAquifer, \ BoundaryStencil const & stencil, \ globalIndex const rankOffset, \ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp index f0c54d476ce..46beab765fd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp @@ -22,12 +22,13 @@ #include "common/DataLayouts.hpp" #include "common/DataTypes.hpp" #include "common/GEOS_RAJA_Interface.hpp" -#include "constitutive/fluid/layouts.hpp" -#include "constitutive/relativePermeability/layouts.hpp" -#include "constitutive/capillaryPressure/layouts.hpp" +#include "constitutive/fluid/MultiFluidBase.hpp" +#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" +#include "constitutive/capillaryPressure/CapillaryPressureBase.hpp" #include "fieldSpecification/AquiferBoundaryCondition.hpp" #include "finiteVolume/BoundaryStencil.hpp" -#include "mesh/ElementRegionManager.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" namespace geosx { @@ -40,68 +41,218 @@ using namespace constitutive; /******************************** PhaseMobilityKernel ********************************/ /** - * @brief Functions to compute phase mobilities and derivatives from density, viscosity and relperm + * @class PhaseMobilityKernel + * @tparam NUM_COMP number of fluid components + * @tparam NUM_PHASE number of fluid phases + * @brief Define the interface for the property kernel in charge of computing the phase mobilities */ -struct PhaseMobilityKernel +template< integer NUM_COMP, integer NUM_PHASE > +class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { - template< localIndex NC, localIndex NP > +public: + + using Base = CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; + using Base::numComp; + + /// Compile time value for the number of phases + static constexpr integer numPhase = NUM_PHASE; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] relperm the relperm model + */ + PhaseMobilityKernel( ObjectManagerBase & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + : Base(), + m_phaseVolFrac( subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFraction >() ), + m_dPhaseVolFrac_dPres( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dPressure >() ), + m_dPhaseVolFrac_dComp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dGlobalCompDensity >() ), + m_dCompFrac_dCompDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity >() ), + m_phaseDens( fluid.phaseDensity() ), + m_dPhaseDens_dPres( fluid.dPhaseDensity_dPressure() ), + m_dPhaseDens_dComp( fluid.dPhaseDensity_dGlobalCompFraction() ), + m_phaseVisc( fluid.phaseViscosity() ), + m_dPhaseVisc_dPres( fluid.dPhaseViscosity_dPressure() ), + m_dPhaseVisc_dComp( fluid.dPhaseViscosity_dGlobalCompFraction() ), + m_phaseRelPerm( relperm.phaseRelPerm() ), + m_dPhaseRelPerm_dPhaseVolFrac( relperm.dPhaseRelPerm_dPhaseVolFraction() ), + m_phaseMob( subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseMobility >() ), + m_dPhaseMob_dPres( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dPressure >() ), + m_dPhaseMob_dComp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dGlobalCompDensity >() ) + {} + + /** + * @brief Compute the phase mobilities in an element + * @tparam FUNC the type of the function that can be used to customize the kernel + * @param[in] ei the element index + * @param[in] phaseMobilityKernelOp the function used to customize the kernel + */ + template< typename FUNC = CompositionalMultiphaseBaseKernels::NoOpFunc > GEOSX_HOST_DEVICE - static void - compute( arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dCompFrac_dCompDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseDens_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseDens_dComp, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseVisc, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseVisc_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseVisc_dComp, - arraySlice1d< real64 const, relperm::USD_RELPERM - 2 > const & phaseRelPerm, - arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > const & dPhaseRelPerm_dPhaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & phaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & dPhaseVolFrac_dPres, - arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const & dPhaseVolFrac_dComp, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & phaseMob, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & dPhaseMob_dPres, - arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const & dPhaseMob_dComp ); - - template< localIndex NC, localIndex NP > - static void - launch( localIndex const size, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ); + void compute( localIndex const ei, + FUNC && phaseMobilityKernelOp = CompositionalMultiphaseBaseKernels::NoOpFunc{} ) const + { + arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const phaseDens = m_phaseDens[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const dPhaseDens_dPres = m_dPhaseDens_dPres[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseDens_dComp = m_dPhaseDens_dComp[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const phaseVisc = m_phaseVisc[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const dPhaseVisc_dPres = m_dPhaseVisc_dPres[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseVisc_dComp = m_dPhaseVisc_dComp[ei][0]; + arraySlice1d< real64 const, relperm::USD_RELPERM - 2 > const phaseRelPerm = m_phaseRelPerm[ei][0]; + arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > const dPhaseRelPerm_dPhaseVolFrac = m_dPhaseRelPerm_dPhaseVolFrac[ei][0]; + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const phaseVolFrac = m_phaseVolFrac[ei]; + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const dPhaseVolFrac_dPres = m_dPhaseVolFrac_dPres[ei]; + arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const dPhaseVolFrac_dComp = m_dPhaseVolFrac_dComp[ei]; + arraySlice1d< real64, compflow::USD_PHASE - 1 > const phaseMob = m_phaseMob[ei]; + arraySlice1d< real64, compflow::USD_PHASE - 1 > const dPhaseMob_dPres = m_dPhaseMob_dPres[ei]; + arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const dPhaseMob_dComp = m_dPhaseMob_dComp[ei]; + + real64 dRelPerm_dC[numComp]{}; + real64 dDens_dC[numComp]{}; + real64 dVisc_dC[numComp]{}; + + for( integer ip = 0; ip < numPhase; ++ip ) + { + + // compute the phase mobility only if the phase is present + bool const phaseExists = (phaseVolFrac[ip] > 0); + if( !phaseExists ) + { + phaseMob[ip] = 0.0; + dPhaseMob_dPres[ip] = 0.0; + for( integer jc = 0; jc < numComp; ++jc ) + { + dPhaseMob_dComp[ip][jc] = 0.0; + } + continue; + } + + real64 const density = phaseDens[ip]; + real64 const dDens_dP = dPhaseDens_dPres[ip]; + applyChainRule( numComp, dCompFrac_dCompDens, dPhaseDens_dComp[ip], dDens_dC ); + + real64 const viscosity = phaseVisc[ip]; + real64 const dVisc_dP = dPhaseVisc_dPres[ip]; + applyChainRule( numComp, dCompFrac_dCompDens, dPhaseVisc_dComp[ip], dVisc_dC ); + + real64 const relPerm = phaseRelPerm[ip]; + real64 dRelPerm_dP = 0.0; + for( integer ic = 0; ic < numComp; ++ic ) + { + dRelPerm_dC[ic] = 0.0; + } + + for( integer jp = 0; jp < numPhase; ++jp ) + { + real64 const dRelPerm_dS = dPhaseRelPerm_dPhaseVolFrac[ip][jp]; + dRelPerm_dP += dRelPerm_dS * dPhaseVolFrac_dPres[jp]; + + for( integer jc = 0; jc < numComp; ++jc ) + { + dRelPerm_dC[jc] += dRelPerm_dS * dPhaseVolFrac_dComp[jp][jc]; + } + } + + real64 const mobility = relPerm * density / viscosity; + + phaseMob[ip] = mobility; + dPhaseMob_dPres[ip] = dRelPerm_dP * density / viscosity + + mobility * (dDens_dP / density - dVisc_dP / viscosity); + + // compositional derivatives + for( integer jc = 0; jc < numComp; ++jc ) + { + dPhaseMob_dComp[ip][jc] = dRelPerm_dC[jc] * density / viscosity + + mobility * (dDens_dC[jc] / density - dVisc_dC[jc] / viscosity); + } + + // call the lambda in the phase loop to allow the reuse of the relperm, density, viscosity, and mobility + // possible use: assemble the derivatives wrt temperature + phaseMobilityKernelOp( ip, phaseMob[ip], dPhaseMob_dPres[ip], dPhaseMob_dComp[ip] ); + } + } + +protected: + + // inputs + + /// Views on the phase volume fractions + arrayView2d< real64 const, compflow::USD_PHASE > m_phaseVolFrac; + arrayView2d< real64 const, compflow::USD_PHASE > m_dPhaseVolFrac_dPres; + arrayView3d< real64 const, compflow::USD_PHASE_DC > m_dPhaseVolFrac_dComp; + arrayView3d< real64 const, compflow::USD_COMP_DC > m_dCompFrac_dCompDens; + + /// Views on the phase densities + arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseDens; + arrayView3d< real64 const, multifluid::USD_PHASE > m_dPhaseDens_dPres; + arrayView4d< real64 const, multifluid::USD_PHASE_DC > m_dPhaseDens_dComp; + + /// Views on the phase viscosities + arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseVisc; + arrayView3d< real64 const, multifluid::USD_PHASE > m_dPhaseVisc_dPres; + arrayView4d< real64 const, multifluid::USD_PHASE_DC > m_dPhaseVisc_dComp; + + /// Views on the phase relative permeabilities + arrayView3d< real64 const, relperm::USD_RELPERM > m_phaseRelPerm; + arrayView4d< real64 const, relperm::USD_RELPERM_DS > m_dPhaseRelPerm_dPhaseVolFrac; + + // outputs + + /// Views on the phase mobilities + arrayView2d< real64, compflow::USD_PHASE > m_phaseMob; + arrayView2d< real64, compflow::USD_PHASE > m_dPhaseMob_dPres; + arrayView3d< real64, compflow::USD_PHASE_DC > m_dPhaseMob_dComp; - template< localIndex NC, localIndex NP > - static void - launch( SortedArrayView< localIndex const > const & targetSet, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ); }; +/** + * @class PhaseMobilityKernelFactory + */ +class PhaseMobilityKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numComp the number of fluid components + * @param[in] numPhase the number of fluid phases + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] relperm the relperm model + */ + template< typename POLICY > + static void + createAndLaunch( integer const numComp, + integer const numPhase, + ObjectManagerBase & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + { + if( numPhase == 2 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + } +}; /******************************** FluxKernel ********************************/ @@ -120,10 +271,10 @@ struct FluxKernel template< typename VIEWTYPE > using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; - template< localIndex NC, localIndex MAX_NUM_ELEMS, localIndex MAX_STENCIL_SIZE > + template< integer NC, localIndex MAX_NUM_ELEMS, localIndex MAX_STENCIL_SIZE > GEOSX_HOST_DEVICE static void - compute( localIndex const numPhases, + compute( integer const numPhases, localIndex const stencilSize, localIndex const numFluxElems, arraySlice1d< localIndex const > const seri, @@ -153,9 +304,9 @@ struct FluxKernel arraySlice1d< real64 > const localFlux, arraySlice2d< real64 > const localFluxJacobian ); - template< localIndex NC, typename STENCILWRAPPER_TYPE > + template< integer NC, typename STENCILWRAPPER_TYPE > static void - launch( localIndex const numPhases, + launch( integer const numPhases, STENCILWRAPPER_TYPE const & stencilWrapper, globalIndex const rankOffset, ElementViewConst< arrayView1d< globalIndex const > > const & dofNumber, @@ -205,10 +356,10 @@ struct CFLFluxKernel template< typename VIEWTYPE > using ElementView = ElementRegionManager::ElementView< VIEWTYPE >; - template< localIndex NC, localIndex NUM_ELEMS, localIndex MAX_STENCIL_SIZE > + template< integer NC, localIndex NUM_ELEMS, localIndex MAX_STENCIL_SIZE > GEOSX_HOST_DEVICE static void - compute( localIndex const numPhases, + compute( integer const numPhases, localIndex const stencilSize, real64 const & dt, arraySlice1d< localIndex const > const seri, @@ -226,9 +377,9 @@ struct CFLFluxKernel ElementView< arrayView2d< real64, compflow::USD_PHASE > > const & phaseOutflux, ElementView< arrayView2d< real64, compflow::USD_COMP > > const & compOutflux ); - template< localIndex NC, typename STENCILWRAPPER_TYPE > + template< integer NC, typename STENCILWRAPPER_TYPE > static void - launch( localIndex const numPhases, + launch( integer const numPhases, real64 const & dt, STENCILWRAPPER_TYPE const & stencil, ElementViewConst< arrayView1d< real64 const > > const & pres, @@ -256,7 +407,7 @@ struct CFLKernel static constexpr real64 minPhaseMobility = 1e-12; static constexpr real64 minComponentFraction = 1e-12; - template< localIndex NP > + template< integer NP > GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE static void @@ -268,7 +419,7 @@ struct CFLKernel arraySlice1d< real64 const, compflow::USD_PHASE- 1 > phaseOutflux, real64 & phaseCFLNumber ); - template< localIndex NC > + template< integer NC > GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE static void @@ -278,7 +429,7 @@ struct CFLKernel arraySlice1d< real64 const, compflow::USD_COMP - 1 > compOutflux, real64 & compCFLNumber ); - template< localIndex NC, localIndex NP > + template< integer NC, integer NP > static void launch( localIndex const size, arrayView1d< real64 const > const & volume, @@ -315,11 +466,11 @@ struct AquiferBCKernel template< typename VIEWTYPE > using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; - template< localIndex NC > + template< integer NC > GEOSX_HOST_DEVICE static void - compute( localIndex const numPhases, - localIndex const ipWater, + compute( integer const numPhases, + integer const ipWater, bool const allowAllPhasesIntoAquifer, real64 const & aquiferVolFlux, real64 const & dAquiferVolFlux_dPres, @@ -339,10 +490,10 @@ struct AquiferBCKernel real64 ( &localFlux )[NC], real64 ( &localFluxJacobian )[NC][NC+1] ); - template< localIndex NC > + template< integer NC > static void - launch( localIndex const numPhases, - localIndex const ipWater, + launch( integer const numPhases, + integer const ipWater, bool const allowAllPhasesIntoAquifer, BoundaryStencil const & stencil, globalIndex const rankOffset, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index 02bb3a858fb..48c717b0a3e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -892,58 +892,16 @@ void CompositionalMultiphaseHybridFVM::updatePhaseMobility( ObjectManagerBase & { GEOSX_MARK_FUNCTION; - // note that the phase mobility computed here does NOT include phase density - - // outputs - - arrayView2d< real64, compflow::USD_PHASE > const phaseMob = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::phaseMobility >(); - - arrayView2d< real64, compflow::USD_PHASE > const dPhaseMob_dPres = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dPressure >(); - - arrayView3d< real64, compflow::USD_PHASE_DC > const dPhaseMob_dComp = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dGlobalCompDensity >(); - - // inputs - - arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFraction >(); - - arrayView2d< real64 const, compflow::USD_PHASE > const dPhaseVolFrac_dPres = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dPressure >(); - - arrayView3d< real64 const, compflow::USD_PHASE_DC > const dPhaseVolFrac_dComp = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dGlobalCompDensity >(); - - arrayView3d< real64 const, compflow::USD_COMP_DC > const dCompFrac_dCompDens = - dataGroup.getExtrinsicData< extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity >(); - MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, m_fluidModelNames[targetIndex] ); - - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc = fluid.phaseViscosity(); - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres = fluid.dPhaseViscosity_dPressure(); - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp = fluid.dPhaseViscosity_dGlobalCompFraction(); - RelativePermeabilityBase const & relperm = getConstitutiveModel< RelativePermeabilityBase >( dataGroup, m_relPermModelNames[targetIndex] ); - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm = relperm.phaseRelPerm(); - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac = relperm.dPhaseRelPerm_dPhaseVolFraction(); - - KernelLaunchSelector2< PhaseMobilityKernel >( m_numComponents, m_numPhases, - dataGroup.size(), - dCompFrac_dCompDens, - phaseVisc, - dPhaseVisc_dPres, - dPhaseVisc_dComp, - phaseRelPerm, - dPhaseRelPerm_dPhaseVolFrac, - phaseVolFrac, - dPhaseVolFrac_dPres, - dPhaseVolFrac_dComp, - phaseMob, - dPhaseMob_dPres, - dPhaseMob_dComp ); + CompositionalMultiphaseHybridFVMKernels:: + PhaseMobilityKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid, + relperm ); } REGISTER_CATALOG_ENTRY( SolverBase, CompositionalMultiphaseHybridFVM, std::string const &, Group * const ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp index 244be1c3922..cd7222d1431 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp @@ -34,7 +34,7 @@ namespace CompositionalMultiphaseHybridFVMKernels /******************************** UpwindingHelper ********************************/ -template< localIndex NC, localIndex NP > +template< integer NC, integer NP > void UpwindingHelper:: upwindViscousCoefficient( localIndex const (&localIds)[ 3 ], @@ -70,24 +70,24 @@ UpwindingHelper:: real64 totalMob = 0; real64 dTotalMob_dPres = 0; real64 dTotalMob_dCompDens[ NC ]{}; - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { totalMob = totalMob + phaseMob[er][esr][ei][ip]; dTotalMob_dPres = dTotalMob_dPres + dPhaseMob_dPres[er][esr][ei][ip]; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dTotalMob_dCompDens[ic] = dTotalMob_dCompDens[ic] + dPhaseMob_dCompDens[er][esr][ei][ip][ic]; } } real64 const totalMobInv = 1.0 / totalMob; - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { // 3) Compute viscous mobility ratio: \frac{\lambda_{\ell}}{\lambda_T} real64 const upwMobRatio = phaseMob[er][esr][ei][ip] * totalMobInv; real64 const dUpwMobRatio_dPres = ( dPhaseMob_dPres[er][esr][ei][ip] - upwMobRatio * dTotalMob_dPres ) * totalMobInv; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dUpwMobRatio_dCompDens[ic] = ( dPhaseMob_dCompDens[er][esr][ei][ip][ic] - upwMobRatio * dTotalMob_dCompDens[ic] ) * totalMobInv; @@ -101,14 +101,14 @@ UpwindingHelper:: real64 const upwDensMobRatio = phaseDens[er][esr][ei][0][ip] * upwMobRatio; real64 const dUpwDensMobRatio_dPres = dPhaseDens_dPres[er][esr][ei][0][ip] * upwMobRatio + phaseDens[er][esr][ei][0][ip] * dUpwMobRatio_dPres; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dUpwDensMobRatio_dCompDens[ic] = dPhaseDens_dC[ic] * upwMobRatio + phaseDens[er][esr][ei][0][ip] * dUpwMobRatio_dCompDens[ic]; } // 5) Multiply density mobility ratio by phase comp fraction: x_{c,\ell} \rho^{up}_{\ell} \frac{\lambda_{\ell}}{\lambda_T} - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { applyChainRule( NC, dCompFrac_dCompDens[er][esr][ei], @@ -117,7 +117,7 @@ UpwindingHelper:: upwPhaseViscCoef[ip][ic] = phaseCompFrac[er][esr][ei][0][ip][ic] * upwDensMobRatio; dUpwPhaseViscCoef_dPres[ip][ic] = dPhaseCompFrac_dPres[er][esr][ei][0][ip][ic] * upwDensMobRatio + phaseCompFrac[er][esr][ei][0][ip][ic] * dUpwDensMobRatio_dPres; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dUpwPhaseViscCoef_dCompDens[ip][ic][jc] = dPhaseCompFrac_dC[jc] * upwDensMobRatio + phaseCompFrac[er][esr][ei][0][ip][ic] * dUpwDensMobRatio_dCompDens[jc]; @@ -129,7 +129,7 @@ UpwindingHelper:: } -template< localIndex NC, localIndex NP > +template< integer NC, integer NP > GEOSX_HOST_DEVICE void UpwindingHelper:: @@ -193,10 +193,10 @@ UpwindingHelper:: real64 dPhaseDens_dC[ NC ]{}; real64 dPhaseCompFrac_dC[ NC ]{}; - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { localIndex k = 0; - for( localIndex jp = 0; jp < NP; ++jp ) + for( integer jp = 0; jp < NP; ++jp ) { if( ip == jp ) { @@ -219,7 +219,7 @@ UpwindingHelper:: dMobRatio_dPres[posd] = ( dPhaseMob_dPres[erd][esrd][eid][jp] * phaseMob[eru][esru][eiu][ip] - mobRatio * dTotalMob_dPres[posd] ) * totalMobInv; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dMobRatio_dCompDens[posu][ic] = ( dPhaseMob_dCompDens[eru][esru][eiu][ip][ic] * phaseMob[erd][esrd][eid][jp] - mobRatio * dTotalMob_dCompDens[posu][ic] ) * totalMobInv; @@ -236,7 +236,7 @@ UpwindingHelper:: dDensMobRatio_dPres[posu] = dPhaseDens_dPres[eru][esru][eiu][0][ip] * mobRatio + phaseDens[eru][esru][eiu][0][ip] * dMobRatio_dPres[posu]; dDensMobRatio_dPres[posd] = phaseDens[eru][esru][eiu][0][ip] * dMobRatio_dPres[posd]; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dDensMobRatio_dCompDens[posu][ic] = dPhaseDens_dC[ic] * mobRatio + phaseDens[eru][esru][eiu][0][ip] * dMobRatio_dCompDens[posu][ic]; @@ -244,7 +244,7 @@ UpwindingHelper:: } // 3.d) Compute the final gravity coefficient \x_{up}_{c,p} \rho_p \frac{\lambda_l \lambda_m}{\lambda_T} - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { applyChainRule( NC, dCompFrac_dCompDens[eru][esru][eiu], @@ -255,7 +255,7 @@ UpwindingHelper:: + phaseCompFrac[eru][esru][eiu][0][ip][ic] * dDensMobRatio_dPres[posu]; dUpwPhaseGravCoef_dPres[ip][k][ic][posd] = phaseCompFrac[eru][esru][eiu][0][ip][ic] * dDensMobRatio_dPres[posd]; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dUpwPhaseGravCoef_dCompDens[ip][k][ic][posu][jc] = dPhaseCompFrac_dC[jc] * densMobRatio + phaseCompFrac[eru][esru][eiu][0][ip][ic] * dDensMobRatio_dCompDens[posu][jc]; @@ -267,7 +267,7 @@ UpwindingHelper:: } } -template< localIndex NC, localIndex NP > +template< integer NC, integer NP > GEOSX_HOST_DEVICE void UpwindingHelper:: @@ -293,7 +293,7 @@ UpwindingHelper:: real64 dPhaseMassDens_dCNeighbor[ NC ]{}; real64 dPhaseMassDens_dC[ NC ]{}; - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { applyChainRule( NC, dCompFrac_dCompDens[er][esr][ei], @@ -305,7 +305,7 @@ UpwindingHelper:: dPhaseMassDens_dCNeighbor ); localIndex k = 0; - for( localIndex jp = 0; jp < NP; ++jp ) + for( integer jp = 0; jp < NP; ++jp ) { if( ip == jp ) { @@ -322,7 +322,7 @@ UpwindingHelper:: dPhaseGravTerm_dPres[ip][k][Pos::NEIGHBOR] = ( -dPhaseMassDens_dPres[ern][esrn][ein][0][ip] + dPhaseMassDens_dPres[ern][esrn][ein][0][jp] ); dPhaseGravTerm_dPres[ip][k][Pos::NEIGHBOR] *= 0.5 * transGravCoef; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dPhaseGravTerm_dCompDens[ip][k][Pos::LOCAL][ic] = -0.5 * transGravCoef * dPhaseMassDens_dCLoc[ic]; dPhaseGravTerm_dCompDens[ip][k][Pos::NEIGHBOR][ic] = -0.5 * transGravCoef * dPhaseMassDens_dCNeighbor[ic]; @@ -331,7 +331,7 @@ UpwindingHelper:: dCompFrac_dCompDens[er][esr][ei], dPhaseMassDens_dCompFrac[er][esr][ei][0][jp], dPhaseMassDens_dC ); - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dPhaseGravTerm_dCompDens[ip][k][Pos::LOCAL][ic] += 0.5 * transGravCoef * dPhaseMassDens_dC[ic]; } @@ -348,7 +348,7 @@ UpwindingHelper:: } } -template< localIndex NC, localIndex NP > +template< integer NC, integer NP > GEOSX_HOST_DEVICE void UpwindingHelper:: @@ -369,7 +369,7 @@ UpwindingHelper:: phaseGravTerm, totalMobIds, totalMobPos ); - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { localIndex const er = totalMobIds[ip][0]; localIndex const esr = totalMobIds[ip][1]; @@ -377,7 +377,7 @@ UpwindingHelper:: localIndex const pos = totalMobPos[ip]; totalMob = totalMob + phaseMob[er][esr][ei][ip]; dTotalMob_dPres[pos] = dTotalMob_dPres[pos] + dPhaseMob_dPres[er][esr][ei][pos]; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dTotalMob_dCompDens[pos][ic] = dTotalMob_dCompDens[pos][ic] + dPhaseMob_dCompDens[er][esr][ei][ip][ic]; } @@ -421,7 +421,7 @@ UpwindingHelper:: } } -template< localIndex NP > +template< integer NP > GEOSX_HOST_DEVICE void UpwindingHelper:: @@ -460,7 +460,7 @@ UpwindingHelper:: { // TODO Francois: this should be improved // currently this implements the algorithm proposed by SH Lee - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { if( ( gravTerm[ip][0] >= 0 && gravTerm[ip][1] >= 0 ) || // includes the no-buoyancy case ( ( fabs( gravTerm[ip][0] ) >= fabs( gravTerm[ip][1] ) ) && gravTerm[ip][1] >= 0 ) || @@ -587,7 +587,7 @@ INST_UpwindingHelperNP( 3 ); /******************************** AssemblerKernelHelper ********************************/ -template< localIndex NF, localIndex NC, localIndex NP > +template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE void AssemblerKernelHelper:: @@ -617,7 +617,7 @@ AssemblerKernelHelper:: real64 dPhaseMobPotDif_dCompDens[ NC ]{}; // 0) precompute dPhaseDens_dC since it is always computed at the element center - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { applyChainRule( NC, dElemCompFrac_dCompDens, @@ -625,11 +625,11 @@ AssemblerKernelHelper:: dPhaseMassDens_dC[ip] ); } - for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) { // now in the following nested loop, // we compute the contribution of face jfaceLoc to the one sided total volumetric flux at face iface - for( localIndex jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) { // depth difference between element center and face center @@ -637,7 +637,7 @@ AssemblerKernelHelper:: real64 const fGravCoef = faceGravCoef[elemToFaces[jfaceLoc]]; real64 const gravCoefDif = ccGravCoef - fGravCoef; - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { // 1) compute the potential diff between the cell center and the face center @@ -648,7 +648,7 @@ AssemblerKernelHelper:: real64 const presDif = ccPres - fPres; real64 const dPresDif_dPres = 1; real64 const dPresDif_dFacePres = -1; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dPresDif_dCompDens[ic] = 0.0; // no capillary pressure } @@ -668,7 +668,7 @@ AssemblerKernelHelper:: real64 const dPhaseMobPotDif_dPres = dElemPhaseMob_dPres[ip] * phasePotDif + elemPhaseMob[ip] * (dPresDif_dPres - dPhaseGravDif_dPres); real64 const dPhaseMobPotDif_dFacePres = elemPhaseMob[ip] * dPresDif_dFacePres; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dPhaseMobPotDif_dCompDens[ic] = dElemPhaseMob_dCompDens[ip][ic] * phasePotDif + elemPhaseMob[ip] * (dPresDif_dCompDens[ic] - dPhaseGravDif_dCompDens[ic]); @@ -691,7 +691,7 @@ AssemblerKernelHelper:: } } -template< localIndex NF, localIndex NC, localIndex NP > +template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE void AssemblerKernelHelper:: @@ -729,12 +729,12 @@ AssemblerKernelHelper:: arrayView1d< real64 > const & localRhs ) { using namespace CompositionalMultiphaseUtilities; - localIndex constexpr NDOF = NC+1; + integer constexpr NDOF = NC+1; // dof numbers globalIndex dofColIndicesElemVars[ NDOF*(NF+1) ]{}; globalIndex dofColIndicesFaceVars[ NF ]{}; - for( localIndex idof = 0; idof < NDOF; ++idof ) + for( integer idof = 0; idof < NDOF; ++idof ) { dofColIndicesElemVars[idof] = elemDofNumber[localIds[0]][localIds[1]][localIds[2]] + idof; } @@ -763,7 +763,7 @@ AssemblerKernelHelper:: real64 dUpwPhaseGravCoef_dCompDens[ NP ][ NP-1 ][ NC ][ 2 ][ NC ]{}; // for each element, loop over the one-sided faces - for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) { // 1) Find if there is a neighbor, and if there is, grab the indices of the neighbor element @@ -874,7 +874,7 @@ AssemblerKernelHelper:: // we are ready to assemble the local flux and its derivatives // no need for atomic adds - each row is assembled by a single thread - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { localIndex const eqnRowLocalIndex = LvArray::integerConversion< localIndex >( elemDofNumber[localIds[0]][localIds[1]][localIds[2]] + ic - rankOffset ); @@ -899,7 +899,7 @@ AssemblerKernelHelper:: } } -template< localIndex NF, localIndex NC, localIndex NP > +template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE void AssemblerKernelHelper:: @@ -922,12 +922,12 @@ AssemblerKernelHelper:: globalIndex ( & dofColIndicesElemVars )[ (NC+1)*(NF+1) ], globalIndex ( & dofColIndicesFaceVars )[ NF ] ) { - localIndex constexpr NDOF = NC+1; + integer constexpr NDOF = NC+1; localIndex const elemVarsOffset = NDOF*(ifaceLoc+1); - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { // compute the mass flux at the one-sided face plus its derivatives // add the newly computed flux to the sum @@ -943,7 +943,7 @@ AssemblerKernelHelper:: dDivMassFluxes_dElemVars[ic][0] = dDivMassFluxes_dElemVars[ic][0] + ( elemDofNumber == upwViscDofNumber ) * dt * dUpwPhaseViscCoef_dPres[ip][ic] * oneSidedVolFlux[ifaceLoc]; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dDivMassFluxes_dElemVars[ic][jc+1] = dDivMassFluxes_dElemVars[ic][jc+1] + dt_upwPhaseViscCoef * dOneSidedVolFlux_dCompDens[ifaceLoc][jc]; @@ -957,14 +957,14 @@ AssemblerKernelHelper:: + ( elemDofNumber != upwViscDofNumber ) * dt * dUpwPhaseViscCoef_dPres[ip][ic] * oneSidedVolFlux[ifaceLoc]; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1] = dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1] + ( elemDofNumber != upwViscDofNumber ) * dt * dUpwPhaseViscCoef_dCompDens[ip][ic][jc] * oneSidedVolFlux[ifaceLoc]; } - for( localIndex jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) { dDivMassFluxes_dFaceVars[ic][jfaceLoc] = dDivMassFluxes_dFaceVars[ic][jfaceLoc] + dt_upwPhaseViscCoef * dOneSidedVolFlux_dFacePres[ifaceLoc][jfaceLoc]; @@ -973,14 +973,14 @@ AssemblerKernelHelper:: } // collect the relevant dof numbers - for( localIndex idof = 0; idof < NDOF; ++idof ) + for( integer idof = 0; idof < NDOF; ++idof ) { dofColIndicesElemVars[elemVarsOffset+idof] = neighborDofNumber + idof; } dofColIndicesFaceVars[ifaceLoc] = faceDofNumber; } -template< localIndex NF, localIndex NC, localIndex NP > +template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE void AssemblerKernelHelper:: @@ -995,14 +995,14 @@ AssemblerKernelHelper:: real64 ( & divMassFluxes )[ NC ], real64 ( & dDivMassFluxes_dElemVars )[ NC ][ (NC+1)*(NF+1) ] ) { - localIndex constexpr NDOF = NC+1; + integer constexpr NDOF = NC+1; localIndex const elemVarsOffset = NDOF*(ifaceLoc+1); - for( localIndex ip = 0; ip < NP; ++ip ) + for( integer ip = 0; ip < NP; ++ip ) { - for( localIndex jp = 0; jp < NP - 1; ++jp ) + for( integer jp = 0; jp < NP - 1; ++jp ) { - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { real64 const dt_upwPhaseGravCoef = dt * upwPhaseGravCoef[ip][jp][ic]; @@ -1015,7 +1015,7 @@ AssemblerKernelHelper:: dDivMassFluxes_dElemVars[ic][0] = dDivMassFluxes_dElemVars[ic][0] + dt_upwPhaseGravCoef * dPhaseGravTerm_dPres[ip][jp][Pos::LOCAL]; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dDivMassFluxes_dElemVars[ic][jc+1] = dDivMassFluxes_dElemVars[ic][jc+1] + dt * dUpwPhaseGravCoef_dCompDens[ip][jp][ic][Pos::LOCAL][jc] * phaseGravTerm[ip][jp]; @@ -1029,7 +1029,7 @@ AssemblerKernelHelper:: dDivMassFluxes_dElemVars[ic][elemVarsOffset] = dDivMassFluxes_dElemVars[ic][elemVarsOffset] + dt_upwPhaseGravCoef * dPhaseGravTerm_dPres[ip][jp][Pos::NEIGHBOR]; - for( localIndex jc = 0; jc < NC; ++jc ) + for( integer jc = 0; jc < NC; ++jc ) { dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1] = dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1] + dt * dUpwPhaseGravCoef_dCompDens[ip][jp][ic][Pos::NEIGHBOR][jc] * phaseGravTerm[ip][jp]; @@ -1041,7 +1041,7 @@ AssemblerKernelHelper:: } } -template< localIndex NF, localIndex NC, localIndex NP > +template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE void AssemblerKernelHelper:: @@ -1057,7 +1057,7 @@ AssemblerKernelHelper:: CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - localIndex constexpr NDOF = NC+1; + integer constexpr NDOF = NC+1; // fluxes real64 dFlux_dElemVars[ NDOF ]{}; @@ -1066,13 +1066,13 @@ AssemblerKernelHelper:: // dof numbers globalIndex dofColIndicesElemVars[ NDOF ]{}; globalIndex dofColIndicesFaceVars[ NF ]{}; - for( localIndex idof = 0; idof < NDOF; ++idof ) + for( integer idof = 0; idof < NDOF; ++idof ) { dofColIndicesElemVars[idof] = elemDofNumber + idof; } // for each element, loop over the local (one-sided) faces - for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) { if( faceGhostRank[elemToFaces[ifaceLoc]] >= 0 ) { @@ -1082,12 +1082,12 @@ AssemblerKernelHelper:: // flux at this face real64 const flux = oneSidedVolFlux[ifaceLoc]; dFlux_dElemVars[0] = dOneSidedVolFlux_dPres[ifaceLoc]; - for( localIndex ic = 0; ic < NC; ++ic ) + for( integer ic = 0; ic < NC; ++ic ) { dFlux_dElemVars[ic+1] = dOneSidedVolFlux_dCompDens[ifaceLoc][ic]; } - for( localIndex jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) { dFlux_dFaceVars[jfaceLoc] = dOneSidedVolFlux_dFacePres[ifaceLoc][jfaceLoc]; dofColIndicesFaceVars[jfaceLoc] = faceDofNumber[elemToFaces[jfaceLoc]]; @@ -1264,7 +1264,7 @@ INST_AssemblerKernelHelper( 6, 5, 3 ); /******************************** AssemblerKernel ********************************/ -template< localIndex NF, localIndex NC, localIndex NP > +template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE void AssemblerKernel:: @@ -1487,7 +1487,7 @@ INST_AssemblerKernel( 6, 5, 3 ); /******************************** FluxKernel ********************************/ -template< localIndex NF, localIndex NC, localIndex NP, typename IP_TYPE > +template< integer NF, integer NC, integer NP, typename IP_TYPE > void FluxKernel:: launch( localIndex er, localIndex esr, @@ -1745,203 +1745,6 @@ INST_FluxKernel( 6, 5, 3, mimeticInnerProduct::BdVLMInnerProduct const ); #undef INST_FluxKernel - -/******************************** PhaseMobilityKernel ********************************/ - -template< localIndex NC, localIndex NP > -GEOSX_HOST_DEVICE -void -PhaseMobilityKernel:: - compute( arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dCompFrac_dCompDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseVisc, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseVisc_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseVisc_dComp, - arraySlice1d< real64 const, relperm::USD_RELPERM - 2 > const & phaseRelPerm, - arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > const & dPhaseRelPerm_dPhaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & phaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & dPhaseVolFrac_dPres, - arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const & dPhaseVolFrac_dComp, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & phaseMob, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & dPhaseMob_dPres, - arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const & dPhaseMob_dComp ) -{ - real64 dRelPerm_dC[NC]; - real64 dVisc_dC[NC]; - - for( localIndex ip = 0; ip < NP; ++ip ) - { - // compute the phase mobility only if the phase is present - bool const phaseExists = (phaseVolFrac[ip] > 0); - if( !phaseExists ) - { - phaseMob[ip] = 0.; - dPhaseMob_dPres[ip] = 0.; - for( localIndex jc = 0; jc < NC; ++jc ) - { - dPhaseMob_dComp[ip][jc] = 0.; - } - continue; - } - - real64 const viscosity = phaseVisc[ip]; - real64 const dVisc_dP = dPhaseVisc_dPres[ip]; - applyChainRule( NC, dCompFrac_dCompDens, dPhaseVisc_dComp[ip], dVisc_dC ); - - real64 const relPerm = phaseRelPerm[ip]; - real64 dRelPerm_dP = 0.0; - for( localIndex ic = 0; ic < NC; ++ic ) - { - dRelPerm_dC[ic] = 0.0; - } - - for( localIndex jp = 0; jp < NP; ++jp ) - { - real64 const dRelPerm_dS = dPhaseRelPerm_dPhaseVolFrac[ip][jp]; - dRelPerm_dP += dRelPerm_dS * dPhaseVolFrac_dPres[jp]; - - for( localIndex jc = 0; jc < NC; ++jc ) - { - dRelPerm_dC[jc] += dRelPerm_dS * dPhaseVolFrac_dComp[jp][jc]; - } - } - - real64 const mobility = relPerm / viscosity; - - phaseMob[ip] = mobility; - dPhaseMob_dPres[ip] = dRelPerm_dP / viscosity - - mobility * dVisc_dP / viscosity; - - // compositional derivatives - for( localIndex jc = 0; jc < NC; ++jc ) - { - dPhaseMob_dComp[ip][jc] = dRelPerm_dC[jc] / viscosity - - mobility * dVisc_dC[jc] / viscosity; - } - } -} - -template< localIndex NC, localIndex NP > -void PhaseMobilityKernel:: - launch( localIndex const size, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ) -{ - forAll< parallelDevicePolicy<> >( size, [=] GEOSX_HOST_DEVICE ( localIndex const a ) - { - compute< NC, NP >( dCompFrac_dCompDens[a], - phaseVisc[a][0], - dPhaseVisc_dPres[a][0], - dPhaseVisc_dComp[a][0], - phaseRelPerm[a][0], - dPhaseRelPerm_dPhaseVolFrac[a][0], - phaseVolFrac[a], - dPhaseVolFrac_dPres[a], - dPhaseVolFrac_dComp[a], - phaseMob[a], - dPhaseMob_dPres[a], - dPhaseMob_dComp[a] ); - } ); -} - -template< localIndex NC, localIndex NP > -void PhaseMobilityKernel:: - launch( SortedArrayView< localIndex const > const & targetSet, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ) -{ - forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const i ) - { - localIndex const a = targetSet[ i ]; - compute< NC, NP >( dCompFrac_dCompDens[a], - phaseVisc[a][0], - dPhaseVisc_dPres[a][0], - dPhaseVisc_dComp[a][0], - phaseRelPerm[a][0], - dPhaseRelPerm_dPhaseVolFrac[a][0], - phaseVolFrac[a], - dPhaseVolFrac_dPres[a], - dPhaseVolFrac_dComp[a], - phaseMob[a], - dPhaseMob_dPres[a], - dPhaseMob_dComp[a] ); - } ); -} - -#define INST_PhaseMobilityKernel( NC, NP ) \ - template \ - void \ - PhaseMobilityKernel:: \ - launch< NC, NP >( localIndex const size, \ - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, \ - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, \ - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, \ - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, \ - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, \ - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, \ - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, \ - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ); \ - template \ - void \ - PhaseMobilityKernel:: \ - launch< NC, NP >( SortedArrayView< localIndex const > const & targetSet, \ - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, \ - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, \ - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, \ - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, \ - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, \ - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, \ - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, \ - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ) - -INST_PhaseMobilityKernel( 1, 1 ); -INST_PhaseMobilityKernel( 2, 1 ); -INST_PhaseMobilityKernel( 3, 1 ); -INST_PhaseMobilityKernel( 4, 1 ); -INST_PhaseMobilityKernel( 5, 1 ); - -INST_PhaseMobilityKernel( 1, 2 ); -INST_PhaseMobilityKernel( 2, 2 ); -INST_PhaseMobilityKernel( 3, 2 ); -INST_PhaseMobilityKernel( 4, 2 ); -INST_PhaseMobilityKernel( 5, 2 ); - -INST_PhaseMobilityKernel( 1, 3 ); -INST_PhaseMobilityKernel( 2, 3 ); -INST_PhaseMobilityKernel( 3, 3 ); -INST_PhaseMobilityKernel( 4, 3 ); -INST_PhaseMobilityKernel( 5, 3 ); - -#undef INST_PhaseMobilityKernel - - } // namespace CompositionalMultiphaseHybridFVMKernels } // namespace geosx diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp index 68b26964e63..9f03002f050 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp @@ -20,12 +20,14 @@ #define GEOSX_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEHYBRIDFVMKERNELS_HPP #include "common/DataTypes.hpp" -#include "constitutive/fluid/layouts.hpp" -#include "constitutive/relativePermeability/layouts.hpp" -#include "constitutive/capillaryPressure/layouts.hpp" -#include "linearAlgebra/interfaces/InterfaceTypes.hpp" -#include "mesh/MeshLevel.hpp" + +#include "constitutive/fluid/MultiFluidBase.hpp" #include "constitutive/permeability/PermeabilityBase.hpp" +#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" +#include "mesh/ElementRegionManager.hpp" +#include "mesh/ObjectManagerBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" + namespace geosx { @@ -78,7 +80,7 @@ struct UpwindingHelper * @param[out] dUpwPhaseViscCoef_dCompDens the derivatives of the upwind viscous transport coef wrt component density at this face * @param[out] upwViscDofNumber the dof number of the upwind cell at this face */ - template< localIndex NC, localIndex NP > + template< integer NC, integer NP > GEOSX_HOST_DEVICE static void upwindViscousCoefficient( localIndex const (&localIds)[ 3 ], @@ -124,7 +126,7 @@ struct UpwindingHelper * @param[inout] dUpwPhaseGravCoef_dPres the derivative of the upwinded buoyancy transport coefficient wrt pressure * @param[inout] dUpwPhaseGravCoef_dCompDens the derivative of the upwinded buoyancy transport coefficient wrt component density */ - template< localIndex NC, localIndex NP > + template< integer NC, integer NP > GEOSX_HOST_DEVICE static void upwindBuoyancyCoefficient( localIndex const (&localIds)[ 3 ], @@ -165,7 +167,7 @@ struct UpwindingHelper * @param[inout] dPhaseGravTerm_dCompDens the derivatives of the gravCoef multiplied by the difference in phase densities wrt component * density */ - template< localIndex NC, localIndex NP > + template< integer NC, integer NP > GEOSX_HOST_DEVICE static void computePhaseGravTerm( localIndex const (&localIds)[ 3 ], @@ -191,7 +193,7 @@ struct UpwindingHelper * @param[inout] dTotalMob_dPres the derivative of the upwinded total mobility wrt pressure * @param[inout] dTotalMob_dCompDens the derivative of the upwinded total mobility wrt component density */ - template< localIndex NC, localIndex NP > + template< integer NC, integer NP > GEOSX_HOST_DEVICE static void computeUpwindedTotalMobility( localIndex const (&localIds)[ 3 ], @@ -235,7 +237,7 @@ struct UpwindingHelper * @param[out] totalMobIds for each phase, triplet of indices of the upwind element * @param[out] totalMobPos for each phase, flag specifying with the upwind element is local or neighbor */ - template< localIndex NP > + template< integer NP > GEOSX_HOST_DEVICE static void setIndicesForTotalMobilityUpwinding( localIndex const (&localIds)[ 3 ], @@ -284,7 +286,7 @@ struct AssemblerKernelHelper * @param[out] dOneSidedVolFlux_dFacePres the derivatives of the vol fluxes wrt to this element's face pressures * @param[out] dOneSidedVolFlux_dCompDens the derivatives of the vol fluxes wrt to this element's component density */ - template< localIndex NF, localIndex NC, localIndex NP > + template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE static void applyGradient( arrayView1d< real64 const > const & facePres, @@ -336,7 +338,7 @@ struct AssemblerKernelHelper * @param[inout] localMatrix the Jacobian matrix * @param[inout] localRhs the residual */ - template< localIndex NF, localIndex NC, localIndex NP > + template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE static void assembleFluxDivergence( localIndex const (&localIds)[ 3 ], @@ -391,7 +393,7 @@ struct AssemblerKernelHelper * @param[inout] dofColIndicesElemVars degrees of freedom of the cells involved in the flux divergence * @param[inout] dofColIndicesFaceVars degrees of freedom of the faces involved in the flux divergence */ - template< localIndex NF, localIndex NC, localIndex NP > + template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE static void assembleViscousFlux( localIndex const ifaceLoc, @@ -429,7 +431,7 @@ struct AssemblerKernelHelper * @param[inout] dDivMassFluxes_dElemVars the derivatives of the flux divergence wrt the element centered vars (pres and comp dens) * @param[inout] dofColIndicesElemVars degrees of freedom of the cells involved in the flux divergence */ - template< localIndex NF, localIndex NC, localIndex NP > + template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE static void assembleBuoyancyFlux( localIndex const ifaceLoc, @@ -457,7 +459,7 @@ struct AssemblerKernelHelper * @param[inout] matrix the jacobian matrix * @param[inout] rhs the residual */ - template< localIndex NF, localIndex NC, localIndex NP > + template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE static void assembleFaceConstraints( arrayView1d< globalIndex const > const & faceDofNumber, @@ -528,7 +530,7 @@ struct AssemblerKernel * @param[inout] matrix the system matrix * @param[inout] rhs the system right-hand side vector */ - template< localIndex NF, localIndex NC, localIndex NP > + template< integer NF, integer NC, integer NP > GEOSX_HOST_DEVICE static void compute( localIndex const er, localIndex const esr, localIndex const ei, @@ -623,7 +625,7 @@ struct FluxKernel * @param[inout] matrix the system matrix * @param[inout] rhs the system right-hand side vector */ - template< localIndex NF, localIndex NC, localIndex NP, typename IP_TYPE > + template< integer NF, integer NC, integer NP, typename IP_TYPE > static void launch( localIndex er, localIndex esr, CellElementSubRegion const & subRegion, @@ -667,59 +669,202 @@ struct FluxKernel /******************************** PhaseMobilityKernel ********************************/ /** - * @brief Functions to compute phase mobilities and derivatives from viscosity and relperm + * @class PhaseMobilityKernel + * @tparam NUM_COMP number of fluid components + * @tparam NUM_PHASE number of fluid phases + * @brief Define the interface for the property kernel in charge of computing the phase mobilities */ -struct PhaseMobilityKernel +template< integer NUM_COMP, integer NUM_PHASE > +class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { - template< localIndex NC, localIndex NP > +public: + + using Base = CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; + using Base::numComp; + + /// Compile time value for the number of phases + static constexpr integer numPhase = NUM_PHASE; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] relperm the relperm model + */ + PhaseMobilityKernel( ObjectManagerBase & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + : Base(), + m_phaseVolFrac( subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFraction >() ), + m_dPhaseVolFrac_dPres( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dPressure >() ), + m_dPhaseVolFrac_dComp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dGlobalCompDensity >() ), + m_dCompFrac_dCompDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity >() ), + m_phaseVisc( fluid.phaseViscosity() ), + m_dPhaseVisc_dPres( fluid.dPhaseViscosity_dPressure() ), + m_dPhaseVisc_dComp( fluid.dPhaseViscosity_dGlobalCompFraction() ), + m_phaseRelPerm( relperm.phaseRelPerm() ), + m_dPhaseRelPerm_dPhaseVolFrac( relperm.dPhaseRelPerm_dPhaseVolFraction() ), + m_phaseMob( subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseMobility >() ), + m_dPhaseMob_dPres( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dPressure >() ), + m_dPhaseMob_dComp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dGlobalCompDensity >() ) + {} + + /** + * @brief Compute the phase mobilities in an element + * @tparam FUNC the type of the function that can be used to customize the kernel + * @param[in] ei the element index + * @param[in] phaseMobilityKernelOp the function used to customize the kernel + */ + template< typename FUNC = CompositionalMultiphaseBaseKernels::NoOpFunc > GEOSX_HOST_DEVICE - static void - compute( arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dCompFrac_dCompDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseVisc, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseVisc_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseVisc_dComp, - arraySlice1d< real64 const, relperm::USD_RELPERM - 2 > const & phaseRelPerm, - arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > const & dPhaseRelPerm_dPhaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & phaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & dPhaseVolFrac_dPres, - arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const & dPhaseVolFrac_dComp, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & phaseMob, - arraySlice1d< real64, compflow::USD_PHASE - 1 > const & dPhaseMob_dPres, - arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const & dPhaseMob_dComp ); - - template< localIndex NC, localIndex NP > - static void - launch( localIndex const size, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ); - - template< localIndex NC, localIndex NP > - static void - launch( SortedArrayView< localIndex const > const & targetSet, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseVisc_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseVisc_dComp, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm, - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dComp, - arrayView2d< real64, compflow::USD_PHASE > const & phaseMob, - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseMob_dPres, - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseMob_dComp ); + void compute( localIndex const ei, + FUNC && phaseMobilityKernelOp = CompositionalMultiphaseBaseKernels::NoOpFunc{} ) const + { + arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const phaseVisc = m_phaseVisc[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const dPhaseVisc_dPres = m_dPhaseVisc_dPres[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseVisc_dComp = m_dPhaseVisc_dComp[ei][0]; + arraySlice1d< real64 const, relperm::USD_RELPERM - 2 > const phaseRelPerm = m_phaseRelPerm[ei][0]; + arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > const dPhaseRelPerm_dPhaseVolFrac = m_dPhaseRelPerm_dPhaseVolFrac[ei][0]; + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const phaseVolFrac = m_phaseVolFrac[ei]; + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const dPhaseVolFrac_dPres = m_dPhaseVolFrac_dPres[ei]; + arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const dPhaseVolFrac_dComp = m_dPhaseVolFrac_dComp[ei]; + arraySlice1d< real64, compflow::USD_PHASE - 1 > const phaseMob = m_phaseMob[ei]; + arraySlice1d< real64, compflow::USD_PHASE - 1 > const dPhaseMob_dPres = m_dPhaseMob_dPres[ei]; + arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const dPhaseMob_dComp = m_dPhaseMob_dComp[ei]; + + real64 dRelPerm_dC[numComp]{}; + real64 dVisc_dC[numComp]{}; + + for( integer ip = 0; ip < numPhase; ++ip ) + { + // compute the phase mobility only if the phase is present + bool const phaseExists = (phaseVolFrac[ip] > 0); + if( !phaseExists ) + { + phaseMob[ip] = 0.; + dPhaseMob_dPres[ip] = 0.; + for( integer jc = 0; jc < numComp; ++jc ) + { + dPhaseMob_dComp[ip][jc] = 0.; + } + continue; + } + + real64 const viscosity = phaseVisc[ip]; + real64 const dVisc_dP = dPhaseVisc_dPres[ip]; + applyChainRule( numComp, dCompFrac_dCompDens, dPhaseVisc_dComp[ip], dVisc_dC ); + + real64 const relPerm = phaseRelPerm[ip]; + real64 dRelPerm_dP = 0.0; + for( integer ic = 0; ic < numComp; ++ic ) + { + dRelPerm_dC[ic] = 0.0; + } + + for( integer jp = 0; jp < numPhase; ++jp ) + { + real64 const dRelPerm_dS = dPhaseRelPerm_dPhaseVolFrac[ip][jp]; + dRelPerm_dP += dRelPerm_dS * dPhaseVolFrac_dPres[jp]; + + for( integer jc = 0; jc < numComp; ++jc ) + { + dRelPerm_dC[jc] += dRelPerm_dS * dPhaseVolFrac_dComp[jp][jc]; + } + } + + real64 const mobility = relPerm / viscosity; + + phaseMob[ip] = mobility; + dPhaseMob_dPres[ip] = dRelPerm_dP / viscosity + - mobility * dVisc_dP / viscosity; + + // compositional derivatives + for( integer jc = 0; jc < numComp; ++jc ) + { + dPhaseMob_dComp[ip][jc] = dRelPerm_dC[jc] / viscosity + - mobility * dVisc_dC[jc] / viscosity; + } + + // call the lambda in the phase loop to allow the reuse of the relperm, viscosity, and mobility + // possible use: assemble the derivatives wrt temperature + phaseMobilityKernelOp( ip, phaseMob[ip], dPhaseMob_dPres[ip], dPhaseMob_dComp[ip] ); + } + } + + +protected: + + // inputs + + /// Views on the phase volume fractions + arrayView2d< real64 const, compflow::USD_PHASE > m_phaseVolFrac; + arrayView2d< real64 const, compflow::USD_PHASE > m_dPhaseVolFrac_dPres; + arrayView3d< real64 const, compflow::USD_PHASE_DC > m_dPhaseVolFrac_dComp; + arrayView3d< real64 const, compflow::USD_COMP_DC > m_dCompFrac_dCompDens; + + /// Views on the phase viscosities + arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseVisc; + arrayView3d< real64 const, multifluid::USD_PHASE > m_dPhaseVisc_dPres; + arrayView4d< real64 const, multifluid::USD_PHASE_DC > m_dPhaseVisc_dComp; + + /// Views on the phase relative permeabilities + arrayView3d< real64 const, relperm::USD_RELPERM > m_phaseRelPerm; + arrayView4d< real64 const, relperm::USD_RELPERM_DS > m_dPhaseRelPerm_dPhaseVolFrac; + + // outputs + + /// Views on the phase mobilities + arrayView2d< real64, compflow::USD_PHASE > m_phaseMob; + arrayView2d< real64, compflow::USD_PHASE > m_dPhaseMob_dPres; + arrayView3d< real64, compflow::USD_PHASE_DC > m_dPhaseMob_dComp; + }; +/** + * @class PhaseMobilityKernelFactory + */ +class PhaseMobilityKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numComp the number of fluid components + * @param[in] numPhase the number of fluid phases + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] relperm the relperm model + */ + template< typename POLICY > + static void + createAndLaunch( integer const numComp, + integer const numPhase, + ObjectManagerBase & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + { + if( numPhase == 2 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + } +}; /******************************** ResidualNormKernel ********************************/ @@ -765,7 +910,7 @@ struct ResidualNormKernel { // compute a normalizer to obtain a dimensionless norm real64 sumMobOld = 0.0; - for( localIndex ip = 0; ip < numPhases; ++ip ) + for( integer ip = 0; ip < numPhases; ++ip ) { sumMobOld += phaseMobOld[er][esr][ei][ip]; } @@ -826,7 +971,7 @@ struct SolutionCheckKernel struct PrecomputeKernel { - template< typename IP_TYPE, localIndex NF > + template< typename IP_TYPE, integer NF > static void launch( localIndex const subRegionSize, localIndex const faceManagerSize, @@ -859,7 +1004,7 @@ struct PrecomputeKernel lengthTolerance, transMatrix ); - for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) { mimFaceGravCoefNumerator[elemToFaces[ei][ifaceLoc]] += elemGravCoef[ei] * transMatrix[ifaceLoc][ifaceLoc]; mimFaceGravCoefDenominator[elemToFaces[ei][ifaceLoc]] += transMatrix[ifaceLoc][ifaceLoc]; @@ -901,7 +1046,7 @@ void KernelLaunchSelectorFaceSwitch( T value, LAMBDA && lambda ) } // namespace internal template< typename KERNELWRAPPER, typename IP_TYPE, typename ... ARGS > -void KernelLaunchSelector( localIndex numFacesInElem, localIndex numComps, localIndex numPhases, ARGS && ... args ) +void KernelLaunchSelector( integer numFacesInElem, integer numComps, integer numPhases, ARGS && ... args ) { // Ideally this would be inside the dispatch, but it breaks on Summit with GCC 9.1.0 and CUDA 11.0.3. if( numPhases == 2 ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.cpp deleted file mode 100644 index d1f01a0c44d..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.cpp +++ /dev/null @@ -1,27 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2018-2020 TotalEnergies - * Copyright (c) 2019- GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file SinglePhaseBaseKernels.cpp - */ - -#include "SinglePhaseBaseKernels.hpp" - -namespace geosx -{ - -namespace SinglePhaseBaseKernels -{} // namespace SinglePhaseBaseKernels - -} // namespace geosx diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index f735a8ed030..07161ab32c2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -443,25 +443,10 @@ void CompositionalMultiphaseWell::updateComponentFraction( WellElementSubRegion { GEOSX_MARK_FUNCTION; - // outputs - arrayView2d< real64, compflow::USD_COMP > const & compFrac = - subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::globalCompFractionString() ); - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens = - subRegion.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( viewKeyStruct::dGlobalCompFraction_dGlobalCompDensityString() ); - - // inputs - arrayView2d< real64 const, compflow::USD_COMP > const & compDens = - subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::globalCompDensityString() ); - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens = - subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::deltaGlobalCompDensityString() ); - - CompositionalMultiphaseBaseKernels::KernelLaunchSelector1< CompositionalMultiphaseBaseKernels::ComponentFractionKernel - >( numFluidComponents(), - subRegion.size(), - compDens, - dCompDens, - compFrac, - dCompFrac_dCompDens ); + CompositionalMultiphaseBaseKernels:: + ComponentFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + subRegion ); } @@ -759,94 +744,28 @@ void CompositionalMultiphaseWell::updatePhaseVolumeFraction( WellElementSubRegio { GEOSX_MARK_FUNCTION; - // outputs - - arrayView2d< real64, compflow::USD_PHASE > const & phaseVolFrac = - subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::phaseVolumeFractionString() ); - arrayView2d< real64, compflow::USD_PHASE > const & dPhaseVolFrac_dPres = - subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::dPhaseVolumeFraction_dPressureString() ); - arrayView3d< real64, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dCompDens = - subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( viewKeyStruct::dPhaseVolumeFraction_dGlobalCompDensityString() ); - - // inputs - - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens = - subRegion.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( viewKeyStruct::dGlobalCompFraction_dGlobalCompDensityString() ); - arrayView2d< real64 const, compflow::USD_COMP > const & compDens = - subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::globalCompDensityString() ); - arrayView2d< real64 const, compflow::USD_COMP > const & dCompDens = - subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::deltaGlobalCompDensityString() ); - MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, m_fluidModelNames[targetIndex] ); - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseFrac = fluid.phaseFraction(); - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseFrac_dPres = fluid.dPhaseFraction_dPressure(); - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseFrac_dComp = fluid.dPhaseFraction_dGlobalCompFraction(); - - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDens = fluid.phaseDensity(); - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseDens_dPres = fluid.dPhaseDensity_dPressure(); - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseDens_dComp = fluid.dPhaseDensity_dGlobalCompFraction(); + CompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + subRegion, + fluid ); - CompositionalMultiphaseBaseKernels::KernelLaunchSelector2< CompositionalMultiphaseBaseKernels::PhaseVolumeFractionKernel - >( numFluidComponents(), numFluidPhases(), - subRegion.size(), - compDens, - dCompDens, - dCompFrac_dCompDens, - phaseDens, - dPhaseDens_dPres, - dPhaseDens_dComp, - phaseFrac, - dPhaseFrac_dPres, - dPhaseFrac_dComp, - phaseVolFrac, - dPhaseVolFrac_dPres, - dPhaseVolFrac_dCompDens ); } void CompositionalMultiphaseWell::updateTotalMassDensity( WellElementSubRegion & subRegion, localIndex const targetIndex ) const { - // outputs - - arrayView1d< real64 > const & totalMassDens = - subRegion.getReference< array1d< real64 > >( viewKeyStruct::totalMassDensityString() ); - arrayView1d< real64 > const & dTotalMassDens_dPres = - subRegion.getReference< array1d< real64 > >( viewKeyStruct::dTotalMassDensity_dPressureString() ); - arrayView2d< real64, compflow::USD_FLUID_DC > const & dTotalMassDens_dCompDens = - subRegion.getReference< array2d< real64, compflow::LAYOUT_FLUID_DC > >( viewKeyStruct::dTotalMassDensity_dGlobalCompDensityString() ); - - // inputs - - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac = - subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::phaseVolumeFractionString() ); - arrayView2d< real64 const, compflow::USD_PHASE > const & dPhaseVolFrac_dPres = - subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::dPhaseVolumeFraction_dPressureString() ); - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dPhaseVolFrac_dCompDens = - subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( viewKeyStruct::dPhaseVolumeFraction_dGlobalCompDensityString() ); - - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dCompFrac_dCompDens = - subRegion.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( viewKeyStruct::dGlobalCompFraction_dGlobalCompDensityString() ); + GEOSX_MARK_FUNCTION; MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, m_fluidModelNames[targetIndex] ); - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseMassDens = fluid.phaseMassDensity(); - arrayView3d< real64 const, multifluid::USD_PHASE > const & dPhaseMassDens_dPres = fluid.dPhaseMassDensity_dPressure(); - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dPhaseMassDens_dComp = fluid.dPhaseMassDensity_dGlobalCompFraction(); - - CompositionalMultiphaseBaseKernels:: - KernelLaunchSelector2< TotalMassDensityKernel >( numFluidComponents(), - numFluidPhases(), - subRegion.size(), - phaseVolFrac, - dPhaseVolFrac_dPres, - dPhaseVolFrac_dCompDens, - dCompFrac_dCompDens, - phaseMassDens, - dPhaseMassDens_dPres, - dPhaseMassDens_dComp, - totalMassDens, - dTotalMassDens_dPres, - dTotalMassDens_dCompDens ); + TotalMassDensityKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + subRegion, + fluid ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp index 9bda15338c5..ebdab02ac9b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp @@ -1755,103 +1755,6 @@ RateInitializationKernel:: } ); } -/******************************** TotalMassDensityKernel ****************************/ - -template< localIndex NC, localIndex NP > -GEOSX_HOST_DEVICE -void -TotalMassDensityKernel:: - compute( arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & phaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & dPhaseVolFrac_dPres, - arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const & dPhaseVolFrac_dCompDens, - arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dCompFrac_dCompDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseMassDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseMassDens_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseMassDens_dComp, - real64 & totalMassDens, - real64 & dTotalMassDens_dPres, - arraySlice1d< real64, compflow::USD_FLUID_DC - 1 > const & dTotalMassDens_dCompDens ) -{ - real64 dMassDens_dC[NC]{}; - - totalMassDens = 0.0; - dTotalMassDens_dPres = 0.0; - for( localIndex ic = 0; ic < NC; ++ic ) - { - dTotalMassDens_dCompDens[ic] = 0.0; - } - - for( localIndex ip = 0; ip < NP; ++ip ) - { - totalMassDens += phaseVolFrac[ip] * phaseMassDens[ip]; - dTotalMassDens_dPres += dPhaseVolFrac_dPres[ip] * phaseMassDens[ip] - + phaseVolFrac[ip] * dPhaseMassDens_dPres[ip]; - - applyChainRule( NC, dCompFrac_dCompDens, dPhaseMassDens_dComp[ip], dMassDens_dC ); - for( localIndex ic = 0; ic < NC; ++ic ) - { - dTotalMassDens_dCompDens[ic] += dPhaseVolFrac_dCompDens[ip][ic] * phaseMassDens[ip] - + phaseVolFrac[ip] * dMassDens_dC[ic]; - } - } -} - -template< localIndex NC, localIndex NP > -void -TotalMassDensityKernel:: - launch( localIndex const size, - arrayView2d< real64 const, compflow::USD_PHASE > const & wellElemPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dWellElemPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dWellElemPhaseVolFrac_dCompDens, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dWellElemCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & wellElemPhaseMassDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dWellElemPhaseMassDens_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dWellElemPhaseMassDens_dComp, - arrayView1d< real64 > const & wellElemTotalMassDens, - arrayView1d< real64 > const & dWellElemTotalMassDens_dPres, - arrayView2d< real64, compflow::USD_FLUID_DC > const & dWellElemTotalMassDens_dCompDens ) -{ - forAll< parallelDevicePolicy<> >( size, [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) - { - compute< NC, NP >( wellElemPhaseVolFrac[iwelem], - dWellElemPhaseVolFrac_dPres[iwelem], - dWellElemPhaseVolFrac_dCompDens[iwelem], - dWellElemCompFrac_dCompDens[iwelem], - wellElemPhaseMassDens[iwelem][0], - dWellElemPhaseMassDens_dPres[iwelem][0], - dWellElemPhaseMassDens_dComp[iwelem][0], - wellElemTotalMassDens[iwelem], - dWellElemTotalMassDens_dPres[iwelem], - dWellElemTotalMassDens_dCompDens[iwelem] ); - } ); -} - -#define INST_TotalMassDensityKernel( NC, NP ) \ - template \ - void TotalMassDensityKernel:: \ - launch< NC, NP >( localIndex const size, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & wellElemPhaseVolFrac, \ - arrayView2d< real64 const, compflow::USD_PHASE > const & dWellElemPhaseVolFrac_dPres, \ - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dWellElemPhaseVolFrac_dCompDens, \ - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dWellElemCompFrac_dCompDens, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & wellElemPhaseMassDens, \ - arrayView3d< real64 const, multifluid::USD_PHASE > const & dWellElemPhaseMassDens_dPres, \ - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dWellElemPhaseMassDens_dComp, \ - arrayView1d< real64 > const & wellElemTotalMassDens, \ - arrayView1d< real64 > const & dWellElemTotalMassDens_dPres, \ - arrayView2d< real64, compflow::USD_FLUID_DC > const & dWellElemTotalMassDens_dCompDens ) - -INST_TotalMassDensityKernel( 1, 2 ); -INST_TotalMassDensityKernel( 2, 2 ); -INST_TotalMassDensityKernel( 3, 2 ); -INST_TotalMassDensityKernel( 4, 2 ); -INST_TotalMassDensityKernel( 5, 2 ); -INST_TotalMassDensityKernel( 1, 3 ); -INST_TotalMassDensityKernel( 2, 3 ); -INST_TotalMassDensityKernel( 3, 3 ); -INST_TotalMassDensityKernel( 4, 3 ); -INST_TotalMassDensityKernel( 5, 3 ); - } // end namespace CompositionalMultiphaseWellKernels diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index 6ad43404883..c43de914887 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -22,6 +22,7 @@ #include "common/DataTypes.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "constitutive/fluid/MultiFluidBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp" #include "physicsSolvers/fluidFlow/wells/WellControls.hpp" @@ -439,39 +440,162 @@ struct RateInitializationKernel /******************************** TotalMassDensityKernel ****************************/ -struct TotalMassDensityKernel +/** + * @class TotalMassDensityKernel + * @tparam NUM_COMP number of fluid components + * @tparam NUM_PHASE number of fluid phases + * @brief Define the interface for the property kernel in charge of computing the total mass density + */ +template< integer NUM_COMP, integer NUM_PHASE > +class TotalMassDensityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { +public: - template< localIndex NC, localIndex NP > + using Base = CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; + // TODO: get rid of this as soon as extrinsicData is available for the wells + using keys = CompositionalMultiphaseWell::viewKeyStruct; + + using Base::numComp; + + /// Compile time value for the number of phases + static constexpr integer numPhase = NUM_PHASE; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + TotalMassDensityKernel( dataRepository::Group & subRegion, + MultiFluidBase const & fluid ) + : Base(), + m_phaseVolFrac( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( keys::phaseVolumeFractionString() ) ), + m_dPhaseVolFrac_dPres( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( keys::dPhaseVolumeFraction_dPressureString() ) ), + m_dPhaseVolFrac_dCompDens( subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( keys::dPhaseVolumeFraction_dGlobalCompDensityString() ) ), + m_dCompFrac_dCompDens( subRegion.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( keys::dGlobalCompFraction_dGlobalCompDensityString() ) ), + m_phaseMassDens( fluid.phaseMassDensity() ), + m_dPhaseMassDens_dPres( fluid.dPhaseMassDensity_dPressure() ), + m_dPhaseMassDens_dComp( fluid.dPhaseMassDensity_dGlobalCompFraction() ), + m_totalMassDens( subRegion.getReference< array1d< real64 > >( keys::totalMassDensityString() ) ), + m_dTotalMassDens_dPres( subRegion.getReference< array1d< real64 > >( keys::dTotalMassDensity_dPressureString() ) ), + m_dTotalMassDens_dCompDens( subRegion.getReference< array2d< real64, compflow::LAYOUT_FLUID_DC > >( keys::dTotalMassDensity_dGlobalCompDensityString() ) ) + {} + + /** + * @brief Compute the total mass density in an element + * @tparam FUNC the type of the function that can be used to customize the kernel + * @param[in] ei the element index + * @param[in] totalMassDensityKernelOp the function used to customize the kernel + */ + template< typename FUNC = CompositionalMultiphaseBaseKernels::NoOpFunc > GEOSX_HOST_DEVICE - static void - compute( arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & phaseVolFrac, - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & dPhaseVolFrac_dPres, - arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const & dPhaseVolFrac_dCompDens, - arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dCompFrac_dCompDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & phaseMassDens, - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & dPhaseMassDens_dPres, - arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dPhaseMassDens_dComp, - real64 & totalMassDens, - real64 & dTotalMassDens_dPres, - arraySlice1d< real64, compflow::USD_FLUID_DC - 1 > const & dTotalMassDens_dCompDens ); + void compute( localIndex const ei, + FUNC && totalMassDensityKernelOp = CompositionalMultiphaseBaseKernels::NoOpFunc{} ) const + { + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac = m_phaseVolFrac[ei]; + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > dPhaseVolFrac_dPres = m_dPhaseVolFrac_dPres[ei]; + arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > dPhaseVolFrac_dCompDens = m_dPhaseVolFrac_dCompDens[ei]; + arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseMassDens = m_phaseMassDens[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > dPhaseMassDens_dPres = m_dPhaseMassDens_dPres[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseMassDens_dComp = m_dPhaseMassDens_dComp[ei][0]; + real64 & totalMassDens = m_totalMassDens[ei]; + real64 & dTotalMassDens_dPres = m_dTotalMassDens_dPres[ei]; + arraySlice1d< real64, compflow::USD_FLUID_DC - 1 > dTotalMassDens_dCompDens = m_dTotalMassDens_dCompDens[ei]; + + real64 dMassDens_dC[numComp]{}; + + totalMassDens = 0.0; + dTotalMassDens_dPres = 0.0; + for( integer ic = 0; ic < numComp; ++ic ) + { + dTotalMassDens_dCompDens[ic] = 0.0; + } - template< localIndex NC, localIndex NP > - static void - launch( localIndex const size, - arrayView2d< real64 const, compflow::USD_PHASE > const & wellElemPhaseVolFrac, - arrayView2d< real64 const, compflow::USD_PHASE > const & dWellElemPhaseVolFrac_dPres, - arrayView3d< real64 const, compflow::USD_PHASE_DC > const & dWellElemPhaseVolFrac_dCompDens, - arrayView3d< real64 const, compflow::USD_COMP_DC > const & dWellElemCompFrac_dCompDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & wellElemPhaseMassDens, - arrayView3d< real64 const, multifluid::USD_PHASE > const & dWellElemPhaseMassDens_dPres, - arrayView4d< real64 const, multifluid::USD_PHASE_DC > const & dWellElemPhaseMassDens_dComp, - arrayView1d< real64 > const & wellElemTotalMassDens, - arrayView1d< real64 > const & dWellElemTotalMassDens_dPres, - arrayView2d< real64, compflow::USD_FLUID_DC > const & dWellElemTotalMassDens_dCompDens ); + for( integer ip = 0; ip < numPhase; ++ip ) + { + totalMassDens += phaseVolFrac[ip] * phaseMassDens[ip]; + dTotalMassDens_dPres += dPhaseVolFrac_dPres[ip] * phaseMassDens[ip] + + phaseVolFrac[ip] * dPhaseMassDens_dPres[ip]; + + applyChainRule( numComp, dCompFrac_dCompDens, dPhaseMassDens_dComp[ip], dMassDens_dC ); + for( integer ic = 0; ic < numComp; ++ic ) + { + dTotalMassDens_dCompDens[ic] += dPhaseVolFrac_dCompDens[ip][ic] * phaseMassDens[ip] + + phaseVolFrac[ip] * dMassDens_dC[ic]; + } + totalMassDensityKernelOp( ip, totalMassDens, dTotalMassDens_dPres, dTotalMassDens_dCompDens ); + } + } + +protected: + + // inputs + + /// Views on phase volume fractions + arrayView2d< real64 const, compflow::USD_PHASE > m_phaseVolFrac; + arrayView2d< real64 const, compflow::USD_PHASE > m_dPhaseVolFrac_dPres; + arrayView3d< real64 const, compflow::USD_PHASE_DC > m_dPhaseVolFrac_dCompDens; + arrayView3d< real64 const, compflow::USD_COMP_DC > m_dCompFrac_dCompDens; + + /// Views on phase mass densities + arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseMassDens; + arrayView3d< real64 const, multifluid::USD_PHASE > m_dPhaseMassDens_dPres; + arrayView4d< real64 const, multifluid::USD_PHASE_DC > m_dPhaseMassDens_dComp; + + // outputs + + /// Views on total mass densities + arrayView1d< real64 > m_totalMassDens; + arrayView1d< real64 > m_dTotalMassDens_dPres; + arrayView2d< real64, compflow::USD_FLUID_DC > m_dTotalMassDens_dCompDens; + +}; + +/** + * @class TotalMassDensityKernelFactory + */ +class TotalMassDensityKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numComp the number of fluid components + * @param[in] numPhase the number of fluid phases + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + template< typename POLICY > + static void + createAndLaunch( integer const numComp, + integer const numPhase, + dataRepository::Group & subRegion, + MultiFluidBase const & fluid ) + { + if( numPhase == 2 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + TotalMassDensityKernel< NUM_COMP, 2 > kernel( subRegion, fluid ); + TotalMassDensityKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + TotalMassDensityKernel< NUM_COMP, 3 > kernel( subRegion, fluid ); + TotalMassDensityKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + } }; + /******************************** ResidualNormKernel ********************************/ struct ResidualNormKernel