Skip to content

feat: Prevent ghost elements from modifying their own data in contact solvers #4000

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bd713 wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

feat: Prevent ghost elements from modifying their own data in contact solvers #4000

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,8 @@ void ContactSolverBase::synchronizeFractureState( DomainPartition & domain ) con
FieldIdentifiers fieldsToBeSync;

fieldsToBeSync.addElementFields( { contact::fractureState::key(),
contact::traction::key() },
contact::traction::key(),
contact::iterativePenalty::key() },
{ getUniqueFractureRegionName() } );

CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync,
Expand Down
104 changes: 82 additions & 22 deletions ...onents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -255,6 +255,9 @@ void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & do
// Create the lists of interface elements that have same type.
createFaceTypeList( domain );

// Synchronize fracture state before reading it in ghost elements
synchronizeFractureState( domain );

// Create the lists of interface elements that have same type and same fracture state.
updateStickSlipList( domain );

Expand Down Expand Up @@ -324,14 +327,13 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepSetup( real64 const &
real64 const & dt,
DomainPartition & domain )
{

SolidMechanicsLagrangianFEM::implicitStepSetup( time_n, dt, domain );

// Compute rotation matrices (only for owned elements)
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
string_array const & )
{

FaceManager & faceManager = mesh.getFaceManager();
ElementRegionManager & elemManager = mesh.getElemManager();

Expand All @@ -341,44 +343,75 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepSetup( real64 const &
arrayView2d< real64 const > const faceNormal = faceManager.faceNormal();
arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst();

arrayView2d< real64 > const incrBubbleDisp =
faceManager.getField< contact::incrementalBubbleDisplacement >();

arrayView3d< real64 > const
rotationMatrix = subRegion.getField< contact::rotationMatrix >().toView();

arrayView3d< real64 > const rotationMatrix = subRegion.getField< contact::rotationMatrix >().toView();
arrayView2d< real64 > const unitNormal = subRegion.getNormalVector();
arrayView2d< real64 > const unitTangent1 = subRegion.getTangentVector1();
arrayView2d< real64 > const unitTangent2 = subRegion.getTangentVector2();

// Compute rotation matrices
arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

solidMechanicsConformingContactKernels::ComputeRotationMatricesKernel::
launch< parallelDevicePolicy<> >( subRegion.size(),
ghostRank,
faceNormal,
elemsToFaces,
rotationMatrix,
unitNormal,
unitTangent1,
unitTangent2 );
} );

// Synchronize rotation matrices to ghost elements
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
string_array const & )
{
FieldIdentifiers fieldsToBeSync;

fieldsToBeSync.addElementFields( { contact::rotationMatrix::key() },
{ getUniqueFractureRegionName() } );

CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync,
mesh,
domain.getNeighbors(),
true );
} );

// Compute tolerances and initialize fields
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
string_array const & )
{
FaceManager & faceManager = mesh.getFaceManager();
ElementRegionManager & elemManager = mesh.getElemManager();

// Set the tollerances
SurfaceElementRegion & region = elemManager.getRegion< SurfaceElementRegion >( getUniqueFractureRegionName() );
FaceElementSubRegion & subRegion = region.getUniqueSubRegion< FaceElementSubRegion >();

arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst();
arrayView2d< real64 > const incrBubbleDisp = faceManager.getField< contact::incrementalBubbleDisplacement >();
arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

// Set the tolerances
computeTolerances( domain );

// Set array to update penalty coefficients
arrayView2d< real64 > const dispJumpUpdPenalty =
subRegion.getReference< array2d< real64 > >( viewKeyStruct::dispJumpUpdPenaltyString() );

arrayView2d< real64 > const
iterativePenalty = subRegion.getField< contact::iterativePenalty >().toView();
arrayView2d< real64 > const iterativePenalty = subRegion.getField< contact::iterativePenalty >().toView();
arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >();

if( m_simultaneous )
{
// Set the iterative penalty coefficients
forAll< parallelDevicePolicy<> >( subRegion.size(),
[=]
GEOS_HOST_DEVICE ( localIndex const k )
[=] GEOS_HOST_DEVICE ( localIndex const k )
{
if( ghostRank[k] >= 0 )
{
return;
}
if( fractureState[k] == contact::FractureState::Stick )
{
iterativePenalty[k][2] = iterativePenalty[k][1];
Expand All @@ -394,15 +427,21 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepSetup( real64 const &
} );
}

// Initialize element fields (only for owned elements)
forAll< parallelDevicePolicy<> >( subRegion.size(),
[ = ]
GEOS_HOST_DEVICE ( localIndex const k )
[ = ] GEOS_HOST_DEVICE ( localIndex const k )
{
if( ghostRank[k] < 0 )
{
LvArray::tensorOps::fill< 3 >( dispJumpUpdPenalty[k], 0.0 );
}
} );

// Initialize face fields directly
forAll< parallelDevicePolicy<> >( faceManager.size(),
[ = ] GEOS_HOST_DEVICE ( localIndex const kf )
{
LvArray::tensorOps::fill< 3 >( dispJumpUpdPenalty[k], 0.0 );
localIndex const kf0 = elemsToFaces[k][0];
localIndex const kf1 = elemsToFaces[k][1];
LvArray::tensorOps::fill< 3 >( incrBubbleDisp[kf0], 0.0 );
LvArray::tensorOps::fill< 3 >( incrBubbleDisp[kf1], 0.0 );
LvArray::tensorOps::fill< 3 >( incrBubbleDisp[kf], 0.0 );
} );
} );

Expand All @@ -421,9 +460,9 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepSetup( real64 const &
domain.getNeighbors(),
true );
} );

}


void SolidMechanicsAugmentedLagrangianContact::assembleSystem( real64 const time,
real64 const dt,
DomainPartition & domain,
Expand Down Expand Up @@ -757,10 +796,16 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepComplete( real64 cons
arrayView1d< real64 > const slip = subRegion.getField< contact::slip >();
arrayView1d< real64 > const tangentialTraction = subRegion.getField< contact::tangentialTraction >();

arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

forAll< parallelDevicePolicy<> >( subRegion.size(),
[ = ]
GEOS_HOST_DEVICE ( localIndex const kfe )
{
if( ghostRank[kfe] >= 0 )
{
return;
}
// Compute the slip
real64 const deltaDisp[2] = { deltaDispJump[kfe][1],
deltaDispJump[kfe][2] };
Expand Down Expand Up @@ -1056,6 +1101,7 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit
{
solidMechanicsALMKernels::ComputeTractionSimultaneousKernel::
launch< parallelDevicePolicy<> >( subRegion.size(),
ghostRank,
iterativePenalty,
traction,
dispJump,
Expand All @@ -1067,6 +1113,7 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit
solidMechanicsALMKernels::ComputeTractionKernel::
launch< parallelDevicePolicy<> >( subRegion.size(),
frictionWrapper,
ghostRank,
iterativePenalty,
traction,
dispJump,
Expand Down Expand Up @@ -1173,11 +1220,17 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit
FaceElementSubRegion & subRegion )
{

arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

arrayView2d< real64 > const traction_new_v = traction_new.toView();
arrayView2d< real64 > const traction = subRegion.getField< contact::traction >();

forAll< parallelDevicePolicy<> >( subRegion.size(), [ = ] GEOS_HOST_DEVICE ( localIndex const kfe )
{
if( ghostRank[kfe] >= 0 )
{
return;
}
LvArray::tensorOps::copy< 3 >( traction[kfe], traction_new_v[kfe] );
} );
} );
Expand All @@ -1198,6 +1251,8 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit
string const & frictionLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() );
FrictionBase const & frictionLaw = getConstitutiveModel< FrictionBase >( subRegion, frictionLawName );

arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

arrayView2d< real64 > const traction = subRegion.getField< contact::traction >();

arrayView1d< real64 const > const normalTractionTolerance =
Expand All @@ -1224,6 +1279,7 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit
solidMechanicsALMKernels::UpdateStateKernel::
launch< parallelDevicePolicy<> >( subRegion.size(),
frictionWrapper,
ghostRank,
oldDispJump,
dispJump,
iterativePenalty,
Expand All @@ -1236,6 +1292,10 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit
forAll< parallelDevicePolicy<> >( subRegion.size(), [ = ]
GEOS_HOST_DEVICE ( localIndex const kfe )
{
if( ghostRank[kfe] >= 0 )
{
return;
}
dispJumpUpdPenalty[kfe][0] = dispJump[kfe][0];
dispJumpUpdPenalty[kfe][1] = deltaDispJump[kfe][1];
dispJumpUpdPenalty[kfe][2] = deltaDispJump[kfe][2];
Expand Down
71 changes: 55 additions & 16 deletions ...ponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -318,11 +318,11 @@ void SolidMechanicsLagrangeContactBubbleStab::setSparsityPattern( DomainPartitio

void SolidMechanicsLagrangeContactBubbleStab::computeRotationMatrices( DomainPartition & domain ) const
{
// Compute rotation matrices (owned only)
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
string_array const & )
{

FaceManager & faceManager = mesh.getFaceManager();
ElementRegionManager & elemManager = mesh.getElemManager();

Expand All @@ -342,27 +342,60 @@ void SolidMechanicsLagrangeContactBubbleStab::computeRotationMatrices( DomainPar
arrayView2d< real64 > const unitTangent1 = subRegion.getTangentVector1();
arrayView2d< real64 > const unitTangent2 = subRegion.getTangentVector2();

// Compute rotation matrices
solidMechanicsConformingContactKernels::ComputeRotationMatricesKernel::launch< parallelDevicePolicy<> >( subRegion.size(),
faceNormal,
elemsToFaces,
rotationMatrix,
unitNormal,
unitTangent1,
unitTangent2 );
arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

solidMechanicsConformingContactKernels::ComputeRotationMatricesKernel::launch< parallelDevicePolicy<> >(
subRegion.size(),
ghostRank,
faceNormal,
elemsToFaces,
rotationMatrix,
unitNormal,
unitTangent1,
unitTangent2 );
} );

forAll< parallelDevicePolicy<> >( subRegion.size(),
[ = ]
GEOS_HOST_DEVICE ( localIndex const k )
// Synchronize rotation matrices to ghost elements
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
string_array const & )
{
FieldIdentifiers fieldsToBeSync;

fieldsToBeSync.addElementFields( { contact::rotationMatrix::key() },
{ getUniqueFractureRegionName() } );

CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync,
mesh,
domain.getNeighbors(),
true );
} );

// Initialize bubble displacement
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
string_array const & )
{
FaceManager & faceManager = mesh.getFaceManager();
ElementRegionManager & elemManager = mesh.getElemManager();

SurfaceElementRegion & region = elemManager.getRegion< SurfaceElementRegion >( getUniqueFractureRegionName() );
FaceElementSubRegion & subRegion = region.getUniqueSubRegion< FaceElementSubRegion >();

arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst();
arrayView2d< real64 > const incrBubbleDisp =
faceManager.getField< contact::incrementalBubbleDisplacement >();

// Initialize face fields directly
forAll< parallelDevicePolicy<> >( faceManager.size(),
[ = ] GEOS_HOST_DEVICE ( localIndex const kf )
{
localIndex const kf0 = elemsToFaces[k][0];
localIndex const kf1 = elemsToFaces[k][1];
LvArray::tensorOps::fill< 3 >( incrBubbleDisp[kf0], 0.0 );
LvArray::tensorOps::fill< 3 >( incrBubbleDisp[kf1], 0.0 );
LvArray::tensorOps::fill< 3 >( incrBubbleDisp[kf], 0.0 );
} );
} );
}


void SolidMechanicsLagrangeContactBubbleStab::implicitStepSetup( real64 const & time_n,
real64 const & dt,
DomainPartition & domain )
Expand Down Expand Up @@ -505,8 +538,14 @@ void SolidMechanicsLagrangeContactBubbleStab::implicitStepComplete( real64 const
arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >();
arrayView2d< real64 > const oldDispJump = subRegion.getField< contact::oldDispJump >();

arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe )
{
if( ghostRank[kfe] >= 0 )
{
return;
}
LvArray::tensorOps::fill< 3 >( deltaDispJump[kfe], 0.0 );
LvArray::tensorOps::fill< 3 >( deltaTraction[kfe], 0.0 );
LvArray::tensorOps::copy< 3 >( oldDispJump[kfe], dispJump[kfe] );
Expand Down
6 changes: 6 additions & 0 deletions ...coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernels.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -401,6 +401,7 @@ struct ComputeTractionKernel
* @tparam POLICY the type of policy used in the kernel launch
* @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates
* @param[in] size the size of the subregion
* @param[in] ghostRank the ghost rank of each element
* @param[in] penalty the array containing the tangential penalty matrix
* @param[in] traction the array containing the current traction
* @param[in] dispJump the array containing the displacement jump
Expand All @@ -411,6 +412,7 @@ struct ComputeTractionKernel
static void
launch( localIndex const size,
CONTACT_WRAPPER const & contactWrapper,
arrayView1d< integer const > const & ghostRank,
arrayView2d< real64 const > const & penalty,
arrayView2d< real64 const > const & traction,
arrayView2d< real64 const > const & dispJump,
Expand All @@ -420,6 +422,10 @@ struct ComputeTractionKernel

forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
{
if( ghostRank[k] >= 0 )
{
return;
}

contactWrapper.updateTractionOnly( dispJump[k], deltaDispJump[k],
penalty[k], traction[k], tractionNew[k] );
Expand Down
Loading
Loading