diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 6d64562d74a..2aadc509e55 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -256,11 +256,26 @@ void CompositionalMultiphaseBase::postInputInitialization() getWrapperDataContext( viewKeyStruct::minScalingFactorString() ) << ": The minumum scaling factor must be smaller or equal to 1.0" ); - if( m_isThermal && m_useSimpleAccumulation == 1 ) // useSimpleAccumulation is not yet compatible with thermal + if( m_isThermal && m_useSimpleAccumulation ) // useSimpleAccumulation is not yet compatible with thermal { - GEOS_LOG_RANK_0( "'useSimpleAccumulation' is not yet implemented for thermal simulation. Switched to phase sum accumulation." ); + GEOS_LOG_RANK_0( GEOS_FMT( "{}: '{}' is not yet implemented for thermal simulation. Switched to phase sum accumulation.", + getDataContext(), viewKeyStruct::useSimpleAccumulationString() ) ); m_useSimpleAccumulation = 0; } + + if( m_useZFormulation ) + { + if( m_isThermal ) // useZFormulation is not yet compatible with thermal + { + GEOS_ERROR( GEOS_FMT( "{}: '{}' is currently not available for thermal simulations", + getDataContext(), viewKeyStruct::useZFormulationFlagString() ) ); + } + if( m_isJumpStabilized ) // useZFormulation is not yet compatible with pressure stabilization + { + GEOS_ERROR( GEOS_FMT( "{}: pressure stabilization is not yet supported by {}", + getDataContext(), viewKeyStruct::useZFormulationFlagString() ) ); + } + } } void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) @@ -323,10 +338,12 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) if( m_useZFormulation ) + { // Z formulation - component densities and pressure are primary unknowns // (n_c-1) overall compositions + one pressure // Testing: let's have sum_c zc = 1 as explicit equation for now m_numDofPerCell = m_numComponents + 1; + } else { // default formulation - component densities and pressure are primary unknowns @@ -406,16 +423,14 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< pressureScalingFactor >( getName() ); subRegion.registerField< temperatureScalingFactor >( getName() ); - subRegion.registerField< globalCompDensityScalingFactor >( getName() ); - subRegion.registerField< globalCompFractionScalingFactor >( getName() ); // The resizing of the arrays needs to happen here, before the call to initializePreSubGroups, // to make sure that the dimensions are properly set before the timeHistoryOutput starts its initialization. + subRegion.registerField< globalCompFraction >( getName() ). + setDimLabels( 1, fluid.componentNames() ). + reference().resizeDimension< 1 >( m_numComponents ); if( m_useZFormulation ) { - subRegion.registerField< globalCompFraction >( getName() ). - setDimLabels( 1, fluid.componentNames() ). - reference().resizeDimension< 1 >( m_numComponents ); subRegion.registerField< globalCompFraction_n >( getName() ). setDimLabels( 1, fluid.componentNames() ). reference().resizeDimension< 1 >( m_numComponents ); @@ -425,6 +440,7 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) setDimLabels( 1, fluid.componentNames() ). reference().resizeDimension< 1 >( m_numComponents ); } + subRegion.registerField< globalCompFractionScalingFactor >( getName() ); } else { @@ -439,9 +455,7 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) setDimLabels( 1, fluid.componentNames() ). reference().resizeDimension< 1 >( m_numComponents ); } - subRegion.registerField< globalCompFraction >( getName() ). - setDimLabels( 1, fluid.componentNames() ). - reference().resizeDimension< 1 >( m_numComponents ); + subRegion.registerField< globalCompDensityScalingFactor >( getName() ); subRegion.registerField< dGlobalCompFraction_dGlobalCompDensity >( getName() ). reference().resizeDimension< 1, 2 >( m_numComponents, m_numComponents ); } @@ -707,43 +721,36 @@ real64 CompositionalMultiphaseBase::updatePhaseVolumeFraction( ObjectManagerBase string const & fluidName = dataGroup.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, fluidName ); - real64 maxDeltaPhaseVolFrac = - m_isThermal ? - thermalCompositionalMultiphaseBaseKernels:: - PhaseVolumeFractionKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dataGroup, - fluid ) -: isothermalCompositionalMultiphaseBaseKernels:: - PhaseVolumeFractionKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dataGroup, - fluid ); - - return maxDeltaPhaseVolFrac; -} - -real64 CompositionalMultiphaseBase::updatePhaseVolumeFractionZFormulation( ObjectManagerBase & dataGroup ) const -{ - GEOS_MARK_FUNCTION; - - string const & fluidName = dataGroup.getReference< string >( viewKeyStruct::fluidNamesString() ); - MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, fluidName ); - - // For now: isothermal only - GEOS_ERROR_IF( m_isThermal, GEOS_FMT( "{}: Z Formulation is currently not available for thermal simulations", getDataContext() ) ); - - real64 maxDeltaPhaseVolFrac = - isothermalCompositionalMultiphaseBaseKernels:: - PhaseVolumeFractionZFormulationKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dataGroup, - fluid ); - - return maxDeltaPhaseVolFrac; + if( m_isThermal ) + { + return thermalCompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid ); + } + else + { + if( m_useZFormulation ) + { + return isothermalCompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionZFormulationKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid ); + } + else + { + return isothermalCompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid ); + } + } } void CompositionalMultiphaseBase::updateFluidModel( ObjectManagerBase & dataGroup ) const @@ -926,9 +933,8 @@ real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & sub if( m_useZFormulation ) { // For p, z_c as the primary unknowns - updateFluidModel( subRegion ); // rho_T is now a function of p, z_c from volume balance + updateFluidModel( subRegion ); // rho_T is now a function of p, z_c from volume balance updateCompAmountZFormulation( subRegion ); - maxDeltaPhaseVolFrac = updatePhaseVolumeFractionZFormulation( subRegion ); } else { @@ -936,10 +942,8 @@ real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & sub updateGlobalComponentFraction( subRegion ); updateFluidModel( subRegion ); updateCompAmount( subRegion ); - maxDeltaPhaseVolFrac = updatePhaseVolumeFraction( subRegion ); - } - + maxDeltaPhaseVolFrac = updatePhaseVolumeFraction( subRegion ); updateRelPermModel( subRegion ); updatePhaseMobility( subRegion ); updateCapPressureModel( subRegion ); @@ -1011,13 +1015,12 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, if( m_useZFormulation ) { updateCompAmountZFormulation( subRegion ); - updatePhaseVolumeFractionZFormulation( subRegion ); } else { updateCompAmount( subRegion ); - updatePhaseVolumeFraction( subRegion ); } + updatePhaseVolumeFraction( subRegion ); // Update the constitutive models that only depend on // - the primary variables @@ -1036,9 +1039,9 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, // - This step depends on phaseVolFraction RelativePermeabilityBase & relPerm = getConstitutiveModel< RelativePermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::relPermNamesString() ) ); - relPerm.saveConvergedPhaseVolFractionState( phaseVolFrac ); // this needs to happen before calling updateRelPermModel + relPerm.saveConvergedPhaseVolFractionState( phaseVolFrac ); // this needs to happen before calling updateRelPermModel updateRelPermModel( subRegion ); - relPerm.saveConvergedState(); // this needs to happen after calling updateRelPermModel + relPerm.saveConvergedState(); // this needs to happen after calling updateRelPermModel string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); @@ -1068,7 +1071,7 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, CapillaryPressureBase const & capPressure = getConstitutiveModel< CapillaryPressureBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::capPressureNamesString() ) ); - capPressure.initializeRockState( porosity, permeability ); // this needs to happen before calling updateCapPressureModel + capPressure.initializeRockState( porosity, permeability ); // this needs to happen before calling updateCapPressureModel updateCapPressureModel( subRegion ); } @@ -1200,7 +1203,7 @@ void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition real64 const equilTolerance = fs.getEquilibrationTolerance(); real64 const datumElevation = fs.getDatumElevation(); real64 const datumPressure = fs.getDatumPressure(); - string const initPhaseName = fs.getInitPhaseName(); // will go away when GOC/WOC are implemented + string const initPhaseName = fs.getInitPhaseName(); // will go away when GOC/WOC are implemented localIndex const equilIndex = equilNameToEquilId.at( fs.getName() ); real64 const minElevation = LvArray::math::min( globalMinElevation[equilIndex], datumElevation ); @@ -1208,7 +1211,7 @@ void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition real64 const elevationIncrement = LvArray::math::min( fs.getElevationIncrement(), maxElevation - minElevation ); localIndex const numPointsInTable = ( elevationIncrement > 0 ) ? std::ceil( (maxElevation - minElevation) / elevationIncrement ) + 1 : 1; - real64 const eps = 0.1 * (maxElevation - minElevation); // we add a small buffer to only log in the pathological cases + real64 const eps = 0.1 * (maxElevation - minElevation); // we add a small buffer to only log in the pathological cases GEOS_LOG_RANK_0_IF( ( (datumElevation > globalMaxElevation[equilIndex]+eps) || (datumElevation < globalMinElevation[equilIndex]-eps) ), getCatalogName() << " " << getDataContext() << ": By looking at the elevation of the cell centers in this model, GEOS found that " << @@ -1247,7 +1250,7 @@ void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition auto itRegionFilter = regionFilter.find( region.getName() ); if( itRegionFilter == regionFilter.end() ) { - return; // the region is not in target, there is nothing to do + return; // the region is not in target, there is nothing to do } string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); @@ -1577,27 +1580,6 @@ void CompositionalMultiphaseBase::assembleZFormulationAccumulation( DomainPartit MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); CoupledSolidBase const & solid = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - // TODO: add the thermal case - /* - if( m_isThermal ) - { - thermalCompositionalMultiphaseBaseKernels:: - AccumulationKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - m_useTotalMassEquation, - dofKey, - subRegion, - fluid, - solid, - localMatrix, - localRhs ); - } - else - { - } - */ // isothermal for now isothermalCompositionalMultiphaseBaseKernels:: AccumulationZFormulationKernelFactory:: @@ -1791,7 +1773,7 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time, return; } - real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor here! + real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor here! massProd += rhsValue; if( useTotalMassEquation > 0 ) { @@ -1800,7 +1782,7 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time, localRhs[totalMassBalanceRow] += rhsValue; if( fluidComponentId < numFluidComponents - 1 ) { - globalIndex const compMassBalanceRow = totalMassBalanceRow + fluidComponentId + 1; // component mass bal equations are shifted + globalIndex const compMassBalanceRow = totalMassBalanceRow + fluidComponentId + 1; // component mass bal equations are shifted localRhs[compMassBalanceRow] += rhsValue; } } @@ -1928,7 +1910,7 @@ bool CompositionalMultiphaseBase::validateDirichletBC( DomainPartition & domain, bcConsistent = false; GEOS_WARNING( BCMessage::invalidComponentIndex( comp, fs.getName(), fields::flow::globalCompFraction::key() ) ); - return; // can't check next part with invalid component id + return; // can't check next part with invalid component id } ComponentMask< MAX_NC > & compMask = subRegionSetMap[setName]; @@ -2244,7 +2226,7 @@ void CompositionalMultiphaseBase::keepVariablesConstantDuringInitStep( real64 co rankOffset, localMatrix, rhsValue, - pres[ei], // freeze the current pressure value + pres[ei], // freeze the current pressure value pres[ei] ); localRhs[localRow] = rhsValue; @@ -2255,7 +2237,7 @@ void CompositionalMultiphaseBase::keepVariablesConstantDuringInitStep( real64 co rankOffset, localMatrix, rhsValue, - temp[ei], // freeze the current temperature value + temp[ei], // freeze the current temperature value temp[ei] ); localRhs[localRow + numComp + 1] = rhsValue; } @@ -2267,7 +2249,7 @@ void CompositionalMultiphaseBase::keepVariablesConstantDuringInitStep( real64 co rankOffset, localMatrix, rhsValue, - compDens[ei][ic], // freeze the current component density values + compDens[ei][ic], // freeze the current component density values compDens[ei][ic] ); localRhs[localRow + ic + 1] = rhsValue; } @@ -2869,11 +2851,11 @@ void CompositionalMultiphaseBase::implicitStepComplete( real64 const & time, CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); if( m_keepVariablesConstantDuringInitStep ) { - porousMaterial.ignoreConvergedState(); // newPorosity <- porosity_n + porousMaterial.ignoreConvergedState(); // newPorosity <- porosity_n } else { - porousMaterial.saveConvergedState(); // porosity_n <- porosity + porousMaterial.saveConvergedState(); // porosity_n <- porosity } // Step 4: save converged state for the relperm model to handle hysteresis @@ -3017,7 +2999,7 @@ void CompositionalMultiphaseBase::saveSequentialIterationState( DomainPartition } ); } ); - m_sequentialCompDensChange = MpiWrapper::max( maxCompDensChange ); // store to be later used for convergence check + m_sequentialCompDensChange = MpiWrapper::max( maxCompDensChange ); // store to be later used for convergence check } void CompositionalMultiphaseBase::updateState( DomainPartition & domain ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index b6fd40cb6e0..0503eaad4c1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -116,12 +116,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ real64 updatePhaseVolumeFraction( ObjectManagerBase & dataGroup ) const; - /** - * @brief Recompute phase volume fractions (saturations) from constitutive and primary variables (pressure and global component fractions) - * @param dataGroup the group storing the required fields - */ - real64 updatePhaseVolumeFractionZFormulation( ObjectManagerBase & dataGroup ) const; - /** * @brief Update all relevant fluid models using current values of pressure and composition * @param dataGroup the group storing the required fields diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 694008ccd14..84fed924935 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -329,6 +329,7 @@ void CompositionalMultiphaseFVM::assembleZFormulationFluxTerms( real64 const dt, elemDofKey, m_hasCapPressure, m_useTotalMassEquation, + m_useNewGravity, fluxApprox.upwindingParams(), getName(), mesh.getElemManager(), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/FluxComputeZFormulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/FluxComputeZFormulationKernel.hpp index 0e8805e8d98..0286d1e6746 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/FluxComputeZFormulationKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/FluxComputeZFormulationKernel.hpp @@ -272,12 +272,11 @@ class FluxComputeZFormulationKernel : public FluxComputeKernelBase localIndex k_up = -1; - - isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFluxZFormulation::compute< numComp, numFluxSupportPoints > ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), + m_kernelFlags.isSet( KernelFlags::NewGravity ), seri, sesri, sei, trans, dTrans_dPres, @@ -472,6 +471,7 @@ class FluxComputeZFormulationKernelFactory string const & dofKey, integer const hasCapPressure, integer const useTotalMassEquation, + integer const useNewGravity, UpwindingParameters upwindingParams, string const & solverName, ElementRegionManager const & elemManager, @@ -494,16 +494,16 @@ class FluxComputeZFormulationKernelFactory kernelFlags.set( KernelFlags::CapPressure ); if( useTotalMassEquation ) kernelFlags.set( KernelFlags::TotalMassEquation ); + if( useNewGravity ) + kernelFlags.set( KernelFlags::NewGravity ); if( upwindingParams.upwindingScheme == UpwindingScheme::C1PPU && isothermalCompositionalMultiphaseFVMKernelUtilities::epsC1PPU > 0 ) { - kernelFlags.set( KernelFlags::C1PPU ); - //GEOS_FMT("CompositionalMultiphaseBase {}: Z Formulation is currently not available for C1PPU ", getDataContext() ); + GEOS_ERROR( "Z Formulation is currently not available for C1PPU" ); } else if( upwindingParams.upwindingScheme == UpwindingScheme::IHU ) { - kernelFlags.set( KernelFlags::IHU ); - //GEOS_FMT("CompositionalMultiphaseBase {}: Z Formulation is currently not available for IHU ", getDataContext() ); + GEOS_ERROR( "Z Formulation is currently not available for IHU" ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PPUPhaseFluxZFormulation.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PPUPhaseFluxZFormulation.hpp index d30d14406a5..cfb4ba9b12f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PPUPhaseFluxZFormulation.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PPUPhaseFluxZFormulation.hpp @@ -74,6 +74,7 @@ struct PPUPhaseFluxZFormulation compute( integer const numPhase, integer const ip, integer const hasCapPressure, + integer const useNewGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -104,7 +105,7 @@ struct PPUPhaseFluxZFormulation real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; real64 dGravHead_dC[numFluxSupportPoints][numComp]{}; - PotGradZFormulation::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, seri, sesri, sei, trans, dTrans_dPres, pres, + PotGradZFormulation::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, useNewGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, phaseVolFrac, dPhaseVolFrac, phaseMassDens, dPhaseMassDens, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PotGradZFormulation.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PotGradZFormulation.hpp index 6aa9be4d161..4cd4a758ea2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PotGradZFormulation.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PotGradZFormulation.hpp @@ -46,6 +46,7 @@ struct PotGradZFormulation compute ( integer const numPhase, integer const ip, integer const hasCapPressure, + integer const useNewGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -86,45 +87,7 @@ struct PotGradZFormulation real64 gravHead = 0.0; real64 dCapPressure_dC[numComp]{}; - // calculate quantities on primary connected cells - integer denom = 0; - for( integer i = 0; i < numFluxSupportPoints; ++i ) - { - localIndex const er = seri[i]; - localIndex const esr = sesri[i]; - localIndex const ei = sei[i]; - - bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( !phaseExists ) - { - continue; - } - - // density - real64 const density = phaseMassDens[er][esr][ei][0][ip]; - real64 const dDens_dP = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP]; - - // average density and derivatives - densMean += density; - dDensMean_dP[i] = dDens_dP; - for( integer jc = 0; jc < numComp; ++jc ) - { - dDensMean_dC[i][jc] = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dC+jc]; - } - denom++; - } - if( denom > 1 ) - { - densMean /= denom; - for( integer i = 0; i < numFluxSupportPoints; ++i ) - { - dDensMean_dP[i] /= denom; - for( integer jc = 0; jc < numComp; ++jc ) - { - dDensMean_dC[i][jc] /= denom; - } - } - } + calculateMeanDensity( useNewGravity, ip, seri, sesri, sei, phaseVolFrac, phaseMassDens, dPhaseMassDens, densMean, dDensMean_dP, dDensMean_dC ); /// compute the TPFA potential difference for( integer i = 0; i < numFluxSupportPoints; i++ ) @@ -190,6 +153,60 @@ struct PotGradZFormulation } + template< integer numComp, integer numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + calculateMeanDensity( integer const useNewGravity, + integer const ip, + localIndex const ( &seri )[numFluxSupportPoints], + localIndex const ( &sesri )[numFluxSupportPoints], + localIndex const ( &sei )[numFluxSupportPoints], + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, + ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, + real64 & densMean, real64 ( & dDensMean_dP)[numFluxSupportPoints], real64 ( & dDensMean_dC )[numFluxSupportPoints][numComp] ) + { + // calculate quantities on primary connected cells + integer denom = 0; + for( integer i = 0; i < numFluxSupportPoints; ++i ) + { + localIndex const er = seri[i]; + localIndex const esr = sesri[i]; + localIndex const ei = sei[i]; + + bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); + if( useNewGravity && !phaseExists ) + { + continue; + } + + // density + real64 const density = phaseMassDens[er][esr][ei][0][ip]; + real64 const dDens_dP = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP]; + + // average density and derivatives + densMean += density; + dDensMean_dP[i] = dDens_dP; + for( integer jc = 0; jc < numComp; ++jc ) + { + dDensMean_dC[i][jc] = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dC+jc]; + } + denom++; + } + if( denom > 1 ) + { + densMean /= denom; + for( integer i = 0; i < numFluxSupportPoints; ++i ) + { + dDensMean_dP[i] /= denom; + for( integer jc = 0; jc < numComp; ++jc ) + { + dDensMean_dC[i][jc] /= denom; + } + } + } + } + }; } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities