From db1718a74bf0754c553d5bf2b15f3d6eb08e47f4 Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Wed, 10 Aug 2022 13:27:36 +0200 Subject: [PATCH] Improved initialization performance for rotational periodic BC by splitting the work among the compute node processors. Fixed bug due to missing tolerance, which now increases the bounding box size by 1 percent of the length in h and r during the iSide/jSide test. This hopefully results in less or no lost particles. --- .../NIG_DSMC/RotPeriodicBC/parameter.ini | 8 +- .../pre-hopr-mortar/hopr_mortar.ini | 13 +- .../NIG_DSMC/RotPeriodicBC/readme.md | 16 + src/io_hdf5/hdf5_output_elemdata.f90 | 4 +- src/mesh/mesh_vars.f90 | 7 +- .../boundary/particle_boundary_condition.f90 | 31 +- .../boundary/particle_boundary_init.f90 | 478 +++++++++++++----- .../boundary/particle_boundary_vars.f90 | 21 +- src/particles/particle_init.f90 | 2 +- 9 files changed, 419 insertions(+), 161 deletions(-) diff --git a/regressioncheck/NIG_DSMC/RotPeriodicBC/parameter.ini b/regressioncheck/NIG_DSMC/RotPeriodicBC/parameter.ini index bea622aa1..be9e737e2 100644 --- a/regressioncheck/NIG_DSMC/RotPeriodicBC/parameter.ini +++ b/regressioncheck/NIG_DSMC/RotPeriodicBC/parameter.ini @@ -6,9 +6,11 @@ IniExactFunc = 0 ! =============================================================================== ! ! DISCRETIZATION ! =============================================================================== ! -N = 1 ! Polynomial degree -NAnalyze = 1 ! Number of analyze points -NVisu = 0 +N = 1 ! Polynomial degree +NAnalyze = 1 ! Number of analyze points +NVisu = 1 +VisuParticles = T +!DisplayLostParticles = T ! =============================================================================== ! ! MESH ! =============================================================================== ! diff --git a/regressioncheck/NIG_DSMC/RotPeriodicBC/pre-hopr-mortar/hopr_mortar.ini b/regressioncheck/NIG_DSMC/RotPeriodicBC/pre-hopr-mortar/hopr_mortar.ini index 56e8cc2f5..cd3334b2d 100644 --- a/regressioncheck/NIG_DSMC/RotPeriodicBC/pre-hopr-mortar/hopr_mortar.ini +++ b/regressioncheck/NIG_DSMC/RotPeriodicBC/pre-hopr-mortar/hopr_mortar.ini @@ -1,14 +1,15 @@ -DEFVAR=(INT): i01 = 6 ! no. elems in left and right block -DEFVAR=(INT): i02 = 6 ! no. elems in upper block (should be twice the value of i01) +DEFVAR=(INT): i01 = 6!2 ! no. elems in left and right block +DEFVAR=(INT): i02 = 6!6!2 ! no. elems in upper block (should be twice the value of i01) -DEFVAR=(INT): ir1 = 25 ! no. elems in r for first ring +DEFVAR=(INT): ir1 = 25!3 ! no. elems in r for first ring DEFVAR=(REAL): r01 = 3.5 ! middle square dim DEFVAR=(REAL): r02 = 7.0 ! middle square dim DEFVAR=(REAL): s0 = 0.2857142857142857 ! middle square dim -DEFVAR=(INT): iz = 50 ! +DEFVAR=(INT): iz1 = 50!3 ! +DEFVAR=(INT): iz2 = 100!6 ! DEFVAR=(REAL): lz = 1 ! length of domain in z DEFVAR=(REAL): f1 = 1. ! stretching factor in first ring @@ -33,7 +34,7 @@ nZones = 2 ! number of boxes !right-lower (x+) Corner =(/r01 , 0. , -lz ,, r02 , 0. , -lz ,, r02 , r02 , -lz ,, r01 , r01 , -lz ,, r01 , 0. , lz ,, r02 , 0. , lz ,, r02 , r02 , lz ,, r01 , r01 , lz /) -nElems =(/ir1,i01,iz/) ! number of elements in each direction +nElems =(/ir1,i01,iz1/) ! number of elements in each direction BCIndex =(/1 , 5 , 4 , 0 , 3 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) ! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions @@ -42,7 +43,7 @@ factor =(/f1,1.,1./) ! stretching !right-upper (y+) Corner =(/0. , r01 , -lz ,, r01 , r01 , -lz ,, r02 , r02 , -lz ,, 0. , r02 , -lz ,, 0. , r01 , lz ,, r01 , r01 , lz ,, r02 , r02 , lz ,, 0. , r02 , lz /) -nElems =(/i02,50,100/) ! number of elements in each direction +nElems =(/i02,iz1,iz2/) ! number of elements in each direction BCIndex =(/1 , 3 , 0 , 4 , 6 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) ! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions elemtype =108 ! element type (108: Hexahedral) diff --git a/regressioncheck/NIG_DSMC/RotPeriodicBC/readme.md b/regressioncheck/NIG_DSMC/RotPeriodicBC/readme.md index a70351413..659ce7d1e 100644 --- a/regressioncheck/NIG_DSMC/RotPeriodicBC/readme.md +++ b/regressioncheck/NIG_DSMC/RotPeriodicBC/readme.md @@ -7,4 +7,20 @@ finding a corresponding rotationally periodic side) - Debugging information should be written to TestRotatingWall_LostRotPeriodicSides_000.00000000000000000.h5 file (cannot be tested automatically unfortunately) when setting CalcMeshInfo=T +- for debugging, use 2 procs and smaller mesh via + + DEFVAR=(INT): i01 = 2 ! no. elems in left and right block + DEFVAR=(INT): i02 = 2 ! no. elems in upper block (should be twice the value of i01) + + DEFVAR=(INT): ir1 = 3 ! no. elems in r for first ring + DEFVAR=(REAL): r01 = 3.5 ! middle square dim + DEFVAR=(REAL): r02 = 7.0 ! middle square dim + DEFVAR=(REAL): s0 = 0.2857142857142857 ! middle square dim + + + DEFVAR=(INT): iz1 = 3 + DEFVAR=(INT): iz2 = 6 + DEFVAR=(REAL): lz = 1 ! length of domain in z + + DEFVAR=(REAL): f1 = 1. ! stretching factor in first ring diff --git a/src/io_hdf5/hdf5_output_elemdata.f90 b/src/io_hdf5/hdf5_output_elemdata.f90 index dfa810a05..4785200ba 100644 --- a/src/io_hdf5/hdf5_output_elemdata.f90 +++ b/src/io_hdf5/hdf5_output_elemdata.f90 @@ -28,10 +28,10 @@ MODULE MOD_HDF5_Output_ElemData #if USE_MPI PUBLIC :: WriteMyInvisibleRankToHDF5 +#endif /*USE_MPI*/ #if defined(PARTICLES) PUBLIC :: WriteLostRotPeriodicSidesToHDF5 #endif /*defined(PARTICLES)*/ -#endif /*USE_MPI*/ PUBLIC :: WriteAdditionalElemData !=================================================================================================================================== @@ -191,6 +191,7 @@ SUBROUTINE WriteMyInvisibleRankToHDF5() WRITE(UNIT_stdOut,'(a)',ADVANCE='YES')'DONE' #endif END SUBROUTINE WriteMyInvisibleRankToHDF5 +#endif /*USE_MPI*/ #if defined(PARTICLES) @@ -247,7 +248,6 @@ SUBROUTINE WriteLostRotPeriodicSidesToHDF5() #endif END SUBROUTINE WriteLostRotPeriodicSidesToHDF5 #endif /*defined(PARTICLES)*/ -#endif /*USE_MPI*/ END MODULE MOD_HDF5_Output_ElemData diff --git a/src/mesh/mesh_vars.f90 b/src/mesh/mesh_vars.f90 index afaa1ce92..7a9ab2604 100644 --- a/src/mesh/mesh_vars.f90 +++ b/src/mesh/mesh_vars.f90 @@ -174,10 +174,11 @@ MODULE MOD_Mesh_Vars !< [1:4] (mortar) neighbors !< [1:6] local sides !< [OffSetElem+1:OffsetElem+PP_nElems] -INTEGER(KIND=8),ALLOCATABLE :: ElemGlobalID(:) !< global element id of each element -INTEGER(KIND=8),ALLOCATABLE :: myInvisibleRank(:) !< global proc ID which the current proc cannot see (particle +INTEGER(KIND=8),ALLOCATABLE :: ElemGlobalID(:) !< global element id of each element +INTEGER(KIND=8),ALLOCATABLE :: myInvisibleRank(:) !< global proc ID which the current proc cannot see (particle !< communication) -INTEGER(KIND=8),ALLOCATABLE :: LostRotPeriodicSides(:) !< Number of lost sides during rotational periodic search +INTEGER(KIND=8),ALLOCATABLE :: LostRotPeriodicSides(:) !< Number of lost sides during rotational periodic search +!LOGICAL :: RotPeriodicReBuild !< Force re-building of mapping (might already exist) !----------------------------------------------------------------------------------------------------------------------------------- CHARACTER(LEN=255),ALLOCATABLE :: BoundaryName(:) CHARACTER(LEN=255) :: MeshFile !< name of hdf5 meshfile (write with ending .h5!) diff --git a/src/particles/boundary/particle_boundary_condition.f90 b/src/particles/boundary/particle_boundary_condition.f90 index 85ba4411d..3c88bd280 100644 --- a/src/particles/boundary/particle_boundary_condition.f90 +++ b/src/particles/boundary/particle_boundary_condition.f90 @@ -398,7 +398,7 @@ SUBROUTINE RotPeriodicBC(PartID,SideID,ElemID) USE MOD_Particle_Mesh_Tools ,ONLY: ParticleInsideQuad3D USE MOD_part_operations ,ONLY: RemoveParticle USE MOD_part_tools ,ONLY: StoreLostParticleProperties -USE MOD_Particle_Tracking_Vars ,ONLY: NbrOfLostParticles, TrackInfo, CountNbrOfLostParts +USE MOD_Particle_Tracking_Vars ,ONLY: NbrOfLostParticles, TrackInfo, CountNbrOfLostParts,DisplayLostParticles USE MOD_DSMC_Vars ,ONLY: DSMC, AmbipolElecVelo #ifdef CODE_ANALYZE USE MOD_Particle_Tracking_Vars ,ONLY: PartOut,MPIRankOut @@ -413,7 +413,7 @@ SUBROUTINE RotPeriodicBC(PartID,SideID,ElemID) INTEGER,INTENT(INOUT),OPTIONAL :: ElemID !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: SideID2, ElemID2, iNeigh, RotSideID +INTEGER :: iNeigh, RotSideID, GlobalElemID REAL :: adaptTimeStep LOGICAL :: FoundInElem REAL :: LastPartPos_old(1:3),Velo_old(1:3), Velo_oldAmbi(1:3) @@ -480,7 +480,7 @@ SUBROUTINE RotPeriodicBC(PartID,SideID,ElemID) END IF END SELECT END ASSOCIATE -! (2) update particle positon after periodic BC +! (2) update particle position after periodic BC PartState(1:3,PartID) = LastPartPos(1:3,PartID) + (1.0 - TrackInfo%alpha/TrackInfo%lengthPartTrajectory) * dt*RKdtFrac & * PartState(4:6,PartID) * adaptTimeStep ! compute moved particle || rest of movement @@ -496,26 +496,31 @@ SUBROUTINE RotPeriodicBC(PartID,SideID,ElemID) IF(PartID.EQ.PARTOUT)THEN IPWRITE(UNIT_stdout,'(I0,A)') ' RotPeriodicBC: ' IPWRITE(UNIT_stdout,'(I0,A,3(1X,G0))') ' ParticlePosition-pp: ',PartState(1:3,PartID) - IPWRITE(UNIT_stdout,'(I0,A,3(1X,G0))') ' LastPartPo-pp: ',LastPartPos(1:3,PartID) + IPWRITE(UNIT_stdout,'(I0,A,3(1X,G0))') ' LastPartPosition-pp: ',LastPartPos(1:3,PartID) END IF END IF #endif /*CODE_ANALYZE*/ ! (3) move particle from old element to new element -RotSideID = SurfSide2RotPeriodicSide((GlobalSide2SurfSide(SURF_SIDEID,SideID))) +RotSideID = SurfSide2RotPeriodicSide(GlobalSide2SurfSide(SURF_SIDEID,SideID)) DO iNeigh=1,NumRotPeriodicNeigh(RotSideID) - SideID2 = RotPeriodicSideMapping(RotSideID,iNeigh) - IF(SideID2.EQ.-1) THEN - CALL abort(__STAMP__,' ERROR: Halo-rot-periodic side has no corresponding side.') - END IF - ElemID2 = SideInfo_Shared(SIDE_ELEMID,SideID2) - ! find rotational periodic SideID2 through localization in all potentional rotational periodic sides - CALL ParticleInsideQuad3D(LastPartPos(1:3,PartID),ElemID2,FoundInElem) + ! find rotational periodic elem through localization in all potential rotational periodic sides + GlobalElemID = RotPeriodicSideMapping(RotSideID,iNeigh) + IF(GlobalElemID.EQ.-1) CALL abort(__STAMP__,' ERROR: Halo-rot-periodic side has no corresponding element.') + CALL ParticleInsideQuad3D(LastPartPos(1:3,PartID),GlobalElemID,FoundInElem) IF(FoundInElem) THEN - ElemID = ElemID2 + ElemID = GlobalElemID EXIT END IF END DO IF(.NOT.FoundInElem) THEN + ! Particle appears to have not crossed any of the checked sides. Deleted! + IF(DisplayLostParticles)THEN + IPWRITE(*,*) 'Error in Particle TriaTracking! Particle Number',PartID,'lost. Element:', ElemID,'(species:',PartSpecies(PartID),')' + IPWRITE(*,*) 'LastPos: ', LastPartPos(1:3,PartID) + IPWRITE(*,*) 'Pos: ', PartState(1:3,PartID) + IPWRITE(*,*) 'Velo: ', PartState(4:6,PartID) + IPWRITE(*,*) 'Particle deleted!' + END IF ! DisplayLostParticles IF(CountNbrOfLostParts)THEN CALL StoreLostParticleProperties(PartID,ElemID) NbrOfLostParticles=NbrOfLostParticles+1 diff --git a/src/particles/boundary/particle_boundary_init.f90 b/src/particles/boundary/particle_boundary_init.f90 index b86b403dd..d8c42e4fc 100644 --- a/src/particles/boundary/particle_boundary_init.f90 +++ b/src/particles/boundary/particle_boundary_init.f90 @@ -116,6 +116,7 @@ SUBROUTINE DefineParametersParticleBoundary() CALL prms%CreateIntOption( 'Part-Boundary[$]-RotPeriodicDir' & , 'Angular degree of rotation periodicity in [deg].' & , '0', numberedmulti=.TRUE.) +!CALL prms%CreateLogicalOption( 'Part-RotPeriodicReBuild', 'Force re-creation of rotational periodic mapping (which might already exist in the mesh file).', '.FALSE.') CALL prms%CreateRealOption( 'Part-Boundary[$]-WallTemp2' & , 'Second wall temperature (in [K]) of reflective particle boundary for a temperature gradient.' & , '0.', numberedmulti=.TRUE.) @@ -528,85 +529,219 @@ SUBROUTINE InitializeVariablesPartBoundary() END SUBROUTINE InitializeVariablesPartBoundary +!=================================================================================================================================== +!> Read mapping for rotational periodicity from mesh file. If it does not yet exist, build the mapping and store in mesh.h5 for +!> faster initialization later on. +!=================================================================================================================================== +SUBROUTINE InitParticleBoundaryRotPeriodic() +! MODULES +USE MOD_Globals +USE MOD_HDF5_Input +USE MOD_HDF5_output ,ONLY: WriteAttributeToHDF5 +USE MOD_Mesh_Vars ,ONLY: MeshFile!,RotPeriodicReBuild +USE MOD_Particle_MPI_Vars ,ONLY: halo_eps_velo +USE MOD_TimeDisc_Vars ,ONLY: ManualTimeStep +USE MOD_ReadInTools ,ONLY: GETLOGICAL +USE MOD_Globals_Vars ,ONLY: ProjectName +#if USE_LOADBALANCE +USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance +#endif /*USE_LOADBALANCE*/ +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------! +! INPUT / OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +LOGICAL :: DatasetFound +CHARACTER(LEN=64) :: DatasetName,hilf +REAL :: StartT,EndT +INTEGER :: notMappedTotal +!=================================================================================================================================== + +LBWRITE(UNIT_stdOut,'(A)') ' INIT ROTATIONAL PERIODIC BOUNDARY...' +!RotPeriodicReBuild = GETLOGICAL('Part-RotPeriodicReBuild') +GETTIME(StartT) + +DatasetFound = .FALSE. +notMappedTotal = 0 +WRITE(UNIT=hilf,FMT='(ES10.4)') halo_eps_velo +DatasetName = 'RotPeriodicMap-v'//TRIM(hilf) +WRITE(UNIT=hilf,FMT='(ES10.4)') ManualTimeStep +DatasetName = TRIM(DatasetName)//'-dt'//TRIM(hilf) +CALL OpenDataFile(MeshFile,create=.FALSE.,single=.FALSE.,readOnly=.FALSE.,communicatorOpt=MPI_COMM_WORLD) +CALL DatasetExists(File_ID,TRIM(DatasetName),DatasetFound) +CALL CloseDataFile() + +IF(DatasetFound)THEN!.AND.(.NOT.RotPeriodicReBuild))THEN + ! Read mapping from mesh file + LBWRITE(Unit_StdOut,'(A)')" Reading ["//TRIM(DatasetName)//"] from mesh file ["//TRIM(MeshFile)//"]" + ! This feature is not implemented yet. If more performance is required, e.g., during load balancing, store the mapping in the + ! mesh file and simply read it during initialization without the need for searching the connections again + !CALL ReadArray(TRIM(DatasetName),2,(/x,x,x/),0_IK,1,RealArray=mapping) + CALL ReadAttribute(File_ID,'notMappedTotal',1,IntScalar=notMappedTotal) + CALL CloseDataFile() +ELSE + ! Create mapping and store in mesh file + CALL BuildParticleBoundaryRotPeriodic(notMappedTotal) + + ! Store number of unmapped elements in mesh file + IF(MPIRoot)THEN + CALL OpenDataFile(MeshFile,create=.FALSE.,single=.TRUE.,readOnly=.FALSE.) + CALL WriteAttributeToHDF5(File_ID,'notMappedTotal',1,IntegerScalar=notMappedTotal,Overwrite=.TRUE.) + CALL CloseDataFile() + END IF ! MPIRoot +END IF ! DatasetFound + + +IF(notMappedTotal.GT.0)THEN + LBWRITE(Unit_StdOut,'(A,I0,A)') ' | Warning: Found ',notMappedTotal,' rot periodic sides that did not find a corresponding side.' + LBWRITE(Unit_StdOut,'(A)')" | The halo region (halo flag 2) merely reaches a rot periodic BC side but not any further." + LBWRITE(Unit_StdOut,'(A)')" | See ElemData container 'LostRotPeriodicSides' for more information on where sides were unmatched." + LBWRITE(Unit_StdOut,'(A)')" | This information is written to "//TRIM(ProjectName)//"_LostRotPeriodicSides.h5 (only when CalcMeshInfo=T)" + !LBWRITE(Unit_StdOut,'(A)')" | If the file does not exist, it can be re-created with RotPeriodicReBuild=T" +END IF ! notMappedTotal.GT.0 + + +GETTIME(EndT) +LBWRITE(UNIT_stdOut,'(A,F0.3,A)') ' INIT ROTATIONAL PERIODIC BOUNDARY DONE! [',EndT-StartT,'s]' +LBWRITE(UNIT_StdOut,'(132("-"))') + +END SUBROUTINE InitParticleBoundaryRotPeriodic + + !=================================================================================================================================== !> Build Mapping for rotational periodicity: RotPeriodicSide -> SideID2 (Side on corresponding BC). !> In RotPeriodicBC (particle_boundary_condition.f90): SideID -> SurfSideID -> RotPeriodicSide !> RotPeriodicSide -> SideID2 !> (1) counting rotational periodic sides and build mapping from SurfSideID -> RotPeriodicSide -!> (2) find Side on corresponding BC and build mapping RotPeriodicSide -> SideID2 (and vice versa) +!> (2) Build bounding boxes (in 2D reference system) for all nRotPeriodicSides +!> (3) find Side on corresponding BC and build mapping RotPeriodicSide -> SideID2 (and vice versa) !> counting potential rotational periodic sides (for not conform meshes) -!> (3) reallocate array due to number of potential rotational periodic sides +!> (4) reallocate array due to number of potential rotational periodic sides !=================================================================================================================================== -SUBROUTINE InitParticleBoundaryRotPeriodic() +SUBROUTINE BuildParticleBoundaryRotPeriodic(notMappedTotal) ! MODULES USE MOD_Globals -USE MOD_Particle_Boundary_Vars ,ONLY: RotPeriodicSide2GlobalSide,nComputeNodeSurfTotalSides,SurfSide2GlobalSide,PartBound +USE MOD_Particle_Boundary_Vars ,ONLY: nComputeNodeSurfTotalSides,SurfSide2GlobalSide,PartBound USE MOD_Particle_Boundary_Vars ,ONLY: RotPeriodicSideMapping, NumRotPeriodicNeigh, SurfSide2RotPeriodicSide USE MOD_Particle_Mesh_Vars ,ONLY: SideInfo_Shared, NodeCoords_Shared, ElemSideNodeID_Shared, GEO, ElemInfo_Shared USE MOD_Mesh_Tools ,ONLY: GetCNElemID -#if USE_MPI -USE MOD_Analyze_Vars ,ONLY: CalcMeshInfo +USE MOD_Particle_Mesh_Vars ,ONLY: SideInfo_Shared USE MOD_Mesh_Vars ,ONLY: LostRotPeriodicSides,nElems +USE MOD_Analyze_Vars ,ONLY: CalcMeshInfo USE MOD_IO_HDF5 ,ONLY: AddToElemData,ElementOut USE MOD_HDF5_Output_ElemData ,ONLY: WriteLostRotPeriodicSidesToHDF5 -USE MOD_Globals_Vars ,ONLY: ProjectName +#if USE_MPI +USE MOD_Particle_Boundary_Vars ,ONLY: SurfSide2RotPeriodicSide_Shared,SurfSide2RotPeriodicSide_Shared_Win +USE MOD_Particle_Boundary_Vars ,ONLY: Rot2Glob_temp_Shared,Rot2Glob_temp_Shared_Win +USE MOD_Particle_Boundary_Vars ,ONLY: RotPeriodicSideMapping_temp_Shared,RotPeriodicSideMapping_temp_Shared_Win +USE MOD_Particle_Boundary_Vars ,ONLY: RotPeriodicSideMapping_Shared,RotPeriodicSideMapping_Shared_Win +USE MOD_Particle_Boundary_Vars ,ONLY: BoundingBox_Shared,BoundingBox_Shared_Win +USE MOD_Particle_Boundary_Vars ,ONLY: NumRotPeriodicNeigh_Shared,NumRotPeriodicNeigh_Shared_Win +USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED +USE MOD_MPI_Shared_Vars ,ONLY: myComputeNodeRank,nComputeNodeProcessors +USE MOD_MPI_Shared +#else +USE MOD_Particle_Boundary_Vars ,ONLY: nSurfTotalSides #endif /*USE_MPI*/ -#if USE_LOADBALANCE -USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance -#endif /*USE_LOADBALANCE*/ !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------! ! OUTPUT VARIABLES +INTEGER,INTENT(OUT) :: notMappedTotal !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: iSide, jSide, nRotPeriodicSides, SideID,SideID2, MaxNumRotPeriodicNeigh, iNode, jNode, iNeigh -INTEGER :: NodeID, CNElemID, CNElemID2, LocSideID, k, l, m, NodeID2, LocSideID2 -INTEGER,ALLOCATABLE :: Rot2Glob_temp(:) -INTEGER,ALLOCATABLE :: RotPeriodicSideMapping_temp(:,:) -LOGICAL :: isMapped,SideIsMapped -REAL :: iNodeVec(1:3), jNodeVec(1:3) -REAL :: iNodeR, iNodeH, jNodeR, jNodeH, Node2Rmin, Node2Rmax, Node2Hmin, Node2Hmax -INTEGER,PARAMETER :: NbrOfRotConnections=1000 -INTEGER :: notMapped -#if USE_MPI -INTEGER :: notMappedTotal -REAL :: StartT,EndT -#endif /*USE_MPI*/ +INTEGER :: iSide, jSide, nRotPeriodicSides, SideID,SideID2, MaxNumRotPeriodicNeigh, iNode, iNeigh +INTEGER :: NodeID, CNElemID, LocSideID, k, l, m +LOGICAL :: mySide +REAL :: iNodeVec(1:3), jNodeVec(1:3) +REAL :: iNodeR, iNodeH, jNodeR, jNodeH, Node2Rmin, Node2Rmax, Node2Hmin, Node2Hmax, dh, dr +INTEGER,PARAMETER :: NbrOfRotConnections=1000 +INTEGER :: notMapped,GlobalElemID +INTEGER :: firstSide,lastSide,offsetSide +INTEGER,ALLOCPOINT,DIMENSION(:) :: Rot2Glob_temp +INTEGER,ALLOCPOINT,DIMENSION(:,:) :: RotPeriodicSideMapping_temp +REAL,ALLOCPOINT,DIMENSION(:,:) :: BoundingBox !=================================================================================================================================== -LBWRITE(UNIT_stdOut,'(A)') ' INIT ROTATIONAL PERIODIC BOUNDARY...' -#if USE_MPI -IF(MPIRoot) StartT=MPI_WTIME() -#endif /*USE_MPI*/ +nRotPeriodicSides = 0 +offsetSide = 0 -ALLOCATE(Rot2Glob_temp(nComputeNodeSurfTotalSides)) -ALLOCATE(SurfSide2RotPeriodicSide(nComputeNodeSurfTotalSides)) +! Surf sides are shared, array calculation can be distributed +#if USE_MPI +CALL Allocate_Shared((/nComputeNodeSurfTotalSides/) , SurfSide2RotPeriodicSide_Shared_Win , SurfSide2RotPeriodicSide_Shared) +CALL MPI_WIN_LOCK_ALL(0 , SurfSide2RotPeriodicSide_Shared_Win , IERROR) +CALL Allocate_Shared((/nComputeNodeSurfTotalSides/) , Rot2Glob_temp_Shared_Win , Rot2Glob_temp_Shared) +CALL MPI_WIN_LOCK_ALL(0 , Rot2Glob_temp_Shared_Win , IERROR) +SurfSide2RotPeriodicSide => SurfSide2RotPeriodicSide_Shared +Rot2Glob_temp => Rot2Glob_temp_Shared +IF (myComputeNodeRank.EQ.0) SurfSide2RotPeriodicSide = -1. +IF (myComputeNodeRank.EQ.0) Rot2Glob_temp = -1. +CALL BARRIER_AND_SYNC(SurfSide2RotPeriodicSide_Shared_Win , MPI_COMM_SHARED) +CALL BARRIER_AND_SYNC(Rot2Glob_temp_Shared_Win , MPI_COMM_SHARED) +firstSide = INT(REAL( myComputeNodeRank *nComputeNodeSurfTotalSides)/REAL(nComputeNodeProcessors))+1 +lastSide = INT(REAL((myComputeNodeRank+1)*nComputeNodeSurfTotalSides)/REAL(nComputeNodeProcessors)) +#else +firstSide = 1 +lastSide = nSurfTotalSides +ALLOCATE(SurfSide2RotPeriodicSide(firstSide:lastSide)) SurfSide2RotPeriodicSide(:) = -1 -nRotPeriodicSides=0 +ALLOCATE(Rot2Glob_temp(firstSide:lastSide)) +#endif /*USE_MPI*/ ! (1) Count rotational periodic sides and build mapping from SurfSideID -> RotPeriodicSide -DO iSide=1, nComputeNodeSurfTotalSides +#if USE_MPI +! Only when more than 1 core per node +IF(nComputeNodeProcessors.GT.1)THEN + ! Loop once and get offset + DO iSide=firstSide, lastSide + SideID = SurfSide2GlobalSide(SURF_SIDEID,iSide) + IF(PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).EQ.PartBound%RotPeriodicBC) & + nRotPeriodicSides = nRotPeriodicSides + 1 + END DO + + CALL MPI_EXSCAN(nRotPeriodicSides,offsetSide,1,MPI_INTEGER,MPI_SUM,MPI_COMM_SHARED,iError) + nRotPeriodicSides = offsetSide +END IF ! nComputeNodeProcessors.GT.1 +#endif /*USE_MPI*/ + +DO iSide=firstSide, lastSide SideID = SurfSide2GlobalSide(SURF_SIDEID,iSide) IF(PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).EQ.PartBound%RotPeriodicBC) THEN nRotPeriodicSides = nRotPeriodicSides + 1 Rot2Glob_temp(nRotPeriodicSides) = SideID - SurfSide2RotPeriodicSide(iSide) = nRotPeriodicSides + SurfSide2RotPeriodicSide(iSide) = nRotPeriodicSides ! Store RotSideID END IF END DO -ALLOCATE(RotPeriodicSide2GlobalSide(nRotPeriodicSides)) +#if USE_MPI +! last proc knows CN total number of rot periodic sides +CALL MPI_BCAST(nRotPeriodicSides,1,MPI_INTEGER,nComputeNodeProcessors-1,MPI_COMM_SHARED,iError) +CALL BARRIER_AND_SYNC(SurfSide2RotPeriodicSide_Shared_Win , MPI_COMM_SHARED) +CALL BARRIER_AND_SYNC(Rot2Glob_temp_Shared_Win , MPI_COMM_SHARED) +CALL Allocate_Shared((/nRotPeriodicSides/) , NumRotPeriodicNeigh_Shared_Win , NumRotPeriodicNeigh_Shared) +CALL MPI_WIN_LOCK_ALL(0 , NumRotPeriodicNeigh_Shared_Win , IERROR) +NumRotPeriodicNeigh => NumRotPeriodicNeigh_Shared +CALL Allocate_Shared((/nRotPeriodicSides,NbrOfRotConnections/) , RotPeriodicSideMapping_temp_Shared_Win , RotPeriodicSideMapping_temp_Shared) +CALL MPI_WIN_LOCK_ALL(0 , RotPeriodicSideMapping_temp_Shared_Win , IERROR) +RotPeriodicSideMapping_temp => RotPeriodicSideMapping_temp_Shared +IF (myComputeNodeRank.EQ.0) NumRotPeriodicNeigh = 0 +IF (myComputeNodeRank.EQ.0) RotPeriodicSideMapping_temp = -1 +CALL BARRIER_AND_SYNC(NumRotPeriodicNeigh_Shared_Win , MPI_COMM_SHARED) +CALL BARRIER_AND_SYNC(RotPeriodicSideMapping_temp_Shared_Win , MPI_COMM_SHARED) +firstSide = INT(REAL( myComputeNodeRank *nRotPeriodicSides)/REAL(nComputeNodeProcessors))+1 +lastSide = INT(REAL((myComputeNodeRank+1)*nRotPeriodicSides)/REAL(nComputeNodeProcessors)) +#else +firstSide = 1 +lastSide = nRotPeriodicSides ALLOCATE(NumRotPeriodicNeigh(nRotPeriodicSides)) +NumRotPeriodicNeigh = 0 ! number of potential rotational periodic sides is unknown => allocate mapping array with fixed number of NbrOfRotConnections ! and reallocate at the end of subroutine ALLOCATE(RotPeriodicSideMapping_temp(nRotPeriodicSides,NbrOfRotConnections)) +RotPeriodicSideMapping_temp = -1 +#endif /*USE_MPI*/ -DO iSide=1, nRotPeriodicSides - RotPeriodicSide2GlobalSide(iSide) = Rot2Glob_temp(iSide) - NumRotPeriodicNeigh(iSide) = 0 - RotPeriodicSideMapping_temp(iSide,1:NbrOfRotConnections) = -1 -END DO notMapped=0 MaxNumRotPeriodicNeigh = 0 ! Defining rotation matrix @@ -615,87 +750,127 @@ SUBROUTINE InitParticleBoundaryRotPeriodic() k = 1 l = 2 m = 3 - CASE(2) ! x-rotation axis + CASE(2) ! y-rotation axis k = 2 l = 3 m = 1 - CASE(3) ! x-rotation axis + CASE(3) ! z-rotation axis k = 3 l = 1 m = 2 END SELECT -! (2) find Side on corresponding BC and build mapping RotPeriodicSide -> SideID2 (and vice versa) -! counting potential rotational periodic sides (for non-conforming meshes) -DO iSide=1, nRotPeriodicSides - SideID = RotPeriodicSide2GlobalSide(iSide) + +! (2) Build bounding boxes (in 2D reference system) for all nRotPeriodicSides +#if USE_MPI +CALL Allocate_Shared((/4,nRotPeriodicSides/) , BoundingBox_Shared_Win , BoundingBox_Shared) +CALL MPI_WIN_LOCK_ALL(0 , BoundingBox_Shared_Win , IERROR) +BoundingBox => BoundingBox_Shared +CALL BARRIER_AND_SYNC(BoundingBox_Shared_Win , MPI_COMM_SHARED) +#else +ALLOCATE(BoundingBox(4,nRotPeriodicSides)) +#endif /*USE_MPI*/ +DO iSide = firstSide, lastSide + + ! Get side information + SideID = Rot2Glob_temp(iSide) CNElemID = GetCNElemID(SideInfo_Shared(SIDE_ELEMID,SideID)) LocSideID = SideInfo_Shared(SIDE_LOCALID,SideID) - SideIsMapped = .FALSE. - ! check if at least one node of iSide is inside bounding box of a Side on corresponding BC + + ! Loop over all 4 nodes DO iNode=1, 4 - NodeID = ElemSideNodeID_Shared(iNode,LocSideID,CNElemID) + 1 - iNodeVec(1:3) = NodeCoords_Shared(1:3,NodeID) - iNodeH = iNodeVec(k) - iNodeR = SQRT(iNodeVec(l)*iNodeVec(l)+iNodeVec(m)*iNodeVec(m)) - DO jSide=1, nRotPeriodicSides - SideID2 = RotPeriodicSide2GlobalSide(jSide) - ! is on same RotPeriodicBC? - IF(PartBound%RotPeriodicDir(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).EQ. & - PartBound%RotPeriodicDir(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID2)))) CYCLE - isMapped = .FALSE. - ! check whether RotPeriodicSides is already mapped - IF(NumRotPeriodicNeigh(jSide).GT. 0) THEN - DO iNeigh=1, NumRotPeriodicNeigh(jSide) - IF(RotPeriodicSideMapping_temp(jSide,iNeigh).EQ.SideID) THEN - isMapped = .TRUE. - SideIsMapped = .TRUE. - EXIT - END IF - END DO - END IF - IF(isMapped) CYCLE - ! get ElemID for node mapping - CNElemID2 = GetCNElemID(SideInfo_Shared(SIDE_ELEMID,SideID2)) - LocSideID2 = SideInfo_Shared(SIDE_LOCALID,SideID2) - ! calc bounding box of jSide - DO jNode=1, 4 - NodeID2 = ElemSideNodeID_Shared(jNode,LocSideID2,CNElemID2) + 1 - jNodeVec(1:3) = NodeCoords_Shared(1:3,NodeID2) - jNodeH = jNodeVec(k) - jNodeR = SQRT(jNodeVec(l)*jNodeVec(l)+jNodeVec(m)*jNodeVec(m)) - IF(jNode.EQ. 1) THEN - Node2Hmin = jNodeH - Node2Hmax = jNodeH - Node2Rmin = jNodeR - Node2Rmax = jNodeR - ELSE - Node2Hmin = MIN(Node2Hmin,jNodeH) - Node2Hmax = MAX(Node2Hmax,jNodeH) - Node2Rmin = MIN(Node2Rmin,jNodeR) - Node2Rmax = MAX(Node2Rmax,jNodeR) - END IF - END DO - IF( ( (Node2Hmin.LE.iNodeH).AND.(iNodeH.LE.Node2Hmax) ) .AND. & - ( (Node2Rmin.LE.iNodeR).AND.(iNodeR.LE.Node2Rmax) ) ) THEN - ! at least one node of iSide is inside bounding box of jSide => - ! 1. increase NumRotPeriodicNeigh - ! 2. map: iSide => SideID2 and - ! vise versa: jSide => SideID s.o. - NumRotPeriodicNeigh(iSide) = NumRotPeriodicNeigh(iSide) + 1 - IF(NumRotPeriodicNeigh(iSide).GT. NbrOfRotConnections) THEN - CALL abort(__STAMP__,' ERROR: Number of rotational periodic side exceed fixed number of ',IntInfoOpt=NbrOfRotConnections) - END IF - RotPeriodicSideMapping_temp(iSide,NumRotPeriodicNeigh(iSide)) = SideID2 + NodeID = ElemSideNodeID_Shared(iNode,LocSideID,CNElemID) + 1 + jNodeVec(1:3) = NodeCoords_Shared(1:3,NodeID) + jNodeH = jNodeVec(k) + jNodeR = SQRT(jNodeVec(l)**2+jNodeVec(m)**2) + IF(iNode.EQ. 1) THEN + Node2Hmin = jNodeH + Node2Hmax = jNodeH + Node2Rmin = jNodeR + Node2Rmax = jNodeR + ELSE + Node2Hmin = MIN(Node2Hmin,jNodeH) + Node2Hmax = MAX(Node2Hmax,jNodeH) + Node2Rmin = MIN(Node2Rmin,jNodeR) + Node2Rmax = MAX(Node2Rmax,jNodeR) + END IF + END DO + ! Add tolerance by increasing the bounding box size by 1 percent of the length in h and r + dh = Node2Hmax-Node2Hmin + dr = Node2Rmax-Node2Rmin + BoundingBox(1,iSide) = Node2Hmin-dh*0.01 + BoundingBox(2,iSide) = Node2Hmax+dh*0.01 + BoundingBox(3,iSide) = Node2Rmin-dr*0.01 + BoundingBox(4,iSide) = Node2Rmax+dr*0.01 + +END DO ! iSide = firstSide, lastSide +#if USE_MPI +CALL BARRIER_AND_SYNC(BoundingBox_Shared_Win , MPI_COMM_SHARED) +#endif /*USE_MPI*/ + +! (3) find Side on corresponding BC and build mapping RotPeriodicSide -> SideID2 (and vice versa) +! counting potential rotational periodic sides (for non-conforming meshes) +! Use named loops: Loop over the assigned iSides and compare against all nRotPeriodicSides +iSideLoop: DO iSide = firstSide, lastSide + SideID = Rot2Glob_temp(iSide) + CNElemID = GetCNElemID(SideInfo_Shared(SIDE_ELEMID,SideID)) + LocSideID = SideInfo_Shared(SIDE_LOCALID,SideID) + + jSideLoop: DO jSide = 1, nRotPeriodicSides + SideID2 = Rot2Glob_temp(jSide) + + ! Check if both sides are on the same boundary, i.e., they cannot be connected + IF(PartBound%RotPeriodicDir(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).EQ. & + PartBound%RotPeriodicDir(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID2)))) CYCLE jSideLoop + + ! Check if jSide is assigned to the proc + mySide = (jSide.GE.firstSide).AND.(jSide.LE.lastSide) + + ! Check if the side was added in a previous step + IF(mySide)THEN + DO iNeigh = 1, NumRotPeriodicNeigh(jSide) + IF(RotPeriodicSideMapping_temp(jSide,iNeigh).EQ.SideID) CYCLE jSideLoop + END DO ! iNeigh = 1, NumRotPeriodicNeigh(jSide) + END IF ! mySide + + ! Loop the 4 nodes of iSide and test against the bounding box of jSide + iNodeLoop: DO iNode = 1, 4 + NodeID = ElemSideNodeID_Shared(iNode,LocSideID,CNElemID) + 1 + ! Calculate node coordinates in reference system + iNodeVec(1:3) = NodeCoords_Shared(1:3,NodeID) + iNodeH = iNodeVec(k) + iNodeR = SQRT(iNodeVec(l)**2+iNodeVec(m)**2) + + ! Cycle if outside of bounding box of jSide + IF(BoundingBox(1,jSide).GT.iNodeH) CYCLE iNodeLoop + IF(BoundingBox(2,jSide).LT.iNodeH) CYCLE iNodeLoop + IF(BoundingBox(3,jSide).GT.iNodeR) CYCLE iNodeLoop + IF(BoundingBox(4,jSide).LT.iNodeR) CYCLE iNodeLoop + + ! Check if the side was added in a previous step + DO iNeigh = 1, NumRotPeriodicNeigh(iSide) + IF(RotPeriodicSideMapping_temp(iSide,iNeigh).EQ.SideID2) CYCLE iNodeLoop + END DO ! iNeigh = 1, NumRotPeriodicNeigh(iSide) + + ! Found connection + NumRotPeriodicNeigh(iSide) = NumRotPeriodicNeigh(iSide) + 1 + IF(NumRotPeriodicNeigh(iSide).GT.NbrOfRotConnections) CALL abort(__STAMP__,& + ' ERROR: Number of rotational periodic side exceed fixed number of ',IntInfoOpt=NbrOfRotConnections) + RotPeriodicSideMapping_temp(iSide,NumRotPeriodicNeigh(iSide)) = SideID2 + + ! Only do vice versa mapping if the processor has been assigned this side + IF(mySide)THEN NumRotPeriodicNeigh(jSide) = NumRotPeriodicNeigh(jSide) + 1 - IF(NumRotPeriodicNeigh(jSide).GT. NbrOfRotConnections) THEN - CALL abort(__STAMP__,' ERROR: Number of rotational periodic side exceed fixed number of ',IntInfoOpt=NbrOfRotConnections) - END IF + IF(NumRotPeriodicNeigh(jSide).GT.NbrOfRotConnections) CALL abort(__STAMP__,& + ' jSide: Number of rotational periodic side exceed fixed number of ',IntInfoOpt=NbrOfRotConnections) RotPeriodicSideMapping_temp(jSide,NumRotPeriodicNeigh(jSide)) = SideID - SideIsMapped = .TRUE. - END IF - END DO ! jSide=1, nRotPeriodicSides - END DO ! iNode=1, 4 - IF(.NOT.SideIsMapped) THEN + END IF ! mySide + + END DO iNodeLoop ! iNode = 1, 4 + + END DO jSideLoop ! jSide = 1, nRotPeriodicSides + + ! Check if iSide could not be mapped to any other side. There should be at least one + IF(NumRotPeriodicNeigh(iSide).EQ.0) THEN IF(ElemInfo_Shared(ELEM_HALOFLAG,SideInfo_Shared(SIDE_ELEMID,SideID)).NE.3) THEN ! Count number of sides that could not be mapped (warning output + info in h5 file when CalcMeshInfo=T) notMapped = notMapped + 1 @@ -704,25 +879,47 @@ SUBROUTINE InitParticleBoundaryRotPeriodic() NumRotPeriodicNeigh(iSide) = 1 RotPeriodicSideMapping_temp(iSide,NumRotPeriodicNeigh(iSide)) = -1 END IF ! ElemInfo_Shared(ELEM_HALOFLAG,SideInfo_Shared(SIDE_ELEMID,SideID)).NE.3 - END IF ! .NOT.SideIsMapped -END DO ! iSide=1, nRotPeriodicSides + END IF ! NumRotPeriodicNeigh(iSide).EQ.0 -! (3) reallocate array due to number of potential rotational periodic sides +END DO iSideLoop ! iSide = firstSide, lastSide + +! (4) reallocate array due to number of potential rotational periodic sides +#if USE_MPI +CALL BARRIER_AND_SYNC(NumRotPeriodicNeigh_Shared_Win, MPI_COMM_SHARED) +CALL BARRIER_AND_SYNC(RotPeriodicSideMapping_temp_Shared_Win, MPI_COMM_SHARED) +! The allreduce is only required when a global array for writing to .h5 is to be used +!CALL MPI_ALLREDUCE(MAXVAL(NumRotPeriodicNeigh) , MaxNumRotPeriodicNeigh , 1 , MPI_INTEGER , MPI_MAX , MPI_COMM_WORLD , iError) +#endif /*USE_MPI*/ MaxNumRotPeriodicNeigh = MAXVAL(NumRotPeriodicNeigh) + +#if USE_MPI +CALL Allocate_Shared((/nRotPeriodicSides,NbrOfRotConnections/) , RotPeriodicSideMapping_Shared_Win , RotPeriodicSideMapping_Shared) +CALL MPI_WIN_LOCK_ALL(0 , RotPeriodicSideMapping_Shared_Win , IERROR) +RotPeriodicSideMapping => RotPeriodicSideMapping_Shared +IF (myComputeNodeRank.EQ.0) RotPeriodicSideMapping = -1 +CALL BARRIER_AND_SYNC(RotPeriodicSideMapping_Shared_Win , MPI_COMM_SHARED) +#else ALLOCATE(RotPeriodicSideMapping(nRotPeriodicSides,MaxNumRotPeriodicNeigh)) +RotPeriodicSideMapping = -1 +#endif /*USE_MPI*/ + DO iSide=1, nRotPeriodicSides DO iNeigh=1, MaxNumRotPeriodicNeigh - RotPeriodicSideMapping(iSide,iNeigh) = RotPeriodicSideMapping_temp(iSide,iNeigh) + SideID = RotPeriodicSideMapping_temp(iSide,iNeigh) + IF(SideID.NE.-1)THEN + GlobalElemID = SideInfo_Shared(SIDE_ELEMID,SideID) + RotPeriodicSideMapping(iSide,iNeigh) = GlobalElemID + END IF ! SideID.NE.-1 END DO END DO #if USE_MPI CALL MPI_ALLREDUCE(notMapped , notMappedTotal , 1 , MPI_INTEGER , MPI_SUM , MPI_COMM_WORLD , IERROR) +#else +notMappedTotal = notMapped +#endif /*USE_MPI*/ + IF(notMappedTotal.GT.0)THEN - SWRITE(Unit_StdOut,'(A,I0,A)') ' | Warning: Found ',notMappedTotal,' rot periodic sides that did not find a corresponding side.' - SWRITE(Unit_StdOut,'(A)')" | The halo region (halo flag 2) merely reaches a rot periodic BC side but not any further." - SWRITE(Unit_StdOut,'(A)')" | See ElemData container 'LostRotPeriodicSides' for more information on where sides were unmatched." - SWRITE(Unit_StdOut,'(A)')" | This information is written to "//TRIM(ProjectName)//"_LostRotPeriodicSides.h5 (only when CalcMeshInfo=T)" IF(CalcMeshInfo)THEN ALLOCATE(LostRotPeriodicSides(1:nElems)) LostRotPeriodicSides=notMapped @@ -731,19 +928,25 @@ SUBROUTINE InitParticleBoundaryRotPeriodic() END IF ! CalcMeshInfo !IF(MPIroot) CALL abort(__STAMP__,'At least one rot periodic side did not find a corresponding side.') END IF ! notMappedTotal.GT.0 -#endif /*USE_MPI*/ #if USE_MPI -IF(MPIRoot)THEN - EndT=MPI_WTIME() - LBWRITE(UNIT_stdOut,'(A,F0.3,A)') ' INIT ROTATIONAL PERIODIC BOUNDARY DONE! [',EndT-StartT,'s]' - LBWRITE(UNIT_StdOut,'(132("-"))') -END IF ! MPIRoot +CALL MPI_BARRIER(MPI_COMM_SHARED,iERROR) +CALL UNLOCK_AND_FREE(Rot2Glob_temp_Shared_Win) +ADEALLOCATE(Rot2Glob_temp_Shared) +ADEALLOCATE(Rot2Glob_temp) +CALL UNLOCK_AND_FREE(RotPeriodicSideMapping_temp_Shared_Win) +ADEALLOCATE(RotPeriodicSideMapping_temp_Shared) +ADEALLOCATE(RotPeriodicSideMapping_temp) +CALL UNLOCK_AND_FREE(BoundingBox_Shared_Win) +ADEALLOCATE(BoundingBox_Shared) +ADEALLOCATE(BoundingBox) #else -WRITE(UNIT_stdOut,'(A)') ' INIT ROTATIONAL PERIODIC BOUNDARY DONE!' +DEALLOCATE(Rot2Glob_temp) +DEALLOCATE(RotPeriodicSideMapping_temp) +DEALLOCATE(BoundingBox) #endif /*USE_MPI*/ -END SUBROUTINE InitParticleBoundaryRotPeriodic +END SUBROUTINE BuildParticleBoundaryRotPeriodic !=================================================================================================================================== @@ -1092,12 +1295,13 @@ SUBROUTINE FinalizeParticleBoundary() USE MOD_Globals USE MOD_Particle_Boundary_Vars #if USE_MPI -USE MOD_MPI_Shared_vars ,ONLY: MPI_COMM_SHARED +USE MOD_MPI_Shared_vars ,ONLY: MPI_COMM_SHARED USE MOD_MPI_Shared #endif #if USE_LOADBALANCE -USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance,UseH5IOLoadBalance +USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance,UseH5IOLoadBalance #endif /*USE_LOADBALANCE*/ +USE MOD_Particle_Mesh_Vars ,ONLY: GEO !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT VARIABLES @@ -1136,11 +1340,27 @@ SUBROUTINE FinalizeParticleBoundary() SDEALLOCATE(PartBound%Dielectric) SDEALLOCATE(PartBound%BoundaryParticleOutputHDF5) SDEALLOCATE(PartBound%RadiativeEmissivity) + ! Rotational periodic boundary condition -SDEALLOCATE(RotPeriodicSide2GlobalSide) -SDEALLOCATE(NumRotPeriodicNeigh) -SDEALLOCATE(RotPeriodicSideMapping) -SDEALLOCATE(SurfSide2RotPeriodicSide) +IF(GEO%RotPeriodicBC)THEN +#if USE_MPI + CALL MPI_BARRIER(MPI_COMM_SHARED,iERROR) + CALL UNLOCK_AND_FREE(SurfSide2RotPeriodicSide_Shared_Win) + CALL UNLOCK_AND_FREE(NumRotPeriodicNeigh_Shared_Win) + ADEALLOCATE(SurfSide2RotPeriodicSide_Shared) + ADEALLOCATE(SurfSide2RotPeriodicSide) + ADEALLOCATE(NumRotPeriodicNeigh_Shared) + ADEALLOCATE(NumRotPeriodicNeigh) + CALL UNLOCK_AND_FREE(RotPeriodicSideMapping_Shared_Win) + ADEALLOCATE(RotPeriodicSideMapping_Shared) + ADEALLOCATE(RotPeriodicSideMapping) +#else + SDEALLOCATE(SurfSide2RotPeriodicSide) + SDEALLOCATE(NumRotPeriodicNeigh) + SDEALLOCATE(RotPeriodicSideMapping) +#endif +END IF ! GEO%RotPeriodicBC + ! Adaptive wall temperature (e.g. calculate from sampled heat flux) IF (ANY(PartBound%UseAdaptedWallTemp)) THEN #if USE_MPI diff --git a/src/particles/boundary/particle_boundary_vars.f90 b/src/particles/boundary/particle_boundary_vars.f90 index 13f1754cd..23fdd744a 100644 --- a/src/particles/boundary/particle_boundary_vars.f90 +++ b/src/particles/boundary/particle_boundary_vars.f90 @@ -107,10 +107,23 @@ MODULE MOD_Particle_Boundary_Vars ! ==================================================================== ! Rotational periodic sides -INTEGER,ALLOCATABLE :: RotPeriodicSide2GlobalSide(:) ! Mapping BC-side with PartBoundCond=6 to Global Side ID -INTEGER,ALLOCATABLE :: NumRotPeriodicNeigh(:) ! Number of adjacent Neigbours sites in rotational periodic BC -INTEGER,ALLOCATABLE :: RotPeriodicSideMapping(:,:) ! Mapping between rotational periodic sides. -INTEGER,ALLOCATABLE :: SurfSide2RotPeriodicSide(:) ! Mapping between surf side and periodic sides. +INTEGER,ALLOCPOINT,DIMENSION(:) :: NumRotPeriodicNeigh ! Number of adjacent Neigbours sites in rotational periodic BC +INTEGER,ALLOCPOINT,DIMENSION(:,:) :: RotPeriodicSideMapping ! Mapping between rotational periodic sides. +INTEGER,ALLOCPOINT,DIMENSION(:) :: SurfSide2RotPeriodicSide ! Mapping between surf side and periodic sides. +#if USE_MPI +INTEGER,POINTER,DIMENSION(:) :: SurfSide2RotPeriodicSide_Shared +INTEGER :: SurfSide2RotPeriodicSide_Shared_Win +INTEGER,POINTER,DIMENSION(:) :: NumRotPeriodicNeigh_Shared +INTEGER :: NumRotPeriodicNeigh_Shared_Win +INTEGER,POINTER,DIMENSION(:) :: Rot2Glob_temp_Shared +INTEGER :: Rot2Glob_temp_Shared_Win +INTEGER,POINTER,DIMENSION(:,:) :: RotPeriodicSideMapping_temp_Shared +INTEGER :: RotPeriodicSideMapping_temp_Shared_Win +INTEGER,POINTER,DIMENSION(:,:) :: RotPeriodicSideMapping_Shared +INTEGER :: RotPeriodicSideMapping_Shared_Win +REAL,POINTER,DIMENSION(:,:) :: BoundingBox_Shared +INTEGER :: BoundingBox_Shared_Win +#endif /*USE_MPI*/ !----------------------------------------------------------------------------------------------------------------------------------- ! required variables diff --git a/src/particles/particle_init.f90 b/src/particles/particle_init.f90 index dd4852fc1..d3d0f0a24 100644 --- a/src/particles/particle_init.f90 +++ b/src/particles/particle_init.f90 @@ -341,7 +341,7 @@ SUBROUTINE InitParticles() ! (the following IF arguments have to be considered in FinalizeParticleBoundarySampling as well) IF (WriteMacroSurfaceValues.OR.DSMC%CalcSurfaceVal.OR.(ANY(PartBound%Reactive)).OR.(nPorousBC.GT.0).OR.GEO%RotPeriodicBC) THEN CALL InitParticleBoundarySampling() - CALL InitParticleBoundaryRotPeriodic() + IF(GEO%RotPeriodicBC) CALL InitParticleBoundaryRotPeriodic() CALL InitAdaptiveWallTemp() END IF