From d5de7f366e55b8a6328e3791cbf7a8e731eb2184 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 11:08:40 +0200 Subject: [PATCH 01/48] Update SolidMechanicsLagrangianFEM.cpp --- .../solidMechanics/SolidMechanicsLagrangianFEM.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 5794878d59e..8f60e9a5d03 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -946,20 +946,28 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >(); + constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, , [&] (auto const solidModel) + { + using SOLID_TYPE = decltype( solidModel ); + finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) { using FE_TYPE = decltype( finiteElement ); - AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< CellElementSubRegion, FE_TYPE, parallelDevicePolicy<> >( nodeManager, + AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< CellElementSubRegion, FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager, mesh.getEdgeManager(), mesh.getFaceManager(), subRegion, finiteElement, + solidModel, disp, strain ); } ); + } ); + + } ); } ); From 662836f3d6d232b1b4d751bc871807a18021de99 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 11:39:35 +0200 Subject: [PATCH 02/48] Update StrainHelper.hpp --- .../solidMechanics/kernels/StrainHelper.hpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index d77b1cb1e26..09fd4c90d2a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -33,9 +33,11 @@ namespace geos * @class AverageStrainOverQuadraturePoints * @tparam SUBREGION_TYPE the subRegion type * @tparam FE_TYPE the finite element type + * @tparam SOLID_TYPE the solid mechanics constitutuve type */ template< typename SUBREGION_TYPE, - typename FE_TYPE > + typename FE_TYPE, + typename SOLID_TYPE > class AverageStrainOverQuadraturePoints : public AverageOverQuadraturePointsBase< SUBREGION_TYPE, FE_TYPE > @@ -65,6 +67,7 @@ class AverageStrainOverQuadraturePoints : FaceManager const & faceManager, SUBREGION_TYPE const & elementSubRegion, FE_TYPE const & finiteElementSpace, + SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, fields::solidMechanics::arrayView2dLayoutStrain const avgStrain ): Base( nodeManager, @@ -72,6 +75,7 @@ class AverageStrainOverQuadraturePoints : faceManager, elementSubRegion, finiteElementSpace ), + m_solidUpdate(solidModel.createKernelUpdates()), m_displacement( displacement ), m_avgStrain( avgStrain ) {} @@ -160,6 +164,9 @@ class AverageStrainOverQuadraturePoints : protected: + /// The material + typename SOLID_TYPE::kernelWrapper const m_solidUpdate; + /// The displacement solution fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const m_displacement; From df505c225022543be10f07b2767ebcb6a9ebbc9b Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 12:19:21 +0200 Subject: [PATCH 03/48] Update SolidMechanicsFields.hpp --- .../solidMechanics/SolidMechanicsFields.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp index d12b9b4e243..5733304fa32 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp @@ -87,6 +87,14 @@ DECLARE_FIELD( strain, WRITE_AND_READ, "Average strain in cell" ); +DECLARE_FIELD( plasticStrain, + "plasticStrain", + array2dLayoutStrain, + 0, + LEVEL_0, + WRITE_AND_READ, + "Average plastic strain in cell" ); + DECLARE_FIELD( incrementalBubbleDisplacement, "incrementalBubbleDisplacement", array2d< real64 >, From 3b0bd7b6a4412b491ecf14be683732d0494c182e Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 12:26:08 +0200 Subject: [PATCH 04/48] Update SolidMechanicsLagrangianFEM.cpp --- .../solidMechanics/SolidMechanicsLagrangianFEM.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 8f60e9a5d03..ae3fd2bc201 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -159,6 +159,7 @@ void SolidMechanicsLagrangianFEM::registerDataOnMesh( Group & meshBodies ) setConstitutiveNamesCallSuper( subRegion ); subRegion.registerField< solidMechanics::strain >( getName() ).reference().resizeDimension< 1 >( 6 ); + subRegion.registerField< solidMechanics::plasticStrain >( getName() ).reference().resizeDimension< 1 >( 6 ); } ); NodeManager & nodes = meshLevel.getNodeManager(); @@ -945,6 +946,7 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS constitutiveRelation.saveConvergedState(); solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >(); + solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >(); constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, , [&] (auto const solidModel) { @@ -961,7 +963,8 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS finiteElement, solidModel, disp, - strain ); + strain, + plasticStrain ); } ); From 78e4db41916cb1c136857e633011bd6f41c20e21 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 12:28:52 +0200 Subject: [PATCH 05/48] Update StrainHelper.hpp --- .../solidMechanics/kernels/StrainHelper.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 09fd4c90d2a..ec97b907f2c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -69,7 +69,8 @@ class AverageStrainOverQuadraturePoints : FE_TYPE const & finiteElementSpace, SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, - fields::solidMechanics::arrayView2dLayoutStrain const avgStrain ): + fields::solidMechanics::arrayView2dLayoutStrain const avgStrain + fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain): Base( nodeManager, edgeManager, faceManager, @@ -77,7 +78,8 @@ class AverageStrainOverQuadraturePoints : finiteElementSpace ), m_solidUpdate(solidModel.createKernelUpdates()), m_displacement( displacement ), - m_avgStrain( avgStrain ) + m_avgStrain( avgStrain ), + m_avgPlasticStrain( avgPlasticStrain ) {} /** @@ -173,6 +175,9 @@ class AverageStrainOverQuadraturePoints : /// The average strain fields::solidMechanics::arrayView2dLayoutStrain const m_avgStrain; + /// The average plastic strain + fields::solidMechanics::arrayView2dLayoutStrain const m_avgPlasticStrain; + }; From a23f26f29283face85dad941bc2300cd315dbfe8 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 13:45:03 +0200 Subject: [PATCH 06/48] Update StrainHelper.hpp --- .../solidMechanics/kernels/StrainHelper.hpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index ec97b907f2c..1887f17fc3d 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -23,6 +23,7 @@ #include "common/DataTypes.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "finiteElement/FiniteElementDispatch.hpp" +#include "constitutive/ConsitutivePassThru.hpp" #include "mesh/CellElementSubRegion.hpp" #include "mesh/utilities/AverageOverQuadraturePointsKernel.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" @@ -111,6 +112,7 @@ class AverageStrainOverQuadraturePoints : for( int icomp = 0; icomp < 6; ++icomp ) { m_avgStrain[k][icomp] = 0.0; + m_avgPlasticStrain[k][icomp] = 0.0; } } @@ -132,9 +134,13 @@ class AverageStrainOverQuadraturePoints : real64 strain[6] = {0.0}; FE_TYPE::symmetricGradient( dNdX, stack.uLocal, strain ); + real64 elasticStrain[6]; + m_solidUpdate.getElasticStrain(k, q, elasticStrain); + for( int icomp = 0; icomp < 6; ++icomp ) { m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k]; + m_avgPlasticStrain[k][icomp] += detJxW*(strain[icomp] - elasticStrain[icomp])/m_elementVolume[k]; } } @@ -205,6 +211,7 @@ class AverageStrainOverQuadraturePointsKernelFactory */ template< typename SUBREGION_TYPE, typename FE_TYPE, + typename SOLID_TYPE, typename POLICY > static void createAndLaunch( NodeManager & nodeManager, @@ -212,14 +219,16 @@ class AverageStrainOverQuadraturePointsKernelFactory FaceManager const & faceManager, SUBREGION_TYPE const & elementSubRegion, FE_TYPE const & finiteElementSpace, + SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, - fields::solidMechanics::arrayView2dLayoutStrain const avgStrain ) + fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, + fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain,) { - AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE > + AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE, SOLID_TYPE > kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace, - displacement, avgStrain ); + solidModel, displacement, avgStrain ); - AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >::template + AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE, SOLID_TYPE >::template kernelLaunch< POLICY >( elementSubRegion.size(), kernel ); } }; From c2325c3a0bbd0c4ee49f36760010a7110e5f8bde Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 13:51:26 +0200 Subject: [PATCH 07/48] Update StrainHelper.hpp --- .../physicsSolvers/solidMechanics/kernels/StrainHelper.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 1887f17fc3d..8806d0891a3 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -23,7 +23,7 @@ #include "common/DataTypes.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "finiteElement/FiniteElementDispatch.hpp" -#include "constitutive/ConsitutivePassThru.hpp" +#include "constitutive/ConstitutivePassThru.hpp" #include "mesh/CellElementSubRegion.hpp" #include "mesh/utilities/AverageOverQuadraturePointsKernel.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" From 27a4801bd2f4eccbc3a3f3cf74fa40166aaa8f7c Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 13:58:14 +0200 Subject: [PATCH 08/48] Update StrainHelper.hpp --- .../physicsSolvers/solidMechanics/kernels/StrainHelper.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 8806d0891a3..073097c0454 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -70,7 +70,7 @@ class AverageStrainOverQuadraturePoints : FE_TYPE const & finiteElementSpace, SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, - fields::solidMechanics::arrayView2dLayoutStrain const avgStrain + fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain): Base( nodeManager, edgeManager, From 18a94b21cc17084cc98fc65b20055b7af57771ac Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 14:03:47 +0200 Subject: [PATCH 09/48] Update StrainHelper.hpp --- .../physicsSolvers/solidMechanics/kernels/StrainHelper.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 073097c0454..46df9fd7938 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -200,6 +200,7 @@ class AverageStrainOverQuadraturePointsKernelFactory * @brief Create a new kernel and launch * @tparam SUBREGION_TYPE the subRegion type * @tparam FE_TYPE the finite element type + * @tparam SOLID_TYPE the constitutive type * @tparam POLICY the kernel policy * @param nodeManager the node manager * @param edgeManager the edge manager @@ -222,11 +223,11 @@ class AverageStrainOverQuadraturePointsKernelFactory SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, - fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain,) + fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain) { AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE, SOLID_TYPE > kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace, - solidModel, displacement, avgStrain ); + solidModel, displacement, avgStrain, avgPlasticStrain ); AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE, SOLID_TYPE >::template kernelLaunch< POLICY >( elementSubRegion.size(), kernel ); From c7e03b916e01a978ee413dbae035a0d5ed612102 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 14:22:49 +0200 Subject: [PATCH 10/48] Update SolidMechanicsLagrangianFEM.cpp --- .../solidMechanics/SolidMechanicsLagrangianFEM.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index ae3fd2bc201..88f051dacd1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -948,7 +948,7 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >(); solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >(); - constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, , [&] (auto const solidModel) + constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto const solidModel) { using SOLID_TYPE = decltype( solidModel ); From 2f4e691cc7e9cfb321baa86a2e38674fc1f17637 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 14:36:28 +0200 Subject: [PATCH 11/48] Update SolidMechanicsLagrangianFEM.cpp --- .../solidMechanics/SolidMechanicsLagrangianFEM.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 88f051dacd1..a6b03ee022c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -956,7 +956,7 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) { using FE_TYPE = decltype( finiteElement ); - AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< CellElementSubRegion, FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager, + AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, parallelDevicePolicy<> >( nodeManager, mesh.getEdgeManager(), mesh.getFaceManager(), subRegion, From 3cf18caa6dbcbb4d44a7a2d225da00d3af84d8da Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 14:38:54 +0200 Subject: [PATCH 12/48] Update StrainHelper.hpp --- .../solidMechanics/kernels/StrainHelper.hpp | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 46df9fd7938..e832c2c4d00 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -36,17 +36,16 @@ namespace geos * @tparam FE_TYPE the finite element type * @tparam SOLID_TYPE the solid mechanics constitutuve type */ -template< typename SUBREGION_TYPE, - typename FE_TYPE, +template< typename FE_TYPE, typename SOLID_TYPE > class AverageStrainOverQuadraturePoints : - public AverageOverQuadraturePointsBase< SUBREGION_TYPE, + public AverageOverQuadraturePointsBase< CellElementSubRegion, FE_TYPE > { public: /// Alias for the base class; - using Base = AverageOverQuadraturePointsBase< SUBREGION_TYPE, + using Base = AverageOverQuadraturePointsBase< CellElementSubRegion, FE_TYPE >; using Base::m_elementVolume; @@ -66,7 +65,7 @@ class AverageStrainOverQuadraturePoints : AverageStrainOverQuadraturePoints( NodeManager & nodeManager, EdgeManager const & edgeManager, FaceManager const & faceManager, - SUBREGION_TYPE const & elementSubRegion, + CellElementSubRegion const & elementSubRegion, FE_TYPE const & finiteElementSpace, SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, @@ -210,26 +209,25 @@ class AverageStrainOverQuadraturePointsKernelFactory * @param property the property at quadrature points * @param averageProperty the property averaged over quadrature points */ - template< typename SUBREGION_TYPE, - typename FE_TYPE, + template< typename FE_TYPE, typename SOLID_TYPE, typename POLICY > static void createAndLaunch( NodeManager & nodeManager, EdgeManager const & edgeManager, FaceManager const & faceManager, - SUBREGION_TYPE const & elementSubRegion, + CellElementSubRegion const & elementSubRegion, FE_TYPE const & finiteElementSpace, SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain) { - AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE, SOLID_TYPE > + AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE > kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace, solidModel, displacement, avgStrain, avgPlasticStrain ); - AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE, SOLID_TYPE >::template + AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >::template kernelLaunch< POLICY >( elementSubRegion.size(), kernel ); } }; From 4284d711d538c6cd3fc4df6c21d20150a8e5d91a Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 14:42:04 +0200 Subject: [PATCH 13/48] Update SolidMechanicsLagrangianFEM.cpp --- .../solidMechanics/SolidMechanicsLagrangianFEM.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index a6b03ee022c..449b6a197af 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -948,9 +948,8 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >(); solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >(); - constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto const solidModel) + constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto solidModel) { - using SOLID_TYPE = decltype( solidModel ); finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) From cd42322e1d06fd707621a0048fdab55cce3e1092 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 15:03:09 +0200 Subject: [PATCH 14/48] Update SolidMechanicsLagrangianFEM.cpp --- .../solidMechanics/SolidMechanicsLagrangianFEM.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 449b6a197af..89da0d72b65 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -948,7 +948,7 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >(); solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >(); - constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto solidModel) + constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto & solidModel) { finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); From 6700051029e43032396deb88fba8b0ba33a3c390 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 15:14:50 +0200 Subject: [PATCH 15/48] Update SolidMechanicsLagrangianFEM.cpp --- .../solidMechanics/SolidMechanicsLagrangianFEM.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 89da0d72b65..6773b54e477 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -951,11 +951,13 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto & solidModel) { + using SOLID_TYPE = decltype( solidModel ); + finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) { using FE_TYPE = decltype( finiteElement ); - AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, parallelDevicePolicy<> >( nodeManager, + AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager, mesh.getEdgeManager(), mesh.getFaceManager(), subRegion, From 5c8f164eb2542b0e26a7eed05ba7014ead056e7d Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 15:32:36 +0200 Subject: [PATCH 16/48] Update SolidMechanicsLagrangianFEM.cpp --- .../solidMechanics/SolidMechanicsLagrangianFEM.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 6773b54e477..9cdfa11c000 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -951,7 +951,7 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto & solidModel) { - using SOLID_TYPE = decltype( solidModel ); + using SOLID_TYPE = TYPEOFREF( solidModel ); finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) From 7a51c307d862d1c0edec81e28c4458c5b2e30cf4 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 15:41:48 +0200 Subject: [PATCH 17/48] Update StrainHelper.hpp --- .../physicsSolvers/solidMechanics/kernels/StrainHelper.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index e832c2c4d00..7bdded4867f 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -172,7 +172,7 @@ class AverageStrainOverQuadraturePoints : protected: /// The material - typename SOLID_TYPE::kernelWrapper const m_solidUpdate; + typename SOLID_TYPE::KernelWrapper const m_solidUpdate; /// The displacement solution fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const m_displacement; From ca89c990b7ab1d3bfd6e9e28816f2cb240b15132 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 3 Oct 2024 16:01:09 +0200 Subject: [PATCH 18/48] Update StrainHelper.hpp --- .../physicsSolvers/solidMechanics/kernels/StrainHelper.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 7bdded4867f..7df079038fd 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -133,7 +133,7 @@ class AverageStrainOverQuadraturePoints : real64 strain[6] = {0.0}; FE_TYPE::symmetricGradient( dNdX, stack.uLocal, strain ); - real64 elasticStrain[6]; + real64 elasticStrain[6] = {0.0}; m_solidUpdate.getElasticStrain(k, q, elasticStrain); for( int icomp = 0; icomp < 6; ++icomp ) From cf6fe9851cdec45ffd422160048c95fd1c1f4d18 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Fri, 11 Oct 2024 15:30:39 -0700 Subject: [PATCH 19/48] initial implementation, behavior is almost right but need to debug initial states --- .../constitutive/solid/ElasticIsotropic.hpp | 22 + .../ElasticIsotropicPressureDependent.hpp | 66 ++ .../constitutive/solid/SolidBase.hpp | 18 + "src/coreComponents/constitutive/solid/\177" | 573 ++++++++++++++++++ .../SolidMechanicsLagrangianFEM.cpp | 7 +- .../solidMechanics/kernels/StrainHelper.hpp | 23 +- 6 files changed, 701 insertions(+), 8 deletions(-) create mode 100644 "src/coreComponents/constitutive/solid/\177" diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index b6388f8e907..12975c80ae8 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -135,6 +135,11 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates localIndex const q, real64 ( &elasticStrain )[6] ) const override final; + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + GEOS_HOST_DEVICE virtual real64 getBulkModulus( localIndex const k ) const override final { @@ -228,6 +233,23 @@ void ElasticIsotropicUpdates::getElasticStrain( localIndex const k, elasticStrain[5] = m_newStress[k][q][5] / m_shearModulus[k]; } +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const +{ + real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] ); + real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] ); + + elasticStrainInc[0] = ( (m_newStress[k][q][0] - m_oldStress[k][q][0]) - nu*(m_newStress[k][q][1] - m_oldStress[k][q][1]) - nu*(m_newStress[k][q][2] - m_oldStress[k][q][2]))/E; + elasticStrainInc[1] = (-nu*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_newStress[k][q][1] - m_oldStress[k][q][1]) - nu*(m_newStress[k][q][2] - m_oldStress[k][q][2]))/E; + elasticStrainInc[2] = (-nu*(m_newStress[k][q][0] - m_oldStress[k][q][0]) - nu*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_newStress[k][q][2] - m_oldStress[k][q][2]))/E; + + elasticStrainInc[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_shearModulus[k]; + elasticStrainInc[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_shearModulus[k]; + elasticStrainInc[5] = (m_newStress[k][q][5] - m_oldStress[k][q][5]) / m_shearModulus[k]; +} GEOS_HOST_DEVICE inline diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 3718765d3b2..671a002421a 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -115,6 +115,11 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates localIndex const q, real64 ( &elasticStrain )[6] ) const override final; + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + GEOS_HOST_DEVICE virtual void viscousStateUpdate( localIndex const k, localIndex const q, @@ -223,6 +228,67 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrain( localIndex cons } +GEOS_HOST_DEVICE +inline +void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const +{ + real64 const mu = m_shearModulus[k]; + real64 const p0 = m_refPressure; + real64 const eps_v0 = m_refStrainVol; + real64 const Cr = m_recompressionIndex[k]; + real64 deviator[6]; + real64 stress[6]; + real64 P; + real64 Q; + real64 elasticStrainVol; + real64 elasticStrainDev; + + for( localIndex i=0; i<6; ++i ) + { + stress[i] = m_newStress[k][q][i]; + } + + twoInvariant::stressDecomposition( stress, + P, + Q, + deviator ); + + elasticStrainVol = std::log( P/p0 ) * Cr * (-1.0) + eps_v0; + elasticStrainDev = Q/3./mu; + + twoInvariant::strainRecomposition( elasticStrainVol, + elasticStrainDev, + deviator, + elasticStrainInc ); + + real64 oldStrain[6]; + for( localIndex i=0; i<6; ++i ) + { + stress[i] = m_oldStress[k][q][i]; + } + + twoInvariant::stressDecomposition( stress, + P, + Q, + deviator ); + + elasticStrainVol = std::log( P/p0 ) * Cr * (-1.0) + eps_v0; + elasticStrainDev = Q/3./mu; + + twoInvariant::strainRecomposition( elasticStrainVol, + elasticStrainDev, + deviator, + oldStrain ); + + for (localIndex i = 0; i<6; ++i) + { + elasticStrainInc[i] -= oldStrain[i]; + } +} + + GEOS_HOST_DEVICE inline void ElasticIsotropicPressureDependentUpdates::smallStrainUpdate( localIndex const k, diff --git a/src/coreComponents/constitutive/solid/SolidBase.hpp b/src/coreComponents/constitutive/solid/SolidBase.hpp index cbcf7ffaca8..f8cd559f645 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.hpp +++ b/src/coreComponents/constitutive/solid/SolidBase.hpp @@ -337,6 +337,24 @@ class SolidBaseUpdates GEOS_ERROR( "getElasticStrain() not implemented for this model" ); } + /** + * @brief Return the current elastic strain increment at a given material point (small-strain interface) + * + * @param k the element inex + * @param q the quadrature index + * @param elasticStrainInc Current elastic strain increment + */ + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc )[6] ) const + { + GEOS_UNUSED_VAR( k ); + GEOS_UNUSED_VAR( q ); + GEOS_UNUSED_VAR( elasticStrainInc ); + GEOS_ERROR( "getElasticStrainInc() not implemented for this model (called when computing plastic strains)" ); + } + /** * @brief Perform a viscous (rate-dependent) state update * diff --git "a/src/coreComponents/constitutive/solid/\177" "b/src/coreComponents/constitutive/solid/\177" new file mode 100644 index 00000000000..af8484f08ae --- /dev/null +++ "b/src/coreComponents/constitutive/solid/\177" @@ -0,0 +1,573 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ElasticIsotropic.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPIC_HPP_ +#define GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPIC_HPP_ + +#include "SolidBase.hpp" +#include "PropertyConversions.hpp" +#include "SolidModelDiscretizationOpsIsotropic.hpp" +#include "constitutive/ExponentialRelation.hpp" +#include "LvArray/src/tensorOps.hpp" + +namespace geos +{ + +namespace constitutive +{ + +/** + * @class ElasticIsotropicUpdates + * + * Class to provide elastic isotropic material updates that may be + * called from a kernel function. + */ +class ElasticIsotropicUpdates : public SolidBaseUpdates +{ +public: + /** + * @brief Constructor + * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. + * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. + * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. + * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. + * @param[in] disableInelasticity Flag to disable plasticity for inelastic models + */ + ElasticIsotropicUpdates( arrayView1d< real64 const > const & bulkModulus, + arrayView1d< real64 const > const & shearModulus, + arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView3d< real64, solid::STRESS_USD > const & newStress, + arrayView3d< real64, solid::STRESS_USD > const & oldStress, + const bool & disableInelasticity ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + m_bulkModulus( bulkModulus ), + m_shearModulus( shearModulus ) + {} + + /// Deleted default constructor + ElasticIsotropicUpdates() = delete; + + /// Default copy constructor + ElasticIsotropicUpdates( ElasticIsotropicUpdates const & ) = default; + + /// Default move constructor + ElasticIsotropicUpdates( ElasticIsotropicUpdates && ) = default; + + /// Deleted copy assignment operator + ElasticIsotropicUpdates & operator=( ElasticIsotropicUpdates const & ) = delete; + + /// Deleted move assignment operator + ElasticIsotropicUpdates & operator=( ElasticIsotropicUpdates && ) = delete; + + /// Use the "isotropic" form of inner product compression + using DiscretizationOps = SolidModelDiscretizationOpsIsotropic; + + /// Use base version of saveConvergedState + using SolidBaseUpdates::saveConvergedState; + + GEOS_HOST_DEVICE + virtual void smallStrainNoStateUpdate_StressOnly( localIndex const k, + localIndex const q, + real64 const ( &totalStrain )[6], + real64 ( &stress )[6] ) const override final; + + GEOS_HOST_DEVICE + virtual void smallStrainNoStateUpdate( localIndex const k, + localIndex const q, + real64 const ( &totalStrain )[6], + real64 ( &stress )[6], + real64 ( &stiffness )[6][6] ) const override final; + + GEOS_HOST_DEVICE + virtual void smallStrainNoStateUpdate( localIndex const k, + localIndex const q, + real64 const ( &totalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness ) const final; + + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_StressOnly( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrement )[6], + real64 ( &stress )[6] ) const override; + + GEOS_HOST_DEVICE + void smallStrainUpdate( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrement )[6], + real64 ( &stress )[6], + real64 ( &stiffness )[6][6] ) const; + + GEOS_HOST_DEVICE + virtual void smallStrainUpdate( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrement )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness ) const; + + GEOS_HOST_DEVICE + virtual void getElasticStiffness( localIndex const k, + localIndex const q, + real64 ( &stiffness )[6][6] ) const override; + + GEOS_HOST_DEVICE + virtual void getElasticStrain( localIndex const k, + localIndex const q, + real64 ( &elasticStrain )[6] ) const override final; + + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + + GEOS_HOST_DEVICE + virtual real64 getBulkModulus( localIndex const k ) const override final + { + return m_bulkModulus[k]; + } + + GEOS_HOST_DEVICE + virtual real64 getShearModulus( localIndex const k ) const override final + { + return m_shearModulus[k]; + } + + GEOS_HOST_DEVICE + virtual void viscousStateUpdate( localIndex const k, + localIndex const q, + real64 beta ) const override; + + // TODO: confirm hyper stress/strain measures before activatiing + + /* + GEOS_HOST_DEVICE + virtual void hyperUpdate( localIndex const k, + localIndex const q, + real64 const ( & FminusI )[3][3], + real64 ( & stress )[6] ) const override final; + + GEOS_HOST_DEVICE + virtual void hyperUpdate( localIndex const k, + localIndex const q, + real64 const ( & FminusI )[3][3], + real64 ( & stress )[6], + real64 ( & stiffness )[6][6] ) const override final; + */ + +protected: + + /// A reference to the ArrayView holding the bulk modulus for each element. + arrayView1d< real64 const > const m_bulkModulus; + + /// A reference to the ArrayView holding the shear modulus for each element. + arrayView1d< real64 const > const m_shearModulus; + +}; + + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::getElasticStiffness( localIndex const k, + localIndex const q, + real64 ( & stiffness )[6][6] ) const +{ + GEOS_UNUSED_VAR( q ); + real64 const G = m_shearModulus[k]; + real64 const lambda = conversions::bulkModAndShearMod::toFirstLame( m_bulkModulus[k], G ); + + LvArray::tensorOps::fill< 6, 6 >( stiffness, 0 ); + + stiffness[0][0] = lambda + 2*G; + stiffness[0][1] = lambda; + stiffness[0][2] = lambda; + + stiffness[1][0] = lambda; + stiffness[1][1] = lambda + 2*G; + stiffness[1][2] = lambda; + + stiffness[2][0] = lambda; + stiffness[2][1] = lambda; + stiffness[2][2] = lambda + 2*G; + + stiffness[3][3] = G; + stiffness[4][4] = G; + stiffness[5][5] = G; +} + + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::getElasticStrain( localIndex const k, + localIndex const q, + real64 ( & elasticStrain)[6] ) const +{ + real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] ); + real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] ); + + elasticStrain[0] = ( m_newStress[k][q][0] - nu*m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E; + elasticStrain[1] = (-nu*m_newStress[k][q][0] + m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E; + elasticStrain[2] = (-nu*m_newStress[k][q][0] - nu*m_newStress[k][q][1] + m_newStress[k][q][2])/E; + + elasticStrain[3] = m_newStress[k][q][3] / m_shearModulus[k]; + elasticStrain[4] = m_newStress[k][q][4] / m_shearModulus[k]; + elasticStrain[5] = m_newStress[k][q][5] / m_shearModulus[k]; +} + + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::smallStrainNoStateUpdate_StressOnly( localIndex const k, + localIndex const q, + real64 const ( &totalStrain )[6], + real64 ( & stress )[6] ) const +{ + GEOS_UNUSED_VAR( q ); + + real64 const twoG = 2 * m_shearModulus[k]; + real64 const lambda = conversions::bulkModAndShearMod::toFirstLame( m_bulkModulus[k], m_shearModulus[k] ); + real64 const vol = lambda * ( totalStrain[0] + totalStrain[1] + totalStrain[2] ); + + stress[0] = vol + twoG * totalStrain[0]; + stress[1] = vol + twoG * totalStrain[1]; + stress[2] = vol + twoG * totalStrain[2]; + + stress[3] = m_shearModulus[k] * totalStrain[3]; + stress[4] = m_shearModulus[k] * totalStrain[4]; + stress[5] = m_shearModulus[k] * totalStrain[5]; +} + + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::smallStrainNoStateUpdate( localIndex const k, + localIndex const q, + real64 const ( &totalStrain )[6], + real64 ( & stress )[6], + real64 ( & stiffness )[6][6] ) const +{ + smallStrainNoStateUpdate_StressOnly( k, q, totalStrain, stress ); + getElasticStiffness( k, q, stiffness ); +} + + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::smallStrainNoStateUpdate( localIndex const k, + localIndex const q, + real64 const ( &totalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness ) const +{ + smallStrainNoStateUpdate_StressOnly( k, q, totalStrain, stress ); + stiffness.m_bulkModulus = m_bulkModulus[k]; + stiffness.m_shearModulus = m_shearModulus[k]; +} + + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::smallStrainUpdate_StressOnly( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrement )[6], + real64 ( & stress )[6] ) const +{ + GEOS_UNUSED_VAR( timeIncrement ); + smallStrainNoStateUpdate_StressOnly( k, q, strainIncrement, stress ); // stress = incrementalStress + LvArray::tensorOps::add< 6 >( stress, m_oldStress[k][q] ); // stress += m_oldStress + saveStress( k, q, stress ); // m_newStress = stress +} + + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::smallStrainUpdate( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrement )[6], + real64 ( & stress )[6], + real64 ( & stiffness )[6][6] ) const +{ + smallStrainUpdate_StressOnly( k, q, timeIncrement, strainIncrement, stress ); + getElasticStiffness( k, q, stiffness ); +} + + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::smallStrainUpdate( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrement )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness ) const +{ + smallStrainUpdate_StressOnly( k, q, timeIncrement, strainIncrement, stress ); + stiffness.m_bulkModulus = m_bulkModulus[k]; + stiffness.m_shearModulus = m_shearModulus[k]; +} + +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ElasticIsotropicUpdates::viscousStateUpdate( localIndex const k, + localIndex const q, + real64 beta ) const +{ + GEOS_UNUSED_VAR( k ); + GEOS_UNUSED_VAR( q ); + GEOS_UNUSED_VAR( beta ); +} + + +// TODO: need to confirm stress / strain measures before activating hyper inferface +/* + GEOS_HOST_DEVICE + inline + void ElasticIsotropicUpdates::hyperUpdate( localIndex const k, + localIndex const q, + real64 const (&FmI)[3][3], + real64 ( & stress )[ 6 ] ) const + { + GEOS_UNUSED_VAR(q); + + real64 const C1 = 0.5 * m_shearModulus[k]; + real64 const D1 = 0.5 * m_bulkModulus[k]; + real64 const detFm1 = FmI[0][0] + FmI[1][1] + FmI[2][2] + - FmI[1][2]*FmI[2][1] + FmI[1][1]*FmI[2][2] + + FmI[0][2]*(-FmI[2][0] - FmI[1][1]*FmI[2][0] + FmI[1][0]*FmI[2][1]) + + FmI[0][1]*(-FmI[1][0] + FmI[1][2]*FmI[2][0] - FmI[1][0]*FmI[2][2]) + + FmI[0][0]*( FmI[1][1] - FmI[1][2]*FmI[2][1] + FmI[2][2] + FmI[1][1]*FmI[2][2]); + + + real64 const p = -2 * D1 * ( detFm1 + 1.0 ) * detFm1; + real64 devB[6] = { 1/3 * (2 * FmI[0][0] * (2 + FmI[0][0]) - FmI[1][1] * (2 + FmI[1][1]) - FmI[2][2] * (2 + FmI[2][2]) + ++ + 2 * FmI[0][1]*FmI[0][1] + 2 * FmI[0][2] * FmI[0][2] - FmI[1][0] * FmI[1][0] - FmI[1][2] * + +FmI[1][2] - + FmI[2][0] * FmI[2][0] - FmI[2][1] * FmI[2][1]), + 1/3 * (-FmI[0][0] * (2 + FmI[0][0]) + 2 * FmI[1][1] * ( 2 + FmI[1][1]) - FmI[2][2] * (2 + + +FmI[2][2]) - + FmI[0][1]*FmI[0][1] - FmI[0][2]*FmI[0][2] + 2 * FmI[1][0]*FmI[1][0] + 2 * + +FmI[1][2]*FmI[1][2] - FmI[2][0]*FmI[2][0] - FmI[2][1]* + FmI[2][1]), + 1/3 *(-FmI[0][0] * (2 + FmI[0][0]) - FmI[1][1] * (2 + FmI[1][1]) + 2 * FmI[2][2] * (2 + FmI[2][2]) + +- + FmI[0][1]*FmI[0][1] - FmI[0][2]*FmI[0][2] - FmI[1][0]*FmI[1][0] - FmI[1][2]*FmI[1][2] + 2 * + +FmI[2][0]*FmI[2][0] + 2 * FmI[2][1]* + FmI[2][1]), + FmI[1][2] + FmI[1][0] * FmI[2][0] + FmI[2][1] + FmI[1][1]*FmI[2][1] + FmI[1][2]*FmI[2][2], + FmI[0][2] + FmI[2][0] + FmI[0][0] * FmI[2][0] + FmI[0][1]*FmI[2][1] + FmI[0][2]*FmI[2][2], + FmI[0][1] + FmI[1][0] + FmI[0][0] * FmI[1][0] + FmI[0][1]*FmI[1][1] + FmI[0][2]*FmI[1][2] + }; + + real64 const C = 2 * C1 / pow( detFm1 + 1, 2.0/3.0 ); + stress[0] = -p + C * devB[0]; + stress[1] = -p + C * devB[1]; + stress[2] = -p + C * devB[2]; + stress[3] = C * devB[3]; + stress[4] = C * devB[4]; + stress[5] = C * devB[5]; + } + */ + + +/** + * @class ElasticIsotropic + * + * Class to provide an elastic isotropic material response. + */ +class ElasticIsotropic : public SolidBase +{ +public: + + /// Alias for ElasticIsotropicUpdates + using KernelWrapper = ElasticIsotropicUpdates; + + /** + * constructor + * @param[in] name name of the instance in the catalog + * @param[in] parent the group which contains this instance + */ + ElasticIsotropic( string const & name, Group * const parent ); + + /** + * Default Destructor + */ + virtual ~ElasticIsotropic() override; + + /** + * @name Static Factory Catalog members and functions + */ + ///@{ + + /// string name to use for this class in the catalog + static constexpr auto m_catalogNameString = "ElasticIsotropic"; + + /** + * @brief Static catalog string + * @return A string that is used to register/lookup this class in the registry + */ + static std::string catalogName() { return m_catalogNameString; } + + /** + * @brief Get catalog name + * @return Name string + */ + virtual string getCatalogName() const override { return catalogName(); } + + ///@} + + /// Keys for data specified in this class. + struct viewKeyStruct : public SolidBase::viewKeyStruct + { + /// string/key for default bulk modulus + static constexpr char const * defaultBulkModulusString() { return "defaultBulkModulus"; } + + /// string/key for default poisson ratio + static constexpr char const * defaultPoissonRatioString() { return "defaultPoissonRatio"; } + + /// string/key for default shear modulus + static constexpr char const * defaultShearModulusString() { return "defaultShearModulus"; } + + /// string/key for default Young's modulus + static constexpr char const * defaultYoungModulusString() { return "defaultYoungModulus"; } + + /// string/key for bulk modulus + static constexpr char const * bulkModulusString() { return "bulkModulus"; } + + /// string/key for shear modulus + static constexpr char const * shearModulusString() { return "shearModulus"; } + + }; + + /** + * @brief Accessor for bulk modulus + * @return A const reference to arrayView1d containing the bulk + * modulus (at every element). + */ + arrayView1d< real64 > const bulkModulus() { return m_bulkModulus; } + + /** + * @brief Const accessor for bulk modulus + * @return A const reference to arrayView1d containing the bulk + * modulus (at every element). + */ + arrayView1d< real64 const > const bulkModulus() const { return m_bulkModulus; } + + /** + * @brief Accessor for shear modulus + * @return A const reference to arrayView1d containing the shear + * modulus (at every element). + */ + arrayView1d< real64 > const shearModulus() { return m_shearModulus; } + + /** + * @brief Const accessor for shear modulus + * @return A const reference to arrayView1d containing the + * shear modulus (at every element). + */ + arrayView1d< real64 const > const shearModulus() const { return m_shearModulus; } + + GEOS_HOST_DEVICE + virtual arrayView1d< real64 const > getBulkModulus() const override final + { + return m_bulkModulus; + } + GEOS_HOST_DEVICE + virtual arrayView1d< real64 const > getShearModulus() const override final + { + return m_shearModulus; + } + + /** + * @brief Create a instantiation of the ElasticIsotropicUpdate class + * that refers to the data in this. + * @param includeState Flag whether to pass state arrays that may not be needed for "no-state" updates + * @return An instantiation of ElasticIsotropicUpdate. + */ + ElasticIsotropicUpdates createKernelUpdates( bool const includeState = true ) const + { + if( includeState ) + { + return ElasticIsotropicUpdates( m_bulkModulus, + m_shearModulus, + m_thermalExpansionCoefficient, + m_newStress, + m_oldStress, + m_disableInelasticity ); + } + else // for "no state" updates, pass empty views to avoid transfer of stress data to device + { + return ElasticIsotropicUpdates( m_bulkModulus, + m_shearModulus, + m_thermalExpansionCoefficient, + arrayView3d< real64, solid::STRESS_USD >(), + arrayView3d< real64, solid::STRESS_USD >(), + m_disableInelasticity ); + } + } + + /** + * @brief Construct an update kernel for a derived type. + * @tparam UPDATE_KERNEL The type of update kernel from the derived type. + * @tparam PARAMS The parameter pack to hold the constructor parameters for + * the derived update kernel. + * @param constructorParams The constructor parameter for the derived type. + * @return An @p UPDATE_KERNEL object. + */ + template< typename UPDATE_KERNEL, typename ... PARAMS > + UPDATE_KERNEL createDerivedKernelUpdates( PARAMS && ... constructorParams ) const + { + return UPDATE_KERNEL( std::forward< PARAMS >( constructorParams )..., + m_bulkModulus, + m_shearModulus, + m_thermalExpansionCoefficient, + m_newStress, + m_oldStress, + m_disableInelasticity ); + } + +protected: + + /// Post-process XML data + virtual void postInputInitialization() override; + + /// The default value of the bulk modulus for any new allocations. + real64 m_defaultBulkModulus; + + /// The default value of the shear modulus for any new allocations. + real64 m_defaultShearModulus; + + /// The bulk modulus for each upper level dimension (i.e. cell) of *this + array1d< real64 > m_bulkModulus; + + /// The shear modulus for each upper level dimension (i.e. cell) of *this + array1d< real64 > m_shearModulus; + +}; + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPIC_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 2d94c190a91..37f3c9d59e1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -160,8 +160,8 @@ void SolidMechanicsLagrangianFEM::registerDataOnMesh( Group & meshBodies ) { setConstitutiveNamesCallSuper( subRegion ); - subRegion.registerField< solidMechanics::strain >( getName() ).reference().resizeDimension< 1 >( 6 ); subRegion.registerField< solidMechanics::strain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 ); + subRegion.registerField< solidMechanics::plasticStrain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 ); } ); @@ -456,6 +456,7 @@ void SolidMechanicsLagrangianFEM::initializePostInitialConditionsPreSubGroups() m_targetNodes.insert( m_nonSendOrReceiveNodes.begin(), m_nonSendOrReceiveNodes.end() ); + } ); } ); @@ -946,8 +947,8 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS { string const & solidMaterialName = subRegion.template getReference< string >( viewKeyStruct::solidMaterialNamesString() ); SolidBase & constitutiveRelation = getConstitutiveModel< SolidBase >( subRegion, solidMaterialName ); - constitutiveRelation.saveConvergedState(); + solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >(); solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >(); @@ -967,6 +968,7 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS finiteElement, solidModel, disp, + uhat, strain, plasticStrain ); } ); @@ -974,6 +976,7 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS } ); + constitutiveRelation.saveConvergedState(); } ); } ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 7df079038fd..3c1e3067c30 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -69,6 +69,7 @@ class AverageStrainOverQuadraturePoints : FE_TYPE const & finiteElementSpace, SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, + fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc, fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain): Base( nodeManager, @@ -78,6 +79,7 @@ class AverageStrainOverQuadraturePoints : finiteElementSpace ), m_solidUpdate(solidModel.createKernelUpdates()), m_displacement( displacement ), + m_displacementInc( displacementInc ), m_avgStrain( avgStrain ), m_avgPlasticStrain( avgPlasticStrain ) {} @@ -86,7 +88,8 @@ class AverageStrainOverQuadraturePoints : * @copydoc finiteElement::KernelBase::StackVariables */ struct StackVariables : Base::StackVariables - {real64 uLocal[FE_TYPE::maxSupportPoints][3]; }; + {real64 uLocal[FE_TYPE::maxSupportPoints][3]; + real64 uHatLocal[FE_TYPE::maxSupportPoints][3]; }; /** * @brief Performs the setup phase for the kernel. @@ -105,13 +108,14 @@ class AverageStrainOverQuadraturePoints : for( int i = 0; i < 3; ++i ) { stack.uLocal[a][i] = m_displacement[localNodeIndex][i]; + stack.uHatLocal[a][i] = m_displacementInc[localNodeIndex][i]; } } for( int icomp = 0; icomp < 6; ++icomp ) { m_avgStrain[k][icomp] = 0.0; - m_avgPlasticStrain[k][icomp] = 0.0; + //m_avgPlasticStrain[k][icomp] = 0.0; } } @@ -131,15 +135,18 @@ class AverageStrainOverQuadraturePoints : real64 dNdX[ FE_TYPE::maxSupportPoints ][3]; real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, stack.feStack, dNdX ); real64 strain[6] = {0.0}; + real64 strainInc[6] = {0.0}; FE_TYPE::symmetricGradient( dNdX, stack.uLocal, strain ); + FE_TYPE::symmetricGradient( dNdX, stack.uHatLocal, strainInc ); - real64 elasticStrain[6] = {0.0}; - m_solidUpdate.getElasticStrain(k, q, elasticStrain); + real64 elasticStrainInc[6] = {0.0}; + m_solidUpdate.getElasticStrainInc(k, q, elasticStrainInc); for( int icomp = 0; icomp < 6; ++icomp ) { m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k]; - m_avgPlasticStrain[k][icomp] += detJxW*(strain[icomp] - elasticStrain[icomp])/m_elementVolume[k]; + m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; + //m_avgPlasticStrain[k][icomp] += detJxW*(elasticStrainInc[icomp])/m_elementVolume[k]; } } @@ -177,6 +184,9 @@ class AverageStrainOverQuadraturePoints : /// The displacement solution fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const m_displacement; + /// The displacement increment + fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const m_displacementInc; + /// The average strain fields::solidMechanics::arrayView2dLayoutStrain const m_avgStrain; @@ -220,12 +230,13 @@ class AverageStrainOverQuadraturePointsKernelFactory FE_TYPE const & finiteElementSpace, SOLID_TYPE const & solidModel, fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, + fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc, fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain) { AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE > kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace, - solidModel, displacement, avgStrain, avgPlasticStrain ); + solidModel, displacement, displacementInc, avgStrain, avgPlasticStrain ); AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >::template kernelLaunch< POLICY >( elementSubRegion.size(), kernel ); From 8fda259d00a4d5a5e6d232289be7ba6adf86b31b Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Fri, 11 Oct 2024 15:32:16 -0700 Subject: [PATCH 20/48] remove extraneous file --- "src/coreComponents/constitutive/solid/\177" | 573 ------------------- 1 file changed, 573 deletions(-) delete mode 100644 "src/coreComponents/constitutive/solid/\177" diff --git "a/src/coreComponents/constitutive/solid/\177" "b/src/coreComponents/constitutive/solid/\177" deleted file mode 100644 index af8484f08ae..00000000000 --- "a/src/coreComponents/constitutive/solid/\177" +++ /dev/null @@ -1,573 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2024 Total, S.A - * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2023-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file ElasticIsotropic.hpp - */ - -#ifndef GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPIC_HPP_ -#define GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPIC_HPP_ - -#include "SolidBase.hpp" -#include "PropertyConversions.hpp" -#include "SolidModelDiscretizationOpsIsotropic.hpp" -#include "constitutive/ExponentialRelation.hpp" -#include "LvArray/src/tensorOps.hpp" - -namespace geos -{ - -namespace constitutive -{ - -/** - * @class ElasticIsotropicUpdates - * - * Class to provide elastic isotropic material updates that may be - * called from a kernel function. - */ -class ElasticIsotropicUpdates : public SolidBaseUpdates -{ -public: - /** - * @brief Constructor - * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. - * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. - * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. - * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. - * @param[in] disableInelasticity Flag to disable plasticity for inelastic models - */ - ElasticIsotropicUpdates( arrayView1d< real64 const > const & bulkModulus, - arrayView1d< real64 const > const & shearModulus, - arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView3d< real64, solid::STRESS_USD > const & newStress, - arrayView3d< real64, solid::STRESS_USD > const & oldStress, - const bool & disableInelasticity ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), - m_bulkModulus( bulkModulus ), - m_shearModulus( shearModulus ) - {} - - /// Deleted default constructor - ElasticIsotropicUpdates() = delete; - - /// Default copy constructor - ElasticIsotropicUpdates( ElasticIsotropicUpdates const & ) = default; - - /// Default move constructor - ElasticIsotropicUpdates( ElasticIsotropicUpdates && ) = default; - - /// Deleted copy assignment operator - ElasticIsotropicUpdates & operator=( ElasticIsotropicUpdates const & ) = delete; - - /// Deleted move assignment operator - ElasticIsotropicUpdates & operator=( ElasticIsotropicUpdates && ) = delete; - - /// Use the "isotropic" form of inner product compression - using DiscretizationOps = SolidModelDiscretizationOpsIsotropic; - - /// Use base version of saveConvergedState - using SolidBaseUpdates::saveConvergedState; - - GEOS_HOST_DEVICE - virtual void smallStrainNoStateUpdate_StressOnly( localIndex const k, - localIndex const q, - real64 const ( &totalStrain )[6], - real64 ( &stress )[6] ) const override final; - - GEOS_HOST_DEVICE - virtual void smallStrainNoStateUpdate( localIndex const k, - localIndex const q, - real64 const ( &totalStrain )[6], - real64 ( &stress )[6], - real64 ( &stiffness )[6][6] ) const override final; - - GEOS_HOST_DEVICE - virtual void smallStrainNoStateUpdate( localIndex const k, - localIndex const q, - real64 const ( &totalStrain )[6], - real64 ( &stress )[6], - DiscretizationOps & stiffness ) const final; - - GEOS_HOST_DEVICE - virtual void smallStrainUpdate_StressOnly( localIndex const k, - localIndex const q, - real64 const & timeIncrement, - real64 const ( &strainIncrement )[6], - real64 ( &stress )[6] ) const override; - - GEOS_HOST_DEVICE - void smallStrainUpdate( localIndex const k, - localIndex const q, - real64 const & timeIncrement, - real64 const ( &strainIncrement )[6], - real64 ( &stress )[6], - real64 ( &stiffness )[6][6] ) const; - - GEOS_HOST_DEVICE - virtual void smallStrainUpdate( localIndex const k, - localIndex const q, - real64 const & timeIncrement, - real64 const ( &strainIncrement )[6], - real64 ( &stress )[6], - DiscretizationOps & stiffness ) const; - - GEOS_HOST_DEVICE - virtual void getElasticStiffness( localIndex const k, - localIndex const q, - real64 ( &stiffness )[6][6] ) const override; - - GEOS_HOST_DEVICE - virtual void getElasticStrain( localIndex const k, - localIndex const q, - real64 ( &elasticStrain )[6] ) const override final; - - GEOS_HOST_DEVICE - virtual void getElasticStrainInc( localIndex const k, - localIndex const q, - real64 ( &elasticStrainInc )[6] ) const override final; - - GEOS_HOST_DEVICE - virtual real64 getBulkModulus( localIndex const k ) const override final - { - return m_bulkModulus[k]; - } - - GEOS_HOST_DEVICE - virtual real64 getShearModulus( localIndex const k ) const override final - { - return m_shearModulus[k]; - } - - GEOS_HOST_DEVICE - virtual void viscousStateUpdate( localIndex const k, - localIndex const q, - real64 beta ) const override; - - // TODO: confirm hyper stress/strain measures before activatiing - - /* - GEOS_HOST_DEVICE - virtual void hyperUpdate( localIndex const k, - localIndex const q, - real64 const ( & FminusI )[3][3], - real64 ( & stress )[6] ) const override final; - - GEOS_HOST_DEVICE - virtual void hyperUpdate( localIndex const k, - localIndex const q, - real64 const ( & FminusI )[3][3], - real64 ( & stress )[6], - real64 ( & stiffness )[6][6] ) const override final; - */ - -protected: - - /// A reference to the ArrayView holding the bulk modulus for each element. - arrayView1d< real64 const > const m_bulkModulus; - - /// A reference to the ArrayView holding the shear modulus for each element. - arrayView1d< real64 const > const m_shearModulus; - -}; - - -GEOS_HOST_DEVICE -inline -void ElasticIsotropicUpdates::getElasticStiffness( localIndex const k, - localIndex const q, - real64 ( & stiffness )[6][6] ) const -{ - GEOS_UNUSED_VAR( q ); - real64 const G = m_shearModulus[k]; - real64 const lambda = conversions::bulkModAndShearMod::toFirstLame( m_bulkModulus[k], G ); - - LvArray::tensorOps::fill< 6, 6 >( stiffness, 0 ); - - stiffness[0][0] = lambda + 2*G; - stiffness[0][1] = lambda; - stiffness[0][2] = lambda; - - stiffness[1][0] = lambda; - stiffness[1][1] = lambda + 2*G; - stiffness[1][2] = lambda; - - stiffness[2][0] = lambda; - stiffness[2][1] = lambda; - stiffness[2][2] = lambda + 2*G; - - stiffness[3][3] = G; - stiffness[4][4] = G; - stiffness[5][5] = G; -} - - -GEOS_HOST_DEVICE -inline -void ElasticIsotropicUpdates::getElasticStrain( localIndex const k, - localIndex const q, - real64 ( & elasticStrain)[6] ) const -{ - real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] ); - real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] ); - - elasticStrain[0] = ( m_newStress[k][q][0] - nu*m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E; - elasticStrain[1] = (-nu*m_newStress[k][q][0] + m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E; - elasticStrain[2] = (-nu*m_newStress[k][q][0] - nu*m_newStress[k][q][1] + m_newStress[k][q][2])/E; - - elasticStrain[3] = m_newStress[k][q][3] / m_shearModulus[k]; - elasticStrain[4] = m_newStress[k][q][4] / m_shearModulus[k]; - elasticStrain[5] = m_newStress[k][q][5] / m_shearModulus[k]; -} - - -GEOS_HOST_DEVICE -inline -void ElasticIsotropicUpdates::smallStrainNoStateUpdate_StressOnly( localIndex const k, - localIndex const q, - real64 const ( &totalStrain )[6], - real64 ( & stress )[6] ) const -{ - GEOS_UNUSED_VAR( q ); - - real64 const twoG = 2 * m_shearModulus[k]; - real64 const lambda = conversions::bulkModAndShearMod::toFirstLame( m_bulkModulus[k], m_shearModulus[k] ); - real64 const vol = lambda * ( totalStrain[0] + totalStrain[1] + totalStrain[2] ); - - stress[0] = vol + twoG * totalStrain[0]; - stress[1] = vol + twoG * totalStrain[1]; - stress[2] = vol + twoG * totalStrain[2]; - - stress[3] = m_shearModulus[k] * totalStrain[3]; - stress[4] = m_shearModulus[k] * totalStrain[4]; - stress[5] = m_shearModulus[k] * totalStrain[5]; -} - - -GEOS_HOST_DEVICE -inline -void ElasticIsotropicUpdates::smallStrainNoStateUpdate( localIndex const k, - localIndex const q, - real64 const ( &totalStrain )[6], - real64 ( & stress )[6], - real64 ( & stiffness )[6][6] ) const -{ - smallStrainNoStateUpdate_StressOnly( k, q, totalStrain, stress ); - getElasticStiffness( k, q, stiffness ); -} - - -GEOS_HOST_DEVICE -inline -void ElasticIsotropicUpdates::smallStrainNoStateUpdate( localIndex const k, - localIndex const q, - real64 const ( &totalStrain )[6], - real64 ( & stress )[6], - DiscretizationOps & stiffness ) const -{ - smallStrainNoStateUpdate_StressOnly( k, q, totalStrain, stress ); - stiffness.m_bulkModulus = m_bulkModulus[k]; - stiffness.m_shearModulus = m_shearModulus[k]; -} - - -GEOS_HOST_DEVICE -inline -void ElasticIsotropicUpdates::smallStrainUpdate_StressOnly( localIndex const k, - localIndex const q, - real64 const & timeIncrement, - real64 const ( &strainIncrement )[6], - real64 ( & stress )[6] ) const -{ - GEOS_UNUSED_VAR( timeIncrement ); - smallStrainNoStateUpdate_StressOnly( k, q, strainIncrement, stress ); // stress = incrementalStress - LvArray::tensorOps::add< 6 >( stress, m_oldStress[k][q] ); // stress += m_oldStress - saveStress( k, q, stress ); // m_newStress = stress -} - - -GEOS_HOST_DEVICE -inline -void ElasticIsotropicUpdates::smallStrainUpdate( localIndex const k, - localIndex const q, - real64 const & timeIncrement, - real64 const ( &strainIncrement )[6], - real64 ( & stress )[6], - real64 ( & stiffness )[6][6] ) const -{ - smallStrainUpdate_StressOnly( k, q, timeIncrement, strainIncrement, stress ); - getElasticStiffness( k, q, stiffness ); -} - - -GEOS_HOST_DEVICE -inline -void ElasticIsotropicUpdates::smallStrainUpdate( localIndex const k, - localIndex const q, - real64 const & timeIncrement, - real64 const ( &strainIncrement )[6], - real64 ( & stress )[6], - DiscretizationOps & stiffness ) const -{ - smallStrainUpdate_StressOnly( k, q, timeIncrement, strainIncrement, stress ); - stiffness.m_bulkModulus = m_bulkModulus[k]; - stiffness.m_shearModulus = m_shearModulus[k]; -} - -GEOS_HOST_DEVICE -GEOS_FORCE_INLINE -void ElasticIsotropicUpdates::viscousStateUpdate( localIndex const k, - localIndex const q, - real64 beta ) const -{ - GEOS_UNUSED_VAR( k ); - GEOS_UNUSED_VAR( q ); - GEOS_UNUSED_VAR( beta ); -} - - -// TODO: need to confirm stress / strain measures before activating hyper inferface -/* - GEOS_HOST_DEVICE - inline - void ElasticIsotropicUpdates::hyperUpdate( localIndex const k, - localIndex const q, - real64 const (&FmI)[3][3], - real64 ( & stress )[ 6 ] ) const - { - GEOS_UNUSED_VAR(q); - - real64 const C1 = 0.5 * m_shearModulus[k]; - real64 const D1 = 0.5 * m_bulkModulus[k]; - real64 const detFm1 = FmI[0][0] + FmI[1][1] + FmI[2][2] - - FmI[1][2]*FmI[2][1] + FmI[1][1]*FmI[2][2] - + FmI[0][2]*(-FmI[2][0] - FmI[1][1]*FmI[2][0] + FmI[1][0]*FmI[2][1]) - + FmI[0][1]*(-FmI[1][0] + FmI[1][2]*FmI[2][0] - FmI[1][0]*FmI[2][2]) - + FmI[0][0]*( FmI[1][1] - FmI[1][2]*FmI[2][1] + FmI[2][2] + FmI[1][1]*FmI[2][2]); - - - real64 const p = -2 * D1 * ( detFm1 + 1.0 ) * detFm1; - real64 devB[6] = { 1/3 * (2 * FmI[0][0] * (2 + FmI[0][0]) - FmI[1][1] * (2 + FmI[1][1]) - FmI[2][2] * (2 + FmI[2][2]) - ++ - 2 * FmI[0][1]*FmI[0][1] + 2 * FmI[0][2] * FmI[0][2] - FmI[1][0] * FmI[1][0] - FmI[1][2] * - +FmI[1][2] - - FmI[2][0] * FmI[2][0] - FmI[2][1] * FmI[2][1]), - 1/3 * (-FmI[0][0] * (2 + FmI[0][0]) + 2 * FmI[1][1] * ( 2 + FmI[1][1]) - FmI[2][2] * (2 + - +FmI[2][2]) - - FmI[0][1]*FmI[0][1] - FmI[0][2]*FmI[0][2] + 2 * FmI[1][0]*FmI[1][0] + 2 * - +FmI[1][2]*FmI[1][2] - FmI[2][0]*FmI[2][0] - FmI[2][1]* - FmI[2][1]), - 1/3 *(-FmI[0][0] * (2 + FmI[0][0]) - FmI[1][1] * (2 + FmI[1][1]) + 2 * FmI[2][2] * (2 + FmI[2][2]) - +- - FmI[0][1]*FmI[0][1] - FmI[0][2]*FmI[0][2] - FmI[1][0]*FmI[1][0] - FmI[1][2]*FmI[1][2] + 2 * - +FmI[2][0]*FmI[2][0] + 2 * FmI[2][1]* - FmI[2][1]), - FmI[1][2] + FmI[1][0] * FmI[2][0] + FmI[2][1] + FmI[1][1]*FmI[2][1] + FmI[1][2]*FmI[2][2], - FmI[0][2] + FmI[2][0] + FmI[0][0] * FmI[2][0] + FmI[0][1]*FmI[2][1] + FmI[0][2]*FmI[2][2], - FmI[0][1] + FmI[1][0] + FmI[0][0] * FmI[1][0] + FmI[0][1]*FmI[1][1] + FmI[0][2]*FmI[1][2] - }; - - real64 const C = 2 * C1 / pow( detFm1 + 1, 2.0/3.0 ); - stress[0] = -p + C * devB[0]; - stress[1] = -p + C * devB[1]; - stress[2] = -p + C * devB[2]; - stress[3] = C * devB[3]; - stress[4] = C * devB[4]; - stress[5] = C * devB[5]; - } - */ - - -/** - * @class ElasticIsotropic - * - * Class to provide an elastic isotropic material response. - */ -class ElasticIsotropic : public SolidBase -{ -public: - - /// Alias for ElasticIsotropicUpdates - using KernelWrapper = ElasticIsotropicUpdates; - - /** - * constructor - * @param[in] name name of the instance in the catalog - * @param[in] parent the group which contains this instance - */ - ElasticIsotropic( string const & name, Group * const parent ); - - /** - * Default Destructor - */ - virtual ~ElasticIsotropic() override; - - /** - * @name Static Factory Catalog members and functions - */ - ///@{ - - /// string name to use for this class in the catalog - static constexpr auto m_catalogNameString = "ElasticIsotropic"; - - /** - * @brief Static catalog string - * @return A string that is used to register/lookup this class in the registry - */ - static std::string catalogName() { return m_catalogNameString; } - - /** - * @brief Get catalog name - * @return Name string - */ - virtual string getCatalogName() const override { return catalogName(); } - - ///@} - - /// Keys for data specified in this class. - struct viewKeyStruct : public SolidBase::viewKeyStruct - { - /// string/key for default bulk modulus - static constexpr char const * defaultBulkModulusString() { return "defaultBulkModulus"; } - - /// string/key for default poisson ratio - static constexpr char const * defaultPoissonRatioString() { return "defaultPoissonRatio"; } - - /// string/key for default shear modulus - static constexpr char const * defaultShearModulusString() { return "defaultShearModulus"; } - - /// string/key for default Young's modulus - static constexpr char const * defaultYoungModulusString() { return "defaultYoungModulus"; } - - /// string/key for bulk modulus - static constexpr char const * bulkModulusString() { return "bulkModulus"; } - - /// string/key for shear modulus - static constexpr char const * shearModulusString() { return "shearModulus"; } - - }; - - /** - * @brief Accessor for bulk modulus - * @return A const reference to arrayView1d containing the bulk - * modulus (at every element). - */ - arrayView1d< real64 > const bulkModulus() { return m_bulkModulus; } - - /** - * @brief Const accessor for bulk modulus - * @return A const reference to arrayView1d containing the bulk - * modulus (at every element). - */ - arrayView1d< real64 const > const bulkModulus() const { return m_bulkModulus; } - - /** - * @brief Accessor for shear modulus - * @return A const reference to arrayView1d containing the shear - * modulus (at every element). - */ - arrayView1d< real64 > const shearModulus() { return m_shearModulus; } - - /** - * @brief Const accessor for shear modulus - * @return A const reference to arrayView1d containing the - * shear modulus (at every element). - */ - arrayView1d< real64 const > const shearModulus() const { return m_shearModulus; } - - GEOS_HOST_DEVICE - virtual arrayView1d< real64 const > getBulkModulus() const override final - { - return m_bulkModulus; - } - GEOS_HOST_DEVICE - virtual arrayView1d< real64 const > getShearModulus() const override final - { - return m_shearModulus; - } - - /** - * @brief Create a instantiation of the ElasticIsotropicUpdate class - * that refers to the data in this. - * @param includeState Flag whether to pass state arrays that may not be needed for "no-state" updates - * @return An instantiation of ElasticIsotropicUpdate. - */ - ElasticIsotropicUpdates createKernelUpdates( bool const includeState = true ) const - { - if( includeState ) - { - return ElasticIsotropicUpdates( m_bulkModulus, - m_shearModulus, - m_thermalExpansionCoefficient, - m_newStress, - m_oldStress, - m_disableInelasticity ); - } - else // for "no state" updates, pass empty views to avoid transfer of stress data to device - { - return ElasticIsotropicUpdates( m_bulkModulus, - m_shearModulus, - m_thermalExpansionCoefficient, - arrayView3d< real64, solid::STRESS_USD >(), - arrayView3d< real64, solid::STRESS_USD >(), - m_disableInelasticity ); - } - } - - /** - * @brief Construct an update kernel for a derived type. - * @tparam UPDATE_KERNEL The type of update kernel from the derived type. - * @tparam PARAMS The parameter pack to hold the constructor parameters for - * the derived update kernel. - * @param constructorParams The constructor parameter for the derived type. - * @return An @p UPDATE_KERNEL object. - */ - template< typename UPDATE_KERNEL, typename ... PARAMS > - UPDATE_KERNEL createDerivedKernelUpdates( PARAMS && ... constructorParams ) const - { - return UPDATE_KERNEL( std::forward< PARAMS >( constructorParams )..., - m_bulkModulus, - m_shearModulus, - m_thermalExpansionCoefficient, - m_newStress, - m_oldStress, - m_disableInelasticity ); - } - -protected: - - /// Post-process XML data - virtual void postInputInitialization() override; - - /// The default value of the bulk modulus for any new allocations. - real64 m_defaultBulkModulus; - - /// The default value of the shear modulus for any new allocations. - real64 m_defaultShearModulus; - - /// The bulk modulus for each upper level dimension (i.e. cell) of *this - array1d< real64 > m_bulkModulus; - - /// The shear modulus for each upper level dimension (i.e. cell) of *this - array1d< real64 > m_shearModulus; - -}; - -} /* namespace constitutive */ - -} /* namespace geos */ - -#endif /* GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPIC_HPP_ */ From f5f1fa790c0e2fc49418485b7b0346cfbe461f41 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Tue, 15 Oct 2024 08:48:53 -0700 Subject: [PATCH 21/48] initial proposal for handling prescribed displacements and sign differences --- .../solidMechanics/kernels/StrainHelper.hpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 3c1e3067c30..36e7badd5c0 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -145,7 +145,14 @@ class AverageStrainOverQuadraturePoints : for( int icomp = 0; icomp < 6; ++icomp ) { m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k]; - m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; + + // This is maybe bad on gpu + // How to hamdle magnitudes? + if ((std::abs(strainInc[icomp]) > std::abs(elasticStrainInc[icomp]) ) && (std::abs(strainInc[icomp]) > 1.0e-8 )) + { + m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; + } + //m_avgStrain[k][icomp] += detJxW*(strainInc[icomp])/m_elementVolume[k]; //m_avgPlasticStrain[k][icomp] += detJxW*(elasticStrainInc[icomp])/m_elementVolume[k]; } } From 48db73defc4e2d98868428b87382b143967e13f5 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 4 Dec 2024 09:54:55 -0800 Subject: [PATCH 22/48] simple method for handling BCs, may want to put in a disclaimer --- .../solidMechanics/kernels/StrainHelper.hpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 36e7badd5c0..a75a2eff9d7 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -146,12 +146,17 @@ class AverageStrainOverQuadraturePoints : { m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k]; - // This is maybe bad on gpu - // How to hamdle magnitudes? - if ((std::abs(strainInc[icomp]) > std::abs(elasticStrainInc[icomp]) ) && (std::abs(strainInc[icomp]) > 1.0e-8 )) + if (std::abs(strainInc[icomp]) > 1.0e-8) { m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; } + + // This is maybe bad on gpu + // How to hamdle magnitudes? + //if ((std::abs(strainInc[icomp]) > std::abs(elasticStrainInc[icomp]) ) && (std::abs(strainInc[icomp]) > 1.0e-8 )) + //{ + // m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; + //} //m_avgStrain[k][icomp] += detJxW*(strainInc[icomp])/m_elementVolume[k]; //m_avgPlasticStrain[k][icomp] += detJxW*(elasticStrainInc[icomp])/m_elementVolume[k]; } From 583e7dea8e8c083fac21a71abb4f6b070cbf734d Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 5 Dec 2024 11:25:20 -0800 Subject: [PATCH 23/48] Update SolidMechanicsStateReset.cpp --- .../solidMechanics/SolidMechanicsStateReset.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index 6c35edca982..b747cae2658 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp @@ -98,6 +98,16 @@ bool SolidMechanicsStateReset::execute( real64 const time_n, } nodeManager.getField< solidMechanics::totalDisplacement >().zero(); nodeManager.getField< solidMechanics::incrementalDisplacement >().zero(); + + ElementRegionManager & elemManager = meshLevel.getElemManager(); + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + subRegion.getField< solidMechanics::strain >().zero(); + subRegion.getField< solidMechanics::plasticStrain >().zero(); + } ); + } // Option 2: enable / disable inelastic behavior From e9bd9536202c37887816e917113e53f387ac2755 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 5 Dec 2024 11:28:00 -0800 Subject: [PATCH 24/48] Update StrainHelper.hpp --- .../solidMechanics/kernels/StrainHelper.hpp | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index a75a2eff9d7..3bd2a096744 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -146,19 +146,12 @@ class AverageStrainOverQuadraturePoints : { m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k]; + // This is a hack to handle boundary conditions such as those seen in plane-strain wellbore problems + // Essentially, if bcs are constraining the strain (and thus total displacement), we do not accumulate any plastic strain (regardless of stresses in material law) if (std::abs(strainInc[icomp]) > 1.0e-8) { m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; } - - // This is maybe bad on gpu - // How to hamdle magnitudes? - //if ((std::abs(strainInc[icomp]) > std::abs(elasticStrainInc[icomp]) ) && (std::abs(strainInc[icomp]) > 1.0e-8 )) - //{ - // m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; - //} - //m_avgStrain[k][icomp] += detJxW*(strainInc[icomp])/m_elementVolume[k]; - //m_avgPlasticStrain[k][icomp] += detJxW*(elasticStrainInc[icomp])/m_elementVolume[k]; } } From 25ddb6a53c9483ba7b97b8d91ff72e26dad59e70 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 5 Dec 2024 12:00:05 -0800 Subject: [PATCH 25/48] Update SolidMechanicsStateReset.cpp --- .../physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index b747cae2658..c83a9f159ac 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp @@ -99,7 +99,7 @@ bool SolidMechanicsStateReset::execute( real64 const time_n, nodeManager.getField< solidMechanics::totalDisplacement >().zero(); nodeManager.getField< solidMechanics::incrementalDisplacement >().zero(); - ElementRegionManager & elemManager = meshLevel.getElemManager(); + ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, ElementSubRegionBase & subRegion ) From ccc1b0d9270d8611832e0330cbf603b23f1910f6 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 5 Dec 2024 13:31:03 -0800 Subject: [PATCH 26/48] uncrustify --- .../constitutive/solid/ElasticIsotropic.hpp | 2 +- .../ElasticIsotropicPressureDependent.hpp | 12 +++---- .../constitutive/solid/SolidBase.hpp | 4 +-- .../SolidMechanicsLagrangianFEM.cpp | 34 +++++++++---------- .../SolidMechanicsStateReset.cpp | 14 ++++---- .../solidMechanics/kernels/StrainHelper.hpp | 13 +++---- 6 files changed, 40 insertions(+), 39 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index 12975c80ae8..16c306ba1c9 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -139,7 +139,7 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates virtual void getElasticStrainInc( localIndex const k, localIndex const q, real64 ( &elasticStrainInc )[6] ) const override final; - + GEOS_HOST_DEVICE virtual real64 getBulkModulus( localIndex const k ) const override final { diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 671a002421a..978f3c5bb0b 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -117,9 +117,9 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates GEOS_HOST_DEVICE virtual void getElasticStrainInc( localIndex const k, - localIndex const q, - real64 ( &elasticStrainInc )[6] ) const override final; - + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + GEOS_HOST_DEVICE virtual void viscousStateUpdate( localIndex const k, localIndex const q, @@ -231,8 +231,8 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrain( localIndex cons GEOS_HOST_DEVICE inline void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex const k, - localIndex const q, - real64 ( & elasticStrainInc)[6] ) const + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const { real64 const mu = m_shearModulus[k]; real64 const p0 = m_refPressure; @@ -282,7 +282,7 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c deviator, oldStrain ); - for (localIndex i = 0; i<6; ++i) + for( localIndex i = 0; i<6; ++i ) { elasticStrainInc[i] -= oldStrain[i]; } diff --git a/src/coreComponents/constitutive/solid/SolidBase.hpp b/src/coreComponents/constitutive/solid/SolidBase.hpp index f8cd559f645..17537bc2fbe 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.hpp +++ b/src/coreComponents/constitutive/solid/SolidBase.hpp @@ -346,8 +346,8 @@ class SolidBaseUpdates */ GEOS_HOST_DEVICE virtual void getElasticStrainInc( localIndex const k, - localIndex const q, - real64 ( & elasticStrainInc )[6] ) const + localIndex const q, + real64 ( & elasticStrainInc )[6] ) const { GEOS_UNUSED_VAR( k ); GEOS_UNUSED_VAR( q ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 345ae737b17..a4d2b4e8945 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -949,30 +949,30 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS string const & solidMaterialName = subRegion.template getReference< string >( viewKeyStruct::solidMaterialNamesString() ); SolidBase & constitutiveRelation = getConstitutiveModel< SolidBase >( subRegion, solidMaterialName ); - + solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >(); solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >(); - constitutive::ConstitutivePassThru< SolidBase >::execute(constitutiveRelation, [&] (auto & solidModel) + constitutive::ConstitutivePassThru< SolidBase >::execute( constitutiveRelation, [&] ( auto & solidModel ) { using SOLID_TYPE = TYPEOFREF( solidModel ); - finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); - finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) - { - using FE_TYPE = decltype( finiteElement ); - AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager, - mesh.getEdgeManager(), - mesh.getFaceManager(), - subRegion, - finiteElement, - solidModel, - disp, - uhat, - strain, - plasticStrain ); - } ); + finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName()); + finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement ) + { + using FE_TYPE = decltype( finiteElement ); + AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager, + mesh.getEdgeManager(), + mesh.getFaceManager(), + subRegion, + finiteElement, + solidModel, + disp, + uhat, + strain, + plasticStrain ); + } ); } ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index c83a9f159ac..c7ff94d5e58 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp @@ -101,13 +101,13 @@ bool SolidMechanicsStateReset::execute( real64 const time_n, ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, - [&]( localIndex const, - ElementSubRegionBase & subRegion ) - { - subRegion.getField< solidMechanics::strain >().zero(); - subRegion.getField< solidMechanics::plasticStrain >().zero(); - } ); - + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + subRegion.getField< solidMechanics::strain >().zero(); + subRegion.getField< solidMechanics::plasticStrain >().zero(); + } ); + } // Option 2: enable / disable inelastic behavior diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 3bd2a096744..bfcd31cdc5f 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -71,13 +71,13 @@ class AverageStrainOverQuadraturePoints : fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc, fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, - fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain): + fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain ): Base( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace ), - m_solidUpdate(solidModel.createKernelUpdates()), + m_solidUpdate( solidModel.createKernelUpdates()), m_displacement( displacement ), m_displacementInc( displacementInc ), m_avgStrain( avgStrain ), @@ -140,15 +140,16 @@ class AverageStrainOverQuadraturePoints : FE_TYPE::symmetricGradient( dNdX, stack.uHatLocal, strainInc ); real64 elasticStrainInc[6] = {0.0}; - m_solidUpdate.getElasticStrainInc(k, q, elasticStrainInc); + m_solidUpdate.getElasticStrainInc( k, q, elasticStrainInc ); for( int icomp = 0; icomp < 6; ++icomp ) { m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k]; // This is a hack to handle boundary conditions such as those seen in plane-strain wellbore problems - // Essentially, if bcs are constraining the strain (and thus total displacement), we do not accumulate any plastic strain (regardless of stresses in material law) - if (std::abs(strainInc[icomp]) > 1.0e-8) + // Essentially, if bcs are constraining the strain (and thus total displacement), we do not accumulate any plastic strain (regardless + // of stresses in material law) + if( std::abs( strainInc[icomp] ) > 1.0e-8 ) { m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; } @@ -237,7 +238,7 @@ class AverageStrainOverQuadraturePointsKernelFactory fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement, fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc, fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, - fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain) + fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain ) { AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE > kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace, From 179e4fd2a1b5e1748c6cb840497819dcb1921f82 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 11 Dec 2024 11:23:24 -0800 Subject: [PATCH 27/48] Apply suggestions from code review Co-authored-by: Matteo Cusini <49037133+CusiniM@users.noreply.github.com> --- .../solid/ElasticIsotropicPressureDependent.hpp | 12 ++++++------ src/coreComponents/constitutive/solid/SolidBase.hpp | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 978f3c5bb0b..5351f4b2a7e 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -238,8 +238,8 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c real64 const p0 = m_refPressure; real64 const eps_v0 = m_refStrainVol; real64 const Cr = m_recompressionIndex[k]; - real64 deviator[6]; - real64 stress[6]; + real64 deviator[6]{}; + real64 stress[6]{}; real64 P; real64 Q; real64 elasticStrainVol; @@ -255,15 +255,15 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c Q, deviator ); - elasticStrainVol = std::log( P/p0 ) * Cr * (-1.0) + eps_v0; - elasticStrainDev = Q/3./mu; + real64 const elasticStrainVol = LvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0; + real64 const elasticStrainDev = Q/3./mu; twoInvariant::strainRecomposition( elasticStrainVol, elasticStrainDev, deviator, elasticStrainInc ); - real64 oldStrain[6]; + real64 oldStrain[6]{}; for( localIndex i=0; i<6; ++i ) { stress[i] = m_oldStress[k][q][i]; @@ -275,7 +275,7 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c deviator ); elasticStrainVol = std::log( P/p0 ) * Cr * (-1.0) + eps_v0; - elasticStrainDev = Q/3./mu; + real64 const elasticStrainDev = Q/3./mu; twoInvariant::strainRecomposition( elasticStrainVol, elasticStrainDev, diff --git a/src/coreComponents/constitutive/solid/SolidBase.hpp b/src/coreComponents/constitutive/solid/SolidBase.hpp index 17537bc2fbe..a1b651f5d94 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.hpp +++ b/src/coreComponents/constitutive/solid/SolidBase.hpp @@ -352,7 +352,7 @@ class SolidBaseUpdates GEOS_UNUSED_VAR( k ); GEOS_UNUSED_VAR( q ); GEOS_UNUSED_VAR( elasticStrainInc ); - GEOS_ERROR( "getElasticStrainInc() not implemented for this model (called when computing plastic strains)" ); + GEOS_ERROR( "getElasticStrainInc() of SolidBase was called." ); } /** From 834bba6dd7f2ab087325e14eae720f49f39268a1 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 11 Dec 2024 11:27:58 -0800 Subject: [PATCH 28/48] Apply suggestions from code review Co-authored-by: Matteo Cusini <49037133+CusiniM@users.noreply.github.com> --- .../constitutive/solid/ElasticIsotropicPressureDependent.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 5351f4b2a7e..74e4ddc92d8 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -242,8 +242,6 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c real64 stress[6]{}; real64 P; real64 Q; - real64 elasticStrainVol; - real64 elasticStrainDev; for( localIndex i=0; i<6; ++i ) { From 3b9372ecddbf24a525c33069adc8477373cf3407 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 11 Dec 2024 11:37:03 -0800 Subject: [PATCH 29/48] Update ElasticIsotropicPressureDependent.hpp --- .../solid/ElasticIsotropicPressureDependent.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 74e4ddc92d8..c007d03e61f 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -253,8 +253,8 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c Q, deviator ); - real64 const elasticStrainVol = LvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0; - real64 const elasticStrainDev = Q/3./mu; + real64 elasticStrainVol = LvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0; + real64 elasticStrainDev = Q/3./mu; twoInvariant::strainRecomposition( elasticStrainVol, elasticStrainDev, @@ -273,7 +273,7 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c deviator ); elasticStrainVol = std::log( P/p0 ) * Cr * (-1.0) + eps_v0; - real64 const elasticStrainDev = Q/3./mu; + elasticStrainDev = Q/3./mu; twoInvariant::strainRecomposition( elasticStrainVol, elasticStrainDev, From 136fe651402d963de8e318d219a71fb1e06aa368 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 11 Dec 2024 13:17:14 -0800 Subject: [PATCH 30/48] Update ElasticIsotropicPressureDependent.hpp --- .../solid/ElasticIsotropicPressureDependent.hpp | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index c007d03e61f..9e8ccba4a97 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -243,10 +243,7 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c real64 P; real64 Q; - for( localIndex i=0; i<6; ++i ) - { - stress[i] = m_newStress[k][q][i]; - } + LvArray::tensorOps::copy<6>( stress, m_newStress[k][q] ); twoInvariant::stressDecomposition( stress, P, @@ -262,17 +259,15 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c elasticStrainInc ); real64 oldStrain[6]{}; - for( localIndex i=0; i<6; ++i ) - { - stress[i] = m_oldStress[k][q][i]; - } + + LvArray::tensorOps::copy<6>( stress, m_oldStress[k][q] ); twoInvariant::stressDecomposition( stress, P, Q, deviator ); - elasticStrainVol = std::log( P/p0 ) * Cr * (-1.0) + eps_v0; + elasticStrainVol = lvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0; elasticStrainDev = Q/3./mu; twoInvariant::strainRecomposition( elasticStrainVol, From 246dbc83466b326d654c2d0d07c56a2598e77b91 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 11 Dec 2024 13:23:13 -0800 Subject: [PATCH 31/48] Update ElasticIsotropicPressureDependent.hpp --- .../constitutive/solid/ElasticIsotropicPressureDependent.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 84e4bd7f5ea..0dee4c5ad97 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -267,7 +267,7 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c Q, deviator ); - elasticStrainVol = lvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0; + elasticStrainVol = LvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0; elasticStrainDev = Q/3./mu; twoInvariant::strainRecomposition( elasticStrainVol, From bb4dfc8cad6d1a6fbb4c6c8d5122cad562c0fbf7 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 11 Dec 2024 14:58:02 -0800 Subject: [PATCH 32/48] uncrustify --- .../constitutive/solid/ElasticIsotropicPressureDependent.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 0dee4c5ad97..01f36bf937f 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -243,7 +243,7 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c real64 P; real64 Q; - LvArray::tensorOps::copy<6>( stress, m_newStress[k][q] ); + LvArray::tensorOps::copy< 6 >( stress, m_newStress[k][q] ); twoInvariant::stressDecomposition( stress, P, @@ -260,7 +260,7 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex c real64 oldStrain[6]{}; - LvArray::tensorOps::copy<6>( stress, m_oldStress[k][q] ); + LvArray::tensorOps::copy< 6 >( stress, m_oldStress[k][q] ); twoInvariant::stressDecomposition( stress, P, From 4e6c411f61b30fea161a7ef9eb27b4e9126661b6 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 11 Dec 2024 15:44:54 -0800 Subject: [PATCH 33/48] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 31ce96119a9..c6feedf2165 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3450-9221-37d940c + baseline: integratedTests/baseline_integratedTests-pr3384-9299-246dbc8 allow_fail: all: '' streak: '' From ce59c5411dbb00611c5e3e7c53640a32834d8129 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Wed, 11 Dec 2024 15:45:44 -0800 Subject: [PATCH 34/48] Update BASELINE_NOTES.md --- BASELINE_NOTES.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 3b097223255..cd5d47c179f 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3384 (2024-12-11) +===================== +Added plastic strain output. + PR #3450 (2024-12-08) ===================== Added test for explicit runge kutta sprinslider. From beae829a76c90f12fc1f0147d54d489b317f578a Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Fri, 13 Dec 2024 14:24:48 -0800 Subject: [PATCH 35/48] add elastic strain functions to ElasticOrthotropic --- .../constitutive/solid/ElasticOrthotropic.hpp | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index 4445b88cc12..d048fec9efd 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -147,6 +147,16 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates real64 ( &stress )[6], DiscretizationOps & stiffness ) const final; + GEOS_HOST_DEVICE + virtual void getElasticStrain( localIndex const k, + localIndex const q, + real64 ( &elasticStrain )[6] ) const override final; + + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + // miscellaneous getters GEOS_HOST_DEVICE @@ -219,6 +229,44 @@ void ElasticOrthotropicUpdates::getElasticStiffness( localIndex const k, stiffness[5][5] = m_c66[k]; } + +GEOS_HOST_DEVICE +inline +void ElasticOrthotropicUpdates::getElasticStrain( localIndex const k, + localIndex const q, + real64 ( & elasticStrain)[6] ) const +{ + + real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]); + + elasticStrain[0] = ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*m_newStress[k][q][0] + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*m_newStress[k][q][1] + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*m_newStress[k][q][2] ) / detC; + elasticStrain[1] = ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*m_newStress[k][q][0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][1] + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*m_newStress[k][q][2] ) / detC; + elasticStrain[2] = ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*m_newStress[k][q][0] + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*m_newStress[k][q][1] + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*m_newStress[k][q][2] ) / detC; + + elasticStrain[3] = m_newStress[k][q][3] / m_c44[k]; + elasticStrain[4] = m_newStress[k][q][4] / m_c55[k]; + elasticStrain[5] = m_newStress[k][q][5] / m_c66[k]; +} + +GEOS_HOST_DEVICE +inline +void ElasticOrthotropicUpdates::getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const +{ + + real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]); + + elasticStrain[0] = ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrain[1] = ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrain[2] = ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + + elasticStrain[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_c44[k]; + elasticStrain[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_c55[k]; + elasticStrain[5] = (m_newStress[k][q][5] - m_oldStress[k][q][5]) / m_c66[k]; +} + + inline GEOS_HOST_DEVICE void ElasticOrthotropicUpdates::smallStrainNoStateUpdate_StressOnly( localIndex const k, From a5164dcb4b571919ff784ad513661a6b92d45a87 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Fri, 13 Dec 2024 14:36:52 -0800 Subject: [PATCH 36/48] Update ElasticOrthotropic.hpp --- .../constitutive/solid/ElasticOrthotropic.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index d048fec9efd..446ac06caee 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -257,13 +257,13 @@ void ElasticOrthotropicUpdates::getElasticStrainInc( localIndex const k, real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]); - elasticStrain[0] = ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrain[1] = ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrain[2] = ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[0] = ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[1] = ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[2] = ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrain[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_c44[k]; - elasticStrain[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_c55[k]; - elasticStrain[5] = (m_newStress[k][q][5] - m_oldStress[k][q][5]) / m_c66[k]; + elasticStrainInc[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_c44[k]; + elasticStrainInc[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_c55[k]; + elasticStrainInc[5] = (m_newStress[k][q][5] - m_oldStress[k][q][5]) / m_c66[k]; } From d4039917e284593cde2cb25f2b8cb186cf8233d7 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Fri, 13 Dec 2024 14:46:40 -0800 Subject: [PATCH 37/48] Add elastic strain functions to ElasticTransverseIsotropic --- .../solid/ElasticTransverseIsotropic.hpp | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index f2c86d4439a..26324dccfad 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -144,6 +144,16 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates GEOS_HOST_DEVICE virtual void getElasticStiffness( localIndex const k, localIndex const q, real64 ( &stiffness )[6][6] ) const override final; + GEOS_HOST_DEVICE + virtual void getElasticStrain( localIndex const k, + localIndex const q, + real64 ( &elasticStrain )[6] ) const override final; + + GEOS_HOST_DEVICE + virtual void getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( &elasticStrainInc )[6] ) const override final; + /** * @brief Getter for apparent shear modulus. * @return reference to shear modulus that will be used for computing stabilization scalling parameter. @@ -200,6 +210,43 @@ void ElasticTransverseIsotropicUpdates::getElasticStiffness( localIndex const k, stiffness[5][5] = m_c66[k]; } +GEOS_HOST_DEVICE +inline +void ElasticTransverseIsotropicUpdates::getElasticStrain( localIndex const k, + localIndex const q, + real64 ( & elasticStrain)[6] ) const +{ + real64 const c12 = ( m_c11[k] - 2.0 * m_c66[k] ); + real64 const detC = m_c11[k]*(m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k]) - c12*(c12*m_c33[k] - m_c13[k]*m_c13[k]) + m_c13[k]*(c12*m_c13[k] - m_c11[k]*m_c13[k]); + + elasticStrain[0] = ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][0] + (m_c13[k]*m_c13[k] - c12*m_c33[k])*m_newStress[k][q][1] + (c12*m_c13[k] - m_c13[k]*m_c11[k])*m_newStress[k][q][2] ) / detC; + elasticStrain[1] = ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*m_newStress[k][q][0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][1] + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*m_newStress[k][q][2] ) / detC; + elasticStrain[2] = ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*m_newStress[k][q][0] + (c12*m_c13[k] - m_c11[k]*m_c13[k])*m_newStress[k][q][1] + (m_c11[k]*m_c11[k] - c12*c12)*m_newStress[k][q][2] ) / detC; + + elasticStrain[3] = m_newStress[k][q][3] / m_c44[k]; + elasticStrain[4] = m_newStress[k][q][4] / m_c44[k]; + elasticStrain[5] = m_newStress[k][q][5] / m_c66[k]; +} + +GEOS_HOST_DEVICE +inline +void ElasticTransverseIsotropicUpdates::getElasticStrainInc( localIndex const k, + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const +{ + + real64 const c12 = ( m_c11[k] - 2.0 * m_c66[k] ); + real64 const detC = m_c11[k]*(m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k]) - c12*(c12*m_c33[k] - m_c13[k]*m_c13[k]) + m_c13[k]*(c12*m_c13[k] - m_c11[k]*m_c13[k]); + + elasticStrainInc[0] = ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c13[k] - c12*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (c12*m_c13[k] - m_c13[k]*m_c11[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[1] = ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[2] = ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (c12*m_c13[k] - m_c11[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c11[k]*m_c11[k] - c12*c12)*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + + elasticStrainInc[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_c44[k]; + elasticStrainInc[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_c44[k]; + elasticStrainInc[5] = (m_newStress[k][q][5] - m_oldStress[k][q][5]) / m_c66[k]; +} + inline GEOS_HOST_DEVICE void ElasticTransverseIsotropicUpdates::smallStrainNoStateUpdate_StressOnly( localIndex const k, From b04133dbda190bfb5c0c528430adc3b019931da2 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Mon, 16 Dec 2024 09:42:50 -0800 Subject: [PATCH 38/48] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index c6feedf2165..c6b7d19b084 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3384-9299-246dbc8 + baseline: integratedTests/baseline_integratedTests-pr3384-9342-d403991 allow_fail: all: '' streak: '' From d8a1680b3e95d81bcc10362a6bc3b1a59b86ebde Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Mon, 16 Dec 2024 09:51:44 -0800 Subject: [PATCH 39/48] uncrustify --- .../constitutive/solid/ElasticOrthotropic.hpp | 32 +++++++++++++------ .../solid/ElasticTransverseIsotropic.hpp | 26 +++++++++------ 2 files changed, 39 insertions(+), 19 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index 446ac06caee..0958b3b89cf 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -233,15 +233,21 @@ void ElasticOrthotropicUpdates::getElasticStiffness( localIndex const k, GEOS_HOST_DEVICE inline void ElasticOrthotropicUpdates::getElasticStrain( localIndex const k, - localIndex const q, - real64 ( & elasticStrain)[6] ) const + localIndex const q, + real64 ( & elasticStrain)[6] ) const { real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]); - elasticStrain[0] = ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*m_newStress[k][q][0] + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*m_newStress[k][q][1] + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*m_newStress[k][q][2] ) / detC; - elasticStrain[1] = ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*m_newStress[k][q][0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][1] + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*m_newStress[k][q][2] ) / detC; - elasticStrain[2] = ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*m_newStress[k][q][0] + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*m_newStress[k][q][1] + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*m_newStress[k][q][2] ) / detC; + elasticStrain[0] = + ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*m_newStress[k][q][0] + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*m_newStress[k][q][1] + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*m_newStress[k][q][2] ) / + detC; + elasticStrain[1] = + ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*m_newStress[k][q][0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][1] + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*m_newStress[k][q][2] ) / + detC; + elasticStrain[2] = + ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*m_newStress[k][q][0] + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*m_newStress[k][q][1] + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*m_newStress[k][q][2] ) / + detC; elasticStrain[3] = m_newStress[k][q][3] / m_c44[k]; elasticStrain[4] = m_newStress[k][q][4] / m_c55[k]; @@ -251,15 +257,21 @@ void ElasticOrthotropicUpdates::getElasticStrain( localIndex const k, GEOS_HOST_DEVICE inline void ElasticOrthotropicUpdates::getElasticStrainInc( localIndex const k, - localIndex const q, - real64 ( & elasticStrainInc)[6] ) const + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const { real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]); - elasticStrainInc[0] = ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrainInc[1] = ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrainInc[2] = ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[0] = + ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[1] = + ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[2] = + ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; elasticStrainInc[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_c44[k]; elasticStrainInc[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_c55[k]; diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index 26324dccfad..759d4cbcf0d 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -213,14 +213,16 @@ void ElasticTransverseIsotropicUpdates::getElasticStiffness( localIndex const k, GEOS_HOST_DEVICE inline void ElasticTransverseIsotropicUpdates::getElasticStrain( localIndex const k, - localIndex const q, - real64 ( & elasticStrain)[6] ) const + localIndex const q, + real64 ( & elasticStrain)[6] ) const { real64 const c12 = ( m_c11[k] - 2.0 * m_c66[k] ); real64 const detC = m_c11[k]*(m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k]) - c12*(c12*m_c33[k] - m_c13[k]*m_c13[k]) + m_c13[k]*(c12*m_c13[k] - m_c11[k]*m_c13[k]); - elasticStrain[0] = ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][0] + (m_c13[k]*m_c13[k] - c12*m_c33[k])*m_newStress[k][q][1] + (c12*m_c13[k] - m_c13[k]*m_c11[k])*m_newStress[k][q][2] ) / detC; - elasticStrain[1] = ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*m_newStress[k][q][0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][1] + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*m_newStress[k][q][2] ) / detC; + elasticStrain[0] = + ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][0] + (m_c13[k]*m_c13[k] - c12*m_c33[k])*m_newStress[k][q][1] + (c12*m_c13[k] - m_c13[k]*m_c11[k])*m_newStress[k][q][2] ) / detC; + elasticStrain[1] = + ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*m_newStress[k][q][0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][1] + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*m_newStress[k][q][2] ) / detC; elasticStrain[2] = ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*m_newStress[k][q][0] + (c12*m_c13[k] - m_c11[k]*m_c13[k])*m_newStress[k][q][1] + (m_c11[k]*m_c11[k] - c12*c12)*m_newStress[k][q][2] ) / detC; elasticStrain[3] = m_newStress[k][q][3] / m_c44[k]; @@ -231,16 +233,22 @@ void ElasticTransverseIsotropicUpdates::getElasticStrain( localIndex const k, GEOS_HOST_DEVICE inline void ElasticTransverseIsotropicUpdates::getElasticStrainInc( localIndex const k, - localIndex const q, - real64 ( & elasticStrainInc)[6] ) const + localIndex const q, + real64 ( & elasticStrainInc)[6] ) const { real64 const c12 = ( m_c11[k] - 2.0 * m_c66[k] ); real64 const detC = m_c11[k]*(m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k]) - c12*(c12*m_c33[k] - m_c13[k]*m_c13[k]) + m_c13[k]*(c12*m_c13[k] - m_c11[k]*m_c13[k]); - elasticStrainInc[0] = ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c13[k] - c12*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (c12*m_c13[k] - m_c13[k]*m_c11[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrainInc[1] = ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrainInc[2] = ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (c12*m_c13[k] - m_c11[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_c11[k]*m_c11[k] - c12*c12)*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[0] = + ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c13[k] - c12*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + + (c12*m_c13[k] - m_c13[k]*m_c11[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[1] = + ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; + elasticStrainInc[2] = + ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (c12*m_c13[k] - m_c11[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + + (m_c11[k]*m_c11[k] - c12*c12)*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; elasticStrainInc[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_c44[k]; elasticStrainInc[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_c44[k]; From 01cce1e59aa429c3e4072eea5cb21fa85b424a77 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Mon, 16 Dec 2024 12:12:10 -0800 Subject: [PATCH 40/48] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index c6b7d19b084..6dfbf981b0c 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3384-9342-d403991 + baseline: integratedTests/baseline_integratedTests-pr3479-9362-cffefcc allow_fail: all: '' streak: '' From caef37b933e8528d9a28be918db3b0821db3615d Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Mon, 16 Dec 2024 14:22:27 -0800 Subject: [PATCH 41/48] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 6dfbf981b0c..63e604a8ca5 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3479-9362-cffefcc + baseline: integratedTests/baseline_integratedTests-pr3384-9386-cfe31b2 allow_fail: all: '' streak: '' From 69ade8dfbdea2ca8849797d238b8c8fde167cb09 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Mon, 16 Dec 2024 15:48:10 -0800 Subject: [PATCH 42/48] restructure elasticIsotropic --- .../constitutive/solid/ElasticIsotropic.hpp | 51 +++++++++++++------ 1 file changed, 35 insertions(+), 16 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index 7edc5e7d5f0..268cf961a58 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -176,6 +176,13 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates protected: + GEOS_HOST_DEVICE + virtual void computeElasticStrain( localIndex const k, + localIndex const q, + real64 const ( &stress )[6], + real64 ( &elasticStrainInc )[6] ) const; + + /// A reference to the ArrayView holding the bulk modulus for each element. arrayView1d< real64 const > const m_bulkModulus; @@ -215,22 +222,37 @@ void ElasticIsotropicUpdates::getElasticStiffness( localIndex const k, } +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::computeElasticStrain( localIndex const k, + localIndex const q, + real64 const ( &stress )[6], + real64 ( & elasticStrain)[6] ) const +{ + GEOS_UNUSED_VAR( q ); + real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] ); + real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] ); + + elasticStrain[0] = ( stress[0] - nu*stress[1] - nu*stress[2])/E; + elasticStrain[1] = (-nu*stress[0] + stress[1] - nu*stress[2])/E; + elasticStrain[2] = (-nu*stress[0] - nu*stress[1] + stress[2])/E; + + elasticStrain[3] = stress[3] / m_shearModulus[k]; + elasticStrain[4] = stress[4] / m_shearModulus[k]; + elasticStrain[5] = stress[5] / m_shearModulus[k]; +} + GEOS_HOST_DEVICE inline void ElasticIsotropicUpdates::getElasticStrain( localIndex const k, localIndex const q, real64 ( & elasticStrain)[6] ) const { - real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] ); - real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] ); - elasticStrain[0] = ( m_newStress[k][q][0] - nu*m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E; - elasticStrain[1] = (-nu*m_newStress[k][q][0] + m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E; - elasticStrain[2] = (-nu*m_newStress[k][q][0] - nu*m_newStress[k][q][1] + m_newStress[k][q][2])/E; + real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrain ); - elasticStrain[3] = m_newStress[k][q][3] / m_shearModulus[k]; - elasticStrain[4] = m_newStress[k][q][4] / m_shearModulus[k]; - elasticStrain[5] = m_newStress[k][q][5] / m_shearModulus[k]; } GEOS_HOST_DEVICE @@ -239,16 +261,13 @@ void ElasticIsotropicUpdates::getElasticStrainInc( localIndex const k, localIndex const q, real64 ( & elasticStrainInc)[6] ) const { - real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] ); - real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] ); - elasticStrainInc[0] = ( (m_newStress[k][q][0] - m_oldStress[k][q][0]) - nu*(m_newStress[k][q][1] - m_oldStress[k][q][1]) - nu*(m_newStress[k][q][2] - m_oldStress[k][q][2]))/E; - elasticStrainInc[1] = (-nu*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_newStress[k][q][1] - m_oldStress[k][q][1]) - nu*(m_newStress[k][q][2] - m_oldStress[k][q][2]))/E; - elasticStrainInc[2] = (-nu*(m_newStress[k][q][0] - m_oldStress[k][q][0]) - nu*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + (m_newStress[k][q][2] - m_oldStress[k][q][2]))/E; + real64 stress[6] = + {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrainInc ); - elasticStrainInc[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_shearModulus[k]; - elasticStrainInc[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_shearModulus[k]; - elasticStrainInc[5] = (m_newStress[k][q][5] - m_oldStress[k][q][5]) / m_shearModulus[k]; } GEOS_HOST_DEVICE From 576e7bb5c3ea782d7e9bb2c900d1758ef1ba7a8d Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Mon, 16 Dec 2024 15:49:02 -0800 Subject: [PATCH 43/48] restructure elasticOrthotropic --- .../constitutive/solid/ElasticOrthotropic.hpp | 61 ++++++++++++------- 1 file changed, 38 insertions(+), 23 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index 0958b3b89cf..d9e6d6a4e83 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -172,6 +172,13 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates return LvArray::math::max( LvArray::math::max( m_c44[k], m_c55[k] ), m_c66[k] ); } +protected: + GEOS_HOST_DEVICE + virtual void computeElasticStrain( localIndex const k, + localIndex const q, + real64 const (&stress)[6], + real64 ( &elasticStrain )[6] ) const; + private: /// A reference to the ArrayView holding c11 for each element. arrayView1d< real64 const > const m_c11; @@ -232,26 +239,43 @@ void ElasticOrthotropicUpdates::getElasticStiffness( localIndex const k, GEOS_HOST_DEVICE inline -void ElasticOrthotropicUpdates::getElasticStrain( localIndex const k, - localIndex const q, - real64 ( & elasticStrain)[6] ) const +void ElasticOrthotropicUpdates::computeElasticStrain( localIndex const k, + localIndex const q, + real64 const (&stress)[6], + real64 ( & elasticStrain)[6] ) const { + GEOS_UNUSED_VAR( q ); + real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]); elasticStrain[0] = - ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*m_newStress[k][q][0] + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*m_newStress[k][q][1] + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*m_newStress[k][q][2] ) / + ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*stress[0] + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*stress[1] + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*stress[2] ) / detC; elasticStrain[1] = - ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*m_newStress[k][q][0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][1] + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*m_newStress[k][q][2] ) / + ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*stress[0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[1] + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*stress[2] ) / detC; elasticStrain[2] = - ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*m_newStress[k][q][0] + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*m_newStress[k][q][1] + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*m_newStress[k][q][2] ) / + ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*stress[0] + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*stress[1] + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*stress[2] ) / detC; - elasticStrain[3] = m_newStress[k][q][3] / m_c44[k]; - elasticStrain[4] = m_newStress[k][q][4] / m_c55[k]; - elasticStrain[5] = m_newStress[k][q][5] / m_c66[k]; + elasticStrain[3] = stress[3] / m_c44[k]; + elasticStrain[4] = stress[4] / m_c55[k]; + elasticStrain[5] = stress[5] / m_c66[k]; +} + + +GEOS_HOST_DEVICE +inline +void ElasticOrthotropicUpdates::getElasticStrain( localIndex const k, + localIndex const q, + real64 ( & elasticStrain)[6] ) const +{ + + real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrain ); + } GEOS_HOST_DEVICE @@ -261,21 +285,12 @@ void ElasticOrthotropicUpdates::getElasticStrainInc( localIndex const k, real64 ( & elasticStrainInc)[6] ) const { - real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]); + real64 stress[6] = + {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrainInc ); - elasticStrainInc[0] = - ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + - (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrainInc[1] = - ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + - (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrainInc[2] = - ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + - (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - - elasticStrainInc[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_c44[k]; - elasticStrainInc[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_c55[k]; - elasticStrainInc[5] = (m_newStress[k][q][5] - m_oldStress[k][q][5]) / m_c66[k]; } From 54f91fc8cc3b52ad225db1a7dfc0c92d6399c77c Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Mon, 16 Dec 2024 15:49:33 -0800 Subject: [PATCH 44/48] restructure elasticTransverseIsotropic --- .../solid/ElasticTransverseIsotropic.hpp | 63 ++++++++++++------- 1 file changed, 39 insertions(+), 24 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index 759d4cbcf0d..d78cc24385d 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -165,6 +165,15 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates } +protected: + + GEOS_HOST_DEVICE + virtual void computeElasticStrain( localIndex const k, + localIndex const q, + real64 const (&stress)[6], + real64 ( &elasticStrain )[6] ) const; + + private: /// A reference to the ArrayView holding c11 for each element. @@ -212,22 +221,38 @@ void ElasticTransverseIsotropicUpdates::getElasticStiffness( localIndex const k, GEOS_HOST_DEVICE inline -void ElasticTransverseIsotropicUpdates::getElasticStrain( localIndex const k, - localIndex const q, - real64 ( & elasticStrain)[6] ) const +void ElasticTransverseIsotropicUpdates::computeElasticStrain( localIndex const k, + localIndex const q, + real64 const (&stress)[6], + real64 ( & elasticStrain)[6] ) const { + GEOS_UNUSED_VAR( q ); real64 const c12 = ( m_c11[k] - 2.0 * m_c66[k] ); real64 const detC = m_c11[k]*(m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k]) - c12*(c12*m_c33[k] - m_c13[k]*m_c13[k]) + m_c13[k]*(c12*m_c13[k] - m_c11[k]*m_c13[k]); elasticStrain[0] = - ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][0] + (m_c13[k]*m_c13[k] - c12*m_c33[k])*m_newStress[k][q][1] + (c12*m_c13[k] - m_c13[k]*m_c11[k])*m_newStress[k][q][2] ) / detC; + ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[0] + (m_c13[k]*m_c13[k] - c12*m_c33[k])*stress[1] + (c12*m_c13[k] - m_c13[k]*m_c11[k])*stress[2] ) / detC; elasticStrain[1] = - ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*m_newStress[k][q][0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*m_newStress[k][q][1] + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*m_newStress[k][q][2] ) / detC; - elasticStrain[2] = ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*m_newStress[k][q][0] + (c12*m_c13[k] - m_c11[k]*m_c13[k])*m_newStress[k][q][1] + (m_c11[k]*m_c11[k] - c12*c12)*m_newStress[k][q][2] ) / detC; + ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*stress[0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[1] + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*stress[2] ) / detC; + elasticStrain[2] = ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*stress[0] + (c12*m_c13[k] - m_c11[k]*m_c13[k])*stress[1] + (m_c11[k]*m_c11[k] - c12*c12)*stress[2] ) / detC; + + elasticStrain[3] = stress[3] / m_c44[k]; + elasticStrain[4] = stress[4] / m_c44[k]; + elasticStrain[5] = stress[5] / m_c66[k]; +} + + +GEOS_HOST_DEVICE +inline +void ElasticTransverseIsotropicUpdates::getElasticStrain( localIndex const k, + localIndex const q, + real64 ( & elasticStrain)[6] ) const +{ + + real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrain ); - elasticStrain[3] = m_newStress[k][q][3] / m_c44[k]; - elasticStrain[4] = m_newStress[k][q][4] / m_c44[k]; - elasticStrain[5] = m_newStress[k][q][5] / m_c66[k]; } GEOS_HOST_DEVICE @@ -237,22 +262,12 @@ void ElasticTransverseIsotropicUpdates::getElasticStrainInc( localIndex const k, real64 ( & elasticStrainInc)[6] ) const { - real64 const c12 = ( m_c11[k] - 2.0 * m_c66[k] ); - real64 const detC = m_c11[k]*(m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k]) - c12*(c12*m_c33[k] - m_c13[k]*m_c13[k]) + m_c13[k]*(c12*m_c13[k] - m_c11[k]*m_c13[k]); + real64 stress[6] = + {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + + computeElasticStrain( k, q, stress, elasticStrainInc ); - elasticStrainInc[0] = - ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c13[k]*m_c13[k] - c12*m_c33[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + - (c12*m_c13[k] - m_c13[k]*m_c11[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrainInc[1] = - ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + - (m_c13[k]*c12 - m_c11[k]*m_c13[k])*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - elasticStrainInc[2] = - ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*(m_newStress[k][q][0] - m_oldStress[k][q][0]) + (c12*m_c13[k] - m_c11[k]*m_c13[k])*(m_newStress[k][q][1] - m_oldStress[k][q][1]) + - (m_c11[k]*m_c11[k] - c12*c12)*(m_newStress[k][q][2] - m_oldStress[k][q][2]) ) / detC; - - elasticStrainInc[3] = (m_newStress[k][q][3] - m_oldStress[k][q][3]) / m_c44[k]; - elasticStrainInc[4] = (m_newStress[k][q][4] - m_oldStress[k][q][4]) / m_c44[k]; - elasticStrainInc[5] = (m_newStress[k][q][5] - m_oldStress[k][q][5]) / m_c66[k]; } inline From 93b5268a41c780d6c858dd7beaf4ebd24ca49b60 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Mon, 16 Dec 2024 15:55:49 -0800 Subject: [PATCH 45/48] uncrustify --- src/coreComponents/constitutive/solid/ElasticIsotropic.hpp | 2 +- src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp | 2 +- .../constitutive/solid/ElasticTransverseIsotropic.hpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index 268cf961a58..97573fce5f3 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -264,7 +264,7 @@ void ElasticIsotropicUpdates::getElasticStrainInc( localIndex const k, real64 stress[6] = {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], - m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; computeElasticStrain( k, q, stress, elasticStrainInc ); diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index d9e6d6a4e83..49179c5a95a 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -287,7 +287,7 @@ void ElasticOrthotropicUpdates::getElasticStrainInc( localIndex const k, real64 stress[6] = {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], - m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; computeElasticStrain( k, q, stress, elasticStrainInc ); diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index d78cc24385d..17fabe5479e 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -264,7 +264,7 @@ void ElasticTransverseIsotropicUpdates::getElasticStrainInc( localIndex const k, real64 stress[6] = {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3], - m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; + m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]}; computeElasticStrain( k, q, stress, elasticStrainInc ); From 592b9913036881dad0244cbba90c67864203fe56 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 19 Dec 2024 09:46:12 -0800 Subject: [PATCH 46/48] Update StrainHelper.hpp to output tensor strain not engineering --- .../physicsSolvers/solidMechanics/kernels/StrainHelper.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp index 26fb8b4d09b..faf51f89351 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StrainHelper.hpp @@ -142,16 +142,18 @@ class AverageStrainOverQuadraturePoints : real64 elasticStrainInc[6] = {0.0}; m_solidUpdate.getElasticStrainInc( k, q, elasticStrainInc ); + real64 conversionFactor[6] = {1.0, 1.0, 1.0, 0.5, 0.5, 0.5}; // used for converting from engineering shear to tensor shear + for( int icomp = 0; icomp < 6; ++icomp ) { - m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k]; + m_avgStrain[k][icomp] += conversionFactor[icomp]*detJxW*strain[icomp]/m_elementVolume[k]; // This is a hack to handle boundary conditions such as those seen in plane-strain wellbore problems // Essentially, if bcs are constraining the strain (and thus total displacement), we do not accumulate any plastic strain (regardless // of stresses in material law) if( std::abs( strainInc[icomp] ) > 1.0e-8 ) { - m_avgPlasticStrain[k][icomp] += detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; + m_avgPlasticStrain[k][icomp] += conversionFactor[icomp]*detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; } } } From 2c13f68e2bc46b9c3fccae4280e2b8a1dcd641d8 Mon Sep 17 00:00:00 2001 From: Ryan Aronson Date: Thu, 19 Dec 2024 20:34:47 -0800 Subject: [PATCH 47/48] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 63e604a8ca5..5f8608272be 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3384-9386-cfe31b2 + baseline: integratedTests/baseline_integratedTests-pr3384-9442-381d6ba allow_fail: all: '' streak: '' From 0ce3b3278c0199f4793baa03bd2dc734e441334a Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 7 Jan 2025 14:34:29 -0600 Subject: [PATCH 48/48] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 6f9b1951cf3..6009de53427 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3486-9492-f0c817c + baseline: integratedTests/baseline_integratedTests-pr3384-9542-4e0f693 allow_fail: all: '' streak: ''