Skip to content

Commit

Permalink
Adding field to output total strain (#3142)
Browse files Browse the repository at this point in the history
  • Loading branch information
ryar9534 authored Jun 21, 2024
1 parent 3bf12d2 commit d6b6e0b
Show file tree
Hide file tree
Showing 7 changed files with 268 additions and 1 deletion.
2 changes: 1 addition & 1 deletion .integrated_tests.yaml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
baselines:
bucket: geosx
baseline: integratedTests/baseline_integratedTests-pr3170-5687-d6e16c7
baseline: integratedTests/baseline_integratedTests-pr3142-5705-ebdd3cb

allow_fail:
all: ''
Expand Down
4 changes: 4 additions & 0 deletions BASELINE_NOTES.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 #3142 (2024-06-20)
======================
Adding output of total strain. Rebaseline because of new inclusion of strain in output.

PR #3170 (2024-06-19)
======================
Fix tutorial example for thermal debonding wellbore problem. Test case modified.
Expand Down
9 changes: 9 additions & 0 deletions src/coreComponents/common/DataLayouts.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -154,16 +154,25 @@ namespace cells
/// Cell node map permutation when using cuda.
using NODE_MAP_PERMUTATION = RAJA::PERM_JI;

/// Cell strain permutation when using cuda
using STRAIN_PERM = RAJA::PERM_JI;

#else

/// Cell node map permutation when not using cuda.
using NODE_MAP_PERMUTATION = RAJA::PERM_IJ;

/// Cell strain permutation when not using cuda
using STRAIN_PERM = RAJA::PERM_IJ;

#endif

/// Cell node map unit stride dimension.
static constexpr int NODE_MAP_USD = LvArray::typeManipulation::getStrideOneDimension( NODE_MAP_PERMUTATION {} );

/// Cell strain unit stride dimension
static constexpr int STRAIN_USD = LvArray::typeManipulation::getStrideOneDimension( STRAIN_PERM {} );

} // namespace cells

namespace solid
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,10 @@ using array2dLayoutIncrDisplacement = array2d< real64, nodes::INCR_DISPLACEMENT_
using arrayView2dLayoutIncrDisplacement = arrayView2d< real64, nodes::INCR_DISPLACEMENT_USD >;
using arrayViewConst2dLayoutIncrDisplacement = arrayView2d< real64 const, nodes::INCR_DISPLACEMENT_USD >;

using array2dLayoutStrain = array2d< real64, cells::STRAIN_PERM >;
using arrayView2dLayoutStrain = arrayView2d< real64, cells::STRAIN_USD >;
using arrayViewConst2dLayoutStrain = arrayView2d< real64 const, cells::STRAIN_USD >;

using array2dLayoutVelocity = array2d< real64, nodes::VELOCITY_PERM >;
using arrayView2dLayoutVelocity = arrayView2d< real64, nodes::VELOCITY_USD >;
using arrayViewConst2dLayoutVelocity = arrayView2d< real64 const, nodes::VELOCITY_USD >;
Expand Down Expand Up @@ -74,6 +78,14 @@ DECLARE_FIELD( incrementalDisplacement,
WRITE_AND_READ,
"Incremental displacements for the current time step on the nodes" );

DECLARE_FIELD( strain,
"strain",
array2dLayoutStrain,
0,
LEVEL_0,
WRITE_AND_READ,
"Average strain in cell" );

DECLARE_FIELD( incrementalBubbleDisplacement,
"incrementalBubbleDisplacement",
array2d< real64 >,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,8 @@ void SolidMechanicsLagrangianFEM::registerDataOnMesh( Group & meshBodies )
ElementSubRegionBase & subRegion )
{
setConstitutiveNamesCallSuper( subRegion );

subRegion.registerField< solidMechanics::strain >( getName() ).reference().resizeDimension< 1 >( 6 );
} );

NodeManager & nodes = meshLevel.getNodeManager();
Expand Down Expand Up @@ -904,6 +906,9 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS
solidMechanics::arrayView2dLayoutIncrDisplacement const uhat =
nodeManager.getField< solidMechanics::incrementalDisplacement >();

solidMechanics::arrayView2dLayoutTotalDisplacement const disp =
nodeManager.getField< solidMechanics::totalDisplacement >();

if( this->m_timeIntegrationOption == TimeIntegrationOption::ImplicitDynamic )
{
solidMechanics::arrayView2dLayoutAcceleration const a_n = nodeManager.getField< solidMechanics::acceleration >();
Expand Down Expand Up @@ -933,6 +938,23 @@ 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 >();

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,
mesh.getEdgeManager(),
mesh.getFaceManager(),
subRegion,
finiteElement,
disp,
strain );
} );


} );
} );

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include "codingUtilities/EnumStrings.hpp"
#include "common/TimingMacros.hpp"
#include "kernels/SolidMechanicsLagrangianFEMKernels.hpp"
#include "kernels/StrainHelper.hpp"
#include "mesh/MeshForLoopInterface.hpp"
#include "mesh/mpiCommunications/CommunicationTools.hpp"
#include "mesh/mpiCommunications/MPI_iCommData.hpp"
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,219 @@
/*
* ------------------------------------------------------------------------------------------------------------
* SPDX-License-Identifier: LGPL-2.1-only
*
* Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
* Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
* Copyright (c) 2018-2020 TotalEnergies
* Copyright (c) 2019- GEOSX Contributors
* All rights reserved
*
* See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
* ------------------------------------------------------------------------------------------------------------
*/

/**
* @file StrainHelper.hpp
*/

#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_
#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_

#include "common/DataTypes.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
#include "finiteElement/FiniteElementDispatch.hpp"
#include "mesh/CellElementSubRegion.hpp"
#include "mesh/utilities/AverageOverQuadraturePointsKernel.hpp"
#include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp"

namespace geos
{
/**
* @class AverageStrainOverQuadraturePoints
* @tparam SUBREGION_TYPE the subRegion type
* @tparam FE_TYPE the finite element type
*/
template< typename SUBREGION_TYPE,
typename FE_TYPE >
class AverageStrainOverQuadraturePoints :
public AverageOverQuadraturePointsBase< SUBREGION_TYPE,
FE_TYPE >
{
public:

/// Alias for the base class;
using Base = AverageOverQuadraturePointsBase< SUBREGION_TYPE,
FE_TYPE >;

using Base::m_elementVolume;
using Base::m_elemsToNodes;
using Base::m_finiteElementSpace;

/**
* @brief Constructor for the class
* @param nodeManager the node manager
* @param edgeManager the edge manager
* @param faceManager the face manager
* @param elementSubRegion the element subRegion
* @param finiteElementSpace the finite element space
* @param displacement the displacement solution field
* @param avgStrain the strain averaged over quadrature points
*/
AverageStrainOverQuadraturePoints( NodeManager & nodeManager,
EdgeManager const & edgeManager,
FaceManager const & faceManager,
SUBREGION_TYPE const & elementSubRegion,
FE_TYPE const & finiteElementSpace,
fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
fields::solidMechanics::arrayView2dLayoutStrain const avgStrain ):
Base( nodeManager,
edgeManager,
faceManager,
elementSubRegion,
finiteElementSpace ),
m_displacement( displacement ),
m_avgStrain( avgStrain )
{}

/**
* @copydoc finiteElement::KernelBase::StackVariables
*/
struct StackVariables : Base::StackVariables
{real64 uLocal[FE_TYPE::maxSupportPoints][3]; };

/**
* @brief Performs the setup phase for the kernel.
* @param k The element index.
* @param stack The StackVariable object that hold the stack variables.
*/
GEOS_HOST_DEVICE
void setup( localIndex const k,
StackVariables & stack ) const
{
Base::setup( k, stack );

for( localIndex a = 0; a < FE_TYPE::maxSupportPoints; ++a )
{
localIndex const localNodeIndex = m_elemsToNodes( k, a );
for( int i = 0; i < 3; ++i )
{
stack.uLocal[a][i] = m_displacement[localNodeIndex][i];
}
}

for( int icomp = 0; icomp < 6; ++icomp )
{
m_avgStrain[k][icomp] = 0.0;
}
}

/**
* @brief Increment the average property with the contribution of the property at this quadrature point
* @param k The element index
* @param q The quadrature point index
* @param stack The StackVariables object that hold the stack variables.
*/
GEOS_HOST_DEVICE
void quadraturePointKernel( localIndex const k,
localIndex const q,
StackVariables & stack ) const
{
//real64 const weight = FE_TYPE::transformedQuadratureWeight( q, stack.xLocal, stack.feStack ) / m_elementVolume[k];

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};
FE_TYPE::symmetricGradient( dNdX, stack.uLocal, strain );

for( int icomp = 0; icomp < 6; ++icomp )
{
m_avgStrain[k][icomp] += detJxW*strain[icomp]/m_elementVolume[k];
}
}

/**
* @brief Launch the kernel over the elements in the subRegion
* @tparam POLICY the kernel policy
* @tparam KERNEL_TYPE the type of kernel
* @param numElems the number of elements in the subRegion
* @param kernelComponent the kernel component
*/
template< typename POLICY,
typename KERNEL_TYPE >
static void
kernelLaunch( localIndex const numElems,
KERNEL_TYPE const & kernelComponent )
{
forAll< POLICY >( numElems,
[=] GEOS_HOST_DEVICE ( localIndex const k )
{
typename KERNEL_TYPE::StackVariables stack;

kernelComponent.setup( k, stack );
for( integer q = 0; q < FE_TYPE::numQuadraturePoints; ++q )
{
kernelComponent.quadraturePointKernel( k, q, stack );
}
} );
}

protected:

/// The displacement solution
fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const m_displacement;

/// The average strain
fields::solidMechanics::arrayView2dLayoutStrain const m_avgStrain;

};



/**
* @class AverageStrainOverQuadraturePointsKernelFactory
* @brief Class to create and launch the kernel
*/
class AverageStrainOverQuadraturePointsKernelFactory
{
public:

/**
* @brief Create a new kernel and launch
* @tparam SUBREGION_TYPE the subRegion type
* @tparam FE_TYPE the finite element type
* @tparam POLICY the kernel policy
* @param nodeManager the node manager
* @param edgeManager the edge manager
* @param faceManager the face manager
* @param elementSubRegion the element subRegion
* @param finiteElementSpace the finite element space
* @param property the property at quadrature points
* @param averageProperty the property averaged over quadrature points
*/
template< typename SUBREGION_TYPE,
typename FE_TYPE,
typename POLICY >
static void
createAndLaunch( NodeManager & nodeManager,
EdgeManager const & edgeManager,
FaceManager const & faceManager,
SUBREGION_TYPE const & elementSubRegion,
FE_TYPE const & finiteElementSpace,
fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
fields::solidMechanics::arrayView2dLayoutStrain const avgStrain )
{
AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >
kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace,
displacement, avgStrain );

AverageStrainOverQuadraturePoints< SUBREGION_TYPE, FE_TYPE >::template
kernelLaunch< POLICY >( elementSubRegion.size(), kernel );
}
};



}


#endif /* GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_ */

0 comments on commit d6b6e0b

Please sign in to comment.