diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 6dfbf981b0c..746f223b649 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3479-9362-cffefcc + baseline: integratedTests/baseline_integratedTests-pr3495-9552-257c87e allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 0bc5c80f561..c5617ba8ed9 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,18 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3495 (2024-01-08) +===================== +Add missing logic to support switching from fixed mass rate injection rate constraint to max injection pressure. + +PR #3384 (2024-01-07) +===================== +Added plastic strain output. + +PR #3486 (2025-01-06) +===================== +useNewGravity became gravityDensityScheme. + PR #3479 (2024-12-15) ===================== Refine inputFiles/compositionalMultiphaseFlow: shift reference pressures to initial pressures, make nonlinear tuning more reasonable, minimize output. diff --git a/inputFiles/compositionalMultiphaseWell/compositionalMultiphaseWell.ats b/inputFiles/compositionalMultiphaseWell/compositionalMultiphaseWell.ats index ea1d810c064..16de5c4d90c 100644 --- a/inputFiles/compositionalMultiphaseWell/compositionalMultiphaseWell.ats +++ b/inputFiles/compositionalMultiphaseWell/compositionalMultiphaseWell.ats @@ -76,6 +76,22 @@ decks = [ partitions=[(1, 1, 1), (2, 2, 1)], restart_step=12, check_step=25, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="isothm_mass_inj_table", + description= + "Isothermal mass injection constraint testing control switching", + partitions=[(1, 1, 1), (1, 2, 1)], + restart_step=24, + check_step=28, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="isothm_vol_inj_table", + description= + "Isothermal vol injection constraint testing control switching, should have same results as isothm_mass_inj_table.xml", + partitions=[(1, 1, 1), (1, 2, 1)], + restart_step=24, + check_step=28, restartcheck_params=RestartcheckParameters(**restartcheck_params)) ] diff --git a/inputFiles/compositionalMultiphaseWell/isothm_mass_inj_table.xml b/inputFiles/compositionalMultiphaseWell/isothm_mass_inj_table.xml new file mode 100644 index 00000000000..14076fab917 --- /dev/null +++ b/inputFiles/compositionalMultiphaseWell/isothm_mass_inj_table.xml @@ -0,0 +1,229 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/compositionalMultiphaseWell/isothm_vol_inj_table.xml b/inputFiles/compositionalMultiphaseWell/isothm_vol_inj_table.xml new file mode 100644 index 00000000000..d39421e36fc --- /dev/null +++ b/inputFiles/compositionalMultiphaseWell/isothm_vol_inj_table.xml @@ -0,0 +1,234 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index bed64cb6e6a..97573fce5f3 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -135,6 +135,11 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates localIndex const q, real64 ( &elasticStrain )[6] ) const override final; + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + GEOS_HOST_DEVICE virtual real64 getBulkModulus( localIndex const k ) const override final { @@ -171,6 +176,13 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates protected: + GEOS_HOST_DEVICE + virtual void computeElasticStrain( localIndex const k, + localIndex const q, + real64 const ( &stress )[6], + real64 ( &elasticStrainInc )[6] ) const; + + /// A reference to the ArrayView holding the bulk modulus for each element. arrayView1d< real64 const > const m_bulkModulus; @@ -210,24 +222,53 @@ void ElasticIsotropicUpdates::getElasticStiffness( localIndex const k, } +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::computeElasticStrain( localIndex const k, + localIndex const q, + real64 const ( &stress )[6], + real64 ( & elasticStrain)[6] ) const +{ + GEOS_UNUSED_VAR( q ); + real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] ); + real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] ); + + elasticStrain[0] = ( stress[0] - nu*stress[1] - nu*stress[2])/E; + elasticStrain[1] = (-nu*stress[0] + stress[1] - nu*stress[2])/E; + elasticStrain[2] = (-nu*stress[0] - nu*stress[1] + stress[2])/E; + + elasticStrain[3] = stress[3] / m_shearModulus[k]; + elasticStrain[4] = stress[4] / m_shearModulus[k]; + elasticStrain[5] = stress[5] / m_shearModulus[k]; +} + GEOS_HOST_DEVICE inline void ElasticIsotropicUpdates::getElasticStrain( localIndex const k, localIndex const q, real64 ( & elasticStrain)[6] ) const { - real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] ); - real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] ); - elasticStrain[0] = ( m_newStress[k][q][0] - nu*m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E; - elasticStrain[1] = (-nu*m_newStress[k][q][0] + m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E; - elasticStrain[2] = (-nu*m_newStress[k][q][0] - nu*m_newStress[k][q][1] + m_newStress[k][q][2])/E; + real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrain ); - elasticStrain[3] = m_newStress[k][q][3] / m_shearModulus[k]; - elasticStrain[4] = m_newStress[k][q][4] / m_shearModulus[k]; - elasticStrain[5] = m_newStress[k][q][5] / m_shearModulus[k]; } +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const +{ + + real64 stress[6] = + {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrainInc ); + +} GEOS_HOST_DEVICE inline diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 4181e8e2013..01f36bf937f 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -115,6 +115,11 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates localIndex const q, real64 ( &elasticStrain )[6] ) const override final; + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + GEOS_HOST_DEVICE virtual void viscousStateUpdate( localIndex const k, localIndex const q, @@ -223,6 +228,60 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrain( localIndex cons } +GEOS_HOST_DEVICE +inline +void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const +{ + real64 const mu = m_shearModulus[k]; + real64 const p0 = m_refPressure; + real64 const eps_v0 = m_refStrainVol; + real64 const Cr = m_recompressionIndex[k]; + real64 deviator[6]{}; + real64 stress[6]{}; + real64 P; + real64 Q; + + LvArray::tensorOps::copy< 6 >( stress, m_newStress[k][q] ); + + twoInvariant::stressDecomposition( stress, + P, + Q, + deviator ); + + real64 elasticStrainVol = LvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0; + real64 elasticStrainDev = Q/3./mu; + + twoInvariant::strainRecomposition( elasticStrainVol, + elasticStrainDev, + deviator, + elasticStrainInc ); + + real64 oldStrain[6]{}; + + LvArray::tensorOps::copy< 6 >( stress, m_oldStress[k][q] ); + + twoInvariant::stressDecomposition( stress, + P, + Q, + deviator ); + + elasticStrainVol = LvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0; + elasticStrainDev = Q/3./mu; + + twoInvariant::strainRecomposition( elasticStrainVol, + elasticStrainDev, + deviator, + oldStrain ); + + for( localIndex i = 0; i<6; ++i ) + { + elasticStrainInc[i] -= oldStrain[i]; + } +} + + GEOS_HOST_DEVICE inline void ElasticIsotropicPressureDependentUpdates::smallStrainUpdate( localIndex const k, diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index 4445b88cc12..49179c5a95a 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -147,6 +147,16 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates real64 ( &stress )[6], DiscretizationOps & stiffness ) const final; + GEOS_HOST_DEVICE + virtual void getElasticStrain( localIndex const k, + localIndex const q, + real64 ( &elasticStrain )[6] ) const override final; + + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + // miscellaneous getters GEOS_HOST_DEVICE @@ -162,6 +172,13 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates return LvArray::math::max( LvArray::math::max( m_c44[k], m_c55[k] ), m_c66[k] ); } +protected: + GEOS_HOST_DEVICE + virtual void computeElasticStrain( localIndex const k, + localIndex const q, + real64 const (&stress)[6], + real64 ( &elasticStrain )[6] ) const; + private: /// A reference to the ArrayView holding c11 for each element. arrayView1d< real64 const > const m_c11; @@ -219,6 +236,64 @@ void ElasticOrthotropicUpdates::getElasticStiffness( localIndex const k, stiffness[5][5] = m_c66[k]; } + +GEOS_HOST_DEVICE +inline +void ElasticOrthotropicUpdates::computeElasticStrain( localIndex const k, + localIndex const q, + real64 const (&stress)[6], + real64 ( & elasticStrain)[6] ) const +{ + + GEOS_UNUSED_VAR( q ); + + real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]); + + elasticStrain[0] = + ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*stress[0] + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*stress[1] + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*stress[2] ) / + detC; + elasticStrain[1] = + ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*stress[0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[1] + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*stress[2] ) / + detC; + elasticStrain[2] = + ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*stress[0] + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*stress[1] + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*stress[2] ) / + detC; + + elasticStrain[3] = stress[3] / m_c44[k]; + elasticStrain[4] = stress[4] / m_c55[k]; + elasticStrain[5] = stress[5] / m_c66[k]; +} + + +GEOS_HOST_DEVICE +inline +void ElasticOrthotropicUpdates::getElasticStrain( localIndex const k, + localIndex const q, + real64 ( & elasticStrain)[6] ) const +{ + + real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrain ); + +} + +GEOS_HOST_DEVICE +inline +void ElasticOrthotropicUpdates::getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const +{ + + real64 stress[6] = + {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrainInc ); + +} + + inline GEOS_HOST_DEVICE void ElasticOrthotropicUpdates::smallStrainNoStateUpdate_StressOnly( localIndex const k, diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index f2c86d4439a..17fabe5479e 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -144,6 +144,16 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates GEOS_HOST_DEVICE virtual void getElasticStiffness( localIndex const k, localIndex const q, real64 ( &stiffness )[6][6] ) const override final; + GEOS_HOST_DEVICE + virtual void getElasticStrain( localIndex const k, + localIndex const q, + real64 ( &elasticStrain )[6] ) const override final; + + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + /** * @brief Getter for apparent shear modulus. * @return reference to shear modulus that will be used for computing stabilization scalling parameter. @@ -155,6 +165,15 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates } +protected: + + GEOS_HOST_DEVICE + virtual void computeElasticStrain( localIndex const k, + localIndex const q, + real64 const (&stress)[6], + real64 ( &elasticStrain )[6] ) const; + + private: /// A reference to the ArrayView holding c11 for each element. @@ -200,6 +219,57 @@ void ElasticTransverseIsotropicUpdates::getElasticStiffness( localIndex const k, stiffness[5][5] = m_c66[k]; } +GEOS_HOST_DEVICE +inline +void ElasticTransverseIsotropicUpdates::computeElasticStrain( localIndex const k, + localIndex const q, + real64 const (&stress)[6], + real64 ( & elasticStrain)[6] ) const +{ + GEOS_UNUSED_VAR( q ); + real64 const c12 = ( m_c11[k] - 2.0 * m_c66[k] ); + real64 const detC = m_c11[k]*(m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k]) - c12*(c12*m_c33[k] - m_c13[k]*m_c13[k]) + m_c13[k]*(c12*m_c13[k] - m_c11[k]*m_c13[k]); + + elasticStrain[0] = + ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[0] + (m_c13[k]*m_c13[k] - c12*m_c33[k])*stress[1] + (c12*m_c13[k] - m_c13[k]*m_c11[k])*stress[2] ) / detC; + elasticStrain[1] = + ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*stress[0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[1] + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*stress[2] ) / detC; + elasticStrain[2] = ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*stress[0] + (c12*m_c13[k] - m_c11[k]*m_c13[k])*stress[1] + (m_c11[k]*m_c11[k] - c12*c12)*stress[2] ) / detC; + + elasticStrain[3] = stress[3] / m_c44[k]; + elasticStrain[4] = stress[4] / m_c44[k]; + elasticStrain[5] = stress[5] / m_c66[k]; +} + + +GEOS_HOST_DEVICE +inline +void ElasticTransverseIsotropicUpdates::getElasticStrain( localIndex const k, + localIndex const q, + real64 ( & elasticStrain)[6] ) const +{ + + real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrain ); + +} + +GEOS_HOST_DEVICE +inline +void ElasticTransverseIsotropicUpdates::getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const +{ + + real64 stress[6] = + {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrainInc ); + +} + inline GEOS_HOST_DEVICE void ElasticTransverseIsotropicUpdates::smallStrainNoStateUpdate_StressOnly( localIndex const k, diff --git a/src/coreComponents/constitutive/solid/SolidBase.hpp b/src/coreComponents/constitutive/solid/SolidBase.hpp index 24b05eea71a..b56bd38b53a 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.hpp +++ b/src/coreComponents/constitutive/solid/SolidBase.hpp @@ -337,6 +337,24 @@ class SolidBaseUpdates GEOS_ERROR( "getElasticStrain() not implemented for this model" ); } + /** + * @brief Return the current elastic strain increment at a given material point (small-strain interface) + * + * @param k the element inex + * @param q the quadrature index + * @param elasticStrainInc Current elastic strain increment + */ + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc )[6] ) const + { + GEOS_UNUSED_VAR( k ); + GEOS_UNUSED_VAR( q ); + GEOS_UNUSED_VAR( elasticStrainInc ); + GEOS_ERROR( "getElasticStrainInc() of SolidBase was called." ); + } + /** * @brief Perform a viscous (rate-dependent) state update * diff --git a/src/coreComponents/finiteVolume/docs/FiniteVolume.rst.txt b/src/coreComponents/finiteVolume/docs/FiniteVolume.rst.txt index 755e215d64e..ad65ea50dd4 100644 --- a/src/coreComponents/finiteVolume/docs/FiniteVolume.rst.txt +++ b/src/coreComponents/finiteVolume/docs/FiniteVolume.rst.txt @@ -16,17 +16,27 @@ The numerical flux is obtained using the following expression for the mass flux .. math:: F_{KL} = \Upsilon_{KL} \frac{\rho^{upw}}{\mu^{upw}} \big( p_K - p_L - \rho^{avg} g ( d_K - d_L ) \big), -where :math:`p_K` is the pressure of cell :math:`K`, :math:`d_K` is the depth of cell :math:`K`, and :math:`\Upsilon_{KL}` is the standard TPFA transmissibility coefficient at the interface. + +where :math:`p_K` is the pressure of cell :math:`K`, :math:`\rho^{avg}` is the average fluid density, :math:`d_K` is the depth of cell :math:`K`, and :math:`\Upsilon_{KL}` is the standard TPFA transmissibility coefficient at the interface. The fluid density, :math:`\rho^{upw}`, and the fluid viscosity, :math:`\mu^{upw}`, are upwinded using the sign of the potential difference at the interface. -This is currently the only available discretization in the :ref:`CompositionalMultiphaseFlow`. +For :ref:`CompositionalMultiphaseFlow` there are two options to compute the average density, :math:`\rho^{avg}`. The desired option can be selected using the `gravityDensityScheme` parameter: + +#. `ArithmeticAverage`: :math:`\rho^{avg}` is computed using simple arithmetic average: :math:`\rho^{avg} = 0.5 \cdot (rho_K + rho_L)`, where :math:`rho_K` and :math:`rho_K` are densities in the two cells. + +#. `PhasePresence`: average phase density is computed using checking for phase presence: + + * :math:`\rho^{avg} = 0.5 \cdot (\rho_K + \rho_L)` if phase is present in both cells :math:`K` and :math:`L` + + * :math:`\rho^{avg} = \rho_K` if phase is present only in cell :math:`K` + + * :math:`\rho^{avg} = \rho_L` if phase is present only in cell :math:`L` Hybrid FVM ~~~~~~~~~~ This discretization scheme overcomes the limitations of the standard TPFA on non K-orthogonal meshes. The hybrid finite-volume scheme--equivalent to the well-known hybrid Mimetic Finite Difference (MFD) scheme--remains consistent with the pressure equation even when the mesh does not satisfy the K-orthogonality condition. -This numerical scheme is currently implemented in the `SinglePhaseHybridFVM` solver. The hybrid FVM scheme uses both cell-centered and face-centered pressure degrees of freedom. The one-sided face flux, :math:`F_{K,f}`, at face :math:`f` of cell :math:`K` is computed as: @@ -60,5 +70,3 @@ For a given interior face :math:`f` between two neighboring cells :math:`K` and We obtain a numerical scheme with :math:`n_{\textit{cells}}` cell-centered degrees of freedom and :math:`n_{\textit{faces}}` face-centered pressure degrees of freedom. The system involves :math:`n_{\textit{cells}}` mass conservation equations and :math:`n_{\textit{faces}}` face-based constraints. The linear systems can be efficiently solved using the MultiGrid Reduction (MGR) preconditioner implemented in the Hypre linear algebra package. - -The implementation of the hybrid FVM scheme for :ref:`CompositionalMultiphaseFlow` is in progress. diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index 14ec3cbfd45..669f721b5f4 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -458,15 +458,6 @@ real64 PhysicsSolverBase::setNextDtBasedOnNewtonIter( real64 const & currentDt ) return nextDt; } - -real64 PhysicsSolverBase::setNextDtBasedOnCFL( const geos::real64 & currentDt, geos::DomainPartition & domain ) -{ - GEOS_UNUSED_VAR( currentDt, domain ); - return LvArray::NumericLimits< real64 >::max; // i.e., not implemented -} - - - real64 PhysicsSolverBase::linearImplicitStep( real64 const & time_n, real64 const & dt, integer const GEOS_UNUSED_PARAM( cycleNumber ), diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index 92b004f8624..95a87f1a87f 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -239,17 +239,6 @@ class PhysicsSolverBase : public ExecutableGroup virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, DomainPartition & domain ); - /** - * @brief function to set the next dt based on state change - * @param [in] currentDt the current time step size - * @param[in] domain the domain object - * @return the prescribed time step size - */ - virtual real64 setNextDtBasedOnCFL( real64 const & currentDt, - DomainPartition & domain ); - - - /** * @brief Entry function for an explicit time integration step * @param time_n time at the beginning of the step diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 66c79677033..5bb8896898d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -53,7 +53,6 @@ #include "physicsSolvers/fluidFlow/kernels/compositional/SolidInternalEnergyUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/HydrostaticPressureKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/StatisticsKernel.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp" #if defined( __INTEL_COMPILER ) #pragma GCC optimize "O0" @@ -78,7 +77,6 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, m_allowCompDensChopping( 1 ), m_useTotalMassEquation( 1 ), m_useSimpleAccumulation( 1 ), - m_useNewGravity( 0 ), m_minCompDens( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision ) { //START_SPHINX_INCLUDE_00 @@ -147,12 +145,6 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, setApplyDefaultValue( 1 ). setDescription( "Flag indicating whether local (cell-wise) chopping of negative compositions is allowed" ); - this->registerWrapper( viewKeyStruct::targetFlowCFLString(), &m_targetFlowCFL ). - setApplyDefaultValue( -1. ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Target CFL condition `CFL condition `_" - " when computing the next timestep." ); - this->registerWrapper( viewKeyStruct::useTotalMassEquationString(), &m_useTotalMassEquation ). setSizedFromParent( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -165,12 +157,6 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, setApplyDefaultValue( 1 ). setDescription( "Flag indicating whether simple accumulation form is used" ); - this->registerWrapper( viewKeyStruct::useNewGravityString(), &m_useNewGravity ). - setSizedFromParent( 0 ). - setInputFlag( InputFlags::OPTIONAL ). - setApplyDefaultValue( 0 ). - setDescription( "Flag indicating whether new gravity treatment is used" ); - this->registerWrapper( viewKeyStruct::minCompDensString(), &m_minCompDens ). setSizedFromParent( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -359,16 +345,6 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) InputError ); } - if( m_targetFlowCFL > 0 ) - { - subRegion.registerField< fields::flow::phaseOutflux >( getName() ). - reference().resizeDimension< 1 >( m_numPhases ); - subRegion.registerField< fields::flow::componentOutflux >( getName() ). - reference().resizeDimension< 1 >( m_numComponents ); - subRegion.registerField< fields::flow::phaseCFLNumber >( getName() ); - subRegion.registerField< fields::flow::componentCFLNumber >( getName() ); - } - string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); @@ -1349,6 +1325,15 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom { GEOS_MARK_FUNCTION; + using namespace isothermalCompositionalMultiphaseBaseKernels; + + BitFlags< KernelFlags > kernelFlags; + if( m_useTotalMassEquation ) + kernelFlags.set( KernelFlags::TotalMassEquation ); + if( m_useSimpleAccumulation ) + kernelFlags.set( KernelFlags::SimpleAccumulation ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, arrayView1d< string const > const & regionNames ) @@ -1371,7 +1356,7 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), - m_useTotalMassEquation, + kernelFlags, dofKey, subRegion, fluid, @@ -1386,8 +1371,7 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), - m_useTotalMassEquation, - m_useSimpleAccumulation, + kernelFlags, dofKey, subRegion, fluid, @@ -2181,172 +2165,6 @@ real64 CompositionalMultiphaseBase::setNextDtBasedOnStateChange( real64 const & return std::min( std::min( nextDtPressure, std::min( nextDtPhaseVolFrac, nextDtCompDens ) ), nextDtTemperature ); } -real64 CompositionalMultiphaseBase::setNextDtBasedOnCFL( const geos::real64 & currentDt, geos::DomainPartition & domain ) -{ - - real64 maxPhaseCFL, maxCompCFL; - - computeCFLNumbers( domain, currentDt, maxPhaseCFL, maxCompCFL ); - - GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max phase CFL number = {}", getName(), maxPhaseCFL ) ); - GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max component CFL number = {} ", getName(), maxCompCFL ) ); - - return std::min( m_targetFlowCFL * currentDt / maxCompCFL, - m_targetFlowCFL * currentDt / maxPhaseCFL ); -} - -void CompositionalMultiphaseBase::computeCFLNumbers( geos::DomainPartition & domain, const geos::real64 & dt, - geos::real64 & maxPhaseCFL, geos::real64 & maxCompCFL ) -{ - GEOS_MARK_FUNCTION; - - integer const numPhases = numFluidPhases(); - integer const numComps = numFluidComponents(); - - // Step 1: reset the arrays involved in the computation of CFL numbers - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - mesh.getElemManager().forElementSubRegions( regionNames, - [&]( localIndex const, - ElementSubRegionBase & subRegion ) - { - arrayView2d< real64, compflow::USD_PHASE > const & phaseOutflux = - subRegion.getField< fields::flow::phaseOutflux >(); - arrayView2d< real64, compflow::USD_COMP > const & compOutflux = - subRegion.getField< fields::flow::componentOutflux >(); - phaseOutflux.zero(); - compOutflux.zero(); - } ); - - // Step 2: compute the total volumetric outflux for each reservoir cell by looping over faces - NumericalMethodsManager & numericalMethodManager = domain.getNumericalMethodManager(); - FiniteVolumeManager & fvManager = numericalMethodManager.getFiniteVolumeManager(); - FluxApproximationBase & fluxApprox = fvManager.getFluxApproximation( getDiscretizationName() ); - - isothermalCompositionalMultiphaseFVMKernels:: - CFLFluxKernel::CompFlowAccessors compFlowAccessors( mesh.getElemManager(), getName() ); - isothermalCompositionalMultiphaseFVMKernels:: - CFLFluxKernel::MultiFluidAccessors multiFluidAccessors( mesh.getElemManager(), getName() ); - isothermalCompositionalMultiphaseFVMKernels:: - CFLFluxKernel::PermeabilityAccessors permeabilityAccessors( mesh.getElemManager(), getName() ); - isothermalCompositionalMultiphaseFVMKernels:: - CFLFluxKernel::RelPermAccessors relPermAccessors( mesh.getElemManager(), getName() ); - - // TODO: find a way to compile with this modifiable accessors in CompFlowAccessors, and remove them from here - ElementRegionManager::ElementViewAccessor< arrayView2d< real64, compflow::USD_PHASE > > const phaseOutfluxAccessor = - mesh.getElemManager().constructViewAccessor< array2d< real64, compflow::LAYOUT_PHASE >, - arrayView2d< real64, compflow::USD_PHASE > >( fields::flow::phaseOutflux::key() ); - - ElementRegionManager::ElementViewAccessor< arrayView2d< real64, compflow::USD_COMP > > const compOutfluxAccessor = - mesh.getElemManager().constructViewAccessor< array2d< real64, compflow::LAYOUT_COMP >, - arrayView2d< real64, compflow::USD_COMP > >( fields::flow::componentOutflux::key() ); - - - fluxApprox.forAllStencils( mesh, [&] ( auto & stencil ) - { - typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); - - // While this kernel is waiting for a factory class, pass all the accessors here - isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector1 - < isothermalCompositionalMultiphaseFVMKernels::CFLFluxKernel >( numComps, - numPhases, - m_useNewGravity, - dt, - stencilWrapper, - compFlowAccessors.get( fields::flow::pressure{} ), - compFlowAccessors.get( fields::flow::gravityCoefficient{} ), - compFlowAccessors.get( fields::flow::phaseVolumeFraction{} ), - permeabilityAccessors.get( fields::permeability::permeability{} ), - permeabilityAccessors.get( fields::permeability::dPerm_dPressure{} ), - relPermAccessors.get( fields::relperm::phaseRelPerm{} ), - multiFluidAccessors.get( fields::multifluid::phaseViscosity{} ), - multiFluidAccessors.get( fields::multifluid::phaseDensity{} ), - multiFluidAccessors.get( fields::multifluid::phaseMassDensity{} ), - multiFluidAccessors.get( fields::multifluid::phaseCompFraction{} ), - phaseOutfluxAccessor.toNestedView(), - compOutfluxAccessor.toNestedView() ); - } ); - } ); - - // Step 3: finalize the (cell-based) computation of the CFL numbers - real64 localMaxPhaseCFLNumber = 0.0; - real64 localMaxCompCFLNumber = 0.0; - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - mesh.getElemManager().forElementSubRegions( regionNames, - [&]( localIndex const, - ElementSubRegionBase & subRegion ) - { - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseOutflux = - subRegion.getField< fields::flow::phaseOutflux >(); - arrayView2d< real64 const, compflow::USD_COMP > const & compOutflux = - subRegion.getField< fields::flow::componentOutflux >(); - - arrayView1d< real64 > const & phaseCFLNumber = subRegion.getField< fields::flow::phaseCFLNumber >(); - arrayView1d< real64 > const & compCFLNumber = subRegion.getField< fields::flow::componentCFLNumber >(); - - arrayView1d< real64 const > const & volume = subRegion.getElementVolume(); - - arrayView2d< real64 const, compflow::USD_COMP > const & compDens = - subRegion.getField< fields::flow::globalCompDensity >(); - arrayView2d< real64 const, compflow::USD_COMP > const compFrac = - subRegion.getField< fields::flow::globalCompFraction >(); - arrayView2d< real64, compflow::USD_PHASE > const phaseVolFrac = - subRegion.getField< fields::flow::phaseVolumeFraction >(); - - Group const & constitutiveModels = subRegion.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() ); - - string const & fluidName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::fluidNamesString() ); - MultiFluidBase const & fluid = constitutiveModels.getGroup< MultiFluidBase >( fluidName ); - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc = fluid.phaseViscosity(); - - string const & relpermName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::relPermNamesString() ); - RelativePermeabilityBase const & relperm = constitutiveModels.getGroup< RelativePermeabilityBase >( relpermName ); - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm = relperm.phaseRelPerm(); - arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac = relperm.dPhaseRelPerm_dPhaseVolFraction(); - - string const & solidName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::solidNamesString() ); - CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); - arrayView2d< real64 const > const & porosity = solid.getPorosity(); - - real64 subRegionMaxPhaseCFLNumber = 0.0; - real64 subRegionMaxCompCFLNumber = 0.0; - - isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector2 - < isothermalCompositionalMultiphaseFVMKernels::CFLKernel >( numComps, numPhases, - subRegion.size(), - volume, - porosity, - compDens, - compFrac, - phaseVolFrac, - phaseRelPerm, - dPhaseRelPerm_dPhaseVolFrac, - phaseVisc, - phaseOutflux, - compOutflux, - phaseCFLNumber, - compCFLNumber, - subRegionMaxPhaseCFLNumber, - subRegionMaxCompCFLNumber ); - - localMaxPhaseCFLNumber = LvArray::math::max( localMaxPhaseCFLNumber, subRegionMaxPhaseCFLNumber ); - localMaxCompCFLNumber = LvArray::math::max( localMaxCompCFLNumber, subRegionMaxCompCFLNumber ); - - } ); - } ); - - maxPhaseCFL = MpiWrapper::max( localMaxPhaseCFLNumber ); - maxCompCFL = MpiWrapper::max( localMaxCompCFLNumber ); - -} - - void CompositionalMultiphaseBase::resetStateToBeginningOfStep( DomainPartition & domain ) { GEOS_MARK_FUNCTION; @@ -2606,13 +2424,4 @@ bool CompositionalMultiphaseBase::checkSequentialSolutionIncrements( DomainParti return isConverged && (m_sequentialCompDensChange < m_maxSequentialCompDensChange); } -real64 CompositionalMultiphaseBase::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain ) -{ - - if( m_targetFlowCFL < 0 ) - return PhysicsSolverBase::setNextDt( currentDt, domain ); - else - return setNextDtBasedOnCFL( currentDt, domain ); -} - } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 1cec9ca685b..318662d00fc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -68,6 +68,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual void registerDataOnMesh( Group & meshBodies ) override; + virtual void registerDataForCFL( Group & meshBodies ) { GEOS_UNUSED_VAR( meshBodies ); } + /** * @defgroup Solver Interface Functions * @@ -268,7 +270,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase 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 * useNewGravityString() { return "useNewGravity"; } static constexpr char const * minCompDensString() { return "minCompDens"; } static constexpr char const * maxSequentialCompDensChangeString() { return "maxSequentialCompDensChange"; } static constexpr char const * minScalingFactorString() { return "minScalingFactor"; } @@ -370,19 +371,12 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, DomainPartition & domain ) override; - void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL ); - - /** - * @brief function to set the next time step size - * @param[in] currentDt the current time step size - * @param[in] domain the domain object - * @return the prescribed time step size - */ - real64 setNextDt( real64 const & currentDt, - DomainPartition & domain ) override; - virtual real64 setNextDtBasedOnCFL( real64 const & currentDt, - DomainPartition & domain ) override; + virtual void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL ) + { + GEOS_UNUSED_VAR( domain, dt, maxPhaseCFL, maxCompCFL ); + GEOS_ERROR( GEOS_FMT( "{}: computeCFLNumbers is not implemented for {}", getDataContext(), getCatalogName() ) ); + } virtual void initializePostInitialConditionsPreSubGroups() override; @@ -487,9 +481,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase /// flag indicating whether simple accumulation form is used integer m_useSimpleAccumulation; - /// flag indicating whether new gravity treatment is used - integer m_useNewGravity; - /// minimum allowed global component density real64 m_minCompDens; @@ -500,9 +491,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase real64 m_sequentialCompDensChange; real64 m_maxSequentialCompDensChange; - /// the targeted CFL for timestep - real64 m_targetFlowCFL; - private: /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 9a64efef23d..d2dbcb2f3a0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -53,6 +53,7 @@ #include "physicsSolvers/fluidFlow/kernels/compositional/PhaseMobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseMobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/AquiferBCKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp" namespace geos { @@ -100,6 +101,18 @@ CompositionalMultiphaseFVM::CompositionalMultiphaseFVM( const string & name, setApplyDefaultValue( ScalingType::Global ). setDescription( "Solution scaling type." "Valid options:\n* " + EnumStrings< ScalingType >::concat( "\n* " ) ); + + this->registerWrapper( viewKeyStruct::gravityDensitySchemeString(), &m_gravityDensityScheme ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( GravityDensityScheme::ArithmeticAverage ). + setDescription( "Scheme for density treatment in gravity" ); + + this->registerWrapper( viewKeyStruct::targetFlowCFLString(), &m_targetFlowCFL ). + setApplyDefaultValue( -1. ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Target CFL condition `CFL condition `_" + " when computing the next timestep." ); } void CompositionalMultiphaseFVM::postInputInitialization() @@ -112,6 +125,36 @@ void CompositionalMultiphaseFVM::postInputInitialization() } } +void CompositionalMultiphaseFVM::registerDataOnMesh( Group & meshBodies ) +{ + using namespace fields::flow; + + CompositionalMultiphaseBase::registerDataOnMesh( meshBodies ); + + if( m_targetFlowCFL > 0 ) + { + registerDataForCFL( meshBodies ); + } +} + +void CompositionalMultiphaseFVM::registerDataForCFL( Group & meshBodies ) +{ + forDiscretizationOnMeshTargets( meshBodies, [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + subRegion.registerField< fields::flow::phaseOutflux >( getName()).reference().resizeDimension< 1 >( m_numPhases ); + subRegion.registerField< fields::flow::componentOutflux >( getName()).reference().resizeDimension< 1 >( m_numComponents ); + subRegion.registerField< fields::flow::phaseCFLNumber >( getName()); + subRegion.registerField< fields::flow::componentCFLNumber >( getName()); + } ); + } ); +} + void CompositionalMultiphaseFVM::initializePreSubGroups() { CompositionalMultiphaseBase::initializePreSubGroups(); @@ -164,6 +207,20 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, { GEOS_MARK_FUNCTION; + using namespace isothermalCompositionalMultiphaseFVMKernels; + + BitFlags< KernelFlags > kernelFlags; + if( m_hasCapPressure ) + kernelFlags.set( KernelFlags::CapPressure ); + if( m_hasDiffusion ) + kernelFlags.set( KernelFlags::Diffusion ); + if( m_hasDispersion ) + kernelFlags.set( KernelFlags::Dispersion ); + if( m_useTotalMassEquation ) + kernelFlags.set( KernelFlags::TotalMassEquation ); + if( m_gravityDensityScheme == GravityDensityScheme::PhasePresence ) + kernelFlags.set( KernelFlags::CheckPhasePresenceInGravity ); + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, arrayView1d< string const > const & ) @@ -172,6 +229,13 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName ); + auto const & upwindingParams = fluxApprox.upwindingParams(); + if( upwindingParams.upwindingScheme == UpwindingScheme::C1PPU && + isothermalCompositionalMultiphaseFVMKernelUtilities::epsC1PPU > 0 ) + kernelFlags.set( KernelFlags::C1PPU ); + else if( upwindingParams.upwindingScheme == UpwindingScheme::IHU ) + kernelFlags.set( KernelFlags::IHU ); + string const & elemDofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); fluxApprox.forAllStencils( mesh, [&] ( auto & stencil ) @@ -187,8 +251,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, m_numPhases, dofManager.rankOffset(), elemDofKey, - m_hasCapPressure, - m_useTotalMassEquation, + kernelFlags, getName(), mesh.getElemManager(), stencilWrapper, @@ -206,8 +269,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, m_numPhases, dofManager.rankOffset(), elemDofKey, - m_hasCapPressure, - m_useTotalMassEquation, + kernelFlags, getName(), mesh.getElemManager(), stencilWrapper, @@ -229,10 +291,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, m_numPhases, dofManager.rankOffset(), elemDofKey, - m_hasCapPressure, - m_useTotalMassEquation, - m_useNewGravity, - fluxApprox.upwindingParams(), + kernelFlags, getName(), mesh.getElemManager(), stencilWrapper, @@ -254,9 +313,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, m_numPhases, dofManager.rankOffset(), elemDofKey, - m_hasDiffusion, - m_hasDispersion, - m_useTotalMassEquation, + kernelFlags, getName(), mesh.getElemManager(), stencilWrapper, @@ -272,9 +329,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, m_numPhases, dofManager.rankOffset(), elemDofKey, - m_hasDiffusion, - m_hasDispersion, - m_useTotalMassEquation, + kernelFlags, getName(), mesh.getElemManager(), stencilWrapper, @@ -296,6 +351,14 @@ void CompositionalMultiphaseFVM::assembleStabilizedFluxTerms( real64 const dt, { GEOS_MARK_FUNCTION; + using namespace isothermalCompositionalMultiphaseFVMKernels; + + BitFlags< KernelFlags > kernelFlags; + if( m_hasCapPressure ) + kernelFlags.set( KernelFlags::CapPressure ); + if( m_useTotalMassEquation ) + kernelFlags.set( KernelFlags::TotalMassEquation ); + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, arrayView1d< string const > const & ) @@ -321,8 +384,7 @@ void CompositionalMultiphaseFVM::assembleStabilizedFluxTerms( real64 const dt, m_numPhases, dofManager.rankOffset(), elemDofKey, - m_hasCapPressure, - m_useTotalMassEquation, + kernelFlags, getName(), mesh.getElemManager(), stencilWrapper, @@ -1008,6 +1070,12 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n, GEOS_ERROR_IF( !bcConsistent, GEOS_FMT( "CompositionalMultiphaseBase {}: inconsistent boundary conditions", getDataContext() ) ); } + using namespace isothermalCompositionalMultiphaseFVMKernels; + + BitFlags< KernelFlags > kernelFlags; + if( m_useTotalMassEquation ) + kernelFlags.set( KernelFlags::TotalMassEquation ); + FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); @@ -1070,7 +1138,7 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n, createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, dofManager.rankOffset(), - m_useTotalMassEquation, + kernelFlags, elemDofKey, getName(), faceManager, @@ -1196,6 +1264,180 @@ void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, } +real64 CompositionalMultiphaseFVM::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain ) +{ + + if( m_targetFlowCFL < 0 ) + return CompositionalMultiphaseBase::setNextDt( currentDt, domain ); + else + return setNextDtBasedOnCFL( currentDt, domain ); +} + +real64 CompositionalMultiphaseFVM::setNextDtBasedOnCFL( const geos::real64 & currentDt, geos::DomainPartition & domain ) +{ + + real64 maxPhaseCFL, maxCompCFL; + + computeCFLNumbers( domain, currentDt, maxPhaseCFL, maxCompCFL ); + + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max phase CFL number = {}", getName(), maxPhaseCFL ) ); + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max component CFL number = {} ", getName(), maxCompCFL ) ); + + return std::min( m_targetFlowCFL * currentDt / maxCompCFL, + m_targetFlowCFL * currentDt / maxPhaseCFL ); +} + +void CompositionalMultiphaseFVM::computeCFLNumbers( geos::DomainPartition & domain, const geos::real64 & dt, + geos::real64 & maxPhaseCFL, geos::real64 & maxCompCFL ) +{ + GEOS_MARK_FUNCTION; + + integer const numPhases = numFluidPhases(); + integer const numComps = numFluidComponents(); + + // Step 1: reset the arrays involved in the computation of CFL numbers + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + arrayView2d< real64, compflow::USD_PHASE > const & phaseOutflux = + subRegion.getField< fields::flow::phaseOutflux >(); + arrayView2d< real64, compflow::USD_COMP > const & compOutflux = + subRegion.getField< fields::flow::componentOutflux >(); + phaseOutflux.zero(); + compOutflux.zero(); + } ); + + // Step 2: compute the total volumetric outflux for each reservoir cell by looping over faces + NumericalMethodsManager & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase & fluxApprox = fvManager.getFluxApproximation( getDiscretizationName() ); + + isothermalCompositionalMultiphaseFVMKernels:: + CFLFluxKernel::CompFlowAccessors compFlowAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + CFLFluxKernel::MultiFluidAccessors multiFluidAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + CFLFluxKernel::PermeabilityAccessors permeabilityAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + CFLFluxKernel::RelPermAccessors relPermAccessors( mesh.getElemManager(), getName() ); + + // TODO: find a way to compile with this modifiable accessors in CompFlowAccessors, and remove them from here + ElementRegionManager::ElementViewAccessor< arrayView2d< real64, compflow::USD_PHASE > > const phaseOutfluxAccessor = + mesh.getElemManager().constructViewAccessor< array2d< real64, compflow::LAYOUT_PHASE >, + arrayView2d< real64, compflow::USD_PHASE > >( fields::flow::phaseOutflux::key() ); + + ElementRegionManager::ElementViewAccessor< arrayView2d< real64, compflow::USD_COMP > > const compOutfluxAccessor = + mesh.getElemManager().constructViewAccessor< array2d< real64, compflow::LAYOUT_COMP >, + arrayView2d< real64, compflow::USD_COMP > >( fields::flow::componentOutflux::key() ); + + + fluxApprox.forAllStencils( mesh, [&] ( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + // While this kernel is waiting for a factory class, pass all the accessors here + isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector1 + < isothermalCompositionalMultiphaseFVMKernels::CFLFluxKernel >( numComps, + numPhases, + m_gravityDensityScheme == GravityDensityScheme::PhasePresence, + dt, + stencilWrapper, + compFlowAccessors.get( fields::flow::pressure{} ), + compFlowAccessors.get( fields::flow::gravityCoefficient{} ), + compFlowAccessors.get( fields::flow::phaseVolumeFraction{} ), + permeabilityAccessors.get( fields::permeability::permeability{} ), + permeabilityAccessors.get( fields::permeability::dPerm_dPressure{} ), + relPermAccessors.get( fields::relperm::phaseRelPerm{} ), + multiFluidAccessors.get( fields::multifluid::phaseViscosity{} ), + multiFluidAccessors.get( fields::multifluid::phaseDensity{} ), + multiFluidAccessors.get( fields::multifluid::phaseMassDensity{} ), + multiFluidAccessors.get( fields::multifluid::phaseCompFraction{} ), + phaseOutfluxAccessor.toNestedView(), + compOutfluxAccessor.toNestedView() ); + } ); + } ); + + // Step 3: finalize the (cell-based) computation of the CFL numbers + real64 localMaxPhaseCFLNumber = 0.0; + real64 localMaxCompCFLNumber = 0.0; + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + arrayView2d< real64 const, compflow::USD_PHASE > const & phaseOutflux = + subRegion.getField< fields::flow::phaseOutflux >(); + arrayView2d< real64 const, compflow::USD_COMP > const & compOutflux = + subRegion.getField< fields::flow::componentOutflux >(); + + arrayView1d< real64 > const & phaseCFLNumber = subRegion.getField< fields::flow::phaseCFLNumber >(); + arrayView1d< real64 > const & compCFLNumber = subRegion.getField< fields::flow::componentCFLNumber >(); + + arrayView1d< real64 const > const & volume = subRegion.getElementVolume(); + + arrayView2d< real64 const, compflow::USD_COMP > const & compDens = + subRegion.getField< fields::flow::globalCompDensity >(); + arrayView2d< real64 const, compflow::USD_COMP > const compFrac = + subRegion.getField< fields::flow::globalCompFraction >(); + arrayView2d< real64, compflow::USD_PHASE > const phaseVolFrac = + subRegion.getField< fields::flow::phaseVolumeFraction >(); + + Group const & constitutiveModels = subRegion.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() ); + + string const & fluidName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::fluidNamesString() ); + MultiFluidBase const & fluid = constitutiveModels.getGroup< MultiFluidBase >( fluidName ); + arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseVisc = fluid.phaseViscosity(); + + string const & relpermName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::relPermNamesString() ); + RelativePermeabilityBase const & relperm = constitutiveModels.getGroup< RelativePermeabilityBase >( relpermName ); + arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelPerm = relperm.phaseRelPerm(); + arrayView4d< real64 const, relperm::USD_RELPERM_DS > const & dPhaseRelPerm_dPhaseVolFrac = relperm.dPhaseRelPerm_dPhaseVolFraction(); + + string const & solidName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::solidNamesString() ); + CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); + arrayView2d< real64 const > const & porosity = solid.getPorosity(); + + real64 subRegionMaxPhaseCFLNumber = 0.0; + real64 subRegionMaxCompCFLNumber = 0.0; + + isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector2 + < isothermalCompositionalMultiphaseFVMKernels::CFLKernel >( numComps, numPhases, + subRegion.size(), + volume, + porosity, + compDens, + compFrac, + phaseVolFrac, + phaseRelPerm, + dPhaseRelPerm_dPhaseVolFrac, + phaseVisc, + phaseOutflux, + compOutflux, + phaseCFLNumber, + compCFLNumber, + subRegionMaxPhaseCFLNumber, + subRegionMaxCompCFLNumber ); + + localMaxPhaseCFLNumber = LvArray::math::max( localMaxPhaseCFLNumber, subRegionMaxPhaseCFLNumber ); + localMaxCompCFLNumber = LvArray::math::max( localMaxCompCFLNumber, subRegionMaxCompCFLNumber ); + + } ); + } ); + + maxPhaseCFL = MpiWrapper::max( localMaxPhaseCFLNumber ); + maxCompCFL = MpiWrapper::max( localMaxCompCFLNumber ); + +} + //START_SPHINX_INCLUDE_01 REGISTER_CATALOG_ENTRY( PhysicsSolverBase, CompositionalMultiphaseFVM, string const &, Group * const ) //END_SPHINX_INCLUDE_01 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp index f14b2436740..896dc87b7f1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp @@ -25,6 +25,26 @@ namespace geos { +/** + * @brief Options for density treatment in gravity + */ +enum class GravityDensityScheme : integer +{ + ArithmeticAverage, ///< average phase density is computed using simple arithmetic average: + /// rho_ave = 0.5 * (rho_i + rho_j) + PhasePresence, ///< average phase density is computed using checking for phase presence: + /// rho_ave = 0.5 * (rho_i + rho_j) if phase is present in both cells i and j + /// = rho_i if phase is present in only cell i + /// = rho_j if phase is present in only cell j +}; + +/** + * @brief Strings for options for density treatment in gravity + */ +ENUM_STRINGS( GravityDensityScheme, + "ArithmeticAverage", + "PhasePresence" ); + /** * @class CompositionalMultiphaseFVM * @@ -85,6 +105,11 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase */ /**@{*/ + virtual void + registerDataOnMesh( Group & MeshBodies ) override; + + virtual void registerDataForCFL( Group & meshBodies ) override; + virtual void setupDofs( DomainPartition const & domain, DofManager & dofManager ) const override; @@ -150,6 +175,20 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const override; + virtual void computeCFLNumbers( DomainPartition & domain, + real64 const & dt, + real64 & maxPhaseCFL, + real64 & maxCompCFL ) override; + + /** + * @brief function to set the next time step size + * @param[in] currentDt the current time step size + * @param[in] domain the domain object + * @return the prescribed time step size + */ + real64 setNextDt( real64 const & currentDt, + DomainPartition & domain ) override; + struct viewKeyStruct : CompositionalMultiphaseBase::viewKeyStruct { // DBC parameters @@ -162,7 +201,8 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase static constexpr char const * contMultiplierDBCString() { return "contMultiplierDBC"; } // nonlinear solver parameters - static constexpr char const * scalingTypeString() { return "scalingType"; } + static constexpr char const * scalingTypeString() { return "scalingType"; } + static constexpr char const * gravityDensitySchemeString() { return "gravityDensityScheme"; } }; /** @@ -191,8 +231,10 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase virtual void postInputInitialization() override; - virtual void - initializePreSubGroups() override; + virtual void initializePreSubGroups() override; + + real64 setNextDtBasedOnCFL( real64 const & currentDt, + DomainPartition & domain ); struct DBCParameters { @@ -215,6 +257,12 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase /// Solution scaling type ScalingType m_scalingType; + /// scheme for density treatment in gravity + GravityDensityScheme m_gravityDensityScheme; + + /// the targeted CFL for timestep + real64 m_targetFlowCFL; + private: /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp index 254c2a72a17..7f44f5f5c16 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp @@ -142,22 +142,13 @@ void CompositionalMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) } } } - - // if we have to compute CFL numbers later, we need to register additional variables - if( m_computeCFLNumbers ) - { - elemManager.forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) - { - subRegion.registerField< fields::flow::phaseOutflux >( getName() ). - reference().resizeDimension< 1 >( numPhases ); - subRegion.registerField< fields::flow::componentOutflux >( getName() ). - reference().resizeDimension< 1 >( numComps ); - subRegion.registerField< fields::flow::phaseCFLNumber >( getName() ); - subRegion.registerField< fields::flow::componentCFLNumber >( getName() ); - } ); - } } ); + + // if we have to compute CFL numbers later, we need to register additional variables + if( m_computeCFLNumbers ) + { + m_solver->registerDataForCFL( meshBodies ); + } } bool CompositionalMultiphaseStatistics::execute( real64 const time_n, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp index c88ec507449..e283affcf43 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp @@ -474,8 +474,7 @@ class AccumulationKernelFactory createAndLaunch( integer const numComps, integer const numPhases, globalIndex const rankOffset, - integer const useTotalMassEquation, - integer const useSimpleAccumulation, + BitFlags< KernelFlags > kernelFlags, string const dofKey, ElementSubRegionBase const & subRegion, constitutive::MultiFluidBase const & fluid, @@ -488,12 +487,6 @@ class AccumulationKernelFactory integer constexpr NUM_COMP = NC(); integer constexpr NUM_DOF = NC()+1; - BitFlags< KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( KernelFlags::TotalMassEquation ); - if( useSimpleAccumulation ) - kernelFlags.set( KernelFlags::SimpleAccumulation ); - AccumulationKernel< NUM_COMP, NUM_DOF > kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs, kernelFlags ); AccumulationKernel< NUM_COMP, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp index 6b31d7b9165..f92335dd726 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp @@ -79,7 +79,7 @@ struct C1PPUPhaseFlux compute( integer const numPhase, integer const ip, integer const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -111,8 +111,9 @@ struct C1PPUPhaseFlux real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; real64 dGravHead_dC[numFluxSupportPoints][numComp]{}; - PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, useNewGravity, seri, sesri, sei, trans, dTrans_dPres, pres, - gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, + PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, checkPhasePresenceInGravity, + seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, + phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.cpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.cpp index 0b34792763f..ab2b2838d7b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.cpp @@ -39,7 +39,7 @@ inline void CFLFluxKernel:: compute( integer const numPhases, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const stencilSize, real64 const dt, arraySlice1d< localIndex const > const seri, @@ -69,7 +69,7 @@ CFLFluxKernel:: real64 gravHead{}; // calculate quantities on primary connected cells - calculateMeanDensity( useNewGravity, ip, stencilSize, seri, sesri, sei, phaseVolFrac, phaseMassDens, densMean ); + calculateMeanDensity( checkPhasePresenceInGravity, ip, stencilSize, seri, sesri, sei, phaseVolFrac, phaseMassDens, densMean ); //***** calculation of phase volumetric flux ***** @@ -123,7 +123,7 @@ CFLFluxKernel:: GEOS_HOST_DEVICE inline void -CFLFluxKernel::calculateMeanDensity( integer const useNewGravity, integer const ip, localIndex const stencilSize, +CFLFluxKernel::calculateMeanDensity( integer const checkPhasePresenceInGravity, integer const ip, localIndex const stencilSize, arraySlice1d< localIndex const > const seri, arraySlice1d< localIndex const > const sesri, arraySlice1d< localIndex const > const sei, @@ -139,7 +139,7 @@ CFLFluxKernel::calculateMeanDensity( integer const useNewGravity, integer const localIndex const ei = sei[i]; bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( useNewGravity && !phaseExists ) + if( checkPhasePresenceInGravity && !phaseExists ) { continue; } @@ -156,7 +156,7 @@ CFLFluxKernel::calculateMeanDensity( integer const useNewGravity, integer const template< integer NC, typename STENCILWRAPPER_TYPE > void CFLFluxKernel:: launch( integer const numPhases, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, real64 const dt, STENCILWRAPPER_TYPE const & stencilWrapper, ElementViewConst< arrayView1d< real64 const > > const & pres, @@ -189,7 +189,7 @@ void CFLFluxKernel:: dTrans_dPres ); CFLFluxKernel::compute< NC >( numPhases, - useNewGravity, + checkPhasePresenceInGravity, sei[iconn].size(), dt, seri[iconn], @@ -213,7 +213,7 @@ void CFLFluxKernel:: template \ void CFLFluxKernel:: \ launch< NC, STENCILWRAPPER_TYPE >( integer const numPhases, \ - integer const useNewGravity, \ + integer const checkPhasePresenceInGravity, \ real64 const dt, \ STENCILWRAPPER_TYPE const & stencil, \ ElementViewConst< arrayView1d< real64 const > > const & pres, \ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp index a13c9170bc6..b4574fe1ece 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp @@ -87,7 +87,7 @@ struct CFLFluxKernel template< integer NC > GEOS_HOST_DEVICE inline static void compute( integer const numPhases, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const stencilSize, real64 const dt, arraySlice1d< localIndex const > const seri, @@ -106,7 +106,8 @@ struct CFLFluxKernel ElementView< arrayView2d< real64, compflow::USD_COMP > > const & compOutflux ); GEOS_HOST_DEVICE inline static void - calculateMeanDensity( integer const useNewGravity, integer const ip, localIndex const stencilSize, + calculateMeanDensity( integer const checkPhasePresenceInGravity, + integer const ip, localIndex const stencilSize, arraySlice1d< localIndex const > const seri, arraySlice1d< localIndex const > const sesri, arraySlice1d< localIndex const > const sei, @@ -117,7 +118,7 @@ struct CFLFluxKernel template< integer NC, typename STENCILWRAPPER_TYPE > static void launch( integer const numPhases, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, real64 const dt, STENCILWRAPPER_TYPE const & stencil, ElementViewConst< arrayView1d< real64 const > > const & pres, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DiffusionDispersionFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DiffusionDispersionFluxComputeKernel.hpp index f4293d5f83d..a078cb402d1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DiffusionDispersionFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DiffusionDispersionFluxComputeKernel.hpp @@ -725,9 +725,7 @@ class DiffusionDispersionFluxComputeKernelFactory integer const numPhases, globalIndex const rankOffset, string const & dofKey, - integer const hasDiffusion, - integer const hasDispersion, - integer const useTotalMassEquation, + BitFlags< KernelFlags > kernelFlags, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, @@ -744,10 +742,6 @@ class DiffusionDispersionFluxComputeKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); - BitFlags< KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( KernelFlags::TotalMassEquation ); - using kernelType = DiffusionDispersionFluxComputeKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); @@ -760,7 +754,8 @@ class DiffusionDispersionFluxComputeKernelFactory diffusionAccessors, dispersionAccessors, porosityAccessors, dt, localMatrix, localRhs, kernelFlags ); kernelType::template launch< POLICY >( stencilWrapper.size(), - hasDiffusion, hasDispersion, + kernelFlags.isSet( KernelFlags::Diffusion ), + kernelFlags.isSet( KernelFlags::Dispersion ), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp index 3d09e878fbf..b3544545b66 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp @@ -185,7 +185,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl // // We use the lambda below (called **inside** the phase loop of the base computeFlux) to compute dissipation terms Base::computeFlux( iconn, stack, [&] ( integer const ip, - integer const GEOS_UNUSED_PARAM( useNewGravity ), + integer const GEOS_UNUSED_PARAM( checkPhasePresenceInGravity ), localIndex const (&k)[2], localIndex const (&seri)[2], localIndex const (&sesri)[2], @@ -349,8 +349,7 @@ class FluxComputeKernelFactory integer const numPhases, globalIndex const rankOffset, string const & dofKey, - integer const hasCapPressure, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, @@ -374,12 +373,6 @@ class FluxComputeKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); - BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags; - if( hasCapPressure ) - kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure ); - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::TotalMassEquation ); - using KERNEL_TYPE = FluxComputeKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename KERNEL_TYPE::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename KERNEL_TYPE::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp index 83af007a79f..78271f54416 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp @@ -277,7 +277,7 @@ class FluxComputeKernel : public FluxComputeKernelBase ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), - m_kernelFlags.isSet( KernelFlags::NewGravity ), + m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ), seri, sesri, sei, trans, dTrans_dPres, @@ -304,7 +304,7 @@ class FluxComputeKernel : public FluxComputeKernelBase ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), - m_kernelFlags.isSet( KernelFlags::NewGravity ), + m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ), seri, sesri, sei, trans, dTrans_dPres, @@ -331,7 +331,7 @@ class FluxComputeKernel : public FluxComputeKernelBase ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), - m_kernelFlags.isSet( KernelFlags::NewGravity ), + m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ), seri, sesri, sei, trans, dTrans_dPres, @@ -355,12 +355,12 @@ class FluxComputeKernel : public FluxComputeKernelBase // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives // possible use: assemble the derivatives wrt temperature, and the flux term of the energy equation for this phase - compFluxKernelOp( ip, m_kernelFlags.isSet( KernelFlags::NewGravity ), + compFluxKernelOp( ip, m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ), k, seri, sesri, sei, connectionIndex, k_up, seri[k_up], sesri[k_up], sei[k_up], potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); - } // loop over phases + } // loop over phases /// populate local flux vector and derivatives for( integer ic = 0; ic < numComp; ++ic ) @@ -527,10 +527,7 @@ class FluxComputeKernelFactory integer const numPhases, globalIndex const rankOffset, string const & dofKey, - integer const hasCapPressure, - integer const useTotalMassEquation, - integer const useNewGravity, - UpwindingParameters upwindingParams, + BitFlags< KernelFlags > kernelFlags, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, @@ -547,20 +544,6 @@ class FluxComputeKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); - BitFlags< KernelFlags > kernelFlags; - if( hasCapPressure ) - kernelFlags.set( KernelFlags::CapPressure ); - if( useTotalMassEquation ) - kernelFlags.set( KernelFlags::TotalMassEquation ); - if( useNewGravity ) - kernelFlags.set( KernelFlags::NewGravity ); - if( upwindingParams.upwindingScheme == UpwindingScheme::C1PPU && - isothermalCompositionalMultiphaseFVMKernelUtilities::epsC1PPU > 0 ) - kernelFlags.set( KernelFlags::C1PPU ); - else if( upwindingParams.upwindingScheme == UpwindingScheme::IHU ) - kernelFlags.set( KernelFlags::IHU ); - - using kernelType = FluxComputeKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp index ff96044c991..c7b85da8d97 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp @@ -45,17 +45,19 @@ enum class KernelFlags { /// Flag to specify whether capillary pressure is used or not CapPressure = 1 << 0, // 1 + /// Flag to specify whether diffusion is used or not + Diffusion = 1 << 1, // 2 + /// Flag to specify whether dispersion is used or not + Dispersion = 1 << 2, // 4 /// Flag indicating whether total mass equation is formed or not - TotalMassEquation = 1 << 1, // 2 - /// Flag indicating whether new gravity treatment is used or not - NewGravity = 1 << 2, // 4 + TotalMassEquation = 1 << 3, // 8 + /// Flag indicating whether gravity treatment is checking phase presence or not + CheckPhasePresenceInGravity = 1 << 4, // 16 /// Flag indicating whether C1-PPU is used or not - C1PPU = 1 << 3, // 8 + C1PPU = 1 << 5, // 32 /// Flag indicating whether IHU is used or not - IHU = 1 << 4 // 16 + IHU = 1 << 6 // 64 /// Add more flags like that if needed: - // Flag6 = 1 << 5, // 32 - // Flag7 = 1 << 6, // 64 // Flag8 = 1 << 7 //128 }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp index 98ee6989a18..9365a96c830 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp @@ -140,7 +140,7 @@ upwindMobilityGravity( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, integer const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex & upwindDir, real64 & mobility, real64( &dMobility_dP), @@ -176,7 +176,7 @@ upwindMobilityGravity( localIndex const numPhase, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, upwindDir ); localIndex const er_up = seri[upwindDir]; @@ -408,7 +408,7 @@ computeFractionalFlowGravity( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, integer const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex & k_up_main, real64 & fractionalFlow, real64 ( & dFractionalFlow_dP)[numFluxSupportPoints], @@ -459,7 +459,7 @@ computeFractionalFlowGravity( localIndex const numPhase, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, k_up, mob, dMob_dP, @@ -662,7 +662,7 @@ struct computePotentialGravity GEOS_HOST_DEVICE static void compute( localIndex const GEOS_UNUSED_PARAM( numPhase ), localIndex const ip, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], @@ -702,7 +702,9 @@ struct computePotentialGravity } } - calculateMeanDensity( useNewGravity, ip, seri, sesri, sei, phaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, dProp_dComp, + calculateMeanDensity( checkPhasePresenceInGravity, ip, seri, sesri, sei, + phaseVolFrac, dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, dProp_dComp, densMean, dDensMean_dPres, dDensMean_dComp ); // compute potential difference MPFA-style @@ -730,7 +732,7 @@ struct computePotentialGravity template< localIndex numComp, localIndex numFluxSupportPoints > GEOS_HOST_DEVICE - static void calculateMeanDensity( integer const useNewGravity, + static void calculateMeanDensity( integer const checkPhasePresenceInGravity, localIndex const ip, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], @@ -750,7 +752,7 @@ struct computePotentialGravity localIndex const ei = sei[i]; bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( useNewGravity && !phaseExists ) + if( checkPhasePresenceInGravity && !phaseExists ) { continue; } @@ -878,7 +880,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, localIndex const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex( &k_up), localIndex (&k_up_o), real64 & phaseFlux, @@ -898,7 +900,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, // UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, ip, - useNewGravity, + checkPhasePresenceInGravity, seri, sesri, sei, @@ -944,7 +946,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, k_up, fflow, dFflow_dP, @@ -964,7 +966,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, //Fetch pot for phase j!=i defined as \rho_j g dz/dx UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, jp, - useNewGravity, + checkPhasePresenceInGravity, seri, sesri, sei, @@ -1012,7 +1014,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, k_up_o, mobOther, dMobOther_dP, @@ -1341,7 +1343,7 @@ class UpwindScheme ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, integer const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex & upwindDir ) { @@ -1368,7 +1370,7 @@ class UpwindScheme phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, pot ); //all definition has been changed to fit pot>0 => first cell is upstream @@ -1567,9 +1569,8 @@ class HybridUpwind : public UpwindScheme ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, integer const GEOS_UNUSED_PARAM( hasCapPressure ), - integer const useNewGravity, - real64 & potential - ) + integer const checkPhasePresenceInGravity, + real64 & potential ) { //Form total velocity @@ -1587,7 +1588,7 @@ class HybridUpwind : public UpwindScheme UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, ipp, - useNewGravity, + checkPhasePresenceInGravity, seri, sesri, sei, @@ -1715,7 +1716,7 @@ struct IHUPhaseFlux compute( integer const numPhase, integer const ip, integer const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -1763,7 +1764,8 @@ struct IHUPhaseFlux for( integer jp = 0; jp < numPhase; ++jp ) { - PPUPhaseFlux::compute( numPhase, jp, hasCapPressure, useNewGravity, + PPUPhaseFlux::compute( numPhase, jp, + hasCapPressure, checkPhasePresenceInGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, @@ -1917,7 +1919,7 @@ struct IHUPhaseFlux phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, - useNewGravity, + checkPhasePresenceInGravity, k_up_g, k_up_og, gravitationalPhaseFlux, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp index 483e8c4c2ba..468adabb4b1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp @@ -75,7 +75,7 @@ struct PPUPhaseFlux compute( integer const numPhase, integer const ip, integer const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -107,10 +107,12 @@ struct PPUPhaseFlux real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; real64 dGravHead_dC[numFluxSupportPoints][numComp]{}; - PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, useNewGravity, seri, sesri, sei, trans, dTrans_dPres, pres, - gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, - phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP, - dPresGrad_dC, dGravHead_dP, dGravHead_dC ); + PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, checkPhasePresenceInGravity, + seri, sesri, sei, trans, dTrans_dPres, pres, + gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, + dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); // *** upwinding *** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp index 2b4522343e3..b783b061712 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp @@ -46,7 +46,7 @@ struct PotGrad compute ( integer const numPhase, integer const ip, integer const hasCapPressure, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -88,7 +88,11 @@ struct PotGrad real64 gravHead = 0.0; real64 dCapPressure_dC[numComp]{}; - calculateMeanDensity( useNewGravity, ip, seri, sesri, sei, phaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, densMean, dDensMean_dP, dDensMean_dC ); + calculateMeanDensity( checkPhasePresenceInGravity, ip, + seri, sesri, sei, + phaseVolFrac, dCompFrac_dCompDens, + phaseMassDens, dPhaseMassDens, + densMean, dDensMean_dP, dDensMean_dC ); /// compute the TPFA potential difference for( integer i = 0; i < numFluxSupportPoints; i++ ) @@ -157,7 +161,7 @@ struct PotGrad template< integer numComp, integer numFluxSupportPoints > GEOS_HOST_DEVICE static void - calculateMeanDensity( integer const useNewGravity, + calculateMeanDensity( integer const checkPhasePresenceInGravity, integer const ip, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], @@ -178,7 +182,7 @@ struct PotGrad localIndex const ei = sei[i]; bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( useNewGravity && !phaseExists ) + if( checkPhasePresenceInGravity && !phaseExists ) { continue; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp index eb679f2942c..74332aaf862 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp @@ -193,7 +193,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl // // We use the lambda below (called **inside** the phase loop of the base computeFlux) to compute stabilization terms Base::computeFlux( iconn, stack, [&] ( integer const ip, - integer const GEOS_UNUSED_PARAM( useNewGravity ), + integer const GEOS_UNUSED_PARAM( checkPhasePresenceInGravity ), localIndex const (&k)[2], localIndex const (&seri)[2], localIndex const (&sesri)[2], @@ -327,8 +327,7 @@ class FluxComputeKernelFactory integer const numPhases, globalIndex const rankOffset, string const & dofKey, - integer const hasCapPressure, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, @@ -346,12 +345,6 @@ class FluxComputeKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); - BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags; - if( hasCapPressure ) - kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure ); - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::TotalMassEquation ); - using KERNEL_TYPE = FluxComputeKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename KERNEL_TYPE::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename KERNEL_TYPE::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp index 58d6aa12cf4..e5839bd023d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp @@ -311,7 +311,7 @@ class AccumulationKernelFactory createAndLaunch( localIndex const numComps, localIndex const numPhases, globalIndex const rankOffset, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, string const dofKey, ElementSubRegionBase const & subRegion, constitutive::MultiFluidBase const & fluid, @@ -325,10 +325,6 @@ class AccumulationKernelFactory localIndex constexpr NUM_COMP = NC(); localIndex constexpr NUM_DOF = NC()+2; - BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); - AccumulationKernel< NUM_COMP, NUM_DOF > kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs, kernelFlags ); AccumulationKernel< NUM_COMP, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalDiffusionDispersionFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalDiffusionDispersionFluxComputeKernel.hpp index 957d98bee43..c5f2d67951d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalDiffusionDispersionFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalDiffusionDispersionFluxComputeKernel.hpp @@ -302,9 +302,7 @@ class DiffusionDispersionFluxComputeKernelFactory integer const numPhases, globalIndex const rankOffset, string const & dofKey, - integer const hasDiffusion, - integer const hasDispersion, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, @@ -322,10 +320,6 @@ class DiffusionDispersionFluxComputeKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); - BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::TotalMassEquation ); - using kernelType = DiffusionDispersionFluxComputeKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); @@ -338,7 +332,8 @@ class DiffusionDispersionFluxComputeKernelFactory diffusionAccessors, dispersionAccessors, porosityAccessors, dt, localMatrix, localRhs, kernelFlags ); kernelType::template launch< POLICY >( stencilWrapper.size(), - hasDiffusion, hasDispersion, + kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::Diffusion ), + kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::Dispersion ), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalDirichletFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalDirichletFluxComputeKernel.hpp index 5df89fc5e3e..009ab3ce63c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalDirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalDirichletFluxComputeKernel.hpp @@ -430,7 +430,7 @@ class DirichletFluxComputeKernelFactory createAndLaunch( integer const numComps, integer const numPhases, globalIndex const rankOffset, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, string const & dofKey, string const & solverName, FaceManager const & faceManager, @@ -453,11 +453,6 @@ class DirichletFluxComputeKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); - // for now, we neglect capillary pressure in the kernel - BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::TotalMassEquation ); - using KernelType = DirichletFluxComputeKernel< NUM_COMP, NUM_DOF, typename FluidType::KernelWrapper >; typename KernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename KernelType::ThermalCompFlowAccessors thermalCompFlowAccessors( elemManager, solverName ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalFluxComputeKernel.hpp index 519b1a538bb..1eb02a44f38 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalFluxComputeKernel.hpp @@ -197,7 +197,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl // such as potGrad, phaseFlux, and the indices of the upwind cell // We use the lambda below (called **inside** the phase loop of the base computeFlux) to access these variables Base::computeFlux( iconn, stack, [&] ( integer const ip, - integer const useNewGravity, + integer const checkPhasePresenceInGravity, localIndex const (&k)[2], localIndex const (&seri)[2], localIndex const (&sesri)[2], @@ -235,7 +235,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl localIndex const ei = sei[i]; bool const phaseExists = (m_phaseVolFrac[er_up][esr_up][ei_up][ip] > 0); - if( useNewGravity && !phaseExists ) + if( checkPhasePresenceInGravity && !phaseExists ) { continue; } @@ -522,8 +522,7 @@ class FluxComputeKernelFactory integer const numPhases, globalIndex const rankOffset, string const & dofKey, - integer const hasCapPressure, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, string const & solverName, ElementRegionManager const & elemManager, STENCILWRAPPER const & stencilWrapper, @@ -541,12 +540,6 @@ class FluxComputeKernelFactory elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); - BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags; - if( hasCapPressure ) - kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure ); - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::TotalMassEquation ); - using KernelType = FluxComputeKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; typename KernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); typename KernelType::ThermalCompFlowAccessors thermalCompFlowAccessors( elemManager, solverName ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 1f13044d87d..b891467863d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -260,6 +260,8 @@ void CompositionalMultiphaseWell::registerDataOnMesh( Group & meshBodies ) wellControls.registerWrapper< real64 >( viewKeyStruct::massDensityString() ); + wellControls.registerWrapper< real64 >( viewKeyStruct::currentMassRateString() ); + // write rates output header // the rank that owns the reference well element is responsible if( m_writeCSV > 0 && subRegion.isLocallyOwned() ) @@ -452,6 +454,7 @@ void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n WellControls::Control const currentControl = wellControls.getControl(); real64 const & targetTotalRate = wellControls.getTargetTotalRate( time_n + dt ); real64 const & targetPhaseRate = wellControls.getTargetPhaseRate( time_n + dt ); + real64 const & targetMassRate = wellControls.getTargetMassRate( time_n + dt ); GEOS_THROW_IF( wellControls.isInjector() && currentControl == WellControls::Control::PHASEVOLRATE, "WellControls " << wellControls.getDataContext() << @@ -470,6 +473,18 @@ void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n "WellControls " << wellControls.getDataContext() << ": Target phase rate cannot be used for injectors", InputError ); + GEOS_THROW_IF( wellControls.isProducer() && !isZero( targetTotalRate ), + "WellControls " << wellControls.getDataContext() << + ": Target total rate cannot be used for producers", + InputError ); + GEOS_THROW_IF( wellControls.isProducer() && !isZero( targetMassRate ), + "WellControls " << wellControls.getDataContext() << + ": Target mass rate cannot be used for producers", + InputError ); + GEOS_THROW_IF( !m_useMass && !isZero( targetMassRate ), + "WellControls " << wellControls.getDataContext() << + ": Target mass rate cannot with useMass=0", + InputError ); // The user always provides positive rates, but these rates are later multiplied by -1 internally for producers GEOS_THROW_IF( wellControls.isProducer() && targetPhaseRate > 0.0, @@ -686,6 +701,10 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubReg real64 & currentTotalVolRate = wellControls.getReference< real64 >( CompositionalMultiphaseWell::viewKeyStruct::currentTotalVolRateString() ); + + real64 & currentMassRate = + wellControls.getReference< real64 >( CompositionalMultiphaseWell::viewKeyStruct::currentMassRateString() ); + arrayView1d< real64 > const & dCurrentTotalVolRate = wellControls.getReference< array1d< real64 > >( CompositionalMultiphaseWell::viewKeyStruct::dCurrentTotalVolRateString() ); @@ -721,6 +740,7 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubReg dCurrentTotalVolRate, currentPhaseVolRate, dCurrentPhaseVolRate, + ¤tMassRate, &iwelemRef, &logLevel, &wellControlsName, @@ -757,7 +777,8 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubReg // Step 2: update the total volume rate real64 const currentTotalRate = connRate[iwelemRef]; - + // Assumes useMass is true + currentMassRate = currentTotalRate; // Step 2.1: compute the inverse of the total density and derivatives massDensity = totalDens[iwelemRef][0]; real64 const totalDensInv = 1.0 / totalDens[iwelemRef][0]; @@ -1090,6 +1111,9 @@ void CompositionalMultiphaseWell::assembleFluxTerms( real64 const & time, { GEOS_MARK_FUNCTION; + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; + if( m_useTotalMassEquation ) + kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); string const wellDofKey = dofManager.getKey( wellElementDofName()); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, @@ -1114,7 +1138,7 @@ void CompositionalMultiphaseWell::assembleFluxTerms( real64 const & time, createAndLaunch< parallelDevicePolicy<> >( numComponents, dt, dofManager.rankOffset(), - m_useTotalMassEquation, + kernelFlags, wellDofKey, well_controls, subRegion, @@ -1129,7 +1153,7 @@ void CompositionalMultiphaseWell::assembleFluxTerms( real64 const & time, createAndLaunch< parallelDevicePolicy<> >( numComponents, dt, dofManager.rankOffset(), - m_useTotalMassEquation, + kernelFlags, wellDofKey, well_controls, subRegion, @@ -1152,6 +1176,11 @@ void CompositionalMultiphaseWell::assembleAccumulationTerms( real64 const & time GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( time ); GEOS_UNUSED_VAR( dt ); + + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; + if( m_useTotalMassEquation ) + kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); + string const wellDofKey = dofManager.getKey( wellElementDofName() ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, @@ -1177,7 +1206,7 @@ void CompositionalMultiphaseWell::assembleAccumulationTerms( real64 const & time numPhases, wellControls.isProducer(), dofManager.rankOffset(), - m_useTotalMassEquation, + kernelFlags, wellDofKey, subRegion, fluid, @@ -1192,7 +1221,7 @@ void CompositionalMultiphaseWell::assembleAccumulationTerms( real64 const & time numPhases, wellControls.isProducer(), dofManager.rankOffset(), - m_useTotalMassEquation, + kernelFlags, wellDofKey, subRegion, fluid, @@ -1951,6 +1980,12 @@ void CompositionalMultiphaseWell::assemblePressureRelations( real64 const & time GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, GEOS_FMT( "Control switch for well {} from BHP constraint to phase volumetric rate constraint", subRegion.getName() ) ); } + else if( wellControls.getInputControl() == WellControls::Control::MASSRATE ) + { + wellControls.switchToMassRateControl( wellControls.getTargetMassRate( timeAtEndOfStep ) ); + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, + GEOS_FMT( "Control switch for well {} from BHP constraint to mass rate constraint", subRegion.getName()) ); + } else { wellControls.switchToTotalRateControl( wellControls.getTargetTotalRate( timeAtEndOfStep ) ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index aa70d4f4507..72614596376 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -296,6 +296,8 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr char const * currentTotalVolRateString() { return "currentTotalVolumetricRate"; } static constexpr char const * dCurrentTotalVolRateString() { return "dCurrentTotalVolumetricRate"; } + static constexpr char const * currentMassRateString() { return "currentMassRate"; } + static constexpr char const * dCurrentTotalVolRate_dPresString() { return "dCurrentTotalVolumetricRate_dPres"; } static constexpr char const * dCurrentTotalVolRate_dCompDensString() { return "dCurrentTotalVolumetricRate_dCompDens"; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index 7ce65b6cbec..47877476c2f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -71,27 +71,32 @@ WellControls::WellControls( string const & name, Group * const parent ) registerWrapper( viewKeyStruct::targetBHPString(), &m_targetBHP ). setDefaultValue( 0.0 ). setInputFlag( InputFlags::OPTIONAL ). + setRestartFlags( RestartFlags::WRITE_AND_READ ). setDescription( "Target bottom-hole pressure [Pa]" ); registerWrapper( viewKeyStruct::targetTotalRateString(), &m_targetTotalRate ). setDefaultValue( 0.0 ). setInputFlag( InputFlags::OPTIONAL ). + setRestartFlags( RestartFlags::WRITE_AND_READ ). setDescription( "Target total volumetric rate (if useSurfaceConditions: [surface m^3/s]; else [reservoir m^3/s])" ); registerWrapper( viewKeyStruct::targetPhaseRateString(), &m_targetPhaseRate ). setDefaultValue( 0.0 ). setInputFlag( InputFlags::OPTIONAL ). + setRestartFlags( RestartFlags::WRITE_AND_READ ). setDescription( "Target phase volumetric rate (if useSurfaceConditions: [surface m^3/s]; else [reservoir m^3/s])" ); registerWrapper( viewKeyStruct::targetMassRateString(), &m_targetMassRate ). setDefaultValue( 0.0 ). setInputFlag( InputFlags::OPTIONAL ). + setRestartFlags( RestartFlags::WRITE_AND_READ ). setDescription( "Target Mass Rate rate ( [kg^3/s])" ); registerWrapper( viewKeyStruct::targetPhaseNameString(), &m_targetPhaseName ). setRTTypeName( rtTypes::CustomTypes::groupNameRef ). setDefaultValue( "" ). setInputFlag( InputFlags::OPTIONAL ). + setRestartFlags( RestartFlags::WRITE_AND_READ ). setDescription( "Name of the target phase" ); registerWrapper( viewKeyStruct::refElevString(), &m_refElevation ). diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index 3335b7d98dc..9e202096385 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -149,6 +149,12 @@ class WellControls : public dataRepository::Group */ Control getControl() const { return m_currentControl; } + /** + * @brief Get the input control type for the well. + * @return the Control enum enforced at the well + */ + Control getInputControl() const { return m_inputControl; } + /** * @brief Getter for the reference elevation where the BHP control is enforced * @return the reference elevation diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp index 42ec6911080..01fb94d2e72 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp @@ -37,6 +37,7 @@ inline void ControlEquationHelper:: switchControl( bool const isProducer, + WellControls::Control const & inputControl, WellControls::Control const & currentControl, integer const phasePhaseIndex, real64 const & targetBHP, @@ -46,6 +47,7 @@ ControlEquationHelper:: real64 const & currentBHP, arrayView1d< real64 const > const & currentPhaseVolRate, real64 const & currentTotalVolRate, + real64 const & currentMassRate, WellControls::Control & newControl ) { // if isViable is true at the end of the following checks, no need to switch @@ -58,7 +60,7 @@ ControlEquationHelper:: // Currently, the available constraints are: // - Producer: BHP, PHASEVOLRATE - // - Injector: BHP, TOTALVOLRATE + // - Injector: BHP, TOTALVOLRATE, MASSRATE // TODO: support GRAT, WRAT, LIQUID for producers and check if any of the active constraint is violated @@ -71,6 +73,10 @@ ControlEquationHelper:: controlIsViable = ( LvArray::math::abs( currentPhaseVolRate[phasePhaseIndex] ) <= LvArray::math::abs( targetPhaseRate ) ); } // the control is viable if the reference total rate is below the max rate for injectors + else if( inputControl == WellControls::Control::MASSRATE ) + { + controlIsViable = ( LvArray::math::abs( currentMassRate ) <= LvArray::math::abs( targetMassRate ) ); + } else { controlIsViable = ( LvArray::math::abs( currentTotalVolRate ) <= LvArray::math::abs( targetTotalRate ) ); @@ -219,17 +225,6 @@ ControlEquationHelper:: if constexpr ( IS_THERMAL ) dControlEqn[COFFSET_WJ::dT] = massDensity*dCurrentTotalVolRate[COFFSET_WJ::dT]; } - // Total mass rate control - else if( currentControl == WellControls::Control::MASSRATE ) - { - controlEqn = massDensity*currentTotalVolRate - targetMassRate; - dControlEqn[COFFSET_WJ::dP] = massDensity*dCurrentTotalVolRate[COFFSET_WJ::dP]; - dControlEqn[COFFSET_WJ::dQ] = massDensity*dCurrentTotalVolRate[COFFSET_WJ::dQ]; - for( integer ic = 0; ic < NC; ++ic ) - { - dControlEqn[COFFSET_WJ::dC+ic] = massDensity*dCurrentTotalVolRate[COFFSET_WJ::dC+ic]; - } - } else { GEOS_ERROR( "This constraint is not supported in CompositionalMultiphaseWell" ); @@ -323,6 +318,7 @@ PressureRelationKernel:: // static well control data bool const isProducer = wellControls.isProducer(); WellControls::Control const currentControl = wellControls.getControl(); + WellControls::Control const inputControl = wellControls.getInputControl(); real64 const targetBHP = wellControls.getTargetBHP( timeAtEndOfStep ); real64 const targetTotalRate = wellControls.getTargetTotalRate( timeAtEndOfStep ); real64 const targetPhaseRate = wellControls.getTargetPhaseRate( timeAtEndOfStep ); @@ -344,6 +340,9 @@ PressureRelationKernel:: arrayView1d< real64 const > const & dCurrentTotalVolRate = wellControls.getReference< array1d< real64 > >( CompositionalMultiphaseWell::viewKeyStruct::dCurrentTotalVolRateString() ); + real64 const & currentMassRate = + wellControls.getReference< real64 >( CompositionalMultiphaseWell::viewKeyStruct::currentMassRateString() ); + real64 const & massDensity = wellControls.getReference< real64 >( CompositionalMultiphaseWell::viewKeyStruct::massDensityString() ); @@ -358,6 +357,7 @@ PressureRelationKernel:: { WellControls::Control newControl = currentControl; ControlEquationHelper::switchControl( isProducer, + inputControl, currentControl, targetPhaseIndex, targetBHP, @@ -367,6 +367,7 @@ PressureRelationKernel:: currentBHP, currentPhaseVolRate, currentTotalVolRate, + currentMassRate, newControl ); if( currentControl != newControl ) { @@ -737,7 +738,6 @@ RateInitializationKernel:: else if( control == WellControls::Control::MASSRATE ) { connRate[iwelem] = targetMassRate; - connRate[iwelem] = targetMassRate* totalDens[iwelem][0]; } else { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp index 41d63148bfb..c2e9078ba47 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp @@ -134,6 +134,7 @@ struct ControlEquationHelper static void switchControl( bool const isProducer, + WellControls::Control const & inputControl, WellControls::Control const & currentControl, integer const phasePhaseIndex, real64 const & targetBHP, @@ -143,6 +144,7 @@ struct ControlEquationHelper real64 const & currentBHP, arrayView1d< real64 const > const & currentPhaseVolRate, real64 const & currentTotalVolRate, + real64 const & currentMassRate, WellControls::Control & newControl ); template< integer NC, integer IS_THERMAL > @@ -1331,7 +1333,7 @@ class ElementBasedAssemblyKernelFactory localIndex const numPhases, integer const isProducer, globalIndex const rankOffset, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, string const dofKey, WellElementSubRegion const & subRegion, constitutive::MultiFluidBase const & fluid, @@ -1344,10 +1346,6 @@ class ElementBasedAssemblyKernelFactory integer constexpr istherm = IS_THERMAL(); - BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); - ElementBasedAssemblyKernel< NUM_COMP, istherm > kernel( numPhases, isProducer, rankOffset, dofKey, subRegion, fluid, localMatrix, localRhs, kernelFlags ); ElementBasedAssemblyKernel< NUM_COMP, istherm >::template @@ -1888,7 +1886,7 @@ class FaceBasedAssemblyKernelFactory createAndLaunch( integer const numComps, real64 const dt, globalIndex const rankOffset, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, string const dofKey, WellControls const & wellControls, WellElementSubRegion const & subRegion, @@ -1899,14 +1897,7 @@ class FaceBasedAssemblyKernelFactory { integer constexpr NUM_COMP = NC(); - - BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); - using kernelType = FaceBasedAssemblyKernel< NUM_COMP, 0 >; - - kernelType kernel( dt, rankOffset, dofKey, wellControls, subRegion, localMatrix, localRhs, kernelFlags ); kernelType::template launch< POLICY >( subRegion.size(), kernel ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp index e7beb347528..6e225367a36 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp @@ -638,7 +638,7 @@ class ElementBasedAssemblyKernelFactory localIndex const numPhases, integer const isProducer, globalIndex const rankOffset, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, string const dofKey, WellElementSubRegion const & subRegion, MultiFluidBase const & fluid, @@ -650,11 +650,6 @@ class ElementBasedAssemblyKernelFactory { localIndex constexpr NUM_COMP = NC(); - - BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); - ElementBasedAssemblyKernel< NUM_COMP > kernel( numPhases, isProducer, rankOffset, dofKey, subRegion, fluid, localMatrix, localRhs, kernelFlags ); ElementBasedAssemblyKernel< NUM_COMP >::template @@ -1089,7 +1084,7 @@ class FaceBasedAssemblyKernelFactory createAndLaunch( integer const numComps, real64 const dt, globalIndex const rankOffset, - integer const useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, string const dofKey, WellControls const & wellControls, WellElementSubRegion const & subRegion, @@ -1101,15 +1096,7 @@ class FaceBasedAssemblyKernelFactory { integer constexpr NUM_COMP = NC(); - - BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); - - using kernelType = FaceBasedAssemblyKernel< NUM_COMP >; - - kernelType kernel( dt, rankOffset, dofKey, wellControls, subRegion, fluid, localMatrix, localRhs, kernelFlags ); kernelType::template launch< POLICY >( subRegion.size(), kernel ); } ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp index ea2c6f34b3e..afa810d945a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp @@ -281,6 +281,10 @@ assembleCouplingTerms( real64 const time_n, this->getCatalogName(), this->getName() ), std::runtime_error ); + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; + if( Base::wellSolver()->useTotalMassEquation() ) + kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); + this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel const & mesh, arrayView1d< string const > const & regionNames ) @@ -322,7 +326,6 @@ assembleCouplingTerms( real64 const time_n, string const wellDofKey = dofManager.getKey( Base::wellSolver()->wellElementDofName() ); areWellsShut = 0; - integer useTotalMassEquation=Base::wellSolver()->useTotalMassEquation(); integer numCrossflowPerforations=0; if( isThermal ( ) ) { @@ -337,7 +340,7 @@ assembleCouplingTerms( real64 const time_n, resDofNumber, perforationData, fluid, - useTotalMassEquation, + kernelFlags, detectCrossflow, numCrossflowPerforations, localRhs, @@ -355,7 +358,7 @@ assembleCouplingTerms( real64 const time_n, resDofNumber, perforationData, fluid, - useTotalMassEquation, + kernelFlags, detectCrossflow, numCrossflowPerforations, localRhs, diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp index 5fdb80fdf62..d1d7a131817 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp @@ -306,7 +306,7 @@ class IsothermalCompositionalMultiPhaseFluxKernelFactory ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const * const perforationData, MultiFluidBase const & fluid, - integer const & useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, bool const & detectCrossflow, integer & numCrossFlowPerforations, arrayView1d< real64 > const & localRhs, @@ -317,16 +317,9 @@ class IsothermalCompositionalMultiPhaseFluxKernelFactory { integer constexpr NUM_COMP = NC(); - - BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); - - using kernelType = IsothermalCompositionalMultiPhaseFluxKernel< NUM_COMP, 0 >; - - - kernelType kernel( dt, rankOffset, wellDofKey, subRegion, resDofNumber, perforationData, fluid, localRhs, localMatrix, detectCrossflow, numCrossFlowPerforations, kernelFlags ); + kernelType kernel( dt, rankOffset, wellDofKey, subRegion, resDofNumber, perforationData, + fluid, localRhs, localMatrix, detectCrossflow, numCrossFlowPerforations, kernelFlags ); kernelType::template launch< POLICY >( perforationData->size(), kernel ); } ); @@ -559,7 +552,7 @@ class ThermalCompositionalMultiPhaseFluxKernelFactory ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const * const perforationData, MultiFluidBase const & fluid, - integer const & useTotalMassEquation, + BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, bool const & detectCrossflow, integer & numCrossFlowPerforations, arrayView1d< real64 > const & localRhs, @@ -570,16 +563,9 @@ class ThermalCompositionalMultiPhaseFluxKernelFactory { integer constexpr NUM_COMP = NC(); - - BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags; - if( useTotalMassEquation ) - kernelFlags.set( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ); - - using kernelType = ThermalCompositionalMultiPhaseFluxKernel< NUM_COMP, 1 >; - - - kernelType kernel( dt, isProducer, rankOffset, wellDofKey, subRegion, resDofNumber, perforationData, fluid, localRhs, localMatrix, detectCrossflow, numCrossFlowPerforations, kernelFlags ); + kernelType kernel( dt, isProducer, rankOffset, wellDofKey, subRegion, resDofNumber, perforationData, + fluid, localRhs, localMatrix, detectCrossflow, numCrossFlowPerforations, kernelFlags ); kernelType::template launch< POLICY >( perforationData->size(), kernel ); } ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp index 68117918b7f..8ffee2ab6f3 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp @@ -87,6 +87,14 @@ DECLARE_FIELD( strain, WRITE_AND_READ, "Average strain in cell" ); +DECLARE_FIELD( plasticStrain, + "plasticStrain", + array2dLayoutStrain, + 0, + LEVEL_0, + WRITE_AND_READ, + "Average plastic strain in cell" ); + DECLARE_FIELD( incrementalBubbleDisplacement, "incrementalBubbleDisplacement", array2d< real64 >, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 70cb8b3c56e..fa5c1685db3 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -162,6 +162,8 @@ void SolidMechanicsLagrangianFEM::registerDataOnMesh( Group & meshBodies ) setConstitutiveNamesCallSuper( subRegion ); subRegion.registerField< solidMechanics::strain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 ); + subRegion.registerField< solidMechanics::plasticStrain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 ); + } ); NodeManager & nodes = meshLevel.getNodeManager(); @@ -455,6 +457,7 @@ void SolidMechanicsLagrangianFEM::initializePostInitialConditionsPreSubGroups() m_targetNodes.insert( m_nonSendOrReceiveNodes.begin(), m_nonSendOrReceiveNodes.end() ); + } ); } ); @@ -946,23 +949,36 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS { string const & solidMaterialName = subRegion.template getReference< string >( viewKeyStruct::solidMaterialNamesString() ); SolidBase & constitutiveRelation = getConstitutiveModel< SolidBase >( subRegion, solidMaterialName ); - constitutiveRelation.saveConvergedState(); + solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >(); + solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >(); - finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); - finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) + constitutive::ConstitutivePassThru< SolidBase >::execute( constitutiveRelation, [&] ( auto & solidModel ) { - using FE_TYPE = decltype( finiteElement ); - AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< CellElementSubRegion, FE_TYPE, parallelDevicePolicy<> >( nodeManager, - mesh.getEdgeManager(), - mesh.getFaceManager(), - subRegion, - finiteElement, - disp, - strain ); + + using SOLID_TYPE = TYPEOFREF( solidModel ); + + finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); + finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) + { + using FE_TYPE = decltype( finiteElement ); + AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager, + mesh.getEdgeManager(), + mesh.getFaceManager(), + subRegion, + finiteElement, + solidModel, + disp, + uhat, + strain, + plasticStrain ); + } ); + + } ); + constitutiveRelation.saveConvergedState(); } ); } ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index 764d86d8f38..083047a880e 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp @@ -98,6 +98,16 @@ bool SolidMechanicsStateReset::execute( real64 const time_n, } nodeManager.getField< solidMechanics::totalDisplacement >().zero(); nodeManager.getField< solidMechanics::incrementalDisplacement >().zero(); + + ElementRegionManager & elemManager = mesh.getElemManager(); + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + subRegion.getField< solidMechanics::strain >().zero(); + subRegion.getField< solidMechanics::plasticStrain >().zero(); + } ); + } // Option 2: enable / disable inelastic behavior diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index d1c5567eff4..faf51f89351 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -23,6 +23,7 @@ #include "common/DataTypes.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "finiteElement/FiniteElementDispatch.hpp" +#include "constitutive/ConstitutivePassThru.hpp" #include "mesh/CellElementSubRegion.hpp" #include "mesh/utilities/AverageOverQuadraturePointsKernel.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" @@ -33,17 +34,18 @@ namespace geos * @class AverageStrainOverQuadraturePoints * @tparam SUBREGION_TYPE the subRegion type * @tparam FE_TYPE the finite element type + * @tparam SOLID_TYPE the solid mechanics constitutuve type */ -template< typename SUBREGION_TYPE, - typename FE_TYPE > +template< typename FE_TYPE, + typename SOLID_TYPE > class AverageStrainOverQuadraturePoints : - public AverageOverQuadraturePointsBase< SUBREGION_TYPE, + public AverageOverQuadraturePointsBase< CellElementSubRegion, FE_TYPE > { public: /// Alias for the base class; - using Base = AverageOverQuadraturePointsBase< SUBREGION_TYPE, + using Base = AverageOverQuadraturePointsBase< CellElementSubRegion, FE_TYPE >; using Base::m_elementVolume; @@ -63,24 +65,31 @@ class AverageStrainOverQuadraturePoints : AverageStrainOverQuadraturePoints( NodeManager & nodeManager, EdgeManager const & edgeManager, FaceManager const & faceManager, - SUBREGION_TYPE const & elementSubRegion, + CellElementSubRegion const & elementSubRegion, FE_TYPE const & finiteElementSpace, + SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, - fields::solidMechanics::arrayView2dLayoutStrain const avgStrain ): + fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc, + fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, + fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain ): Base( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace ), + m_solidUpdate( solidModel.createKernelUpdates()), m_displacement( displacement ), - m_avgStrain( avgStrain ) + m_displacementInc( displacementInc ), + m_avgStrain( avgStrain ), + m_avgPlasticStrain( avgPlasticStrain ) {} /** * @copydoc finiteElement::KernelBase::StackVariables */ struct StackVariables : Base::StackVariables - {real64 uLocal[FE_TYPE::maxSupportPoints][3]; }; + {real64 uLocal[FE_TYPE::maxSupportPoints][3]; + real64 uHatLocal[FE_TYPE::maxSupportPoints][3]; }; /** * @brief Performs the setup phase for the kernel. @@ -99,12 +108,14 @@ class AverageStrainOverQuadraturePoints : for( int i = 0; i < 3; ++i ) { stack.uLocal[a][i] = m_displacement[localNodeIndex][i]; + stack.uHatLocal[a][i] = m_displacementInc[localNodeIndex][i]; } } for( int icomp = 0; icomp < 6; ++icomp ) { m_avgStrain[k][icomp] = 0.0; + //m_avgPlasticStrain[k][icomp] = 0.0; } } @@ -124,11 +135,26 @@ class AverageStrainOverQuadraturePoints : real64 dNdX[ FE_TYPE::maxSupportPoints ][3]; real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, stack.feStack, dNdX ); real64 strain[6] = {0.0}; + real64 strainInc[6] = {0.0}; FE_TYPE::symmetricGradient( dNdX, stack.uLocal, strain ); + FE_TYPE::symmetricGradient( dNdX, stack.uHatLocal, strainInc ); + + real64 elasticStrainInc[6] = {0.0}; + m_solidUpdate.getElasticStrainInc( k, q, elasticStrainInc ); + + real64 conversionFactor[6] = {1.0, 1.0, 1.0, 0.5, 0.5, 0.5}; // used for converting from engineering shear to tensor shear for( int icomp = 0; icomp < 6; ++icomp ) { - m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k]; + m_avgStrain[k][icomp] += conversionFactor[icomp]*detJxW*strain[icomp]/m_elementVolume[k]; + + // This is a hack to handle boundary conditions such as those seen in plane-strain wellbore problems + // Essentially, if bcs are constraining the strain (and thus total displacement), we do not accumulate any plastic strain (regardless + // of stresses in material law) + if( std::abs( strainInc[icomp] ) > 1.0e-8 ) + { + m_avgPlasticStrain[k][icomp] += conversionFactor[icomp]*detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; + } } } @@ -160,12 +186,21 @@ class AverageStrainOverQuadraturePoints : protected: + /// The material + typename SOLID_TYPE::KernelWrapper const m_solidUpdate; + /// The displacement solution fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const m_displacement; + /// The displacement increment + fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const m_displacementInc; + /// The average strain fields::solidMechanics::arrayView2dLayoutStrain const m_avgStrain; + /// The average plastic strain + fields::solidMechanics::arrayView2dLayoutStrain const m_avgPlasticStrain; + }; @@ -182,6 +217,7 @@ class AverageStrainOverQuadraturePointsKernelFactory * @brief Create a new kernel and launch * @tparam SUBREGION_TYPE the subRegion type * @tparam FE_TYPE the finite element type + * @tparam SOLID_TYPE the constitutive type * @tparam POLICY the kernel policy * @param nodeManager the node manager * @param edgeManager the edge manager @@ -191,23 +227,26 @@ class AverageStrainOverQuadraturePointsKernelFactory * @param property the property at quadrature points * @param averageProperty the property averaged over quadrature points */ - template< typename SUBREGION_TYPE, - typename FE_TYPE, + template< typename FE_TYPE, + typename SOLID_TYPE, typename POLICY > static void createAndLaunch( NodeManager & nodeManager, EdgeManager const & edgeManager, FaceManager const & faceManager, - SUBREGION_TYPE const & elementSubRegion, + CellElementSubRegion const & elementSubRegion, FE_TYPE const & finiteElementSpace, + SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, - fields::solidMechanics::arrayView2dLayoutStrain const avgStrain ) + fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc, + fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, + fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain ) { - AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE > + AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE > kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace, - displacement, avgStrain ); + solidModel, displacement, displacementInc, avgStrain, avgPlasticStrain ); - AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >::template + AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >::template kernelLaunch< POLICY >( elementSubRegion.size(), kernel ); } };