diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 6dfbf981b0..5f8608272b 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-pr3384-9442-381d6ba allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 0bc5c80f56..fa00433a6e 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ 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 #3384 (2024-12-16) +===================== +Added plastic strain output. + PR #3479 (2024-12-15) ===================== Refine inputFiles/compositionalMultiphaseFlow: shift reference pressures to initial pressures, make nonlinear tuning more reasonable, minimize output. diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index bed64cb6e6..97573fce5f 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 4181e8e201..01f36bf937 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 4445b88cc1..49179c5a95 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 f2c86d4439..17fabe5479 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 24b05eea71..b56bd38b53 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/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp index 68117918b7..8ffee2ab6f 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 70cb8b3c56..fa5c1685db 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 764d86d8f3..083047a880 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 d1c5567eff..faf51f8935 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 ); } };