Skip to content

Commit

Permalink
Further trials to fix the parallel numbering issue with mesh split
Browse files Browse the repository at this point in the history
  • Loading branch information
@raback
raback committed Jan 31, 2025
1 parent fc663ee commit 29022a3
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 46 deletions.
128 changes: 83 additions & 45 deletions fem/src/MeshUtils.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2498,10 +2498,12 @@ FUNCTION LoadMesh2( Model, MeshDirPar, MeshNamePar,&
RETURN
END IF
END IF

! Prepare the mesh for next steps.
! For example, create non-nodal mesh structures, periodic projectors etc.
CALL PrepareMesh(Model,Mesh,Parallel,Def_Dofs,mySolver)
i = ListGetInteger( VList, 'Mesh Levels', GotIt )
CALL PrepareMesh(Model,Mesh,Parallel,Def_Dofs,mySolver,InitOnly=(i>1))

CALL Info(Caller,'Preparing mesh done',Level=8)


Expand Down Expand Up @@ -2934,16 +2936,18 @@ END FUNCTION LoadMesh2
!> Enlarge the coordinate vectors for p-elements.
!> Generate static projector for periodic BCS.
!-------------------------------------------------------------------
SUBROUTINE PrepareMesh( Model, Mesh, Parallel, Def_Dofs, mySolver )

SUBROUTINE PrepareMesh( Model, Mesh, Parallel, Def_Dofs, mySolver, InitOnly)
TYPE(Model_t) :: Model
TYPE(Mesh_t), POINTER :: Mesh
LOGICAL :: Parallel
INTEGER, OPTIONAL :: Def_Dofs(:,:), mySolver
TYPE(ValueList_t), POINTER :: Vlist
LOGICAL, OPTIONAL :: InitOnly

LOGICAL :: Found
CHARACTER(*),PARAMETER :: Caller='PrepareMesh'


IF( PRESENT( mySolver ) ) THEN
VList => Model % Solvers(mySolver) % Values
ELSE
Expand All @@ -2960,12 +2964,24 @@ SUBROUTINE PrepareMesh( Model, Mesh, Parallel, Def_Dofs, mySolver )
IF( ListGetLogical( Vlist,'Constant Stencil', Found ) ) THEN
CALL SetEqualElementIndeces( Mesh )
END IF

! CALL CreateIntersectionBCs(Model,Mesh)

! We know that not all structures are needed on this level.
! So we can do just partial initializations and leave early.
IF(PRESENT(InitOnly)) THEN
IF(InitOnly) THEN
CALL NonNodalElements()
!CALL EnlargeCoordinates( Mesh )
CALL Info(Caller,'Skipping rest of mesh preparation!',Found )
RETURN
END IF
END IF


IF( ListGetLogical( Vlist,'Increase Element Order',Found ) ) THEN
! We need to follow the boundary also for the new nodes of the quadratic mesh.
CALL FollowCurvedBoundary( Model, Mesh, .FALSE. )
CALL EnlargeCoordinates( Mesh )
CALL FollowCurvedBoundary( Model, Mesh, .FALSE. )
CALL IncreaseElementOrder( Model, Mesh )
END IF

Expand All @@ -2981,7 +2997,7 @@ SUBROUTINE PrepareMesh( Model, Mesh, Parallel, Def_Dofs, mySolver )
CALL EnlargeCoordinates( Mesh )

CALL FollowCurvedBoundary( Model, Mesh, .FALSE. )

CALL GeneratePeriodicProjectors( Model, Mesh )

IF( ListGetLogical( Vlist,'Inspect Quadratic Mesh', Found ) ) THEN
Expand Down Expand Up @@ -3171,8 +3187,8 @@ SUBROUTINE NonNodalElements()
!------------------------------------------------------------
body_id0 = -1; FoundDef=.FALSE.; FoundEq=.FALSE.
ElementDef = ' '


!
! Check whether face DOFs have been generated by "-quad_face b: ..." or
! "-tri_face b: ..."
!
Expand Down Expand Up @@ -3528,10 +3544,14 @@ SUBROUTINE ParallelNonNodalElements()
END IF

! Create parallel numbering of faces
CALL ResetTimer('SParFaceNumbering')
CALL SParFaceNumbering(Mesh, .TRUE. )
CALL CheckTimer('SParFaceNumbering',Level=7,Delete=.TRUE.)

! Create parallel numbering for edges
CALL ResetTimer('SParEdgeNumbering')
CALL SParEdgeNumbering(Mesh, .TRUE.)
CALL CheckTimer('SParEdgeNumbering',Level=7,Delete=.TRUE.)

! There are mainly implemented for parallel debugging.
! The whole sequence is only activated when "Max Output Level >= 10".
Expand Down Expand Up @@ -14845,7 +14865,7 @@ END SUBROUTINE PolynomBoundaryFit
SUBROUTINE FollowCurvedBoundary(Model, Mesh, SetP )
TYPE(Model_t) :: Model
TYPE(Mesh_t), POINTER :: Mesh
LOGICAL, OPTIONAL :: SetP
LOGICAL :: SetP

LOGICAL :: Found
REAL(KIND=dp) :: FitParams(12)
Expand Down Expand Up @@ -14974,7 +14994,7 @@ SUBROUTINE SetCurvedBoundary()
!------------------------------------------------------------------------------
REAL(KIND=dp) :: R, Rminor, r1, rat, f, gradf(3)
REAL(KIND=dp) :: Nrm(3), Tngt1(3), Tngt2(3), Orig(3), Coord(3), NtCoord(3), PlaneCoord(3)
INTEGER :: i,j,k,l,t,n
INTEGER :: i,j,k,l,t,n,elem
LOGICAL, POINTER :: DoneNode(:)
TYPE(Element_t), POINTER :: Element
LOGICAL :: Parallel
Expand Down Expand Up @@ -15010,14 +15030,16 @@ SUBROUTINE SetCurvedBoundary()
END IF

Parallel = ( ParEnv % PEs > 1 .AND. .NOT. Mesh % SingleMesh )

PRINT *,'SetP:',SetP

IF(.NOT. SetP) THEN
ALLOCATE( DoneNode(Mesh % NumberOfNodes))
DoneNode = .FALSE.

DO t=Mesh % NumberOfBulkElements+1, &
DO elem=Mesh % NumberOfBulkElements+1, &
Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
Element => Mesh % Elements(t)
Element => Mesh % Elements(elem)
IF ( Element % BoundaryInfo % Constraint &
== Model % BCs(bc_ind) % Tag ) THEN
n = Element % TYPE % NumberOfNodes
Expand Down Expand Up @@ -15084,9 +15106,9 @@ SUBROUTINE SetCurvedBoundary()
END IF

IF( SetP ) THEN
DO t=Mesh % NumberOfBulkElements+1, &
DO elem=Mesh % NumberOfBulkElements+1, &
Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
Element => Mesh % Elements(t)
Element => Mesh % Elements(elem)
IF ( Element % BoundaryInfo % Constraint &
/= Model % BCs(bc_ind) % Tag ) CYCLE
n = Element % TYPE % NumberOfNodes
Expand Down Expand Up @@ -15168,7 +15190,7 @@ SUBROUTINE SetCurvedBoundary()
Coord = Coord + Orig
! Solve for desired coordinate displacement rather than absolute coordinate value
Coord = Coord - Coord0

! Create equation involving mass matrix that solves for the coordinates at the p-dofs
DO q=1,nd
MASS(1:nd,q) = MASS(1:nd,q) + Weight * Basis(1:nd) * Basis(q)
Expand All @@ -15187,8 +15209,17 @@ SUBROUTINE SetCurvedBoundary()
END DO

CALL LUdecomp(MASS,nd,pivot,Erroneous)
IF (Erroneous) CALL Fatal('SetCurvedBoundary', 'LU-decomposition fails')

IF (Erroneous) THEN
PRINT *,'Element:',elem,ip % n,nd,n,dim
DO i=1,dim
PRINT *,'FORCE',i,FORCE(i,1:nd)
END DO
DO i=n+1,nd
PRINT *,'MASS',i,MASS(i,1:nd)
END DO
CALL Fatal('SetCurvedBoundary', 'LU-decomposition fails')
END IF

DO i=1,dim
x(1:nd) = FORCE(i,1:nd)
CALL LUSolve(nd,MASS,x,pivot)
Expand Down Expand Up @@ -21934,10 +21965,9 @@ FUNCTION SplitMeshEqual(Mesh,h) RESULT( NewMesh )
FaceCnt = FaceCnt+1
FacePerm(i) = NodeCnt
END IF
END DO

WRITE( Message, * ) 'Added nodes in the center of faces : ', FaceCnt
CALL Info( 'SplitMeshEqual', Message, Level=10 )
END DO
IF(FaceCnt>0) CALL Info( 'SplitMeshEqual','Added '//I2S(FaceCnt)//' nodes in the center of faces',Level=10)

!
! For quads and bricks, count centerpoints:
! -----------------------------------------
Expand All @@ -21949,10 +21979,9 @@ FUNCTION SplitMeshEqual(Mesh,h) RESULT( NewMesh )
NodeCnt = NodeCnt + 1
NodeIt = NodeIt + 1
END SELECT
END DO

WRITE( Message, * ) 'Added nodes in the center of bulks : ', NodeIt
CALL Info( 'SplitMeshEqual', Message, Level=10 )
END DO
IF(NodeIt>0) CALL Info( 'SplitMeshEqual','Added '//I2S(NodeIt)//' nodes in the center of bulks',Level=10)

!
! new mesh nodecoordinate arrays:
! -------------------------------
Expand Down Expand Up @@ -22003,10 +22032,9 @@ FUNCTION SplitMeshEqual(Mesh,h) RESULT( NewMesh )
y(j) = SUM(v(Edge % NodeIndexes))/k
z(j) = SUM(w(Edge % NodeIndexes))/k
END IF
END DO

CALL Info('SplitMeshEqual','Added edge centers to the nodes list.', Level=10 )
!
END DO
CALL Info('SplitMeshEqual','Added edge centers to the nodes list.', Level=15 )

! add quad face centers for bricks and prisms(wedges):
! ----------------------------
j = Mesh % NumberOfNodes + Mesh % NumberOfEdges
Expand All @@ -22033,11 +22061,9 @@ FUNCTION SplitMeshEqual(Mesh,h) RESULT( NewMesh )
z(j) = SUM(w(Face % NodeIndexes))/k
END IF
END IF
END DO
END DO
CALL Info('SplitMeshEqual','Added face centers to the nodes list.', Level=15 )


CALL Info('SplitMeshEqual','Added face centers to the nodes list.', Level=10 )
!
! add centerpoint for quads & bricks:
! -----------------------------------
DO i=1,Mesh % NumberOfBulkElements
Expand Down Expand Up @@ -22088,7 +22114,9 @@ FUNCTION SplitMeshEqual(Mesh,h) RESULT( NewMesh )
END IF
END SELECT
END DO
!
CALL Info('SplitMeshEqual','Added quad and brick centers to the nodes list.', Level=15 )


! Update new mesh node count:
! ---------------------------
NewMesh % NumberOfEdges = 0
Expand Down Expand Up @@ -23508,7 +23536,7 @@ FUNCTION SplitMeshEqual(Mesh,h) RESULT( NewMesh )
! Our boundary may be a circle, cylider or sphere surface.
! Honor those shapes when splitting the mesh!
CALL FollowCurvedBoundary( CurrentModel, NewMesh, .FALSE. )


!call writemeshtodisk( NewMesh, "." )
!stop
Expand All @@ -23529,7 +23557,7 @@ SUBROUTINE UpdateParallelMesh( Mesh, NewMesh )
LOGICAL :: Found
!------------------------------------------------------------------------------

!

! Update mesh interfaces for parallel execution.
! ==============================================
!
Expand Down Expand Up @@ -23781,11 +23809,15 @@ SUBROUTINE UpdateParallelMesh( Mesh, NewMesh )
END IF

DO i=1,Mesh % NumberOfNodes
CALL AllocateVector( NewMesh % ParallelInfo % NeighbourList(i) % Neighbours, &
SIZE( Mesh % ParallelInfo % Neighbourlist(i) % Neighbours) )

NewMesh % ParallelInfo % NeighbourList(i) % Neighbours = &
Mesh % ParallelInfo % NeighbourList(i) % Neighbours
IF(.NOT. ASSOCIATED( Mesh % ParallelInfo % NeighbourList(i) % Neighbours ) ) THEN
CALL AllocateVector( NewMesh % ParallelInfo % NeighbourList(i) % Neighbours, 1 )
NewMesh % ParallelInfo % NeighbourList(i) % Neighbours = ParEnv % MyPe
ELSE
CALL AllocateVector( NewMesh % ParallelInfo % NeighbourList(i) % Neighbours, &
SIZE( Mesh % ParallelInfo % Neighbourlist(i) % Neighbours) )
NewMesh % ParallelInfo % NeighbourList(i) % Neighbours = &
Mesh % ParallelInfo % NeighbourList(i) % Neighbours
END IF
END DO
!
! Take care of the new mesh internal nodes.
Expand Down Expand Up @@ -23848,18 +23880,24 @@ SUBROUTINE UpdateParallelMesh( Mesh, NewMesh )
Reorder = [ (i, i=1,NewMesh % NumberOfNodes) ]

k = NewMesh % Nodes % NumberOfNodes - Mesh % Nodes % NumberOfNodes

CALL ResetTimer('ParallelGlobalNumbering')
CALL ParallelGlobalNumbering( NewMesh, Mesh, k, Reorder )
CALL CheckTimer('ParallelGlobalNumbering',Level=7,Delete=.TRUE.)


! Account for the reordering of the nodes:
! ----------------------------------------
DO i=1,NewMesh % NumberOfBulkElements + &
NewMesh % NumberOfBoundaryElements
NewMesh % Elements(i) % NodeIndexes = &
Reorder( NewMesh % Elements(i) % NodeIndexes )
NewMesh % NumberOfBoundaryElements
NewMesh % Elements(i) % NodeIndexes = &
Reorder( NewMesh % Elements(i) % NodeIndexes )
END DO

! DEALLOCATE( IntCnts, IntArray, Reorder )
! DEALLOCATE( Reorder )
! DEALLOCATE( Reorder )


!------------------------------------------------------------------------------
END SUBROUTINE UpdateParallelMesh
END FUNCTION SplitMeshEqual
Expand Down
2 changes: 1 addition & 1 deletion fem/src/ModelDescription.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2985,7 +2985,7 @@ FUNCTION LoadModel( ModelName,BoundariesOnly,numprocs,mype,MeshIndex) RESULT( Mo
IF( MeshLevels > 1 ) THEN
! This has been commented out, but is needed. There may be some issues in parallel still...
IF(ListGetLogical(Model % Simulation,'Prepare Mesh Before Split',GotIt) ) THEN
CALL PrepareMesh( Model, NewMesh, ParEnv % PEs > 1 )
CALL PrepareMesh( Model, NewMesh, ParEnv % PEs > 1, InitOnly = (i<MeshLevels) )
END IF
END IF

Expand Down

0 comments on commit 29022a3

Please sign in to comment.