From 736cc7ccf4192f5088b641ec1ddb9ce0df534099 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 20 Dec 2024 14:04:38 -0600 Subject: [PATCH] build fix, code simplification, input checks --- .../fluidFlow/CompositionalMultiphaseBase.cpp | 301 +++++++----------- .../fluidFlow/CompositionalMultiphaseBase.hpp | 45 +-- .../CompositionalMultiphaseBaseFields.hpp | 15 +- .../fluidFlow/CompositionalMultiphaseFVM.cpp | 218 ++++++------- .../fluidFlow/CompositionalMultiphaseFVM.hpp | 8 - .../CompositionalMultiphaseHybridFVM.cpp | 11 - .../CompositionalMultiphaseHybridFVM.hpp | 7 - .../ReactiveCompositionalMultiphaseOBL.cpp | 4 + 8 files changed, 241 insertions(+), 368 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 2aadc509e5..1c7ba200bd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -270,6 +270,11 @@ void CompositionalMultiphaseBase::postInputInitialization() GEOS_ERROR( GEOS_FMT( "{}: '{}' is currently not available for thermal simulations", getDataContext(), viewKeyStruct::useZFormulationFlagString() ) ); } + if( m_hasDiffusion || m_hasDispersion ) + { + GEOS_ERROR( GEOS_FMT( "{}: {} is currently not available for diffusion or dispersion", + getDataContext(), viewKeyStruct::useZFormulationFlagString() ) ); + } if( m_isJumpStabilized ) // useZFormulation is not yet compatible with pressure stabilization { GEOS_ERROR( GEOS_FMT( "{}: pressure stabilization is not yet supported by {}", @@ -434,12 +439,13 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< globalCompFraction_n >( getName() ). setDimLabels( 1, fluid.componentNames() ). reference().resizeDimension< 1 >( m_numComponents ); - if( m_isFixedStressPoromechanicsUpdate ) - { - subRegion.registerField< globalCompFraction_k >( getName() ). - setDimLabels( 1, fluid.componentNames() ). - reference().resizeDimension< 1 >( m_numComponents ); - } + // may be needed later for sequential poromechanics implementation + //if( m_isFixedStressPoromechanicsUpdate ) + //{ + // subRegion.registerField< globalCompFraction_k >( getName() ). + // setDimLabels( 1, fluid.componentNames() ). + // reference().resizeDimension< 1 >( m_numComponents ); + //} subRegion.registerField< globalCompFractionScalingFactor >( getName() ); } else @@ -721,10 +727,11 @@ real64 CompositionalMultiphaseBase::updatePhaseVolumeFraction( ObjectManagerBase string const & fluidName = dataGroup.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, fluidName ); - if( m_isThermal ) + if( m_useZFormulation ) { - return thermalCompositionalMultiphaseBaseKernels:: - PhaseVolumeFractionKernelFactory:: + // isothermal for now + return isothermalCompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionZFormulationKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dataGroup, @@ -732,10 +739,10 @@ real64 CompositionalMultiphaseBase::updatePhaseVolumeFraction( ObjectManagerBase } else { - if( m_useZFormulation ) + if( m_isThermal ) { - return isothermalCompositionalMultiphaseBaseKernels:: - PhaseVolumeFractionZFormulationKernelFactory:: + return thermalCompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dataGroup, @@ -832,49 +839,43 @@ void CompositionalMultiphaseBase::updateCompAmount( ElementSubRegionBase & subRe { GEOS_MARK_FUNCTION; + arrayView2d< real64, compflow::USD_COMP > const compAmount = subRegion.getField< fields::flow::compAmount >(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ); CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); - arrayView1d< real64 const > const volume = subRegion.getElementVolume(); - arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >(); - arrayView2d< real64, compflow::USD_COMP > const compAmount = subRegion.getField< fields::flow::compAmount >(); integer const numComp = m_numComponents; - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + if( m_useZFormulation ) { - for( integer ic = 0; ic < numComp; ++ic ) - { - compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compDens[ei][ic]; - } - } ); -} - -void CompositionalMultiphaseBase::updateCompAmountZFormulation( ElementSubRegionBase & subRegion ) const -{ - GEOS_MARK_FUNCTION; - - string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ); - CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); - arrayView1d< real64 const > const volume = subRegion.getElementVolume(); - arrayView2d< real64 const, compflow::USD_COMP > const compFrac = subRegion.getField< fields::flow::globalCompFraction >(); - arrayView2d< real64, compflow::USD_COMP > const compAmount = subRegion.getField< fields::flow::compAmount >(); - // access total density stored in the fluid - string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); - MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); - arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity(); - - integer const numComp = m_numComponents; + arrayView2d< real64 const, compflow::USD_COMP > const compFrac = subRegion.getField< fields::flow::globalCompFraction >(); + // access total density stored in the fluid + string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); + MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); + arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity(); - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + for( integer ic = 0; ic < numComp; ++ic ) + { + // m = phi*V*z_c*rho_T + compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compFrac[ei][ic] * totalDens[ei][0]; + } + } ); + } + else { - for( integer ic = 0; ic < numComp; ++ic ) + arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) { - // m = phi*V*z_c*rho_T - compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compFrac[ei][ic] * totalDens[ei][0]; - } - } ); + for( integer ic = 0; ic < numComp; ++ic ) + { + compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compDens[ei][ic]; + } + } ); + } } void CompositionalMultiphaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const @@ -930,19 +931,13 @@ real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & sub real64 maxDeltaPhaseVolFrac; - if( m_useZFormulation ) - { - // For p, z_c as the primary unknowns - updateFluidModel( subRegion ); // rho_T is now a function of p, z_c from volume balance - updateCompAmountZFormulation( subRegion ); - } - else + if( !m_useZFormulation ) { // For p, rho_c as the primary unknowns updateGlobalComponentFraction( subRegion ); - updateFluidModel( subRegion ); - updateCompAmount( subRegion ); } + updateFluidModel( subRegion ); + updateCompAmount( subRegion ); maxDeltaPhaseVolFrac = updatePhaseVolumeFraction( subRegion ); updateRelPermModel( subRegion ); updatePhaseMobility( subRegion ); @@ -988,16 +983,16 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, compDens[ei][ic] = totalDens[ei][0] * compFrac[ei][ic]; } } ); - } - // with initial component densities defined - check if they need to be corrected to avoid zero diags etc - if( m_allowCompDensChopping ) - { - chopNegativeDensities( subRegion ); + // with initial component densities defined - check if they need to be corrected to avoid zero diags etc + if( m_allowCompDensChopping ) + { + chopNegativeDensities( subRegion ); + } } } ); - // with initial component densities defined - check if they need to be corrected to avoid zero diags etc + // check if comp fractions need to be corrected to avoid zero diags etc if( m_useZFormulation ) { DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); @@ -1012,14 +1007,7 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, { // Initialize/update dependent state quantities - if( m_useZFormulation ) - { - updateCompAmountZFormulation( subRegion ); - } - else - { - updateCompAmount( subRegion ); - } + updateCompAmount( subRegion ); updatePhaseVolumeFraction( subRegion ); // Update the constitutive models that only depend on @@ -1388,16 +1376,8 @@ void CompositionalMultiphaseBase::initializePostInitialConditionsPreSubGroups() arrayView1d< string const > const & regionNames ) { FieldIdentifiers fieldsToBeSync; - if( m_useZFormulation ) - { - fieldsToBeSync.addElementFields( { fields::flow::globalCompFraction::key() }, - regionNames ); - } - else - { - fieldsToBeSync.addElementFields( { fields::flow::globalCompDensity::key() }, - regionNames ); - } + fieldsToBeSync.addElementFields( {m_useZFormulation ? fields::flow::globalCompFraction::key() : fields::flow::globalCompDensity::key() }, + regionNames ); CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); @@ -1461,50 +1441,34 @@ void CompositionalMultiphaseBase::assembleSystem( real64 const GEOS_UNUSED_PARAM { GEOS_MARK_FUNCTION; - // Accumulation term - if( m_useZFormulation ) + // Accumulation and other local terms + assembleLocalTerms( domain, + dofManager, + localMatrix, + localRhs ); + + if( m_isJumpStabilized ) { - assembleZFormulationAccumulation( domain, - dofManager, - localMatrix, - localRhs ); - - assembleZFormulationFluxTerms( dt, - domain, - dofManager, - localMatrix, - localRhs ); + assembleStabilizedFluxTerms( dt, + domain, + dofManager, + localMatrix, + localRhs ); } else { - assembleAccumulationAndVolumeBalanceTerms( domain, - dofManager, - localMatrix, - localRhs ); - - if( m_isJumpStabilized ) - { - assembleStabilizedFluxTerms( dt, - domain, - dofManager, - localMatrix, - localRhs ); - } - else - { - assembleFluxTerms( dt, - domain, - dofManager, - localMatrix, - localRhs ); - } + assembleFluxTerms( dt, + domain, + dofManager, + localMatrix, + localRhs ); } } -void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const +void CompositionalMultiphaseBase::assembleLocalTerms( DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const { GEOS_MARK_FUNCTION; @@ -1523,14 +1487,16 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); CoupledSolidBase const & solid = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - if( m_isThermal ) + if( m_useZFormulation ) { - thermalCompositionalMultiphaseBaseKernels:: - AccumulationKernelFactory:: + // isothermal for now + isothermalCompositionalMultiphaseBaseKernels:: + AccumulationZFormulationKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), m_useTotalMassEquation, + m_useSimpleAccumulation, dofKey, subRegion, fluid, @@ -1540,64 +1506,42 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom } else { - isothermalCompositionalMultiphaseBaseKernels:: - AccumulationKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - m_useTotalMassEquation, - m_useSimpleAccumulation, - dofKey, - subRegion, - fluid, - solid, - localMatrix, - localRhs ); + if( m_isThermal ) + { + thermalCompositionalMultiphaseBaseKernels:: + AccumulationKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + m_useTotalMassEquation, + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); + } + else + { + isothermalCompositionalMultiphaseBaseKernels:: + AccumulationKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + m_useTotalMassEquation, + m_useSimpleAccumulation, + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); + } } } ); } ); } -void CompositionalMultiphaseBase::assembleZFormulationAccumulation( DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const -{ - GEOS_MARK_FUNCTION; - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel const & mesh, - arrayView1d< string const > const & regionNames ) - { - mesh.getElemManager().forElementSubRegions( regionNames, - [&]( localIndex const, - ElementSubRegionBase const & subRegion ) - { - string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); - string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); - string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() ); - - MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); - CoupledSolidBase const & solid = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - - // isothermal for now - isothermalCompositionalMultiphaseBaseKernels:: - AccumulationZFormulationKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - m_useTotalMassEquation, - m_useSimpleAccumulation, - dofKey, - subRegion, - fluid, - solid, - localMatrix, - localRhs ); - } ); - } ); -} - void CompositionalMultiphaseBase::applyBoundaryConditions( real64 const time_n, real64 const dt, DomainPartition & domain, @@ -2926,12 +2870,13 @@ void CompositionalMultiphaseBase::saveConvergedState( ElementSubRegionBase & sub arrayView2d< real64, compflow::USD_COMP > const & compFrac_n = subRegion.template getField< fields::flow::globalCompFraction_n >(); compFrac_n.setValues< parallelDevicePolicy<> >( compFrac ); - if( m_isFixedStressPoromechanicsUpdate ) - { - arrayView2d< real64, compflow::USD_COMP > const & compFrac_k = - subRegion.template getField< fields::flow::globalCompFraction_k >(); - compFrac_k.setValues< parallelDevicePolicy<> >( compFrac ); - } + // may be useful later for sequential poromechanics implementation + //if( m_isFixedStressPoromechanicsUpdate ) + //{ + // arrayView2d< real64, compflow::USD_COMP > const & compFrac_k = + // subRegion.template getField< fields::flow::globalCompFraction_k >(); + // compFrac_k.setValues< parallelDevicePolicy<> >( compFrac ); + //} } else { @@ -2948,8 +2893,6 @@ void CompositionalMultiphaseBase::saveConvergedState( ElementSubRegionBase & sub } } - - arrayView2d< real64 const, compflow::USD_COMP > const & compAmount = subRegion.template getField< fields::flow::compAmount >(); arrayView2d< real64, compflow::USD_COMP > const & compAmount_n = diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 0503eaad4c..174e2c856d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -140,12 +140,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ void updateCompAmount( ElementSubRegionBase & subRegion ) const; - /** - * @brief Update components mass/moles if Z formulation is used - * @param subRegion the subregion storing the required fields - */ - void updateCompAmountZFormulation( ElementSubRegionBase & subRegion ) const; - /** * @brief Update energy * @param subRegion the subregion storing the required fields @@ -197,7 +191,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase { return m_useMass ? units::Unit::Mass : units::Unit::Mole; } /** - * @brief assembles the accumulation and volume balance terms for all cells + * @brief assembles the accumulation other local terms for all cells * @param time_n previous time value * @param dt time step * @param domain the physical domain object @@ -205,10 +199,10 @@ class CompositionalMultiphaseBase : public FlowSolverBase * @param localMatrix the system matrix * @param localRhs the system right-hand side vector */ - void assembleAccumulationAndVolumeBalanceTerms( DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; + void assembleLocalTerms( DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; /** * @brief assembles the flux terms for all cells @@ -241,35 +235,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const = 0; - /** - * @brief assembles the accumulation term (Z formulation) for all cells - * @param time_n previous time value - * @param dt time step - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localMatrix the system matrix - * @param localRhs the system right-hand side vector - */ - void assembleZFormulationAccumulation( DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; - - /** - * @brief assembles the flux terms (Z formulation) for all cells - * @param time_n previous time value - * @param dt time step - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param matrix the system matrix - * @param rhs the system right-hand side vector - */ - virtual void - assembleZFormulationFluxTerms( real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const = 0; /**@}*/ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp index d3ff384165..3f520747ae 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp @@ -80,13 +80,14 @@ DECLARE_FIELD( globalCompFraction_n, NO_WRITE, "Global component fraction at the previous converged time step" ); -DECLARE_FIELD( globalCompFraction_k, - "globalCompFraction_k", - array2dLayoutComp, - 0, - NOPLOT, - NO_WRITE, - "Global component fraction updates at the previous sequential iteration" ); +// may be needed later for sequential poromechanics implementation +//DECLARE_FIELD( globalCompFraction_k, +// "globalCompFraction_k", +// array2dLayoutComp, +// 0, +// NOPLOT, +// NO_WRITE, +// "Global component fraction updates at the previous sequential iteration" ); DECLARE_FIELD( faceGlobalCompFraction, "faceGlobalCompFraction", diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 84fed92493..e5636f6920 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -111,9 +111,18 @@ void CompositionalMultiphaseFVM::postInputInitialization() { CompositionalMultiphaseBase::postInputInitialization(); - if( m_scalingType == ScalingType::Local && m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None ) + if( m_scalingType == ScalingType::Local && + m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None ) { - GEOS_ERROR( GEOS_FMT( "{}: line search is not supported for {} = {}", getName(), viewKeyStruct::scalingTypeString(), EnumStrings< ScalingType >::toString( ScalingType::Local )) ); + GEOS_ERROR( GEOS_FMT( "{}: line search is not supported for {} = {}", + getName(), viewKeyStruct::scalingTypeString(), + EnumStrings< ScalingType >::toString( ScalingType::Local )) ); + } + + if( m_useZFormulation && m_dbcParams.useDBC ) // useZFormulation is not compatible with DBC + { + GEOS_ERROR( GEOS_FMT( "{}: '{}' is not compatible with {}", + getDataContext(), viewKeyStruct::useZFormulationFlagString(), viewKeyStruct::useDBCString() ) ); } } @@ -184,16 +193,20 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); // Convective flux - if( m_isThermal ) + + if( m_useZFormulation ) { - thermalCompositionalMultiphaseFVMKernels:: - FluxComputeKernelFactory:: + // isothermal only for now + isothermalCompositionalMultiphaseFVMKernels:: + FluxComputeZFormulationKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), elemDofKey, m_hasCapPressure, m_useTotalMassEquation, + m_useNewGravity, + fluxApprox.upwindingParams(), getName(), mesh.getElemManager(), stencilWrapper, @@ -203,9 +216,9 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, } else { - if( m_dbcParams.useDBC ) + if( m_isThermal ) { - dissipationCompositionalMultiphaseFVMKernels:: + thermalCompositionalMultiphaseFVMKernels:: FluxComputeKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, @@ -218,74 +231,95 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, stencilWrapper, dt, localMatrix.toViewConstSizes(), - localRhs.toView(), - m_dbcParams.omega, - getNonlinearSolverParameters().m_numNewtonIterations, - m_dbcParams.continuation, - m_dbcParams.miscible, - m_dbcParams.kappamin, - m_dbcParams.contMultiplier ); + localRhs.toView() ); } else { - isothermalCompositionalMultiphaseFVMKernels:: - FluxComputeKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - elemDofKey, - m_hasCapPressure, - m_useTotalMassEquation, - m_useNewGravity, - fluxApprox.upwindingParams(), - getName(), - mesh.getElemManager(), - stencilWrapper, - dt, - localMatrix.toViewConstSizes(), - localRhs.toView() ); + if( m_dbcParams.useDBC ) + { + dissipationCompositionalMultiphaseFVMKernels:: + FluxComputeKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + m_hasCapPressure, + m_useTotalMassEquation, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView(), + m_dbcParams.omega, + getNonlinearSolverParameters().m_numNewtonIterations, + m_dbcParams.continuation, + m_dbcParams.miscible, + m_dbcParams.kappamin, + m_dbcParams.contMultiplier ); + } + else + { + isothermalCompositionalMultiphaseFVMKernels:: + FluxComputeKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + m_hasCapPressure, + m_useTotalMassEquation, + m_useNewGravity, + fluxApprox.upwindingParams(), + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } } - } - // Diffusive and dispersive flux - if( m_hasDiffusion || m_hasDispersion ) - { + // Diffusive and dispersive flux - if( m_isThermal ) + if( m_hasDiffusion || m_hasDispersion ) { - thermalCompositionalMultiphaseFVMKernels:: - DiffusionDispersionFluxComputeKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - elemDofKey, - m_hasDiffusion, - m_hasDispersion, - m_useTotalMassEquation, - getName(), - mesh.getElemManager(), - stencilWrapper, - dt, - localMatrix.toViewConstSizes(), - localRhs.toView() ); - } - else - { - isothermalCompositionalMultiphaseFVMKernels:: - DiffusionDispersionFluxComputeKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - elemDofKey, - m_hasDiffusion, - m_hasDispersion, - m_useTotalMassEquation, - getName(), - mesh.getElemManager(), - stencilWrapper, - dt, - localMatrix.toViewConstSizes(), - localRhs.toView() ); + + if( m_isThermal ) + { + thermalCompositionalMultiphaseFVMKernels:: + DiffusionDispersionFluxComputeKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + m_hasDiffusion, + m_hasDispersion, + m_useTotalMassEquation, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } + else + { + isothermalCompositionalMultiphaseFVMKernels:: + DiffusionDispersionFluxComputeKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + m_hasDiffusion, + m_hasDispersion, + m_useTotalMassEquation, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } } } @@ -293,54 +327,6 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, } ); } -void CompositionalMultiphaseFVM::assembleZFormulationFluxTerms( real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const -{ - GEOS_MARK_FUNCTION; - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel const & mesh, - arrayView1d< string const > const & ) - { - NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); - FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); - FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName ); - - string const & elemDofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); - - fluxApprox.forAllStencils( mesh, [&] ( auto & stencil ) - { - typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); - - GEOS_ERROR_IF( m_isThermal, GEOS_FMT( - "{}: Z Formulation is currently not available for thermal simulations", getDataContext() ) ); - GEOS_ERROR_IF( m_hasDiffusion || m_hasDispersion, GEOS_FMT( - "{}: Z Formulation is currently not available for Diffusion or Dispersion", getDataContext() ) ); - - // isothermal only for now - isothermalCompositionalMultiphaseFVMKernels:: - FluxComputeZFormulationKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - elemDofKey, - m_hasCapPressure, - m_useTotalMassEquation, - m_useNewGravity, - fluxApprox.upwindingParams(), - getName(), - mesh.getElemManager(), - stencilWrapper, - dt, - localMatrix.toViewConstSizes(), - localRhs.toView() ); - } ); - } ); -} - void CompositionalMultiphaseFVM::assembleStabilizedFluxTerms( real64 const dt, DomainPartition const & domain, DofManager const & dofManager, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp index 7e77160184..2798752e69 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp @@ -143,14 +143,6 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const override; - virtual void - assembleZFormulationFluxTerms( real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const override; - - virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index a4e379681e..909997e7d6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -419,17 +419,6 @@ void CompositionalMultiphaseHybridFVM::assembleFluxTerms( real64 const dt, } ); } -void CompositionalMultiphaseHybridFVM::assembleZFormulationFluxTerms( real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const -{ - // z formulation not implemented - GEOS_UNUSED_VAR( dt, domain, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "Z formulation not yet available for this flow solver" ); -} - void CompositionalMultiphaseHybridFVM::assembleStabilizedFluxTerms( real64 const dt, DomainPartition const & domain, DofManager const & dofManager, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp index bb85c071ae..5fa78d20ce 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp @@ -139,13 +139,6 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const override; - virtual void - assembleZFormulationFluxTerms( real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const override; - virtual void assembleStabilizedFluxTerms( real64 const dt, DomainPartition const & domain, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp index 8d28b7bfcd..07c38cd272 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp @@ -272,6 +272,10 @@ void ReactiveCompositionalMultiphaseOBL::registerDataOnMesh( Group & meshBodies subRegion.registerField< bcGlobalCompFraction >( solverName ). reference().resizeDimension< 1 >( m_numComponents ); + // we need to register this fiels in any case (if there is a single component or not) + // to be able to pass the view to OBLOperatorsKernel + subRegion.registerField< globalCompFraction_n >( solverName ). + reference().resizeDimension< 1 >( m_numComponents ); // in principle, referencePorosity could be used directly from solid model, // but was duplicated to remove dependency on solid subRegion.registerField< referencePorosity >( solverName );