diff --git a/inputFiles/phaseField/PhaseFieldFracture.ats b/inputFiles/phaseField/PhaseFieldFracture.ats index 88b731578d3..44087dbea67 100644 --- a/inputFiles/phaseField/PhaseFieldFracture.ats +++ b/inputFiles/phaseField/PhaseFieldFracture.ats @@ -28,7 +28,7 @@ decks = [ check_step=10, restartcheck_params=RestartcheckParameters(**restartcheck_params)), TestDeck( - name="PhaseFieldFracture_SpectralSplit", + name="PhaseFieldFracture_SpectralSplit_Cohesive", description= 'Testing the spectral split of the stress tensor in PF Fracture', partitions=((1, 1, 1), (2, 2, 1)), @@ -36,13 +36,37 @@ decks = [ check_step=30, restartcheck_params=RestartcheckParameters(**restartcheck_params)), TestDeck( - name="PhaseFieldFracture_VolDevSplit", + name="PhaseFieldFracture_SpectralSplit_Brittle", + description= + 'Testing the spectral split with the brittle degradation function in PF Fracture', + partitions=((1, 1, 1), (2, 2, 1)), + restart_step=15, + check_step=30, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="PhaseFieldFracture_VolDevSplit_Brittle", description= 'Testing the volumetric-deviatoric split of the stress tensor in PF Fracture', partitions=((1, 1, 1), (2, 2, 1)), restart_step=15, check_step=30, restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="PhaseFieldFracture_VolDevSplit_Cohesive", + description= + 'Testing the volumetric-deviatoric split with the cohesive degradation function in PF Fracture', + partitions=((1, 1, 1), (2, 2, 1)), + restart_step=15, + check_step=30, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="PhaseFieldFracture_MultiRegion", + description= + 'Testing PF Fracture on one region of a multi-region problem, with mechanics only elsewhere', + partitions=((1, 1, 1), (2, 2, 1)), + restart_step=10, + check_step=20, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), TestDeck(name="PhaseFieldFracture_Nucleation_smoke", description='Testing the nucleation model for PF Fracture', partitions=((1, 1, 1), (2, 2, 1)), diff --git a/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml b/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml index 4a6716a4385..5d36f307c94 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml @@ -32,7 +32,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - localDissipation="Linear" logLevel="1" targetRegions="{ Region1 }"> --> - - + criticalStrainEnergy="0" + localDissipationOption="Quadratic"/> - - - - @@ -33,8 +33,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - - localDissipation="Quadratic" logLevel="1" targetRegions="{ Region1 }"> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + \ No newline at end of file diff --git a/inputFiles/phaseField/PhaseFieldFracture_MultiRegion.xml b/inputFiles/phaseField/PhaseFieldFracture_MultiRegion.xml new file mode 100644 index 00000000000..f72ab40b5a3 --- /dev/null +++ b/inputFiles/phaseField/PhaseFieldFracture_MultiRegion.xml @@ -0,0 +1,191 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/PhaseFieldFracture_Nucleation_base.xml b/inputFiles/phaseField/PhaseFieldFracture_Nucleation_base.xml index 791cf4b0f52..529b380c615 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_Nucleation_base.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_Nucleation_base.xml @@ -35,8 +35,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - - localDissipation="Linear" irreversibilityFlag="1" damageUpperBound="1.0" logLevel="1" diff --git a/inputFiles/phaseField/PhaseFieldFracture_Nucleation_benchmark.xml b/inputFiles/phaseField/PhaseFieldFracture_Nucleation_benchmark.xml index 815c858c769..691317aa71f 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_Nucleation_benchmark.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_Nucleation_benchmark.xml @@ -58,45 +58,13 @@ lengthScale="0.2" defaultCriticalFractureEnergy="9.1e-2" criticalStrainEnergy="0" - extDrivingForceFlag="1" + fractureModelType="Nucleation" defaultTensileStrength="27.0" defaultCompressiveStrength="77.0" defaultDeltaCoefficient="3.22"/> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit_Brittle.xml b/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit_Brittle.xml new file mode 100644 index 00000000000..550f11d18c6 --- /dev/null +++ b/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit_Brittle.xml @@ -0,0 +1,61 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit_Cohesive.xml b/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit_Cohesive.xml new file mode 100644 index 00000000000..b2a6c77a889 --- /dev/null +++ b/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit_Cohesive.xml @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/PhaseFieldFracture_VolDevSplit_Brittle.xml b/inputFiles/phaseField/PhaseFieldFracture_VolDevSplit_Brittle.xml new file mode 100644 index 00000000000..0b14aa0f029 --- /dev/null +++ b/inputFiles/phaseField/PhaseFieldFracture_VolDevSplit_Brittle.xml @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/PhaseFieldFracture_VolDevSplit_Cohesive.xml b/inputFiles/phaseField/PhaseFieldFracture_VolDevSplit_Cohesive.xml new file mode 100644 index 00000000000..441abc768c7 --- /dev/null +++ b/inputFiles/phaseField/PhaseFieldFracture_VolDevSplit_Cohesive.xml @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Injection.xml b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Injection.xml index 97d75c706b4..3f9297ff800 100644 --- a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Injection.xml +++ b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Injection.xml @@ -41,7 +41,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - localDissipation="Linear" irreversibilityFlag="1" damageUpperBound="1.0" fracturePressureTermFlag="1" @@ -109,38 +108,6 @@ - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotchTension_main.xml b/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotchTension_main.xml new file mode 100644 index 00000000000..b59dcee244a --- /dev/null +++ b/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotchTension_main.xml @@ -0,0 +1,110 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotch_base.xml b/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotch_base.xml new file mode 100644 index 00000000000..b9ba13e4e14 --- /dev/null +++ b/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotch_base.xml @@ -0,0 +1,99 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/benchmark/singleEdgeNotch/meshes/singleEdgeCrackSquareShear.vtu b/inputFiles/phaseField/benchmark/singleEdgeNotch/meshes/singleEdgeCrackSquareShear.vtu new file mode 100644 index 00000000000..8869d567651 --- /dev/null +++ b/inputFiles/phaseField/benchmark/singleEdgeNotch/meshes/singleEdgeCrackSquareShear.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:954ac69fe9d18539c63eaef49c76270daefc561d0bc30fc626abf3a0ff75c5e7 +size 8556918 diff --git a/inputFiles/phaseField/benchmark/singleEdgeNotch/meshes/singleEdgeCrackSquareTension.vtu b/inputFiles/phaseField/benchmark/singleEdgeNotch/meshes/singleEdgeCrackSquareTension.vtu new file mode 100644 index 00000000000..d81b83593e6 --- /dev/null +++ b/inputFiles/phaseField/benchmark/singleEdgeNotch/meshes/singleEdgeCrackSquareTension.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fc5cd9f6016ec7c97635f540367a031b561cdc2e463314272325042a56e67e6b +size 4130730 diff --git a/inputFiles/phaseField/benchmark/threePointsBending/PhaseFieldFracture_ThreePointsBending_base.xml b/inputFiles/phaseField/benchmark/threePointsBending/PhaseFieldFracture_ThreePointsBending_base.xml new file mode 100644 index 00000000000..af89bfa5837 --- /dev/null +++ b/inputFiles/phaseField/benchmark/threePointsBending/PhaseFieldFracture_ThreePointsBending_base.xml @@ -0,0 +1,95 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/benchmark/threePointsBending/PhaseFieldFracture_ThreePointsBending_main.xml b/inputFiles/phaseField/benchmark/threePointsBending/PhaseFieldFracture_ThreePointsBending_main.xml new file mode 100644 index 00000000000..a1022a3e852 --- /dev/null +++ b/inputFiles/phaseField/benchmark/threePointsBending/PhaseFieldFracture_ThreePointsBending_main.xml @@ -0,0 +1,112 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/benchmark/threePointsBending/mesh/threePointsBendingSingleNotch.vtu b/inputFiles/phaseField/benchmark/threePointsBending/mesh/threePointsBendingSingleNotch.vtu new file mode 100644 index 00000000000..b13d4a24854 --- /dev/null +++ b/inputFiles/phaseField/benchmark/threePointsBending/mesh/threePointsBendingSingleNotch.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2c27bd8a7f93400a56a2e1c746b1e120b2b37bfa9ae5b539db044e3eff4cf436 +size 10562784 diff --git a/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/PhaseFieldFracture_ThreePointsBendWithHoles_base.xml b/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/PhaseFieldFracture_ThreePointsBendWithHoles_base.xml new file mode 100644 index 00000000000..c4c1354cece --- /dev/null +++ b/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/PhaseFieldFracture_ThreePointsBendWithHoles_base.xml @@ -0,0 +1,106 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/PhaseFieldFracture_ThreePointsBendWithHoles_main.xml b/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/PhaseFieldFracture_ThreePointsBendWithHoles_main.xml new file mode 100644 index 00000000000..450bdae5833 --- /dev/null +++ b/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/PhaseFieldFracture_ThreePointsBendWithHoles_main.xml @@ -0,0 +1,139 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/mesh/threePointsBendingWithHoles.vtu b/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/mesh/threePointsBendingWithHoles.vtu new file mode 100644 index 00000000000..9cb1821bf0b --- /dev/null +++ b/inputFiles/phaseField/benchmark/threePointsBendingWithHoles/mesh/threePointsBendingWithHoles.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e1c1f36ee6c041755fe55b5b7f4c08e611f3b7ef2891d68f3cc1095db10c58e6 +size 49152688 diff --git a/src/coreComponents/constitutive/solid/Damage.cpp b/src/coreComponents/constitutive/solid/Damage.cpp index 525d6e0f0bf..acd9b3a7a9f 100644 --- a/src/coreComponents/constitutive/solid/Damage.cpp +++ b/src/coreComponents/constitutive/solid/Damage.cpp @@ -48,10 +48,16 @@ Damage< BASE >::Damage( string const & name, Group * const parent ): setInputFlag( InputFlags::OPTIONAL ). setDescription( "The lower limit of the degradation function" ); - this->registerWrapper( viewKeyStruct::extDrivingForceFlagString(), &m_extDrivingForceFlag ). - setApplyDefaultValue( 0 ). + this->registerWrapper( viewKeyStruct::fractureModelTypeString(), &m_fractureModelType ). + setApplyDefaultValue( FractureModelType::Brittle ). setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Whether to have external driving force. Can be 0 or 1" ); + setDescription( "Type of crack/fracture model. Can be Brittle, Cohesive, or Nucleation" ); + + this->registerWrapper( viewKeyStruct::localDissipationOptionString(), &m_localDissipationOption ). + setApplyDefaultValue( LocalDissipationOption::Linear ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Type of local dissipation function. Must match the damage solver's " + "localDissipation option. Can be Linear or Quadratic" ); this->registerWrapper( viewKeyStruct::defaultTensileStrengthString(), &m_defaultTensileStrength ). setApplyDefaultValue( 0.0 ). @@ -76,7 +82,9 @@ Damage< BASE >::Damage( string const & name, Group * const parent ): this->template registerField< fields::solid::damageGrad >( &m_damageGrad ); - this->template registerField< fields::solid::strainEnergyDensity >( &m_strainEnergyDensity ); + this->template registerField< fields::solid::crackDrivingForce >( &m_crackDrivingForce ); + + this->template registerField< fields::solid::oldCrackDrivingForce >( &m_oldCrackDrivingForce ); this->template registerField< fields::solid::volStrain >( &m_volStrain ); @@ -99,35 +107,25 @@ void Damage< BASE >::postInputInitialization() { BASE::postInputInitialization(); - GEOS_ERROR_IF( m_extDrivingForceFlag != 0 && m_extDrivingForceFlag!= 1, - "invalid external driving force flag option - must" - " be 0 or 1", + GEOS_ERROR_IF( m_fractureModelType != FractureModelType::Brittle && m_localDissipationOption != LocalDissipationOption::Linear, + "the Cohesive and Nucleation crack models are only supported with the Linear " + "local dissipation option", + this->getDataContext() ); + GEOS_ERROR_IF( m_fractureModelType == FractureModelType::Cohesive && m_criticalStrainEnergy <= 0.0, + "criticalStrainEnergy must be positive when the Cohesive crack model is used", this->getDataContext() ); - GEOS_ERROR_IF( m_extDrivingForceFlag == 1 && m_defaultTensileStrength <= 0.0, + GEOS_ERROR_IF( m_fractureModelType == FractureModelType::Nucleation && m_defaultTensileStrength <= 0.0, "tensile strength must be input and positive when the" - " external driving force flag is turned on", + " Nucleation crack model is used", this->getDataContext() ); - GEOS_ERROR_IF( m_extDrivingForceFlag == 1 && m_defaultCompressiveStrength <= 0.0, + GEOS_ERROR_IF( m_fractureModelType == FractureModelType::Nucleation && m_defaultCompressiveStrength <= 0.0, "compressive strength must be input and positive when the" - " external driving force flag is turned on", + " Nucleation crack model is used", this->getDataContext() ); - GEOS_ERROR_IF( m_extDrivingForceFlag == 1 && m_defaultDeltaCoefficient < 0.0, + GEOS_ERROR_IF( m_fractureModelType == FractureModelType::Nucleation && m_defaultDeltaCoefficient < 0.0, "delta coefficient must be input and non-negative when the" - " external driving force flag is turned on", + " Nucleation crack model is used", this->getDataContext() ); - - // set results as array default values - this->template getField< fields::solid::criticalFractureEnergy >(). - setApplyDefaultValue( m_defaultCriticalFractureEnergy ); - - this->template getField< fields::solid::tensileStrength >(). - setApplyDefaultValue( m_defaultTensileStrength ); - - this->template getField< fields::solid::compressiveStrength >(). - setApplyDefaultValue( m_defaultCompressiveStrength ); - - this->template getField< fields::solid::deltaCoefficient >(). - setApplyDefaultValue( m_defaultDeltaCoefficient ); } template< typename BASE > @@ -136,7 +134,8 @@ void Damage< BASE >::allocateConstitutiveData( Group & parent, localIndex const m_newDamage.resize( 0, numPts ); m_oldDamage.resize( 0, numPts ); m_damageGrad.resize( 0, numPts, 3 ); - m_strainEnergyDensity.resize( 0, numPts ); + m_crackDrivingForce.resize( 0, numPts ); + m_oldCrackDrivingForce.resize( 0, numPts ); m_volStrain.resize( 0, numPts ); m_extDrivingForce.resize( 0, numPts ); m_biotCoefficient.resize( parent.size() ); @@ -145,6 +144,20 @@ void Damage< BASE >::allocateConstitutiveData( Group & parent, localIndex const m_compressiveStrength.resize( parent.size() ); m_deltaCoefficient.resize( parent.size() ); + // apply the defaults here, after resizing: setApplyDefaultValue() in postInputInitialization() + // runs before the arrays exist, so it has nothing to fill in. + this->template getField< fields::solid::criticalFractureEnergy >(). + setApplyDefaultValue( m_defaultCriticalFractureEnergy ); + + this->template getField< fields::solid::tensileStrength >(). + setApplyDefaultValue( m_defaultTensileStrength ); + + this->template getField< fields::solid::compressiveStrength >(). + setApplyDefaultValue( m_defaultCompressiveStrength ); + + this->template getField< fields::solid::deltaCoefficient >(). + setApplyDefaultValue( m_defaultDeltaCoefficient ); + BASE::allocateConstitutiveData( parent, numPts ); } diff --git a/src/coreComponents/constitutive/solid/Damage.hpp b/src/coreComponents/constitutive/solid/Damage.hpp index b54bb30f372..b222285c77a 100644 --- a/src/coreComponents/constitutive/solid/Damage.hpp +++ b/src/coreComponents/constitutive/solid/Damage.hpp @@ -20,11 +20,9 @@ * * In a phase-field for fracture model, the damage variable affects the Elasticity equation * with the degradation of the stresses. Instead of sigma = C : epsilon, we have sigma = g(d)*C:epsilon, - * where g(d) is the degradation function. This degradation function can either be a quadratic one - * (set LORENTZ 0) or a quasi-quadratic one (set LORENTZ 1). In general, the quadratic one will give you - * brittle fracture behaviour. The quasi-quadratic one, combined with linear dissipation, will give you - * cohesive fracture behaviour, with a user-defined critical stress. If you use quadratic dissipation in - * your damage solver, set QUADRATIC_DISSIPATION to 1. + * where g(d) is the degradation function. The base Damage model uses a quadratic degradation + * function. DamageSpectral uses a quasi-quadratic one which, combined with linear dissipation, + * gives cohesive fracture behaviour with a user-defined critical stress. * * References: * @@ -43,6 +41,7 @@ #ifndef GEOS_CONSTITUTIVE_SOLID_DAMAGE_HPP_ #define GEOS_CONSTITUTIVE_SOLID_DAMAGE_HPP_ +#include "common/format/EnumStrings.hpp" #include "constitutive/solid/SolidBase.hpp" #include "InvariantDecompositions.hpp" #include "ElasticIsotropic.hpp" @@ -52,6 +51,30 @@ namespace geos namespace constitutive { +/// Type of local dissipation function used by the phase-field damage solver. +enum class LocalDissipationOption : integer +{ + Linear, + Quadratic, +}; + +ENUM_STRINGS( LocalDissipationOption, + "Linear", + "Quadratic" ); + +/// Type of crack/fracture model used by the phase-field damage model. +enum class FractureModelType : integer +{ + Brittle, + Cohesive, + Nucleation, +}; + +ENUM_STRINGS( FractureModelType, + "Brittle", + "Cohesive", + "Nucleation" ); + // DAMAGE MODEL UPDATES // // NOTE: This model uses the m_newStress array to represent the stress in an @@ -73,14 +96,16 @@ class DamageUpdates : public UPDATE_BASE DamageUpdates( arrayView2d< real64 > const & inputNewDamage, arrayView2d< real64 > const & inputOldDamage, arrayView3d< real64 > const & inputDamageGrad, - arrayView2d< real64 > const & inputStrainEnergyDensity, + arrayView2d< real64 > const & inputCrackDrivingForce, + arrayView2d< real64 > const & inputOldCrackDrivingForce, arrayView2d< real64 > const & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, arrayView1d< real64 > const & inputCriticalFractureEnergy, real64 const & inputcriticalStrainEnergy, real64 const & inputDegradationLowerLimit, - integer const & inputExtDrivingForceFlag, + FractureModelType const & inputFractureModelType, + LocalDissipationOption const & inputLocalDissipationOption, arrayView1d< real64 > const & inputTensileStrength, arrayView1d< real64 > const & inputCompressiveStrength, arrayView1d< real64 > const & inputDeltaCoefficient, @@ -90,14 +115,16 @@ class DamageUpdates : public UPDATE_BASE m_newDamage( inputNewDamage ), m_oldDamage( inputOldDamage ), m_damageGrad( inputDamageGrad ), - m_strainEnergyDensity( inputStrainEnergyDensity ), + m_crackDrivingForce( inputCrackDrivingForce ), + m_oldCrackDrivingForce( inputOldCrackDrivingForce ), m_volStrain( inputVolumetricStrain ), m_extDrivingForce ( inputExtDrivingForce ), m_lengthScale( inputLengthScale ), m_criticalFractureEnergy( inputCriticalFractureEnergy ), m_criticalStrainEnergy( inputcriticalStrainEnergy ), m_degradationLowerLimit( inputDegradationLowerLimit ), - m_extDrivingForceFlag( inputExtDrivingForceFlag ), + m_fractureModelType( inputFractureModelType ), + m_localDissipationOption( inputLocalDissipationOption ), m_tensileStrength( inputTensileStrength ), m_compressiveStrength( inputCompressiveStrength ), m_deltaCoefficient( inputDeltaCoefficient ), @@ -114,28 +141,33 @@ class DamageUpdates : public UPDATE_BASE using UPDATE_BASE::m_disableInelasticity; - //Standard quadratic degradation functions + //Degradation functions: quadratic (Brittle/Nucleation) or Lorentz-type rational (Cohesive, + //Geelen et al., 2019, CMAME; AT1-only, validated at input-parsing time). inline GEOS_HOST_DEVICE virtual real64 getDegradationValue( localIndex const k, localIndex const q ) const { - real64 pf; + real64 const pf = LvArray::math::max( LvArray::math::min( 1.0, m_newDamage( k, q )), 0.0 ); - if( m_extDrivingForceFlag ) - { - pf = LvArray::math::max( LvArray::math::min( 1.0, m_newDamage( k, q )), 0.0 ); - } - else + switch( m_fractureModelType ) { - pf = m_newDamage( k, q ); + case FractureModelType::Cohesive: + { + real64 const m = 3*m_criticalFractureEnergy[k]/(8*m_lengthScale*m_criticalStrainEnergy); + real64 const p = 1; + return pow( 1 - pf, 2 ) / ( pow( 1 - pf, 2 ) + m*pf*(1 + p*pf) ); + } + case FractureModelType::Brittle: + case FractureModelType::Nucleation: + default: + { + // Set a lower bound tolerance for the degradation + real64 const eps = m_degradationLowerLimit; + return ((1 - eps)*(1 - pf)*(1 - pf) + eps); + } } - - // Set a lower bound tolerance for the degradation - real64 const eps = m_degradationLowerLimit; - - return ((1 - eps)*(1 - pf)*(1 - pf) + eps); } @@ -143,9 +175,19 @@ class DamageUpdates : public UPDATE_BASE GEOS_HOST_DEVICE virtual real64 getDegradationDerivative( localIndex const k, real64 const d ) const { - GEOS_UNUSED_VAR( k ); - - return -2*(1 - d); + switch( m_fractureModelType ) + { + case FractureModelType::Cohesive: + { + real64 const m = 3*m_criticalFractureEnergy[k]/(8*m_lengthScale*m_criticalStrainEnergy); + real64 const p = 1; + return -m*(1 - d)*(1 + (2*p + 1)*d) / pow( pow( 1-d, 2 ) + m*d*(1+p*d), 2 ); + } + case FractureModelType::Brittle: + case FractureModelType::Nucleation: + default: + return -2*(1 - d); + } } @@ -153,9 +195,19 @@ class DamageUpdates : public UPDATE_BASE GEOS_HOST_DEVICE virtual real64 getDegradationSecondDerivative( localIndex const k, real64 const d ) const { - GEOS_UNUSED_VAR( k, d ); - - return 2.0; + switch( m_fractureModelType ) + { + case FractureModelType::Cohesive: + { + real64 const m = 3*m_criticalFractureEnergy[k]/(8*m_lengthScale*m_criticalStrainEnergy); + real64 const p = 1; + return -2*m*( pow( d, 3 )*(2*m*p*p + m*p + 2*p + 1) + pow( d, 2 )*(-3*m*p*p -3*p) + d*(-3*m*p - 3) + (-m+p+2) )/pow( pow( 1-d, 2 ) + m*d*(1+p*d), 3 ); + } + case FractureModelType::Brittle: + case FractureModelType::Nucleation: + default: + return 2.0; + } } //Damage dependence function on fluid pressure terms and its derivatives @@ -249,7 +301,12 @@ class DamageUpdates : public UPDATE_BASE m_volStrain( k, q ) = traceOfStrain; - if( m_extDrivingForceFlag ) + // update crack driving force history variable before stress is degraded below + real64 const sed = SolidBaseUpdates::getStrainEnergyDensity( k, q ); + + m_crackDrivingForce( k, q ) = fmax( sed, m_oldCrackDrivingForce( k, q ) ); + + if( m_fractureModelType == FractureModelType::Nucleation ) { real64 stressP; real64 stressQ; @@ -294,21 +351,11 @@ class DamageUpdates : public UPDATE_BASE } - // TODO: The code below assumes the strain energy density will never be - // evaluated in a non-converged / garbage configuration. - GEOS_HOST_DEVICE - virtual real64 getStrainEnergyDensity( localIndex const k, - localIndex const q ) const override + virtual real64 getCrackDrivingForce( localIndex const k, + localIndex const q ) const { - real64 const sed = SolidBaseUpdates::getStrainEnergyDensity( k, q ); - - if( sed > m_strainEnergyDensity( k, q ) ) - { - m_strainEnergyDensity( k, q ) = sed; - } - - return m_strainEnergyDensity( k, q ); + return m_crackDrivingForce( k, q ); } GEOS_HOST_DEVICE @@ -335,17 +382,25 @@ class DamageUpdates : public UPDATE_BASE virtual real64 getEnergyThreshold( localIndex const k, localIndex const q ) const { - #if LORENTZ - return m_criticalStrainEnergy; - #else - if( m_extDrivingForceFlag ) - return 3*m_criticalFractureEnergy[k]/(16 * m_lengthScale) + 0.5 * m_extDrivingForce( k, q ); - else - return 3*m_criticalFractureEnergy[k]/(16 * m_lengthScale); - - #endif - - + switch( m_fractureModelType ) + { + case FractureModelType::Cohesive: + return m_criticalStrainEnergy; + + case FractureModelType::Nucleation: + return 3*m_criticalFractureEnergy[k]/(16 * m_lengthScale) + 0.5 * m_extDrivingForce( k, q ); + + case FractureModelType::Brittle: + default: + switch( m_localDissipationOption ) + { + case LocalDissipationOption::Linear: + return 3*m_criticalFractureEnergy[k]/(16 * m_lengthScale); + case LocalDissipationOption::Quadratic: + default: + return 0.0; // unused: AT2 (Quadratic) dissipation never reads the energy threshold + } + } } GEOS_HOST_DEVICE @@ -360,6 +415,7 @@ class DamageUpdates : public UPDATE_BASE { ElasticIsotropicUpdates::saveConvergedState( k, q ); m_oldDamage[k][q] = m_newDamage[k][q]; + m_oldCrackDrivingForce[k][q] = m_crackDrivingForce[k][q]; } /// The new damage value on all quadrature points @@ -372,7 +428,10 @@ class DamageUpdates : public UPDATE_BASE arrayView3d< real64 > const m_damageGrad; /// The strain energy density to drive fracture at the quadrature point - arrayView2d< real64 > const m_strainEnergyDensity; + arrayView2d< real64 > const m_crackDrivingForce; + + /// The crack driving force at the last converged state + arrayView2d< real64 > const m_oldCrackDrivingForce; /// The volumetric strain at the quadrature point arrayView2d< real64 > const m_volStrain; @@ -392,8 +451,11 @@ class DamageUpdates : public UPDATE_BASE /// The lower limit of the degradation function real64 const m_degradationLowerLimit; - /// The flag to indicate if the external driving force is used for fracture nucleation - integer const m_extDrivingForceFlag; + /// The type of crack/fracture model (Brittle, Cohesive, or Nucleation) + FractureModelType const m_fractureModelType; + + /// The type of local dissipation function used by the phase-field solver + LocalDissipationOption const m_localDissipationOption; /// A reference view to the tensile strength for each element arrayView1d< real64 > const m_tensileStrength; @@ -437,20 +499,26 @@ class Damage : public BASE arrayView2d< real64 const > getExtDrivingForce() const { return m_extDrivingForce; } + FractureModelType getFractureModelType() const { return m_fractureModelType; } + + LocalDissipationOption getLocalDissipationOption() const { return m_localDissipationOption; } + KernelWrapper createKernelUpdates() const { return BASE::template createDerivedKernelUpdates< KernelWrapper >( m_newDamage.toView(), m_oldDamage.toView(), m_damageGrad.toView(), - m_strainEnergyDensity.toView(), + m_crackDrivingForce.toView(), + m_oldCrackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, m_criticalFractureEnergy.toView(), m_criticalStrainEnergy, m_degradationLowerLimit, - m_extDrivingForceFlag, + m_fractureModelType, + m_localDissipationOption, m_tensileStrength.toView(), m_compressiveStrength.toView(), m_deltaCoefficient.toView(), @@ -467,8 +535,10 @@ class Damage : public BASE static constexpr char const * criticalStrainEnergyString() { return "criticalStrainEnergy"; } /// string/key for degradation lower limit static constexpr char const * degradationLowerLimitString() { return "degradationLowerLimit"; } - // string/key for c_e switch - static constexpr char const * extDrivingForceFlagString() { return "extDrivingForceFlag"; } + /// string/key for the crack/fracture model type + static constexpr char const * fractureModelTypeString() { return "fractureModelType"; } + /// string/key for the local dissipation option + static constexpr char const * localDissipationOptionString() { return "localDissipationOption"; } /// string/key for the default tensile strength static constexpr char const * defaultTensileStrengthString() { return "defaultTensileStrength"; } /// string/key for the default compressive strength @@ -490,7 +560,10 @@ class Damage : public BASE array3d< real64 > m_damageGrad; /// The strain energy density to drive fracture at the quadrature point - array2d< real64 > m_strainEnergyDensity; + array2d< real64 > m_crackDrivingForce; + + /// The crack driving force at the last converged state + array2d< real64 > m_oldCrackDrivingForce; /// The volumetric strain at the quadrature point array2d< real64 > m_volStrain; @@ -510,8 +583,11 @@ class Damage : public BASE /// The lower limit of the degradation function real64 m_degradationLowerLimit; - /// The flag to indicate if the external driving force is used for fracture nucleation - integer m_extDrivingForceFlag; + /// The type of crack/fracture model (Brittle, Cohesive, or Nucleation) + FractureModelType m_fractureModelType; + + /// The type of local dissipation function used by the phase-field solver + LocalDissipationOption m_localDissipationOption; /// The default value of the tensile strength real64 m_defaultTensileStrength; diff --git a/src/coreComponents/constitutive/solid/DamageSpectral.cpp b/src/coreComponents/constitutive/solid/DamageSpectral.cpp index 45d6631ebb8..7eefac31fce 100644 --- a/src/coreComponents/constitutive/solid/DamageSpectral.cpp +++ b/src/coreComponents/constitutive/solid/DamageSpectral.cpp @@ -32,7 +32,21 @@ namespace constitutive template< typename BASE > DamageSpectral< BASE >::DamageSpectral( string const & name, Group * const parent ): Damage< BASE >( name, parent ) -{} +{ + // Preserve historical behavior: DamageSpectral defaults to the Cohesive degradation function. + this->template getWrapper< FractureModelType >( Damage< BASE >::viewKeyStruct::fractureModelTypeString() ). + setApplyDefaultValue( FractureModelType::Cohesive ); +} + +template< typename BASE > +void DamageSpectral< BASE >::postInputInitialization() +{ + Damage< BASE >::postInputInitialization(); + + GEOS_ERROR_IF( this->getFractureModelType() == FractureModelType::Nucleation, + "the Nucleation crack model is not supported with the Spectral split", + this->getDataContext() ); +} typedef DamageSpectral< ElasticIsotropic > DamageSpectralElasticIsotropic; diff --git a/src/coreComponents/constitutive/solid/DamageSpectral.hpp b/src/coreComponents/constitutive/solid/DamageSpectral.hpp index 03a5e3bff73..c6861d40da9 100644 --- a/src/coreComponents/constitutive/solid/DamageSpectral.hpp +++ b/src/coreComponents/constitutive/solid/DamageSpectral.hpp @@ -27,8 +27,6 @@ #include "SolidBase.hpp" #include "SolidModelDiscretizationOpsFullyAnisotropic.hpp" -#define QUADRATIC_DISSIPATION 0 - namespace geos { namespace constitutive @@ -42,21 +40,24 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > DamageSpectralUpdates( arrayView2d< real64 > const & inputNewDamage, arrayView2d< real64 > const & inputOldDamage, arrayView3d< real64 > const & inputDamageGrad, - arrayView2d< real64 > const & inputStrainEnergyDensity, + arrayView2d< real64 > const & inputCrackDrivingForce, + arrayView2d< real64 > const & inputOldCrackDrivingForce, arrayView2d< real64 > const & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, arrayView1d< real64 > const & inputCriticalFractureEnergy, real64 const & inputcriticalStrainEnergy, real64 const & inputDegradationLowerLimit, - int const & inputExtDrivingForceFlag, + FractureModelType const & inputFractureModelType, + LocalDissipationOption const & inputLocalDissipationOption, arrayView1d< real64 > const & inputTensileStrength, arrayView1d< real64 > const & inputCompressiveStrength, arrayView1d< real64 > const & inputDeltaCoefficient, arrayView1d< real64 > const & inputBiotCoefficient, PARAMS && ... baseParams ): - DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputStrainEnergyDensity, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, - inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputExtDrivingForceFlag, + DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputCrackDrivingForce, inputOldCrackDrivingForce, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, + inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputFractureModelType, + inputLocalDissipationOption, inputTensileStrength, inputCompressiveStrength, inputDeltaCoefficient, inputBiotCoefficient, std::forward< PARAMS >( baseParams )... ) {} @@ -65,13 +66,10 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > using DamageUpdates< UPDATE_BASE >::smallStrainUpdate; using DamageUpdates< UPDATE_BASE >::saveConvergedState; - using DamageUpdates< UPDATE_BASE >::getDegradationValue; - using DamageUpdates< UPDATE_BASE >::getDegradationDerivative; - using DamageUpdates< UPDATE_BASE >::getDegradationSecondDerivative; - using DamageUpdates< UPDATE_BASE >::getEnergyThreshold; - using DamageUpdates< UPDATE_BASE >::m_strainEnergyDensity; + using DamageUpdates< UPDATE_BASE >::m_crackDrivingForce; + using DamageUpdates< UPDATE_BASE >::m_oldCrackDrivingForce; using DamageUpdates< UPDATE_BASE >::m_volStrain; using DamageUpdates< UPDATE_BASE >::m_criticalStrainEnergy; using DamageUpdates< UPDATE_BASE >::m_extDrivingForce; @@ -80,7 +78,8 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > using DamageUpdates< UPDATE_BASE >::m_newDamage; using DamageUpdates< UPDATE_BASE >::m_oldDamage; using DamageUpdates< UPDATE_BASE >::m_damageGrad; - using DamageUpdates< UPDATE_BASE >::m_extDrivingForceFlag; + using DamageUpdates< UPDATE_BASE >::m_fractureModelType; + using DamageUpdates< UPDATE_BASE >::m_localDissipationOption; using DamageUpdates< UPDATE_BASE >::m_tensileStrength; using DamageUpdates< UPDATE_BASE >::m_compressiveStrength; using DamageUpdates< UPDATE_BASE >::m_deltaCoefficient; @@ -91,51 +90,6 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > using UPDATE_BASE::m_bulkModulus; // TODO: model below strongly assumes iso elasticity, templating not so useful using UPDATE_BASE::m_shearModulus; - // Lorentz type degradation functions - - inline - GEOS_HOST_DEVICE - virtual real64 getDegradationValue( localIndex const k, - localIndex const q ) const override - { - #if QUADRATIC_DISSIPATION - real64 m = m_criticalFractureEnergy[k]/(2*m_lengthScale*m_criticalStrainEnergy); - #else - real64 m = 3*m_criticalFractureEnergy[k]/(8*m_lengthScale*m_criticalStrainEnergy); - #endif - real64 p = 1; - return pow( 1 - m_newDamage( k, q ), 2 ) /( pow( 1 - m_newDamage( k, q ), 2 ) + m * m_newDamage( k, q ) * (1 + p*m_newDamage( k, q )) ); - } - - - inline - GEOS_HOST_DEVICE - virtual real64 getDegradationDerivative( localIndex const k, real64 const d ) const override - { - #if QUADRATIC_DISSIPATION - real64 m = m_criticalFractureEnergy[k]/(2*m_lengthScale*m_criticalStrainEnergy); - #else - real64 m = 3*m_criticalFractureEnergy[k]/(8*m_lengthScale*m_criticalStrainEnergy); - #endif - real64 p = 1; - return -m*(1 - d)*(1 + (2*p + 1)*d) / pow( pow( 1-d, 2 ) + m*d*(1+p*d), 2 ); - } - - - inline - GEOS_HOST_DEVICE - virtual real64 getDegradationSecondDerivative( localIndex const k, real64 const d ) const override - { - #if QUADRATIC_DISSIPATION - real64 m = m_criticalFractureEnergy[k]/(2*m_lengthScale*m_criticalStrainEnergy); - #else - real64 m = 3*m_criticalFractureEnergy[k]/(8*m_lengthScale*m_criticalStrainEnergy); - #endif - real64 p = 1; - return -2*m*( pow( d, 3 )*(2*m*p*p + m*p + 2*p + 1) + pow( d, 2 )*(-3*m*p*p -3*p) + d*(-3*m*p - 3) + (-m+p+2) )/pow( pow( 1-d, 2 ) + m*d*(1+p*d), 3 ); - } - - GEOS_HOST_DEVICE void smallStrainUpdate( localIndex const k, localIndex const q, @@ -250,10 +204,7 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > real64 const sed = 0.5 * lambda * tracePlus * tracePlus + mu * doubleContraction( positivePartOfStrain, positivePartOfStrain ); - if( sed > m_strainEnergyDensity( k, q ) ) - { - m_strainEnergyDensity( k, q ) = sed; - } + m_crackDrivingForce( k, q ) = fmax( sed, m_oldCrackDrivingForce( k, q ) ); } @@ -268,25 +219,6 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } - - GEOS_HOST_DEVICE - virtual real64 getStrainEnergyDensity( localIndex const k, - localIndex const q ) const override final - { - return m_strainEnergyDensity( k, q ); - } - - - GEOS_HOST_DEVICE - virtual real64 getEnergyThreshold( localIndex const k, - localIndex const q ) const override final - { - GEOS_UNUSED_VAR( k ); - GEOS_UNUSED_VAR( q ); - - return m_criticalStrainEnergy; - } - }; @@ -300,14 +232,16 @@ class DamageSpectral : public Damage< BASE > using Damage< BASE >::m_newDamage; using Damage< BASE >::m_oldDamage; using Damage< BASE >::m_damageGrad; - using Damage< BASE >::m_strainEnergyDensity; + using Damage< BASE >::m_crackDrivingForce; + using Damage< BASE >::m_oldCrackDrivingForce; using Damage< BASE >::m_volStrain; using Damage< BASE >::m_extDrivingForce; using Damage< BASE >::m_criticalFractureEnergy; using Damage< BASE >::m_lengthScale; using Damage< BASE >::m_criticalStrainEnergy; using Damage< BASE >::m_degradationLowerLimit; - using Damage< BASE >::m_extDrivingForceFlag; + using Damage< BASE >::m_fractureModelType; + using Damage< BASE >::m_localDissipationOption; using Damage< BASE >::m_tensileStrength; using Damage< BASE >::m_compressiveStrength; using Damage< BASE >::m_deltaCoefficient; @@ -319,20 +253,23 @@ class DamageSpectral : public Damage< BASE > static string catalogName() { return string( "DamageSpectral" ) + BASE::catalogName(); } virtual string getCatalogName() const override { return catalogName(); } + virtual void postInputInitialization() override; KernelWrapper createKernelUpdates() const { return BASE::template createDerivedKernelUpdates< KernelWrapper >( m_newDamage.toView(), m_oldDamage.toView(), m_damageGrad.toView(), - m_strainEnergyDensity.toView(), + m_crackDrivingForce.toView(), + m_oldCrackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, m_criticalFractureEnergy.toView(), m_criticalStrainEnergy, m_degradationLowerLimit, - m_extDrivingForceFlag, + m_fractureModelType, + m_localDissipationOption, m_tensileStrength.toView(), m_compressiveStrength.toView(), m_deltaCoefficient.toView(), diff --git a/src/coreComponents/constitutive/solid/DamageVolDev.cpp b/src/coreComponents/constitutive/solid/DamageVolDev.cpp index 0aa0e6983c8..09dd714f45a 100644 --- a/src/coreComponents/constitutive/solid/DamageVolDev.cpp +++ b/src/coreComponents/constitutive/solid/DamageVolDev.cpp @@ -34,6 +34,16 @@ DamageVolDev< BASE >::DamageVolDev( string const & name, Group * const parent ): Damage< BASE >( name, parent ) {} +template< typename BASE > +void DamageVolDev< BASE >::postInputInitialization() +{ + Damage< BASE >::postInputInitialization(); + + GEOS_ERROR_IF( this->getFractureModelType() == FractureModelType::Nucleation, + "the Nucleation crack model is not supported with the Vol-dev split", + this->getDataContext() ); +} + typedef DamageVolDev< ElasticIsotropic > DamageVolDevElasticIsotropic; REGISTER_CATALOG_ENTRY( ConstitutiveBase, DamageVolDevElasticIsotropic, string const &, Group * const ) diff --git a/src/coreComponents/constitutive/solid/DamageVolDev.hpp b/src/coreComponents/constitutive/solid/DamageVolDev.hpp index 951eded451f..1c990aee9b2 100644 --- a/src/coreComponents/constitutive/solid/DamageVolDev.hpp +++ b/src/coreComponents/constitutive/solid/DamageVolDev.hpp @@ -38,21 +38,24 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > DamageVolDevUpdates( arrayView2d< real64 > const & inputNewDamage, arrayView2d< real64 > const & inputOldDamage, arrayView3d< real64 > const & inputDamageGrad, - arrayView2d< real64 > const & inputStrainEnergyDensity, + arrayView2d< real64 > const & inputCrackDrivingForce, + arrayView2d< real64 > const & inputOldCrackDrivingForce, arrayView2d< real64 > const & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, arrayView1d< real64 > const & inputCriticalFractureEnergy, real64 const & inputcriticalStrainEnergy, real64 const & inputDegradationLowerLimit, - int const & inputExtDrivingForceFlag, + FractureModelType const & inputFractureModelType, + LocalDissipationOption const & inputLocalDissipationOption, arrayView1d< real64 > const & inputTensileStrength, arrayView1d< real64 > const & inputCompressiveStrength, arrayView1d< real64 > const & inputDeltaCoefficient, arrayView1d< real64 > const & inputBiotCoefficient, PARAMS && ... baseParams ): - DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputStrainEnergyDensity, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, - inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputExtDrivingForceFlag, + DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputCrackDrivingForce, inputOldCrackDrivingForce, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, + inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputFractureModelType, + inputLocalDissipationOption, inputTensileStrength, inputCompressiveStrength, inputDeltaCoefficient, inputBiotCoefficient, std::forward< PARAMS >( baseParams )... ) {} @@ -63,11 +66,9 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > using DamageUpdates< UPDATE_BASE >::saveConvergedState; using DamageUpdates< UPDATE_BASE >::getDegradationValue; - using DamageUpdates< UPDATE_BASE >::getDegradationDerivative; - using DamageUpdates< UPDATE_BASE >::getDegradationSecondDerivative; - using DamageUpdates< UPDATE_BASE >::getEnergyThreshold; - using DamageUpdates< UPDATE_BASE >::m_strainEnergyDensity; + using DamageUpdates< UPDATE_BASE >::m_crackDrivingForce; + using DamageUpdates< UPDATE_BASE >::m_oldCrackDrivingForce; using DamageUpdates< UPDATE_BASE >::m_volStrain; using DamageUpdates< UPDATE_BASE >::m_criticalStrainEnergy; using DamageUpdates< UPDATE_BASE >::m_extDrivingForce; @@ -76,7 +77,8 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > using DamageUpdates< UPDATE_BASE >::m_newDamage; using DamageUpdates< UPDATE_BASE >::m_oldDamage; using DamageUpdates< UPDATE_BASE >::m_damageGrad; - using DamageUpdates< UPDATE_BASE >::m_extDrivingForceFlag; + using DamageUpdates< UPDATE_BASE >::m_fractureModelType; + using DamageUpdates< UPDATE_BASE >::m_localDissipationOption; using DamageUpdates< UPDATE_BASE >::m_tensileStrength; using DamageUpdates< UPDATE_BASE >::m_compressiveStrength; using DamageUpdates< UPDATE_BASE >::m_deltaCoefficient; @@ -138,7 +140,6 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > stress ); // update strain energy density - // TODO: refactor as a proper history variable update. the code below doesn't allow for rewinds. real64 sed = 0.5 * (stressQ * devStrain) / factor; @@ -147,18 +148,7 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > sed += 0.5 * (stressP * volStrain) / factor; } - if( sed > m_strainEnergyDensity( k, q ) ) - { - m_strainEnergyDensity( k, q ) = sed; - } - } - - - GEOS_HOST_DEVICE - virtual real64 getStrainEnergyDensity( localIndex const k, - localIndex const q ) const override final - { - return m_strainEnergyDensity( k, q ); + m_crackDrivingForce( k, q ) = fmax( sed, m_oldCrackDrivingForce( k, q ) ); } }; @@ -174,14 +164,16 @@ class DamageVolDev : public Damage< BASE > using Damage< BASE >::m_newDamage; using Damage< BASE >::m_oldDamage; using Damage< BASE >::m_damageGrad; - using Damage< BASE >::m_strainEnergyDensity; + using Damage< BASE >::m_crackDrivingForce; + using Damage< BASE >::m_oldCrackDrivingForce; using Damage< BASE >::m_volStrain; using Damage< BASE >::m_extDrivingForce; using Damage< BASE >::m_criticalFractureEnergy; using Damage< BASE >::m_lengthScale; using Damage< BASE >::m_criticalStrainEnergy; using Damage< BASE >::m_degradationLowerLimit; - using Damage< BASE >::m_extDrivingForceFlag; + using Damage< BASE >::m_fractureModelType; + using Damage< BASE >::m_localDissipationOption; using Damage< BASE >::m_tensileStrength; using Damage< BASE >::m_compressiveStrength; using Damage< BASE >::m_deltaCoefficient; @@ -193,20 +185,23 @@ class DamageVolDev : public Damage< BASE > static string catalogName() { return string( "DamageVolDev" ) + BASE::catalogName(); } virtual string getCatalogName() const override { return catalogName(); } + virtual void postInputInitialization() override; KernelWrapper createKernelUpdates() const { return BASE::template createDerivedKernelUpdates< KernelWrapper >( m_newDamage.toView(), m_oldDamage.toView(), m_damageGrad.toView(), - m_strainEnergyDensity.toView(), + m_crackDrivingForce.toView(), + m_oldCrackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, m_criticalFractureEnergy.toView(), m_criticalStrainEnergy, m_degradationLowerLimit, - m_extDrivingForceFlag, + m_fractureModelType, + m_localDissipationOption, m_tensileStrength.toView(), m_compressiveStrength.toView(), m_deltaCoefficient.toView(), diff --git a/src/coreComponents/constitutive/solid/SolidFields.hpp b/src/coreComponents/constitutive/solid/SolidFields.hpp index 8ed466eadc7..bfb1eda4c51 100644 --- a/src/coreComponents/constitutive/solid/SolidFields.hpp +++ b/src/coreComponents/constitutive/solid/SolidFields.hpp @@ -154,13 +154,21 @@ DECLARE_FIELD( damageGrad, WRITE_AND_READ, "Material damage gradient" ); -DECLARE_FIELD( strainEnergyDensity, - "strainEnergyDensity", +DECLARE_FIELD( crackDrivingForce, + "crackDrivingForce", array2d< real64 >, 0, LEVEL_0, WRITE_AND_READ, - "Material strain energy density" ); + "Material crack driving force" ); + +DECLARE_FIELD( oldCrackDrivingForce, + "oldCrackDrivingForce", + array2d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Material crack driving force at the last converged state" ); DECLARE_FIELD( volStrain, "volStrain", diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp index 751dd7dea77..9ca96a93c85 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp @@ -20,11 +20,9 @@ #include "PhaseFieldFractureSolver.hpp" -#include "discretizationMethods/NumericalMethodsManager.hpp" -#include "fieldSpecification/TractionBoundaryCondition.hpp" +#include "constitutive/ConstitutivePassThru.hpp" #include "finiteElement/FiniteElementDispatch.hpp" #include "mesh/DomainPartition.hpp" -#include "mesh/utilities/ComputationalGeometry.hpp" namespace geos { @@ -37,11 +35,6 @@ PhaseFieldFractureSolver::PhaseFieldFractureSolver( const string & name, Base( name, parent ) {} -PhaseFieldFractureSolver::~PhaseFieldFractureSolver() -{ - // TODO Auto-generated destructor stub -} - void PhaseFieldFractureSolver::postInputInitialization() { Base::postInputInitialization(); @@ -52,57 +45,61 @@ void PhaseFieldFractureSolver::postInputInitialization() void PhaseFieldFractureSolver::mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) { + GEOS_MARK_FUNCTION; + + if( solverType == static_cast< integer >( SolverType::Damage ) ) + { + mapDamageToQuadrature( domain ); + } +} +void PhaseFieldFractureSolver::mapDamageToQuadrature( DomainPartition & domain ) +{ GEOS_MARK_FUNCTION; - if( solverType == static_cast< integer >( SolverType::Damage ) ) + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) { - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { - NodeManager & nodeManager = mesh.getNodeManager(); + NodeManager & nodeManager = mesh.getNodeManager(); - string const & damageFieldName = damageSolver()->getFieldName(); + string const & damageFieldName = damageSolver()->getFieldName(); + string const & discretizationName = damageSolver()->getDiscretizationName(); - string const & discretizationName = damageSolver()->getDiscretizationName(); + arrayView1d< real64 const > const nodalDamage = nodeManager.getReference< array1d< real64 > >( damageFieldName ); - //should get reference to damage field here. - arrayView1d< real64 const > const nodalDamage = nodeManager.getReference< array1d< real64 > >( damageFieldName ); + ElementRegionManager & elemManager = mesh.getElemManager(); - ElementRegionManager & elemManager = mesh.getElemManager(); + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [discretizationName, nodalDamage] + ( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + string const & solidModelName = elementSubRegion.getReference< string >( SolidMechanicsLagrangianFEM::viewKeyStruct::solidMaterialNamesString()); + constitutive::SolidBase & + solidModel = elementSubRegion.getConstitutiveModel< constitutive::SolidBase >( solidModelName ); - // begin region loop - elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [discretizationName, nodalDamage] - ( localIndex const, - CellElementSubRegion & elementSubRegion ) + ConstitutivePassThru< DamageBase >::execute( solidModel, [&elementSubRegion, discretizationName, nodalDamage]( auto & damageModel ) { - string const & solidModelName = elementSubRegion.getReference< string >( SolidMechanicsLagrangianFEM::viewKeyStruct::solidMaterialNamesString()); - constitutive::SolidBase & - solidModel = elementSubRegion.getConstitutiveModel< constitutive::SolidBase >( solidModelName ); + using CONSTITUTIVE_TYPE = TYPEOFREF( damageModel ); + typename CONSTITUTIVE_TYPE::KernelWrapper constitutiveUpdate = damageModel.createKernelUpdates(); - ConstitutivePassThru< DamageBase >::execute( solidModel, [&elementSubRegion, discretizationName, nodalDamage]( auto & damageModel ) - { - using CONSTITUTIVE_TYPE = TYPEOFREF( damageModel ); - typename CONSTITUTIVE_TYPE::KernelWrapper constitutiveUpdate = damageModel.createKernelUpdates(); - - arrayView2d< real64 > const damageFieldOnMaterial = constitutiveUpdate.m_newDamage; - arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemToNodes = elementSubRegion.nodeList(); + arrayView2d< real64 > const damageFieldOnMaterial = constitutiveUpdate.m_newDamage; + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemToNodes = elementSubRegion.nodeList(); - finiteElement::FiniteElementBase const & - fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( discretizationName ); + finiteElement::FiniteElementBase const & + fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( discretizationName ); - finiteElement::FiniteElementDispatchHandler< ALL_FE_TYPES >::dispatch3D( fe, [=, &elementSubRegion] ( auto & finiteElement ) - { - using FE_TYPE = TYPEOFREF( finiteElement ); + finiteElement::FiniteElementDispatchHandler< ALL_FE_TYPES >::dispatch3D( fe, [=, &elementSubRegion] ( auto & finiteElement ) + { + using FE_TYPE = TYPEOFREF( finiteElement ); - DamageInterpolationKernel< FE_TYPE > interpolationKernel( elementSubRegion ); + DamageInterpolationKernel< FE_TYPE > interpolationKernel( elementSubRegion ); - interpolationKernel.interpolateDamage( elemToNodes, nodalDamage, damageFieldOnMaterial ); - } ); + interpolationKernel.interpolateDamage( elemToNodes, nodalDamage, damageFieldOnMaterial ); } ); } ); } ); - } + } ); } REGISTER_CATALOG_ENTRY( PhysicsSolverBase, PhaseFieldFractureSolver, string const &, Group * const ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.hpp index a1aae598694..e493d213381 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.hpp @@ -39,14 +39,20 @@ class PhaseFieldFractureSolver : public CoupledSolver< SolidMechanicsLagrangianF using Base::m_rhs; using Base::m_solution; + /** + * @brief main constructor for PhaseFieldFractureSolver objects + * @param name the name of this instantiation of PhaseFieldFractureSolver in the repository + * @param parent the parent group of this instantiation of PhaseFieldFractureSolver + */ PhaseFieldFractureSolver( const string & name, Group * const parent ); - ~PhaseFieldFractureSolver() override; + /// Destructor for the class + ~PhaseFieldFractureSolver() override = default; /** - * @brief name of the node manager in the object catalog - * @return string that contains the catalog name to generate a new NodeManager object through the object catalog. + * @brief name of the solver in the object catalog + * @return string that contains the catalog name to generate a new PhaseFieldFractureSolver object through the object catalog. */ static string catalogName() { @@ -78,20 +84,28 @@ class PhaseFieldFractureSolver : public CoupledSolver< SolidMechanicsLagrangianF } /** - * @brief accessor for the pointer to the flow solver - * @return a pointer to the flow solver + * @brief accessor for the pointer to the damage solver + * @return a pointer to the damage solver */ PhaseFieldDamageFEM * damageSolver() const { return std::get< toUnderlying( SolverType::Damage ) >( m_solvers ); } - virtual void mapSolutionBetweenSolvers( DomainPartition & Domain, integer const idx ) override final; + virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override final; protected: virtual void initializePostInitialConditionsPreSubGroups() override final {} +private: + + /** + * @brief interpolate the nodal damage field onto the quadrature points of the solid constitutive model + * @param[in] domain the domain partition + */ + void mapDamageToQuadrature( DomainPartition & domain ); + }; @@ -120,9 +134,6 @@ struct DamageInterpolationKernel for( localIndex a = 0; a < numNodesPerElement; ++a ) { damageFieldOnMaterial( k, q ) += N[a] * nodalDamage[elemToNodes( k, a )]; - //solution is probably not going to work because the solution of the coupled solver - //has both damage and displacements. Using the damageResult field from the Damage solver - //is probably better } } diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.cpp index 06bb5259f06..c3e44cfc1ad 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.cpp @@ -22,9 +22,6 @@ #include "fieldSpecification/TractionBoundaryCondition.hpp" #include "mesh/DomainPartition.hpp" -#include "mesh/MeshForLoopInterface.hpp" -#include "mesh/utilities/ComputationalGeometry.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" namespace geos { @@ -52,11 +49,6 @@ void PhaseFieldPoromechanicsSolver::postInputInitialization() getNonlinearSolverParameters().m_couplingType = NonlinearSolverParameters::CouplingType::Sequential; } -PhaseFieldPoromechanicsSolver::~PhaseFieldPoromechanicsSolver() -{ - // TODO Auto-generated destructor stub -} - void PhaseFieldPoromechanicsSolver::mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) { if( solverType == static_cast< integer >( SolverType::Damage ) ) @@ -74,12 +66,10 @@ void PhaseFieldPoromechanicsSolver::mapSolutionBetweenSolvers( DomainPartition & string const & discretizationName = damageSolver()->getDiscretizationName(); - //should get reference to damage field here. arrayView1d< real64 const > const nodalDamage = nodeManager.getReference< array1d< real64 > >( damageFieldName ); ElementRegionManager & elemManager = mesh.getElemManager(); - // begin region loop elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [discretizationName, xNodes, nodalDamage] ( localIndex const, CellElementSubRegion & elementSubRegion ) @@ -114,7 +104,7 @@ void PhaseFieldPoromechanicsSolver::mapSolutionBetweenSolvers( DomainPartition & } else if( solverType == static_cast< integer >( SolverType::Poromechanics ) ) { - poromechancisSolver()->flowSolver()->updatePressureGradient( domain ); + poromechanicsSolver()->flowSolver()->updatePressureGradient( domain ); } } diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp index 75de5e22b49..cdaa953aa0f 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp @@ -24,7 +24,6 @@ #include "physicsSolvers/multiphysics/CoupledSolver.hpp" #include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" #include "physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" namespace geos { @@ -43,7 +42,7 @@ class PhaseFieldPoromechanicsSolver : public CoupledSolver< SinglePhasePoromecha PhaseFieldPoromechanicsSolver( const string & name, Group * const parent ); - ~PhaseFieldPoromechanicsSolver() override; + ~PhaseFieldPoromechanicsSolver() override = default; /** * @brief name of the node manager in the object catalog @@ -73,7 +72,7 @@ class PhaseFieldPoromechanicsSolver : public CoupledSolver< SinglePhasePoromecha * @brief accessor for the pointer to the poromechanics solver * @return a pointer to the poromechanics solver */ - SinglePhasePoromechanics< SinglePhaseBase > * poromechancisSolver() const + SinglePhasePoromechanics< SinglePhaseBase > * poromechanicsSolver() const { return std::get< toUnderlying( SolverType::Poromechanics ) >( m_solvers ); } @@ -89,8 +88,6 @@ class PhaseFieldPoromechanicsSolver : public CoupledSolver< SinglePhasePoromecha virtual void mapSolutionBetweenSolvers( DomainPartition & Domain, integer const idx ) override final; - void mapDamageAndGradientToQuadrature( DomainPartition & domain ); - void applyDamageOnTractionBC( DomainPartition & domain ); protected: @@ -162,4 +159,4 @@ struct DamageAndDamageGradientInterpolationKernel } /* namespace geos */ -#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_PhaseFieldPoromechanicsSOLVER_HPP_ */ +#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_PHASEFIELDPOROMECHANICSSOLVER_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp index bfe09b66e6d..bdd978ad0e0 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp @@ -14,40 +14,26 @@ */ /** - * @file PhaseFieldDamageFEM.hpp + * @file PhaseFieldDamageFEM.cpp */ #include "PhaseFieldDamageFEM.hpp" #include "PhaseFieldDamageFEMKernels.hpp" #include "PhaseFieldPressurizedDamageFEMKernels.hpp" -#include -#include - -#include "common/TimingMacros.hpp" -#include "dataRepository/Group.hpp" -#include "mesh/mpiCommunications/CommunicationTools.hpp" -#include "mesh/mpiCommunications/NeighborCommunicator.hpp" #include "common/DataTypes.hpp" -#include "constitutive/ConstitutivePassThru.hpp" +#include "common/TimingMacros.hpp" #include "constitutive/solid/Damage.hpp" -#include "physicsSolvers/LogLevelsInfo.hpp" #include "constitutive/solid/SolidBase.hpp" #include "fieldSpecification/FieldSpecificationImpl.hpp" #include "finiteElement/FiniteElementDiscretization.hpp" -#include "finiteElement/Kinematics.h" - #include "mesh/DomainPartition.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" +#include "physicsSolvers/LogLevelsInfo.hpp" namespace geos { -namespace dataRepository -{ -namespace keys -{} -} // namespace dataRepository - using namespace dataRepository; using namespace constitutive; @@ -67,10 +53,6 @@ PhaseFieldDamageFEM::PhaseFieldDamageFEM( const string & name, setInputFlag( InputFlags::REQUIRED ). setDescription( "name of field variable" ); - registerWrapper( viewKeyStruct::localDissipationOptionString(), &m_localDissipationOption ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "Type of local dissipation function. Can be Linear or Quadratic" ); - registerWrapper( viewKeyStruct::irreversibilityFlagString(), &m_irreversibilityFlag ). setApplyDefaultValue( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -89,10 +71,7 @@ PhaseFieldDamageFEM::PhaseFieldDamageFEM( const string & name, addLogLevel< logInfo::ResidualNorm >(); } -PhaseFieldDamageFEM::~PhaseFieldDamageFEM() -{ - // TODO Auto-generated destructor stub -} +PhaseFieldDamageFEM::~PhaseFieldDamageFEM() = default; void PhaseFieldDamageFEM::registerDataOnMesh( Group & meshBodies ) { @@ -122,21 +101,6 @@ void PhaseFieldDamageFEM::registerDataOnMesh( Group & meshBodies ) } ); } -void PhaseFieldDamageFEM::postInputInitialization() -{ - PhysicsSolverBase::postInputInitialization(); - - // Set basic parameters for solver - // m_linearSolverParameters.logLevel = 0; - // m_linearSolverParameters.solverType = "gmres"; - // m_linearSolverParameters.krylov.tolerance = 1e-8; - // m_linearSolverParameters.krylov.maxIterations = 250; - // m_linearSolverParameters.krylov.maxRestart = 250; - // m_linearSolverParameters.preconditionerType = "amg"; - // m_linearSolverParameters.amg.smootherType = "gaussSeidel"; - // m_linearSolverParameters.amg.coarseType = "direct"; -} - real64 PhaseFieldDamageFEM::solverStep( real64 const & time_n, real64 const & dt, const int cycleNumber, @@ -197,52 +161,21 @@ void PhaseFieldDamageFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM( time_n arrayView1d< real64 > const & localRhs ) { GEOS_MARK_FUNCTION; + + localMatrix.zero(); + localRhs.zero(); + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, string_array const & regionNames ) { NodeManager & nodeManager = mesh.getNodeManager(); - arrayView1d< globalIndex const > const & dofIndex = nodeManager.getReference< array1d< globalIndex > >( dofManager.getKey( m_fieldName ) ); - - // Initialize all entries to zero -#if 1 // Andre...this is the new code - localMatrix.zero(); - localRhs.zero(); + arrayView1d< globalIndex const > const & dofIndex = + nodeManager.getReference< array1d< globalIndex > >( dofManager.getKey( m_fieldName ) ); - auto const localDissipation = m_localDissipationOption == LocalDissipation::Linear ? - PhaseFieldDamageKernelLocalDissipation::Linear : - PhaseFieldDamageKernelLocalDissipation::Quadratic; - - if( m_fracturePressureTermFlag ) - { - PhaseFieldPressurizedDamageKernelFactory kernelFactory( dofIndex, - dofManager.rankOffset(), - localMatrix, - localRhs, - dt, - m_fieldName, - localDissipation ); - - finiteElement:: - regionBasedKernelApplication< parallelDevicePolicy<>, - constitutive::DamageBase, - CellElementSubRegion >( mesh, - regionNames, - this->getDiscretizationName(), - viewKeyStruct::solidModelNamesString(), - kernelFactory ); - } - else + auto const launchKernels = [&]( auto kernelFactory ) { - PhaseFieldDamageKernelFactory kernelFactory( dofIndex, - dofManager.rankOffset(), - localMatrix, - localRhs, - dt, - m_fieldName, - localDissipation ); - finiteElement:: regionBasedKernelApplication< parallelDevicePolicy<>, constitutive::DamageBase, @@ -251,177 +184,26 @@ void PhaseFieldDamageFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM( time_n this->getDiscretizationName(), viewKeyStruct::solidModelNamesString(), kernelFactory ); - } - -#else // this has your changes to the old base code - matrix.zero(); - rhs.zero(); + }; - matrix.open(); - rhs.open(); - - // begin region loop - for( localIndex er = 0; er < elemManager->numRegions(); ++er ) + if( m_fracturePressureTermFlag ) { - ElementRegionBase & elementRegion = elemManager->GetRegion( er ); - - elementRegion.forElementSubRegionsIndex< CellElementSubRegion >( [&]( localIndex const GEOS_UNUSED_PARAM( esr ), - CellElementSubRegion & elementSubRegion ) - { - - constitutive::ConstitutiveBase & - solidModel = elementSubRegion.getConstitutiveModel< constitutive::ConstitutiveBase >( m_solidModelName ); - - constitutive::ConstitutivePassThru< constitutive::DamageBase >::execute( solidModel, - [&]( auto & damageModel ) - { - using CONSTITUTIVE_TYPE = TYPEOFREF( damageModel ); - typename CONSTITUTIVE_TYPE::KernelWrapper constitutiveUpdate = damageModel.createKernelUpdates(); - - arrayView4d< real64 const > const & - dNdX = elementSubRegion.dNdX(); - - arrayView2d< real64 const > const & - detJ = elementSubRegion.detJ(); - - localIndex const numNodesPerElement = elementSubRegion.numNodesPerElement(); - arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemNodes = elementSubRegion.nodeList(); - - // arrayView1d const & - // coeff = elementSubRegion.getReference >(viewKeyStruct::coeffNameString()); - - globalIndex_array elemDofIndex( numNodesPerElement ); - real64_array element_rhs( numNodesPerElement ); - real64_array2d element_matrix( numNodesPerElement, numNodesPerElement ); - - arrayView1d< integer const > const & elemGhostRank = elementSubRegion.ghostRank(); - std::unique_ptr< FiniteElementBase > finiteElement = feDiscretization->getFiniteElement( elementSubRegion.getElementType() ); - localIndex const n_q_points = finiteElement->n_quadrature_points(); - - //real64 ell = m_lengthScale; //phase-field length scale - real64 ell = constitutiveUpdate.getRegularizationLength(); - //real64 Gc = m_criticalFractureEnergy; //energy release rate - real64 Gc = constitutiveUpdate.getCriticalFractureEnergy(); - - real64 threshold = constitutiveUpdate.getEnergyThreshold();//elastic energy threshold - use when Local Dissipation is linear - - arrayView1d< real64 > const & nodalDamage = nodeManager.getReference< array1d< real64 > >( m_fieldName ); - //real64 diffusion = 1.0; - // begin element loop, skipping ghost elements - for( localIndex k = 0; k < elementSubRegion.size(); ++k ) - { - if( elemGhostRank[k] < 0 ) - { - element_rhs = 0.0; - element_matrix = 0.0; - for( localIndex q = 0; q < n_q_points; ++q ) - { - real64 const strainEnergyDensity = constitutiveUpdate.calculateStrainEnergyDensity( k, q ); - real64 D = 0; //max between threshold and - // Elastic energy - if( m_localDissipationOption == LocalDissipation::Linear ) - { - D = std::max( threshold, strainEnergyDensity ); - //D = max(strainEnergy(k,q), strainEnergy(k,q));//debbuging line - remove after testing - } - //Interpolate d and grad_d - - real64 qp_damage = 0.0; - R1Tensor qp_grad_damage; - R1Tensor temp; - for( localIndex a = 0; a < numNodesPerElement; ++a ) - { - qp_damage += finiteElement->value( a, q ) * nodalDamage[elemNodes( k, a )]; - temp = dNdX[k][q][a]; - temp *= nodalDamage[elemNodes( k, a )]; - qp_grad_damage += temp; - - } - //std::cout << "Damage: " << qp_damage <value( a, q ); - //element_rhs(a) += detJ[k][q] * Na * myFunc(Xq, Yq, Zq); //older reaction diffusion solver - if( m_localDissipationOption == LocalDissipation::Linear ) - { - // element_rhs( a ) += detJ[k][q] * (Na * (ell * D - 3 * Gc / 16 )/ Gc - - // 0.375*pow( ell, 2 ) * LvArray::tensorOps::AiBi<3>( qp_grad_damage, dNdX[k][q][a] ) - // - - // (ell * D/Gc) * Na * qp_damage); - - element_rhs( a ) += detJ[k][q] * ( -3 * Na / 16 - - 0.375*pow( ell, 2 ) * LvArray::tensorOps::AiBi< 3 >( qp_grad_damage, dNdX[k][q][a] ) - - (0.5 * ell * D/Gc) * Na * constitutiveUpdate.GetDegradationDerivative( qp_damage )); - - } - else - { - // element_rhs( a ) += detJ[k][q] * (Na * (2 * ell) * strainEnergyDensity / Gc - - // (pow( ell, 2 ) * LvArray::tensorOps::AiBi<3>( qp_grad_damage, dNdX[k][q][a] ) + - // Na * qp_damage * (1 + 2 * ell*strainEnergyDensity/Gc)) ); - - - element_rhs( a ) -= detJ[k][q] * (Na * qp_damage + - (pow( ell, 2 ) * LvArray::tensorOps::AiBi< 3 >( qp_grad_damage, dNdX[k][q][a] ) + - Na * constitutiveUpdate.GetDegradationDerivative( qp_damage ) * (ell*strainEnergyDensity/Gc)) ); - } - - for( localIndex b = 0; b < numNodesPerElement; ++b ) - { - real64 Nb = finiteElement->value( b, q ); - if( m_localDissipationOption == LocalDissipation::Linear ) - { - // element_matrix( a, b ) -= detJ[k][q] * - // (0.375*pow( ell, 2 ) * LvArray::tensorOps::AiBi<3>( dNdX[k][q][a], dNdX[k][q][b] ) + - // (ell * D/Gc) * Na * Nb); - // - element_matrix( a, b ) -= detJ[k][q] * - (0.375*pow( ell, 2 ) * LvArray::tensorOps::AiBi< 3 >( dNdX[k][q][a], dNdX[k][q][b] ) + - (0.5 * ell * D/Gc) * constitutiveUpdate.GetDegradationSecondDerivative( qp_damage ) * Na * Nb); - - } - else - { - // element_matrix( a, b ) -= detJ[k][q] * - // ( pow( ell, 2 ) * LvArray::tensorOps::AiBi<3>( dNdX[k][q][a], dNdX[k][q][b] ) + - // Na * Nb * (1 + 2 * ell*strainEnergyDensity/Gc ) - // ); - - element_matrix( a, b ) -= detJ[k][q] * - ( pow( ell, 2 ) * LvArray::tensorOps::AiBi< 3 >( dNdX[k][q][a], dNdX[k][q][b] ) + - Na * Nb * (1 + constitutiveUpdate.GetDegradationSecondDerivative( qp_damage ) * ell*strainEnergyDensity/Gc ) - ); - } - } - } - } - matrix.add( elemDofIndex, elemDofIndex, element_matrix ); - rhs.add( elemDofIndex, element_rhs ); - } - } - } ); - } ); + launchKernels( PhaseFieldPressurizedDamageKernelFactory( dofIndex, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + m_fieldName ) ); } - matrix.close(); - rhs.close(); - - if( getLogLevel() == 2 ) + else { - GEOS_LOG_RANK_0( "After PhaseFieldDamageFEM::AssembleSystem" ); - GEOS_LOG_RANK_0( "\nJacobian:\n" ); - std::cout << matrix; - GEOS_LOG_RANK_0( "\nResidual:\n" ); - std::cout << rhs; + launchKernels( PhaseFieldDamageKernelFactory( dofIndex, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + m_fieldName ) ); } - - if( getLogLevel() >= 3 ) - { - NonlinearSolverParameters & solverParams = getNonlinearSolverParameters(); - integer newtonIter = solverParams.m_numNewtonIterations; -#endif } ); } @@ -481,24 +263,6 @@ void PhaseFieldDamageFEM::applyBoundaryConditions( GEOS_LOG_RANK_0( "\nResidual:\n" ); std::cout << localRhs; } -// -// if( getLogLevel() >= 3 ) -// { -// NonlinearSolverParameters & solverParams = getNonlinearSolverParameters(); -// integer newtonIter = solverParams.m_numNewtonIterations; -// -// string filename_mat = "matrix_bc_" + std::to_string( time_n ) + "_" + -// std::to_string( newtonIter ) + ".mtx"; -// matrix.write( filename_mat ); -// -// string filename_rhs = "rhs_bc_" + std::to_string( time_n ) + "_" + -// std::to_string( newtonIter ) + ".mtx"; -// rhs.write( filename_rhs ); -// -// GEOS_LOG_RANK_0( "After PhaseFieldDamageFEM::applyBoundaryConditions" ); -// GEOS_LOG_RANK_0( "Jacobian: written to " << filename_mat ); -// GEOS_LOG_RANK_0( "Residual: written to " << filename_rhs ); -// } } real64 @@ -527,7 +291,8 @@ PhaseFieldDamageFEM::calculateResidualNorm( real64 const & GEOS_UNUSED_PARAM( ti forAll< parallelDevicePolicy<> >( nodeManager.size(), [localRhs, localSum, dofNumber, rankOffset, ghostRank] GEOS_HOST_DEVICE ( localIndex const k ) { - if( ghostRank[k] < 0 ) + // nodes outside the target regions carry no damage dof + if( ghostRank[k] < 0 && dofNumber[k] >= 0 ) { localIndex const localRow = LvArray::integerConversion< localIndex >( dofNumber[k] - rankOffset ); localSum += localRhs[localRow] * localRhs[localRow]; @@ -632,7 +397,7 @@ void PhaseFieldDamageFEM::applyIrreversibilityConstraint( DofManager const & dof globalIndex const rankOffSet = dofManager.rankOffset(); - real64 const damangeUpperBound = m_damageUpperBound; + real64 const damageUpperBound = m_damageUpperBound; forAll< parallelDevicePolicy<> >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const nodeIndex ) { @@ -642,7 +407,7 @@ void PhaseFieldDamageFEM::applyIrreversibilityConstraint( DofManager const & dof { real64 const damageAtNode = nodalDamage[nodeIndex]; - if( damageAtNode >= damangeUpperBound ) + if( damageAtNode >= damageUpperBound ) { // Specify the contribution to rhs diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp index d4d0bb5ad93..3208e6f9501 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp @@ -25,11 +25,6 @@ #include "fieldSpecification/FieldSpecificationManager.hpp" #include "physicsSolvers/PhysicsSolverBase.hpp" -struct stabledt -{ - double m_maxdt; -}; - namespace geos { namespace dataRepository @@ -141,16 +136,9 @@ class PhaseFieldDamageFEM : public PhysicsSolverBase ExplicitTransient }; - enum class LocalDissipation - { - Linear, - Quadratic, - }; - struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct { static constexpr char const * coeffNameString() { return "coeffField"; } - static constexpr char const * localDissipationOptionString() { return "localDissipation"; } static constexpr char const * irreversibilityFlagString() { return "irreversibilityFlag"; } static constexpr char const * damageUpperBoundString() { return "damageUpperBound"; } static constexpr char const * fracturePressureTermFlagString() { return "fracturePressureTermFlag"; } @@ -160,42 +148,22 @@ class PhaseFieldDamageFEM : public PhysicsSolverBase dataRepository::ViewKey fieldVarName = { "fieldName" }; } PhaseFieldDamageFEMViewKeys; - inline ParallelVector const * getSolution() const - { - return &m_solution; - } - - inline globalIndex getSize() const - { - return m_matrix.numGlobalRows(); - } - string const & getFieldName() const { return m_fieldName; } -protected: - virtual void postInputInitialization() override final; - private: string m_fieldName; - stabledt m_stabledt; TimeIntegrationOption m_timeIntegrationOption; - LocalDissipation m_localDissipationOption; integer m_irreversibilityFlag; real64 m_damageUpperBound; integer m_fracturePressureTermFlag; array1d< real64 > m_coeff; - - PhaseFieldDamageFEM(); }; /// Declare strings associated with enumeration values. -ENUM_STRINGS( PhaseFieldDamageFEM::LocalDissipation, - "Linear", - "Quadratic" ); ENUM_STRINGS( PhaseFieldDamageFEM::TimeIntegrationOption, "SteadyState", "ImplicitTransient", diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp index 8bca1266398..109f8561778 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp @@ -20,27 +20,22 @@ #ifndef GEOS_PHYSICSSOLVERS_SIMPLEPDE_PHASEFIELDDAMAGEKERNELS_HPP_ #define GEOS_PHYSICSSOLVERS_SIMPLEPDE_PHASEFIELDDAMAGEKERNELS_HPP_ +#include "constitutive/solid/Damage.hpp" #include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" #include "finiteElement/elementFormulations/FiniteElementOperators.hpp" namespace geos { -enum class PhaseFieldDamageKernelLocalDissipation -{ - Linear, - Quadratic, -}; +/// Type of local dissipation function used in the phase-field damage model, shared with the +/// constitutive layer since it also affects degradation-function calibration there. +using PhaseFieldDamageKernelLocalDissipation = constitutive::LocalDissipationOption; //***************************************************************************** /** * @brief Implements kernels for solving the Damage(or phase-field) equation * in a phase-field fracture problem. * @copydoc geos::finiteElement::KernelBase - * @tparam NUM_NODES_PER_ELEM The number of nodes per element for the - * @p SUBREGION_TYPE. - * @tparam UNUSED An unused parameter since we are assuming that the test and - * trial space have the same number of support points. * * ### PhaseFieldDamageKernel Description * Implements the KernelBase interface functions required for solving the @@ -48,13 +43,7 @@ enum class PhaseFieldDamageKernelLocalDissipation * It uses the finite element kernel application functions such as * geos::finiteElement::RegionBasedKernelApplication. * - * In this implementation, the template parameter @p NUM_NODES_PER_ELEM is used - * in place of both @p NUM_TEST_SUPPORT_POINTS_PER_ELEM and - * @p NUM_TRIAL_SUPPORT_POINTS_PER_ELEM, which are assumed to be equal. This - * results in the @p UNUSED template parameter as only the NUM_NODES_PER_ELEM - * is passed to the ImplicitKernelBase template to form the base class. - * - * Additionally, the number of degrees of freedom per support point for both + * The number of degrees of freedom per support point for both * the test and trial spaces are specified as `1` when specifying the base * class. */ @@ -111,8 +100,7 @@ class PhaseFieldDamageKernel : CRSMatrixView< real64, globalIndex const > const inputMatrix, arrayView1d< real64 > const inputRhs, real64 const inputDt, - string const fieldName, - LocalDissipation localDissipationOption ): + string const fieldName ): Base( nodeManager, edgeManager, faceManager, @@ -127,9 +115,7 @@ class PhaseFieldDamageKernel : inputDt ), m_X( nodeManager.referencePosition()), m_nodalDamage( nodeManager.template getReference< array1d< real64 > >( fieldName )), - m_quadDamage( inputConstitutiveType.getNewDamage() ), - m_quadExtDrivingForce( inputConstitutiveType.getExtDrivingForce() ), - m_localDissipationOption( localDissipationOption ) + m_quadExtDrivingForce( inputConstitutiveType.getExtDrivingForce() ) {} //*************************************************************************** @@ -187,7 +173,7 @@ class PhaseFieldDamageKernel : } /** - * @copydoc geos::finiteElement::ImplicitKernelBase::quadraturePointJacobianContribution + * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel */ GEOS_HOST_DEVICE inline @@ -196,8 +182,8 @@ class PhaseFieldDamageKernel : StackVariables & stack ) const { - real64 const strainEnergyDensity = m_constitutiveUpdate.getStrainEnergyDensity( k, q ); - real64 const ell = m_constitutiveUpdate.getRegularizationLength(); + real64 const crackDrivingForce = m_constitutiveUpdate.getCrackDrivingForce( k, q ); + real64 const regularizationLength = m_constitutiveUpdate.getRegularizationLength(); real64 const Gc = m_constitutiveUpdate.getCriticalFractureEnergy( k ); real64 const threshold = m_constitutiveUpdate.getEnergyThreshold( k, q ); @@ -211,46 +197,54 @@ class PhaseFieldDamageKernel : real64 qp_grad_damage[3] = {0, 0, 0}; finiteElement::feOps::valueAndGradient( N, dNdX, stack.nodalDamageLocal, qp_damage, qp_grad_damage ); - real64 D = 0; //max between threshold and - // Elastic energy - if( m_localDissipationOption == LocalDissipation::Linear ) + LocalDissipation const localDissipationOption = m_constitutiveUpdate.m_localDissipationOption; + + // Crack driving force, floored by the threshold for Linear dissipation. + real64 const effectiveCrackDrivingForce = localDissipationOption == LocalDissipation::Linear ? + fmax( threshold, crackDrivingForce ) : + crackDrivingForce; + + // Coefficients that differ between the Linear and Quadratic dissipation models. + real64 localDissipation; // local dissipation contribution to the residual + real64 localDissipationDeriv; // its damage-derivative, for the Jacobian + real64 nonlocalDissipationGradientCoeff; // scaling of the damage-gradient term + real64 drivingForceCoeff; // scaling of the driving force terms + if( localDissipationOption == LocalDissipation::Linear ) + { + localDissipation = 3.0 / 16.0; + localDissipationDeriv = 0.0; + nonlocalDissipationGradientCoeff = 0.375; + drivingForceCoeff = 0.5; + } + else { - D = fmax( threshold, strainEnergyDensity ); + localDissipation = qp_damage; + localDissipationDeriv = 1.0; + nonlocalDissipationGradientCoeff = 1.0; + drivingForceCoeff = 1.0; } + real64 const degradationDeriv = m_constitutiveUpdate.getDegradationDerivative( k, qp_damage ); + real64 const degradationSecondDeriv = m_constitutiveUpdate.getDegradationSecondDerivative( k, qp_damage ); + real64 const scaledDrivingForce = drivingForceCoeff * regularizationLength * effectiveCrackDrivingForce / Gc; + // The external driving force only enters the Linear (nucleation) model. + real64 const scaledExtDrivingForce = localDissipationOption == LocalDissipation::Linear ? + 0.5 * regularizationLength * m_quadExtDrivingForce[k][q] / Gc : + 0.0; + for( localIndex a = 0; a < numNodesPerElem; ++a ) { - if( m_localDissipationOption == LocalDissipation::Linear ) - { - stack.localResidual[ a ] -= detJ * ( 3 * N[a] / 16 - + 0.375* ell * ell * LvArray::tensorOps::AiBi< 3 >( qp_grad_damage, dNdX[a] ) - + (0.5 * ell * D/Gc) * m_constitutiveUpdate.getDegradationDerivative( k, qp_damage ) * N[a] - + 0.5 * ell * m_quadExtDrivingForce[k][q]/Gc * N[a] ); - } - else - { - stack.localResidual[ a ] -= detJ * ( N[a] * qp_damage - + ( ell * ell * LvArray::tensorOps::AiBi< 3 >( qp_grad_damage, dNdX[a] ) - + N[a] * (ell*strainEnergyDensity/Gc) * m_constitutiveUpdate.getDegradationDerivative( k, qp_damage ) ) ); + stack.localResidual[ a ] -= detJ * ( localDissipation * N[a] + + nonlocalDissipationGradientCoeff * regularizationLength * regularizationLength * LvArray::tensorOps::AiBi< 3 >( qp_grad_damage, dNdX[a] ) + + scaledDrivingForce * degradationDeriv * N[a] + + scaledExtDrivingForce * N[a] ); - } for( localIndex b = 0; b < numNodesPerElem; ++b ) { - if( m_localDissipationOption == LocalDissipation::Linear ) - { - stack.localJacobian[ a ][ b ] -= detJ * ( 0.375* ell * ell * LvArray::tensorOps::AiBi< 3 >( dNdX[a], dNdX[b] ) - + (0.5 * ell * D/Gc) * m_constitutiveUpdate.getDegradationSecondDerivative( k, qp_damage ) * N[a] * N[b] ); - - } - else - { - stack.localJacobian[ a ][ b ] -= detJ * ( pow( ell, 2 ) * LvArray::tensorOps::AiBi< 3 >( dNdX[a], dNdX[b] ) - + N[a] * N[b] * (1 + m_constitutiveUpdate.getDegradationSecondDerivative( k, qp_damage ) * ell * strainEnergyDensity/Gc ) ); - } + stack.localJacobian[ a ][ b ] -= detJ * ( nonlocalDissipationGradientCoeff * regularizationLength * regularizationLength * LvArray::tensorOps::AiBi< 3 >( dNdX[a], dNdX[b] ) + + ( localDissipationDeriv + scaledDrivingForce * degradationSecondDeriv ) * N[a] * N[b] ); } } - - } /** @@ -293,14 +287,9 @@ class PhaseFieldDamageKernel : /// The global primary field array. arrayView1d< real64 const > const m_nodalDamage; - /// The array containing the damage on each quadrature point of all elements - arrayView2d< real64 const > const m_quadDamage; - /// The array containing the external driving force on each quadrature point of all elements arrayView2d< real64 const > const m_quadExtDrivingForce; - PhaseFieldDamageKernelLocalDissipation m_localDissipationOption; - }; using PhaseFieldDamageKernelFactory = finiteElement::KernelFactory< PhaseFieldDamageKernel, @@ -309,8 +298,7 @@ using PhaseFieldDamageKernelFactory = finiteElement::KernelFactory< PhaseFieldDa CRSMatrixView< real64, globalIndex const > const, arrayView1d< real64 > const, real64 const, - string const, - PhaseFieldDamageKernelLocalDissipation >; + string const >; } // namespace geos diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp index 0623667d56c..072eaf03963 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp @@ -53,22 +53,8 @@ class PhaseFieldPressurizedDamageKernel : CONSTITUTIVE_TYPE, FE_TYPE >; - using Base::numDofPerTestSupportPoint; - using Base::numDofPerTrialSupportPoint; - using Base::m_dofNumber; - using Base::m_dofRankOffset; - using Base::m_matrix; - using Base::m_rhs; using Base::m_elemsToNodes; using Base::m_constitutiveUpdate; - using Base::m_finiteElementSpace; - using Base::m_dt; - using Base::m_X; - using Base::m_nodalDamage; - using Base::m_quadDamage; - using Base::m_quadExtDrivingForce; - using Base::m_localDissipationOption; - using LocalDissipation = PhaseFieldDamageKernelLocalDissipation; /// Maximum number of nodes per element, which is equal to the maxNumTestSupportPointPerElem and /// maxNumTrialSupportPointPerElem by definition. When the FE_TYPE is not a Virtual Element, this @@ -77,7 +63,7 @@ class PhaseFieldPressurizedDamageKernel : /** * @brief Constructor - * @copydoc geosx::finiteElement::ImplicitKernelBase::ImplicitKernelBase + * @copydoc geos::finiteElement::ImplicitKernelBase::ImplicitKernelBase * @param fieldName The name of the primary field * (i.e. Temperature, Pressure, etc.) */ @@ -93,8 +79,7 @@ class PhaseFieldPressurizedDamageKernel : CRSMatrixView< real64, globalIndex const > const inputMatrix, arrayView1d< real64 > const inputRhs, real64 const inputDt, - string const fieldName, - LocalDissipation localDissipationOption ): + string const fieldName ): Base( nodeManager, edgeManager, faceManager, @@ -107,8 +92,7 @@ class PhaseFieldPressurizedDamageKernel : inputMatrix, inputRhs, inputDt, - fieldName, - localDissipationOption ), + fieldName ), m_disp( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ), m_fluidPressure( elementSubRegion.template getField< fields::flow::pressure >() ), m_fluidPressureGradient( elementSubRegion.template getReference< array2d< real64 > >( "pressureGradient" ) ) @@ -117,9 +101,9 @@ class PhaseFieldPressurizedDamageKernel : //*************************************************************************** /** * @class StackVariables - * @copydoc geosx::finiteElement::ImplicitKernelBase::StackVariables + * @copydoc geos::finiteElement::ImplicitKernelBase::StackVariables * - * Adds a stack array for the primary field. + * Adds a stack array for the element local nodal displacement. */ struct StackVariables : Base::StackVariables { @@ -140,7 +124,7 @@ class PhaseFieldPressurizedDamageKernel : /** * @brief Copy global values from primary field to a local stack array. - * @copydoc geosx::finiteElement::ImplicitKernelBase::setup + * @copydoc geos::finiteElement::ImplicitKernelBase::setup */ GEOS_HOST_DEVICE GEOS_FORCE_INLINE @@ -156,6 +140,12 @@ class PhaseFieldPressurizedDamageKernel : } } + /** + * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel + * + * Adds the fracture pressure contribution of Fei et al. (2023, IJNAMG) + * on top of the base damage kernel. + */ GEOS_HOST_DEVICE GEOS_FORCE_INLINE void quadraturePointKernel( localIndex const k, @@ -164,12 +154,12 @@ class PhaseFieldPressurizedDamageKernel : { Base::quadraturePointKernel( k, q, stack ); - real64 const ell = m_constitutiveUpdate.getRegularizationLength(); + real64 const regularizationLength = m_constitutiveUpdate.getRegularizationLength(); real64 const Gc = m_constitutiveUpdate.getCriticalFractureEnergy( k ); real64 const volStrain = m_constitutiveUpdate.getVolStrain( k, q ); real64 const biotCoeff = m_constitutiveUpdate.getBiotCoefficient( k ); - //Interpolate d and grad_d + //Interpolate d and u real64 N[ numNodesPerElem ]; real64 dNdX[ numNodesPerElem ][ 3 ]; real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, dNdX ); @@ -187,17 +177,21 @@ class PhaseFieldPressurizedDamageKernel : elemPresGradient[i] = m_fluidPressureGradient( k, i ); } + real64 const pressureDamageDeriv = m_constitutiveUpdate.pressureDamageFunctionDerivative( qp_damage ); + real64 const pressureDamageSecondDeriv = m_constitutiveUpdate.pressureDamageFunctionSecondDerivative( qp_damage ); + + // Pressure work density scaled by the dissipation normalization. + real64 const scaledPressureWork = 0.5 * regularizationLength/Gc * + ( ( 1.0 - biotCoeff ) * volStrain * m_fluidPressure( k ) + + LvArray::tensorOps::AiBi< 3 >( qp_disp, elemPresGradient ) ); + for( localIndex a = 0; a < numNodesPerElem; ++a ) { - /// Add pressure effects - stack.localResidual[ a ] -= detJ * 0.5 * ell/Gc * ( ( 1.0 - biotCoeff ) * volStrain * m_fluidPressure( k ) * m_constitutiveUpdate.pressureDamageFunctionDerivative( qp_damage ) * N[a] - + LvArray::tensorOps::AiBi< 3 >( qp_disp, elemPresGradient ) * m_constitutiveUpdate.pressureDamageFunctionDerivative( qp_damage ) * N[a] ); + stack.localResidual[ a ] -= detJ * scaledPressureWork * pressureDamageDeriv * N[a]; for( localIndex b = 0; b < numNodesPerElem; ++b ) { - stack.localJacobian[ a ][ b ] -= detJ * 0.5 * ell/Gc * - ( ( 1.0 - biotCoeff ) * volStrain * m_fluidPressure( k ) * m_constitutiveUpdate.pressureDamageFunctionSecondDerivative( qp_damage ) * N[a] * N[b] - + LvArray::tensorOps::AiBi< 3 >( qp_disp, elemPresGradient ) * m_constitutiveUpdate.pressureDamageFunctionSecondDerivative( qp_damage ) * N[a] * N[b] ); + stack.localJacobian[ a ][ b ] -= detJ * scaledPressureWork * pressureDamageSecondDeriv * N[a] * N[b]; } } } @@ -219,8 +213,7 @@ using PhaseFieldPressurizedDamageKernelFactory = finiteElement::KernelFactory< P CRSMatrixView< real64, globalIndex const > const, arrayView1d< real64 > const, real64 const, - string const, - PhaseFieldDamageKernelLocalDissipation >; + string const >; } // namespace geos diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index e78160e9ce0..7175090413d 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -4707,7 +4707,7 @@ When set to `all` output both convergence & iteration information to a csv.--> - + - +