From 0cb5390e0ebaf40e22dce463b6766f6a5d890e1b Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Wed, 15 Jul 2026 22:16:59 -0700 Subject: [PATCH 01/12] PhaseFieldFractureSolver cleanup --- .../multiphysics/PhaseFieldFractureSolver.cpp | 81 +++++++++---------- .../multiphysics/PhaseFieldFractureSolver.hpp | 29 ++++--- 2 files changed, 59 insertions(+), 51 deletions(-) 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 } } From baf40f8a1aea093c2455e685390fb6780ad01894 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Thu, 16 Jul 2026 13:49:22 -0700 Subject: [PATCH 02/12] PhaseFieldDamageFEM cleanup, fix typo, condense, rename --- .../simplePDE/PhaseFieldDamageFEM.cpp | 292 ++---------------- .../simplePDE/PhaseFieldDamageFEM.hpp | 34 +- .../simplePDE/PhaseFieldDamageFEMKernels.hpp | 96 +++--- .../PhaseFieldPressurizedDamageFEMKernels.hpp | 47 ++- src/coreComponents/schema/schema.xsd | 4 +- 5 files changed, 108 insertions(+), 365 deletions(-) diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp index bfe09b66e6d..feb00753893 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; @@ -89,10 +75,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 +105,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 +165,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(); - - 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 ); + arrayView1d< globalIndex const > const & dofIndex = + nodeManager.getReference< array1d< globalIndex > >( dofManager.getKey( m_fieldName ) ); - 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 +188,28 @@ 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, + m_localDissipationOption ) ); } - 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, + m_localDissipationOption ) ); } - - if( getLogLevel() >= 3 ) - { - NonlinearSolverParameters & solverParams = getNonlinearSolverParameters(); - integer newtonIter = solverParams.m_numNewtonIterations; -#endif } ); } @@ -481,24 +269,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 @@ -632,7 +402,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 +412,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..7ed3952d603 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 @@ -40,6 +35,9 @@ class FieldSpecification; class FiniteElementBase; class DomainPartition; +/// Forward declaration of the local dissipation enum defined in PhaseFieldDamageFEMKernels.hpp. +enum class PhaseFieldDamageKernelLocalDissipation : integer; + class PhaseFieldDamageFEM : public PhysicsSolverBase { public: @@ -141,11 +139,8 @@ class PhaseFieldDamageFEM : public PhysicsSolverBase ExplicitTransient }; - enum class LocalDissipation - { - Linear, - Quadratic, - }; + /// The type of local dissipation function, shared with the damage kernels. + using LocalDissipation = PhaseFieldDamageKernelLocalDissipation; struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct { @@ -160,27 +155,13 @@ 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; @@ -188,14 +169,9 @@ class PhaseFieldDamageFEM : public PhysicsSolverBase 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..b501c539a71 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp @@ -20,27 +20,30 @@ #ifndef GEOS_PHYSICSSOLVERS_SIMPLEPDE_PHASEFIELDDAMAGEKERNELS_HPP_ #define GEOS_PHYSICSSOLVERS_SIMPLEPDE_PHASEFIELDDAMAGEKERNELS_HPP_ +#include "common/format/EnumStrings.hpp" #include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" #include "finiteElement/elementFormulations/FiniteElementOperators.hpp" namespace geos { -enum class PhaseFieldDamageKernelLocalDissipation +/// Type of local dissipation function used in the phase-field damage model. +enum class PhaseFieldDamageKernelLocalDissipation : integer { Linear, Quadratic, }; +/// Declare strings associated with enumeration values. +ENUM_STRINGS( PhaseFieldDamageKernelLocalDissipation, + "Linear", + "Quadratic" ); + //***************************************************************************** /** * @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 +51,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. */ @@ -127,7 +124,6 @@ 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 ) {} @@ -187,7 +183,7 @@ class PhaseFieldDamageKernel : } /** - * @copydoc geos::finiteElement::ImplicitKernelBase::quadraturePointJacobianContribution + * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel */ GEOS_HOST_DEVICE inline @@ -197,7 +193,7 @@ class PhaseFieldDamageKernel : { real64 const strainEnergyDensity = m_constitutiveUpdate.getStrainEnergyDensity( k, q ); - real64 const ell = m_constitutiveUpdate.getRegularizationLength(); + real64 const regularizationLength = m_constitutiveUpdate.getRegularizationLength(); real64 const Gc = m_constitutiveUpdate.getCriticalFractureEnergy( k ); real64 const threshold = m_constitutiveUpdate.getEnergyThreshold( k, q ); @@ -211,46 +207,52 @@ 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 + // Elastic energy, floored by the threshold for Linear dissipation. + real64 const crackDrivingForce = m_localDissipationOption == LocalDissipation::Linear ? + fmax( threshold, strainEnergyDensity ) : + strainEnergyDensity; + + // 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( m_localDissipationOption == LocalDissipation::Linear ) { - D = fmax( threshold, strainEnergyDensity ); + localDissipation = 3.0 / 16.0; + localDissipationDeriv = 0.0; + nonlocalDissipationGradientCoeff = 0.375; + drivingForceCoeff = 0.5; } + else + { + 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 * crackDrivingForce / Gc; + // The external driving force only enters the Linear (nucleation) model. + real64 const scaledExtDrivingForce = m_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,13 +295,11 @@ 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; + /// The type of local dissipation function (Linear or Quadratic). + PhaseFieldDamageKernelLocalDissipation const m_localDissipationOption; }; diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp index 0623667d56c..548d0e732ea 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp @@ -53,21 +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 @@ -77,7 +64,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.) */ @@ -117,9 +104,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 +127,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 +143,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 +157,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 +180,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]; } } } 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.--> - + - + From 24f37b0ca1c7f6293213ff2bce07e8ab6e509c16 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Thu, 16 Jul 2026 15:00:05 -0700 Subject: [PATCH 03/12] added a check for nucleation model to see if the dissipaiton option is linear --- .../constitutive/solid/Damage.hpp | 2 ++ .../simplePDE/PhaseFieldDamageFEM.cpp | 28 +++++++++++++++++++ .../simplePDE/PhaseFieldDamageFEM.hpp | 2 ++ 3 files changed, 32 insertions(+) diff --git a/src/coreComponents/constitutive/solid/Damage.hpp b/src/coreComponents/constitutive/solid/Damage.hpp index b54bb30f372..d720cb3741b 100644 --- a/src/coreComponents/constitutive/solid/Damage.hpp +++ b/src/coreComponents/constitutive/solid/Damage.hpp @@ -437,6 +437,8 @@ class Damage : public BASE arrayView2d< real64 const > getExtDrivingForce() const { return m_extDrivingForce; } + integer getExtDrivingForceFlag() const { return m_extDrivingForceFlag; } + KernelWrapper createKernelUpdates() const { diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp index feb00753893..99fdef89afe 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp @@ -23,6 +23,7 @@ #include "common/DataTypes.hpp" #include "common/TimingMacros.hpp" +#include "constitutive/ConstitutivePassThru.hpp" #include "constitutive/solid/Damage.hpp" #include "constitutive/solid/SolidBase.hpp" #include "fieldSpecification/FieldSpecificationImpl.hpp" @@ -105,6 +106,33 @@ void PhaseFieldDamageFEM::registerDataOnMesh( Group & meshBodies ) } ); } +void PhaseFieldDamageFEM::initializePreSubGroups() +{ + PhysicsSolverBase::initializePreSubGroups(); + + DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & subRegion ) + { + string const & solidModelName = subRegion.getReference< string >( viewKeyStruct::solidModelNamesString() ); + SolidBase & solidModel = subRegion.getConstitutiveModel< SolidBase >( solidModelName ); + + ConstitutivePassThru< DamageBase >::execute( solidModel, [&]( auto & damageModel ) + { + GEOS_ERROR_IF( damageModel.getExtDrivingForceFlag() && m_localDissipationOption != LocalDissipation::Linear, + GEOS_FMT( "{}: the external driving force (extDrivingForceFlag) is only supported " + "with the Linear local dissipation function.", getName() ) ); + } ); + } ); + } ); +} + real64 PhaseFieldDamageFEM::solverStep( real64 const & time_n, real64 const & dt, const int cycleNumber, diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp index 7ed3952d603..aed89ae6806 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp @@ -58,6 +58,8 @@ class PhaseFieldDamageFEM : public PhysicsSolverBase virtual void registerDataOnMesh( Group & meshBodies ) override final; + virtual void initializePreSubGroups() override; + /** * @defgroup Solver Interface Functions * From 9b33c2c3f555b17fb4e4180a6256d6aa3427b997 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Fri, 17 Jul 2026 19:56:39 -0700 Subject: [PATCH 04/12] move localDissipationoption to constitutive model and added FractureModelType option --- .../constitutive/solid/Damage.cpp | 33 ++-- .../constitutive/solid/Damage.hpp | 161 +++++++++++++----- .../constitutive/solid/DamageSpectral.cpp | 16 +- .../constitutive/solid/DamageSpectral.hpp | 78 ++------- .../constitutive/solid/DamageVolDev.cpp | 10 ++ .../constitutive/solid/DamageVolDev.hpp | 19 ++- .../simplePDE/PhaseFieldDamageFEM.cpp | 38 +---- .../simplePDE/PhaseFieldDamageFEM.hpp | 10 -- .../simplePDE/PhaseFieldDamageFEMKernels.hpp | 36 ++-- .../PhaseFieldPressurizedDamageFEMKernels.hpp | 10 +- 10 files changed, 200 insertions(+), 211 deletions(-) diff --git a/src/coreComponents/constitutive/solid/Damage.cpp b/src/coreComponents/constitutive/solid/Damage.cpp index 525d6e0f0bf..d051cdb9e3e 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 ). @@ -99,21 +105,24 @@ 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 diff --git a/src/coreComponents/constitutive/solid/Damage.hpp b/src/coreComponents/constitutive/solid/Damage.hpp index d720cb3741b..407db05c7b9 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 @@ -80,7 +103,8 @@ class DamageUpdates : public UPDATE_BASE 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, @@ -97,7 +121,8 @@ class DamageUpdates : public UPDATE_BASE 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 +139,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 +173,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 +193,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 +299,7 @@ class DamageUpdates : public UPDATE_BASE m_volStrain( k, q ) = traceOfStrain; - if( m_extDrivingForceFlag ) + if( m_fractureModelType == FractureModelType::Nucleation ) { real64 stressP; real64 stressQ; @@ -335,17 +385,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 @@ -392,8 +450,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,7 +498,9 @@ class Damage : public BASE arrayView2d< real64 const > getExtDrivingForce() const { return m_extDrivingForce; } - integer getExtDrivingForceFlag() const { return m_extDrivingForceFlag; } + FractureModelType getFractureModelType() const { return m_fractureModelType; } + + LocalDissipationOption getLocalDissipationOption() const { return m_localDissipationOption; } KernelWrapper createKernelUpdates() const @@ -452,7 +515,8 @@ class Damage : public BASE m_criticalFractureEnergy.toView(), m_criticalStrainEnergy, m_degradationLowerLimit, - m_extDrivingForceFlag, + m_fractureModelType, + m_localDissipationOption, m_tensileStrength.toView(), m_compressiveStrength.toView(), m_deltaCoefficient.toView(), @@ -469,8 +533,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 @@ -512,8 +578,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..2b29f64fbae 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 @@ -49,14 +47,16 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > 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, + inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputFractureModelType, + inputLocalDissipationOption, inputTensileStrength, inputCompressiveStrength, inputDeltaCoefficient, inputBiotCoefficient, std::forward< PARAMS >( baseParams )... ) {} @@ -65,11 +65,7 @@ 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_volStrain; @@ -80,7 +76,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 +88,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, @@ -276,17 +228,6 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > 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; - } - }; @@ -307,7 +248,8 @@ class DamageSpectral : public Damage< BASE > 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,6 +261,7 @@ 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 { @@ -332,7 +275,8 @@ class DamageSpectral : public Damage< BASE > 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..ed62c22258d 100644 --- a/src/coreComponents/constitutive/solid/DamageVolDev.hpp +++ b/src/coreComponents/constitutive/solid/DamageVolDev.hpp @@ -45,14 +45,16 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > 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, + inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputFractureModelType, + inputLocalDissipationOption, inputTensileStrength, inputCompressiveStrength, inputDeltaCoefficient, inputBiotCoefficient, std::forward< PARAMS >( baseParams )... ) {} @@ -63,9 +65,6 @@ 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_volStrain; @@ -76,7 +75,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; @@ -181,7 +181,8 @@ class DamageVolDev : public Damage< BASE > 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,6 +194,7 @@ 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 { @@ -206,7 +208,8 @@ class DamageVolDev : public Damage< BASE > 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/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp index 99fdef89afe..1f735fc070b 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp @@ -23,7 +23,6 @@ #include "common/DataTypes.hpp" #include "common/TimingMacros.hpp" -#include "constitutive/ConstitutivePassThru.hpp" #include "constitutive/solid/Damage.hpp" #include "constitutive/solid/SolidBase.hpp" #include "fieldSpecification/FieldSpecificationImpl.hpp" @@ -54,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 ). @@ -106,33 +101,6 @@ void PhaseFieldDamageFEM::registerDataOnMesh( Group & meshBodies ) } ); } -void PhaseFieldDamageFEM::initializePreSubGroups() -{ - PhysicsSolverBase::initializePreSubGroups(); - - DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { - ElementRegionManager & elemManager = mesh.getElemManager(); - - elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & subRegion ) - { - string const & solidModelName = subRegion.getReference< string >( viewKeyStruct::solidModelNamesString() ); - SolidBase & solidModel = subRegion.getConstitutiveModel< SolidBase >( solidModelName ); - - ConstitutivePassThru< DamageBase >::execute( solidModel, [&]( auto & damageModel ) - { - GEOS_ERROR_IF( damageModel.getExtDrivingForceFlag() && m_localDissipationOption != LocalDissipation::Linear, - GEOS_FMT( "{}: the external driving force (extDrivingForceFlag) is only supported " - "with the Linear local dissipation function.", getName() ) ); - } ); - } ); - } ); -} - real64 PhaseFieldDamageFEM::solverStep( real64 const & time_n, real64 const & dt, const int cycleNumber, @@ -225,8 +193,7 @@ void PhaseFieldDamageFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM( time_n localMatrix, localRhs, dt, - m_fieldName, - m_localDissipationOption ) ); + m_fieldName ) ); } else { @@ -235,8 +202,7 @@ void PhaseFieldDamageFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM( time_n localMatrix, localRhs, dt, - m_fieldName, - m_localDissipationOption ) ); + m_fieldName ) ); } } ); } diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp index aed89ae6806..3208e6f9501 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.hpp @@ -35,9 +35,6 @@ class FieldSpecification; class FiniteElementBase; class DomainPartition; -/// Forward declaration of the local dissipation enum defined in PhaseFieldDamageFEMKernels.hpp. -enum class PhaseFieldDamageKernelLocalDissipation : integer; - class PhaseFieldDamageFEM : public PhysicsSolverBase { public: @@ -58,8 +55,6 @@ class PhaseFieldDamageFEM : public PhysicsSolverBase virtual void registerDataOnMesh( Group & meshBodies ) override final; - virtual void initializePreSubGroups() override; - /** * @defgroup Solver Interface Functions * @@ -141,13 +136,9 @@ class PhaseFieldDamageFEM : public PhysicsSolverBase ExplicitTransient }; - /// The type of local dissipation function, shared with the damage kernels. - using LocalDissipation = PhaseFieldDamageKernelLocalDissipation; - 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"; } @@ -165,7 +156,6 @@ class PhaseFieldDamageFEM : public PhysicsSolverBase private: string m_fieldName; TimeIntegrationOption m_timeIntegrationOption; - LocalDissipation m_localDissipationOption; integer m_irreversibilityFlag; real64 m_damageUpperBound; integer m_fracturePressureTermFlag; diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp index b501c539a71..2189dbe2f05 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp @@ -20,24 +20,16 @@ #ifndef GEOS_PHYSICSSOLVERS_SIMPLEPDE_PHASEFIELDDAMAGEKERNELS_HPP_ #define GEOS_PHYSICSSOLVERS_SIMPLEPDE_PHASEFIELDDAMAGEKERNELS_HPP_ -#include "common/format/EnumStrings.hpp" +#include "constitutive/solid/Damage.hpp" #include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" #include "finiteElement/elementFormulations/FiniteElementOperators.hpp" namespace geos { -/// Type of local dissipation function used in the phase-field damage model. -enum class PhaseFieldDamageKernelLocalDissipation : integer -{ - Linear, - Quadratic, -}; - -/// Declare strings associated with enumeration values. -ENUM_STRINGS( 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; //***************************************************************************** /** @@ -108,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, @@ -124,8 +115,7 @@ class PhaseFieldDamageKernel : inputDt ), m_X( nodeManager.referencePosition()), m_nodalDamage( nodeManager.template getReference< array1d< real64 > >( fieldName )), - m_quadExtDrivingForce( inputConstitutiveType.getExtDrivingForce() ), - m_localDissipationOption( localDissipationOption ) + m_quadExtDrivingForce( inputConstitutiveType.getExtDrivingForce() ) {} //*************************************************************************** @@ -207,8 +197,10 @@ class PhaseFieldDamageKernel : real64 qp_grad_damage[3] = {0, 0, 0}; finiteElement::feOps::valueAndGradient( N, dNdX, stack.nodalDamageLocal, qp_damage, qp_grad_damage ); + LocalDissipation const localDissipationOption = m_constitutiveUpdate.m_localDissipationOption; + // Elastic energy, floored by the threshold for Linear dissipation. - real64 const crackDrivingForce = m_localDissipationOption == LocalDissipation::Linear ? + real64 const crackDrivingForce = localDissipationOption == LocalDissipation::Linear ? fmax( threshold, strainEnergyDensity ) : strainEnergyDensity; @@ -217,7 +209,7 @@ class PhaseFieldDamageKernel : 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( m_localDissipationOption == LocalDissipation::Linear ) + if( localDissipationOption == LocalDissipation::Linear ) { localDissipation = 3.0 / 16.0; localDissipationDeriv = 0.0; @@ -236,7 +228,7 @@ class PhaseFieldDamageKernel : real64 const degradationSecondDeriv = m_constitutiveUpdate.getDegradationSecondDerivative( k, qp_damage ); real64 const scaledDrivingForce = drivingForceCoeff * regularizationLength * crackDrivingForce / Gc; // The external driving force only enters the Linear (nucleation) model. - real64 const scaledExtDrivingForce = m_localDissipationOption == LocalDissipation::Linear ? + real64 const scaledExtDrivingForce = localDissipationOption == LocalDissipation::Linear ? 0.5 * regularizationLength * m_quadExtDrivingForce[k][q] / Gc : 0.0; @@ -298,9 +290,6 @@ class PhaseFieldDamageKernel : /// The array containing the external driving force on each quadrature point of all elements arrayView2d< real64 const > const m_quadExtDrivingForce; - /// The type of local dissipation function (Linear or Quadratic). - PhaseFieldDamageKernelLocalDissipation const 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 548d0e732ea..072eaf03963 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldPressurizedDamageFEMKernels.hpp @@ -55,7 +55,6 @@ class PhaseFieldPressurizedDamageKernel : using Base::m_elemsToNodes; using Base::m_constitutiveUpdate; - 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 @@ -80,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, @@ -94,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" ) ) @@ -216,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 From f2818e795430e272775af0f4b95c681db2621bdc Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Fri, 17 Jul 2026 19:57:28 -0700 Subject: [PATCH 05/12] updated xmls --- inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml | 1 - inputFiles/phaseField/PhaseFieldFracture_DamageAndLoad.xml | 5 ++--- inputFiles/phaseField/PhaseFieldFracture_DamageOnly.xml | 2 -- inputFiles/phaseField/PhaseFieldFracture_Nucleation_base.xml | 2 -- .../phaseField/PhaseFieldFracture_Nucleation_benchmark.xml | 2 +- .../phaseField/PhaseFieldFracture_Nucleation_calibration.xml | 2 +- .../phaseField/PhaseFieldFracture_Nucleation_smoke.xml | 4 +--- inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml | 2 -- inputFiles/phaseField/PhaseFieldFracture_VolDevSplit.xml | 5 ++--- .../PhaseFieldPoromechanics_Nucleation_Injection.xml | 1 - .../PhaseFieldPoromechanics_Nucleation_Wellbore.xml | 1 - .../phaseField/PhaseFieldPoromechanics_Nucleation_base.xml | 2 +- .../phaseField/PhaseFieldPoromechanics_Nucleation_smoke.xml | 1 - .../phaseField/PhaseFieldPoromechanics_multiCracks.xml | 1 - 14 files changed, 8 insertions(+), 23 deletions(-) diff --git a/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml b/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml index 4a6716a4385..2888ddfe3a9 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"/> diff --git a/inputFiles/phaseField/PhaseFieldFracture_DamageOnly.xml b/inputFiles/phaseField/PhaseFieldFracture_DamageOnly.xml index b22c9a2bec1..54c7d63bbed 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_DamageOnly.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_DamageOnly.xml @@ -35,8 +35,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - - localDissipation="Linear" logLevel="1" targetRegions="{ Region1 }"> diff --git a/inputFiles/phaseField/PhaseFieldFracture_Nucleation_calibration.xml b/inputFiles/phaseField/PhaseFieldFracture_Nucleation_calibration.xml index 43c39532410..0d4dce08eba 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_Nucleation_calibration.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_Nucleation_calibration.xml @@ -58,7 +58,7 @@ lengthScale="0.4" defaultCriticalFractureEnergy="4e-3" criticalStrainEnergy="0" - extDrivingForceFlag="1" + fractureModelType="Nucleation" defaultTensileStrength="5.0" defaultCompressiveStrength="40.0" defaultDeltaCoefficient="0.0"/> diff --git a/inputFiles/phaseField/PhaseFieldFracture_Nucleation_smoke.xml b/inputFiles/phaseField/PhaseFieldFracture_Nucleation_smoke.xml index 083819e314a..42027245959 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_Nucleation_smoke.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_Nucleation_smoke.xml @@ -35,8 +35,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - - localDissipation="Linear" irreversibilityFlag="1" damageUpperBound="1.0" logLevel="1" @@ -111,7 +109,7 @@ lengthScale="0.001" defaultCriticalFractureEnergy="91" criticalStrainEnergy="0" - extDrivingForceFlag="1" + fractureModelType="Nucleation" defaultTensileStrength="27e6" defaultCompressiveStrength="77e6" defaultDeltaCoefficient="0.0"/> diff --git a/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml b/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml index e0efbc3ba00..4c65946802a 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml @@ -33,8 +33,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - - localDissipation="Linear" logLevel="1" targetRegions="{ Region1 }"> + criticalStrainEnergy="14.88" + localDissipationOption="Quadratic"/> diff --git a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Injection.xml b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Injection.xml index 97d75c706b4..4cf29770ca8 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" diff --git a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Wellbore.xml b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Wellbore.xml index 1f3b5d12e0d..c08cec8eb5c 100644 --- a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Wellbore.xml +++ b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_Wellbore.xml @@ -40,7 +40,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - localDissipation="Linear" irreversibilityFlag="1" damageUpperBound="0.95" fracturePressureTermFlag="1" diff --git a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_base.xml b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_base.xml index 491ca55d9db..9486fae7b2b 100644 --- a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_base.xml +++ b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_base.xml @@ -29,7 +29,7 @@ lengthScale="0.5" defaultCriticalFractureEnergy="4e-3" criticalStrainEnergy="0" - extDrivingForceFlag="1" + fractureModelType="Nucleation" degradationLowerLimit="1e-6" defaultTensileStrength="5.5" defaultCompressiveStrength="40.0" diff --git a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_smoke.xml b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_smoke.xml index dff245983f4..83c206f642b 100644 --- a/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_smoke.xml +++ b/inputFiles/phaseField/PhaseFieldPoromechanics_Nucleation_smoke.xml @@ -39,7 +39,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - localDissipation="Linear" irreversibilityFlag="1" damageUpperBound="1.0" fracturePressureTermFlag="1" diff --git a/inputFiles/phaseField/PhaseFieldPoromechanics_multiCracks.xml b/inputFiles/phaseField/PhaseFieldPoromechanics_multiCracks.xml index 08f64cd5dd7..f557ff6e48d 100644 --- a/inputFiles/phaseField/PhaseFieldPoromechanics_multiCracks.xml +++ b/inputFiles/phaseField/PhaseFieldPoromechanics_multiCracks.xml @@ -36,7 +36,6 @@ discretization="FE1" timeIntegrationOption="SteadyState" fieldName="Damage" - localDissipation="Linear" irreversibilityFlag="1" damageUpperBound="1.0" fracturePressureTermFlag="1" From 7dad7b3b74c6c16b963bbdf3694dbab0b89e9227 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sat, 18 Jul 2026 17:41:03 -0700 Subject: [PATCH 06/12] renamed sed to crackDrivingForce --- .../constitutive/solid/Damage.cpp | 4 ++-- .../constitutive/solid/Damage.hpp | 22 +++++++++---------- .../constitutive/solid/DamageSpectral.hpp | 20 ++++++++--------- .../constitutive/solid/DamageVolDev.hpp | 20 ++++++++--------- .../constitutive/solid/SolidFields.hpp | 6 ++--- .../simplePDE/PhaseFieldDamageFEMKernels.hpp | 12 +++++----- 6 files changed, 42 insertions(+), 42 deletions(-) diff --git a/src/coreComponents/constitutive/solid/Damage.cpp b/src/coreComponents/constitutive/solid/Damage.cpp index d051cdb9e3e..314756e8b88 100644 --- a/src/coreComponents/constitutive/solid/Damage.cpp +++ b/src/coreComponents/constitutive/solid/Damage.cpp @@ -82,7 +82,7 @@ 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::volStrain >( &m_volStrain ); @@ -145,7 +145,7 @@ 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_volStrain.resize( 0, numPts ); m_extDrivingForce.resize( 0, numPts ); m_biotCoefficient.resize( parent.size() ); diff --git a/src/coreComponents/constitutive/solid/Damage.hpp b/src/coreComponents/constitutive/solid/Damage.hpp index 407db05c7b9..92f7a8e3a74 100644 --- a/src/coreComponents/constitutive/solid/Damage.hpp +++ b/src/coreComponents/constitutive/solid/Damage.hpp @@ -96,7 +96,7 @@ 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 & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, @@ -114,7 +114,7 @@ class DamageUpdates : public UPDATE_BASE m_newDamage( inputNewDamage ), m_oldDamage( inputOldDamage ), m_damageGrad( inputDamageGrad ), - m_strainEnergyDensity( inputStrainEnergyDensity ), + m_crackDrivingForce( inputCrackDrivingForce ), m_volStrain( inputVolumetricStrain ), m_extDrivingForce ( inputExtDrivingForce ), m_lengthScale( inputLengthScale ), @@ -344,21 +344,21 @@ class DamageUpdates : public UPDATE_BASE } - // TODO: The code below assumes the strain energy density will never be + // TODO: The code below assumes the crack driving force 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 ) ) + if( sed > m_crackDrivingForce( k, q ) ) { - m_strainEnergyDensity( k, q ) = sed; + m_crackDrivingForce( k, q ) = sed; } - return m_strainEnergyDensity( k, q ); + return m_crackDrivingForce( k, q ); } GEOS_HOST_DEVICE @@ -430,7 +430,7 @@ 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 volumetric strain at the quadrature point arrayView2d< real64 > const m_volStrain; @@ -508,7 +508,7 @@ class Damage : public BASE return BASE::template createDerivedKernelUpdates< KernelWrapper >( m_newDamage.toView(), m_oldDamage.toView(), m_damageGrad.toView(), - m_strainEnergyDensity.toView(), + m_crackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, @@ -558,7 +558,7 @@ 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 volumetric strain at the quadrature point array2d< real64 > m_volStrain; diff --git a/src/coreComponents/constitutive/solid/DamageSpectral.hpp b/src/coreComponents/constitutive/solid/DamageSpectral.hpp index 2b29f64fbae..13be1893cab 100644 --- a/src/coreComponents/constitutive/solid/DamageSpectral.hpp +++ b/src/coreComponents/constitutive/solid/DamageSpectral.hpp @@ -40,7 +40,7 @@ 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 & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, @@ -54,7 +54,7 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > arrayView1d< real64 > const & inputDeltaCoefficient, arrayView1d< real64 > const & inputBiotCoefficient, PARAMS && ... baseParams ): - DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputStrainEnergyDensity, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, + DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputCrackDrivingForce, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputFractureModelType, inputLocalDissipationOption, inputTensileStrength, inputCompressiveStrength, inputDeltaCoefficient, inputBiotCoefficient, @@ -67,7 +67,7 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > using DamageUpdates< UPDATE_BASE >::saveConvergedState; using DamageUpdates< UPDATE_BASE >::getDegradationValue; - using DamageUpdates< UPDATE_BASE >::m_strainEnergyDensity; + using DamageUpdates< UPDATE_BASE >::m_crackDrivingForce; using DamageUpdates< UPDATE_BASE >::m_volStrain; using DamageUpdates< UPDATE_BASE >::m_criticalStrainEnergy; using DamageUpdates< UPDATE_BASE >::m_extDrivingForce; @@ -202,9 +202,9 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > real64 const sed = 0.5 * lambda * tracePlus * tracePlus + mu * doubleContraction( positivePartOfStrain, positivePartOfStrain ); - if( sed > m_strainEnergyDensity( k, q ) ) + if( sed > m_crackDrivingForce( k, q ) ) { - m_strainEnergyDensity( k, q ) = sed; + m_crackDrivingForce( k, q ) = sed; } } @@ -222,10 +222,10 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > GEOS_HOST_DEVICE - virtual real64 getStrainEnergyDensity( localIndex const k, - localIndex const q ) const override final + virtual real64 getCrackDrivingForce( localIndex const k, + localIndex const q ) const override final { - return m_strainEnergyDensity( k, q ); + return m_crackDrivingForce( k, q ); } }; @@ -241,7 +241,7 @@ 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_volStrain; using Damage< BASE >::m_extDrivingForce; using Damage< BASE >::m_criticalFractureEnergy; @@ -268,7 +268,7 @@ class DamageSpectral : public Damage< BASE > return BASE::template createDerivedKernelUpdates< KernelWrapper >( m_newDamage.toView(), m_oldDamage.toView(), m_damageGrad.toView(), - m_strainEnergyDensity.toView(), + m_crackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, diff --git a/src/coreComponents/constitutive/solid/DamageVolDev.hpp b/src/coreComponents/constitutive/solid/DamageVolDev.hpp index ed62c22258d..9a0486722d9 100644 --- a/src/coreComponents/constitutive/solid/DamageVolDev.hpp +++ b/src/coreComponents/constitutive/solid/DamageVolDev.hpp @@ -38,7 +38,7 @@ 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 & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, @@ -52,7 +52,7 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > arrayView1d< real64 > const & inputDeltaCoefficient, arrayView1d< real64 > const & inputBiotCoefficient, PARAMS && ... baseParams ): - DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputStrainEnergyDensity, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, + DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputCrackDrivingForce, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputFractureModelType, inputLocalDissipationOption, inputTensileStrength, inputCompressiveStrength, inputDeltaCoefficient, inputBiotCoefficient, @@ -66,7 +66,7 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > using DamageUpdates< UPDATE_BASE >::getDegradationValue; - using DamageUpdates< UPDATE_BASE >::m_strainEnergyDensity; + using DamageUpdates< UPDATE_BASE >::m_crackDrivingForce; using DamageUpdates< UPDATE_BASE >::m_volStrain; using DamageUpdates< UPDATE_BASE >::m_criticalStrainEnergy; using DamageUpdates< UPDATE_BASE >::m_extDrivingForce; @@ -147,18 +147,18 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > sed += 0.5 * (stressP * volStrain) / factor; } - if( sed > m_strainEnergyDensity( k, q ) ) + if( sed > m_crackDrivingForce( k, q ) ) { - m_strainEnergyDensity( k, q ) = sed; + m_crackDrivingForce( k, q ) = sed; } } GEOS_HOST_DEVICE - virtual real64 getStrainEnergyDensity( localIndex const k, - localIndex const q ) const override final + virtual real64 getCrackDrivingForce( localIndex const k, + localIndex const q ) const override final { - return m_strainEnergyDensity( k, q ); + return m_crackDrivingForce( k, q ); } }; @@ -174,7 +174,7 @@ 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_volStrain; using Damage< BASE >::m_extDrivingForce; using Damage< BASE >::m_criticalFractureEnergy; @@ -201,7 +201,7 @@ class DamageVolDev : public Damage< BASE > return BASE::template createDerivedKernelUpdates< KernelWrapper >( m_newDamage.toView(), m_oldDamage.toView(), m_damageGrad.toView(), - m_strainEnergyDensity.toView(), + m_crackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, diff --git a/src/coreComponents/constitutive/solid/SolidFields.hpp b/src/coreComponents/constitutive/solid/SolidFields.hpp index 8ed466eadc7..3593e857117 100644 --- a/src/coreComponents/constitutive/solid/SolidFields.hpp +++ b/src/coreComponents/constitutive/solid/SolidFields.hpp @@ -154,13 +154,13 @@ 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( volStrain, "volStrain", diff --git a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp index 2189dbe2f05..109f8561778 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEMKernels.hpp @@ -182,7 +182,7 @@ class PhaseFieldDamageKernel : StackVariables & stack ) const { - real64 const strainEnergyDensity = m_constitutiveUpdate.getStrainEnergyDensity( k, q ); + 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 ); @@ -199,10 +199,10 @@ class PhaseFieldDamageKernel : LocalDissipation const localDissipationOption = m_constitutiveUpdate.m_localDissipationOption; - // Elastic energy, floored by the threshold for Linear dissipation. - real64 const crackDrivingForce = localDissipationOption == LocalDissipation::Linear ? - fmax( threshold, strainEnergyDensity ) : - strainEnergyDensity; + // 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 @@ -226,7 +226,7 @@ class PhaseFieldDamageKernel : real64 const degradationDeriv = m_constitutiveUpdate.getDegradationDerivative( k, qp_damage ); real64 const degradationSecondDeriv = m_constitutiveUpdate.getDegradationSecondDerivative( k, qp_damage ); - real64 const scaledDrivingForce = drivingForceCoeff * regularizationLength * crackDrivingForce / Gc; + 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 : From 0188a026712e645eefdbc1b82f29445c7674fdd5 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sat, 18 Jul 2026 18:46:28 -0700 Subject: [PATCH 07/12] cleaned up PhaseFieldPoromechanicsSolver --- .../multiphysics/PhaseFieldPoromechanicsSolver.cpp | 12 +----------- .../multiphysics/PhaseFieldPoromechanicsSolver.hpp | 9 +++------ 2 files changed, 4 insertions(+), 17 deletions(-) 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_ */ From f492395ec533a7973adfc6923ea9eca06dd1dc02 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sat, 18 Jul 2026 21:21:19 -0700 Subject: [PATCH 08/12] apply Damage default field values after array allocation --- .../PhaseFieldFracture_CohesiveModel.xml | 8 ----- .../PhaseFieldFracture_DamageAndLoad.xml | 8 ----- .../PhaseFieldFracture_DamageOnly.xml | 8 ----- ...haseFieldFracture_Nucleation_benchmark.xml | 32 ------------------- ...seFieldFracture_Nucleation_calibration.xml | 32 ------------------- .../PhaseFieldFracture_Nucleation_smoke.xml | 32 ------------------- .../PhaseFieldFracture_SpectralSplit.xml | 7 ---- .../PhaseFieldFracture_VolDevSplit.xml | 7 ---- ...ieldPoromechanics_Nucleation_Injection.xml | 32 ------------------- ...FieldPoromechanics_Nucleation_Wellbore.xml | 31 ------------------ ...aseFieldPoromechanics_Nucleation_smoke.xml | 32 ------------------- .../PhaseFieldPoromechanics_multiCracks.xml | 7 ---- .../constitutive/solid/Damage.cpp | 27 ++++++++-------- 13 files changed, 14 insertions(+), 249 deletions(-) diff --git a/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml b/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml index 2888ddfe3a9..5d36f307c94 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_CohesiveModel.xml @@ -118,14 +118,6 @@ objectPath="nodeManager" scale="0" setNames="{xneg, xpos, yneg, ypos, zneg, zpos}"/> --> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ::postInputInitialization() "delta coefficient must be input and non-negative when the" " 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 > @@ -154,6 +141,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 ); } From 50c014ccf8ce3145cc10671e45b2a5a3992c6450 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Wed, 22 Jul 2026 11:55:53 -0700 Subject: [PATCH 09/12] refactored and added spectral/vol-dev split test decks for all fracture-model types --- inputFiles/phaseField/PhaseFieldFracture.ats | 20 +- ...FieldFracture_DecompositionSplit_base.xml} | 68 +----- .../PhaseFieldFracture_SpectralSplit.xml | 198 ------------------ ...aseFieldFracture_SpectralSplit_Brittle.xml | 61 ++++++ ...seFieldFracture_SpectralSplit_Cohesive.xml | 59 ++++++ ...PhaseFieldFracture_VolDevSplit_Brittle.xml | 60 ++++++ ...haseFieldFracture_VolDevSplit_Cohesive.xml | 60 ++++++ 7 files changed, 261 insertions(+), 265 deletions(-) rename inputFiles/phaseField/{PhaseFieldFracture_VolDevSplit.xml => PhaseFieldFracture_DecompositionSplit_base.xml} (60%) delete mode 100644 inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml create mode 100644 inputFiles/phaseField/PhaseFieldFracture_SpectralSplit_Brittle.xml create mode 100644 inputFiles/phaseField/PhaseFieldFracture_SpectralSplit_Cohesive.xml create mode 100644 inputFiles/phaseField/PhaseFieldFracture_VolDevSplit_Brittle.xml create mode 100644 inputFiles/phaseField/PhaseFieldFracture_VolDevSplit_Cohesive.xml diff --git a/inputFiles/phaseField/PhaseFieldFracture.ats b/inputFiles/phaseField/PhaseFieldFracture.ats index 88b731578d3..f22f302d110 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,29 @@ 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_Nucleation_smoke", description='Testing the nucleation model for PF Fracture', partitions=((1, 1, 1), (2, 2, 1)), diff --git a/inputFiles/phaseField/PhaseFieldFracture_VolDevSplit.xml b/inputFiles/phaseField/PhaseFieldFracture_DecompositionSplit_base.xml similarity index 60% rename from inputFiles/phaseField/PhaseFieldFracture_VolDevSplit.xml rename to inputFiles/phaseField/PhaseFieldFracture_DecompositionSplit_base.xml index 6b74828c8d2..ef65044e078 100644 --- a/inputFiles/phaseField/PhaseFieldFracture_VolDevSplit.xml +++ b/inputFiles/phaseField/PhaseFieldFracture_DecompositionSplit_base.xml @@ -7,11 +7,11 @@ name="PhaseFieldSolve" solidSolverName="SolidMechSolve" damageSolverName="DamageSolve" - logLevel="2" + logLevel="1" targetRegions="{ Region1 }"> @@ -71,37 +71,6 @@ xMax="{ 1.01, 1.01, 0.11 }"/> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + \ No newline at end of file diff --git a/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml b/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml deleted file mode 100644 index a5c31a244b5..00000000000 --- a/inputFiles/phaseField/PhaseFieldFracture_SpectralSplit.xml +++ /dev/null @@ -1,198 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 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 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 488e537950c48635d5b31ec6cdcbc6ad8cbb6a87 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Wed, 29 Jul 2026 08:52:40 -0700 Subject: [PATCH 10/12] fix save history of crack driving force --- .../constitutive/solid/Damage.cpp | 3 +++ .../constitutive/solid/Damage.hpp | 25 +++++++++++-------- .../constitutive/solid/DamageSpectral.hpp | 19 +++++--------- .../constitutive/solid/DamageVolDev.hpp | 20 +++++---------- .../constitutive/solid/SolidFields.hpp | 8 ++++++ 5 files changed, 38 insertions(+), 37 deletions(-) diff --git a/src/coreComponents/constitutive/solid/Damage.cpp b/src/coreComponents/constitutive/solid/Damage.cpp index 50133f15dfe..acd9b3a7a9f 100644 --- a/src/coreComponents/constitutive/solid/Damage.cpp +++ b/src/coreComponents/constitutive/solid/Damage.cpp @@ -84,6 +84,8 @@ Damage< BASE >::Damage( string const & name, Group * const parent ): this->template registerField< fields::solid::crackDrivingForce >( &m_crackDrivingForce ); + this->template registerField< fields::solid::oldCrackDrivingForce >( &m_oldCrackDrivingForce ); + this->template registerField< fields::solid::volStrain >( &m_volStrain ); this->template registerField< fields::solid::extDrivingForce >( &m_extDrivingForce ); @@ -133,6 +135,7 @@ void Damage< BASE >::allocateConstitutiveData( Group & parent, localIndex const m_oldDamage.resize( 0, numPts ); m_damageGrad.resize( 0, numPts, 3 ); 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() ); diff --git a/src/coreComponents/constitutive/solid/Damage.hpp b/src/coreComponents/constitutive/solid/Damage.hpp index 92f7a8e3a74..b222285c77a 100644 --- a/src/coreComponents/constitutive/solid/Damage.hpp +++ b/src/coreComponents/constitutive/solid/Damage.hpp @@ -97,6 +97,7 @@ class DamageUpdates : public UPDATE_BASE arrayView2d< real64 > const & inputOldDamage, arrayView3d< real64 > const & inputDamageGrad, arrayView2d< real64 > const & inputCrackDrivingForce, + arrayView2d< real64 > const & inputOldCrackDrivingForce, arrayView2d< real64 > const & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, @@ -115,6 +116,7 @@ class DamageUpdates : public UPDATE_BASE m_oldDamage( inputOldDamage ), m_damageGrad( inputDamageGrad ), m_crackDrivingForce( inputCrackDrivingForce ), + m_oldCrackDrivingForce( inputOldCrackDrivingForce ), m_volStrain( inputVolumetricStrain ), m_extDrivingForce ( inputExtDrivingForce ), m_lengthScale( inputLengthScale ), @@ -299,6 +301,11 @@ class DamageUpdates : public UPDATE_BASE m_volStrain( k, q ) = traceOfStrain; + // 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; @@ -344,20 +351,10 @@ class DamageUpdates : public UPDATE_BASE } - // TODO: The code below assumes the crack driving force will never be - // evaluated in a non-converged / garbage configuration. - GEOS_HOST_DEVICE virtual real64 getCrackDrivingForce( localIndex const k, localIndex const q ) const { - real64 const sed = SolidBaseUpdates::getStrainEnergyDensity( k, q ); - - if( sed > m_crackDrivingForce( k, q ) ) - { - m_crackDrivingForce( k, q ) = sed; - } - return m_crackDrivingForce( k, q ); } @@ -418,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 @@ -432,6 +430,9 @@ class DamageUpdates : public UPDATE_BASE /// The strain energy density to drive fracture at the quadrature point 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; @@ -509,6 +510,7 @@ class Damage : public BASE m_oldDamage.toView(), m_damageGrad.toView(), m_crackDrivingForce.toView(), + m_oldCrackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, @@ -560,6 +562,9 @@ class Damage : public BASE /// The strain energy density to drive fracture at the quadrature point 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; diff --git a/src/coreComponents/constitutive/solid/DamageSpectral.hpp b/src/coreComponents/constitutive/solid/DamageSpectral.hpp index 13be1893cab..c6861d40da9 100644 --- a/src/coreComponents/constitutive/solid/DamageSpectral.hpp +++ b/src/coreComponents/constitutive/solid/DamageSpectral.hpp @@ -41,6 +41,7 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > arrayView2d< real64 > const & inputOldDamage, arrayView3d< real64 > const & inputDamageGrad, arrayView2d< real64 > const & inputCrackDrivingForce, + arrayView2d< real64 > const & inputOldCrackDrivingForce, arrayView2d< real64 > const & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, @@ -54,7 +55,7 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > arrayView1d< real64 > const & inputDeltaCoefficient, arrayView1d< real64 > const & inputBiotCoefficient, PARAMS && ... baseParams ): - DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputCrackDrivingForce, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, + DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputCrackDrivingForce, inputOldCrackDrivingForce, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputFractureModelType, inputLocalDissipationOption, inputTensileStrength, inputCompressiveStrength, inputDeltaCoefficient, inputBiotCoefficient, @@ -68,6 +69,7 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > using DamageUpdates< UPDATE_BASE >::getDegradationValue; 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; @@ -202,10 +204,7 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > real64 const sed = 0.5 * lambda * tracePlus * tracePlus + mu * doubleContraction( positivePartOfStrain, positivePartOfStrain ); - if( sed > m_crackDrivingForce( k, q ) ) - { - m_crackDrivingForce( k, q ) = sed; - } + m_crackDrivingForce( k, q ) = fmax( sed, m_oldCrackDrivingForce( k, q ) ); } @@ -220,14 +219,6 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } - - GEOS_HOST_DEVICE - virtual real64 getCrackDrivingForce( localIndex const k, - localIndex const q ) const override final - { - return m_crackDrivingForce( k, q ); - } - }; @@ -242,6 +233,7 @@ class DamageSpectral : public Damage< BASE > using Damage< BASE >::m_oldDamage; using Damage< BASE >::m_damageGrad; 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; @@ -269,6 +261,7 @@ class DamageSpectral : public Damage< BASE > m_oldDamage.toView(), m_damageGrad.toView(), m_crackDrivingForce.toView(), + m_oldCrackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, diff --git a/src/coreComponents/constitutive/solid/DamageVolDev.hpp b/src/coreComponents/constitutive/solid/DamageVolDev.hpp index 9a0486722d9..1c990aee9b2 100644 --- a/src/coreComponents/constitutive/solid/DamageVolDev.hpp +++ b/src/coreComponents/constitutive/solid/DamageVolDev.hpp @@ -39,6 +39,7 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > arrayView2d< real64 > const & inputOldDamage, arrayView3d< real64 > const & inputDamageGrad, arrayView2d< real64 > const & inputCrackDrivingForce, + arrayView2d< real64 > const & inputOldCrackDrivingForce, arrayView2d< real64 > const & inputVolumetricStrain, arrayView2d< real64 > const & inputExtDrivingForce, real64 const & inputLengthScale, @@ -52,7 +53,7 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > arrayView1d< real64 > const & inputDeltaCoefficient, arrayView1d< real64 > const & inputBiotCoefficient, PARAMS && ... baseParams ): - DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputCrackDrivingForce, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, + DamageUpdates< UPDATE_BASE >( inputNewDamage, inputOldDamage, inputDamageGrad, inputCrackDrivingForce, inputOldCrackDrivingForce, inputVolumetricStrain, inputExtDrivingForce, inputLengthScale, inputCriticalFractureEnergy, inputcriticalStrainEnergy, inputDegradationLowerLimit, inputFractureModelType, inputLocalDissipationOption, inputTensileStrength, inputCompressiveStrength, inputDeltaCoefficient, inputBiotCoefficient, @@ -67,6 +68,7 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > using DamageUpdates< UPDATE_BASE >::getDegradationValue; 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; @@ -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_crackDrivingForce( k, q ) ) - { - m_crackDrivingForce( k, q ) = sed; - } - } - - - GEOS_HOST_DEVICE - virtual real64 getCrackDrivingForce( localIndex const k, - localIndex const q ) const override final - { - return m_crackDrivingForce( k, q ); + m_crackDrivingForce( k, q ) = fmax( sed, m_oldCrackDrivingForce( k, q ) ); } }; @@ -175,6 +165,7 @@ class DamageVolDev : public Damage< BASE > using Damage< BASE >::m_oldDamage; using Damage< BASE >::m_damageGrad; 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; @@ -202,6 +193,7 @@ class DamageVolDev : public Damage< BASE > m_oldDamage.toView(), m_damageGrad.toView(), m_crackDrivingForce.toView(), + m_oldCrackDrivingForce.toView(), m_volStrain.toView(), m_extDrivingForce.toView(), m_lengthScale, diff --git a/src/coreComponents/constitutive/solid/SolidFields.hpp b/src/coreComponents/constitutive/solid/SolidFields.hpp index 3593e857117..bfb1eda4c51 100644 --- a/src/coreComponents/constitutive/solid/SolidFields.hpp +++ b/src/coreComponents/constitutive/solid/SolidFields.hpp @@ -162,6 +162,14 @@ DECLARE_FIELD( crackDrivingForce, WRITE_AND_READ, "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", array2d< real64 >, From 801b32bae8f9452b258e9caada177c912d72dbfe Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Tue, 4 Aug 2026 16:14:33 -0700 Subject: [PATCH 11/12] enabled multi-region run with certain defined regions have no damage dof --- inputFiles/phaseField/PhaseFieldFracture.ats | 8 + .../PhaseFieldFracture_MultiRegion.xml | 191 ++++++++++++++++++ .../simplePDE/PhaseFieldDamageFEM.cpp | 3 +- 3 files changed, 201 insertions(+), 1 deletion(-) create mode 100644 inputFiles/phaseField/PhaseFieldFracture_MultiRegion.xml diff --git a/inputFiles/phaseField/PhaseFieldFracture.ats b/inputFiles/phaseField/PhaseFieldFracture.ats index f22f302d110..44087dbea67 100644 --- a/inputFiles/phaseField/PhaseFieldFracture.ats +++ b/inputFiles/phaseField/PhaseFieldFracture.ats @@ -59,6 +59,14 @@ decks = [ 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_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/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp index 1f735fc070b..bdd978ad0e0 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp +++ b/src/coreComponents/physicsSolvers/simplePDE/PhaseFieldDamageFEM.cpp @@ -291,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]; From 41ea53a8751ee606091062c8de25d2aee9dc3445 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Thu, 6 Aug 2026 16:10:48 -0700 Subject: [PATCH 12/12] Add phase-field benchmark --- .../phaseField/benchmark/.gitattributes | 4 + ...ieldFracture_SingleEdgeNotchShear_main.xml | 103 +++++++++++++ ...ldFracture_SingleEdgeNotchTension_main.xml | 110 ++++++++++++++ ...haseFieldFracture_SingleEdgeNotch_base.xml | 99 +++++++++++++ .../meshes/singleEdgeCrackSquareShear.vtu | 3 + .../meshes/singleEdgeCrackSquareTension.vtu | 3 + ...eFieldFracture_ThreePointsBending_base.xml | 95 ++++++++++++ ...eFieldFracture_ThreePointsBending_main.xml | 112 ++++++++++++++ .../mesh/threePointsBendingSingleNotch.vtu | 3 + ...Fracture_ThreePointsBendWithHoles_base.xml | 106 +++++++++++++ ...Fracture_ThreePointsBendWithHoles_main.xml | 139 ++++++++++++++++++ .../mesh/threePointsBendingWithHoles.vtu | 3 + 12 files changed, 780 insertions(+) create mode 100644 inputFiles/phaseField/benchmark/.gitattributes create mode 100644 inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotchShear_main.xml create mode 100644 inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotchTension_main.xml create mode 100644 inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotch_base.xml create mode 100644 inputFiles/phaseField/benchmark/singleEdgeNotch/meshes/singleEdgeCrackSquareShear.vtu create mode 100644 inputFiles/phaseField/benchmark/singleEdgeNotch/meshes/singleEdgeCrackSquareTension.vtu create mode 100644 inputFiles/phaseField/benchmark/threePointsBending/PhaseFieldFracture_ThreePointsBending_base.xml create mode 100644 inputFiles/phaseField/benchmark/threePointsBending/PhaseFieldFracture_ThreePointsBending_main.xml create mode 100644 inputFiles/phaseField/benchmark/threePointsBending/mesh/threePointsBendingSingleNotch.vtu create mode 100644 inputFiles/phaseField/benchmark/threePointsBendingWithHoles/PhaseFieldFracture_ThreePointsBendWithHoles_base.xml create mode 100644 inputFiles/phaseField/benchmark/threePointsBendingWithHoles/PhaseFieldFracture_ThreePointsBendWithHoles_main.xml create mode 100644 inputFiles/phaseField/benchmark/threePointsBendingWithHoles/mesh/threePointsBendingWithHoles.vtu diff --git a/inputFiles/phaseField/benchmark/.gitattributes b/inputFiles/phaseField/benchmark/.gitattributes new file mode 100644 index 00000000000..95ec89ceb79 --- /dev/null +++ b/inputFiles/phaseField/benchmark/.gitattributes @@ -0,0 +1,4 @@ +singleEdgeNotch/meshes/singleEdgeCrackSquareShear.vtu filter=lfs diff=lfs merge=lfs -text +singleEdgeNotch/meshes/singleEdgeCrackSquareTension.vtu filter=lfs diff=lfs merge=lfs -text +threePointsBending/mesh/threePointsBendingSingleNotch.vtu filter=lfs diff=lfs merge=lfs -text +threePointsBendingWithHoles/mesh/threePointsBendingWithHoles.vtu filter=lfs diff=lfs merge=lfs -text diff --git a/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotchShear_main.xml b/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotchShear_main.xml new file mode 100644 index 00000000000..5925672a3fe --- /dev/null +++ b/inputFiles/phaseField/benchmark/singleEdgeNotch/PhaseFieldFracture_SingleEdgeNotchShear_main.xml @@ -0,0 +1,103 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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