diff --git a/inputFiles/compositionalMultiphaseFlow/co2_flux_3d.xml b/inputFiles/compositionalMultiphaseFlow/co2_flux_3d.xml index 714ac6337af..f3334b3e9d4 100644 --- a/inputFiles/compositionalMultiphaseFlow/co2_flux_3d.xml +++ b/inputFiles/compositionalMultiphaseFlow/co2_flux_3d.xml @@ -8,6 +8,8 @@ discretization="fluidTPFA" temperature="368.15" useMass="1" + useTotalMassEquation="0" + useSimpleAccumulation="1" targetRegions="{ region }"> +struct BitFlags +{ + GEOS_HOST_DEVICE + void set( FLAGS_ENUM flag ) + { + m_FlagValue |= (integer)flag; + } + + GEOS_HOST_DEVICE + bool isSet( FLAGS_ENUM flag ) const + { + return (m_FlagValue & (integer)flag) == (integer)flag; + } + +private: + + using FlagsType = std::underlying_type_t< FLAGS_ENUM >; + FlagsType m_FlagValue{}; + +}; } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 5b3e6998d0b..a1abd0b004d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -67,7 +67,10 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, m_hasDispersion( 0 ), m_keepFlowVariablesConstantDuringInitStep( 0 ), m_minScalingFactor( 0.01 ), - m_allowCompDensChopping( 1 ) + m_allowCompDensChopping( 1 ), + m_useTotalMassEquation( 1 ), + m_useSimpleAccumulation( 0 ), + m_minCompDens( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision ) { //START_SPHINX_INCLUDE_00 this->registerWrapper( viewKeyStruct::inputTemperatureString(), &m_inputTemperature ). @@ -122,6 +125,24 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, setApplyDefaultValue( 1 ). setDescription( "Flag indicating whether local (cell-wise) chopping of negative compositions is allowed" ); + this->registerWrapper( viewKeyStruct::useTotalMassEquationString(), &m_useTotalMassEquation ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1 ). + setDescription( "Flag indicating whether total mass equation is used" ); + + this->registerWrapper( viewKeyStruct::useSimpleAccumulationString(), &m_useSimpleAccumulation ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "Flag indicating whether simple accumulation form is used" ); + + this->registerWrapper( viewKeyStruct::minCompDensString(), &m_minCompDens ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision ). + setDescription( "Minimum allowed global component density" ); + } void CompositionalMultiphaseBase::postProcessInput() @@ -222,6 +243,7 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) m_numPhases = referenceFluid.numFluidPhases(); m_numComponents = referenceFluid.numFluidComponents(); } + // n_c components + one pressure ( + one temperature if needed ) m_numDofPerCell = m_isThermal ? m_numComponents + 2 : m_numComponents + 1; @@ -1300,6 +1322,8 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), + m_useTotalMassEquation, + m_useSimpleAccumulation, dofKey, subRegion, fluid, @@ -1314,6 +1338,8 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), + m_useTotalMassEquation, + m_useSimpleAccumulation, dofKey, subRegion, fluid, @@ -1469,12 +1495,14 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time, integer const fluidComponentId = fs.getComponent(); integer const numFluidComponents = m_numComponents; + integer const useTotalMassEquation = m_useTotalMassEquation; forAll< parallelDevicePolicy<> >( targetSet.size(), [sizeScalingFactor, targetSet, rankOffset, ghostRank, fluidComponentId, numFluidComponents, + useTotalMassEquation, dofNumber, rhsContributionArrayView, localRhs] GEOS_HOST_DEVICE ( localIndex const a ) @@ -1486,15 +1514,23 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time, return; } - // for all "fluid components", we add the value to the total mass balance equation - globalIndex const totalMassBalanceRow = dofNumber[ei] - rankOffset; - localRhs[totalMassBalanceRow] += rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor - // here - - // for all "fluid components" except the last one, we add the value to the component mass balance equation (shifted appropriately) - if( fluidComponentId < numFluidComponents - 1 ) + if( useTotalMassEquation > 0 ) { - globalIndex const compMassBalanceRow = totalMassBalanceRow + fluidComponentId + 1; // component mass bal equations are shifted + // for all "fluid components", we add the value to the total mass balance equation + globalIndex const totalMassBalanceRow = dofNumber[ei] - rankOffset; + localRhs[totalMassBalanceRow] += rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the + // sizeScalingFactor + // here + if( fluidComponentId < numFluidComponents - 1 ) + { + globalIndex const compMassBalanceRow = totalMassBalanceRow + fluidComponentId + 1; // component mass bal equations are shifted + localRhs[compMassBalanceRow] += rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the + // sizeScalingFactor here + } + } + else + { + globalIndex const compMassBalanceRow = dofNumber[ei] - rankOffset + fluidComponentId; localRhs[compMassBalanceRow] += rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the // sizeScalingFactor here } @@ -1919,9 +1955,9 @@ void CompositionalMultiphaseBase::chopNegativeDensities( DomainPartition & domai { for( integer ic = 0; ic < numComp; ++ic ) { - if( compDens[ei][ic] < minDensForDivision ) + if( compDens[ei][ic] < m_minCompDens ) { - compDens[ei][ic] = minDensForDivision; + compDens[ei][ic] = m_minCompDens; } } } @@ -2191,6 +2227,7 @@ void CompositionalMultiphaseBase::saveIterationState( ElementSubRegionBase & sub void CompositionalMultiphaseBase::updateState( DomainPartition & domain ) { + real64 maxDeltaPhaseVolFrac = 0.0; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) @@ -2202,7 +2239,8 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain ) // update porosity, permeability, and solid internal energy updatePorosityAndPermeability( subRegion ); // update all fluid properties - updateFluidState( subRegion ); + real64 const deltaPhaseVolFrac = updateFluidState( subRegion ); + maxDeltaPhaseVolFrac = LvArray::math::max( maxDeltaPhaseVolFrac, deltaPhaseVolFrac ); // for thermal, update solid internal energy if( m_isThermal ) { @@ -2210,6 +2248,8 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain ) } } ); } ); + + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "{}: Max deltaPhaseVolFrac = {}", getName(), maxDeltaPhaseVolFrac ) ); } } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 49e6510f193..a90c23fe6b4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -221,8 +221,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const = 0; - - /**@}*/ struct viewKeyStruct : FlowSolverBase::viewKeyStruct @@ -253,6 +251,9 @@ class CompositionalMultiphaseBase : public FlowSolverBase static constexpr char const * maxRelativePresChangeString() { return "maxRelativePressureChange"; } static constexpr char const * maxRelativeTempChangeString() { return "maxRelativeTemperatureChange"; } static constexpr char const * allowLocalCompDensChoppingString() { return "allowLocalCompDensityChopping"; } + static constexpr char const * useTotalMassEquationString() { return "useTotalMassEquation"; } + static constexpr char const * useSimpleAccumulationString() { return "useSimpleAccumulation"; } + static constexpr char const * minCompDensString() { return "minCompDens"; } }; @@ -357,6 +358,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual void initializePostInitialConditionsPreSubGroups() override; + integer useTotalMassEquation() const { return m_useTotalMassEquation; } + protected: virtual void postProcessInput() override; @@ -445,6 +448,15 @@ class CompositionalMultiphaseBase : public FlowSolverBase /// flag indicating whether local (cell-wise) chopping of negative compositions is allowed integer m_allowCompDensChopping; + /// flag indicating whether total mass equation is used + integer m_useTotalMassEquation; + + /// flag indicating whether simple accumulation form is used + integer m_useSimpleAccumulation; + + /// minimum allowed global component density + real64 m_minCompDens; + /// name of the fluid constitutive model used as a reference for component/phase description string m_referenceFluidModelName; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 68afcffb5bb..227dea7ef2e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -74,8 +74,8 @@ void CompositionalMultiphaseFVM::initializePreSubGroups() CompositionalMultiphaseBase::initializePreSubGroups(); m_linearSolverParameters.get().mgr.strategy = m_isThermal - ? LinearSolverParameters::MGR::StrategyType::thermalCompositionalMultiphaseFVM - : LinearSolverParameters::MGR::StrategyType::compositionalMultiphaseFVM; + ? LinearSolverParameters::MGR::StrategyType::thermalCompositionalMultiphaseFVM + : LinearSolverParameters::MGR::StrategyType::compositionalMultiphaseFVM; DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); @@ -145,6 +145,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, dofManager.rankOffset(), elemDofKey, m_hasCapPressure, + m_useTotalMassEquation, getName(), mesh.getElemManager(), stencilWrapper, @@ -161,6 +162,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, dofManager.rankOffset(), elemDofKey, m_hasCapPressure, + m_useTotalMassEquation, fluxApprox.upwindingParams(), getName(), mesh.getElemManager(), @@ -184,6 +186,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, elemDofKey, m_hasDiffusion, m_hasDispersion, + m_useTotalMassEquation, getName(), mesh.getElemManager(), stencilWrapper, @@ -201,6 +204,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, elemDofKey, m_hasDiffusion, m_hasDispersion, + m_useTotalMassEquation, getName(), mesh.getElemManager(), stencilWrapper, @@ -248,6 +252,7 @@ void CompositionalMultiphaseFVM::assembleStabilizedFluxTerms( real64 const dt, dofManager.rankOffset(), elemDofKey, m_hasCapPressure, + m_useTotalMassEquation, getName(), mesh.getElemManager(), stencilWrapper, @@ -895,6 +900,7 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n, createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), + m_useTotalMassEquation, elemDofKey, getName(), faceManager, @@ -912,6 +918,7 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n, createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), + m_useTotalMassEquation, elemDofKey, getName(), faceManager, @@ -992,6 +999,7 @@ void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, m_numPhases, waterPhaseIndex, allowAllPhasesIntoAquifer, + m_useTotalMassEquation, stencil, dofManager.rankOffset(), elemDofNumber.toNestedViewConst(), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index fab788b2844..a42d5a40b93 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -407,6 +407,7 @@ void CompositionalMultiphaseHybridFVM::assembleFluxTerms( real64 const dt, dofManager.rankOffset(), lengthTolerance, dt, + m_useTotalMassEquation, localMatrix, localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp index 4709da59850..3405e5cdd6d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp @@ -715,6 +715,7 @@ AssemblerKernelHelper:: arrayView1d< real64 const > const & mimFaceGravCoef, arraySlice1d< localIndex const > const & elemToFaces, real64 const & elemGravCoef, + integer const useTotalMassEquation, ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, @@ -864,11 +865,14 @@ AssemblerKernelHelper:: } - // Apply equation/variable change transformation(s) - real64 work[NDOF*(NF+1)]; - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NDOF * ( NF + 1 ), dDivMassFluxes_dElemVars, work ); - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NF, dDivMassFluxes_dFaceVars, work ); - shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, divMassFluxes ); + if( useTotalMassEquation > 0 ) + { + // Apply equation/variable change transformation(s) + real64 work[NDOF * ( NF + 1 )]; + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NDOF * ( NF + 1 ), dDivMassFluxes_dElemVars, work ); + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NF, dDivMassFluxes_dFaceVars, work ); + shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, divMassFluxes ); + } // we are ready to assemble the local flux and its derivatives // no need for atomic adds - each row is assembled by a single thread @@ -1153,6 +1157,7 @@ AssemblerKernelHelper:: arrayView1d< real64 const > const & mimFaceGravCoef, \ arraySlice1d< localIndex const > const & elemToFaces, \ real64 const & elemGravCoef, \ + integer const useTotalMassEquation, \ ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, \ ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, \ ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, \ @@ -1281,6 +1286,7 @@ AssemblerKernel:: arraySlice1d< localIndex const > const & elemToFaces, real64 const & elemPres, real64 const & elemGravCoef, + integer const useTotalMassEquation, ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, @@ -1354,6 +1360,7 @@ AssemblerKernel:: mimFaceGravCoef, elemToFaces, elemGravCoef, + useTotalMassEquation, phaseDens, dPhaseDens, phaseMassDens, @@ -1407,6 +1414,7 @@ AssemblerKernel:: arraySlice1d< localIndex const > const & elemToFaces, \ real64 const & elemPres, \ real64 const & elemGravCoef, \ + integer const useTotalMassEquation, \ ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, \ ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, \ ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, \ @@ -1498,6 +1506,7 @@ FluxKernel:: globalIndex const rankOffset, real64 const lengthTolerance, real64 const dt, + integer const useTotalMassEquation, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -1575,6 +1584,7 @@ FluxKernel:: elemToFaces[ei], elemPres[ei], elemGravCoef[ei], + useTotalMassEquation, phaseDens, dPhaseDens, phaseMassDens, @@ -1627,6 +1637,7 @@ FluxKernel:: globalIndex const rankOffset, \ real64 const lengthTolerance, \ real64 const dt, \ + integer const useTotalMassEquation, \ CRSMatrixView< real64, globalIndex const > const & localMatrix, \ arrayView1d< real64 > const & localRhs ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp index b934ce3699d..3c07c9a99bf 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp @@ -335,6 +335,7 @@ struct AssemblerKernelHelper arrayView1d< real64 const > const & mimFaceGravCoef, arraySlice1d< localIndex const > const & elemToFaces, real64 const & elemGravCoef, + integer const useTotalMassEquation, ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, @@ -522,6 +523,7 @@ struct AssemblerKernel arraySlice1d< localIndex const > const & elemToFaces, real64 const & elemPres, real64 const & elemGravCoef, + const integer useTotalMassEquation, ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, @@ -636,6 +638,7 @@ struct FluxKernel globalIndex const rankOffset, real64 const lengthTolerance, real64 const dt, + integer const useTotalMassEquation, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp index acf8e11e879..20f896eaf52 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp @@ -44,6 +44,19 @@ using namespace constitutive; static constexpr real64 minDensForDivision = 1e-10; +enum class ElementBasedAssemblyKernelFlags +{ + SimpleAccumulation = 1 << 0, // 1 + TotalMassEquation = 1 << 1, // 2 + /// Add more flags like that if needed: + // Flag3 = 1 << 2, // 4 + // Flag4 = 1 << 3, // 8 + // Flag5 = 1 << 4, // 16 + // Flag6 = 1 << 5, // 32 + // Flag7 = 1 << 6, // 64 + // Flag8 = 1 << 7 //128 +}; + /******************************** PropertyKernelBase ********************************/ /** @@ -567,7 +580,8 @@ class ElementBasedAssemblyKernel MultiFluidBase const & fluid, CoupledSolidBase const & solid, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< ElementBasedAssemblyKernelFlags > const kernelFlags ) : m_numPhases( numPhases ), m_rankOffset( rankOffset ), m_dofNumber( subRegion.getReference< array1d< globalIndex > >( dofKey ) ), @@ -586,8 +600,11 @@ class ElementBasedAssemblyKernel m_phaseCompFrac_n( fluid.phaseCompFraction_n() ), m_phaseCompFrac( fluid.phaseCompFraction() ), m_dPhaseCompFrac( fluid.dPhaseCompFraction() ), + m_compDens( subRegion.getField< fields::flow::globalCompDensity >() ), + m_compDens_n( subRegion.getField< fields::flow::globalCompDensity_n >() ), m_localMatrix( localMatrix ), - m_localRhs( localRhs ) + m_localRhs( localRhs ), + m_kernelFlags( kernelFlags ) {} /** @@ -670,6 +687,12 @@ class ElementBasedAssemblyKernel StackVariables & stack, FUNC && phaseAmountKernelOp = NoOpFunc{} ) const { + if( m_kernelFlags.isSet( ElementBasedAssemblyKernelFlags::SimpleAccumulation ) ) + { + computeAccumulationSimple( ei, stack ); + return; + } + using Deriv = multifluid::DerivativeOffset; // construct the slices for variables accessed multiple times @@ -732,9 +755,9 @@ class ElementBasedAssemblyKernel { real64 const dPhaseCompAmount_dC = dPhaseCompFrac_dC[jc] * phaseAmount + phaseCompFrac[ip][ic] * dPhaseAmount_dC[jc]; + stack.localJacobian[ic][jc + 1] += dPhaseCompAmount_dC; } - } // call the lambda in the phase loop to allow the reuse of the phase amounts and their derivatives @@ -742,6 +765,41 @@ class ElementBasedAssemblyKernel phaseAmountKernelOp( ip, phaseAmount, phaseAmount_n, dPhaseAmount_dP, dPhaseAmount_dC ); } + + // check zero diagonal (works only in debug) + for( integer ic = 0; ic < numComp; ++ic ) + { + GEOS_ASSERT_MSG ( LvArray::math::abs( stack.localJacobian[ic][ic] ) > minDensForDivision, + GEOS_FMT( "Zero diagonal in Jacobian: equation {}, value = {}", ic, stack.localJacobian[ic][ic] ) ); + } + } + + template< typename FUNC = NoOpFunc > + GEOS_HOST_DEVICE + void computeAccumulationSimple( localIndex const ei, + StackVariables & stack ) const + { + // ic - index of component whose conservation equation is assembled + // (i.e. row number in local matrix) + for( integer ic = 0; ic < numComp; ++ic ) + { + real64 const compAmount = stack.poreVolume * m_compDens[ei][ic]; + real64 const compAmount_n = stack.poreVolume_n * m_compDens_n[ei][ic]; + + stack.localResidual[ic] += compAmount - compAmount_n; + + // Pavel: commented below is some experiment, needs to be re-tested + //real64 const compDens = (ic == 0 && m_compDens[ei][ic] < 1e-6) ? 1e-3 : m_compDens[ei][ic]; + real64 const dCompAmount_dP = stack.dPoreVolume_dPres * m_compDens[ei][ic]; + GEOS_ASSERT_MSG( !(ic == 0 && LvArray::math::abs( dCompAmount_dP ) < minDensForDivision), + GEOS_FMT( "Zero diagonal in Jacobian: equation {}, value = {}", ic, dCompAmount_dP ) ); + stack.localJacobian[ic][0] += dCompAmount_dP; + + real64 const dCompAmount_dC = stack.poreVolume; + GEOS_ASSERT_MSG( !(ic == ic + 1 && LvArray::math::abs( dCompAmount_dC ) < minDensForDivision), + GEOS_FMT( "Zero diagonal in Jacobian: equation {}, value = {}", ic, dCompAmount_dC ) ); + stack.localJacobian[ic][ic + 1] += dCompAmount_dC; + } } /** @@ -800,16 +858,20 @@ class ElementBasedAssemblyKernel { using namespace compositionalMultiphaseUtilities; - // apply equation/variable change transformation to the component mass balance equations - real64 work[numDof]{}; - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof, stack.localJacobian, work ); - shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localResidual ); + if( m_kernelFlags.isSet( ElementBasedAssemblyKernelFlags::TotalMassEquation ) ) + { + // apply equation/variable change transformation to the component mass balance equations + real64 work[numDof]{}; + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof, stack.localJacobian, work ); + shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localResidual ); + } // add contribution to residual and jacobian into: // - the component mass balance equations (i = 0 to i = numComp-1) // - the volume balance equations (i = numComp) // note that numDof includes derivatives wrt temperature if this class is derived in ThermalKernels - for( integer i = 0; i < numComp+1; ++i ) + integer const numRows = numComp+1; + for( integer i = 0; i < numRows; ++i ) { m_localRhs[stack.localRow + i] += stack.localResidual[i]; m_localMatrix.addToRow< serialAtomic >( stack.localRow + i, @@ -889,11 +951,16 @@ class ElementBasedAssemblyKernel arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const m_phaseCompFrac; arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > const m_dPhaseCompFrac; + // Views on component densities + arrayView2d< real64 const, compflow::USD_COMP > m_compDens; + arrayView2d< real64 const, compflow::USD_COMP > m_compDens_n; + /// View on the local CRS matrix CRSMatrixView< real64, globalIndex const > const m_localMatrix; /// View on the local RHS arrayView1d< real64 > const m_localRhs; + BitFlags< ElementBasedAssemblyKernelFlags > const m_kernelFlags; }; /** @@ -921,6 +988,8 @@ class ElementBasedAssemblyKernelFactory createAndLaunch( integer const numComps, integer const numPhases, globalIndex const rankOffset, + integer const useTotalMassEquation, + integer const useSimpleAccumulation, string const dofKey, ElementSubRegionBase const & subRegion, MultiFluidBase const & fluid, @@ -932,8 +1001,15 @@ class ElementBasedAssemblyKernelFactory { integer constexpr NUM_COMP = NC(); integer constexpr NUM_DOF = NC()+1; + + BitFlags< ElementBasedAssemblyKernelFlags > kernelFlags; + if( useTotalMassEquation ) + kernelFlags.set( ElementBasedAssemblyKernelFlags::TotalMassEquation ); + if( useSimpleAccumulation ) + kernelFlags.set( ElementBasedAssemblyKernelFlags::SimpleAccumulation ); + ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > - kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs, kernelFlags ); ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp index a4fcb82bace..f63126c5cbb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernelUtilities.hpp @@ -218,7 +218,6 @@ struct PPUPhaseFlux compute( integer const numPhase, integer const ip, integer const hasCapPressure, - //real64 const GEOS_UNUSED_PARAM( epsC1PPU ), localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp index 601c42cd589..e6636d930fa 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp @@ -38,9 +38,10 @@ FaceBasedAssemblyKernelBase::FaceBasedAssemblyKernelBase( integer const numPhase DofNumberAccessor const & dofNumberAccessor, CompFlowAccessors const & compFlowAccessors, MultiFluidAccessors const & multiFluidAccessors, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< FaceBasedAssemblyKernelFlags > kernelFlags ) : m_numPhases( numPhases ), m_rankOffset( rankOffset ), m_dt( dt ), @@ -53,7 +54,8 @@ FaceBasedAssemblyKernelBase::FaceBasedAssemblyKernelBase( integer const numPhase m_phaseCompFrac( multiFluidAccessors.get( fields::multifluid::phaseCompFraction {} ) ), m_dPhaseCompFrac( multiFluidAccessors.get( fields::multifluid::dPhaseCompFraction {} ) ), m_localMatrix( localMatrix ), - m_localRhs( localRhs ) + m_localRhs( localRhs ), + m_kernelFlags( kernelFlags ) {} /******************************** CFLFluxKernel ********************************/ @@ -65,7 +67,7 @@ void CFLFluxKernel:: compute( integer const numPhases, localIndex const stencilSize, - real64 const & dt, + real64 const dt, arraySlice1d< localIndex const > const seri, arraySlice1d< localIndex const > const sesri, arraySlice1d< localIndex const > const sei, @@ -156,7 +158,7 @@ template< integer NC, typename STENCILWRAPPER_TYPE > void CFLFluxKernel:: launch( integer const numPhases, - real64 const & dt, + real64 const dt, STENCILWRAPPER_TYPE const & stencilWrapper, ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, @@ -214,7 +216,7 @@ CFLFluxKernel:: template \ void CFLFluxKernel:: \ launch< NC, STENCILWRAPPER_TYPE >( integer const numPhases, \ - real64 const & dt, \ + real64 const dt, \ STENCILWRAPPER_TYPE const & stencil, \ ElementViewConst< arrayView1d< real64 const > > const & pres, \ ElementViewConst< arrayView1d< real64 const > > const & gravCoef, \ @@ -261,7 +263,7 @@ template< integer NP > GEOS_HOST_DEVICE void CFLKernel:: - computePhaseCFL( real64 const & poreVol, + computePhaseCFL( real64 const poreVol, arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac, arraySlice1d< real64 const, relperm::USD_RELPERM - 2 > phaseRelPerm, arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > dPhaseRelPerm_dPhaseVolFrac, @@ -345,7 +347,7 @@ template< integer NC > GEOS_HOST_DEVICE void CFLKernel:: - computeCompCFL( real64 const & poreVol, + computeCompCFL( real64 const poreVol, arraySlice1d< real64 const, compflow::USD_COMP - 1 > compDens, arraySlice1d< real64 const, compflow::USD_COMP - 1 > compFrac, arraySlice1d< real64 const, compflow::USD_COMP - 1 > compOutflux, @@ -462,9 +464,9 @@ AquiferBCKernel:: compute( integer const numPhases, integer const ipWater, bool const allowAllPhasesIntoAquifer, - real64 const & aquiferVolFlux, - real64 const & dAquiferVolFlux_dPres, - real64 const & aquiferWaterPhaseDens, + real64 const aquiferVolFlux, + real64 const dAquiferVolFlux_dPres, + real64 const aquiferWaterPhaseDens, arrayView1d< real64 const > const & aquiferWaterPhaseCompFrac, arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseDens, arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseDens, @@ -473,7 +475,7 @@ AquiferBCKernel:: arraySlice2d< real64 const, multifluid::USD_PHASE_COMP - 2 > phaseCompFrac, arraySlice3d< real64 const, multifluid::USD_PHASE_COMP_DC - 2 > dPhaseCompFrac, arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dCompFrac_dCompDens, - real64 const & dt, + real64 const dt, real64 (& localFlux)[NC], real64 (& localFluxJacobian)[NC][NC+1] ) { @@ -541,11 +543,12 @@ AquiferBCKernel:: launch( integer const numPhases, integer const ipWater, bool const allowAllPhasesIntoAquifer, + integer const useTotalMassEquation, BoundaryStencil const & stencil, globalIndex const rankOffset, ElementViewConst< arrayView1d< globalIndex const > > const & dofNumber, AquiferBoundaryCondition::KernelWrapper const & aquiferBCWrapper, - real64 const & aquiferWaterPhaseDens, + real64 const aquiferWaterPhaseDens, arrayView1d< real64 const > const & aquiferWaterPhaseCompFrac, ElementViewConst< arrayView1d< integer const > > const & ghostRank, ElementViewConst< arrayView1d< real64 const > > const & pres, @@ -558,8 +561,8 @@ AquiferBCKernel:: ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac, ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, - real64 const & timeAtBeginningOfStep, - real64 const & dt, + real64 const timeAtBeginningOfStep, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -622,11 +625,13 @@ AquiferBCKernel:: dofColIndices[jdof] = offset + jdof; } - // Apply equation/variable change transformation(s) - real64 work[NDOF]; - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NDOF, localFluxJacobian, work ); - shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, localFlux ); - + if( useTotalMassEquation > 0 ) + { + // Apply equation/variable change transformation(s) + real64 work[NDOF]; + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NDOF, localFluxJacobian, work ); + shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, localFlux ); + } // Add to residual/jacobian if( ghostRank[er][esr][ei] < 0 ) @@ -654,11 +659,12 @@ AquiferBCKernel:: launch< NC >( integer const numPhases, \ integer const ipWater, \ bool const allowAllPhasesIntoAquifer, \ + integer const useTotalMassEquation, \ BoundaryStencil const & stencil, \ globalIndex const rankOffset, \ ElementViewConst< arrayView1d< globalIndex const > > const & dofNumber, \ AquiferBoundaryCondition::KernelWrapper const & aquiferBCWrapper, \ - real64 const & aquiferWaterPhaseDens, \ + real64 const aquiferWaterPhaseDens, \ arrayView1d< real64 const > const & aquiferWaterPhaseCompFrac, \ ElementViewConst< arrayView1d< integer const > > const & ghostRank, \ ElementViewConst< arrayView1d< real64 const > > const & pres, \ @@ -671,8 +677,8 @@ AquiferBCKernel:: ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, \ ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac, \ ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, \ - real64 const & timeAtBeginningOfStep, \ - real64 const & dt, \ + real64 const timeAtBeginningOfStep, \ + real64 const dt, \ CRSMatrixView< real64, globalIndex const > const & localMatrix, \ arrayView1d< real64 > const & localRhs ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp index 69e9e06c818..a22a30a9736 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp @@ -57,6 +57,22 @@ namespace isothermalCompositionalMultiphaseFVMKernels using namespace constitutive; +enum class FaceBasedAssemblyKernelFlags +{ + /// Flag to specify whether capillary pressure is used or not + CapPressure = 1 << 0, // 1 + /// Flag indicating whether total mass equation is formed or not + TotalMassEquation = 1 << 1, // 2 + /// Flag indicating whether C1-PPU is used or not + C1PPU = 1 << 2, // 4 + /// Add more flags like that if needed: + // Flag4 = 1 << 3, // 8 + // Flag5 = 1 << 4, // 16 + // Flag6 = 1 << 5, // 32 + // Flag7 = 1 << 6, // 64 + // Flag8 = 1 << 7 //128 +}; + /******************************** PhaseMobilityKernel ********************************/ /** @@ -322,15 +338,17 @@ class FaceBasedAssemblyKernelBase * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector + * @param[in] kernelFlags flags packed all together */ FaceBasedAssemblyKernelBase( integer const numPhases, globalIndex const rankOffset, DofNumberAccessor const & dofNumberAccessor, CompFlowAccessors const & compFlowAccessors, MultiFluidAccessors const & multiFluidAccessors, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); + arrayView1d< real64 > const & localRhs, + BitFlags< FaceBasedAssemblyKernelFlags > kernelFlags ); protected: @@ -369,6 +387,8 @@ class FaceBasedAssemblyKernelBase CRSMatrixView< real64, globalIndex const > const m_localMatrix; /// View on the local RHS arrayView1d< real64 > const m_localRhs; + + BitFlags< FaceBasedAssemblyKernelFlags > const m_kernelFlags; }; /** @@ -376,10 +396,9 @@ class FaceBasedAssemblyKernelBase * @tparam NUM_COMP number of fluid components * @tparam NUM_DOF number of degrees of freedom * @tparam STENCILWRAPPER the type of the stencil wrapper - * @tparam PHASE_FLUX_COMPUTE the method employed to compute the flux * @brief Define the interface for the assembly kernel in charge of flux terms */ -template< integer NUM_COMP, integer NUM_DOF, typename STENCILWRAPPER, typename PHASE_FLUX_COMPUTE = isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux > +template< integer NUM_COMP, integer NUM_DOF, typename STENCILWRAPPER > class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase { public: @@ -409,7 +428,6 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase * @brief Constructor for the kernel interface * @param[in] numPhases the number of fluid phases * @param[in] rankOffset the offset of my MPI rank - * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not * @param[in] stencilWrapper reference to the stencil wrapper * @param[in] dofNumberAccessor * @param[in] compFlowAccessors @@ -419,19 +437,20 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector + * @param[in] kernelFlags flags packed together */ FaceBasedAssemblyKernel( integer const numPhases, globalIndex const rankOffset, - integer const hasCapPressure, STENCILWRAPPER const & stencilWrapper, DofNumberAccessor const & dofNumberAccessor, CompFlowAccessors const & compFlowAccessors, MultiFluidAccessors const & multiFluidAccessors, CapPressureAccessors const & capPressureAccessors, PermeabilityAccessors const & permeabilityAccessors, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< FaceBasedAssemblyKernelFlags > kernelFlags ) : FaceBasedAssemblyKernelBase( numPhases, rankOffset, dofNumberAccessor, @@ -439,8 +458,8 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase multiFluidAccessors, dt, localMatrix, - localRhs ), - m_hasCapPressure( hasCapPressure ), + localRhs, + kernelFlags ), m_permeability( permeabilityAccessors.get( fields::permeability::permeability {} ) ), m_dPerm_dPres( permeabilityAccessors.get( fields::permeability::dPerm_dPressure {} ) ), m_phaseMob( compFlowAccessors.get( fields::flow::phaseMobility {} ) ), @@ -602,25 +621,50 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase localIndex k_up = -1; - PHASE_FLUX_COMPUTE::template compute< numComp, numFluxSupportPoints > - ( m_numPhases, - ip, - m_hasCapPressure, - seri, sesri, sei, - trans, - dTrans_dPres, - m_pres, - m_gravCoef, - m_phaseMob, m_dPhaseMob, - m_dPhaseVolFrac, - m_dCompFrac_dCompDens, - m_phaseMassDens, m_dPhaseMassDens, - m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac, - k_up, - potGrad, - phaseFlux, - dPhaseFlux_dP, - dPhaseFlux_dC ); + if( m_kernelFlags.isSet( FaceBasedAssemblyKernelFlags::C1PPU )) + { + isothermalCompositionalMultiphaseFVMKernelUtilities::C1PPUPhaseFlux::compute< numComp, numFluxSupportPoints > + ( m_numPhases, + ip, + m_kernelFlags.isSet( FaceBasedAssemblyKernelFlags::CapPressure ), + seri, sesri, sei, + trans, + dTrans_dPres, + m_pres, + m_gravCoef, + m_phaseMob, m_dPhaseMob, + m_dPhaseVolFrac, + m_dCompFrac_dCompDens, + m_phaseMassDens, m_dPhaseMassDens, + m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac, + k_up, + potGrad, + phaseFlux, + dPhaseFlux_dP, + dPhaseFlux_dC ); + } + else + { + isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux::compute< numComp, numFluxSupportPoints > + ( m_numPhases, + ip, + m_kernelFlags.isSet( FaceBasedAssemblyKernelFlags::CapPressure ), + seri, sesri, sei, + trans, + dTrans_dPres, + m_pres, + m_gravCoef, + m_phaseMob, m_dPhaseMob, + m_dPhaseVolFrac, + m_dCompFrac_dCompDens, + m_phaseMassDens, m_dPhaseMassDens, + m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac, + k_up, + potGrad, + phaseFlux, + dPhaseFlux_dP, + dPhaseFlux_dC ); + } isothermalCompositionalMultiphaseFVMKernelUtilities:: PhaseComponentFlux::compute< numComp, numFluxSupportPoints > @@ -683,12 +727,15 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase { using namespace compositionalMultiphaseUtilities; - // Apply equation/variable change transformation(s) - stackArray1d< real64, maxStencilSize * numDof > work( stack.stencilSize * numDof ); - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numEqn, numDof*stack.stencilSize, stack.numConnectedElems, - stack.localFluxJacobian, work ); - shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, numEqn, stack.numConnectedElems, - stack.localFlux ); + if( m_kernelFlags.isSet( FaceBasedAssemblyKernelFlags::TotalMassEquation ) ) + { + // Apply equation/variable change transformation(s) + stackArray1d< real64, maxStencilSize * numDof > work( stack.stencilSize * numDof ); + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numEqn, numDof * stack.stencilSize, stack.numConnectedElems, + stack.localFluxJacobian, work ); + shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, numEqn, stack.numConnectedElems, + stack.localFlux ); + } // add contribution to residual and jacobian into: // - the component mass balance equations (i = 0 to i = numComp-1) @@ -744,9 +791,6 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase protected: - /// Flag to specify whether capillary pressure is used or not - integer const m_hasCapPressure; - /// Views on permeability ElementViewConst< arrayView3d< real64 const > > const m_permeability; ElementViewConst< arrayView3d< real64 const > > const m_dPerm_dPres; @@ -805,11 +849,12 @@ class FaceBasedAssemblyKernelFactory globalIndex const rankOffset, string const & dofKey, integer const hasCapPressure, + integer const useTotalMassEquation, UpwindingParameters upwindingParams, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -822,34 +867,25 @@ class FaceBasedAssemblyKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + BitFlags< FaceBasedAssemblyKernelFlags > kernelFlags; + if( hasCapPressure ) + kernelFlags.set( FaceBasedAssemblyKernelFlags::CapPressure ); + if( useTotalMassEquation ) + kernelFlags.set( FaceBasedAssemblyKernelFlags::TotalMassEquation ); if( upwindingParams.upwindingScheme == UpwindingScheme::C1PPU && isothermalCompositionalMultiphaseFVMKernelUtilities::epsC1PPU > 0 ) - { - using kernelType = - FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER, isothermalCompositionalMultiphaseFVMKernelUtilities::C1PPUPhaseFlux >; - typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); - typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); - typename kernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName ); - typename kernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); + kernelFlags.set( FaceBasedAssemblyKernelFlags::C1PPU ); - kernelType kernel( numPhases, rankOffset, hasCapPressure, stencilWrapper, dofNumberAccessor, - compFlowAccessors, multiFluidAccessors, capPressureAccessors, permeabilityAccessors, - dt, localMatrix, localRhs ); - kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); - } - else - { - using kernelType = FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; - typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); - typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); - typename kernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName ); - typename kernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); + using kernelType = FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; + typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); + typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); + typename kernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName ); + typename kernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); - kernelType kernel( numPhases, rankOffset, hasCapPressure, stencilWrapper, dofNumberAccessor, - compFlowAccessors, multiFluidAccessors, capPressureAccessors, permeabilityAccessors, - dt, localMatrix, localRhs ); - kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); - } + kernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor, + compFlowAccessors, multiFluidAccessors, capPressureAccessors, permeabilityAccessors, + dt, localMatrix, localRhs, kernelFlags ); + kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } }; @@ -891,6 +927,7 @@ class DiffusionDispersionFaceBasedAssemblyKernel : public FaceBasedAssemblyKerne using AbstractBase = isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelBase; using AbstractBase::m_dPhaseVolFrac; + using AbstractBase::m_kernelFlags; using DiffusionAccessors = StencilMaterialAccessors< DiffusionBase, @@ -920,6 +957,7 @@ class DiffusionDispersionFaceBasedAssemblyKernel : public FaceBasedAssemblyKerne * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector + * @param[in] kernelFlags flags packed together */ DiffusionDispersionFaceBasedAssemblyKernel( integer const numPhases, globalIndex const rankOffset, @@ -930,9 +968,10 @@ class DiffusionDispersionFaceBasedAssemblyKernel : public FaceBasedAssemblyKerne DiffusionAccessors const & diffusionAccessors, DispersionAccessors const & dispersionAccessors, PorosityAccessors const & porosityAccessors, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< FaceBasedAssemblyKernelFlags > kernelFlags ) : FaceBasedAssemblyKernelBase( numPhases, rankOffset, dofNumberAccessor, @@ -940,7 +979,8 @@ class DiffusionDispersionFaceBasedAssemblyKernel : public FaceBasedAssemblyKerne multiFluidAccessors, dt, localMatrix, - localRhs ), + localRhs, + kernelFlags ), m_phaseVolFrac( compFlowAccessors.get( fields::flow::phaseVolumeFraction {} ) ), m_phaseDens( multiFluidAccessors.get( fields::multifluid::phaseDensity {} ) ), m_dPhaseDens( multiFluidAccessors.get( fields::multifluid::dPhaseDensity {} ) ), @@ -1402,12 +1442,15 @@ class DiffusionDispersionFaceBasedAssemblyKernel : public FaceBasedAssemblyKerne { using namespace compositionalMultiphaseUtilities; - // Apply equation/variable change transformation(s) - stackArray1d< real64, maxStencilSize * numDof > work( stack.stencilSize * numDof ); - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numEqn, numDof*stack.stencilSize, stack.numConnectedElems, - stack.localFluxJacobian, work ); - shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, numEqn, stack.numConnectedElems, - stack.localFlux ); + if( m_kernelFlags.isSet( FaceBasedAssemblyKernelFlags::TotalMassEquation ) ) + { + // Apply equation/variable change transformation(s) + stackArray1d< real64, maxStencilSize * numDof > work( stack.stencilSize * numDof ); + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numEqn, numDof*stack.stencilSize, stack.numConnectedElems, + stack.localFluxJacobian, work ); + shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, numEqn, stack.numConnectedElems, + stack.localFlux ); + } // add contribution to residual and jacobian into: // - the component mass balance equations (i = 0 to i = numComp-1) @@ -1447,8 +1490,8 @@ class DiffusionDispersionFaceBasedAssemblyKernel : public FaceBasedAssemblyKerne template< typename POLICY, typename KERNEL_TYPE > static void launch( localIndex const numConnections, - integer const & hasDiffusion, - integer const & hasDispersion, + integer const hasDiffusion, + integer const hasDispersion, KERNEL_TYPE const & kernelComponent ) { GEOS_MARK_FUNCTION; @@ -1532,12 +1575,13 @@ class DiffusionDispersionFaceBasedAssemblyKernelFactory integer const numPhases, globalIndex const rankOffset, string const & dofKey, - integer const & hasDiffusion, - integer const & hasDispersion, + integer const hasDiffusion, + integer const hasDispersion, + integer const useTotalMassEquation, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -1550,6 +1594,10 @@ class DiffusionDispersionFaceBasedAssemblyKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + BitFlags< FaceBasedAssemblyKernelFlags > kernelFlags; + if( useTotalMassEquation ) + kernelFlags.set( FaceBasedAssemblyKernelFlags::TotalMassEquation ); + using kernelType = DiffusionDispersionFaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); @@ -1560,7 +1608,7 @@ class DiffusionDispersionFaceBasedAssemblyKernelFactory kernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor, compFlowAccessors, multiFluidAccessors, diffusionAccessors, dispersionAccessors, porosityAccessors, - dt, localMatrix, localRhs ); + dt, localMatrix, localRhs, kernelFlags ); kernelType::template launch< POLICY >( stencilWrapper.size(), hasDiffusion, hasDispersion, kernel ); @@ -1613,6 +1661,7 @@ class DirichletFaceBasedAssemblyKernel : public FaceBasedAssemblyKernel< NUM_COM using AbstractBase::m_dCompFrac_dCompDens; using AbstractBase::m_localMatrix; using AbstractBase::m_localRhs; + using AbstractBase::m_kernelFlags; using Base = isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, BoundaryStencilWrapper >; using Base::numComp; @@ -1633,7 +1682,6 @@ class DirichletFaceBasedAssemblyKernel : public FaceBasedAssemblyKernel< NUM_COM * @brief Constructor for the kernel interface * @param[in] numPhases the number of fluid phases * @param[in] rankOffset the offset of my MPI rank - * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not * @param[in] faceManager the face manager * @param[in] stencilWrapper reference to the stencil wrapper * @param[in] fluidWrapper reference to the fluid wrapper @@ -1645,10 +1693,10 @@ class DirichletFaceBasedAssemblyKernel : public FaceBasedAssemblyKernel< NUM_COM * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector + * @param[in] kernelFlags flags packed together */ DirichletFaceBasedAssemblyKernel( integer const numPhases, globalIndex const rankOffset, - integer const hasCapPressure, FaceManager const & faceManager, BoundaryStencilWrapper const & stencilWrapper, FLUIDWRAPPER const & fluidWrapper, @@ -1657,12 +1705,12 @@ class DirichletFaceBasedAssemblyKernel : public FaceBasedAssemblyKernel< NUM_COM MultiFluidAccessors const & multiFluidAccessors, CapPressureAccessors const & capPressureAccessors, PermeabilityAccessors const & permeabilityAccessors, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< FaceBasedAssemblyKernelFlags > kernelFlags ) : Base( numPhases, rankOffset, - hasCapPressure, stencilWrapper, dofNumberAccessor, compFlowAccessors, @@ -1671,7 +1719,8 @@ class DirichletFaceBasedAssemblyKernel : public FaceBasedAssemblyKernel< NUM_COM permeabilityAccessors, dt, localMatrix, - localRhs ), + localRhs, + kernelFlags ), m_facePres( faceManager.getField< fields::flow::facePressure >() ), m_faceTemp( faceManager.getField< fields::flow::faceTemperature >() ), m_faceCompFrac( faceManager.getField< fields::flow::faceGlobalCompFraction >() ), @@ -1781,7 +1830,7 @@ class DirichletFaceBasedAssemblyKernel : public FaceBasedAssemblyKernel< NUM_COM StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, multifluid::LAYOUT_PHASE > facePhaseVisc( 1, 1, m_numPhases ); StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, multifluid::LAYOUT_PHASE > facePhaseEnthalpy( 1, 1, m_numPhases ); StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, multifluid::LAYOUT_PHASE > facePhaseInternalEnergy( 1, 1, m_numPhases ); - StackArray< real64, 4, constitutive::MultiFluidBase::MAX_NUM_PHASES *NUM_COMP, + StackArray< real64, 4, constitutive::MultiFluidBase::MAX_NUM_PHASES * NUM_COMP, multifluid::LAYOUT_PHASE_COMP > facePhaseCompFrac( 1, 1, m_numPhases, NUM_COMP ); real64 faceTotalDens = 0.0; @@ -1959,10 +2008,13 @@ class DirichletFaceBasedAssemblyKernel : public FaceBasedAssemblyKernel< NUM_COM using namespace compositionalMultiphaseUtilities; using Order = BoundaryStencil::Order; - // Apply equation/variable change transformation(s) - real64 work[numDof]{}; - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof, stack.localFluxJacobian, work ); - shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localFlux ); + if( AbstractBase::m_kernelFlags.isSet( FaceBasedAssemblyKernelFlags::TotalMassEquation ) ) + { + // Apply equation/variable change transformation(s) + real64 work[numDof]{}; + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof, stack.localFluxJacobian, work ); + shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localFlux ); + } // add contribution to residual and jacobian into: // - the component mass balance equations (i = 0 to i = numComp-1) @@ -2033,13 +2085,14 @@ class DirichletFaceBasedAssemblyKernelFactory createAndLaunch( integer const numComps, integer const numPhases, globalIndex const rankOffset, + integer const useTotalMassEquation, string const & dofKey, string const & solverName, FaceManager const & faceManager, ElementRegionManager const & elemManager, BoundaryStencilWrapper const & stencilWrapper, MultiFluidBase & fluidBase, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -2057,18 +2110,20 @@ class DirichletFaceBasedAssemblyKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + // for now, we neglect capillary pressure in the kernel + BitFlags< FaceBasedAssemblyKernelFlags > kernelFlags; + if( useTotalMassEquation ) + kernelFlags.set( FaceBasedAssemblyKernelFlags::TotalMassEquation ); + using kernelType = DirichletFaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, typename FluidType::KernelWrapper >; typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); typename kernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName ); typename kernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); - // for now, we neglect capillary pressure in the kernel - bool const hasCapPressure = false; - - kernelType kernel( numPhases, rankOffset, hasCapPressure, faceManager, stencilWrapper, fluidWrapper, + kernelType kernel( numPhases, rankOffset, faceManager, stencilWrapper, fluidWrapper, dofNumberAccessor, compFlowAccessors, multiFluidAccessors, capPressureAccessors, permeabilityAccessors, - dt, localMatrix, localRhs ); + dt, localMatrix, localRhs, kernelFlags ); kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } ); @@ -2125,7 +2180,7 @@ struct CFLFluxKernel static void compute( integer const numPhases, localIndex const stencilSize, - real64 const & dt, + real64 const dt, arraySlice1d< localIndex const > const seri, arraySlice1d< localIndex const > const sesri, arraySlice1d< localIndex const > const sei, @@ -2144,7 +2199,7 @@ struct CFLFluxKernel template< integer NC, typename STENCILWRAPPER_TYPE > static void launch( integer const numPhases, - real64 const & dt, + real64 const dt, STENCILWRAPPER_TYPE const & stencil, ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, @@ -2175,7 +2230,7 @@ struct CFLKernel GEOS_HOST_DEVICE inline static void - computePhaseCFL( real64 const & poreVol, + computePhaseCFL( real64 const poreVol, arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac, arraySlice1d< real64 const, relperm::USD_RELPERM - 2 > phaseRelPerm, arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > dPhaseRelPerm_dPhaseVolFrac, @@ -2187,7 +2242,7 @@ struct CFLKernel GEOS_HOST_DEVICE inline static void - computeCompCFL( real64 const & poreVol, + computeCompCFL( real64 const poreVol, arraySlice1d< real64 const, compflow::USD_COMP - 1 > compDens, arraySlice1d< real64 const, compflow::USD_COMP - 1 > compFrac, arraySlice1d< real64 const, compflow::USD_COMP - 1 > compOutflux, @@ -2253,9 +2308,9 @@ struct AquiferBCKernel compute( integer const numPhases, integer const ipWater, bool const allowAllPhasesIntoAquifer, - real64 const & aquiferVolFlux, - real64 const & dAquiferVolFlux_dPres, - real64 const & aquiferWaterPhaseDens, + real64 const aquiferVolFlux, + real64 const dAquiferVolFlux_dPres, + real64 const aquiferWaterPhaseDens, arrayView1d< real64 const > const & aquiferWaterPhaseCompFrac, arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseDens, arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseDens, @@ -2264,7 +2319,7 @@ struct AquiferBCKernel arraySlice2d< real64 const, multifluid::USD_PHASE_COMP - 2 > phaseCompFrac, arraySlice3d< real64 const, multifluid::USD_PHASE_COMP_DC - 2 > dPhaseCompFrac, arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dCompFrac_dCompDens, - real64 const & dt, + real64 const dt, real64 ( &localFlux )[NC], real64 ( &localFluxJacobian )[NC][NC+1] ); @@ -2273,11 +2328,12 @@ struct AquiferBCKernel launch( integer const numPhases, integer const ipWater, bool const allowAllPhasesIntoAquifer, + integer const useTotalMassEquation, BoundaryStencil const & stencil, globalIndex const rankOffset, ElementViewConst< arrayView1d< globalIndex const > > const & dofNumber, AquiferBoundaryCondition::KernelWrapper const & aquiferBCWrapper, - real64 const & aquiferWaterPhaseDens, + real64 const aquiferWaterPhaseDens, arrayView1d< real64 const > const & aquiferWaterPhaseCompFrac, ElementViewConst< arrayView1d< integer const > > const & ghostRank, ElementViewConst< arrayView1d< real64 const > > const & pres, @@ -2290,8 +2346,8 @@ struct AquiferBCKernel ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac, ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, - real64 const & timeAtBeginningOfStep, - real64 const & dt, + real64 const timeAtBeginningOfStep, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/StabilizedCompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/StabilizedCompositionalMultiphaseFVMKernels.hpp index a17b30883d3..b0d6328732f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/StabilizedCompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/StabilizedCompositionalMultiphaseFVMKernels.hpp @@ -75,6 +75,7 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne using AbstractBase::m_dt; using AbstractBase::m_numPhases; + using AbstractBase::m_kernelFlags; using AbstractBase::m_rankOffset; using AbstractBase::m_ghostRank; using AbstractBase::m_dofNumber; @@ -102,7 +103,6 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne * @brief Constructor for the kernel interface * @param[in] numPhases the number of fluid phases * @param[in] rankOffset the offset of my MPI rank - * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not * @param[in] stencilWrapper reference to the stencil wrapper * @param[in] dofNumberAccessor accessor for the dofs numbers * @param[in] compFlowAccessor accessor for wrappers registered by the solver @@ -115,10 +115,10 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector + * @param[in] kernelFlags flags packed together */ FaceBasedAssemblyKernel( integer const numPhases, globalIndex const rankOffset, - integer const hasCapPressure, STENCILWRAPPER const & stencilWrapper, DofNumberAccessor const & dofNumberAccessor, CompFlowAccessors const & compFlowAccessors, @@ -130,10 +130,10 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne RelPermAccessors const & relPermAccessors, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags > kernelFlags ) : Base( numPhases, rankOffset, - hasCapPressure, stencilWrapper, dofNumberAccessor, compFlowAccessors, @@ -142,7 +142,8 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne permeabilityAccessors, dt, localMatrix, - localRhs ), + localRhs, + kernelFlags ), m_pres_n( stabCompFlowAccessors.get( fields::flow::pressure_n {} ) ), m_phaseDens_n( stabMultiFluidAccessors.get( fields::multifluid::phaseDensity_n {} ) ), m_phaseCompFrac_n( stabMultiFluidAccessors.get( fields::multifluid::phaseCompFraction_n {} ) ), @@ -201,8 +202,8 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne localIndex const er_up, localIndex const esr_up, localIndex const ei_up, - real64 const & potGrad, - real64 const & phaseFlux, + real64 const potGrad, + real64 const phaseFlux, real64 const (&dPhaseFlux_dP)[2], real64 const (&dPhaseFlux_dC)[2][numComp] ) { @@ -326,10 +327,11 @@ class FaceBasedAssemblyKernelFactory globalIndex const rankOffset, string const & dofKey, integer const hasCapPressure, + integer const useTotalMassEquation, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -343,6 +345,12 @@ class FaceBasedAssemblyKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + BitFlags< isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags > kernelFlags; + if( hasCapPressure ) + kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags::CapPressure ); + if( useTotalMassEquation ) + kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags::TotalMassEquation ); + using KERNEL_TYPE = FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename KERNEL_TYPE::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename KERNEL_TYPE::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); @@ -352,10 +360,10 @@ class FaceBasedAssemblyKernelFactory typename KERNEL_TYPE::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); typename KERNEL_TYPE::RelPermAccessors relPermAccessors( elemManager, solverName ); - KERNEL_TYPE kernel( numPhases, rankOffset, hasCapPressure, stencilWrapper, dofNumberAccessor, + KERNEL_TYPE kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor, compFlowAccessors, stabCompFlowAccessors, multiFluidAccessors, stabMultiFluidAccessors, capPressureAccessors, permeabilityAccessors, relPermAccessors, - dt, localMatrix, localRhs ); + dt, localMatrix, localRhs, kernelFlags ); KERNEL_TYPE::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp index 42a9886eb7d..7ba1eb09d9e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp @@ -193,8 +193,9 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK MultiFluidBase const & fluid, CoupledSolidBase const & solid, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) - : Base( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ), + arrayView1d< real64 > const & localRhs, + BitFlags< isothermalCompositionalMultiphaseBaseKernels::ElementBasedAssemblyKernelFlags > const kernelFlags ) + : Base( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs, kernelFlags ), m_dPoro_dTemp( solid.getDporosity_dTemperature() ), m_phaseInternalEnergy_n( fluid.phaseInternalEnergy_n() ), m_phaseInternalEnergy( fluid.phaseInternalEnergy() ), @@ -431,6 +432,8 @@ class ElementBasedAssemblyKernelFactory createAndLaunch( localIndex const numComps, localIndex const numPhases, globalIndex const rankOffset, + integer const useTotalMassEquation, + integer const useSimpleAccumulation, string const dofKey, ElementSubRegionBase const & subRegion, MultiFluidBase const & fluid, @@ -443,8 +446,15 @@ class ElementBasedAssemblyKernelFactory { localIndex constexpr NUM_COMP = NC(); localIndex constexpr NUM_DOF = NC()+2; + + BitFlags< isothermalCompositionalMultiphaseBaseKernels::ElementBasedAssemblyKernelFlags > kernelFlags; + if( useTotalMassEquation ) + kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::ElementBasedAssemblyKernelFlags::TotalMassEquation ); + if( useSimpleAccumulation ) + kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::ElementBasedAssemblyKernelFlags::SimpleAccumulation ); + ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > - kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs, kernelFlags ); ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF >::template launch< POLICY, ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > >( subRegion.size(), kernel ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp index 4beb40cb608..28bd1d26eb6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp @@ -199,7 +199,6 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne using Base::maxNumConns; using Base::maxStencilSize; using Base::numFluxSupportPoints; - using Base::m_hasCapPressure; using Base::m_phaseMob; using Base::m_dPhaseMob; using Base::m_dPhaseMassDens; @@ -226,7 +225,6 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne * @brief Constructor for the kernel interface * @param[in] numPhases the number of fluid phases * @param[in] rankOffset the offset of my MPI rank - * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not * @param[in] stencilWrapper reference to the stencil wrapper * @param[in] dofNumberAccessor accessor for the dofs numbers * @param[in] compFlowAccessor accessor for wrappers registered by the solver @@ -239,10 +237,10 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector + * @param[in] kernelFlags flags packed all together */ FaceBasedAssemblyKernel( integer const numPhases, globalIndex const rankOffset, - integer const hasCapPressure, STENCILWRAPPER const & stencilWrapper, DofNumberAccessor const & dofNumberAccessor, CompFlowAccessors const & compFlowAccessors, @@ -252,12 +250,12 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne CapPressureAccessors const & capPressureAccessors, PermeabilityAccessors const & permeabilityAccessors, ThermalConductivityAccessors const & thermalConductivityAccessors, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags > kernelFlags ) : Base( numPhases, rankOffset, - hasCapPressure, stencilWrapper, dofNumberAccessor, compFlowAccessors, @@ -266,7 +264,8 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne permeabilityAccessors, dt, localMatrix, - localRhs ), + localRhs, + kernelFlags ), m_temp( thermalCompFlowAccessors.get( fields::flow::temperature {} ) ), m_phaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::phaseEnthalpy {} ) ), m_dPhaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::dPhaseEnthalpy {} ) ), @@ -325,8 +324,8 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne localIndex const er_up, localIndex const esr_up, localIndex const ei_up, - real64 const & potGrad, - real64 const & phaseFlux, + real64 const potGrad, + real64 const phaseFlux, real64 const (&dPhaseFlux_dP)[2], real64 const (&dPhaseFlux_dC)[2][numComp] ) { @@ -370,7 +369,7 @@ class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKerne // Step 2.1: compute derivative of capillary pressure wrt temperature real64 dCapPressure_dT = 0.0; - if( m_hasCapPressure ) + if( AbstractBase::m_kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags::CapPressure ) ) { for( integer jp = 0; jp < m_numPhases; ++jp ) { @@ -627,10 +626,11 @@ class FaceBasedAssemblyKernelFactory globalIndex const rankOffset, string const & dofKey, integer const hasCapPressure, + integer const useTotalMassEquation, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -644,6 +644,12 @@ class FaceBasedAssemblyKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + BitFlags< isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags > kernelFlags; + if( hasCapPressure ) + kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags::CapPressure ); + if( useTotalMassEquation ) + kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags::TotalMassEquation ); + using KernelType = FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename KernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename KernelType::ThermalCompFlowAccessors thermalCompFlowAccessors( elemManager, solverName ); @@ -653,10 +659,10 @@ class FaceBasedAssemblyKernelFactory typename KernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); typename KernelType::ThermalConductivityAccessors thermalConductivityAccessors( elemManager, solverName ); - KernelType kernel( numPhases, rankOffset, hasCapPressure, stencilWrapper, dofNumberAccessor, + KernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor, compFlowAccessors, thermalCompFlowAccessors, multiFluidAccessors, thermalMultiFluidAccessors, capPressureAccessors, permeabilityAccessors, thermalConductivityAccessors, - dt, localMatrix, localRhs ); + dt, localMatrix, localRhs, kernelFlags ); KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } @@ -722,6 +728,7 @@ class DiffusionDispersionFaceBasedAssemblyKernel : * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector + * @param[in] kernelFlags flags packed together */ DiffusionDispersionFaceBasedAssemblyKernel( integer const numPhases, globalIndex const rankOffset, @@ -732,9 +739,10 @@ class DiffusionDispersionFaceBasedAssemblyKernel : DiffusionAccessors const & diffusionAccessors, DispersionAccessors const & dispersionAccessors, PorosityAccessors const & porosityAccessors, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags > kernelFlags ) : Base( numPhases, rankOffset, stencilWrapper, @@ -746,7 +754,8 @@ class DiffusionDispersionFaceBasedAssemblyKernel : porosityAccessors, dt, localMatrix, - localRhs ) + localRhs, + kernelFlags ) {} struct StackVariables : public Base::StackVariables @@ -791,8 +800,8 @@ class DiffusionDispersionFaceBasedAssemblyKernel : localIndex const er_up, localIndex const esr_up, localIndex const ei_up, - real64 const & compFracGrad, - real64 const & upwindCoefficient ) + real64 const compFracGrad, + real64 const upwindCoefficient ) { // We are in the loop over phases and components, ip provides the current phase index. @@ -863,7 +872,7 @@ class DiffusionDispersionFaceBasedAssemblyKernel : localIndex const er_up, localIndex const esr_up, localIndex const ei_up, - real64 const & compFracGrad ) + real64 const compFracGrad ) { // We are in the loop over phases and components, ip provides the current phase index. @@ -933,12 +942,13 @@ class DiffusionDispersionFaceBasedAssemblyKernelFactory integer const numPhases, globalIndex const rankOffset, string const & dofKey, - integer const & hasDiffusion, - integer const & hasDispersion, + integer const hasDiffusion, + integer const hasDispersion, + integer const useTotalMassEquation, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -952,6 +962,10 @@ class DiffusionDispersionFaceBasedAssemblyKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + BitFlags< isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags > kernelFlags; + if( useTotalMassEquation ) + kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags::TotalMassEquation ); + using kernelType = DiffusionDispersionFaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); @@ -962,7 +976,7 @@ class DiffusionDispersionFaceBasedAssemblyKernelFactory kernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor, compFlowAccessors, multiFluidAccessors, diffusionAccessors, dispersionAccessors, porosityAccessors, - dt, localMatrix, localRhs ); + dt, localMatrix, localRhs, kernelFlags ); kernelType::template launch< POLICY >( stencilWrapper.size(), hasDiffusion, hasDispersion, kernel ); @@ -1015,7 +1029,6 @@ class DirichletFaceBasedAssemblyKernel : public isothermalCompositionalMultiphas using Base::numComp; using Base::numDof; using Base::numEqn; - using Base::m_hasCapPressure; using Base::m_phaseMob; using Base::m_dPhaseMob; using Base::m_dPhaseMassDens; @@ -1045,7 +1058,6 @@ class DirichletFaceBasedAssemblyKernel : public isothermalCompositionalMultiphas * @brief Constructor for the kernel interface * @param[in] numPhases the number of fluid phases * @param[in] rankOffset the offset of my MPI rank - * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not * @param[in] faceManager the face manager * @param[in] stencilWrapper reference to the stencil wrapper * @param[in] fluidWrapper reference to the fluid wrapper @@ -1060,10 +1072,10 @@ class DirichletFaceBasedAssemblyKernel : public isothermalCompositionalMultiphas * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector + * @param[in] kernelFlags flags packed together */ DirichletFaceBasedAssemblyKernel( integer const numPhases, globalIndex const rankOffset, - integer const hasCapPressure, FaceManager const & faceManager, BoundaryStencilWrapper const & stencilWrapper, FLUIDWRAPPER const & fluidWrapper, @@ -1075,12 +1087,12 @@ class DirichletFaceBasedAssemblyKernel : public isothermalCompositionalMultiphas CapPressureAccessors const & capPressureAccessors, PermeabilityAccessors const & permeabilityAccessors, ThermalConductivityAccessors const & thermalConductivityAccessors, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + arrayView1d< real64 > const & localRhs, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags > kernelFlags ) : Base( numPhases, rankOffset, - hasCapPressure, faceManager, stencilWrapper, fluidWrapper, @@ -1091,7 +1103,8 @@ class DirichletFaceBasedAssemblyKernel : public isothermalCompositionalMultiphas permeabilityAccessors, dt, localMatrix, - localRhs ), + localRhs, + kernelFlags ), m_temp( thermalCompFlowAccessors.get( fields::flow::temperature {} ) ), m_phaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::phaseEnthalpy {} ) ), m_dPhaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::dPhaseEnthalpy {} ) ), @@ -1161,12 +1174,12 @@ class DirichletFaceBasedAssemblyKernel : public isothermalCompositionalMultiphas localIndex const esr, localIndex const ei, localIndex const kf, - real64 const & f, // potGrad times trans - real64 const & facePhaseMob, + real64 const f, // potGrad times trans + real64 const facePhaseMob, arraySlice1d< const real64, multifluid::USD_PHASE - 2 > const & facePhaseEnthalpy, arraySlice2d< const real64, multifluid::USD_PHASE_COMP-2 > const & facePhaseCompFrac, - real64 const & phaseFlux, - real64 const & dPhaseFlux_dP, + real64 const phaseFlux, + real64 const dPhaseFlux_dP, real64 const (&dPhaseFlux_dC)[numComp] ) { // We are in the loop over phases, ip provides the current phase index. @@ -1370,13 +1383,14 @@ class DirichletFaceBasedAssemblyKernelFactory createAndLaunch( integer const numComps, integer const numPhases, globalIndex const rankOffset, + integer const useTotalMassEquation, string const & dofKey, string const & solverName, FaceManager const & faceManager, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, MultiFluidBase & fluidBase, - real64 const & dt, + real64 const dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -1395,6 +1409,11 @@ class DirichletFaceBasedAssemblyKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + // for now, we neglect capillary pressure in the kernel + BitFlags< isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags > kernelFlags; + if( useTotalMassEquation ) + kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernelFlags::TotalMassEquation ); + using KernelType = DirichletFaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, typename FluidType::KernelWrapper >; typename KernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename KernelType::ThermalCompFlowAccessors thermalCompFlowAccessors( elemManager, solverName ); @@ -1404,13 +1423,10 @@ class DirichletFaceBasedAssemblyKernelFactory typename KernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); typename KernelType::ThermalConductivityAccessors thermalConductivityAccessors( elemManager, solverName ); - // for now, we neglect capillary pressure in the kernel - bool const hasCapPressure = false; - - KernelType kernel( numPhases, rankOffset, hasCapPressure, faceManager, stencilWrapper, fluidWrapper, + KernelType kernel( numPhases, rankOffset, faceManager, stencilWrapper, fluidWrapper, dofNumberAccessor, compFlowAccessors, thermalCompFlowAccessors, multiFluidAccessors, thermalMultiFluidAccessors, capPressureAccessors, permeabilityAccessors, thermalConductivityAccessors, - dt, localMatrix, localRhs ); + dt, localMatrix, localRhs, kernelFlags ); KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 1b5f63aff78..9e3c978cb50 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -60,6 +60,7 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, m_numPhases( 0 ), m_numComponents( 0 ), m_useMass( false ), + m_useTotalMassEquation( 1 ), m_maxCompFracChange( 1.0 ), m_maxRelativePresChange( 0.2 ), m_minScalingFactor( 0.01 ), @@ -71,6 +72,11 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Use mass formulation instead of molar" ); + this->registerWrapper( viewKeyStruct::useMassFlagString(), &m_useTotalMassEquation ). + setApplyDefaultValue( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Use total mass equation" ); + this->registerWrapper( viewKeyStruct::maxCompFracChangeString(), &m_maxCompFracChange ). setSizedFromParent( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -1036,6 +1042,7 @@ void CompositionalMultiphaseWell::assembleFluxTerms( real64 const GEOS_UNUSED_PA KernelLaunchSelector1< FluxKernel >( numFluidComponents(), subRegion.size(), dofManager.rankOffset(), + m_useTotalMassEquation, wellControls, wellElemDofNumber, nextWellElemIndex, @@ -1103,6 +1110,7 @@ void CompositionalMultiphaseWell::assembleAccumulationTerms( DomainPartition con subRegion.size(), numFluidPhases(), dofManager.rankOffset(), + m_useTotalMassEquation, wellElemDofNumber, wellElemGhostRank, wellElemVolume, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index 9154e28d0b7..b3bce32b232 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -276,6 +276,8 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr char const * useMassFlagString() { return CompositionalMultiphaseBase::viewKeyStruct::useMassFlagString(); } + static constexpr char const * useTotalMassEquationString() { return CompositionalMultiphaseBase::viewKeyStruct::useTotalMassEquationString(); } + static constexpr char const * relPermNamesString() { return CompositionalMultiphaseBase::viewKeyStruct::relPermNamesString(); } static constexpr char const * maxCompFracChangeString() { return CompositionalMultiphaseBase::viewKeyStruct::maxCompFracChangeString(); } @@ -367,6 +369,9 @@ class CompositionalMultiphaseWell : public WellSolverBase /// flag indicating whether mass or molar formulation should be used integer m_useMass; + /// flag indicating whether total mass equation should be used + integer m_useTotalMassEquation; + /// list of relative permeability model names per target region array1d< string > m_relPermModelNames; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp index 51b2f64e73c..70a8e51721d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp @@ -282,6 +282,7 @@ void FluxKernel:: launch( localIndex const size, globalIndex const rankOffset, + integer const useTotalMassEquation, WellControls const & wellControls, arrayView1d< globalIndex const > const & wellElemDofNumber, arrayView1d< localIndex const > const & nextWellElemIndex, @@ -416,11 +417,14 @@ FluxKernel:: oneSidedDofColIndices_dPresCompUp[jdof] = offsetUp + COFFSET::DPRES + jdof; } - // Apply equation/variable change transformation(s) - real64 work[NC+1]{}; - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, 1, oneSidedFluxJacobian_dRate, work ); - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC + 1, oneSidedFluxJacobian_dPresCompUp, work ); - shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, oneSidedFlux ); + if( useTotalMassEquation > 0 ) + { + // Apply equation/variable change transformation(s) + real64 work[NC + 1]{}; + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, 1, oneSidedFluxJacobian_dRate, work ); + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC + 1, oneSidedFluxJacobian_dPresCompUp, work ); + shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, oneSidedFlux ); + } for( integer i = 0; i < NC; ++i ) { @@ -479,11 +483,14 @@ FluxKernel:: dofColIndices_dPresCompUp[jdof] = offsetUp + COFFSET::DPRES + jdof; } - // Apply equation/variable change transformation(s) - real64 work[NC+1]{}; - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC, 1, 2, localFluxJacobian_dRate, work ); - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC, NC + 1, 2, localFluxJacobian_dPresCompUp, work ); - shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( NC, NC, 2, localFlux ); + if( useTotalMassEquation > 0 ) + { + // Apply equation/variable change transformation(s) + real64 work[NC + 1]{}; + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC, 1, 2, localFluxJacobian_dRate, work ); + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC, NC + 1, 2, localFluxJacobian_dPresCompUp, work ); + shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( NC, NC, 2, localFlux ); + } for( integer i = 0; i < 2*NC; ++i ) { @@ -509,6 +516,7 @@ FluxKernel:: void FluxKernel:: \ launch< NC >( localIndex const size, \ globalIndex const rankOffset, \ + integer const useTotalMassEquation, \ WellControls const & wellControls, \ arrayView1d< globalIndex const > const & wellElemDofNumber, \ arrayView1d< localIndex const > const & nextWellElemIndex, \ @@ -1282,6 +1290,7 @@ AccumulationKernel:: launch( localIndex const size, integer const numPhases, globalIndex const rankOffset, + integer const useTotalMassEquation, arrayView1d< globalIndex const > const & wellElemDofNumber, arrayView1d< integer const > const & wellElemGhostRank, arrayView1d< real64 const > const & wellElemVolume, @@ -1341,10 +1350,13 @@ AccumulationKernel:: dofColIndices[idof] = wellElemDofNumber[iwelem] + COFFSET::DPRES + idof; } - // Apply equation/variable change transformation(s) - real64 work[NC+1]; - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC + 1, localAccumJacobian, work ); - shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, localAccum ); + if( useTotalMassEquation > 0 ) + { + // Apply equation/variable change transformation(s) + real64 work[NC + 1]; + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC + 1, localAccumJacobian, work ); + shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, localAccum ); + } // add contribution to residual and jacobian for( integer ic = 0; ic < NC; ++ic ) @@ -1364,6 +1376,7 @@ AccumulationKernel:: launch< NC >( localIndex const size, \ integer const numPhases, \ globalIndex const rankOffset, \ + integer const useTotalMassEquation, \ arrayView1d< globalIndex const > const & wellElemDofNumber, \ arrayView1d< integer const > const & wellElemGhostRank, \ arrayView1d< real64 const > const & wellElemVolume, \ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index d8a8d9b4ba4..d18c512fed0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -163,6 +163,7 @@ struct FluxKernel static void launch( localIndex const size, globalIndex const rankOffset, + integer const useTotalMassEquation, WellControls const & wellControls, arrayView1d< globalIndex const > const & wellElemDofNumber, arrayView1d< localIndex const > const & nextWellElemIndex, @@ -361,6 +362,7 @@ struct AccumulationKernel launch( localIndex const size, integer const numPhases, globalIndex const rankOffset, + integer const useTotalMassEquation, arrayView1d< globalIndex const > const & wellElemDofNumber, arrayView1d< integer const > const & wellElemGhostRank, arrayView1d< real64 const > const & wellElemVolume, diff --git a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp index 3c080074745..480b1b2be92 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp @@ -349,6 +349,8 @@ assembleCouplingTerms( real64 const time_n, arrayView1d< localIndex const > const & resElementIndex = perforationData->getField< fields::perforation::reservoirElementIndex >(); + bool const useTotalMassEquation = this->flowSolver()->useTotalMassEquation() > 0; + RAJA::ReduceSum< parallelDeviceReduce, integer > numCrossflowPerforations( 0 ); // loop over the perforations and add the rates to the residual and jacobian @@ -411,10 +413,13 @@ assembleCouplingTerms( real64 const time_n, } } - // Apply equation/variable change transformation(s) - stackArray1d< real64, 2 * MAX_NUM_DOF > work( 2 * resNumDofs ); - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComps, numComps, resNumDofs*2, 2, localPerfJacobian, work ); - shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComps, numComps, 2, localPerf ); + if( useTotalMassEquation ) + { + // Apply equation/variable change transformation(s) + stackArray1d< real64, 2 * MAX_NUM_DOF > work( 2 * resNumDofs ); + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComps, numComps, resNumDofs * 2, 2, localPerfJacobian, work ); + shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComps, numComps, 2, localPerf ); + } for( localIndex i = 0; i < localPerf.size(); ++i ) { diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp index d07ab7d09ca..81881cbc936 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp @@ -178,6 +178,7 @@ void MultiphasePoromechanics::assembleSystem( real64 const GEOS_UNUSED_PARAM( ti flowDofKey, flowSolver()->numFluidComponents(), flowSolver()->numFluidPhases(), + flowSolver()->useTotalMassEquation(), FlowSolverBase::viewKeyStruct::fluidNamesString() ); } else @@ -194,6 +195,7 @@ void MultiphasePoromechanics::assembleSystem( real64 const GEOS_UNUSED_PARAM( ti flowDofKey, flowSolver()->numFluidComponents(), flowSolver()->numFluidPhases(), + flowSolver()->useTotalMassEquation(), FlowSolverBase::viewKeyStruct::fluidNamesString() ); } } ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp index fbaba63e303..d690038b45c 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp @@ -100,6 +100,7 @@ class MultiphasePoromechanics : string const inputFlowDofKey, localIndex const numComponents, localIndex const numPhases, + integer const useTotalMassEquation, string const fluidModelKey ); //***************************************************************************** @@ -343,6 +344,8 @@ class MultiphasePoromechanics : /// Number of phases localIndex const m_numPhases; + /// Use total mass equation flag + integer const m_useTotalMassEquation; }; using MultiphasePoromechanicsKernelFactory = @@ -356,6 +359,7 @@ using MultiphasePoromechanicsKernelFactory = string const, localIndex const, localIndex const, + integer const, string const >; /** diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp index 06ba0062589..441bd7c0a6d 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp @@ -52,6 +52,7 @@ MultiphasePoromechanics( NodeManager const & nodeManager, string const inputFlowDofKey, localIndex const numComponents, localIndex const numPhases, + integer const useTotalMassEquation, string const fluidModelKey ): Base( nodeManager, edgeManager, @@ -73,7 +74,8 @@ MultiphasePoromechanics( NodeManager const & nodeManager, m_dFluidPhaseVolFrac( elementSubRegion.template getField< fields::flow::dPhaseVolumeFraction >() ), m_dGlobalCompFraction_dGlobalCompDensity( elementSubRegion.template getField< fields::flow::dGlobalCompFraction_dGlobalCompDensity >() ), m_numComponents( numComponents ), - m_numPhases( numPhases ) + m_numPhases( numPhases ), + m_useTotalMassEquation( useTotalMassEquation ) { GEOS_ERROR_IF_GT_MSG( m_numComponents, maxNumComponents, "MultiphasePoromechanics solver allows at most " << @@ -641,12 +643,15 @@ complete( localIndex const k, integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint; constexpr integer maxNumDisplacementDofs = FE_TYPE::maxSupportPoints * numDofPerTestSupportPoint; - // Apply equation/variable change transformation(s) - real64 work[maxNumDisplacementDofs > ( maxNumComponents + 1 ) ? maxNumDisplacementDofs : maxNumComponents + 1]; - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( m_numComponents, numDisplacementDofs, stack.dLocalResidualMass_dDisplacement, work ); - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( m_numComponents, 1, stack.dLocalResidualMass_dPressure, work ); - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( m_numComponents, m_numComponents, stack.dLocalResidualMass_dComponents, work ); - shiftElementsAheadByOneAndReplaceFirstElementWithSum( m_numComponents, stack.localResidualMass ); + if( m_useTotalMassEquation > 0 ) + { + // Apply equation/variable change transformation(s) + real64 work[maxNumDisplacementDofs > ( maxNumComponents + 1 ) ? maxNumDisplacementDofs : maxNumComponents + 1]; + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( m_numComponents, numDisplacementDofs, stack.dLocalResidualMass_dDisplacement, work ); + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( m_numComponents, 1, stack.dLocalResidualMass_dPressure, work ); + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( m_numComponents, m_numComponents, stack.dLocalResidualMass_dComponents, work ); + shiftElementsAheadByOneAndReplaceFirstElementWithSum( m_numComponents, stack.localResidualMass ); + } for( int localNode = 0; localNode < numSupportPoints; ++localNode ) { diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp index 00e979187cb..5c60e6820da 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp @@ -75,6 +75,7 @@ class ThermalMultiphasePoromechanics : using Base::m_pressure_n; using Base::m_numComponents; using Base::m_numPhases; + using Base::m_useTotalMassEquation; using Base::m_solidDensity; using Base::m_fluidPhaseMassDensity; using Base::m_dFluidPhaseMassDensity; @@ -108,6 +109,7 @@ class ThermalMultiphasePoromechanics : string const inputFlowDofKey, localIndex const numComponents, localIndex const numPhases, + integer const useTotalMassEquation, string const fluidModelKey ); /** @@ -338,6 +340,7 @@ using ThermalMultiphasePoromechanicsKernelFactory = string const, localIndex const, localIndex const, + integer const, string const >; } // namespace thermalporomechanicsKernels diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp index 55fcb820bd5..85dd0ec4da6 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp @@ -51,6 +51,7 @@ ThermalMultiphasePoromechanics( NodeManager const & nodeManager, string const inputFlowDofKey, localIndex const numComponents, localIndex const numPhases, + integer const useTotalMassEquation, string const fluidModelKey ): Base( nodeManager, edgeManager, @@ -68,6 +69,7 @@ ThermalMultiphasePoromechanics( NodeManager const & nodeManager, inputFlowDofKey, numComponents, numPhases, + useTotalMassEquation, fluidModelKey ), m_rockInternalEnergy_n( inputConstitutiveType.getInternalEnergy_n() ), m_rockInternalEnergy( inputConstitutiveType.getInternalEnergy() ), @@ -542,9 +544,12 @@ complete( localIndex const k, m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack ); integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint; - // Apply equation/variable change transformation(s) - real64 work[maxNumComponents + 1]{}; - shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( m_numComponents, 1, stack.dLocalResidualMass_dTemperature, work ); + if( m_useTotalMassEquation > 0 ) + { + // Apply equation/variable change transformation(s) + real64 work[maxNumComponents + 1]{}; + shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( m_numComponents, 1, stack.dLocalResidualMass_dTemperature, work ); + } // Step 1: assemble the derivatives of linear momentum balance wrt temperature into the global matrix