diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp index 3faa6bdff59..0dc71c4edec 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp @@ -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, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 17b720bbd64..10213eef099 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -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 ); @@ -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(); @@ -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]; @@ -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 ); } ); } ); @@ -421,9 +460,9 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepSetup( real64 const & domain.getNeighbors(), true ); } ); - } + void SolidMechanicsAugmentedLagrangianContact::assembleSystem( real64 const time, real64 const dt, DomainPartition & domain, @@ -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] }; @@ -1056,6 +1101,7 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit { solidMechanicsALMKernels::ComputeTractionSimultaneousKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), + ghostRank, iterativePenalty, traction, dispJump, @@ -1067,6 +1113,7 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit solidMechanicsALMKernels::ComputeTractionKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), frictionWrapper, + ghostRank, iterativePenalty, traction, dispJump, @@ -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] ); } ); } ); @@ -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 = @@ -1224,6 +1279,7 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit solidMechanicsALMKernels::UpdateStateKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), frictionWrapper, + ghostRank, oldDispJump, dispJump, iterativePenalty, @@ -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]; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp index d73fb84619a..8f988fedca9 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp @@ -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(); @@ -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 ) @@ -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] ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernels.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernels.hpp index 796570528a2..7af616df223 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernels.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernels.hpp @@ -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 @@ -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, @@ -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] ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp index 213877afb03..05e493c0a06 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp @@ -102,6 +102,7 @@ struct UpdateStateKernel * @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] oldDispJump the array containing the old displacement jump (previous time step) * @param[in] dispJump the array containing the displacement jump * @param[in] penalty the array containing the penalty coefficients @@ -114,6 +115,7 @@ struct UpdateStateKernel static void launch( localIndex const size, CONTACT_WRAPPER const & contactWrapper, + arrayView1d< integer const > const & ghostRank, arrayView2d< real64 const > const & oldDispJump, arrayView2d< real64 const > const & dispJump, arrayView2d< real64 > const & penalty, @@ -125,6 +127,10 @@ struct UpdateStateKernel { forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { + if( ghostRank[k] >= 0 ) + { + return; + } real64 const zero = LvArray::NumericLimits< real64 >::epsilon; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp index 07f526cdece..60fb3132e13 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp @@ -415,6 +415,7 @@ struct ComputeTractionSimultaneousKernel * @brief Launch the kernel function to comute traction * @tparam POLICY the type of policy used in the kernel launch * @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 @@ -424,6 +425,7 @@ struct ComputeTractionSimultaneousKernel template< typename POLICY > static void launch( localIndex const size, + arrayView1d< integer const > const & ghostRank, arrayView2d< real64 const > const & penalty, arrayView2d< real64 const > const & traction, arrayView2d< real64 const > const & dispJump, @@ -433,6 +435,10 @@ struct ComputeTractionSimultaneousKernel forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const kfe ) { + if( ghostRank[kfe] >= 0 ) + { + return; + } tractionNew[kfe][0] = traction[kfe][0] + penalty[kfe][0] * dispJump[kfe][0]; tractionNew[kfe][1] = traction[kfe][1] + ( penalty[kfe][2] * deltaDispJump[kfe][1]+ penalty[kfe][4] * deltaDispJump[kfe][2] ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp index 84308200120..46674b7d115 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp @@ -270,9 +270,10 @@ struct ComputeRotationMatricesKernel { /** - * @brief Launch the kernel function to comute rotation matrices + * @brief Launch the kernel function to compute rotation matrices * @tparam POLICY the type of policy used in the kernel launch * @param[in] size the size of the subregion + * @param[in] ghostRank the ghost rank of each element * @param[in] faceNormal the array of array containing the face to nodes map * @param[in] elemsToFaces the array of array containing the element to faces map * @param[out] rotationMatrix the array containing the rotation matrices @@ -280,6 +281,7 @@ struct ComputeRotationMatricesKernel template< typename POLICY > static void launch( localIndex const size, + arrayView1d< integer const > const & ghostRank, arrayView2d< real64 const > const & faceNormal, arrayView2d< localIndex const > const & elemsToFaces, arrayView3d< real64 > const & rotationMatrix, @@ -290,6 +292,10 @@ struct ComputeRotationMatricesKernel forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { + if( ghostRank[k] >= 0 ) + { + return; + } localIndex const f0 = elemsToFaces[k][0]; localIndex const f1 = elemsToFaces[k][1];