diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 8819a54e041..8bc92055db7 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3970-15479-074f42a + baseline: integratedTests/baseline_integratedTests-pr3986-15734-7487221 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 0e0ad277a6d..3efe0e4b57c 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 #3970 (2026-02-27) +===================== +Corrected traction boundary conditions + PR #3970 (2026-02-11) ===================== Bypass well residual calculation for closed wells diff --git a/src/coreComponents/fieldSpecification/TractionBoundaryCondition.cpp b/src/coreComponents/fieldSpecification/TractionBoundaryCondition.cpp index 2dab7ebfa4f..5a98edb67c2 100644 --- a/src/coreComponents/fieldSpecification/TractionBoundaryCondition.cpp +++ b/src/coreComponents/fieldSpecification/TractionBoundaryCondition.cpp @@ -19,6 +19,8 @@ #include "TractionBoundaryCondition.hpp" +#include "finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp" +#include "finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp" #include "functions/TableFunction.hpp" namespace geos @@ -132,15 +134,191 @@ void TractionBoundaryCondition::initializePreSubGroups() } +/** + * @brief Integrate a uniform traction over a face using a standard FE face formulation. + * + * @tparam FE_TYPE The finite element face type (e.g. H1_TriangleFace_Lagrange1_Gauss1_impl, + * H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl). + * Must provide: numNodes, numQuadraturePoints, calcN(), transformedQuadratureWeight(). + * @param traction The traction vector [3] to integrate. + * @param xFace The face node coordinates in FE ordering [numNodes][3]. + * @param permutation Mapping from FE node index to mesh face-local node index. + * Accounts for ordering differences (e.g. tensor-product vs CCW). + * @param faceToNodeMap The face-to-node connectivity map. + * @param kf The face index. + * @param blockLocalDofNumber Block-local DOF numbering. + * @param dofRankOffset The rank offset for DOFs. + * @param localRhs The local RHS vector to assemble into. + */ +template< typename FE_TYPE > +GEOS_HOST_DEVICE +inline +void integrateFaceTraction( real64 const ( &traction )[3], + real64 const ( &xFace )[FE_TYPE::numNodes][3], + int const ( &permutation )[FE_TYPE::numNodes], + ArrayOfArraysView< localIndex const > const & faceToNodeMap, + localIndex const kf, + arrayView1d< globalIndex const > const & blockLocalDofNumber, + globalIndex const dofRankOffset, + arrayView1d< real64 > const & localRhs ) +{ + constexpr localIndex numNodes = FE_TYPE::numNodes; + constexpr localIndex numQuadraturePoints = FE_TYPE::numQuadraturePoints; + + // Accumulate per-node integration weight: w_a = sum_q N_a(q) * |J(q)| * wq + real64 nodalWeight[numNodes] {}; + + for( localIndex q = 0; q < numQuadraturePoints; ++q ) + { + real64 N[numNodes]; + FE_TYPE::calcN( q, N ); + + real64 const detJxW = FE_TYPE::transformedQuadratureWeight( q, xFace ); + + for( localIndex a = 0; a < numNodes; ++a ) + { + nodalWeight[a] += N[a] * detJxW; + } + } + + // Assemble traction * nodalWeight into the RHS. + // Use the permutation to map FE node 'a' back to the mesh face-local node index. + for( localIndex a = 0; a < numNodes; ++a ) + { + localIndex const meshNodeIdx = permutation[a]; + localIndex const dof = blockLocalDofNumber[ faceToNodeMap( kf, meshNodeIdx ) ] - dofRankOffset; + if( dof < 0 || dof >= localRhs.size() ) + continue; + RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+0], traction[0] * nodalWeight[a] ); + RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+1], traction[1] * nodalWeight[a] ); + RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+2], traction[2] * nodalWeight[a] ); + } +} + +/** + * @brief Integrate a uniform traction over a general polygon face. + * + * This is a fallback for faces that are not triangles or quadrilaterals. + * + * @param traction The traction vector [3] to integrate. + * @param faceArea The area of the face. + * @param numNodes The number of nodes of the face. + * @param faceToNodeMap The face-to-node connectivity map. + * @param kf The face index. + * @param blockLocalDofNumber Block-local DOF numbering. + * @param dofRankOffset The rank offset for DOFs. + * @param localRhs The local RHS vector to assemble into. + */ +GEOS_HOST_DEVICE +inline +void integrateFaceTractionPolygon( real64 const ( &traction )[3], + real64 const faceArea, + localIndex const numNodes, + ArrayOfArraysView< localIndex const > const & faceToNodeMap, + localIndex const kf, + arrayView1d< globalIndex const > const & blockLocalDofNumber, + globalIndex const dofRankOffset, + arrayView1d< real64 > const & localRhs ) +{ + real64 const w = faceArea / numNodes; + for( localIndex a = 0; a < numNodes; ++a ) + { + localIndex const dof = blockLocalDofNumber[ faceToNodeMap( kf, a ) ] - dofRankOffset; + if( dof < 0 || dof >= localRhs.size() ) + continue; + RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+0], traction[0] * w ); + RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+1], traction[1] * w ); + RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+2], traction[2] * w ); + } +} + +/** + * @brief Dispatch face traction integration to the appropriate FE formulation. + * + * Selects the FE face type based on numNodes: + * - 3 nodes: H1_TriangleFace_Lagrange1_Gauss1_impl + * - 4 nodes: H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl + * - Other: Equal distribution fallback + * + * @param traction The traction vector [3]. + * @param numNodes The number of nodes of the face. + * @param faceArea The area of the face (only used for polygon fallback). + * @param faceToNodeMap The face-to-node connectivity map. + * @param kf The face index. + * @param nodePositions The reference positions of all nodes. + * @param blockLocalDofNumber Block-local DOF numbering. + * @param dofRankOffset The rank offset for DOFs. + * @param localRhs The local RHS vector to assemble into. + */ +GEOS_HOST_DEVICE +inline +void assembleFaceTraction( real64 const ( &traction )[3], + localIndex const numNodes, + real64 const faceArea, + ArrayOfArraysView< localIndex const > const & faceToNodeMap, + localIndex const kf, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePositions, + arrayView1d< globalIndex const > const & blockLocalDofNumber, + globalIndex const dofRankOffset, + arrayView1d< real64 > const & localRhs ) +{ + using TriFace = finiteElement::H1_TriangleFace_Lagrange1_Gauss1_impl; + using QuadFace = finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl; + + if( numNodes == 3 ) + { + // Get the permutation mapping from FE node ordering to mesh face-local node ordering. + int permutation[TriFace::numNodes]; + TriFace::getPermutation( permutation ); + + // Gather coordinates in FE ordering: xFace[feNode] = coords of mesh node permutation[feNode] + real64 xFace[TriFace::numNodes][3]; + for( localIndex a = 0; a < TriFace::numNodes; ++a ) + { + localIndex const nodeIdx = faceToNodeMap( kf, permutation[a] ); + xFace[a][0] = nodePositions( nodeIdx, 0 ); + xFace[a][1] = nodePositions( nodeIdx, 1 ); + xFace[a][2] = nodePositions( nodeIdx, 2 ); + } + integrateFaceTraction< TriFace >( traction, xFace, permutation, faceToNodeMap, kf, + blockLocalDofNumber, dofRankOffset, localRhs ); + } + else if( numNodes == 4 ) + { + // Get the permutation mapping from FE node ordering to mesh face-local node ordering. + // For quads: FE uses Z-ordering (tensor product), mesh uses CCW ordering. + // getPermutation returns {0, 1, 3, 2} meaning FE node 2 → mesh node 3 and vice versa. + int permutation[QuadFace::numNodes]; + QuadFace::getPermutation( permutation ); + + // Gather coordinates in FE ordering: xFace[feNode] = coords of mesh node permutation[feNode] + real64 xFace[QuadFace::numNodes][3]; + for( localIndex a = 0; a < QuadFace::numNodes; ++a ) + { + localIndex const nodeIdx = faceToNodeMap( kf, permutation[a] ); + xFace[a][0] = nodePositions( nodeIdx, 0 ); + xFace[a][1] = nodePositions( nodeIdx, 1 ); + xFace[a][2] = nodePositions( nodeIdx, 2 ); + } + integrateFaceTraction< QuadFace >( traction, xFace, permutation, faceToNodeMap, kf, + blockLocalDofNumber, dofRankOffset, localRhs ); + } + else + { + integrateFaceTractionPolygon( traction, faceArea, numNodes, faceToNodeMap, kf, + blockLocalDofNumber, dofRankOffset, localRhs ); + } +} + + void TractionBoundaryCondition::launch( real64 const time, arrayView1d< globalIndex const > const blockLocalDofNumber, globalIndex const dofRankOffset, FaceManager const & faceManager, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePositions, SortedArrayView< localIndex const > const & targetSet, arrayView1d< real64 > const & localRhs ) const { - - arrayView1d< real64 const > const faceArea = faceManager.faceArea(); arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); @@ -193,11 +371,9 @@ void TractionBoundaryCondition::launch( real64 const time, real64 const nodalScale = ( nodalScaleFlag == 1 ) ? scaleSet[i] : 1.0; - // TODO consider dispatch if appropriate real64 const tractionMagnitude = spatialFunction ? tractionMagnitudeArrayView[i] : (tractionMagnitude0 * nodalScale); real64 traction[3] = { 0 }; - // TODO consider dispatch if appropriate if( tractionType == TractionType::vector ) { traction[0] = tractionMagnitude * direction[0]; @@ -222,23 +398,11 @@ void TractionBoundaryCondition::launch( real64 const time, traction[1] = inputStress[5] * temp[0] + inputStress[1] * temp[1] + inputStress[3] * temp[2]; traction[2] = inputStress[4] * temp[0] + inputStress[3] * temp[1] + inputStress[2] * temp[2]; } - } - // TODO replace with proper FEM integration - traction[0] *= faceArea[kf] / numNodes; - traction[1] *= faceArea[kf] / numNodes; - traction[2] *= faceArea[kf] / numNodes; - - for( localIndex a=0; a= localRhs.size() ) - continue; - RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+0], traction[0] ); - RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+1], traction[1] ); - RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+2], traction[2] ); - } + // Dispatch face traction integration to the appropriate FE formulation + assembleFaceTraction( traction, numNodes, faceArea[kf], faceToNodeMap, kf, + nodePositions, blockLocalDofNumber, dofRankOffset, localRhs ); } ); } } diff --git a/src/coreComponents/fieldSpecification/TractionBoundaryCondition.hpp b/src/coreComponents/fieldSpecification/TractionBoundaryCondition.hpp index 466f8cd4d2c..9d32b2becc0 100644 --- a/src/coreComponents/fieldSpecification/TractionBoundaryCondition.hpp +++ b/src/coreComponents/fieldSpecification/TractionBoundaryCondition.hpp @@ -23,6 +23,7 @@ #include "FieldSpecificationBase.hpp" #include "mesh/FaceManager.hpp" +#include "mesh/NodeManager.hpp" namespace geos { @@ -71,6 +72,7 @@ class TractionBoundaryCondition : public FieldSpecificationBase * @param blockLocalDofNumber Array of block local DOF numbers for the displacement. * @param dofRankOffset The rank offset for the DOF. * @param faceManager Reference to the face manager (Tractions are applied on faces) + * @param nodePositions The reference position of all nodes (used for proper FEM integration on faces) * @param targetSet The set of faces to apply the BC to. * @param localRhs The RHS of the system to add contributions to. */ @@ -78,6 +80,7 @@ class TractionBoundaryCondition : public FieldSpecificationBase arrayView1d< globalIndex const > const blockLocalDofNumber, globalIndex const dofRankOffset, FaceManager const & faceManager, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePositions, SortedArrayView< localIndex const > const & targetSet, arrayView1d< real64 > const & localRhs ) const; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index dcf51ae7c19..525be734d47 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -499,6 +499,18 @@ real64 SolidMechanicsLagrangianFEM::solverStep( real64 const & time_n, { int const maxNumResolves = m_maxNumResolves; int globallyFractured = 0; + + // Check if mesh was modified by an external event (e.g., SurfaceGen) before this solver step + // This ensures setupSystem is called before implicitStepSetup when topology changes occur + Timestamp const initialMeshModificationTimestamp = getMeshModificationTimestamp( domain ); + MeshLevel & meshLevel = domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ); + if( initialMeshModificationTimestamp > getSystemSetupTimestamp() || meshLevel.getModificationTimestamp() > getSystemSetupTimestamp() ) + { + m_dofManager.clear(); + setupSystem( domain, m_dofManager, m_localMatrix, m_rhs, m_solution, true ); + setSystemSetupTimestamp( LvArray::math::max( initialMeshModificationTimestamp, meshLevel.getModificationTimestamp() ) ); + } + implicitStepSetup( time_n, dt, domain ); for( int solveIter=0; solveIter getSystemSetupTimestamp() || globallyFractured ) + if( meshModificationTimestamp > getSystemSetupTimestamp() || meshLevel.getModificationTimestamp() > getSystemSetupTimestamp() || globallyFractured ) { setupSystem( domain, m_dofManager, m_localMatrix, m_rhs, m_solution ); - setSystemSetupTimestamp( meshModificationTimestamp ); + setSystemSetupTimestamp( LvArray::math::max( meshModificationTimestamp, meshLevel.getModificationTimestamp() ) ); } dtReturn = nonlinearImplicitStep( time_n, @@ -798,6 +810,7 @@ void SolidMechanicsLagrangianFEM::applyTractionBC( real64 const time, blockLocalDofNumber, dofRankOffset, faceManager, + nodeManager.referencePosition(), targetSet, localRhs ); } );