From 89f5c8d6793d6713c99013bae44951e324d442e2 Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Tue, 30 Nov 2021 00:23:23 -0800 Subject: [PATCH 1/3] implemented new kernel interface for properties --- .../PVTFunctions/EzrokhiBrineDensity.hpp | 2 +- .../PVTFunctions/EzrokhiBrineViscosity.hpp | 2 +- .../fluidFlow/CompositionalMultiphaseBase.cpp | 77 +-- .../CompositionalMultiphaseBaseKernels.cpp | 293 +--------- .../CompositionalMultiphaseBaseKernels.hpp | 543 ++++++++++++++---- .../fluidFlow/CompositionalMultiphaseFVM.cpp | 65 +-- .../CompositionalMultiphaseFVMKernels.cpp | 223 ------- .../CompositionalMultiphaseFVMKernels.hpp | 274 +++++++-- .../CompositionalMultiphaseHybridFVM.cpp | 59 +- ...ompositionalMultiphaseHybridFVMKernels.cpp | 197 ------- ...ompositionalMultiphaseHybridFVMKernels.hpp | 252 ++++++-- .../wells/CompositionalMultiphaseWell.cpp | 120 +--- .../CompositionalMultiphaseWellKernels.cpp | 97 ---- .../CompositionalMultiphaseWellKernels.hpp | 178 +++++- 14 files changed, 1050 insertions(+), 1332 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/PVTFunctions/EzrokhiBrineDensity.hpp b/src/coreComponents/constitutive/fluid/PVTFunctions/EzrokhiBrineDensity.hpp index e84cf70a9dc..c2f5e692149 100644 --- a/src/coreComponents/constitutive/fluid/PVTFunctions/EzrokhiBrineDensity.hpp +++ b/src/coreComponents/constitutive/fluid/PVTFunctions/EzrokhiBrineDensity.hpp @@ -238,7 +238,7 @@ void EzrokhiBrineDensityUpdate::compute( real64 const & pressure, value = waterDensity * exponentPowered; - real64 const dValueCoef = log( 10 ) * value; + real64 const dValueCoef = LvArray::math::log( 10 ) * value; real64 const dValue_dPhaseComp = dValueCoef * exponent_dPhaseComp; dValue_dPressure = dValueCoef * exponent_dPressure + waterDensity_dPressure * exponentPowered; diff --git a/src/coreComponents/constitutive/fluid/PVTFunctions/EzrokhiBrineViscosity.hpp b/src/coreComponents/constitutive/fluid/PVTFunctions/EzrokhiBrineViscosity.hpp index 5189c4870c9..cbfa6ee1151 100644 --- a/src/coreComponents/constitutive/fluid/PVTFunctions/EzrokhiBrineViscosity.hpp +++ b/src/coreComponents/constitutive/fluid/PVTFunctions/EzrokhiBrineViscosity.hpp @@ -210,7 +210,7 @@ void EzrokhiBrineViscosityUpdate::compute( real64 const & pressure, real64 const exponentPowered = pow( 10, exponent ); value = waterVisc * exponentPowered; - real64 const dValueCoef = log( 10 ) * value; + real64 const dValueCoef = LvArray::math::log( 10 ) * value; real64 const dValue_dPhaseComp = dValueCoef * exponent_dPhaseComp; dValue_dPressure = dValueCoef * exponent_dPressure; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index f7f2f661dc0..3fc7c779a23 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -371,27 +371,11 @@ void CompositionalMultiphaseBase::updateComponentFraction( Group & dataGroup ) c { GEOSX_MARK_FUNCTION; - // outputs + ComponentFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> > + ( m_numComponents, + dataGroup ); - arrayView2d< real64, compflow::USD_COMP > const & compFrac = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::globalCompFractionString() ); - - arrayView3d< real64, compflow::USD_COMP_DC > const & dCompFrac_dCompDens = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( viewKeyStruct::dGlobalCompFraction_dGlobalCompDensityString() ); - - // inputs - arrayView2d< real64 const, compflow::USD_COMP > const compDens = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::globalCompDensityString() ); - - arrayView2d< real64 const, compflow::USD_COMP > const dCompDens = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::deltaGlobalCompDensityString() ); - - KernelLaunchSelector1< ComponentFractionKernel >( m_numComponents, - dataGroup.size(), - compDens, - dCompDens, - compFrac, - dCompFrac_dCompDens ); } void CompositionalMultiphaseBase::updatePhaseVolumeFraction( Group & dataGroup, @@ -399,52 +383,17 @@ void CompositionalMultiphaseBase::updatePhaseVolumeFraction( Group & dataGroup, { GEOSX_MARK_FUNCTION; - // outputs - - arrayView2d< real64, compflow::USD_PHASE > const phaseVolFrac = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::phaseVolumeFractionString() ); - - arrayView2d< real64, compflow::USD_PHASE > const dPhaseVolFrac_dPres = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::dPhaseVolumeFraction_dPressureString() ); - - arrayView3d< real64, compflow::USD_PHASE_DC > const dPhaseVolFrac_dComp = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( viewKeyStruct::dPhaseVolumeFraction_dGlobalCompDensityString() ); - - // inputs - - arrayView3d< real64 const, compflow::USD_COMP_DC > const dCompFrac_dCompDens = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( viewKeyStruct::dGlobalCompFraction_dGlobalCompDensityString() ); - - arrayView2d< real64 const, compflow::USD_COMP > const compDens = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::globalCompDensityString() ); - - arrayView2d< real64 const, compflow::USD_COMP > const dCompDens = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_COMP > >( viewKeyStruct::deltaGlobalCompDensityString() ); - 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 ); + bool const isIsothermal = true; + PhaseVolumeFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> > + ( isIsothermal, + m_numComponents, + m_numPhases, + dataGroup, + fluid ); + } void CompositionalMultiphaseBase::updateFluidModel( Group & dataGroup, localIndex const targetIndex ) const diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.cpp index 26a0215019e..348646e15ce 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.cpp @@ -17,302 +17,11 @@ */ #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 CompositionalMultiphaseBaseKernels } // namespace geosx diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp index bc4e1cafc18..adc32284a45 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp @@ -24,11 +24,9 @@ #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 "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp" -#include "mesh/ElementRegionManager.hpp" namespace geosx @@ -41,94 +39,440 @@ 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 {} +}; + +/** + * @class PropertyKernelBase + * @tparam NUM_COMP number of fluid components + * @brief Define the base interface for the property update kernels + */ +template< localIndex NUM_COMP > +class PropertyKernelBase +{ +public: + + /// Compile time value for the number of components + static constexpr localIndex numComp = NUM_COMP; - template< localIndex NC > + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + PropertyKernelBase( dataRepository::Group & subRegion ) + { GEOSX_UNUSED_VAR( subRegion ); } + + /** + * @brief Compute the property 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 && kernelOp = NoOpFunc{} ) const + { GEOSX_UNUSED_VAR( ei, kernelOp ); } + + /** + * @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 ) + { + GEOSX_MARK_FUNCTION; - template< localIndex NC > + forAll< POLICY >( numElems, [=] GEOSX_HOST_DEVICE ( localIndex const ei ) + { + kernelComponent.compute( ei ); + } ); + } + + /** + * @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 ) + { + GEOSX_MARK_FUNCTION; + + forAll< POLICY >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const i ) + { + localIndex const ei = targetSet[ i ]; + kernelComponent.compute( ei ); + } ); + } + }; -/******************************** PhaseVolumeFractionKernel ********************************/ +/******************************** 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< localIndex NUM_COMP > +class ComponentFractionKernel : public PropertyKernelBase< NUM_COMP > { - template< localIndex NC, localIndex NP > +public: + + using Base = PropertyKernelBase< NUM_COMP >; + using keys = CompositionalMultiphaseBase::viewKeyStruct; + + using Base::numComp; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + ComponentFractionKernel( dataRepository::Group & subRegion ) + : Base( subRegion ), + m_compDens( subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( keys::globalCompDensityString() ) ), + m_dCompDens( subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( keys::deltaGlobalCompDensityString() ) ), + m_compFrac( subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( keys::globalCompFractionString() ) ), + m_dCompFrac_dCompDens( subRegion.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( keys::dGlobalCompFraction_dGlobalCompDensityString() ) ) + {} + + /** + * @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; + +}; + +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 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( localIndex const numComp, + dataRepository::Group & subRegion ) + { + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + ComponentFractionKernel< NUM_COMP > kernel( subRegion ); + ComponentFractionKernel< NUM_COMP >::template launch< POLICY, ComponentFractionKernel< NUM_COMP > >( 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< localIndex NUM_COMP, localIndex NUM_PHASE > +class PhaseVolumeFractionKernel : public PropertyKernelBase< NUM_COMP > +{ +public: + + using Base = PropertyKernelBase< NUM_COMP >; + using keys = CompositionalMultiphaseBase::viewKeyStruct; + + using Base::numComp; + + /// Compile time value for the number of phases + static constexpr localIndex numPhase = NUM_PHASE; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + PhaseVolumeFractionKernel( dataRepository::Group & subRegion, + MultiFluidBase const & fluid ) + : Base( subRegion ), + 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_dComp( subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( keys::dPhaseVolumeFraction_dGlobalCompDensityString() ) ), + m_compDens( subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( keys::globalCompDensityString() ) ), + m_dCompDens( subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( keys::deltaGlobalCompDensityString() ) ), + m_dCompFrac_dCompDens( subRegion.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( keys::dGlobalCompFraction_dGlobalCompDensityString() ) ), + 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] isIsothermal flag specifying whether the assembly is isothermal or non-isothermal + * @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( bool const isIsothermal, + localIndex const numComp, + localIndex const numPhase, + dataRepository::Group & subRegion, + MultiFluidBase const & fluid ) + { + if( !isIsothermal ) + { GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); } + + if( numPhase == 2 ) + { + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + PhaseVolumeFractionKernel< NUM_COMP, 2 > kernel( subRegion, fluid ); + PhaseVolumeFractionKernel< NUM_COMP, 2 >::template launch< POLICY, PhaseVolumeFractionKernel< NUM_COMP, 2 > >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + PhaseVolumeFractionKernel< NUM_COMP, 3 > kernel( subRegion, fluid ); + PhaseVolumeFractionKernel< NUM_COMP, 3 >::template launch< POLICY, PhaseVolumeFractionKernel< NUM_COMP, 3 > >( subRegion.size(), kernel ); + } ); + } + } +}; /******************************** FluidUpdateKernel ********************************/ @@ -281,21 +625,6 @@ 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 @@ -307,6 +636,8 @@ class ElementBasedAssemblyKernel { public: + using keys = CompositionalMultiphaseBase::viewKeyStruct; + /// Compile time value for the number of components static constexpr localIndex numComp = NUM_COMP; @@ -343,16 +674,16 @@ class ElementBasedAssemblyKernel m_porosityOld( solid.getOldPorosity() ), m_porosityNew( solid.getPorosity() ), m_dPoro_dPres( solid.getDporosity_dPressure() ), - m_dCompFrac_dCompDens( subRegion.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( CompositionalMultiphaseBase::viewKeyStruct::dGlobalCompFraction_dGlobalCompDensityString() ) ), - m_phaseVolFracOld( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( CompositionalMultiphaseBase::viewKeyStruct::phaseVolumeFractionOldString() ) ), - m_phaseVolFrac( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( CompositionalMultiphaseBase::viewKeyStruct::phaseVolumeFractionString() ) ), - m_dPhaseVolFrac_dPres( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( CompositionalMultiphaseBase::viewKeyStruct::dPhaseVolumeFraction_dPressureString() ) ), - m_dPhaseVolFrac_dCompDens( subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( CompositionalMultiphaseBase::viewKeyStruct::dPhaseVolumeFraction_dGlobalCompDensityString() ) ), - m_phaseDensOld( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( CompositionalMultiphaseBase::viewKeyStruct::phaseDensityOldString() ) ), + m_dCompFrac_dCompDens( subRegion.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( keys::dGlobalCompFraction_dGlobalCompDensityString() ) ), + m_phaseVolFracOld( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( keys::phaseVolumeFractionOldString() ) ), + 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_phaseDensOld( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( keys::phaseDensityOldString() ) ), m_phaseDens( fluid.phaseDensity() ), m_dPhaseDens_dPres( fluid.dPhaseDensity_dPressure() ), m_dPhaseDens_dComp( fluid.dPhaseDensity_dGlobalCompFraction() ), - m_phaseCompFracOld( subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_COMP > >( CompositionalMultiphaseBase::viewKeyStruct::phaseComponentFractionOldString() ) ), + m_phaseCompFracOld( subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_COMP > >( keys::phaseComponentFractionOldString() ) ), m_phaseCompFrac( fluid.phaseCompFraction() ), m_dPhaseCompFrac_dPres( fluid.dPhaseCompFraction_dPressure() ), m_dPhaseCompFrac_dComp( fluid.dPhaseCompFraction_dGlobalCompFraction() ), @@ -667,34 +998,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 */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 2ba42d3c61e..9d9efb6a0f8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -515,65 +515,20 @@ void CompositionalMultiphaseFVM::updatePhaseMobility( Group & dataGroup, localIn { 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.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::phaseMobilityString() ); - - arrayView2d< real64, compflow::USD_PHASE > const dPhaseMob_dPres = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::dPhaseMobility_dPressureString() ); - - arrayView3d< real64, compflow::USD_PHASE_DC > const dPhaseMob_dComp = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( viewKeyStruct::dPhaseMobility_dGlobalCompDensityString() ); - - // inputs - - arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::phaseVolumeFractionString() ); - - arrayView2d< real64 const, compflow::USD_PHASE > const dPhaseVolFrac_dPres = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::dPhaseVolumeFraction_dPressureString() ); - - arrayView3d< real64 const, compflow::USD_PHASE_DC > const dPhaseVolFrac_dComp = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( viewKeyStruct::dPhaseVolumeFraction_dGlobalCompDensityString() ); - - arrayView3d< real64 const, compflow::USD_COMP_DC > const dCompFrac_dCompDens = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( viewKeyStruct::dGlobalCompFraction_dGlobalCompDensityString() ); + // 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 ); + bool const isIsothermal = true; + PhaseMobilityKernelFactory:: + createAndLaunch< parallelDevicePolicy<> > + ( isIsothermal, + 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..1e0900af38c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp @@ -33,229 +33,6 @@ 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 > diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp index f0c54d476ce..a91348b839b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp @@ -22,12 +22,11 @@ #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 "fieldSpecification/AquiferBoundaryCondition.hpp" #include "finiteVolume/BoundaryStencil.hpp" -#include "mesh/ElementRegionManager.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" namespace geosx { @@ -40,69 +39,226 @@ 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< localIndex NUM_COMP, localIndex NUM_PHASE > +class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { - template< localIndex NC, localIndex NP > +public: + + using Base = CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; + using keys = CompositionalMultiphaseBase::viewKeyStruct; + + using Base::numComp; + + /// Compile time value for the number of phases + static constexpr localIndex numPhase = NUM_PHASE; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] relperm the relperm model + */ + PhaseMobilityKernel( dataRepository::Group & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + : Base( subRegion ), + 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_dComp( 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_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.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( keys::phaseMobilityString() ) ), + m_dPhaseMob_dPres( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( keys::dPhaseMobility_dPressureString() ) ), + m_dPhaseMob_dComp( subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( keys::dPhaseMobility_dGlobalCompDensityString() ) ) + {} + + /** + * @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 ); + 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( 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< localIndex NC, localIndex NP > +/** + * @class PhaseMobilityKernelFactory + */ +class PhaseMobilityKernelFactory +{ +public: + + /** + * @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] 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 - 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 ); + createAndLaunch( bool const isIsothermal, + localIndex const numComp, + localIndex const numPhase, + dataRepository::Group & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + { + if( !isIsothermal ) + { GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); } + + if( numPhase == 2 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY, PhaseMobilityKernel< NUM_COMP, 2 > >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY, PhaseMobilityKernel< NUM_COMP, 3 > >( subRegion.size(), kernel ); + } ); + } + } }; - /******************************** FluxKernel ********************************/ /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index cad3f9b26b7..93b74b0d7ed 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -890,58 +890,19 @@ void CompositionalMultiphaseHybridFVM::updatePhaseMobility( Group & dataGroup, l { GEOSX_MARK_FUNCTION; - // note that the phase mobility computed here does NOT include phase density - - // outputs - - arrayView2d< real64, compflow::USD_PHASE > const phaseMob = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::phaseMobilityString() ); - - arrayView2d< real64, compflow::USD_PHASE > const dPhaseMob_dPres = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::dPhaseMobility_dPressureString() ); - - arrayView3d< real64, compflow::USD_PHASE_DC > const dPhaseMob_dComp = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( viewKeyStruct::dPhaseMobility_dGlobalCompDensityString() ); - - // inputs - - arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::phaseVolumeFractionString() ); - - arrayView2d< real64 const, compflow::USD_PHASE > const dPhaseVolFrac_dPres = - dataGroup.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( viewKeyStruct::dPhaseVolumeFraction_dPressureString() ); - - arrayView3d< real64 const, compflow::USD_PHASE_DC > const dPhaseVolFrac_dComp = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( viewKeyStruct::dPhaseVolumeFraction_dGlobalCompDensityString() ); - - arrayView3d< real64 const, compflow::USD_COMP_DC > const dCompFrac_dCompDens = - dataGroup.getReference< array3d< real64, compflow::LAYOUT_COMP_DC > >( viewKeyStruct::dGlobalCompFraction_dGlobalCompDensityString() ); - 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 ); + bool const isIsothermal = true; + CompositionalMultiphaseHybridFVMKernels:: + PhaseMobilityKernelFactory:: + createAndLaunch< parallelDevicePolicy<> > + ( isIsothermal, + 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 e03f5a88e67..13d681f6309 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp @@ -1744,203 +1744,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 535e1d358db..26d3081fb97 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp @@ -20,13 +20,11 @@ #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 "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" +#include "constitutive/fluid/MultiFluidBase.hpp" #include "constitutive/permeability/PermeabilityBase.hpp" +#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" + namespace geosx { @@ -668,59 +666,209 @@ 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< localIndex NUM_COMP, localIndex NUM_PHASE > +class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { - template< localIndex NC, localIndex NP > +public: + + using Base = CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; + using keys = CompositionalMultiphaseBase::viewKeyStruct; + + using Base::numComp; + + /// Compile time value for the number of phases + static constexpr localIndex numPhase = NUM_PHASE; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] relperm the relperm model + */ + PhaseMobilityKernel( dataRepository::Group & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + : Base( subRegion ), + 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_dComp( 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_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.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( keys::phaseMobilityString() ) ), + m_dPhaseMob_dPres( subRegion.getReference< array2d< real64, compflow::LAYOUT_PHASE > >( keys::dPhaseMobility_dPressureString() ) ), + m_dPhaseMob_dComp( subRegion.getReference< array3d< real64, compflow::LAYOUT_PHASE_DC > >( keys::dPhaseMobility_dGlobalCompDensityString() ) ) + {} + + /** + * @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 ); + 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; + } - 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 ); + 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; - 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 ); }; +/** + * @class PhaseMobilityKernelFactory + */ +class PhaseMobilityKernelFactory +{ +public: + + /** + * @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] 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( bool const isIsothermal, + localIndex const numComp, + localIndex const numPhase, + dataRepository::Group & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + { + if( !isIsothermal ) + { GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); } + + if( numPhase == 2 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY, PhaseMobilityKernel< NUM_COMP, 2 > >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY, PhaseMobilityKernel< NUM_COMP, 3 > >( subRegion.size(), kernel ); + } ); + } + } +}; /******************************** ResidualNormKernel ********************************/ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 1eb72b0e55c..3f438629e1b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -443,25 +443,11 @@ 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 +745,34 @@ 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(); + bool const isIsothermal = true; + CompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> > + ( isIsothermal, + 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 ); + bool const isIsothermal = true; + TotalMassDensityKernelFactory:: + createAndLaunch< parallelDevicePolicy<> > + ( isIsothermal, + 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..54486562746 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,166 @@ 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< localIndex NUM_COMP, localIndex NUM_PHASE > +class TotalMassDensityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { +public: - template< localIndex NC, localIndex NP > + using Base = CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; + using keys = CompositionalMultiphaseWell::viewKeyStruct; + + using Base::numComp; + + /// Compile time value for the number of phases + static constexpr localIndex 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( subRegion ), + 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 > + 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] isIsothermal flag specifying whether the assembly is isothermal or non-isothermal + * @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 - 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 ); + createAndLaunch( bool const isIsothermal, + localIndex const numComp, + localIndex const numPhase, + dataRepository::Group & subRegion, + MultiFluidBase const & fluid ) + { + if( !isIsothermal ) + { GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); } + if( numPhase == 2 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + TotalMassDensityKernel< NUM_COMP, 2 > kernel( subRegion, fluid ); + TotalMassDensityKernel< NUM_COMP, 2 >::template launch< POLICY, TotalMassDensityKernel< NUM_COMP, 2 > >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + TotalMassDensityKernel< NUM_COMP, 3 > kernel( subRegion, fluid ); + TotalMassDensityKernel< NUM_COMP, 3 >::template launch< POLICY, TotalMassDensityKernel< NUM_COMP, 3 > >( subRegion.size(), kernel ); + } ); + } + } }; + /******************************** ResidualNormKernel ********************************/ struct ResidualNormKernel From db9d26455212e25c0e9b52e89b3059a5fecb6f5a Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Tue, 30 Nov 2021 15:17:58 -0800 Subject: [PATCH 2/3] fixed compilation error --- src/coreComponents/constitutive/fluid/BlackOilFluid.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp b/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp index f8edd50b6b4..9d12043d155 100644 --- a/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp +++ b/src/coreComponents/constitutive/fluid/BlackOilFluid.hpp @@ -292,8 +292,8 @@ BlackOilFluid::KernelWrapper:: { GEOSX_UNUSED_VAR( temperature ); - real64 compMoleFrac[NC_BO]; - real64 dCompMoleFrac_dCompMassFrac[NC_BO][NC_BO]; + real64 compMoleFrac[NC_BO]{}; + real64 dCompMoleFrac_dCompMassFrac[NC_BO][NC_BO]{}; real64 phaseMW[NP_BO]{}; real64 dPhaseMW_dPressure[NP_BO]{}; real64 dPhaseMW_dGlobalCompFraction[NP_BO*NC_BO]{}; From 2ae37eed2e17c251e35bca817c427d563dd7b85f Mon Sep 17 00:00:00 2001 From: Francois Hamon Date: Tue, 14 Dec 2021 14:27:40 -0800 Subject: [PATCH 3/3] addressed Sergey's comments --- .../physicsSolvers/CMakeLists.txt | 2 - .../fluidFlow/CompositionalMultiphaseBase.cpp | 37 ++-- .../CompositionalMultiphaseBaseKernels.cpp | 27 --- .../CompositionalMultiphaseBaseKernels.hpp | 192 +++++++----------- .../fluidFlow/CompositionalMultiphaseFVM.cpp | 13 +- .../CompositionalMultiphaseFVMKernels.cpp | 106 +++++----- .../CompositionalMultiphaseFVMKernels.hpp | 57 +++--- .../CompositionalMultiphaseHybridFVM.cpp | 13 +- ...ompositionalMultiphaseHybridFVMKernels.cpp | 120 +++++------ ...ompositionalMultiphaseHybridFVMKernels.hpp | 55 +++-- .../fluidFlow/SinglePhaseBaseKernels.cpp | 27 --- .../wells/CompositionalMultiphaseWell.cpp | 27 +-- .../CompositionalMultiphaseWellKernels.hpp | 23 +-- 13 files changed, 280 insertions(+), 419 deletions(-) delete mode 100644 src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.cpp delete mode 100644 src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.cpp 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 53a1a6d9f2a..80548de9d38 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -354,9 +354,8 @@ void CompositionalMultiphaseBase::updateComponentFraction( ObjectManagerBase & d GEOSX_MARK_FUNCTION; ComponentFractionKernelFactory:: - createAndLaunch< parallelDevicePolicy<> > - ( m_numComponents, - dataGroup ); + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + dataGroup ); } @@ -367,14 +366,11 @@ void CompositionalMultiphaseBase::updatePhaseVolumeFraction( ObjectManagerBase & MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, m_fluidModelNames[targetIndex] ); - bool const isIsothermal = true; PhaseVolumeFractionKernelFactory:: - createAndLaunch< parallelDevicePolicy<> > - ( isIsothermal, - m_numComponents, - m_numPhases, - dataGroup, - fluid ); + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid ); } @@ -948,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 348646e15ce..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.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 CompositionalMultiphaseBaseKernels.cpp - */ - -#include "CompositionalMultiphaseBaseKernels.hpp" - -namespace geosx -{ - -namespace CompositionalMultiphaseBaseKernels -{} // namespace CompositionalMultiphaseBaseKernels - -} // namespace geosx diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp index 4aa418c427a..84105c61bb1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp @@ -60,33 +60,13 @@ struct NoOpFunc * @tparam NUM_COMP number of fluid components * @brief Define the base interface for the property update kernels */ -template< localIndex NUM_COMP > +template< integer NUM_COMP > class PropertyKernelBase { public: /// Compile time value for the number of components - static constexpr localIndex numComp = NUM_COMP; - - /** - * @brief Constructor - * @param[in] subRegion the element subregion - * @param[in] fluid the fluid model - */ - PropertyKernelBase( dataRepository::Group & subRegion ) - { GEOSX_UNUSED_VAR( subRegion ); } - - /** - * @brief Compute the property 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 && kernelOp = NoOpFunc{} ) const - { GEOSX_UNUSED_VAR( ei, kernelOp ); } + static constexpr integer numComp = NUM_COMP; /** * @brief Performs the kernel launch @@ -100,8 +80,6 @@ class PropertyKernelBase launch( localIndex const numElems, KERNEL_TYPE const & kernelComponent ) { - GEOSX_MARK_FUNCTION; - forAll< POLICY >( numElems, [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { kernelComponent.compute( ei ); @@ -120,8 +98,6 @@ class PropertyKernelBase launch( SortedArrayView< localIndex const > const & targetSet, KERNEL_TYPE const & kernelComponent ) { - GEOSX_MARK_FUNCTION; - forAll< POLICY >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const i ) { localIndex const ei = targetSet[ i ]; @@ -131,6 +107,34 @@ class PropertyKernelBase }; +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 ********************************/ /** @@ -138,7 +142,7 @@ class PropertyKernelBase * @tparam NUM_COMP number of fluid components * @brief Define the interface for the update kernel in charge of computing the phase volume fractions */ -template< localIndex NUM_COMP > +template< integer NUM_COMP > class ComponentFractionKernel : public PropertyKernelBase< NUM_COMP > { public: @@ -152,7 +156,7 @@ class ComponentFractionKernel : public PropertyKernelBase< NUM_COMP > * @param[in] fluid the fluid model */ ComponentFractionKernel( ObjectManagerBase & subRegion ) - : Base( subRegion ), + : Base(), m_compDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::globalCompDensity >() ), m_dCompDens( subRegion.getExtrinsicData< extrinsicMeshData::flow::deltaGlobalCompDensity >() ), m_compFrac( subRegion.getExtrinsicData< extrinsicMeshData::flow::globalCompFraction >() ), @@ -213,34 +217,6 @@ class ComponentFractionKernel : public PropertyKernelBase< NUM_COMP > }; -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 ComponentFractionKernelFactory */ @@ -257,14 +233,14 @@ class ComponentFractionKernelFactory */ template< typename POLICY > static void - createAndLaunch( localIndex const numComp, + createAndLaunch( integer const numComp, ObjectManagerBase & subRegion ) { internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); ComponentFractionKernel< NUM_COMP > kernel( subRegion ); - ComponentFractionKernel< NUM_COMP >::template launch< POLICY, ComponentFractionKernel< NUM_COMP > >( subRegion.size(), kernel ); + ComponentFractionKernel< NUM_COMP >::template launch< POLICY >( subRegion.size(), kernel ); } ); } @@ -278,7 +254,7 @@ class ComponentFractionKernelFactory * @tparam NUM_PHASE number of fluid phases * @brief Define the interface for the property kernel in charge of computing the phase volume fractions */ -template< localIndex NUM_COMP, localIndex NUM_PHASE > +template< integer NUM_COMP, integer NUM_PHASE > class PhaseVolumeFractionKernel : public PropertyKernelBase< NUM_COMP > { public: @@ -287,7 +263,7 @@ class PhaseVolumeFractionKernel : public PropertyKernelBase< NUM_COMP > using Base::numComp; /// Compile time value for the number of phases - static constexpr localIndex numPhase = NUM_PHASE; + static constexpr integer numPhase = NUM_PHASE; /** * @brief Constructor @@ -296,7 +272,7 @@ class PhaseVolumeFractionKernel : public PropertyKernelBase< NUM_COMP > */ PhaseVolumeFractionKernel( ObjectManagerBase & subRegion, MultiFluidBase const & fluid ) - : Base( subRegion ), + : 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 >() ), @@ -433,7 +409,6 @@ class PhaseVolumeFractionKernelFactory /** * @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] numComp the number of fluid components * @param[in] numPhase the number of fluid phases * @param[in] subRegion the element subregion @@ -441,31 +416,27 @@ class PhaseVolumeFractionKernelFactory */ template< typename POLICY > static void - createAndLaunch( bool const isIsothermal, - localIndex const numComp, - localIndex const numPhase, + createAndLaunch( integer const numComp, + integer const numPhase, ObjectManagerBase & subRegion, MultiFluidBase const & fluid ) { - if( !isIsothermal ) - { GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); } - if( numPhase == 2 ) { internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); PhaseVolumeFractionKernel< NUM_COMP, 2 > kernel( subRegion, fluid ); - PhaseVolumeFractionKernel< NUM_COMP, 2 >::template launch< POLICY, PhaseVolumeFractionKernel< NUM_COMP, 2 > >( subRegion.size(), kernel ); + PhaseVolumeFractionKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); } ); } else if( numPhase == 3 ) { internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); PhaseVolumeFractionKernel< NUM_COMP, 3 > kernel( subRegion, fluid ); - PhaseVolumeFractionKernel< NUM_COMP, 3 >::template launch< POLICY, PhaseVolumeFractionKernel< NUM_COMP, 3 > >( subRegion.size(), kernel ); + PhaseVolumeFractionKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); } ); } } @@ -628,19 +599,19 @@ struct CapillaryPressureUpdateKernel * @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 @@ -946,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; @@ -1003,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 @@ -1016,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, @@ -1029,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 ); } ); } @@ -1064,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, @@ -1081,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; @@ -1131,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 ); @@ -1140,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; @@ -1171,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, @@ -1203,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 ); } @@ -1229,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 @@ -1242,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 ) { @@ -1273,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]; } @@ -1295,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, @@ -1326,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 ); } @@ -1356,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]; } @@ -1461,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 ) { @@ -1470,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 2deb1d6a25c..3c607a9767f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -524,15 +524,12 @@ void CompositionalMultiphaseFVM::updatePhaseMobility( ObjectManagerBase & dataGr MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, m_fluidModelNames[targetIndex] ); RelativePermeabilityBase const & relperm = getConstitutiveModel< RelativePermeabilityBase >( dataGroup, m_relPermModelNames[targetIndex] ); - bool const isIsothermal = true; PhaseMobilityKernelFactory:: - createAndLaunch< parallelDevicePolicy<> > - ( isIsothermal, - m_numComponents, - m_numPhases, - dataGroup, - fluid, - relperm ); + 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 1e0900af38c..46339fed0c4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp @@ -35,11 +35,11 @@ namespace CompositionalMultiphaseFVMKernels /******************************** 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, @@ -69,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{}; @@ -120,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]; } @@ -139,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; } @@ -148,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]; } @@ -162,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]; } @@ -179,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; } @@ -197,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; } @@ -213,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]; } @@ -235,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; } @@ -246,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; } @@ -257,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; @@ -266,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; } @@ -278,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]; } @@ -288,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]; @@ -299,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]; @@ -309,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, @@ -409,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; } @@ -430,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, @@ -446,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, \ @@ -504,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, @@ -527,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 @@ -587,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] @@ -597,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, @@ -660,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, \ @@ -704,7 +704,7 @@ INST_CFLFluxKernel( 5, FaceElementToCellStencilWrapper ); /******************************** CFLKernel ********************************/ -template< localIndex NP > +template< integer NP > GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE void @@ -761,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]; @@ -789,7 +789,7 @@ CFLKernel:: } -template< localIndex NC > +template< integer NC > GEOSX_HOST_DEVICE GEOSX_FORCE_INLINE void @@ -803,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 ) { @@ -817,7 +817,7 @@ CFLKernel:: } } -template< localIndex NC, localIndex NP > +template< integer NC, integer NP > void CFLKernel:: launch( localIndex const size, @@ -904,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, @@ -985,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, @@ -1027,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]{}; @@ -1073,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; } @@ -1092,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, @@ -1107,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 af52c4efa6f..46beab765fd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp @@ -46,7 +46,7 @@ using namespace constitutive; * @tparam NUM_PHASE number of fluid phases * @brief Define the interface for the property kernel in charge of computing the phase mobilities */ -template< localIndex NUM_COMP, localIndex NUM_PHASE > +template< integer NUM_COMP, integer NUM_PHASE > class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { public: @@ -55,7 +55,7 @@ class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyK using Base::numComp; /// Compile time value for the number of phases - static constexpr localIndex numPhase = NUM_PHASE; + static constexpr integer numPhase = NUM_PHASE; /** * @brief Constructor @@ -66,7 +66,7 @@ class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyK PhaseMobilityKernel( ObjectManagerBase & subRegion, MultiFluidBase const & fluid, RelativePermeabilityBase const & relperm ) - : Base( subRegion ), + : 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 >() ), @@ -219,7 +219,6 @@ class PhaseMobilityKernelFactory /** * @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] numComp the number of fluid components * @param[in] numPhase the number of fluid phases * @param[in] subRegion the element subregion @@ -228,32 +227,28 @@ class PhaseMobilityKernelFactory */ template< typename POLICY > static void - createAndLaunch( bool const isIsothermal, - localIndex const numComp, - localIndex const numPhase, + createAndLaunch( integer const numComp, + integer const numPhase, ObjectManagerBase & subRegion, MultiFluidBase const & fluid, RelativePermeabilityBase const & relperm ) { - if( !isIsothermal ) - { GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); } - if( numPhase == 2 ) { CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); - PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY, PhaseMobilityKernel< NUM_COMP, 2 > >( subRegion.size(), kernel ); + PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); } ); } else if( numPhase == 3 ) { CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); - PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY, PhaseMobilityKernel< NUM_COMP, 3 > >( subRegion.size(), kernel ); + PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); } ); } } @@ -276,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, @@ -309,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, @@ -361,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, @@ -382,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, @@ -412,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 @@ -424,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 @@ -434,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, @@ -471,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, @@ -495,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 dd124197f41..48c717b0a3e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -895,16 +895,13 @@ void CompositionalMultiphaseHybridFVM::updatePhaseMobility( ObjectManagerBase & MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, m_fluidModelNames[targetIndex] ); RelativePermeabilityBase const & relperm = getConstitutiveModel< RelativePermeabilityBase >( dataGroup, m_relPermModelNames[targetIndex] ); - bool const isIsothermal = true; CompositionalMultiphaseHybridFVMKernels:: PhaseMobilityKernelFactory:: - createAndLaunch< parallelDevicePolicy<> > - ( isIsothermal, - m_numComponents, - m_numPhases, - dataGroup, - fluid, - relperm ); + 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 37b4d16da98..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, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp index ad9b42bd889..9f03002f050 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp @@ -80,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 ], @@ -126,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 ], @@ -167,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 ], @@ -193,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 ], @@ -237,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 ], @@ -286,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, @@ -338,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 ], @@ -393,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, @@ -431,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, @@ -459,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, @@ -530,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, @@ -625,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, @@ -674,7 +674,7 @@ struct FluxKernel * @tparam NUM_PHASE number of fluid phases * @brief Define the interface for the property kernel in charge of computing the phase mobilities */ -template< localIndex NUM_COMP, localIndex NUM_PHASE > +template< integer NUM_COMP, integer NUM_PHASE > class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { public: @@ -683,7 +683,7 @@ class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyK using Base::numComp; /// Compile time value for the number of phases - static constexpr localIndex numPhase = NUM_PHASE; + static constexpr integer numPhase = NUM_PHASE; /** * @brief Constructor @@ -694,7 +694,7 @@ class PhaseMobilityKernel : public CompositionalMultiphaseBaseKernels::PropertyK PhaseMobilityKernel( ObjectManagerBase & subRegion, MultiFluidBase const & fluid, RelativePermeabilityBase const & relperm ) - : Base( subRegion ), + : 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 >() ), @@ -831,7 +831,6 @@ class PhaseMobilityKernelFactory /** * @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] numComp the number of fluid components * @param[in] numPhase the number of fluid phases * @param[in] subRegion the element subregion @@ -840,32 +839,28 @@ class PhaseMobilityKernelFactory */ template< typename POLICY > static void - createAndLaunch( bool const isIsothermal, - localIndex const numComp, - localIndex const numPhase, + createAndLaunch( integer const numComp, + integer const numPhase, ObjectManagerBase & subRegion, MultiFluidBase const & fluid, RelativePermeabilityBase const & relperm ) { - if( !isIsothermal ) - { GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); } - if( numPhase == 2 ) { CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); - PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY, PhaseMobilityKernel< NUM_COMP, 2 > >( subRegion.size(), kernel ); + PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); } ); } else if( numPhase == 3 ) { CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); - PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY, PhaseMobilityKernel< NUM_COMP, 3 > >( subRegion.size(), kernel ); + PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); } ); } } @@ -915,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]; } @@ -976,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, @@ -1009,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]; @@ -1051,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 f53e92b75df..07161ab32c2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -445,9 +445,8 @@ void CompositionalMultiphaseWell::updateComponentFraction( WellElementSubRegion CompositionalMultiphaseBaseKernels:: ComponentFractionKernelFactory:: - createAndLaunch< parallelDevicePolicy<> > - ( m_numComponents, - subRegion ); + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + subRegion ); } @@ -747,15 +746,12 @@ void CompositionalMultiphaseWell::updatePhaseVolumeFraction( WellElementSubRegio MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, m_fluidModelNames[targetIndex] ); - bool const isIsothermal = true; CompositionalMultiphaseBaseKernels:: PhaseVolumeFractionKernelFactory:: - createAndLaunch< parallelDevicePolicy<> > - ( isIsothermal, - m_numComponents, - m_numPhases, - subRegion, - fluid ); + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + subRegion, + fluid ); } @@ -765,14 +761,11 @@ void CompositionalMultiphaseWell::updateTotalMassDensity( WellElementSubRegion & MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, m_fluidModelNames[targetIndex] ); - bool const isIsothermal = true; TotalMassDensityKernelFactory:: - createAndLaunch< parallelDevicePolicy<> > - ( isIsothermal, - m_numComponents, - m_numPhases, - subRegion, - fluid ); + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + subRegion, + fluid ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index d97f03ac3d9..c43de914887 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -446,7 +446,7 @@ struct RateInitializationKernel * @tparam NUM_PHASE number of fluid phases * @brief Define the interface for the property kernel in charge of computing the total mass density */ -template< localIndex NUM_COMP, localIndex NUM_PHASE > +template< integer NUM_COMP, integer NUM_PHASE > class TotalMassDensityKernel : public CompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { public: @@ -458,7 +458,7 @@ class TotalMassDensityKernel : public CompositionalMultiphaseBaseKernels::Proper using Base::numComp; /// Compile time value for the number of phases - static constexpr localIndex numPhase = NUM_PHASE; + static constexpr integer numPhase = NUM_PHASE; /** * @brief Constructor @@ -467,7 +467,7 @@ class TotalMassDensityKernel : public CompositionalMultiphaseBaseKernels::Proper */ TotalMassDensityKernel( dataRepository::Group & subRegion, MultiFluidBase const & fluid ) - : Base( subRegion ), + : 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() ) ), @@ -562,7 +562,6 @@ class TotalMassDensityKernelFactory /** * @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] numComp the number of fluid components * @param[in] numPhase the number of fluid phases * @param[in] subRegion the element subregion @@ -570,31 +569,27 @@ class TotalMassDensityKernelFactory */ template< typename POLICY > static void - createAndLaunch( bool const isIsothermal, - localIndex const numComp, - localIndex const numPhase, + createAndLaunch( integer const numComp, + integer const numPhase, dataRepository::Group & subRegion, MultiFluidBase const & fluid ) { - if( !isIsothermal ) - { GEOSX_ERROR( "CompositionalMultiphaseBase: Thermal simulation is not supported yet: " ); } - if( numPhase == 2 ) { CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); TotalMassDensityKernel< NUM_COMP, 2 > kernel( subRegion, fluid ); - TotalMassDensityKernel< NUM_COMP, 2 >::template launch< POLICY, TotalMassDensityKernel< NUM_COMP, 2 > >( subRegion.size(), kernel ); + TotalMassDensityKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); } ); } else if( numPhase == 3 ) { CompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { - localIndex constexpr NUM_COMP = NC(); + integer constexpr NUM_COMP = NC(); TotalMassDensityKernel< NUM_COMP, 3 > kernel( subRegion, fluid ); - TotalMassDensityKernel< NUM_COMP, 3 >::template launch< POLICY, TotalMassDensityKernel< NUM_COMP, 3 > >( subRegion.size(), kernel ); + TotalMassDensityKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); } ); } }