diff --git a/docs/documentation/developerguide/code_extension.md b/docs/documentation/developerguide/code_extension.md index 9424271d9..b2de53fa1 100644 --- a/docs/documentation/developerguide/code_extension.md +++ b/docs/documentation/developerguide/code_extension.md @@ -34,4 +34,43 @@ Finally, the `WriteSurfSampleToHDF5` routine writes the prepared `MacroSurfaceVa IF (ANY(PartBound%SurfaceModel.EQ.1)) CALL AddVarName(Str2DVarNames,nVar2D_Total,nVarCount,'Sticking_Coefficient') -The order of the variable names and their position in the `MacroSurfaceVal` array has to be the same. Thus, make sure to place the `AddVarName` call at the same position, where you placed the calculation and writing into the `MacroSurfaceVal` array, otherwise the names and values will be mixed up. \ No newline at end of file +The order of the variable names and their position in the `MacroSurfaceVal` array has to be the same. Thus, make sure to place the `AddVarName` call at the same position, where you placed the calculation and writing into the `MacroSurfaceVal` array, otherwise the names and values will be mixed up. + +## Arrays with size of PDM%maxParticleNumber + +If an array is to store particle information, it is usually allocated with + + ALLOCATE(ParticleInformation(1:PDM%maxParticleNumber)) + +But since PDM%maxParticleNumber is dynamic, this behavior must also be captured by this array. Therefore its size has to be changed in the routines `IncreaseMaxParticleNumber` and `ReduceMaxParticleNumber` in `src/particles/particle_tools.f90`. + + IF(ALLOCATED(ParticleInformation)) CALL ChangeSizeArray(ParticleInformation,PDM%maxParticleNumber,NewSize, Default) + +Default is an optional parameter if the new array memory is to be initialized with a specific value. The same must be done for TYPES of size PDM%maxParticleNumber. Please check both routines to see how to do it. + +## Insert new particles + +To add new particles, first create a new particle ID using the GetNextFreePosition function contained in `src/particles/particle_tools.f90` + + NewParticleID = GetNextFreePosition() + +This directly increments the variable PDM%CurrentNextFreePosition by 1 and if necessary adjusts PDM%ParticleVecLength by 1. If this is not desired, it is possible to pass an offset. Then the two variables will not be incremented, which must be done later by the developer. This can happen if the particle generation process is divided into several functions, where each function contains a loop over all new particles (e.g. `src/particles/emission/particle_emission.f90`). + + DO iPart=1,nPart + NewParticleID = GetNextFreePosition(iPart) + END DO + PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + nPart + PDM%ParticleVecLength = MAX(PDM%ParticleVecLength,GetNextFreePosition(0)) + +For the new particle to become a valid particle, the inside flag must be set to true and various other arrays must be filled with meaningful data. See SUBROUTINE CreateParticle in `src/particles/particle_operations.f90`. A basic example of the most important variables is given below: + + newParticleID = GetNextFreePosition() + PDM%ParticleInside(newParticleID) = .TRUE. + PDM%FracPush(newParticleID) = .FALSE. + PDM%IsNewPart(newParticleID) = .TRUE. + PEM%GlobalElemID(newParticleID) = GlobElemID + PEM%LastGlobalElemID(newParticleID) = GlobElemID + PartSpecies(newParticleID) = SpecID + LastPartPos(1:3,newParticleID) = Pos(1:3) + PartState(1:3,newParticleID) = Pos(1:3) + PartState(4:6,newParticleID) = Velocity(1:3) \ No newline at end of file diff --git a/docs/documentation/userguide/features-and-models/particle-initialization-and-emission.md b/docs/documentation/userguide/features-and-models/particle-initialization-and-emission.md index 182fa03e3..de205571b 100644 --- a/docs/documentation/userguide/features-and-models/particle-initialization-and-emission.md +++ b/docs/documentation/userguide/features-and-models/particle-initialization-and-emission.md @@ -1,18 +1,21 @@ (sec:particle-initialization-and-emission)= # Particle Initialization & Emission +The RAM to store the particles is dynamically allocated. However, it is possible to restrict the number of particles per MPI process by setting + + Part-MaxParticleNumber=1000000 + +New memory is allocated in separate chunks because allocating memory for the particle data and copying it to the new memory area is expensive. The chunksize is relative to the particles used and can be set with + + Part-MaxPartNumIncrease=0.1 + +A higher value increases the amount of unnecessary RAM allocated to particles, while a lower value increases the number of memory adjustment operations. The optimal trade-off depends on the simulation and the machine, but it only affects the performance of the simulations, not the quality of the results. + The following section gives an overview of the available options regarding the definition of species and particle initialization and emission. Simulation particles can be inserted initially within the computational domain and/or emitted at every time step. First of all, the number of species is defined by Part-nSpecies=1 - Part-MaxParticleNumber=1000000 - -The maximum particle number is defined per core and should be chosen according to the number of simulation particles you expect, -including a margin to account for imbalances due transient flow features and/or the occurrence of new particles due to chemical -reactions. Example: A node of a HPC cluster has 2 CPUs, each has 12 cores. Thus, the node has 24 cores that share a total of -128GB RAM. Allocating 1000000 particles per core means, you can simulate up to 24 Million particles on a single node in this -example (assuming an even particle distribution). The limiting factor here is the amount of RAM available per node. Regardless whether a standalone PIC, DSMC, or a coupled simulation is performed, the atomic mass [kg], the charge [C] and the weighting factor $w$ [-], sometimes referred to as macro-particle factor (MPF), are required for each species. diff --git a/regressioncheck/NIG_tracking_DSMC/curved/parameter.ini b/regressioncheck/NIG_tracking_DSMC/curved/parameter.ini index 4d9d4567c..ec94670d7 100644 --- a/regressioncheck/NIG_tracking_DSMC/curved/parameter.ini +++ b/regressioncheck/NIG_tracking_DSMC/curved/parameter.ini @@ -21,6 +21,7 @@ RefMappingEps = 1e-5 RefMappingGuess = 3 BezierNewtonTolerance = 1e-4 BezierSplitLimit = 0.8 +Part-MPI-UNFP-afterPartSend = true ! =============================================================================== ! ! OUTPUT / VISUALIZATION ! =============================================================================== ! diff --git a/src/globals/array_operations.f90 b/src/globals/array_operations.f90 new file mode 100644 index 000000000..9bb99dc50 --- /dev/null +++ b/src/globals/array_operations.f90 @@ -0,0 +1,328 @@ +!================================================================================================================================== +! Copyright (c) 2010 - 2018 Prof. Claus-Dieter Munz and Prof. Stefanos Fasoulas +! +! This file is part of PICLas (piclas.boltzplatz.eu/piclas/piclas). PICLas is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 +! of the License, or (at your option) any later version. +! +! PICLas is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License v3.0 for more details. +! +! You should have received a copy of the GNU General Public License along with PICLas. If not, see . +!================================================================================================================================== +#include "piclas.h" + +MODULE MOD_Array_Operations +!=================================================================================================================================== +! Contains tools for array related operations. +!=================================================================================================================================== +! MODULES +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +PRIVATE + +INTERFACE ChangeSizeArray + MODULE PROCEDURE ChangeSizeArrayLOG1, ChangeSizeArrayLOG2, & + ChangeSizeArrayINT1, ChangeSizeArrayINT2, & + ChangeSizeArrayREAL1, ChangeSizeArrayREAL2, ChangeSizeArrayREAL3 +END INTERFACE + +!----------------------------------------------------------------------------------------------------------------------------------- +! GLOBAL VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! Private Part --------------------------------------------------------------------------------------------------------------------- +! Public Part ---------------------------------------------------------------------------------------------------------------------- +PUBLIC :: ChangeSizeArray +!=================================================================================================================================== + +CONTAINS + + +SUBROUTINE ChangeSizeArrayLOG1(Vec,OldSize,NewSize,Default) +!=================================================================================================================================== +! Change size of a vector, Logical 1D +! If NewSize>OldSize and PRESENT(Default): New entries are initialized with default +!=================================================================================================================================== +! MODULES +USE MOD_Globals +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +LOGICAL,ALLOCATABLE,INTENT(INOUT) :: Vec(:) +INTEGER,INTENT(IN) :: OldSize, NewSize +LOGICAL,INTENT(IN),OPTIONAL :: Default +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +LOGICAL,ALLOCATABLE :: TempVec(:) +INTEGER :: ALLOCSTAT +!=================================================================================================================================== +! Allocate new memory space +ALLOCATE(TempVec(NewSize),STAT=ALLOCSTAT) +IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate new Array in ChangeSizeArray') + +! Write old data to new memory space +IF(NewSize.GT.OldSize) THEN + TempVec(1:OldSize)=Vec + IF(PRESENT(Default)) TempVec(OldSize+1:NewSize) = Default +ELSE + TempVec=Vec(1:NewSize) +END IF + +! Switch array pointer +CALL MOVE_ALLOC(TempVec,Vec) + +END SUBROUTINE ChangeSizeArrayLOG1 + + +SUBROUTINE ChangeSizeArrayLOG2(Vec,OldSize,NewSize,Default) +!=================================================================================================================================== +! Change size of a vector, Logical 2D, last dimension +! If NewSize>OldSize and PRESENT(Default): New entries are initialized with default +!=================================================================================================================================== +! MODULES +USE MOD_Globals +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +LOGICAL,ALLOCATABLE,INTENT(INOUT) :: Vec(:,:) +INTEGER,INTENT(IN) :: OldSize, NewSize +LOGICAL,INTENT(IN),OPTIONAL :: Default +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +LOGICAL,ALLOCATABLE :: TempVec(:,:) +INTEGER :: ALLOCSTAT +!=================================================================================================================================== +! Allocate new memory space +ALLOCATE(TempVec(SIZE(Vec,1),NewSize),STAT=ALLOCSTAT) +IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate new Array in ChangeSizeArray') + +! Write old data to new memory space +IF(NewSize.GT.OldSize) THEN + TempVec(:,1:OldSize)=Vec + IF(PRESENT(Default)) TempVec(:,OldSize+1:NewSize) = Default +ELSE + TempVec=Vec(:,1:NewSize) +END IF + +! Switch array pointer +CALL MOVE_ALLOC(TempVec,Vec) + +END SUBROUTINE ChangeSizeArrayLOG2 + + +SUBROUTINE ChangeSizeArrayINT1(Vec,OldSize,NewSize,Default) +!=================================================================================================================================== +! Change size of a vector, Integer 1D +! If NewSize>OldSize and PRESENT(Default): New entries are initialized with default +!=================================================================================================================================== +! MODULES +USE MOD_Globals +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,ALLOCATABLE,INTENT(INOUT) :: Vec(:) +INTEGER,INTENT(IN) :: OldSize, NewSize +INTEGER,INTENT(IN),OPTIONAL :: Default +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER,ALLOCATABLE :: TempVec(:) +INTEGER :: ALLOCSTAT +!=================================================================================================================================== +! Allocate new memory space +ALLOCATE(TempVec(NewSize),STAT=ALLOCSTAT) +IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate new Array in ChangeSizeArray') + +! Write old data to new memory space +IF(NewSize.GT.OldSize) THEN + TempVec(1:OldSize)=Vec + IF(PRESENT(Default)) TempVec(OldSize+1:NewSize) = Default +ELSE + TempVec=Vec(1:NewSize) +END IF + +! Switch array pointer +CALL MOVE_ALLOC(TempVec,Vec) + +END SUBROUTINE ChangeSizeArrayINT1 + + +SUBROUTINE ChangeSizeArrayINT2(Vec,OldSize,NewSize,Default) +!=================================================================================================================================== +! Change size of a vector, Integer 2D, last dimension +! If NewSize>OldSize and PRESENT(Default): New entries are initialized with default +!=================================================================================================================================== +! MODULES +USE MOD_Globals +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,ALLOCATABLE,INTENT(INOUT) :: Vec(:,:) +INTEGER,INTENT(IN) :: OldSize, NewSize +INTEGER,INTENT(IN),OPTIONAL :: Default +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER,ALLOCATABLE :: TempVec(:,:) +INTEGER :: ALLOCSTAT +!=================================================================================================================================== +! Allocate new memory space +ALLOCATE(TempVec(SIZE(Vec,1),NewSize),STAT=ALLOCSTAT) +IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate new Array in ChangeSizeArray') + +! Write old data to new memory space +IF(NewSize.GT.OldSize) THEN + TempVec(:,1:OldSize)=Vec + IF(PRESENT(Default)) TempVec(:,OldSize+1:NewSize) = Default +ELSE + TempVec=Vec(:,1:NewSize) +END IF + +! Switch array pointer +CALL MOVE_ALLOC(TempVec,Vec) + +END SUBROUTINE ChangeSizeArrayINT2 + + +SUBROUTINE ChangeSizeArrayREAL1(Vec,OldSize,NewSize,Default) +!=================================================================================================================================== +! Change size of a vector, Real 1D +! If NewSize>OldSize and PRESENT(Default): New entries are initialized with default +!=================================================================================================================================== +! MODULES +USE MOD_Globals +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +REAL,ALLOCATABLE,INTENT(INOUT) :: Vec(:) +INTEGER,INTENT(IN) :: OldSize, NewSize +REAL,INTENT(IN),OPTIONAL :: Default +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +REAL,ALLOCATABLE :: TempVec(:) +INTEGER :: ALLOCSTAT +!=================================================================================================================================== +! Allocate new memory space +ALLOCATE(TempVec(NewSize),STAT=ALLOCSTAT) +IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate new Array in ChangeSizeArray') + +! Write old data to new memory space +IF(NewSize.GT.OldSize) THEN + TempVec(1:OldSize)=Vec + IF(PRESENT(Default)) TempVec(OldSize+1:NewSize) = Default +ELSE + TempVec=Vec(1:NewSize) +END IF + +! Switch array pointer +CALL MOVE_ALLOC(TempVec,Vec) + +END SUBROUTINE ChangeSizeArrayREAL1 + + +SUBROUTINE ChangeSizeArrayREAL2(Vec,OldSize,NewSize,Default) +!=================================================================================================================================== +! Change size of a vector, Real 2D, last dimension +! If NewSize>OldSize and PRESENT(Default): New entries are initialized with default +!=================================================================================================================================== +! MODULES +USE MOD_Globals +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +REAL,ALLOCATABLE,INTENT(INOUT) :: Vec(:,:) +INTEGER,INTENT(IN) :: OldSize, NewSize +REAL,INTENT(IN),OPTIONAL :: Default +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +REAL,ALLOCATABLE :: TempVec(:,:) +INTEGER :: ALLOCSTAT +!=================================================================================================================================== +! Allocate new memory space +ALLOCATE(TempVec(SIZE(Vec,1),NewSize),STAT=ALLOCSTAT) +IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate new Array in ChangeSizeArray') + +! Write old data to new memory space +IF(NewSize.GT.OldSize) THEN + TempVec(:,1:OldSize)=Vec + IF(PRESENT(Default)) TempVec(:,OldSize+1:NewSize) = Default +ELSE + TempVec=Vec(:,1:NewSize) +END IF + +! Switch array pointer +CALL MOVE_ALLOC(TempVec,Vec) + +END SUBROUTINE ChangeSizeArrayREAL2 + + +SUBROUTINE ChangeSizeArrayREAL3(Vec,OldSize,NewSize,Default) +!=================================================================================================================================== +! Change size of a vector, Real 3D, last dimension +! If NewSize>OldSize and PRESENT(Default): New entries are initialized with default +!=================================================================================================================================== +! MODULES +USE MOD_Globals +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +REAL,ALLOCATABLE,INTENT(INOUT) :: Vec(:,:,:) +INTEGER,INTENT(IN) :: OldSize, NewSize +REAL,INTENT(IN),OPTIONAL :: Default +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +REAL,ALLOCATABLE :: TempVec(:,:,:) +INTEGER :: ALLOCSTAT +!=================================================================================================================================== +! Allocate new memory space +ALLOCATE(TempVec(SIZE(Vec,1),SIZE(Vec,2),NewSize),STAT=ALLOCSTAT) +IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate new Array in ChangeSizeArray') + +! Write old data to new memory space +IF(NewSize.GT.OldSize) THEN + TempVec(:,:,1:OldSize)=Vec + IF(PRESENT(Default)) TempVec(:,:,OldSize+1:NewSize) = Default +ELSE + TempVec=Vec(:,:,1:NewSize) +END IF + +! Switch array pointer +CALL MOVE_ALLOC(TempVec,Vec) + +END SUBROUTINE ChangeSizeArrayREAL3 + + +END MODULE MOD_Array_Operations diff --git a/src/particles/dsmc/dsmc_ambipolardiff.f90 b/src/particles/dsmc/dsmc_ambipolardiff.f90 index bc0bd4754..3262a9c1a 100644 --- a/src/particles/dsmc/dsmc_ambipolardiff.f90 +++ b/src/particles/dsmc/dsmc_ambipolardiff.f90 @@ -80,7 +80,7 @@ SUBROUTINE AD_SetInitElectronVelo(FractNbr,iInit,NbrOfParticle) USE MOD_Particle_Vars USE MOD_Globals_Vars ,ONLY: BoltzmannConst USE MOD_part_emission_tools ,ONLY: CalcVelocity_maxwell_lpn -USE MOD_part_tools ,ONLY: BuildTransGaussNums +USE MOD_part_tools ,ONLY: BuildTransGaussNums, GetNextFreePosition USE MOD_DSMC_Vars ,ONLY: DSMC, AmbipolElecVelo ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -117,33 +117,27 @@ SUBROUTINE AD_SetInitElectronVelo(FractNbr,iInit,NbrOfParticle) SELECT CASE(TRIM(velocityDistribution)) CASE('constant') DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr.GT.0) THEN - IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) - ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) - AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = VeloVecIC(1:3) * VeloIC - END IF + PositionNbr = GetNextFreePosition(i) + IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) + ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) + AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = VeloVecIC(1:3) * VeloIC END DO CASE('maxwell_lpn') DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr.GT.0) THEN - CALL CalcVelocity_maxwell_lpn(DSMC%AmbiDiffElecSpec, Vec3D, Temperature=Species(FractNbr)%Init(iInit)%MWTemperatureIC) - IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) - ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) - AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = VeloIC *VeloVecIC(1:3) + Vec3D(1:3) - END IF + PositionNbr = GetNextFreePosition(i) + CALL CalcVelocity_maxwell_lpn(DSMC%AmbiDiffElecSpec, Vec3D, Temperature=Species(FractNbr)%Init(iInit)%MWTemperatureIC) + IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) + ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) + AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = VeloIC *VeloVecIC(1:3) + Vec3D(1:3) END DO CASE('maxwell') CALL BuildTransGaussNums(NbrOfParticle, iRanPart) maxwellfac = SQRT(BoltzmannConst*Species(FractNbr)%Init(iInit)%MWTemperatureIC/Species(DSMC%AmbiDiffElecSpec)%MassIC) DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr.GT.0) THEN - IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) - ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) - AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = VeloIC *VeloVecIC(1:3) + iRanPart(1:3,i)*maxwellfac - END IF + PositionNbr = GetNextFreePosition(i) + IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) + ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) + AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = VeloIC *VeloVecIC(1:3) + iRanPart(1:3,i)*maxwellfac END DO CASE DEFAULT CALL abort(& @@ -166,6 +160,7 @@ SUBROUTINE AD_SetSFElectronVelo(iSpec,iSFIon,iSample,jSample,iSide,BCSideID,Side USE MOD_Particle_Surfaces, ONLY : CalcNormAndTangBezier USE MOD_Particle_Sampling_Vars ,ONLY: AdaptBCMapElemToSample, AdaptBCMacroVal USE MOD_DSMC_Vars ,ONLY: AmbiPolarSFMapping, AmbipolElecVelo, DSMC +USE MOD_Part_Tools ,ONLY: GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -256,31 +251,29 @@ SUBROUTINE AD_SetSFElectronVelo(iSpec,iSFIon,iSample,jSample,iSide,BCSideID,Side iPart = 0 DO i = NbrOfParticle-PartIns+1,NbrOfParticle iPart = iPart + 1 - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr .NE. 0) THEN - ! In case of side-normal velocities: calc n-vector at particle position, xi was saved in PartState(4:5) - IF (Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. TriaSurfaceFlux) THEN - vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) - vec_t1(1:3) = 0. !dummy - vec_t2(1:3) = 0. !dummy - ELSE IF (Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN - ! CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),xi=PartState(4,PositionNbr),eta=PartState(5,PositionNbr),SideID=SideID ) - CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),xi=particle_xis(2*(iPart-1)+1),eta=particle_xis(2*(iPart-1)+2),SideID=SideID ) - vec_nIn(1:3) = -vec_nIn(1:3) - vec_t1(1:3) = 0. !dummy - vec_t2(1:3) = 0. !dummy - ELSE - vec_nIn(1:3) = VeloVecIC(1:3) - END IF !VeloIsNormal - ! Build complete velo-vector - Vec3D(1:3) = vec_nIn(1:3) * Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIC - ! PartState(4:6,PositionNbr) = Vec3D(1:3) - IF (PositionNbr.GT.0) THEN - IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) - ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) - AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = Vec3D(1:3) - END IF - END IF !PositionNbr .NE. 0 + PositionNbr = GetNextFreePosition(i) + ! In case of side-normal velocities: calc n-vector at particle position, xi was saved in PartState(4:5) + IF (Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. TriaSurfaceFlux) THEN + vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) + vec_t1(1:3) = 0. !dummy + vec_t2(1:3) = 0. !dummy + ELSE IF (Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN + ! CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),xi=PartState(4,PositionNbr),eta=PartState(5,PositionNbr),SideID=SideID ) + CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),xi=particle_xis(2*(iPart-1)+1),eta=particle_xis(2*(iPart-1)+2),SideID=SideID ) + vec_nIn(1:3) = -vec_nIn(1:3) + vec_t1(1:3) = 0. !dummy + vec_t2(1:3) = 0. !dummy + ELSE + vec_nIn(1:3) = VeloVecIC(1:3) + END IF !VeloIsNormal + ! Build complete velo-vector + Vec3D(1:3) = vec_nIn(1:3) * Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIC + ! PartState(4:6,PositionNbr) = Vec3D(1:3) + IF (PositionNbr.GT.0) THEN + IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) + ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) + AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = Vec3D(1:3) + END IF END DO !i = ...NbrOfParticle CASE('maxwell','maxwell_lpn') !-- determine envelope for most efficient ARM [Garcia and Wagner 2006, JCP217-2] @@ -306,151 +299,145 @@ SUBROUTINE AD_SetSFElectronVelo(iSpec,iSFIon,iSample,jSample,iSide,BCSideID,Side iPart = 0 DO i = NbrOfParticle-PartIns+1,NbrOfParticle iPart = iPart + 1 - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr .NE. 0) THEN - !-- 0a.: In case of side-normal velocities: calc n-/t-vectors at particle position, xi was saved in PartState(4:5) - IF (Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. TriaSurfaceFlux) THEN - vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) - vec_t1(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_t1(1:3) - vec_t2(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_t2(1:3) - ELSE IF (Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN - ! CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),tang1=vec_t1(1:3),tang2=vec_t2(1:3) & - ! ,xi=PartState(4,PositionNbr),eta=PartState(5,PositionNbr),SideID=SideID ) - CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),tang1=vec_t1(1:3),tang2=vec_t2(1:3) & - ,xi=particle_xis(2*(iPart-1)+1),eta=particle_xis(2*(iPart-1)+2),SideID=SideID ) - vec_nIn(1:3) = -vec_nIn(1:3) - END IF !VeloIsNormal - !-- 1.: determine zstar (initial generation of potentially too many RVu is for needed indentities of RVu used multiple times! - SELECT CASE(envelope) - CASE(0) - CALL RANDOM_NUMBER(RandVal1) - zstar = -SQRT(-LOG(RandVal1)) - CASE(1) - DO - CALL RANDOM_NUMBER(RandVal2) - zstar = -SQRT(a*a-LOG(RandVal2(1))) - IF ( -(a-zstar)/zstar .GT. RandVal2(2)) THEN + PositionNbr = GetNextFreePosition(i) + !-- 0a.: In case of side-normal velocities: calc n-/t-vectors at particle position, xi was saved in PartState(4:5) + IF (Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. TriaSurfaceFlux) THEN + vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) + vec_t1(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_t1(1:3) + vec_t2(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_t2(1:3) + ELSE IF (Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN + ! CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),tang1=vec_t1(1:3),tang2=vec_t2(1:3) & + ! ,xi=PartState(4,PositionNbr),eta=PartState(5,PositionNbr),SideID=SideID ) + CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),tang1=vec_t1(1:3),tang2=vec_t2(1:3) & + ,xi=particle_xis(2*(iPart-1)+1),eta=particle_xis(2*(iPart-1)+2),SideID=SideID ) + vec_nIn(1:3) = -vec_nIn(1:3) + END IF !VeloIsNormal + !-- 1.: determine zstar (initial generation of potentially too many RVu is for needed indentities of RVu used multiple times! + SELECT CASE(envelope) + CASE(0) + CALL RANDOM_NUMBER(RandVal1) + zstar = -SQRT(-LOG(RandVal1)) + CASE(1) + DO + CALL RANDOM_NUMBER(RandVal2) + zstar = -SQRT(a*a-LOG(RandVal2(1))) + IF ( -(a-zstar)/zstar .GT. RandVal2(2)) THEN + EXIT + END IF + END DO + CASE(2) + z = 0.5*(a-SQRT(a*a+2.)) + beta = a-(1.0-a)*(a-z) + DO + CALL RANDOM_NUMBER(RandVal3) + IF (EXP(-(beta*beta))/(EXP(-(beta*beta))+2.0*(a-z)*(a-beta)*EXP(-(z*z))).GT.RandVal3(1)) THEN + zstar=-SQRT(beta*beta-LOG(RandVal3(2))) + IF ( -(a-zstar)/zstar .GT. RandVal3(3)) THEN EXIT END IF - END DO - CASE(2) - z = 0.5*(a-SQRT(a*a+2.)) - beta = a-(1.0-a)*(a-z) - DO - CALL RANDOM_NUMBER(RandVal3) - IF (EXP(-(beta*beta))/(EXP(-(beta*beta))+2.0*(a-z)*(a-beta)*EXP(-(z*z))).GT.RandVal3(1)) THEN - zstar=-SQRT(beta*beta-LOG(RandVal3(2))) - IF ( -(a-zstar)/zstar .GT. RandVal3(3)) THEN - EXIT - END IF - ELSE - zstar=beta+(a-beta)*RandVal3(2) - IF ( (a-zstar)/(a-z)*EXP(z*z-(zstar*zstar)) .GT. RandVal3(3)) THEN - EXIT - END IF + ELSE + zstar=beta+(a-beta)*RandVal3(2) + IF ( (a-zstar)/(a-z)*EXP(z*z-(zstar*zstar)) .GT. RandVal3(3)) THEN + EXIT END IF - END DO - CASE(3) - DO - CALL RANDOM_NUMBER(RandVal3) - u = RandVal3(1) - IF ( a*SQRT(PI)/(a*SQRT(PI)+1+a*a) .GT. u) THEN + END IF + END DO + CASE(3) + DO + CALL RANDOM_NUMBER(RandVal3) + u = RandVal3(1) + IF ( a*SQRT(PI)/(a*SQRT(PI)+1+a*a) .GT. u) THEN ! IF (.NOT.DoZigguratSampling) THEN !polar method - IF (RandN_in_Mem) THEN !reusing second RandN form previous polar method - RandN = RandN_save - RandN_in_Mem=.FALSE. - ELSE - Velosq = 2 - DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) - CALL RANDOM_NUMBER(RandVal2) - Velo1 = 2.*RandVal2(1) - 1. - Velo2 = 2.*RandVal2(2) - 1. - Velosq = Velo1**2 + Velo2**2 - END DO - RandN = Velo1*SQRT(-2*LOG(Velosq)/Velosq) - RandN_save = Velo2*SQRT(-2*LOG(Velosq)/Velosq) - RandN_in_Mem=.TRUE. - END IF + IF (RandN_in_Mem) THEN !reusing second RandN form previous polar method + RandN = RandN_save + RandN_in_Mem=.FALSE. + ELSE + Velosq = 2 + DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) + CALL RANDOM_NUMBER(RandVal2) + Velo1 = 2.*RandVal2(1) - 1. + Velo2 = 2.*RandVal2(2) - 1. + Velosq = Velo1**2 + Velo2**2 + END DO + RandN = Velo1*SQRT(-2*LOG(Velosq)/Velosq) + RandN_save = Velo2*SQRT(-2*LOG(Velosq)/Velosq) + RandN_in_Mem=.TRUE. + END IF ! ELSE !ziggurat method ! RandN=rnor() ! END IF - zstar = -1./SQRT(2.)*ABS(RandN) - EXIT - ELSE IF ( (a*SQRT(PI)+1.)/(a*SQRT(PI)+1+a*a) .GT. u) THEN - zstar = -SQRT(-LOG(RandVal3(2))) + zstar = -1./SQRT(2.)*ABS(RandN) + EXIT + ELSE IF ( (a*SQRT(PI)+1.)/(a*SQRT(PI)+1+a*a) .GT. u) THEN + zstar = -SQRT(-LOG(RandVal3(2))) + EXIT + ELSE + zstar = (1.0-SQRT(RandVal3(2)))*a + IF (EXP(-(zstar*zstar)).GT.RandVal3(3)) THEN EXIT - ELSE - zstar = (1.0-SQRT(RandVal3(2)))*a - IF (EXP(-(zstar*zstar)).GT.RandVal3(3)) THEN - EXIT - END IF END IF - END DO - CASE(4) - DO - CALL RANDOM_NUMBER(RandVal3) - IF (1.0/(2.0*a*SQRT(PI)+1.0).GT.RandVal3(1)) THEN - zstar=-SQRT(-LOG(RandVal3(2))) - ELSE + END IF + END DO + CASE(4) + DO + CALL RANDOM_NUMBER(RandVal3) + IF (1.0/(2.0*a*SQRT(PI)+1.0).GT.RandVal3(1)) THEN + zstar=-SQRT(-LOG(RandVal3(2))) + ELSE ! IF (.NOT.DoZigguratSampling) THEN !polar method - IF (RandN_in_Mem) THEN !reusing second RandN form previous polar method - RandN = RandN_save - RandN_in_Mem=.FALSE. - ELSE - Velosq = 2 - DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) - CALL RANDOM_NUMBER(RandVal2) - Velo1 = 2.*RandVal2(1) - 1. - Velo2 = 2.*RandVal2(2) - 1. - Velosq = Velo1**2 + Velo2**2 - END DO - RandN = Velo1*SQRT(-2*LOG(Velosq)/Velosq) - RandN_save = Velo2*SQRT(-2*LOG(Velosq)/Velosq) - RandN_in_Mem=.TRUE. - END IF + IF (RandN_in_Mem) THEN !reusing second RandN form previous polar method + RandN = RandN_save + RandN_in_Mem=.FALSE. + ELSE + Velosq = 2 + DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) + CALL RANDOM_NUMBER(RandVal2) + Velo1 = 2.*RandVal2(1) - 1. + Velo2 = 2.*RandVal2(2) - 1. + Velosq = Velo1**2 + Velo2**2 + END DO + RandN = Velo1*SQRT(-2*LOG(Velosq)/Velosq) + RandN_save = Velo2*SQRT(-2*LOG(Velosq)/Velosq) + RandN_in_Mem=.TRUE. + END IF ! ELSE !ziggurat method ! RandN=rnor() ! END IF - zstar = 1./SQRT(2.)*RandN - END IF - IF ( (a-zstar)/a .GT. RandVal3(3)) THEN - EXIT - END IF - END DO - CASE DEFAULT - CALL abort(__STAMP__,'ERROR in SurfaceFlux: Wrong envelope in SetSurfacefluxVelocities!') - END SELECT - !-- 2.: sample normal directions and build complete velo-vector - Vec3D(1:3) = vec_nIn(1:3) * SQRT(2.*BoltzmannConst*T/Species(DSMC%AmbiDiffElecSpec)%MassIC)*(a-zstar) + zstar = 1./SQRT(2.)*RandN + END IF + IF ( (a-zstar)/a .GT. RandVal3(3)) THEN + EXIT + END IF + END DO + CASE DEFAULT + CALL abort(__STAMP__,'ERROR in SurfaceFlux: Wrong envelope in SetSurfacefluxVelocities!') + END SELECT + !-- 2.: sample normal directions and build complete velo-vector + Vec3D(1:3) = vec_nIn(1:3) * SQRT(2.*BoltzmannConst*T/Species(DSMC%AmbiDiffElecSpec)%MassIC)*(a-zstar) ! IF (.NOT.DoZigguratSampling) THEN !polar method - Velosq = 2 - DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) - CALL RANDOM_NUMBER(RandVal2) - Velo1 = 2.*RandVal2(1) - 1. - Velo2 = 2.*RandVal2(2) - 1. - Velosq = Velo1**2 + Velo2**2 - END DO - Velo1 = Velo1*SQRT(-2*LOG(Velosq)/Velosq) - Velo2 = Velo2*SQRT(-2*LOG(Velosq)/Velosq) + Velosq = 2 + DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) + CALL RANDOM_NUMBER(RandVal2) + Velo1 = 2.*RandVal2(1) - 1. + Velo2 = 2.*RandVal2(2) - 1. + Velosq = Velo1**2 + Velo2**2 + END DO + Velo1 = Velo1*SQRT(-2*LOG(Velosq)/Velosq) + Velo2 = Velo2*SQRT(-2*LOG(Velosq)/Velosq) ! ELSE !ziggurat method ! Velo1=rnor() ! Velo2=rnor() ! END IF - Vec3D(1:3) = Vec3D(1:3) + vec_t1(1:3) * ( Velo_t1+Velo1*SQRT(BoltzmannConst*T/Species(DSMC%AmbiDiffElecSpec)%MassIC) ) - Vec3D(1:3) = Vec3D(1:3) + vec_t2(1:3) * ( Velo_t2+Velo2*SQRT(BoltzmannConst*T/Species(DSMC%AmbiDiffElecSpec)%MassIC) ) - IF (PositionNbr.GT.0) THEN - IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) - ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) - AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = Vec3D(1:3) - END IF - ELSE !PositionNbr .EQ. 0 - CALL abort(__STAMP__,'PositionNbr .EQ. 0!') - END IF !PositionNbr .NE. 0 + Vec3D(1:3) = Vec3D(1:3) + vec_t1(1:3) * ( Velo_t1+Velo1*SQRT(BoltzmannConst*T/Species(DSMC%AmbiDiffElecSpec)%MassIC) ) + Vec3D(1:3) = Vec3D(1:3) + vec_t2(1:3) * ( Velo_t2+Velo2*SQRT(BoltzmannConst*T/Species(DSMC%AmbiDiffElecSpec)%MassIC) ) + IF (ALLOCATED(AmbipolElecVelo(PositionNbr)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo) + ALLOCATE(AmbipolElecVelo(PositionNbr)%ElecVelo(3)) + AmbipolElecVelo(PositionNbr)%ElecVelo(1:3) = Vec3D(1:3) END DO !i = ...NbrOfParticle CASE DEFAULT CALL abort(__STAMP__,'ERROR in SurfaceFlux: Wrong velocity distribution!') END SELECT - + END SUBROUTINE AD_SetSFElectronVelo @@ -459,12 +446,13 @@ SUBROUTINE AD_InsertParticles(iPartIndx_Node, nPart, iPartIndx_NodeTotalAmbi, To !> Creating electrons for each actual ion simulation particle, using the stored velocity vector for the electron !=================================================================================================================================== ! MODULES -USE MOD_Globals +USE MOD_Globals USE MOD_DSMC_Vars ,ONLY: BGGas, CollisMode, DSMC, PartStateIntEn, AmbipolElecVelo, RadialWeighting USE MOD_DSMC_Vars ,ONLY: DSMCSumOfFormedParticles, newAmbiParts, iPartIndx_NodeNewAmbi USE MOD_PARTICLE_Vars ,ONLY: PDM, PartSpecies, PartState, PEM, Species, PartMPF, Symmetry, usevMPF USE MOD_PARTICLE_Vars ,ONLY: UseVarTimeStep, PartTimeStep USE MOD_Particle_Tracking ,ONLY: ParticleInsideCheck +USE MOD_Part_Tools ,ONLY: GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -506,12 +494,7 @@ SUBROUTINE AD_InsertParticles(iPartIndx_Node, nPart, iPartIndx_NodeTotalAmbi, To DO iLoop = 1, nNewElectrons DSMCSumOfFormedParticles = DSMCSumOfFormedParticles + 1 - PositionNbr = PDM%nextFreePosition(DSMCSumOfFormedParticles+PDM%CurrentNextFreePosition) - IF (PositionNbr.EQ.0) THEN - CALL Abort(& -__STAMP__& -,'ERROR in Ambipolar Diffusion: MaxParticleNumber too small!') - END IF + PositionNbr = GetNextFreePosition() InsideFlag=.FALSE. iElem = PEM%GlobalElemID(iPartIndx_Node(1)) DO WHILE(.NOT.InsideFlag) @@ -633,7 +616,7 @@ SUBROUTINE AD_DeleteParticles(iPartIndx_Node, nPart_opt) DO iLoop = 1, nElectron IF (ALLOCATED(AmbipolElecVelo(IonIndX(iLoop))%ElecVelo)) DEALLOCATE(AmbipolElecVelo(IonIndX(iLoop))%ElecVelo) ALLOCATE(AmbipolElecVelo(IonIndX(iLoop))%ElecVelo(3)) - AmbipolElecVelo(IonIndX(iLoop))%ElecVelo(1:3) = PartState(4:6,ElecIndx(iLoop)) + AmbipolElecVelo(IonIndX(iLoop))%ElecVelo(1:3) = PartState(4:6,ElecIndx(iLoop)) PDM%ParticleInside(ElecIndx(iLoop)) = .FALSE. END DO diff --git a/src/particles/dsmc/dsmc_bg_gas.f90 b/src/particles/dsmc/dsmc_bg_gas.f90 index 644042a53..893209ebc 100644 --- a/src/particles/dsmc/dsmc_bg_gas.f90 +++ b/src/particles/dsmc/dsmc_bg_gas.f90 @@ -223,6 +223,7 @@ SUBROUTINE BGGas_InsertParticles() USE MOD_Globals ,ONLY: Abort USE MOD_DSMC_Vars ,ONLY: BGGas USE MOD_PARTICLE_Vars ,ONLY: PDM, PartSpecies, PEM +USE MOD_Part_Tools ,ONLY: GetNextFreePosition #if USE_LOADBALANCE USE MOD_LoadBalance_Timers ,ONLY: LBStartTime,LBPauseTime #endif /*USE_LOADBALANCE*/ @@ -255,10 +256,7 @@ SUBROUTINE BGGas_InsertParticles() END IF ! Get a free particle index iNewPart = iNewPart + 1 - PositionNbr = PDM%nextFreePosition(iNewPart+PDM%CurrentNextFreePosition) - IF (PositionNbr.EQ.0) THEN - CALL Abort(__STAMP__,'ERROR in BGGas: MaxParticleNumber should be increased to account for the BGG particles!') - END IF + PositionNbr = GetNextFreePosition() ! Get the background gas species iSpec = BGGas_GetSpecies(PEM%LocalElemID(iPart)) ! Assign particle properties @@ -272,9 +270,6 @@ SUBROUTINE BGGas_InsertParticles() PEM%pNumber(LocalElemID) = PEM%pNumber(LocalElemID) + 1 END IF END DO -! Increase the particle vector length and update the linked list -PDM%ParticleVecLength = MAX(PDM%ParticleVecLength,PositionNbr) -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + iNewPart #if USE_LOADBALANCE CALL LBPauseTime(LB_DSMC,tLBStart) @@ -592,11 +587,13 @@ SUBROUTINE BGGas_PhotoIonization(iSpec,iInit,TotalNbrOfReactions) USE MOD_part_tools ,ONLY: CalcVelocity_maxwell_particle USE MOD_MCC_Vars ,ONLY: PhotoIonFirstLine,PhotoIonLastLine,PhotoReacToReac,PhotonEnergies USE MOD_MCC_Vars ,ONLY: NbrOfPhotonXsecReactions,SpecPhotonXSecInterpolated,MaxPhotonXSec +USE MOD_Part_Tools ,ONLY: GetNextFreePosition, IncreaseMaxParticleNumber ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -INTEGER, INTENT(IN) :: iSpec,iInit,TotalNbrOfReactions +INTEGER, INTENT(IN) :: iSpec,iInit +INTEGER, INTENT(INOUT) :: TotalNbrOfReactions !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- @@ -706,8 +703,9 @@ SUBROUTINE BGGas_PhotoIonization(iSpec,iInit,TotalNbrOfReactions) !> 2.) Delete left-over inserted particles IF(TotalNbrOfReactions.GT.NbrOfParticle) THEN DO iPart = NbrOfParticle+1,TotalNbrOfReactions - PDM%ParticleInside(PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition)) = .FALSE. + PDM%ParticleInside(GetNextFreePosition(iPart)) = .FALSE. END DO + TotalNbrOfReactions = NbrOfParticle ELSE IF(TotalNbrOfReactions.LT.NbrOfParticle) THEN CALL Abort(__STAMP__,'PhotoIonization: Something is wrong, trying to perform more reactions than anticipated!') END IF @@ -723,14 +721,14 @@ SUBROUTINE BGGas_PhotoIonization(iSpec,iInit,TotalNbrOfReactions) DO iPart = 1, NbrOfParticle ! Loop over the particles with a set position (from SetParticlePosition) - ParticleIndex = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) + ParticleIndex = GetNextFreePosition(iPart) IF (DSMC%DoAmbipolarDiff) THEN newAmbiParts = newAmbiParts + 1 iPartIndx_NodeNewAmbi(newAmbiParts) = ParticleIndex END IF iNewPart = iNewPart + 1 ! Get a new index for the second product - NewParticleIndex = PDM%nextFreePosition(iNewPart+PDM%CurrentNextFreePosition+NbrOfParticle) + NewParticleIndex = GetNextFreePosition(iNewPart+NbrOfParticle) IF (NewParticleIndex.EQ.0) THEN CALL Abort(__STAMP__,'ERROR in PhotoIonization: MaxParticleNumber should be increased!') END IF @@ -788,14 +786,23 @@ SUBROUTINE BGGas_PhotoIonization(iSpec,iInit,TotalNbrOfReactions) IF(DSMC%ElectronicModel.GT.0) PartStateIntEn(3,ParticleIndex) = 0. END DO + ! Add the particles initialized through the emission and the background particles -PDM%ParticleVecLength = PDM%ParticleVecLength + NbrOfParticle + iNewPart ! Update the current next free position -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + NbrOfParticle + iNewPart +IF(iNewPart+NbrOfParticle.GT.0) THEN + PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + NbrOfParticle + iNewPart + PDM%ParticleVecLength = MAX(PDM%ParticleVecLength,GetNextFreePosition(0)) +END IF +#ifdef CODE_ANALYZE +IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL Abort(__STAMP__,'PDM%ParticleVeclength exceeds PDM%maxParticleNumber, Difference:',IntInfoOpt=PDM%ParticleVeclength-PDM%maxParticleNumber) +DO iPart=PDM%ParticleVecLength+1,PDM%maxParticleNumber + IF (PDM%ParticleInside(iPart)) THEN + IPWRITE(*,*) iPart,PDM%ParticleVecLength,PDM%maxParticleNumber + CALL Abort(__STAMP__,'Particle outside PDM%ParticleVeclength',IntInfoOpt=iPart) + END IF +END DO +#endif -IF(PDM%ParticleVecLength.GT.PDM%MaxParticleNumber) CALL Abort(__STAMP__& - ,'ERROR in PhotoIonization: ParticleVecLength greater than MaxParticleNumber! Increase the MaxParticleNumber to at least: ' & - , IntInfoOpt=PDM%ParticleVecLength) !> 4.) Perform the reaction, distribute the collision energy (including photon energy) and emit electrons perpendicular !> to the photon's path @@ -830,10 +837,6 @@ SUBROUTINE BGGas_PhotoIonization(iSpec,iInit,TotalNbrOfReactions) END DO END IF ! NbrOfPhotonXsecReactions.GT.0 -! Advance particle vector length and the current next free position with newly created particles -PDM%ParticleVecLength = PDM%ParticleVecLength + DSMCSumOfFormedParticles -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + DSMCSumOfFormedParticles - DSMCSumOfFormedParticles = 0 DEALLOCATE(Coll_pData) @@ -929,7 +932,8 @@ SUBROUTINE BGGas_TraceSpeciesSplit(iElem, nPart, nPair) USE MOD_Globals USE MOD_DSMC_Vars ,ONLY: BGGas, CollisMode, PartStateIntEn, DSMC USE MOD_DSMC_Vars ,ONLY: DSMC, SpecDSMC, VibQuantsPar, PolyatomMolDSMC -USE MOD_Particle_Vars ,ONLY: PDM,PEM,PartSpecies,PartState,PartMPF,Species +USE MOD_Particle_Vars ,ONLY: PEM,PartSpecies,PartState,PartMPF,Species +USE MOD_Part_Tools ,ONLY: GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -978,10 +982,7 @@ SUBROUTINE BGGas_TraceSpeciesSplit(iElem, nPart, nPair) ! --- Create clone of test particle iNewPart = iNewPart + 1 iSplitPart = iSplitPart + 1 - PartIndex = PDM%nextFreePosition(iNewPart+PDM%CurrentNextFreePosition) - IF (PartIndex.EQ.0) THEN - CALL Abort(__STAMP__,'ERROR in BGGas: MaxParticleNumber should be increased to account for the BGG particles!') - END IF + PartIndex = GetNextFreePosition() ! Assign properties but do not use the velocity and energy of the background gas CALL BGGas_AssignParticleProperties(iSpec,iPart,PartIndex,GetVelocity_opt=.FALSE.,GetInternalEnergy_opt=.TRUE.) ! Copy properties from the particle species @@ -1004,10 +1005,7 @@ SUBROUTINE BGGas_TraceSpeciesSplit(iElem, nPart, nPair) iNewPart = iNewPart + 1 iSplitPart = iSplitPart + 1 ! Get a free particle index - bggPartIndex = PDM%nextFreePosition(iNewPart+PDM%CurrentNextFreePosition) - IF (bggPartIndex.EQ.0) THEN - CALL Abort(__STAMP__,'ERROR in BGGas: MaxParticleNumber should be increased to account for the BGG particles!') - END IF + bggPartIndex = GetNextFreePosition() ! Set the pairing partner BGGas%PairingPartner(PartIndex) = bggPartIndex ! Assign properties of the background gas @@ -1021,9 +1019,6 @@ SUBROUTINE BGGas_TraceSpeciesSplit(iElem, nPart, nPair) END IF iPart = PEM%pNext(iPart) END DO -! Increase the particle vector length and the position in the linked list -PDM%ParticleVecLength = MAX(PDM%ParticleVecLength,bggPartIndex) -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + iNewPart ! Set the new number of particles nPart = PEM%pNumber(iElem) diff --git a/src/particles/dsmc/dsmc_chemical_reactions.f90 b/src/particles/dsmc/dsmc_chemical_reactions.f90 index 0bbbf0003..7ab0abccb 100644 --- a/src/particles/dsmc/dsmc_chemical_reactions.f90 +++ b/src/particles/dsmc/dsmc_chemical_reactions.f90 @@ -359,7 +359,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac) USE MOD_DSMC_CollisVec ,ONLY: PostCollVec USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod USE MOD_Particle_Analyze_Vars ,ONLY: ChemEnergySum -USE MOD_part_tools ,ONLY: GetParticleWeight +USE MOD_part_tools ,ONLY: GetParticleWeight, GetNextFreePosition USE MOD_part_operations ,ONLY: RemoveParticle #ifdef CODE_ANALYZE USE MOD_Globals ,ONLY: unit_stdout,myrank @@ -503,8 +503,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac) IF(ProductReac(3).NE.0) THEN ! === Get free particle index for the 3rd product DSMCSumOfFormedParticles = DSMCSumOfFormedParticles + 1 - ReactInx(3) = PDM%nextFreePosition(DSMCSumOfFormedParticles+PDM%CurrentNextFreePosition) - IF (ReactInx(3).EQ.0) CALL abort(__STAMP__,'New Particle Number greater max Part Num in DSMC_Chemistry. Reaction: ',iReac) + ReactInx(3) = GetNextFreePosition() PDM%ParticleInside(ReactInx(3)) = .true. PDM%IsNewPart(ReactInx(3)) = .true. PDM%dtFracPush(ReactInx(3)) = .FALSE. @@ -543,8 +542,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac) IF(ProductReac(4).NE.0) THEN ! === Get free particle index for the 4th product DSMCSumOfFormedParticles = DSMCSumOfFormedParticles + 1 - ReactInx(4) = PDM%nextFreePosition(DSMCSumOfFormedParticles+PDM%CurrentNextFreePosition) - IF (ReactInx(4).EQ.0) CALL abort(__STAMP__,'New Particle Number greater max Part Num in DSMC_Chemistry. Reaction: ',iReac) + ReactInx(4) = GetNextFreePosition() PDM%ParticleInside(ReactInx(4)) = .true. PDM%IsNewPart(ReactInx(4)) = .true. PDM%dtFracPush(ReactInx(4)) = .FALSE. @@ -1576,6 +1574,7 @@ SUBROUTINE PhotoIonization_InsertProducts(iPair, iReac, iInit, InitSpec, iLineOp USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod USE MOD_Particle_Analyze_Vars ,ONLY: ChemEnergySum USE MOD_part_tools ,ONLY: GetParticleWeight, DiceUnitVector, CalcERot_particle, CalcEVib_particle, CalcEElec_particle +USE MOD_Part_Tools ,ONLY: GetNextFreePosition USE MOD_part_emission_tools ,ONLY: CalcVelocity_maxwell_lpn USE MOD_Particle_Analyze_Vars ,ONLY: CalcPartBalance,nPartIn,PartEkinIn USE MOD_Particle_Analyze_Tools ,ONLY: CalcEkinPart @@ -1639,8 +1638,7 @@ SUBROUTINE PhotoIonization_InsertProducts(iPair, iReac, iInit, InitSpec, iLineOp IF(ProductReac(3).NE.0) THEN ! === Get free particle index for the 3rd product DSMCSumOfFormedParticles = DSMCSumOfFormedParticles + 1 - ReactInx(3) = PDM%nextFreePosition(DSMCSumOfFormedParticles+PDM%CurrentNextFreePosition) - IF (ReactInx(3).EQ.0) CALL abort(__STAMP__,'New Particle Number greater max Part Num in DSMC_Chemistry. Reaction: ',iReac) + ReactInx(3) = GetNextFreePosition() PDM%ParticleInside(ReactInx(3)) = .true. PDM%IsNewPart(ReactInx(3)) = .true. PDM%dtFracPush(ReactInx(3)) = .FALSE. @@ -1670,8 +1668,7 @@ SUBROUTINE PhotoIonization_InsertProducts(iPair, iReac, iInit, InitSpec, iLineOp IF(ProductReac(4).NE.0) THEN ! === Get free particle index for the 4th product DSMCSumOfFormedParticles = DSMCSumOfFormedParticles + 1 - ReactInx(4) = PDM%nextFreePosition(DSMCSumOfFormedParticles+PDM%CurrentNextFreePosition) - IF (ReactInx(4).EQ.0) CALL abort(__STAMP__,'New Particle Number greater max Part Num in DSMC_Chemistry. Reaction: ',iReac) + ReactInx(4) = GetNextFreePosition() PDM%ParticleInside(ReactInx(4)) = .true. PDM%IsNewPart(ReactInx(4)) = .true. PDM%dtFracPush(ReactInx(4)) = .FALSE. diff --git a/src/particles/dsmc/dsmc_init.f90 b/src/particles/dsmc/dsmc_init.f90 index 12c3ad141..e9da13502 100644 --- a/src/particles/dsmc/dsmc_init.f90 +++ b/src/particles/dsmc/dsmc_init.f90 @@ -1456,6 +1456,8 @@ SUBROUTINE FinalizeDSMC() SDEALLOCATE(CollInf%omega) SDEALLOCATE(CollInf%dref) SDEALLOCATE(CollInf%Tref) +SDEALLOCATE(CollInf%OldCollPartner) +CollInf%ProhibitDoubleColl=.FALSE. !SDEALLOCATE(VibQuantsPar) ! SDEALLOCATE(XiEq_Surf) SDEALLOCATE(DSMC_Solution) diff --git a/src/particles/dsmc/dsmc_main.f90 b/src/particles/dsmc/dsmc_main.f90 index ad7b99b70..158da45fe 100644 --- a/src/particles/dsmc/dsmc_main.f90 +++ b/src/particles/dsmc/dsmc_main.f90 @@ -129,10 +129,6 @@ SUBROUTINE DSMC_main(DoElement) END DO ! iElem Loop END IF ! CollisMode.NE.0 -! Advance particle vector length and the current next free position with newly created particles -PDM%ParticleVecLength = PDM%ParticleVecLength + DSMCSumOfFormedParticles -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + DSMCSumOfFormedParticles - IF(PDM%ParticleVecLength.GT.PDM%MaxParticleNumber) THEN CALL Abort(__STAMP__& ,'ERROR in DSMC: ParticleVecLength greater than MaxParticleNumber! Increase the MaxParticleNumber to at least: ' & diff --git a/src/particles/dsmc/dsmc_symmetry.f90 b/src/particles/dsmc/dsmc_symmetry.f90 index 8a68d896d..2aa125fd5 100644 --- a/src/particles/dsmc/dsmc_symmetry.f90 +++ b/src/particles/dsmc/dsmc_symmetry.f90 @@ -574,6 +574,7 @@ SUBROUTINE DSMC_2D_SetInClones() USE MOD_Particle_TimeStep ,ONLY: GetParticleTimeStep USE MOD_TimeDisc_Vars ,ONLY: iter USE MOD_Particle_Analyze_Vars ,ONLY: CalcPartBalance, nPartIn +USE MOD_Part_Tools ,ONLY: GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -614,14 +615,7 @@ SUBROUTINE DSMC_2D_SetInClones() ! 2.) Insert the clones at the position they were created DO iPart = 1, RadialWeighting%ClonePartNum(DelayCounter) - PDM%ParticleVecLength = PDM%ParticleVecLength + 1 - PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + 1 - PositionNbr = PDM%nextFreePosition(PDM%CurrentNextFreePosition) - IF (PDM%ParticleVecLength.GT.PDM%maxParticleNumber) THEN - CALL Abort(& - __STAMP__,& - 'ERROR in 2D axisymmetric simulation: New Particle Number greater max Part Num!') - END IF + PositionNbr = GetNextFreePosition() ! Copy particle parameters PDM%ParticleInside(PositionNbr) = .TRUE. PDM%IsNewPart(PositionNbr) = .TRUE. diff --git a/src/particles/emission/particle_br_electron_fluid.f90 b/src/particles/emission/particle_br_electron_fluid.f90 index 1abef0f12..1f6fbc1db 100644 --- a/src/particles/emission/particle_br_electron_fluid.f90 +++ b/src/particles/emission/particle_br_electron_fluid.f90 @@ -920,7 +920,7 @@ SUBROUTINE CreateElectronsFromBRFluid(CreateFromRestartFile) USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared USE MOD_HDG_Vars ,ONLY: ElemToBRRegion,RegionElectronRef USE MOD_Mesh_Tools ,ONLY: GetCNElemID -USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition +USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition, GetNextFreePosition USE MOD_TimeDisc_Vars ,ONLY: time !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE @@ -1022,13 +1022,7 @@ SUBROUTINE CreateElectronsFromBRFluid(CreateFromRestartFile) DO iPart=1,ElemCharge BRNbrOfElectronsCreated = BRNbrOfElectronsCreated + 1 - ! Set the next free position in the particle vector list - PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + 1 - ParticleIndexNbr = PDM%nextFreePosition(PDM%CurrentNextFreePosition) - IF (ParticleIndexNbr.EQ.0) THEN - CALL abort(__STAMP__,'ERROR in CreateElectronsFromBRFluid(): New Particle Number greater max Part Num!') - END IF - PDM%ParticleVecLength = PDM%ParticleVecLength + 1 + ParticleIndexNbr = GetNextFreePosition() !Set new SpeciesID of new particle (electron) PDM%ParticleInside(ParticleIndexNbr) = .TRUE. diff --git a/src/particles/emission/particle_emission.f90 b/src/particles/emission/particle_emission.f90 index 7f772b69c..f573ed720 100644 --- a/src/particles/emission/particle_emission.f90 +++ b/src/particles/emission/particle_emission.f90 @@ -42,7 +42,7 @@ SUBROUTINE ParticleInserting() USE MOD_Timedisc_Vars ,ONLY: RKdtFrac,RKdtFracTotal USE MOD_Particle_Vars USE MOD_PIC_Vars -USE MOD_part_tools ,ONLY: UpdateNextFreePosition +USE MOD_part_tools ,ONLY: UpdateNextFreePosition, GetNextFreePosition, IncreaseMaxParticleNumber USE MOD_DSMC_Vars ,ONLY: useDSMC, CollisMode, SpecDSMC USE MOD_part_emission_tools ,ONLY: DSMC_SetInternalEnr_LauxVFD USE MOD_DSMC_PolyAtomicModel ,ONLY: DSMC_SetInternalEnr_Poly @@ -333,13 +333,11 @@ SUBROUTINE ParticleInserting() IF (useDSMC.AND.(CollisMode.GT.1)) THEN iPart = 1 DO WHILE (iPart.LE.NbrOfParticle) - PositionNbr = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) - IF (PositionNbr.NE.0) THEN - IF (SpecDSMC(i)%PolyatomicMol) THEN - CALL DSMC_SetInternalEnr_Poly(i,iInit,PositionNbr,1) - ELSE - CALL DSMC_SetInternalEnr_LauxVFD(i,iInit,PositionNbr,1) - END IF + PositionNbr = GetNextFreePosition(iPart) + IF (SpecDSMC(i)%PolyatomicMol) THEN + CALL DSMC_SetInternalEnr_Poly(i,iInit,PositionNbr,1) + ELSE + CALL DSMC_SetInternalEnr_LauxVFD(i,iInit,PositionNbr,1) END IF iPart = iPart + 1 END DO @@ -348,13 +346,24 @@ SUBROUTINE ParticleInserting() IF(CalcPartBalance.AND.(NbrOfParticle.GT.0)) THEN nPartIn(i)=nPartIn(i) + NbrOfparticle DO iPart=1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) - IF (PositionNbr .NE. 0) PartEkinIn(i) = PartEkinIn(i) + CalcEkinPart(PositionNbr) + PositionNbr = GetNextFreePosition(iPart) + PartEkinIn(i) = PartEkinIn(i) + CalcEkinPart(PositionNbr) END DO ! iPart END IF ! CalcPartBalance ! Update the current next free position and increase the particle vector length - PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + NbrOfParticle - PDM%ParticleVecLength = PDM%ParticleVecLength + NbrOfParticle + IF(NbrOfParticle.GT.0) THEN + PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + NbrOfParticle + PDM%ParticleVecLength = MAX(PDM%ParticleVecLength,GetNextFreePosition(0)) + END IF +#ifdef CODE_ANALYZE + IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL Abort(__STAMP__,'PDM%ParticleVeclength exceeds PDM%maxParticleNumber, Difference:',IntInfoOpt=PDM%ParticleVeclength-PDM%maxParticleNumber) + DO iPart=PDM%ParticleVecLength+1,PDM%maxParticleNumber + IF (PDM%ParticleInside(iPart)) THEN + IPWRITE(*,*) iPart,PDM%ParticleVecLength,PDM%maxParticleNumber + CALL Abort(__STAMP__,'Particle outside PDM%ParticleVeclength',IntInfoOpt=iPart) + END IF + END DO +#endif ! Complete check if all particles were emitted successfully #if USE_MPI diff --git a/src/particles/emission/particle_emission_init.f90 b/src/particles/emission/particle_emission_init.f90 index a8dd2dcc6..60e27b945 100644 --- a/src/particles/emission/particle_emission_init.f90 +++ b/src/particles/emission/particle_emission_init.f90 @@ -450,9 +450,10 @@ SUBROUTINE InitialParticleInserting() USE MOD_Part_Emission_Tools ,ONLY: SetParticleChargeAndMass,SetParticleMPF,SetParticleTimeStep USE MOD_Part_Pos_and_Velo ,ONLY: SetParticlePosition,SetParticleVelocity,SetPartPosAndVeloEmissionDistribution USE MOD_DSMC_AmbipolarDiffusion ,ONLY: AD_SetInitElectronVelo -USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition +USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition, IncreaseMaxParticleNumber USE MOD_Particle_Vars ,ONLY: Species,nSpecies,PDM,PEM, usevMPF, SpecReset, UseVarTimeStep USE MOD_Restart_Vars ,ONLY: DoRestart +USE MOD_Part_Tools ,ONLY: GetNextFreePosition #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_LOADBALANCE*/ @@ -500,7 +501,7 @@ SUBROUTINE InitialParticleInserting() IF (useDSMC) THEN IF (DSMC%DoAmbipolarDiff) CALL AD_SetInitElectronVelo(iSpec,iInit,NbrOfParticle) DO iPart = 1, NbrOfParticle - PositionNbr = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr .NE. 0) THEN PDM%PartInit(PositionNbr) = iInit ELSE @@ -509,7 +510,16 @@ SUBROUTINE InitialParticleInserting() END DO END IF ! Add new particles to particle vector length - PDM%ParticleVecLength = PDM%ParticleVecLength + NbrOfParticle + IF(NbrOfParticle.GT.0) PDM%ParticleVecLength = MAX(PDM%ParticleVecLength,GetNextFreePosition(NbrOfParticle)) +#ifdef CODE_ANALYZE + IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL Abort(__STAMP__,'PDM%ParticleVeclength exceeds PDM%maxParticleNumber, Difference:',IntInfoOpt=PDM%ParticleVeclength-PDM%maxParticleNumber) + DO iPart=PDM%ParticleVecLength+1,PDM%maxParticleNumber + IF (PDM%ParticleInside(iPart)) THEN + IPWRITE(*,*) iPart,PDM%ParticleVecLength,PDM%maxParticleNumber + CALL Abort(__STAMP__,'Particle outside PDM%ParticleVeclength',IntInfoOpt=iPart) + END IF + END DO +#endif ! Update CALL UpdateNextFreePosition() END IF ! Species(iSpec)%Init(iInit)%ParticleEmissionType.EQ.0 @@ -751,7 +761,7 @@ SUBROUTINE DetermineInitialParticleNumber() USE MOD_Globals_Vars ,ONLY: PI USE MOD_DSMC_Vars ,ONLY: RadialWeighting, DSMC USE MOD_Particle_Mesh_Vars ,ONLY: LocalVolume -USE MOD_Particle_Vars ,ONLY: PDM,Species,nSpecies,SpecReset,Symmetry +USE MOD_Particle_Vars ,ONLY: Species,nSpecies,SpecReset,Symmetry USE MOD_ReadInTools USE MOD_Restart_Vars ,ONLY: DoRestart ! IMPLICIT VARIABLE HANDLING @@ -837,14 +847,6 @@ SUBROUTINE DetermineInitialParticleNumber() END DO ! iInit = 1, Species(iSpec)%NumberOfInits END DO ! iSpec=1,nSpecies -IF(.NOT.RadialWeighting%DoRadialWeighting) THEN - IF (insertParticles.GT.PDM%maxParticleNumber) THEN - IPWRITE(UNIT_stdOut,*)' Maximum particle number : ',PDM%maxParticleNumber - IPWRITE(UNIT_stdOut,*)' To be inserted particles: ',INT(insertParticles,4) - CALL abort(__STAMP__,'Number of to be inserted particles per init-proc exceeds max. particle number! ') - END IF -END IF - END SUBROUTINE DetermineInitialParticleNumber diff --git a/src/particles/emission/particle_emission_tools.f90 b/src/particles/emission/particle_emission_tools.f90 index f73b56709..8467606a9 100644 --- a/src/particles/emission/particle_emission_tools.f90 +++ b/src/particles/emission/particle_emission_tools.f90 @@ -210,8 +210,9 @@ SUBROUTINE SetParticleTimeStep(NbrOfParticle) !> the particle vector, loops over the total number of particles and the indices in the nextFreePosition array. !=================================================================================================================================== ! MODULES -USE MOD_Particle_Vars ,ONLY: PDM, PartTimeStep, PEM, PartState +USE MOD_Particle_Vars ,ONLY: PartTimeStep, PEM, PartState USE MOD_Particle_TimeStep ,ONLY: GetParticleTimeStep +USE MOD_Part_Tools ,ONLY: GetNextFreePosition !---------------------------------------------------------------------------------------------------------------------------------- ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -225,7 +226,7 @@ SUBROUTINE SetParticleTimeStep(NbrOfParticle) INTEGER :: iPart, PositionNbr !=================================================================================================================================== DO iPart=1, NbrOfParticle - PositionNbr = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) + PositionNbr = GetNextFreePosition(iPart) PartTimeStep(PositionNbr) = GetParticleTimeStep(PartState(1,PositionNbr), PartState(2,PositionNbr),PEM%LocalElemID(PositionNbr)) END DO @@ -238,7 +239,8 @@ SUBROUTINE SetParticleChargeAndMass(FractNbr,NbrOfParticle) !=================================================================================================================================== ! MODULES USE MOD_Globals -USE MOD_Particle_Vars, ONLY : PDM, PartSpecies +USE MOD_Particle_Vars ,ONLY: PartSpecies +USE MOD_Part_Tools ,ONLY: GetNextFreePosition !---------------------------------------------------------------------------------------------------------------------------------- ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -250,15 +252,11 @@ SUBROUTINE SetParticleChargeAndMass(FractNbr,NbrOfParticle) INTEGER,INTENT(INOUT) :: NbrOfParticle !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: i,PositionNbr +INTEGER :: iPart,PositionNbr !=================================================================================================================================== -DO i=1, NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr .NE. 0) THEN - PartSpecies(PositionNbr) = FractNbr - ELSE - CALL abort(__STAMP__,'ERROR in SetParticlePosition:ParticleIndexNbr.EQ.0 - maximum nbr of particles reached?') - END IF +DO iPart=1, NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) + PartSpecies(PositionNbr) = FractNbr END DO END SUBROUTINE SetParticleChargeAndMass @@ -270,9 +268,9 @@ SUBROUTINE SetParticleMPF(FractNbr,iInit,NbrOfParticle) !=================================================================================================================================== ! MODULES USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: PDM, PartMPF, Species, PartState +USE MOD_Particle_Vars ,ONLY: PartMPF, Species, PartState USE MOD_DSMC_Vars ,ONLY: RadialWeighting -USE MOD_part_tools ,ONLY: CalcRadWeightMPF +USE MOD_part_tools ,ONLY: CalcRadWeightMPF, GetNextFreePosition !=================================================================================================================================== ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -285,25 +283,19 @@ SUBROUTINE SetParticleMPF(FractNbr,iInit,NbrOfParticle) INTEGER,INTENT(INOUT) :: NbrOfParticle !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: i,PositionNbr -!=================================================================================================================================== -i = 1 -DO WHILE (i .le. NbrOfParticle) - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr .NE. 0) THEN - IF(RadialWeighting%DoRadialWeighting) THEN - PartMPF(PositionNbr) = CalcRadWeightMPF(PartState(2,PositionNbr),FractNbr,PositionNbr) - ELSE - IF(iInit.EQ.-1)THEN - PartMPF(PositionNbr) = Species(FractNbr)%MacroParticleFactor - ELSE - PartMPF(PositionNbr) = Species(FractNbr)%Init(iInit)%MacroParticleFactor ! Use emission-specific MPF (default is species MPF) - END IF ! iInit.EQ.-1 - END IF +INTEGER :: iPart,PositionNbr +!=================================================================================================================================== +DO iPart=1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) + IF(RadialWeighting%DoRadialWeighting) THEN + PartMPF(PositionNbr) = CalcRadWeightMPF(PartState(2,PositionNbr),FractNbr,PositionNbr) ELSE - CALL abort(__STAMP__,'ERROR in SetParticlePosition:ParticleIndexNbr.EQ.0 - maximum nbr of particles reached?') + IF(iInit.EQ.-1)THEN + PartMPF(PositionNbr) = Species(FractNbr)%MacroParticleFactor + ELSE + PartMPF(PositionNbr) = Species(FractNbr)%Init(iInit)%MacroParticleFactor ! Use emission-specific MPF (default is species MPF) + END IF ! iInit.EQ.-1 END IF - i = i + 1 END DO END SUBROUTINE SetParticleMPF @@ -919,6 +911,7 @@ SUBROUTINE SetCellLocalParticlePosition(chunkSize,iSpec,iInit,UseExactPartNum) USE MOD_Particle_Tracking ,ONLY: ParticleInsideCheck USE MOD_Particle_Vars ,ONLY: Species, PDM, PartState, PEM, Symmetry, UseVarTimeStep, PartTimeStep, PartMPF USE MOD_Particle_TimeStep ,ONLY: GetParticleTimeStep +USE MOD_Part_Tools ,ONLY: IncreaseMaxParticleNumber, GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -944,9 +937,7 @@ SUBROUTINE SetCellLocalParticlePosition(chunkSize,iSpec,iInit,UseExactPartNum) INTEGER :: CNElemID !----------------------------------------------------------------------------------------------------------------------------------- IF (UseExactPartNum) THEN - IF(chunkSize.GE.PDM%maxParticleNumber) THEN - CALL abort(__STAMP__,'SetCellLocalParticlePosition: Maximum particle number reached! max. particles needed: ',chunksize) - END IF + IF (Species(iSpec)%Init(iInit)%ParticleEmissionType.EQ.0) CALL IncreaseMaxParticleNumber(chunkSize) CellChunkSize(:)=0 ASSOCIATE( start => GetCNElemID(1+offsetElem),& end => GetCNElemID(nElems+offsetElem)) @@ -956,11 +947,7 @@ SUBROUTINE SetCellLocalParticlePosition(chunkSize,iSpec,iInit,UseExactPartNum) PartDens = Species(iSpec)%Init(iInit)%PartDensity / Species(iSpec)%MacroParticleFactor ! numerical Partdensity is needed IF(RadialWeighting%DoRadialWeighting) PartDens = PartDens * 2. / (RadialWeighting%PartScaleFactor) chunkSize_tmp = INT(PartDens * LocalVolume) - IF(chunkSize_tmp.GE.PDM%maxParticleNumber) THEN - CALL abort(__STAMP__,& - 'ERROR in SetCellLocalParticlePosition: Maximum particle number during sanity check! max. particles needed: ',& - IntInfoOpt=chunkSize_tmp) - END IF + IF (Species(iSpec)%Init(iInit)%ParticleEmissionType.EQ.0) CALL IncreaseMaxParticleNumber(chunkSize_tmp) END IF ichunkSize = 1 @@ -984,42 +971,33 @@ SUBROUTINE SetCellLocalParticlePosition(chunkSize,iSpec,iInit,UseExactPartNum) END IF END IF DO iPart = 1, nPart - ParticleIndexNbr = PDM%nextFreePosition(iChunksize + PDM%CurrentNextFreePosition) - IF (ParticleIndexNbr .ne. 0) THEN - InsideFlag=.FALSE. - DO WHILE(.NOT.InsideFlag) - CALL RANDOM_NUMBER(RandomPos) - IF(Symmetry%Axisymmetric.AND.(.NOT.RadialWeighting%DoRadialWeighting)) THEN - ! Treatment of axisymmetry without weighting - RandomPos(1) = Bounds(1,1) + RandomPos(1)*(Bounds(2,1)-Bounds(1,1)) - RandomPos(2) = SQRT(RandomPos(2)*(Bounds(2,2)**2-Bounds(1,2)**2)+Bounds(1,2)**2) - ELSE - RandomPos = Bounds(1,:) + RandomPos*(Bounds(2,:)-Bounds(1,:)) - END IF - IF(Symmetry%Order.LE.2) RandomPos(3) = 0. - IF(Symmetry%Order.LE.1) RandomPos(2) = 0. - InsideFlag = ParticleInsideCheck(RandomPos,iPart,iGlobalElem) - END DO - PartState(1:3,ParticleIndexNbr) = RandomPos(1:3) - PDM%ParticleInside(ParticleIndexNbr) = .TRUE. - PDM%IsNewPart(ParticleIndexNbr)=.TRUE. - PDM%dtFracPush(ParticleIndexNbr) = .FALSE. - PEM%GlobalElemID(ParticleIndexNbr) = iGlobalElem - ichunkSize = ichunkSize + 1 - IF (UseVarTimeStep) THEN - PartTimeStep(ParticleIndexNbr) = & - GetParticleTimeStep(PartState(1,ParticleIndexNbr), PartState(2,ParticleIndexNbr),iElem) - END IF - IF(RadialWeighting%DoRadialWeighting) THEN - PartMPF(ParticleIndexNbr) = CalcRadWeightMPF(PartState(2,ParticleIndexNbr),iSpec,ParticleIndexNbr) + ParticleIndexNbr = GetNextFreePosition(ichunkSize) + InsideFlag=.FALSE. + DO WHILE(.NOT.InsideFlag) + CALL RANDOM_NUMBER(RandomPos) + IF(Symmetry%Axisymmetric.AND.(.NOT.RadialWeighting%DoRadialWeighting)) THEN + ! Treatment of axisymmetry without weighting + RandomPos(1) = Bounds(1,1) + RandomPos(1)*(Bounds(2,1)-Bounds(1,1)) + RandomPos(2) = SQRT(RandomPos(2)*(Bounds(2,2)**2-Bounds(1,2)**2)+Bounds(1,2)**2) + ELSE + RandomPos = Bounds(1,:) + RandomPos*(Bounds(2,:)-Bounds(1,:)) END IF - ELSE - WRITE(UNIT_stdOut,*) "" - IPWRITE(UNIT_stdOut,*) "ERROR:" - IPWRITE(UNIT_stdOut,*) " iPart :", iPart - IPWRITE(UNIT_stdOut,*) "PDM%maxParticleNumber :", PDM%maxParticleNumber - CALL abort(__STAMP__& - ,'ERROR in SetCellLocalParticlePosition: Maximum particle number reached during inserting! --> ParticleIndexNbr.EQ.0') + IF(Symmetry%Order.LE.2) RandomPos(3) = 0. + IF(Symmetry%Order.LE.1) RandomPos(2) = 0. + InsideFlag = ParticleInsideCheck(RandomPos,iPart,iGlobalElem) + END DO + PartState(1:3,ParticleIndexNbr) = RandomPos(1:3) + PDM%ParticleInside(ParticleIndexNbr) = .TRUE. + PDM%IsNewPart(ParticleIndexNbr)=.TRUE. + PDM%dtFracPush(ParticleIndexNbr) = .FALSE. + PEM%GlobalElemID(ParticleIndexNbr) = iGlobalElem + ichunkSize = ichunkSize + 1 + IF (UseVarTimeStep) THEN + PartTimeStep(ParticleIndexNbr) = & + GetParticleTimeStep(PartState(1,ParticleIndexNbr), PartState(2,ParticleIndexNbr),iElem) + END IF + IF(RadialWeighting%DoRadialWeighting) THEN + PartMPF(ParticleIndexNbr) = CalcRadWeightMPF(PartState(2,ParticleIndexNbr),iSpec,ParticleIndexNbr) END IF END DO END ASSOCIATE @@ -1307,7 +1285,7 @@ SUBROUTINE SetParticlePositionCuboidCylinder(FractNbr,iInit,chunkSize,particle_p !=================================================================================================================================== ! modules USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Species, Symmetry, PDM +USE MOD_Particle_Vars ,ONLY: Species, Symmetry USE MOD_Part_Tools ,ONLY: CalcPartSymmetryPos, CalcRadWeightMPF USE MOD_DSMC_Vars ,ONLY: RadialWeighting !---------------------------------------------------------------------------------------------------------------------------------- @@ -1372,10 +1350,6 @@ SUBROUTINE SetParticlePositionCuboidCylinder(FractNbr,iInit,chunkSize,particle_p IF(Species(FractNbr)%MacroParticleFactor/RadWeightMPF.LT.iRan) THEN i=i+1 CYCLE - ELSE IF(chunkSize2.GT.PDM%maxParticleNumber) THEN - IPWRITE(UNIT_stdOut,*)'Inserted percentage of particles',REAL(i)/REAL(chunkSize)*100 - CALL CollectiveStop(__STAMP__,& - 'Number of to be inserted particles per init-proc exceeds max. particle number! ') END IF END IF END IF @@ -1395,7 +1369,7 @@ SUBROUTINE SetParticlePositionSphere(FractNbr,iInit,chunkSize,particle_positions !=================================================================================================================================== ! modules USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Species, Symmetry, PDM +USE MOD_Particle_Vars ,ONLY: Species, Symmetry USE MOD_Part_tools ,ONLY: DICEUNITVECTOR, CalcPartSymmetryPos, CalcRadWeightMPF USE MOD_DSMC_Vars ,ONLY: RadialWeighting !---------------------------------------------------------------------------------------------------------------------------------- @@ -1434,10 +1408,6 @@ SUBROUTINE SetParticlePositionSphere(FractNbr,iInit,chunkSize,particle_positions IF(Species(FractNbr)%MacroParticleFactor/RadWeightMPF.LT.iRan) THEN i=i+1 CYCLE - ELSE IF(chunkSize2.GT.PDM%maxParticleNumber) THEN - IPWRITE(UNIT_stdOut,*)'Inserted percentage of particles',REAL(i)/REAL(chunkSize)*100 - CALL CollectiveStop(__STAMP__,& - 'Number of to be inserted particles per init-proc exceeds max. particle number! ') END IF END IF END IF @@ -2323,7 +2293,7 @@ SUBROUTINE SetParticlePositionLiu2010SzaboNeutralization(chunkSize,particle_posi ! Count number of emitted particles to compare with chunkSize later on emittedParticles = emittedParticles + 1 ! Emit at random position in element (assume tri-linear element geometry, if position is outside discard the position) - ASSOCIATE( Bounds => BoundsOfElem_Shared(1:2,1:3,GlobalElemID) ) ! 1-2: Min, Max value; 1-3: x,y,z + ASSOCIATE( Bounds => BoundsOfElem_Shared(1:2,1:3,GlobalElemID) ) ! 1-2: Min, Max value; 1-3: x,y,z InsideFlag = .FALSE. DO WHILE(.NOT.InsideFlag) CALL RANDOM_NUMBER(RandomPos) diff --git a/src/particles/emission/particle_macroscopic_restart.f90 b/src/particles/emission/particle_macroscopic_restart.f90 index 4830fc115..fb4c0d41a 100644 --- a/src/particles/emission/particle_macroscopic_restart.f90 +++ b/src/particles/emission/particle_macroscopic_restart.f90 @@ -38,12 +38,12 @@ SUBROUTINE MacroRestart_InsertParticles() USE MOD_Globals USE MOD_Globals_Vars ,ONLY: Pi USE MOD_DSMC_Vars ,ONLY: RadialWeighting, DSMC -USE MOD_part_tools ,ONLY: CalcRadWeightMPF,InitializeParticleMaxwell +USE MOD_part_tools ,ONLY: CalcRadWeightMPF,InitializeParticleMaxwell, IncreaseMaxParticleNumber USE MOD_Mesh_Vars ,ONLY: nElems,offsetElem USE MOD_Particle_TimeStep ,ONLY: GetParticleTimeStep USE MOD_Particle_Vars ,ONLY: Species, PDM, nSpecies, PartState, Symmetry, UseVarTimeStep USE MOD_Restart_Vars ,ONLY: MacroRestartValues -USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared,BoundsOfElem_Shared +USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared,BoundsOfElem_Shared, ElemMidPoint_Shared USE MOD_Particle_Tracking ,ONLY: ParticleInsideCheck !----------------------------------------------------------------------------------------------------------------------------------- ! IMPLICIT VARIABLE HANDLING @@ -57,7 +57,7 @@ SUBROUTINE MacroRestart_InsertParticles() !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: iElem,iSpec,iPart,nPart,locnPart,iHeight,yPartitions,GlobalElemID -REAL :: iRan, RandomPos(3), PartDens, TempMPF, MaxPosTemp, MinPosTemp +REAL :: iRan, RandomPos(3), PartDens, MaxPosTemp, MinPosTemp REAL :: TempVol, Volume LOGICAL :: InsideFlag !=================================================================================================================================== @@ -77,26 +77,26 @@ SUBROUTINE MacroRestart_InsertParticles() IF (iSpec.EQ.DSMC%AmbiDiffElecSpec) CYCLE END IF yPartitions = 6 - PartDens = MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) ! Particle weighting DO iHeight = 1, yPartitions MinPosTemp = Bounds(1,2) + (Bounds(2,2) - Bounds(1,2))/ yPartitions *(iHeight-1.) MaxPosTemp = Bounds(1,2) + (Bounds(2,2) - Bounds(1,2))/ yPartitions *iHeight TempVol = (MaxPosTemp-MinPosTemp)*(Bounds(2,1)-Bounds(1,1)) * Pi * (MaxPosTemp+MinPosTemp) - TempMPF = CalcRadWeightMPF((MaxPosTemp+MinPosTemp)*0.5,iSpec) + PartDens = MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / CalcRadWeightMPF((MaxPosTemp+MinPosTemp)*0.5,iSpec) IF(UseVarTimeStep) THEN - TempMPF = TempMPF * GetParticleTimeStep((Bounds(2,1)+Bounds(1,1))*0.5, (MaxPosTemp+MinPosTemp)*0.5, iElem) + PartDens = PartDens * GetParticleTimeStep((Bounds(2,1)+Bounds(1,1))*0.5, (MaxPosTemp+MinPosTemp)*0.5, iElem) END IF CALL RANDOM_NUMBER(iRan) - nPart = INT(PartDens / TempMPF * TempVol + iRan) + nPart = INT(PartDens * TempVol + iRan) DO iPart = 1, nPart InsideFlag=.FALSE. CALL RANDOM_NUMBER(RandomPos) RandomPos(1) = Bounds(1,1) + RandomPos(1)*(Bounds(2,1)-Bounds(1,1)) RandomPos(2) = MinPosTemp + RandomPos(2)*(MaxPosTemp-MinPosTemp) RandomPos(3) = 0.0 - InsideFlag = ParticleInsideCheck(RandomPos,iPart,GlobalElemID) + InsideFlag = ParticleInsideCheck(RandomPos,1,GlobalElemID) IF (InsideFlag) THEN + IF (locnPart.GE.PDM%maxParticleNumber) CALL IncreaseMaxParticleNumber() PartState(1:3,locnPart) = RandomPos(1:3) CALL InitializeParticleMaxwell(locnPart,iSpec,iElem,Mode=1) locnPart = locnPart + 1 @@ -109,12 +109,13 @@ SUBROUTINE MacroRestart_InsertParticles() IF (DSMC%DoAmbipolarDiff) THEN IF (iSpec.EQ.DSMC%AmbiDiffElecSpec) CYCLE END IF + PartDens = MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / Species(iSpec)%MacroParticleFactor CALL RANDOM_NUMBER(iRan) - TempMPF = Species(iSpec)%MacroParticleFactor IF(UseVarTimeStep) THEN - TempMPF = TempMPF * GetParticleTimeStep((Bounds(2,1)+Bounds(1,1))*0.5, (Bounds(2,2)+Bounds(1,2))*0.5, iElem) + PartDens = PartDens / GetParticleTimeStep(ElemMidPoint_Shared(1,GlobalElemID), ElemMidPoint_Shared(2,GlobalElemID), iElem) END IF - nPart = INT(MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / TempMPF * ElemVolume_Shared(GlobalElemID) + iRan) + nPart = INT(PartDens * ElemVolume_Shared(GlobalElemID) + iRan) + CALL IncreaseMaxParticleNumber(nPart) DO iPart = 1, nPart InsideFlag=.FALSE. DO WHILE (.NOT.InsideFlag) @@ -136,23 +137,25 @@ SUBROUTINE MacroRestart_InsertParticles() IF (DSMC%DoAmbipolarDiff) THEN IF (iSpec.EQ.DSMC%AmbiDiffElecSpec) CYCLE END IF + PartDens = MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / Species(iSpec)%MacroParticleFactor CALL RANDOM_NUMBER(iRan) - TempMPF = Species(iSpec)%MacroParticleFactor IF(UseVarTimeStep) THEN - TempMPF = TempMPF * GetParticleTimeStep((Bounds(2,1)+Bounds(1,1))*0.5, (Bounds(2,2)+Bounds(1,2))*0.5, iElem) + PartDens = PartDens / GetParticleTimeStep(ElemMidPoint_Shared(1,GlobalElemID), ElemMidPoint_Shared(2,GlobalElemID), iElem) END IF - nPart = INT(MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / TempMPF * Volume + iRan) + nPart = INT(PartDens * ElemVolume_Shared(GlobalElemID) + iRan) + CALL IncreaseMaxParticleNumber(nPart) + DO iPart = 1, nPart InsideFlag=.FALSE. - CALL RANDOM_NUMBER(RandomPos(1:2)) - RandomPos(1:2) = Bounds(1,1:2) + RandomPos(1:2)*(Bounds(2,1:2)-Bounds(1,1:2)) - RandomPos(3) = 0.0 - InsideFlag = ParticleInsideCheck(RandomPos,iPart,GlobalElemID) - IF (InsideFlag) THEN - PartState(1:3,locnPart) = RandomPos(1:3) - CALL InitializeParticleMaxwell(locnPart,iSpec,iElem,Mode=1) - locnPart = locnPart + 1 - END IF + DO WHILE(.NOT.InsideFlag) + CALL RANDOM_NUMBER(RandomPos(1:2)) + RandomPos(1:2) = Bounds(1,1:2) + RandomPos(1:2)*(Bounds(2,1:2)-Bounds(1,1:2)) + RandomPos(3) = 0.0 + InsideFlag = ParticleInsideCheck(RandomPos,iPart,GlobalElemID) + END DO + PartState(1:3,locnPart) = RandomPos(1:3) + CALL InitializeParticleMaxwell(locnPart,iSpec,iElem,Mode=1) + locnPart = locnPart + 1 END DO ! nPart END DO ! nSpecies ELSE IF(Symmetry%Order.EQ.1) THEN @@ -161,24 +164,25 @@ SUBROUTINE MacroRestart_InsertParticles() IF (DSMC%DoAmbipolarDiff) THEN IF (iSpec.EQ.DSMC%AmbiDiffElecSpec) CYCLE END IF + PartDens = MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / Species(iSpec)%MacroParticleFactor CALL RANDOM_NUMBER(iRan) - TempMPF = Species(iSpec)%MacroParticleFactor IF(UseVarTimeStep) THEN - TempMPF = TempMPF * GetParticleTimeStep((Bounds(2,1)+Bounds(1,1))*0.5, (Bounds(2,2)+Bounds(1,2))*0.5, iElem) + PartDens = PartDens / GetParticleTimeStep(ElemMidPoint_Shared(1,GlobalElemID), ElemMidPoint_Shared(2,GlobalElemID), iElem) END IF - nPart = INT(MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / TempMPF * Volume + iRan) + nPart = INT(PartDens * ElemVolume_Shared(GlobalElemID) + iRan) + CALL IncreaseMaxParticleNumber(nPart) DO iPart = 1, nPart InsideFlag=.FALSE. - CALL RANDOM_NUMBER(RandomPos(1)) - RandomPos(1:2) = Bounds(1,1) + RandomPos(1)*(Bounds(2,1)-Bounds(1,1)) - RandomPos(2) = 0.0 - RandomPos(3) = 0.0 - InsideFlag = ParticleInsideCheck(RandomPos,iPart,GlobalElemID) - IF (InsideFlag) THEN - PartState(1:3,locnPart) = RandomPos(1:3) - CALL InitializeParticleMaxwell(locnPart,iSpec,iElem,Mode=1) - locnPart = locnPart + 1 - END IF + DO WHILE(.NOT.InsideFlag) + CALL RANDOM_NUMBER(RandomPos(1)) + RandomPos(1:2) = Bounds(1,1) + RandomPos(1)*(Bounds(2,1)-Bounds(1,1)) + RandomPos(2) = 0.0 + RandomPos(3) = 0.0 + InsideFlag = ParticleInsideCheck(RandomPos,iPart,GlobalElemID) + END DO + PartState(1:3,locnPart) = RandomPos(1:3) + CALL InitializeParticleMaxwell(locnPart,iSpec,iElem,Mode=1) + locnPart = locnPart + 1 END DO ! nPart END DO ! nSpecies ELSE @@ -188,31 +192,40 @@ SUBROUTINE MacroRestart_InsertParticles() IF (DSMC%DoAmbipolarDiff) THEN IF (iSpec.EQ.DSMC%AmbiDiffElecSpec) CYCLE END IF + PartDens = MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / Species(iSpec)%MacroParticleFactor CALL RANDOM_NUMBER(iRan) - TempMPF = Species(iSpec)%MacroParticleFactor IF(UseVarTimeStep) THEN - TempMPF = TempMPF * GetParticleTimeStep(iElem=iElem) + PartDens = PartDens / GetParticleTimeStep(ElemMidPoint_Shared(1,GlobalElemID), ElemMidPoint_Shared(2,GlobalElemID), iElem) END IF - nPart = INT(MacroRestartValues(iElem,iSpec,DSMC_NUMDENS) / TempMPF * Volume + iRan) + nPart = INT(PartDens * ElemVolume_Shared(GlobalElemID) + iRan) + CALL IncreaseMaxParticleNumber(nPart) DO iPart = 1, nPart InsideFlag=.FALSE. - CALL RANDOM_NUMBER(RandomPos) - RandomPos(1:3) = Bounds(1,1:3) + RandomPos(1:3)*(Bounds(2,1:3)-Bounds(1,1:3)) - InsideFlag = ParticleInsideCheck(RandomPos,iPart,GlobalElemID) - IF (InsideFlag) THEN - PartState(1:3,locnPart) = RandomPos(1:3) - CALL InitializeParticleMaxwell(locnPart,iSpec,iElem,Mode=1) - locnPart = locnPart + 1 - END IF + DO WHILE(.NOT.InsideFlag) + CALL RANDOM_NUMBER(RandomPos) + RandomPos(1:3) = Bounds(1,1:3) + RandomPos(1:3)*(Bounds(2,1:3)-Bounds(1,1:3)) + InsideFlag = ParticleInsideCheck(RandomPos,iPart,GlobalElemID) + END DO + PartState(1:3,locnPart) = RandomPos(1:3) + CALL InitializeParticleMaxwell(locnPart,iSpec,iElem,Mode=1) + locnPart = locnPart + 1 END DO ! nPart END DO ! nSpecies END IF ! 1D/2D/Axisymmetric/3D END ASSOCIATE END DO ! nElems -IF(locnPart.GE.PDM%maxParticleNumber) CALL abort(__STAMP__,'ERROR in MacroRestart: Increase maxParticleNumber!', locnPart) - +! IF(locnPart.GE.PDM%maxParticleNumber) CALL abort(__STAMP__,'ERROR in MacroRestart: Increase maxParticleNumber!', locnPart) PDM%ParticleVecLength = PDM%ParticleVecLength + locnPart +#ifdef CODE_ANALYZE +IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL Abort(__STAMP__,'PDM%ParticleVeclength exceeds PDM%maxParticleNumber, Difference:',IntInfoOpt=PDM%ParticleVeclength-PDM%maxParticleNumber) +DO iPart=PDM%ParticleVecLength+1,PDM%maxParticleNumber + IF (PDM%ParticleInside(iPart)) THEN + IPWRITE(*,*) iPart,PDM%ParticleVecLength,PDM%maxParticleNumber + CALL Abort(__STAMP__,'Particle outside PDM%ParticleVeclength',IntInfoOpt=iPart) + END IF +END DO +#endif END SUBROUTINE MacroRestart_InsertParticles diff --git a/src/particles/emission/particle_position_and_velocity.f90 b/src/particles/emission/particle_position_and_velocity.f90 index a67d4a9cc..ab9d7ebd3 100644 --- a/src/particles/emission/particle_position_and_velocity.f90 +++ b/src/particles/emission/particle_position_and_velocity.f90 @@ -134,8 +134,8 @@ SUBROUTINE SetParticlePosition(FractNbr,iInit,NbrOfParticle) !=================================================================================================================================== ! modules USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Species,PDM,PartState,FractNbrOld,chunkSizeOld,NeutralizationBalance -USE MOD_Particle_Localization ,ONLY: LocateParticleInElement +USE MOD_Particle_Vars ,ONLY: Species,PDM,PartState,FractNbrOld,chunkSizeOld,NeutralizationBalance, PartPosRef, PEM +USE MOD_Particle_Localization ,ONLY: SinglePointToElement USE MOD_part_emission_tools ,ONLY: IntegerDivide,SetCellLocalParticlePosition,SetParticlePositionPoint USE MOD_part_emission_tools ,ONLY: SetParticlePositionEquidistLine, SetParticlePositionLine, SetParticlePositionDisk USE MOD_part_emission_tools ,ONLY: SetParticlePositionCuboidCylinder, SetParticlePositionGyrotronCircle,SetParticlePositionCircle @@ -146,7 +146,9 @@ SUBROUTINE SetParticlePosition(FractNbr,iInit,NbrOfParticle) USE MOD_part_emission_tools ,ONLY: SetParticlePositionLandmark,SetParticlePositionLandmarkNeutralization USE MOD_part_emission_tools ,ONLY: SetParticlePositionLiu2010Neutralization,SetParticlePositionLiu2010Neutralization3D USE MOD_part_emission_tools ,ONLY: SetParticlePositionLiu2010SzaboNeutralization -USE MOD_DSMC_Vars ,ONLY: RadialWeighting +USE MOD_Eval_xyz ,ONLY: GetPositionInRefElem +USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod +USE MOD_Part_Tools ,ONLY: IncreaseMaxParticleNumber, GetNextFreePosition #if USE_MPI USE MOD_Particle_MPI_Emission ,ONLY: SendEmissionParticlesToProcs USE MOD_Particle_MPI_Vars ,ONLY: PartMPIInitGroup @@ -165,6 +167,7 @@ SUBROUTINE SetParticlePosition(FractNbr,iInit,NbrOfParticle) REAL,ALLOCATABLE :: particle_positions(:) INTEGER :: i,ParticleIndexNbr,allocStat,nChunks, chunkSize INTEGER :: DimSend +INTEGER, ALLOCATABLE :: AcceptedParts(:) #if USE_MPI INTEGER :: InitGroup #endif @@ -221,15 +224,9 @@ SUBROUTINE SetParticlePosition(FractNbr,iInit,NbrOfParticle) IF (PartMPIInitGroup(InitGroup)%MPIROOT.OR.nChunks.GT.1) THEN #endif ! Allocate part pos buffer - IF(RadialWeighting%DoRadialWeighting.AND.(chunkSize.GT.PDM%maxParticleNumber)) THEN - ALLOCATE( particle_positions(1:PDM%maxParticleNumber*DimSend), STAT=allocStat ) - IF (allocStat .NE. 0) & - CALL abort(__STAMP__,'ERROR in SetParticlePosition: cannot allocate particle_positions!') - ELSE - ALLOCATE( particle_positions(1:chunkSize*DimSend), STAT=allocStat ) - IF (allocStat .NE. 0) & - CALL abort(__STAMP__,'ERROR in SetParticlePosition: cannot allocate particle_positions!') - END IF + ALLOCATE( particle_positions(1:chunkSize*DimSend), STAT=allocStat ) + IF (allocStat .NE. 0) & + CALL abort(__STAMP__,'ERROR in SetParticlePosition: cannot allocate particle_positions!') ! Sanity check IF (allocStat .NE. 0) CALL abort(__STAMP__,'ERROR in SetParticlePosition: cannot allocate particle_positions!') @@ -301,26 +298,30 @@ SUBROUTINE SetParticlePosition(FractNbr,iInit,NbrOfParticle) ! Finish emission on local proc ELSE #endif /*USE_MPI*/ - Species(FractNbr)%Init(iInit)%mySumOfMatchedParticles = 0 ParticleIndexNbr = 1 + ALLOCATE(AcceptedParts(0:chunkSize)) + AcceptedParts=-1 + AcceptedParts(0)=0 + DO i = 1,chunkSize + AcceptedParts(i) = SinglePointToElement(particle_positions(DimSend*(i-1)+1:DimSend*(i-1)+DimSend),doHALO=.FALSE.) + IF(AcceptedParts(i).NE.-1) AcceptedParts(0) = AcceptedParts(0) + 1 + END DO + Species(FractNbr)%Init(iInit)%mySumOfMatchedParticles = 0 + IF (Species(FractNbr)%Init(iInit)%ParticleEmissionType.EQ.0) CALL IncreaseMaxParticleNumber(AcceptedParts(0)) DO i = 1,chunkSize ! Find a free position in the PDM array - IF ((i.EQ.1).OR.PDM%ParticleInside(ParticleIndexNbr)) THEN - ParticleIndexNbr = PDM%nextFreePosition(Species(FractNbr)%Init(iInit)%mySumOfMatchedParticles + 1 + PDM%CurrentNextFreePosition) - END IF - IF (ParticleIndexNbr.NE.0) THEN + IF(AcceptedParts(i).NE.-1) THEN + Species(FractNbr)%Init(iInit)%mySumOfMatchedParticles = Species(FractNbr)%Init(iInit)%mySumOfMatchedParticles + 1 + ParticleIndexNbr = GetNextFreePosition(Species(FractNbr)%Init(iInit)%mySumOfMatchedParticles) PartState(1:DimSend,ParticleIndexNbr) = particle_positions(DimSend*(i-1)+1:DimSend*(i-1)+DimSend) - PDM%ParticleInside( ParticleIndexNbr) = .TRUE. - CALL LocateParticleInElement(ParticleIndexNbr,doHALO=.FALSE.) - IF (PDM%ParticleInside(ParticleIndexNbr)) THEN - Species(FractNbr)%Init(iInit)%mySumOfMatchedParticles = Species(FractNbr)%Init(iInit)%mySumOfMatchedParticles + 1 - PDM%IsNewPart(ParticleIndexNbr) = .TRUE. - PDM%dtFracPush(ParticleIndexNbr) = .FALSE. - END IF - ELSE - CALL ABORT(__STAMP__,'ERROR in SetParticlePosition:ParticleIndexNbr.EQ.0 - maximum nbr of particles reached?') + PDM%ParticleInside(ParticleIndexNbr)=.TRUE. + PEM%GlobalElemID(ParticleIndexNbr) = AcceptedParts(i) + IF(TrackingMethod.EQ.REFMAPPING) CALL GetPositionInRefElem(PartState(1:DimSend,ParticleIndexNbr),PartPosRef(1:3,ParticleIndexNbr),AcceptedParts(i)) + PDM%IsNewPart(ParticleIndexNbr) = .TRUE. + PDM%dtFracPush(ParticleIndexNbr) = .FALSE. END IF END DO + DEALLOCATE(AcceptedParts) #if USE_MPI END IF @@ -363,7 +364,7 @@ SUBROUTINE SetParticleVelocity(FractNbr,iInit,NbrOfParticle) USE MOD_part_emission_tools ,ONLY: CalcVelocity_gyrotroncircle USE MOD_Particle_Boundary_Vars ,ONLY: DoBoundaryParticleOutputHDF5 USE MOD_Particle_Boundary_Tools ,ONLY: StoreBoundaryParticleProperties -USE MOD_part_tools ,ONLY: BuildTransGaussNums, InRotRefFrameCheck +USE MOD_part_tools ,ONLY: BuildTransGaussNums, InRotRefFrameCheck, GetNextFreePosition USE MOD_Particle_Vars ,ONLY: CalcBulkElectronTemp,BulkElectronTemp #if USE_HDG USE MOD_HDG_Vars ,ONLY: UseFPC,FPC,UseEPC,EPC @@ -378,7 +379,7 @@ SUBROUTINE SetParticleVelocity(FractNbr,iInit,NbrOfParticle) INTEGER,INTENT(INOUT) :: NbrOfParticle !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: i, PositionNbr +INTEGER :: iPart, PositionNbr CHARACTER(30) :: velocityDistribution REAL :: VeloIC, VeloVecIC(3), maxwellfac, VeloVecNorm REAL :: iRanPart(3, NbrOfParticle), Vec3D(3),MPF @@ -401,15 +402,15 @@ SUBROUTINE SetParticleVelocity(FractNbr,iInit,NbrOfParticle) SELECT CASE(TRIM(velocityDistribution)) CASE('constant') - DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) + DO iPart = 1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr.GT.0) THEN PartState(4:6,PositionNbr) = VeloVecIC(1:3) * VeloIC END IF END DO CASE('gyrotron_circle') - DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) + DO iPart = 1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr.GT.0) THEN CALL CalcVelocity_gyrotroncircle(FractNbr, Vec3D, iInit, PositionNbr) PartState(4:6,PositionNbr) = Vec3D(1:3) @@ -419,8 +420,8 @@ SUBROUTINE SetParticleVelocity(FractNbr,iInit,NbrOfParticle) ! maxwell_lpn: Maxwell low particle number ! 2D_landmark: Ionization profile from T. Charoy, 2D axial-azimuthal particle-in-cell benchmark for low-temperature partially ! magnetized plasmas (2019) - DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) + DO iPart = 1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr.GT.0) THEN CALL CalcVelocity_maxwell_lpn(FractNbr, Vec3D, iInit=iInit) PartState(4:6,PositionNbr) = Vec3D(1:3) @@ -430,8 +431,8 @@ SUBROUTINE SetParticleVelocity(FractNbr,iInit,NbrOfParticle) IF(.NOT.CalcBulkElectronTemp) CALL abort(__STAMP__,& 'Velocity distribution 2D_Liu2010_neutralization requires CalcBulkElectronTemp=T') ! Use the global electron temperature if available - DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) + DO iPart = 1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr.GT.0) THEN CALL CalcVelocity_maxwell_lpn(FractNbr, Vec3D, Temperature=BulkElectronTemp*eV2Kelvin) PartState(4:6,PositionNbr) = Vec3D(1:3) @@ -440,8 +441,8 @@ SUBROUTINE SetParticleVelocity(FractNbr,iInit,NbrOfParticle) END IF END DO CASE('taylorgreenvortex') - DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) + DO iPart = 1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr.GT.0) THEN CALL CalcVelocity_taylorgreenvortex(FractNbr, Vec3D, iInit=iInit, Element=PEM%GlobalElemID(PositionNbr)) PartState(4:6,PositionNbr) = Vec3D(1:3) @@ -450,17 +451,17 @@ SUBROUTINE SetParticleVelocity(FractNbr,iInit,NbrOfParticle) CASE('maxwell') CALL BuildTransGaussNums(NbrOfParticle, iRanPart) maxwellfac = SQRT(BoltzmannConst*Species(FractNbr)%Init(iInit)%MWTemperatureIC/Species(FractNbr)%MassIC) - DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) + DO iPart = 1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr.GT.0) THEN - PartState(4:6,PositionNbr) = VeloIC *VeloVecIC(1:3) + iRanPart(1:3,i)*maxwellfac + PartState(4:6,PositionNbr) = VeloIC *VeloVecIC(1:3) + iRanPart(1:3,iPart)*maxwellfac END IF END DO CASE('IMD') ! read IMD particle velocity from *.chkpt file -> velocity space has already been read when particles position was done ! do nothing CASE('photon_SEE_energy') - DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) + DO iPart = 1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr .NE. 0) THEN CALL CalcVelocity_FromWorkFuncSEE(FractNbr, Vec3D, iInit=iInit) PartState(4:6,PositionNbr) = Vec3D(1:3) @@ -520,8 +521,8 @@ SUBROUTINE SetParticleVelocity(FractNbr,iInit,NbrOfParticle) END SELECT IF(UseRotRefFrame) THEN - DO i = 1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) + DO iPart = 1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (PositionNbr.GT.0) THEN PDM%InRotRefFrame(PositionNbr) = InRotRefFrameCheck(PositionNbr) ! Initialize velocity in the rotational frame of reference @@ -547,9 +548,9 @@ SUBROUTINE SetPartPosAndVeloEmissionDistribution(iSpec,iInit,NbrOfParticle) !USE MOD_Globals USE MOD_PreProc USE MOD_Globals ,ONLY: myrank,UNIT_StdOut,abort -USE MOD_part_tools ,ONLY: InitializeParticleMaxwell,InterpolateEmissionDistribution2D +USE MOD_part_tools ,ONLY: InitializeParticleMaxwell,InterpolateEmissionDistribution2D, GetNextFreePosition USE MOD_Mesh_Vars ,ONLY: nElems,offsetElem -USE MOD_Particle_Vars ,ONLY: Species, PDM, PartState, PEM, LastPartPos +USE MOD_Particle_Vars ,ONLY: Species, PDM, PartState, PEM, LastPartPos, PartPosRef USE MOD_Particle_Tracking ,ONLY: ParticleInsideCheck USE MOD_Mesh_Tools ,ONLY: GetCNElemID USE MOD_Particle_Emission_Vars ,ONLY: EmissionDistributionDim, EmissionDistributionN @@ -559,7 +560,7 @@ SUBROUTINE SetPartPosAndVeloEmissionDistribution(iSpec,iInit,NbrOfParticle) USE MOD_Equation ,ONLY: ExactFunc USE MOD_Mesh_Vars ,ONLY: Elem_xGP,sJ USE MOD_Interpolation_Vars ,ONLY: NodeTypeVISU,NodeType -USE MOD_Eval_xyz ,ONLY: TensorProductInterpolation +USE MOD_Eval_xyz ,ONLY: TensorProductInterpolation, GetPositionInRefElem USE MOD_Mesh_Vars ,ONLY: NGeo,XCL_NGeo,XiCL_NGeo,wBaryCL_NGeo,offsetElem USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared,BoundsOfElem_Shared USE MOD_Mesh_Tools ,ONLY: GetCNElemID @@ -738,12 +739,12 @@ SUBROUTINE SetPartPosAndVeloEmissionDistribution(iSpec,iInit,NbrOfParticle) ! Exclude particles outside of the element IF (InsideFlag) THEN NbrOfParticle = NbrOfParticle + 1 - PositionNbr = PDM%nextFreePosition(NbrOfParticle+PDM%CurrentNextFreePosition) + PositionNbr = GetNextFreePosition(NbrOfParticle) PEM%GlobalElemID(PositionNbr) = GlobalElemID PDM%ParticleInside(PositionNbr) = .TRUE. - IF((PositionNbr.GE.PDM%maxParticleNumber).OR.& - (PositionNbr.EQ.0)) CALL abort(__STAMP__,'Emission: Increase maxParticleNumber!',PositionNbr) PartState(1:3,PositionNbr) = RandomPos(1:3) + IF(TrackingMethod.EQ.REFMAPPING) & + CALL GetPositionInRefElem(PartState(1:3,PositionNbr),PartPosRef(1:3,PositionNbr),GlobalElemID) CALL InitializeParticleMaxwell(PositionNbr,iSpec,iElem,Mode=2,iInit=iInit) END IF END DO ! nPart @@ -775,7 +776,7 @@ SUBROUTINE SetPartPosAndVeloEmissionDistribution(iSpec,iInit,NbrOfParticle) SELECT CASE(TrackingMethod) CASE(REFMAPPING,TRACING) ! Attention: NbrOfParticle+PDM%CurrentNextFreePosition + 1 - PositionNbr = PDM%nextFreePosition(NbrOfParticle+PDM%CurrentNextFreePosition + 1) + PositionNbr = GetNextFreePosition(NbrOfParticle+1) PEM%GlobalElemID(PositionNbr) = GlobalElemID LastPartPos(1:3,PositionNbr) = RandomPos(1:3) InsideFlag = ParticleInsideCheck(RandomPos,PositionNbr,GlobalElemID) @@ -788,11 +789,9 @@ SUBROUTINE SetPartPosAndVeloEmissionDistribution(iSpec,iInit,NbrOfParticle) ! Exclude particles outside of the element IF (InsideFlag) THEN NbrOfParticle = NbrOfParticle + 1 - PositionNbr = PDM%nextFreePosition(NbrOfParticle+PDM%CurrentNextFreePosition) + PositionNbr = GetNextFreePosition(NbrOfParticle) PEM%GlobalElemID(PositionNbr) = GlobalElemID PDM%ParticleInside(PositionNbr) = .TRUE. - IF((PositionNbr.GE.PDM%maxParticleNumber).OR.& - (PositionNbr.EQ.0)) CALL abort(__STAMP__,'Emission: Increase maxParticleNumber!',PositionNbr) PartState(1:3,PositionNbr) = RandomPos(1:3) CALL InitializeParticleMaxwell(PositionNbr,iSpec,iElem,Mode=2,iInit=iInit) END IF diff --git a/src/particles/emission/particle_surface_flux.f90 b/src/particles/emission/particle_surface_flux.f90 index a1aa239e7..08fbb12f1 100644 --- a/src/particles/emission/particle_surface_flux.f90 +++ b/src/particles/emission/particle_surface_flux.f90 @@ -35,11 +35,11 @@ SUBROUTINE ParticleSurfaceflux() ! Modules USE MOD_Globals USE MOD_Particle_Vars -USE MOD_part_tools ,ONLY: CalcRadWeightMPF +USE MOD_part_tools ,ONLY: CalcRadWeightMPF, IncreaseMaxParticleNumber USE MOD_DSMC_Vars ,ONLY: useDSMC, CollisMode, RadialWeighting, DSMC USE MOD_Eval_xyz ,ONLY: GetPositionInRefElem USE MOD_Mesh_Vars ,ONLY: SideToElem, offsetElem -USE MOD_Part_Tools ,ONLY: GetParticleWeight +USE MOD_Part_Tools ,ONLY: GetParticleWeight, GetNextFreePosition USE MOD_Part_Emission_Tools ,ONLY: SetParticleChargeAndMass, SetParticleMPF USE MOD_Particle_Analyze_Vars ,ONLY: CalcPartBalance, CalcSurfFluxInfo, nPartIn, PartEkinIn USE MOD_Particle_Analyze_Tools ,ONLY: CalcEkinPart @@ -218,39 +218,34 @@ SUBROUTINE ParticleSurfaceflux() ParticleIndexNbr = 1 DO iPart=1,PartInsSubSide IF ((iPart.EQ.1).OR.PDM%ParticleInside(ParticleIndexNbr)) & - ParticleIndexNbr = PDM%nextFreePosition(iPartTotal + 1 + PDM%CurrentNextFreePosition) - IF (ParticleIndexNbr .ne. 0) THEN - PartState(1:3,ParticleIndexNbr) = particle_positions(3*(iPart-1)+1:3*(iPart-1)+3) - IF (SF%VeloIsNormal.AND.(.NOT.TriaSurfaceFlux)) THEN - PartState(4:5,ParticleIndexNbr) = particle_xis(2*(iPart-1)+1:2*(iPart-1)+2) !use velo as dummy-storage for xi! - END IF - LastPartPos(1:3,ParticleIndexNbr)=PartState(1:3,ParticleIndexNbr) + ParticleIndexNbr = GetNextFreePosition(iPartTotal+1) + PartState(1:3,ParticleIndexNbr) = particle_positions(3*(iPart-1)+1:3*(iPart-1)+3) + IF (SF%VeloIsNormal.AND.(.NOT.TriaSurfaceFlux)) THEN + PartState(4:5,ParticleIndexNbr) = particle_xis(2*(iPart-1)+1:2*(iPart-1)+2) !use velo as dummy-storage for xi! + END IF + LastPartPos(1:3,ParticleIndexNbr)=PartState(1:3,ParticleIndexNbr) #if defined(IMPA) || defined(ROS) - IF(TrackingMethod.EQ.REFMAPPING) CALL GetPositionInRefElem(PartState(1:3,ParticleIndexNbr),PartPosRef(1:3,ParticleIndexNbr),globElemId) + IF(TrackingMethod.EQ.REFMAPPING) CALL GetPositionInRefElem(PartState(1:3,ParticleIndexNbr),PartPosRef(1:3,ParticleIndexNbr),globElemId) #endif /*IMPA*/ - PDM%ParticleInside(ParticleIndexNbr) = .TRUE. - PDM%dtFracPush(ParticleIndexNbr) = .TRUE. - PDM%IsNewPart(ParticleIndexNbr) = .TRUE. - PEM%GlobalElemID(ParticleIndexNbr) = globElemId - PEM%LastGlobalElemID(ParticleIndexNbr) = globElemId !needed when ParticlePush is not executed, e.g. "delay" - iPartTotal = iPartTotal + 1 - IF (UseVarTimeStep) THEN - PartTimeStep(ParticleIndexNbr) = GetParticleTimeStep(PartState(1,ParticleIndexNbr),PartState(2,ParticleIndexNbr), & - PEM%LocalElemID(ParticleIndexNbr)) - END IF - IF (RadialWeighting%DoRadialWeighting) THEN - PartMPF(ParticleIndexNbr) = CalcRadWeightMPF(PartState(2,ParticleIndexNbr), iSpec,ParticleIndexNbr) - END IF - IF(CalcSurfFluxInfo) THEN - SF%SampledMassflow = SF%SampledMassflow + GetParticleWeight(ParticleIndexNbr) - END IF + PDM%ParticleInside(ParticleIndexNbr) = .TRUE. + PDM%dtFracPush(ParticleIndexNbr) = .TRUE. + PDM%IsNewPart(ParticleIndexNbr) = .TRUE. + PEM%GlobalElemID(ParticleIndexNbr) = globElemId + PEM%LastGlobalElemID(ParticleIndexNbr) = globElemId !needed when ParticlePush is not executed, e.g. "delay" + iPartTotal = iPartTotal + 1 + IF (UseVarTimeStep) THEN + PartTimeStep(ParticleIndexNbr) = GetParticleTimeStep(PartState(1,ParticleIndexNbr),PartState(2,ParticleIndexNbr), & + PEM%LocalElemID(ParticleIndexNbr)) + END IF + IF (RadialWeighting%DoRadialWeighting) THEN + PartMPF(ParticleIndexNbr) = CalcRadWeightMPF(PartState(2,ParticleIndexNbr), iSpec,ParticleIndexNbr) + END IF + IF(CalcSurfFluxInfo) THEN + SF%SampledMassflow = SF%SampledMassflow + GetParticleWeight(ParticleIndexNbr) + END IF #ifdef CODE_ANALYZE - CALL AnalyzePartPos(ParticleIndexNbr) + CALL AnalyzePartPos(ParticleIndexNbr) #endif /*CODE_ANALYZE*/ - ELSE - IPWRITE(*,*) 'Total number of particles emitted: ', PartsEmitted, ' Current number of particles: ', PartInsSubSide - CALL abort(__STAMP__,'ERROR in ParticleSurfaceflux: ParticleIndexNbr.EQ.0 - maximum nbr of particles reached?') - END IF END DO !----- 2a.: set velocities if special for each subside CALL SetSurfacefluxVelocities(iSpec,iSF,iSample,jSample,iSide,BCSideID,SideID,ElemID,NbrOfParticle,PartInsSubSide) @@ -259,12 +254,6 @@ SUBROUTINE ParticleSurfaceflux() IF (useDSMC) THEN IF (DSMC%DoAmbipolarDiff) CALL AD_SetSFElectronVelo(iSpec,iSF,iSample,jSample,iSide,BCSideID,SideID,ElemID,NbrOfParticle,PartInsSubSide,particle_xis) - DO iPart = 1, NbrOfParticle - PositionNbr = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) - IF (PositionNbr .EQ. 0) THEN - CALL abort(__STAMP__,'ERROR in InitialParticleInserting: No free particle index - maximum nbr of particles reached?') - END IF - END DO END IF IF (SF%VeloIsNormal .AND. .NOT.TriaSurfaceFlux) DEALLOCATE(particle_xis) @@ -293,14 +282,24 @@ SUBROUTINE ParticleSurfaceflux() ! Compute number of input particles and energy nPartIn(iSpec)=nPartIn(iSpec) + NbrOfParticle DO iPart=1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) - IF (PositionNbr.NE.0) PartEkinIn(iSpec) = PartEkinIn(iSpec)+CalcEkinPart(PositionNbr) + PositionNbr = GetNextFreePosition(iPart) + PartEkinIn(iSpec) = PartEkinIn(iSpec)+CalcEkinPart(PositionNbr) END DO ! iPart END IF ! CalcPartBalance ! instead of an UpdateNextfreePosition we update the particleVecLength only - enough ?!? - PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + NbrOfParticle - PDM%ParticleVecLength = PDM%ParticleVecLength + NbrOfParticle - IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL abort(__STAMP__, 'Error ParticleSurfaceflux: Maximum number of particles reached!') + IF(iPartTotal.GT.0) THEN + PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + NbrOfParticle + PDM%ParticleVecLength = MAX(PDM%ParticleVecLength,GetNextFreePosition(0)) + END IF +#ifdef CODE_ANALYZE + IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL Abort(__STAMP__,'PDM%ParticleVeclength exceeds PDM%maxParticleNumber, Difference:',IntInfoOpt=PDM%ParticleVeclength-PDM%maxParticleNumber) + DO iPart=PDM%ParticleVecLength+1,PDM%maxParticleNumber + IF (PDM%ParticleInside(iPart)) THEN + IPWRITE(*,*) iPart,PDM%ParticleVecLength,PDM%maxParticleNumber + CALL Abort(__STAMP__,'Particle outside PDM%ParticleVeclength',IntInfoOpt=iPart) + END IF + END DO +#endif #if USE_LOADBALANCE CALL LBPauseTime(LB_SURFFLUX,tLBStart) #endif /*USE_LOADBALANCE*/ @@ -526,9 +525,9 @@ SUBROUTINE SetInnerEnergies(iSpec, iSF, NbrOfParticle) ! MODULES USE MOD_Globals USE MOD_DSMC_Vars ,ONLY: SpecDSMC -USE MOD_Particle_Vars ,ONLY: PDM USE MOD_DSMC_PolyAtomicModel ,ONLY: DSMC_SetInternalEnr_Poly USE MOD_part_emission_tools ,ONLY: DSMC_SetInternalEnr_LauxVFD +USE MOD_Part_Tools ,ONLY: GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -541,16 +540,13 @@ SUBROUTINE SetInnerEnergies(iSpec, iSF, NbrOfParticle) INTEGER :: iPart, PositionNbr !=================================================================================================================================== iPart = 1 -DO WHILE (iPart .le. NbrOfParticle) - PositionNbr = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) - IF (PositionNbr .ne. 0) THEN +DO iPart=1,NbrOfParticle + PositionNbr = GetNextFreePosition(iPart) IF (SpecDSMC(iSpec)%PolyatomicMol) THEN CALL DSMC_SetInternalEnr_Poly(iSpec,iSF,PositionNbr,2) ELSE CALL DSMC_SetInternalEnr_LauxVFD(iSpec, iSF, PositionNbr,2) END IF - END IF - iPart = iPart + 1 END DO END SUBROUTINE SetInnerEnergies @@ -1184,7 +1180,7 @@ SUBROUTINE SetSurfacefluxVelocities(iSpec,iSF,iSample,jSample,iSide,BCSideID,Sid USE MOD_Particle_Surfaces_Vars ,ONLY: SurfMeshSubSideData, TriaSurfaceFlux USE MOD_Particle_Surfaces ,ONLY: CalcNormAndTangBezier USE MOD_Particle_Sampling_Vars ,ONLY: AdaptBCMapElemToSample, AdaptBCMacroVal -USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck +USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck, GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -1265,25 +1261,23 @@ SUBROUTINE SetSurfacefluxVelocities(iSpec,iSF,iSample,jSample,iSide,BCSideID,Sid VeloVecIC(1:3) = VeloVecIC(1:3) / VECNORM(VeloVecIC(1:3)) END IF DO i = NbrOfParticle-PartIns+1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr .NE. 0) THEN - ! In case of side-normal velocities: calc n-vector at particle position, xi was saved in PartState(4:5) - IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. TriaSurfaceFlux) THEN - vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) - vec_t1(1:3) = 0. !dummy - vec_t2(1:3) = 0. !dummy - ELSE IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN - CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),xi=PartState(4,PositionNbr),eta=PartState(5,PositionNbr),SideID=SideID ) - vec_nIn(1:3) = -vec_nIn(1:3) - vec_t1(1:3) = 0. !dummy - vec_t2(1:3) = 0. !dummy - ELSE - vec_nIn(1:3) = VeloVecIC(1:3) - END IF !VeloIsNormal - ! Build complete velo-vector - Vec3D(1:3) = vec_nIn(1:3) * Species(iSpec)%Surfaceflux(iSF)%VeloIC - PartState(4:6,PositionNbr) = Vec3D(1:3) - END IF !PositionNbr .NE. 0 + PositionNbr = GetNextFreePosition(i) + ! In case of side-normal velocities: calc n-vector at particle position, xi was saved in PartState(4:5) + IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. TriaSurfaceFlux) THEN + vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) + vec_t1(1:3) = 0. !dummy + vec_t2(1:3) = 0. !dummy + ELSE IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN + CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),xi=PartState(4,PositionNbr),eta=PartState(5,PositionNbr),SideID=SideID ) + vec_nIn(1:3) = -vec_nIn(1:3) + vec_t1(1:3) = 0. !dummy + vec_t2(1:3) = 0. !dummy + ELSE + vec_nIn(1:3) = VeloVecIC(1:3) + END IF !VeloIsNormal + ! Build complete velo-vector + Vec3D(1:3) = vec_nIn(1:3) * Species(iSpec)%Surfaceflux(iSF)%VeloIC + PartState(4:6,PositionNbr) = Vec3D(1:3) END DO !i = ...NbrOfParticle CASE('maxwell','maxwell_lpn') !-- determine envelope for most efficient ARM [Garcia and Wagner 2006, JCP217-2] @@ -1307,140 +1301,136 @@ SUBROUTINE SetSurfacefluxVelocities(iSpec,iSF,iSample,jSample,iSide,BCSideID,Sid END IF !low speed / high speed / rayleigh flow DO i = NbrOfParticle-PartIns+1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr .NE. 0) THEN - !-- 0a.: In case of side-normal velocities: calc n-/t-vectors at particle position, xi was saved in PartState(4:5) - IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. TriaSurfaceFlux) THEN - vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) - vec_t1(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_t1(1:3) - vec_t2(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_t2(1:3) - ELSE IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN - CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),tang1=vec_t1(1:3),tang2=vec_t2(1:3) & - ,xi=PartState(4,PositionNbr),eta=PartState(5,PositionNbr),SideID=SideID ) - vec_nIn(1:3) = -vec_nIn(1:3) - END IF !VeloIsNormal - !-- 1.: determine zstar (initial generation of potentially too many RVu is for needed indentities of RVu used multiple times! - SELECT CASE(envelope) - CASE(0) - CALL RANDOM_NUMBER(RandVal1) - zstar = -SQRT(-LOG(RandVal1)) - CASE(1) - DO - CALL RANDOM_NUMBER(RandVal2) - zstar = -SQRT(a*a-LOG(RandVal2(1))) - IF ( -(a-zstar)/zstar .GT. RandVal2(2)) THEN + PositionNbr = GetNextFreePosition(i) + !-- 0a.: In case of side-normal velocities: calc n-/t-vectors at particle position, xi was saved in PartState(4:5) + IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. TriaSurfaceFlux) THEN + vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) + vec_t1(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_t1(1:3) + vec_t2(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_t2(1:3) + ELSE IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN + CALL CalcNormAndTangBezier( nVec=vec_nIn(1:3),tang1=vec_t1(1:3),tang2=vec_t2(1:3) & + ,xi=PartState(4,PositionNbr),eta=PartState(5,PositionNbr),SideID=SideID ) + vec_nIn(1:3) = -vec_nIn(1:3) + END IF !VeloIsNormal + !-- 1.: determine zstar (initial generation of potentially too many RVu is for needed indentities of RVu used multiple times! + SELECT CASE(envelope) + CASE(0) + CALL RANDOM_NUMBER(RandVal1) + zstar = -SQRT(-LOG(RandVal1)) + CASE(1) + DO + CALL RANDOM_NUMBER(RandVal2) + zstar = -SQRT(a*a-LOG(RandVal2(1))) + IF ( -(a-zstar)/zstar .GT. RandVal2(2)) THEN + EXIT + END IF + END DO + CASE(2) + z = 0.5*(a-SQRT(a*a+2.)) + beta = a-(1.0-a)*(a-z) + DO + CALL RANDOM_NUMBER(RandVal3) + IF (EXP(-(beta*beta))/(EXP(-(beta*beta))+2.0*(a-z)*(a-beta)*EXP(-(z*z))).GT.RandVal3(1)) THEN + zstar=-SQRT(beta*beta-LOG(RandVal3(2))) + IF ( -(a-zstar)/zstar .GT. RandVal3(3)) THEN EXIT END IF - END DO - CASE(2) - z = 0.5*(a-SQRT(a*a+2.)) - beta = a-(1.0-a)*(a-z) - DO - CALL RANDOM_NUMBER(RandVal3) - IF (EXP(-(beta*beta))/(EXP(-(beta*beta))+2.0*(a-z)*(a-beta)*EXP(-(z*z))).GT.RandVal3(1)) THEN - zstar=-SQRT(beta*beta-LOG(RandVal3(2))) - IF ( -(a-zstar)/zstar .GT. RandVal3(3)) THEN - EXIT - END IF - ELSE - zstar=beta+(a-beta)*RandVal3(2) - IF ( (a-zstar)/(a-z)*EXP(z*z-(zstar*zstar)) .GT. RandVal3(3)) THEN - EXIT - END IF + ELSE + zstar=beta+(a-beta)*RandVal3(2) + IF ( (a-zstar)/(a-z)*EXP(z*z-(zstar*zstar)) .GT. RandVal3(3)) THEN + EXIT END IF - END DO - CASE(3) - DO - CALL RANDOM_NUMBER(RandVal3) - u = RandVal3(1) - IF ( a*SQRT(PI)/(a*SQRT(PI)+1+a*a) .GT. u) THEN + END IF + END DO + CASE(3) + DO + CALL RANDOM_NUMBER(RandVal3) + u = RandVal3(1) + IF ( a*SQRT(PI)/(a*SQRT(PI)+1+a*a) .GT. u) THEN ! IF (.NOT.DoZigguratSampling) THEN !polar method - IF (RandN_in_Mem) THEN !reusing second RandN form previous polar method - RandN = RandN_save - RandN_in_Mem=.FALSE. - ELSE - Velosq = 2 - DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) - CALL RANDOM_NUMBER(RandVal2) - Velo1 = 2.*RandVal2(1) - 1. - Velo2 = 2.*RandVal2(2) - 1. - Velosq = Velo1**2 + Velo2**2 - END DO - RandN = Velo1*SQRT(-2*LOG(Velosq)/Velosq) - RandN_save = Velo2*SQRT(-2*LOG(Velosq)/Velosq) - RandN_in_Mem=.TRUE. - END IF + IF (RandN_in_Mem) THEN !reusing second RandN form previous polar method + RandN = RandN_save + RandN_in_Mem=.FALSE. + ELSE + Velosq = 2 + DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) + CALL RANDOM_NUMBER(RandVal2) + Velo1 = 2.*RandVal2(1) - 1. + Velo2 = 2.*RandVal2(2) - 1. + Velosq = Velo1**2 + Velo2**2 + END DO + RandN = Velo1*SQRT(-2*LOG(Velosq)/Velosq) + RandN_save = Velo2*SQRT(-2*LOG(Velosq)/Velosq) + RandN_in_Mem=.TRUE. + END IF ! ELSE !ziggurat method ! RandN=rnor() ! END IF - zstar = -1./SQRT(2.)*ABS(RandN) - EXIT - ELSE IF ( (a*SQRT(PI)+1.)/(a*SQRT(PI)+1+a*a) .GT. u) THEN - zstar = -SQRT(-LOG(RandVal3(2))) + zstar = -1./SQRT(2.)*ABS(RandN) + EXIT + ELSE IF ( (a*SQRT(PI)+1.)/(a*SQRT(PI)+1+a*a) .GT. u) THEN + zstar = -SQRT(-LOG(RandVal3(2))) + EXIT + ELSE + zstar = (1.0-SQRT(RandVal3(2)))*a + IF (EXP(-(zstar*zstar)).GT.RandVal3(3)) THEN EXIT - ELSE - zstar = (1.0-SQRT(RandVal3(2)))*a - IF (EXP(-(zstar*zstar)).GT.RandVal3(3)) THEN - EXIT - END IF END IF - END DO - CASE(4) - DO - CALL RANDOM_NUMBER(RandVal3) - IF (1.0/(2.0*a*SQRT(PI)+1.0).GT.RandVal3(1)) THEN - zstar=-SQRT(-LOG(RandVal3(2))) - ELSE + END IF + END DO + CASE(4) + DO + CALL RANDOM_NUMBER(RandVal3) + IF (1.0/(2.0*a*SQRT(PI)+1.0).GT.RandVal3(1)) THEN + zstar=-SQRT(-LOG(RandVal3(2))) + ELSE ! IF (.NOT.DoZigguratSampling) THEN !polar method - IF (RandN_in_Mem) THEN !reusing second RandN form previous polar method - RandN = RandN_save - RandN_in_Mem=.FALSE. - ELSE - Velosq = 2 - DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) - CALL RANDOM_NUMBER(RandVal2) - Velo1 = 2.*RandVal2(1) - 1. - Velo2 = 2.*RandVal2(2) - 1. - Velosq = Velo1**2 + Velo2**2 - END DO - RandN = Velo1*SQRT(-2*LOG(Velosq)/Velosq) - RandN_save = Velo2*SQRT(-2*LOG(Velosq)/Velosq) - RandN_in_Mem=.TRUE. - END IF + IF (RandN_in_Mem) THEN !reusing second RandN form previous polar method + RandN = RandN_save + RandN_in_Mem=.FALSE. + ELSE + Velosq = 2 + DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) + CALL RANDOM_NUMBER(RandVal2) + Velo1 = 2.*RandVal2(1) - 1. + Velo2 = 2.*RandVal2(2) - 1. + Velosq = Velo1**2 + Velo2**2 + END DO + RandN = Velo1*SQRT(-2*LOG(Velosq)/Velosq) + RandN_save = Velo2*SQRT(-2*LOG(Velosq)/Velosq) + RandN_in_Mem=.TRUE. + END IF ! ELSE !ziggurat method ! RandN=rnor() ! END IF - zstar = 1./SQRT(2.)*RandN - END IF - IF ( (a-zstar)/a .GT. RandVal3(3)) THEN - EXIT - END IF - END DO - CASE DEFAULT - CALL abort(__STAMP__,'ERROR in SurfaceFlux: Wrong envelope in SetSurfacefluxVelocities!') - END SELECT - !-- 2.: sample normal directions and build complete velo-vector - Vec3D(1:3) = vec_nIn(1:3) * SQRT(2.*BoltzmannConst*T/Species(iSpec)%MassIC)*(a-zstar) + zstar = 1./SQRT(2.)*RandN + END IF + IF ( (a-zstar)/a .GT. RandVal3(3)) THEN + EXIT + END IF + END DO + CASE DEFAULT + CALL abort(__STAMP__,'ERROR in SurfaceFlux: Wrong envelope in SetSurfacefluxVelocities!') + END SELECT + !-- 2.: sample normal directions and build complete velo-vector + Vec3D(1:3) = vec_nIn(1:3) * SQRT(2.*BoltzmannConst*T/Species(iSpec)%MassIC)*(a-zstar) ! IF (.NOT.DoZigguratSampling) THEN !polar method - Velosq = 2 - DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) - CALL RANDOM_NUMBER(RandVal2) - Velo1 = 2.*RandVal2(1) - 1. - Velo2 = 2.*RandVal2(2) - 1. - Velosq = Velo1**2 + Velo2**2 - END DO - Velo1 = Velo1*SQRT(-2*LOG(Velosq)/Velosq) - Velo2 = Velo2*SQRT(-2*LOG(Velosq)/Velosq) + Velosq = 2 + DO WHILE ((Velosq .GE. 1.) .OR. (Velosq .EQ. 0.)) + CALL RANDOM_NUMBER(RandVal2) + Velo1 = 2.*RandVal2(1) - 1. + Velo2 = 2.*RandVal2(2) - 1. + Velosq = Velo1**2 + Velo2**2 + END DO + Velo1 = Velo1*SQRT(-2*LOG(Velosq)/Velosq) + Velo2 = Velo2*SQRT(-2*LOG(Velosq)/Velosq) ! ELSE !ziggurat method ! Velo1=rnor() ! Velo2=rnor() ! END IF - Vec3D(1:3) = Vec3D(1:3) + vec_t1(1:3) * ( Velo_t1+Velo1*SQRT(BoltzmannConst*T/Species(iSpec)%MassIC) ) - Vec3D(1:3) = Vec3D(1:3) + vec_t2(1:3) * ( Velo_t2+Velo2*SQRT(BoltzmannConst*T/Species(iSpec)%MassIC) ) - PartState(4:6,PositionNbr) = Vec3D(1:3) - ELSE !PositionNbr .EQ. 0 - CALL abort(__STAMP__,'PositionNbr .EQ. 0!') - END IF !PositionNbr .NE. 0 + Vec3D(1:3) = Vec3D(1:3) + vec_t1(1:3) * ( Velo_t1+Velo1*SQRT(BoltzmannConst*T/Species(iSpec)%MassIC) ) + Vec3D(1:3) = Vec3D(1:3) + vec_t2(1:3) * ( Velo_t2+Velo2*SQRT(BoltzmannConst*T/Species(iSpec)%MassIC) ) + PartState(4:6,PositionNbr) = Vec3D(1:3) END DO !i = ...NbrOfParticle CASE DEFAULT CALL abort(__STAMP__,'ERROR in SurfaceFlux: Wrong velocity distribution!') @@ -1448,14 +1438,12 @@ SUBROUTINE SetSurfacefluxVelocities(iSpec,iSF,iSample,jSample,iSide,BCSideID,Sid IF(UseRotRefFrame) THEN DO i = NbrOfParticle-PartIns+1,NbrOfParticle - PositionNbr = PDM%nextFreePosition(i+PDM%CurrentNextFreePosition) - IF (PositionNbr.GT.0) THEN - ! Detect if particle is within a RotRefDomain - PDM%InRotRefFrame(PositionNbr) = InRotRefFrameCheck(PositionNbr) - ! Initialize velocity in the rotational frame of reference - IF(PDM%InRotRefFrame(PositionNbr)) THEN - PartVeloRotRef(1:3,PositionNbr) = PartState(4:6,PositionNbr) - CROSS(RotRefFrameOmega(1:3),PartState(1:3,PositionNbr)) - END IF + PositionNbr = GetNextFreePosition(i) + ! Detect if particle is within a RotRefDomain + PDM%InRotRefFrame(PositionNbr) = InRotRefFrameCheck(PositionNbr) + ! Initialize velocity in the rotational frame of reference + IF(PDM%InRotRefFrame(PositionNbr)) THEN + PartVeloRotRef(1:3,PositionNbr) = PartState(4:6,PositionNbr) - CROSS(RotRefFrameOmega(1:3),PartState(1:3,PositionNbr)) END IF END DO END IF diff --git a/src/particles/mcc/mcc.f90 b/src/particles/mcc/mcc.f90 index 71197edb5..47193d11c 100644 --- a/src/particles/mcc/mcc.f90 +++ b/src/particles/mcc/mcc.f90 @@ -63,7 +63,7 @@ SUBROUTINE MonteCarloCollision(iElem) ! ROUTINES USE MOD_DSMC_Analyze ,ONLY: CalcMeanFreePath USE MOD_DSMC_BGGas ,ONLY: BGGas_AssignParticleProperties -USE MOD_part_tools ,ONLY: GetParticleWeight, CalcVelocity_maxwell_particle +USE MOD_part_tools ,ONLY: GetParticleWeight, CalcVelocity_maxwell_particle, GetNextFreePosition USE MOD_Part_Emission_Tools ,ONLY: CalcVelocity_maxwell_lpn USE MOD_DSMC_Collis ,ONLY: DSMC_perform_collision USE MOD_Mesh_Tools ,ONLY: GetCNElemID @@ -274,7 +274,7 @@ SUBROUTINE MonteCarloCollision(iElem) ! Clone the regular particle (re-using the index of the previous particle if it didn't collide) IF(PartIndex.EQ.0) THEN DSMCSumOfFormedParticles = DSMCSumOfFormedParticles + 1 - PartIndex = PDM%nextFreePosition(DSMCSumOfFormedParticles+PDM%CurrentNextFreePosition) + PartIndex = GetNextFreePosition() IF (PartIndex.EQ.0) THEN CALL Abort(__STAMP__,'ERROR in MCC: MaxParticleNumber should be increased!') END IF @@ -410,7 +410,7 @@ SUBROUTINE MonteCarloCollision(iElem) IF(ChemReac%CollCaseInfo(iCase)%HasXSecReaction) THEN IF(bggPartIndex.EQ.0) THEN DSMCSumOfFormedParticles = DSMCSumOfFormedParticles + 1 - bggPartIndex = PDM%nextFreePosition(DSMCSumOfFormedParticles+PDM%CurrentNextFreePosition) + bggPartIndex = GetNextFreePosition() IF (bggPartIndex.EQ.0) THEN CALL Abort(__STAMP__,'ERROR in MCC: MaxParticleNumber should be increased!') END IF @@ -466,7 +466,7 @@ SUBROUTINE MonteCarloCollision(iElem) ! Creating a new background gas particle IF(bggPartIndex.EQ.0) THEN DSMCSumOfFormedParticles = DSMCSumOfFormedParticles + 1 - bggPartIndex = PDM%nextFreePosition(DSMCSumOfFormedParticles+PDM%CurrentNextFreePosition) + bggPartIndex = GetNextFreePosition() IF (bggPartIndex.EQ.0) THEN CALL Abort(__STAMP__,'ERROR in MCC: MaxParticleNumber should be increased!') END IF diff --git a/src/particles/particle_init.f90 b/src/particles/particle_init.f90 index 338c0e887..b475a6239 100644 --- a/src/particles/particle_init.f90 +++ b/src/particles/particle_init.f90 @@ -91,8 +91,16 @@ SUBROUTINE DefineParametersParticles() CALL prms%CreateRealOption( 'InitialIonizationChargeAverage' , 'Average charge for each atom/molecule in the cell '//& '(corresponds to the ionization degree)') -CALL prms%CreateIntOption( 'Part-MaxParticleNumber', 'Maximum number of Particles per proc (used for array init)'& - , '1') +CALL prms%CreateIntOption( 'Part-MaxParticleNumber', 'Maximum allowed particles per core, 0=no limit') +CALL prms%CreateRealOption( 'Part-MaxPartNumIncrease', 'How much shall the PDM%MaxParticleNumber be incresed if it is full'& + , '0.1') +CALL prms%CreateLogicalOption( 'Part-RearrangePartIDs', 'Rearrange PartIDs in the process of reducing maxPartNum to allow lower memory usage'& + , '.TRUE.') +#if USE_MPI +CALL prms%CreateLogicalOption( 'Part-MPI-UNFP-afterPartSend', 'UpdateNextFreePosition after MPIParticleSend to reduce '//& + 'PDM%maxParticleNummber increase and decreace operations'& + , '.FALSE.') +#endif CALL prms%CreateIntOption( 'Part-NumberOfRandomSeeds' , 'Number of Seeds for Random Number Generator'//& 'Choose nRandomSeeds \n'//& '=-1 Random \n'//& @@ -427,9 +435,18 @@ SUBROUTINE InitializeVariables() !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES REAL :: ManualTimeStepParticle ! temporary variable +CHARACTER(32) :: hilf !=================================================================================================================================== ! Read basic particle parameter -PDM%maxParticleNumber = GETINT('Part-maxParticleNumber','1') +WRITE(UNIT=hilf,FMT='(I0)') HUGE(PDM%maxAllowedParticleNumber) +PDM%maxAllowedParticleNumber = GETINT('Part-maxParticleNumber',TRIM(hilf)) +PDM%MaxPartNumIncrease = GETREAL('Part-MaxPartNumIncrease','0.1') +PDM%RearrangePartIDs = GETLOGICAL('Part-RearrangePartIDs','.TRUE.') +#if USE_MPI +PDM%UNFPafterMPIPartSend = GETLOGICAL('Part-MPI-UNFP-afterPartSend','.FALSE.') +#endif +PDM%maxParticleNumber=1 +PDM%ParticleVecLength=0 CALL AllocateParticleArrays() CALL InitializeVariablesRandomNumbers() @@ -1174,6 +1191,7 @@ SUBROUTINE InitialIonization() USE MOD_part_emission_tools ,ONLY: CalcVelocity_maxwell_lpn USE MOD_DSMC_Vars ,ONLY: useDSMC USE MOD_Eval_xyz ,ONLY: TensorProductInterpolation +USE MOD_Part_Tools ,ONLY: GetNextFreePosition #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_LOADBALANCE*/ @@ -1282,9 +1300,7 @@ SUBROUTINE InitialIonization() DO iPart=1,ElemCharge(iElem) ! 1 electron for each charge of each element ! Set the next free position in the particle vector list - PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + 1 - ParticleIndexNbr = PDM%nextFreePosition(PDM%CurrentNextFreePosition) - PDM%ParticleVecLength = PDM%ParticleVecLength + 1 + ParticleIndexNbr = GetNextFreePosition() !Set new SpeciesID of new particle (electron) PDM%ParticleInside(ParticleIndexNbr) = .true. diff --git a/src/particles/particle_mpi/particle_mpi.f90 b/src/particles/particle_mpi/particle_mpi.f90 index 3ca2c5959..0c854afe7 100644 --- a/src/particles/particle_mpi/particle_mpi.f90 +++ b/src/particles/particle_mpi/particle_mpi.f90 @@ -431,6 +431,7 @@ SUBROUTINE MPIParticleSend(UseOldVecLength) USE MOD_Particle_Tracking_Vars, ONLY:TrackingMethod USE MOD_Particle_Vars, ONLY:PartState,PartSpecies,usevMPF,PartMPF,PEM,PDM,PartPosRef,Species USE MOD_Particle_Vars, ONLY:UseRotRefFrame,PartVeloRotRef +USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition #if defined(LSERK) USE MOD_Particle_Vars, ONLY:Pt_temp #endif @@ -864,6 +865,8 @@ SUBROUTINE MPIParticleSend(UseOldVecLength) ! Deallocate sendBuffer after send was successful, see MPIParticleRecv END DO ! iProc +IF(PDM%UNFPafterMPIPartSend) CALL UpdateNextFreePosition() + END SUBROUTINE MPIParticleSend @@ -890,6 +893,7 @@ SUBROUTINE MPIParticleRecv(DoMPIUpdateNextFreePos) USE MOD_Particle_Mesh_Vars ,ONLY: IsExchangeElem USE MOD_Particle_MPI_Vars ,ONLY: ExchangeProcToGlobalProc,DoParticleLatencyHiding USE MOD_Eval_xyz ,ONLY: GetPositionInRefElem +USE MOD_Part_Tools ,ONLY: GetNextFreePosition #if defined(LSERK) USE MOD_Particle_Vars ,ONLY: Pt_temp #endif @@ -926,7 +930,7 @@ SUBROUTINE MPIParticleRecv(DoMPIUpdateNextFreePos) ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: iProc, iPos, nRecv, PartID,jPos, iPart, TempNextFreePosition, ElemID +INTEGER :: iProc, iPos, nRecv, PartID,jPos, iPart, ElemID INTEGER :: recv_status_list(1:MPI_STATUS_SIZE,0:nExchangeProcessors-1) INTEGER :: MessageSize, nRecvParticles #if defined(ROS) || defined(IMPA) @@ -1020,9 +1024,7 @@ SUBROUTINE MPIParticleRecv(DoMPIUpdateNextFreePos) DO iPos=0,MessageSize-1-MsgLengthPoly - MsgLengthElec - MsgLengthAmbi,PartCommSize ! find free position in particle array nRecv = nRecv+1 - PartID = PDM%nextFreePosition(nRecv+PDM%CurrentNextFreePosition) - IF(PartID.EQ.0) CALL ABORT(__STAMP__,& - ' Error in ParticleExchange_parallel. PDM%nextFreePosition=0. Increase Part-MaxParticleNumber! ', nRecv) + PartID = GetNextFreePosition(nRecv) !>> particle position in physical space PartState(1:6,PartID) = PartRecvBuf(iProc)%content(1+iPos: 6+iPos) @@ -1304,21 +1306,28 @@ SUBROUTINE MPIParticleRecv(DoMPIUpdateNextFreePos) END DO ! iProc -TempNextFreePosition = PDM%CurrentNextFreePosition -PDM%ParticleVecLength = PDM%ParticleVecLength + PartMPIExchange%nMPIParticles -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + PartMPIExchange%nMPIParticles -PartMPIExchange%nMPIParticles = 0 -IF(PDM%ParticleVecLength.GT.PDM%MaxParticleNumber) CALL ABORT(__STAMP__& - ,' ParticleVecLegnth>MaxParticleNumber due to MPI-communication! Increase Part-maxParticleNumber or use more processors.') +IF(PartMPIExchange%nMPIParticles.GT.0) THEN + PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + PartMPIExchange%nMPIParticles + PDM%ParticleVecLength = MAX(PDM%ParticleVecLength,GetNextFreePosition(0)) +END IF +#ifdef CODE_ANALYZE +IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL Abort(__STAMP__,'PDM%ParticleVeclength exceeds PDM%maxParticleNumber, Difference:',IntInfoOpt=PDM%ParticleVeclength-PDM%maxParticleNumber) +DO PartID=PDM%ParticleVecLength+1,PDM%maxParticleNumber + IF (PDM%ParticleInside(PartID)) THEN + IPWRITE(*,*) PartID,PDM%ParticleVecLength,PDM%maxParticleNumber + CALL Abort(__STAMP__,'Particle outside PDM%ParticleVeclength',IntInfoOpt=PartID) + END IF +END DO +#endif IF(RadialWeighting%PerformCloning) THEN ! Checking whether received particles have to be cloned or deleted DO iPart = 1,nrecv - PartID = PDM%nextFreePosition(iPart+TempNextFreePosition) + PartID = GetNextFreePosition(iPart-PartMPIExchange%nMPIParticles) IF(ParticleOnProc(PartID)) CALL DSMC_2D_RadialWeighting(PartID,PEM%GlobalElemID(PartID)) END DO END IF - +PartMPIExchange%nMPIParticles = 0 ! deallocate send,receive buffer DO iProc=0,nExchangeProcessors-1 SDEALLOCATE(PartRecvBuf(iProc)%content) diff --git a/src/particles/particle_mpi/particle_mpi_emission.f90 b/src/particles/particle_mpi/particle_mpi_emission.f90 index 04e438d43..ce750a9b6 100644 --- a/src/particles/particle_mpi/particle_mpi_emission.f90 +++ b/src/particles/particle_mpi/particle_mpi_emission.f90 @@ -586,6 +586,7 @@ SUBROUTINE SendEmissionParticlesToProcs(chunkSize,DimSend,FractNbr,iInit,mySumOf USE MOD_Particle_MPI_Vars ,ONLY: EmissionSendBuf,EmissionRecvBuf USE MOD_Particle_Vars ,ONLY: PDM,PEM,PartState,PartPosRef,Species USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod +USE MOD_Part_Tools ,ONLY: GetNextFreePosition #if defined(MEASURE_MPI_WAIT) USE MOD_Particle_MPI_Vars ,ONLY: MPIW8TimePart,MPIW8CountPart #endif /*defined(MEASURE_MPI_WAIT)*/ @@ -911,21 +912,14 @@ SUBROUTINE SendEmissionParticlesToProcs(chunkSize,DimSend,FractNbr,iInit,mySumOf ! Located particle on local proc. ELSE ! Get the next free position in the PDM array - ParticleIndexNbr = PDM%nextFreePosition(mySumOfMatchedParticles + 1 + PDM%CurrentNextFreePosition) - IF (ParticleIndexNbr.NE.0) THEN - ! Fill the PartState manually to avoid a second localization - PartState(1:DimSend,ParticleIndexNbr) = particle_positions(DimSend*(i-1)+1:DimSend*(i-1)+DimSend) - PDM%ParticleInside( ParticleIndexNbr) = .TRUE. - IF (TrackingMethod.EQ.REFMAPPING) THEN - CALL GetPositionInRefElem(PartState(1:3,ParticleIndexNbr),PartPosRef(1:3,ParticleIndexNbr),ElemID) - END IF ! TrackingMethod.EQ.REFMAPPING - PEM%GlobalElemID(ParticleIndexNbr) = ElemID - ELSE - IPWRITE(UNIT_StdOut,'(I0,A,I0,A,I0,A)') " PDM%MaxParticleNumber = ", PDM%MaxParticleNumber," for each processor (",& - PDM%MaxParticleNumber*nProcessors," in total)" - IPWRITE(UNIT_StdOut,'(I0,A)') " Increase value for [Part-maxParticleNumber]!" - CALL ABORT(__STAMP__,'ERROR in ParticleMPIEmission:ParticleIndexNbr.EQ.0 - maximum nbr of particles reached?') - END IF + ParticleIndexNbr = GetNextFreePosition(mySumOfMatchedParticles+1) + ! Fill the PartState manually to avoid a second localization + PartState(1:DimSend,ParticleIndexNbr) = particle_positions(DimSend*(i-1)+1:DimSend*(i-1)+DimSend) + PDM%ParticleInside( ParticleIndexNbr) = .TRUE. + IF (TrackingMethod.EQ.REFMAPPING) THEN + CALL GetPositionInRefElem(PartState(1:3,ParticleIndexNbr),PartPosRef(1:3,ParticleIndexNbr),ElemID) + END IF ! TrackingMethod.EQ.REFMAPPING + PEM%GlobalElemID(ParticleIndexNbr) = ElemID mySumOfMatchedParticles = mySumOfMatchedParticles + 1 END IF ! ElemID.EQ.-1 END IF ! InsideMyBGM(i) @@ -1069,21 +1063,14 @@ SUBROUTINE SendEmissionParticlesToProcs(chunkSize,DimSend,FractNbr,iInit,mySumOf IF (ElemInfo_Shared(ELEM_RANK,ElemID).NE.myRank) CYCLE ! Find a free position in the PDM array - ParticleIndexNbr = PDM%nextFreePosition(mySumOfMatchedParticles + 1 + PDM%CurrentNextFreePosition) - IF (ParticleIndexNbr.NE.0) THEN - ! Fill the PartState manually to avoid a second localization - PartState(1:3,ParticleIndexNbr) = recvPartPos(DimSend*(i-1)+1:DimSend*(i-1)+3) - PDM%ParticleInside( ParticleIndexNbr) = .TRUE. - IF (TrackingMethod.EQ.REFMAPPING) THEN - CALL GetPositionInRefElem(PartState(1:3,ParticleIndexNbr),PartPosRef(1:3,ParticleIndexNbr),ElemID) - END IF ! TrackingMethod.EQ.REFMAPPING - PEM%GlobalElemID(ParticleIndexNbr) = ElemID - ELSE - IPWRITE(UNIT_StdOut,'(I0,A,I0,A,I0,A)') " PDM%MaxParticleNumber = ", PDM%MaxParticleNumber," for each processor (",& - PDM%MaxParticleNumber*nProcessors," in total)" - IPWRITE(UNIT_StdOut,'(I0,A)') " Increase value for [Part-maxParticleNumber]!" - CALL ABORT(__STAMP__,'ERROR in ParticleMPIEmission:ParticleIndexNbr.EQ.0 - maximum nbr of particles reached?') - END IF + ParticleIndexNbr = GetNextFreePosition(mySumOfMatchedParticles+1) + ! Fill the PartState manually to avoid a second localization + PartState(1:3,ParticleIndexNbr) = recvPartPos(DimSend*(i-1)+1:DimSend*(i-1)+3) + PDM%ParticleInside( ParticleIndexNbr) = .TRUE. + IF (TrackingMethod.EQ.REFMAPPING) THEN + CALL GetPositionInRefElem(PartState(1:3,ParticleIndexNbr),PartPosRef(1:3,ParticleIndexNbr),ElemID) + END IF ! TrackingMethod.EQ.REFMAPPING + PEM%GlobalElemID(ParticleIndexNbr) = ElemID mySumOfMatchedParticles = mySumOfMatchedParticles + 1 END DO @@ -1116,21 +1103,14 @@ SUBROUTINE SendEmissionParticlesToProcs(chunkSize,DimSend,FractNbr,iInit,mySumOf DO i = 1,PartMPILocate%nPartsRecv(1,iProc) ! Find a free position in the PDM array - ParticleIndexNbr = PDM%nextFreePosition(mySumOfMatchedParticles + 1 + PDM%CurrentNextFreePosition) - IF (ParticleIndexNbr.NE.0) THEN - ! Fill the PartState manually to avoid a second localization - PartState(1:3,ParticleIndexNbr) = EmissionRecvBuf(iProc)%content(PartCommSize*(i-1)+1:PartCommSize*(i-1)+3) - IF (TrackingMethod.EQ.REFMAPPING) THEN - PartPosRef(1:3,ParticleIndexNbr) = EmissionRecvBuf(iProc)%content(PartCommSize*(i-1)+4:PartCommSize*(i-1)+6) - END IF ! TrackingMethod.EQ.REFMAPPING - PEM%GlobalElemID(ParticleIndexNbr) = INT(EmissionRecvBuf(iProc)%content(PartCommSize*(i)),KIND=4) - PDM%ParticleInside( ParticleIndexNbr) = .TRUE. - ELSE - IPWRITE(UNIT_StdOut,'(I0,A,I0,A,I0,A)') " PDM%MaxParticleNumber = ", PDM%MaxParticleNumber," for each processor (",& - PDM%MaxParticleNumber*nProcessors," in total)" - IPWRITE(UNIT_StdOut,'(I0,A)') " Increase value for [Part-maxParticleNumber]!" - CALL ABORT(__STAMP__,'ERROR in ParticleMPIEmission:ParticleIndexNbr.EQ.0 - maximum nbr of particles reached?') - END IF + ParticleIndexNbr = GetNextFreePosition(mySumOfMatchedParticles+1) + ! Fill the PartState manually to avoid a second localization + PartState(1:3,ParticleIndexNbr) = EmissionRecvBuf(iProc)%content(PartCommSize*(i-1)+1:PartCommSize*(i-1)+3) + IF (TrackingMethod.EQ.REFMAPPING) THEN + PartPosRef(1:3,ParticleIndexNbr) = EmissionRecvBuf(iProc)%content(PartCommSize*(i-1)+4:PartCommSize*(i-1)+6) + END IF ! TrackingMethod.EQ.REFMAPPING + PEM%GlobalElemID(ParticleIndexNbr) = INT(EmissionRecvBuf(iProc)%content(PartCommSize*(i)),KIND=4) + PDM%ParticleInside( ParticleIndexNbr) = .TRUE. mySumOfMatchedParticles = mySumOfMatchedParticles + 1 END DO END DO diff --git a/src/particles/particle_operations.f90 b/src/particles/particle_operations.f90 index 41ab0bda3..1b1f74543 100644 --- a/src/particles/particle_operations.f90 +++ b/src/particles/particle_operations.f90 @@ -46,7 +46,7 @@ SUBROUTINE CreateParticle(SpecID,Pos,GlobElemID,Velocity,RotEnergy,VibEnergy,Ele USE MOD_Eval_xyz ,ONLY: GetPositionInRefElem USE MOD_part_tools ,ONLY: CalcRadWeightMPF USE MOD_Particle_TimeStep ,ONLY: GetParticleTimeStep -USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck +USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck, GetNextFreePosition !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT / OUTPUT VARIABLES @@ -64,15 +64,7 @@ SUBROUTINE CreateParticle(SpecID,Pos,GlobElemID,Velocity,RotEnergy,VibEnergy,Ele INTEGER :: newParticleID !=================================================================================================================================== -! Do not increase the ParticleVecLength for Phantom particles! -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + 1 -newParticleID = PDM%nextFreePosition(PDM%CurrentNextFreePosition) -IF(newParticleID.GT.PDM%ParticleVecLength) PDM%ParticleVecLength = PDM%ParticleVecLength + 1 - -IF(newParticleID.GT.PDM%MaxParticleNumber)THEN - CALL abort(__STAMP__,'CreateParticle: newParticleID.GT.PDM%MaxParticleNumber. '//& - 'Increase Part-maxParticleNumber or use more processors. newParticleID=',IntInfoOpt=newParticleID) -END IF +newParticleID = GetNextFreePosition() PartSpecies(newParticleID) = SpecID LastPartPos(1:3,newParticleID) = Pos(1:3) diff --git a/src/particles/particle_tools.f90 b/src/particles/particle_tools.f90 index 94b31e45a..5be240497 100644 --- a/src/particles/particle_tools.f90 +++ b/src/particles/particle_tools.f90 @@ -78,6 +78,7 @@ MODULE MOD_part_tools PUBLIC :: MergeCells,InRotRefFrameCheck PUBLIC :: CalcPartSymmetryPos PUBLIC :: RotateVectorAroundAxis +PUBLIC :: IncreaseMaxParticleNumber, GetNextFreePosition, ReduceMaxParticleNumber !=================================================================================================================================== CONTAINS @@ -214,12 +215,17 @@ SUBROUTINE UpdateNextFreePosition(WithOutMPIParts) IF (CollInf%ProhibitDoubleColl) CollInf%OldCollPartner(i) = 0 counter = counter + 1 PDM%nextFreePosition(counter) = i +#ifdef CODE_ANALYZE + IF(PDM%ParticleInside(i)) CALL ABORT(& + __STAMP__& + ,'Particle Inside is true but outside of PDM%ParticleVecLength',IntInfoOpt=i) +#endif END DO ! Set nextFreePosition for occupied slots to zero PDM%nextFreePosition(counter+1:PDM%maxParticleNumber) = 0 ! If maxParticleNumber are inside, counter is greater than maxParticleNumber -IF (counter+1.GT.PDM%MaxParticleNumber) PDM%nextFreePosition(PDM%MaxParticleNumber) = 0 +! IF (counter+1.GT.PDM%MaxParticleNumber) PDM%nextFreePosition(PDM%MaxParticleNumber) = 0 #if USE_LOADBALANCE CALL LBPauseTime(LB_UNFP,tLBStart) @@ -1038,7 +1044,7 @@ END FUNCTION CalcEElec_particle SUBROUTINE MergeCells() !=================================================================================================================================== -!> Routine for virtual merging of neighbouring cells. +!> Routine for virtual merging of neighbouring cells. !> Currently, the merging is only done via the number of particles within the cells. !=================================================================================================================================== ! MODULES @@ -1062,7 +1068,7 @@ SUBROUTINE MergeCells() !Nullify every value DO iElem = 1, nElems VirtMergedCells(iElem)%isMerged = .FALSE. - VirtMergedCells(iElem)%MasterCell = 0 + VirtMergedCells(iElem)%MasterCell = 0 VirtMergedCells(iElem)%MergedVolume = 0.0 IF (VirtMergedCells(iElem)%NumOfMergedCells.GT.0) THEN DEALLOCATE(VirtMergedCells(iElem)%MergedCellID) @@ -1088,7 +1094,7 @@ SUBROUTINE MergeCells() IF(VirtualCellMergeSpread.GT.1) THEN IF (VirtMergedCells(iElem)%NumOfMergedCells.GT.0) THEN MasterCellID = VirtMergedCells(LocNBElem)%MasterCell-offSetElem - IF(VirtMergedCells(MasterCellID)%NumOfMergedCells.GE.(MaxNumOfMergedCells-1)) CYCLE NBElemLoop + IF(VirtMergedCells(MasterCellID)%NumOfMergedCells.GE.(MaxNumOfMergedCells-1)) CYCLE NBElemLoop ALLOCATE(tempCellID(VirtMergedCells(MasterCellID)%NumOfMergedCells)) tempCellID = VirtMergedCells(MasterCellID)%MergedCellID DEALLOCATE(VirtMergedCells(MasterCellID)%MergedCellID) @@ -1121,9 +1127,9 @@ SUBROUTINE MergeCells() VirtMergedCells(MasterCellID)%MergedCellID(1:VirtMergedCells(MasterCellID)%NumOfMergedCells-1) = & tempCellID(1:VirtMergedCells(MasterCellID)%NumOfMergedCells-1) VirtMergedCells(MasterCellID)%MergedCellID(VirtMergedCells(MasterCellID)%NumOfMergedCells) = iElem - VirtMergedCells(MasterCellID)%MergedVolume=VirtMergedCells(MasterCellID)%MergedVolume+ElemVolume_Shared(CNElemID) + VirtMergedCells(MasterCellID)%MergedVolume=VirtMergedCells(MasterCellID)%MergedVolume+ElemVolume_Shared(CNElemID) VirtMergedCells(iElem)%MasterCell = VirtMergedCells(LocNBElem)%MasterCell - VirtMergedCells(iElem)%isMerged = .TRUE. + VirtMergedCells(iElem)%isMerged = .TRUE. DEALLOCATE(tempCellID) CYCLE ElemLoop END IF @@ -1167,7 +1173,7 @@ SUBROUTINE MergeCells() VirtMergedCells(MasterCellID)%MergedCellID(1:VirtMergedCells(MasterCellID)%NumOfMergedCells-1) = & tempCellID(1:VirtMergedCells(MasterCellID)%NumOfMergedCells-1) VirtMergedCells(MasterCellID)%MergedCellID(VirtMergedCells(MasterCellID)%NumOfMergedCells) = iElem - VirtMergedCells(MasterCellID)%MergedVolume=VirtMergedCells(MasterCellID)%MergedVolume+ElemVolume_Shared(CNElemID) + VirtMergedCells(MasterCellID)%MergedVolume=VirtMergedCells(MasterCellID)%MergedVolume+ElemVolume_Shared(CNElemID) VirtMergedCells(iElem)%MasterCell = MasterCellID + offSetElem VirtMergedCells(iElem)%isMerged = .TRUE. DEALLOCATE(tempCellID) @@ -1188,8 +1194,8 @@ SUBROUTINE MergeCells() ALLOCATE(VirtMergedCells(iElem)%MergedCellID(VirtMergedCells(iElem)%NumOfMergedCells)) VirtMergedCells(iElem)%MergedCellID(VirtMergedCells(iElem)%NumOfMergedCells) = LocNBElem VirtMergedCells(iElem)%MergedVolume = VirtMergedCells(iElem)%MergedVolume + ElemVolume_Shared(CNNbElem) - VirtMergedCells(iElem)%MasterCell = iElem + offSetElem - VirtMergedCells(LocNBElem)%MasterCell = iElem + offSetElem + VirtMergedCells(iElem)%MasterCell = iElem + offSetElem + VirtMergedCells(LocNBElem)%MasterCell = iElem + offSetElem VirtMergedCells(LocNBElem)%isMerged = .TRUE. ELSE IF(VirtMergedCells(iElem)%NumOfMergedCells.GE.(MaxNumOfMergedCells-1)) CYCLE ElemLoop @@ -1202,14 +1208,14 @@ SUBROUTINE MergeCells() tempCellID(1:VirtMergedCells(iElem)%NumOfMergedCells-1) VirtMergedCells(iElem)%MergedCellID(VirtMergedCells(iElem)%NumOfMergedCells) = LocNBElem VirtMergedCells(iElem)%MergedVolume = VirtMergedCells(iElem)%MergedVolume + ElemVolume_Shared(CNNbElem) - VirtMergedCells(LocNBElem)%MasterCell = iElem + offSetElem + VirtMergedCells(LocNBElem)%MasterCell = iElem + offSetElem VirtMergedCells(LocNBElem)%isMerged = .TRUE. DEALLOCATE(tempCellID) END IF nPartMerged = nPartMerged + PEM%pNumber(LocNBElem) IF (nPartMerged.GT.MinPartNumCellMerge) CYCLE ElemLoop END IF - END DO NBElemLoop + END DO NBElemLoop END IF END DO ElemLoop @@ -1608,4 +1614,792 @@ PPURE FUNCTION RotateVectorAroundAxis(VecIn,Axis,Angle) END FUNCTION RotateVectorAroundAxis + +FUNCTION GetNextFreePosition(Offset) +!=================================================================================================================================== +!> Returns the next free position in the particle vector, if no space is available it increses the maximum particle number +!> ATTENTION: If optional argument is used, the PDM%CurrentNextFreePosition will not be updated +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Particle_Vars ,ONLY: PDM +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,OPTIONAL,INTENT(IN) :: Offset +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +INTEGER :: GetNextFreePosition +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: i +!=================================================================================================================================== +IF(PRESENT(Offset)) THEN + ! IF(PDM%CurrentNextFreePosition+Offset.GT.PDM%MaxParticleNumber) CALL IncreaseMaxParticleNumber(CEILING((PDM%CurrentNextFreePosition+Offset)*(1+PDM%MaxPartNumIncrease)-PDM%MaxParticleNumber)) + IF(PDM%CurrentNextFreePosition+Offset.GT.PDM%MaxParticleNumber) THEN + CALL IncreaseMaxParticleNumber() + IF(PDM%CurrentNextFreePosition.GT.PDM%MaxParticleNumber) THEN + ! This only happens if PDM%CurrentNextFreePosition+Offset is way off (which shouldn't happen) + IPWRITE(UNIT_stdOut,*) "WARNING: PDM%CurrentNextFreePosition+Offset is way off in particle_tools.f90 GetNextFreePosition(Offset), 1" + CALL IncreaseMaxParticleNumber(CEILING((PDM%CurrentNextFreePosition+Offset)*(1+PDM%MaxPartNumIncrease)-PDM%MaxParticleNumber)) + END IF + END IF + + GetNextFreePosition = PDM%nextFreePosition(PDM%CurrentNextFreePosition+Offset) + ! If next free position is equal 0, determine how much more particles are needed to get a position within the particle vector + IF(GetNextFreePosition.EQ.0) THEN + CALL IncreaseMaxParticleNumber() + GetNextFreePosition = PDM%nextFreePosition(PDM%CurrentNextFreePosition+Offset) + IF(GetNextFreePosition.EQ.0) THEN + ! This only happens if PDM%CurrentNextFreePosition+Offset is way off (which shouldn't happen) + IPWRITE(UNIT_stdOut,*) "WARNING: PDM%CurrentNextFreePosition+Offset is way off in particle_tools.f90 GetNextFreePosition(Offset), 2" + IF(PDM%nextFreePosition(1).EQ.0) THEN + i = 0 + ELSE + i = PDM%CurrentNextFreePosition+Offset + DO WHILE(PDM%nextFreePosition(i).EQ.0.AND.i.GT.0) + i = i - 1 + END DO + END IF + ! Increase the maxpartnum + margin + CALL IncreaseMaxParticleNumber(CEILING((PDM%CurrentNextFreePosition+Offset-i)*(1+PDM%MaxPartNumIncrease)+PDM%maxParticleNumber*PDM%MaxPartNumIncrease)) + GetNextFreePosition = PDM%nextFreePosition(PDM%CurrentNextFreePosition+Offset) + END IF + END IF +ELSE + PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + 1 + ! IF(PDM%CurrentNextFreePosition.GT.PDM%MaxParticleNumber) CALL IncreaseMaxParticleNumber(CEILING((PDM%CurrentNextFreePosition)*(1+PDM%MaxPartNumIncrease)-PDM%MaxParticleNumber)) + IF(PDM%CurrentNextFreePosition.GT.PDM%MaxParticleNumber) THEN + CALL IncreaseMaxParticleNumber() + IF(PDM%CurrentNextFreePosition.GT.PDM%MaxParticleNumber) THEN + ! This only happens if PDM%CurrentNextFreePosition is way off (which shouldn't happen) + IPWRITE(UNIT_stdOut,*) "WARNING: PDM%CurrentNextFreePosition is way off in particle_tools.f90 GetNextFreePosition(), 1" + CALL IncreaseMaxParticleNumber(CEILING((PDM%CurrentNextFreePosition)*(1+PDM%MaxPartNumIncrease)-PDM%MaxParticleNumber)) + END IF + END IF + + GetNextFreePosition = PDM%nextFreePosition(PDM%CurrentNextFreePosition) + ! If next free position is equal 0, determine how much more particles are needed to get a position within the particle vector + IF(GetNextFreePosition.EQ.0) THEN + CALL IncreaseMaxParticleNumber() + GetNextFreePosition = PDM%nextFreePosition(PDM%CurrentNextFreePosition) + IF(GetNextFreePosition.EQ.0) THEN + ! This only happens if PDM%CurrentNextFreePosition is way off (which shouldn't happen) + IPWRITE(UNIT_stdOut,*) "WARNING: PDM%CurrentNextFreePosition is way off in particle_tools.f90 GetNextFreePosition(), 2" + IF(PDM%nextFreePosition(1).EQ.0) THEN + i = 0 + ELSE + i = PDM%CurrentNextFreePosition + DO WHILE(PDM%nextFreePosition(i).EQ.0.AND.i.GT.0) + i = i - 1 + END DO + END IF + ! Increase the maxpartnum + margin + CALL IncreaseMaxParticleNumber(CEILING((PDM%CurrentNextFreePosition-i)*(1+PDM%MaxPartNumIncrease)+PDM%maxParticleNumber*PDM%MaxPartNumIncrease)) + GetNextFreePosition = PDM%nextFreePosition(PDM%CurrentNextFreePosition) + END IF + END IF + + IF(PDM%ParticleInside(GetNextFreePosition)) THEN + CALL ABORT(& + __STAMP__& + ,'This Particle is already in use',IntInfoOpt=GetNextFreePosition) + END IF + IF(GetNextFreePosition.GT.PDM%ParticleVecLength) PDM%ParticleVecLength = GetNextFreePosition +END IF +IF(GetNextFreePosition.EQ.0) THEN + CALL ABORT(& +__STAMP__& +,'This should not happen, PDM%MaxParticleNumber reached',IntInfoOpt=PDM%MaxParticleNumber) +END IF + +END FUNCTION GetNextFreePosition + + +SUBROUTINE IncreaseMaxParticleNumber(Amount) +!=================================================================================================================================== +! Increases MaxParticleNumber and increases size of all depended arrays +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Array_Operations ,ONLY: ChangeSizeArray +USE MOD_Particle_Vars +USE MOD_DSMC_Vars +#if USE_MPI +USE MOD_Particle_MPI_Vars ,ONLY: PartShiftVector, PartTargetProc +#endif +USE MOD_PICInterpolation_Vars ,ONLY: FieldAtParticle +#if defined(IMPA) || defined(ROS) +USE MOD_LinearSolver_Vars ,ONLY: PartXK, R_PartXK +USE MOD_TimeDisc_Vars ,ONLY: nRKStages +#endif +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,INTENT(IN),OPTIONAL :: Amount +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: NewSize, i, ii, ALLOCSTAT +TYPE (tAmbipolElecVelo), ALLOCATABLE :: AmbipolElecVelo_New(:) +TYPE (tElectronicDistriPart), ALLOCATABLE :: ElectronicDistriPart_New(:) +TYPE (tPolyatomMolVibQuant), ALLOCATABLE :: VibQuantsPar_New(:) +TYPE (tClonedParticles), ALLOCATABLE :: ClonedParticles_New(:,:) +! REAL :: +!=================================================================================================================================== +IF(PRESENT(Amount)) THEN + IF(Amount.EQ.0) RETURN + NewSize=PDM%MaxParticleNumber+Amount + ! IPWRITE(*,*) "Increase by amount",PDM%MaxParticleNumber,NewSize + IF(NewSize.GT.PDM%maxAllowedParticleNumber)CALL ABORT(& + __STAMP__& + ,'More Particles needed than allowed in PDM%maxAllowedParticleNumber',IntInfoOpt=NewSize) +ELSE + NewSize=MAX(CEILING(PDM%MaxParticleNumber*(1+PDM%MaxPartNumIncrease)),PDM%MaxParticleNumber+1) + IF(PDM%MaxParticleNumber.GE.PDM%maxAllowedParticleNumber) CALL ABORT(& + __STAMP__& + ,'More Particles needed than allowed in PDM%maxAllowedParticleNumber',IntInfoOpt=NewSize) + NewSize=MIN(NewSize,PDM%maxAllowedParticleNumber) + ! IPWRITE(*,*) "Increase by percent",PDM%MaxParticleNumber,NewSize +END IF + +IF(ALLOCATED(PEM%GlobalElemID)) CALL ChangeSizeArray(PEM%GlobalElemID,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PEM%pNext)) CALL ChangeSizeArray(PEM%pNext,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PEM%LastGlobalElemID)) CALL ChangeSizeArray(PEM%LastGlobalElemID,PDM%maxParticleNumber,NewSize) + +IF(ALLOCATED(PDM%ParticleInside)) CALL ChangeSizeArray(PDM%ParticleInside,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PDM%IsNewPart)) CALL ChangeSizeArray(PDM%IsNewPart,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PDM%dtFracPush)) CALL ChangeSizeArray(PDM%dtFracPush,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PDM%InRotRefFrame)) CALL ChangeSizeArray(PDM%InRotRefFrame,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PDM%PartInit)) CALL ChangeSizeArray(PDM%PartInit,PDM%maxParticleNumber,NewSize) + +IF(ALLOCATED(PartState)) CALL ChangeSizeArray(PartState,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(LastPartPos)) CALL ChangeSizeArray(LastPartPos,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartPosRef)) CALL ChangeSizeArray(PartPosRef,PDM%maxParticleNumber,NewSize,-888.) +IF(ALLOCATED(PartSpecies)) CALL ChangeSizeArray(PartSpecies,PDM%maxParticleNumber,NewSize,0) +IF(ALLOCATED(PartTimeStep)) CALL ChangeSizeArray(PartTimeStep,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartMPF)) CALL ChangeSizeArray(PartMPF,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartVeloRotRef)) CALL ChangeSizeArray(PartVeloRotRef,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(PartStateIntEn)) CALL ChangeSizeArray(PartStateIntEn,PDM%maxParticleNumber,NewSize) + +IF(ALLOCATED(Pt_temp)) CALL ChangeSizeArray(Pt_temp,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(Pt)) CALL ChangeSizeArray(Pt,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(FieldAtParticle)) CALL ChangeSizeArray(FieldAtParticle,PDM%maxParticleNumber,NewSize) + +IF(ALLOCATED(InterPlanePartIndx)) CALL ChangeSizeArray(InterPlanePartIndx,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(BGGas%PairingPartner)) CALL ChangeSizeArray(BGGas%PairingPartner,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(CollInf%OldCollPartner)) CALL ChangeSizeArray(CollInf%OldCollPartner,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(ElecRelaxPart)) CALL ChangeSizeArray(ElecRelaxPart,PDM%maxParticleNumber,NewSize) + +#if (PP_TimeDiscMethod==508) || (PP_TimeDiscMethod==509) +IF(ALLOCATED(velocityAtTime)) CALL ChangeSizeArray(velocityAtTime,PDM%maxParticleNumber,NewSize) +#endif + +#if USE_MPI +IF(ALLOCATED(PartTargetProc)) CALL ChangeSizeArray(PartTargetProc,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartShiftVector)) CALL ChangeSizeArray(PartShiftVector,PDM%maxParticleNumber,NewSize) +#endif + +#if defined(IMPA) || defined(ROS) +IF(ALLOCATED(PartXK)) CALL ChangeSizeArray(PartXK,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(R_PartXK)) CALL ChangeSizeArray(R_PartXK,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartStage)) CALL ChangeSizeArray(PartStage,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartStateN)) CALL ChangeSizeArray(PartStateN,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartQ)) CALL ChangeSizeArray(PartQ,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PEM%NormVec)) CALL ChangeSizeArray(PEM%NormVec,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(PEM%PeriodicMoved)) CALL ChangeSizeArray(PEM%PeriodicMoved,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PartDtFrac)) CALL ChangeSizeArray(PartDtFrac,PDM%maxParticleNumber,NewSize,1.) +#endif + +#ifdef IMPA +IF(ALLOCATED(F_PartX0)) CALL ChangeSizeArray(F_PartX0,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(F_PartXk)) CALL ChangeSizeArray(F_PartXk,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(Norm_F_PartX0)) CALL ChangeSizeArray(Norm_F_PartX0,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(Norm_F_PartXk)) CALL ChangeSizeArray(Norm_F_PartXk,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(Norm_F_PartXk_old)) CALL ChangeSizeArray(Norm_F_PartXk_old,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartDeltaX)) CALL ChangeSizeArray(PartDeltaX,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartLambdaAccept)) CALL ChangeSizeArray(PartLambdaAccept,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(DoPartInNewton)) CALL ChangeSizeArray(DoPartInNewton,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartIsImplicit)) CALL ChangeSizeArray(PartIsImplicit,PDM%maxParticleNumber,NewSize,.FALSE.) +#endif /* IMPA */ + +! __ __ __ __ ________ ______ ___________ ___ +! / / / /___ ____/ /___ _/ /____ /_ __/\ \/ / __ \/ ____/ ___/ / | ______________ ___ _______ +! / / / / __ \/ __ / __ `/ __/ _ \ / / \ / /_/ / __/ \__ \ / /| | / ___/ ___/ __ `/ / / / ___/ +! / /_/ / /_/ / /_/ / /_/ / /_/ __/ / / / / ____/ /___ ___/ / / ___ |/ / / / / /_/ / /_/ (__ ) +! \____/ .___/\__,_/\__,_/\__/\___/ /_/ /_/_/ /_____//____/ /_/ |_/_/ /_/ \__,_/\__, /____/ +! /_/ /____/ + +IF(ALLOCATED(AmbipolElecVelo)) THEN + ALLOCATE(AmbipolElecVelo_New(NewSize),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate increased Array in IncreaseMaxParticleNumber') + DO i=1,PDM%maxParticleNumber + CALL MOVE_ALLOC(AmbipolElecVelo(i)%ElecVelo,AmbipolElecVelo_New(i)%ElecVelo) + END DO + DEALLOCATE(AmbipolElecVelo) + CALL MOVE_ALLOC(AmbipolElecVelo_New,AmbipolElecVelo) +END IF + +IF(ALLOCATED(ElectronicDistriPart)) THEN + ALLOCATE(ElectronicDistriPart_New(NewSize),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate increased Array in IncreaseMaxParticleNumber') + DO i=1,PDM%maxParticleNumber + CALL MOVE_ALLOC(ElectronicDistriPart(i)%DistriFunc,ElectronicDistriPart_New(i)%DistriFunc) + END DO + DEALLOCATE(ElectronicDistriPart) + CALL MOVE_ALLOC(ElectronicDistriPart_New,ElectronicDistriPart) +END IF + +IF(ALLOCATED(VibQuantsPar)) THEN + ALLOCATE(VibQuantsPar_New(NewSize),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate increased Array in IncreaseMaxParticleNumber') + DO i=1,PDM%maxParticleNumber + CALL MOVE_ALLOC(VibQuantsPar(i)%Quants,VibQuantsPar_New(i)%Quants) + END DO + DEALLOCATE(VibQuantsPar) + CALL MOVE_ALLOC(VibQuantsPar_New,VibQuantsPar) +END IF + +IF(ALLOCATED(ClonedParticles)) THEN + SELECT CASE(RadialWeighting%CloneMode) + CASE(1) + ALLOCATE(ClonedParticles_new(1:INT(NewSize/RadialWeighting%CloneInputDelay),0:(RadialWeighting%CloneInputDelay-1)),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& + __STAMP__& + ,'Cannot allocate increased Array in IncreaseMaxParticleNumber') + DO ii=0,RadialWeighting%CloneInputDelay-1 + DO i=1,RadialWeighting%ClonePartNum(ii) + ClonedParticles_new(i,ii)%Species=ClonedParticles(i,ii)%Species + ClonedParticles_new(i,ii)%PartState(1:6)=ClonedParticles(i,ii)%PartState(1:6) + ClonedParticles_new(i,ii)%PartStateIntEn(1:3)=ClonedParticles(i,ii)%PartStateIntEn(1:3) + ClonedParticles_new(i,ii)%Element=ClonedParticles(i,ii)%Element + ClonedParticles_new(i,ii)%LastPartPos(1:3)=ClonedParticles(i,ii)%LastPartPos(1:3) + ClonedParticles_new(i,ii)%WeightingFactor=ClonedParticles(i,ii)%WeightingFactor + CALL MOVE_ALLOC(ClonedParticles(i,ii)%VibQuants,ClonedParticles_new(i,ii)%VibQuants) + CALL MOVE_ALLOC(ClonedParticles(i,ii)%DistriFunc,ClonedParticles_new(i,ii)%DistriFunc) + CALL MOVE_ALLOC(ClonedParticles(i,ii)%AmbiPolVelo,ClonedParticles_new(i,ii)%AmbiPolVelo) + END DO + END DO + DEALLOCATE(ClonedParticles) + CALL MOVE_ALLOC(ClonedParticles_New,ClonedParticles) + CASE(2) + ALLOCATE(ClonedParticles_new(1:INT(NewSize/RadialWeighting%CloneInputDelay),0:RadialWeighting%CloneInputDelay),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& + __STAMP__& + ,'Cannot allocate increased Array in IncreaseMaxParticleNumber') + DO ii=0,RadialWeighting%CloneInputDelay + DO i=1,RadialWeighting%ClonePartNum(ii) + ClonedParticles_new(i,ii)%Species=ClonedParticles(i,ii)%Species + ClonedParticles_new(i,ii)%PartState(1:6)=ClonedParticles(i,ii)%PartState(1:6) + ClonedParticles_new(i,ii)%PartStateIntEn(1:3)=ClonedParticles(i,ii)%PartStateIntEn(1:3) + ClonedParticles_new(i,ii)%Element=ClonedParticles(i,ii)%Element + ClonedParticles_new(i,ii)%LastPartPos(1:3)=ClonedParticles(i,ii)%LastPartPos(1:3) + ClonedParticles_new(i,ii)%WeightingFactor=ClonedParticles(i,ii)%WeightingFactor + CALL MOVE_ALLOC(ClonedParticles(i,ii)%VibQuants,ClonedParticles_new(i,ii)%VibQuants) + CALL MOVE_ALLOC(ClonedParticles(i,ii)%DistriFunc,ClonedParticles_new(i,ii)%DistriFunc) + CALL MOVE_ALLOC(ClonedParticles(i,ii)%AmbiPolVelo,ClonedParticles_new(i,ii)%AmbiPolVelo) + END DO + END DO + DEALLOCATE(ClonedParticles) + CALL MOVE_ALLOC(ClonedParticles_New,ClonedParticles) + CASE DEFAULT + CALL Abort(& + __STAMP__,& + 'ERROR in Radial Weighting of 2D/Axisymmetric: The selected cloning mode is not available! Choose between 1 and 2.'//& + ' CloneMode=1: Delayed insertion of clones; CloneMode=2: Delayed randomized insertion of clones') +END SELECT +END IF + + +IF(ALLOCATED(PDM%nextFreePosition)) THEN + CALL ChangeSizeArray(PDM%nextFreePosition,PDM%maxParticleNumber,NewSize,0) + + !Search for first entry where new poition is available + i=1 + DO WHILE(PDM%nextFreePosition(i).NE.0) + i=i+1 + END DO + i=i-1 + ! Fill the free spots with the new entrys + DO ii=1,NewSize-PDM%MaxParticleNumber + PDM%nextFreePosition(i+ii)=ii+PDM%MaxParticleNumber + END DO +END IF + +PDM%MaxParticleNumber=NewSize + +END SUBROUTINE IncreaseMaxParticleNumber + + +SUBROUTINE ReduceMaxParticleNumber() +!=================================================================================================================================== +! Reduces MaxParticleNumber and increases size of all depended arrays +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Array_Operations ,ONLY: ChangeSizeArray +USE MOD_Particle_Vars +USE MOD_DSMC_Vars +#if USE_MPI +USE MOD_Particle_MPI_Vars ,ONLY: PartShiftVector, PartTargetProc +#endif +USE MOD_PICInterpolation_Vars ,ONLY: FieldAtParticle +#if defined(IMPA) || defined(ROS) +USE MOD_LinearSolver_Vars ,ONLY: PartXK, R_PartXK +USE MOD_TimeDisc_Vars ,ONLY: nRKStages +#endif +USE MOD_MCC_Vars ,ONLY: UseMCC +USE MOD_DSMC_Vars ,ONLY: BGGas +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES + +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: NewSize, i, ii, ALLOCSTAT, nPart +TYPE (tAmbipolElecVelo), ALLOCATABLE :: AmbipolElecVelo_New(:) +TYPE (tElectronicDistriPart), ALLOCATABLE :: ElectronicDistriPart_New(:) +TYPE (tPolyatomMolVibQuant), ALLOCATABLE :: VibQuantsPar_New(:) +TYPE (tClonedParticles), ALLOCATABLE :: ClonedParticles_New(:,:) +! REAL :: +!=================================================================================================================================== + +nPart=0 +DO i=1,PDM%ParticleVecLength + IF(PDM%ParticleInside(i)) nPart = nPart + 1 +END DO + +IF(DSMC%DoAmbipolarDiff) THEN + DO i=1,PDM%ParticleVecLength + IF(PDM%ParticleInside(i).AND.ALLOCATED(AmbipolElecVelo(i)%ElecVelo)) nPart = nPart + 1 + END DO +END IF + +IF(BGGas%NumberOfSpecies.GT.0.AND..NOT.UseMCC) nPart=nPart*2 + +! Reduce Arrays only for at least PDM%maxParticleNumber*PDM%MaxPartNumIncrease free spots +IF (nPart.GE.PDM%maxParticleNumber/(1.+PDM%MaxPartNumIncrease)**2) RETURN + +! Maintain nPart*PDM%MaxPartNumIncrease free spots +Newsize=MAX(CEILING(nPart*(1.+PDM%MaxPartNumIncrease)),1) +IF (Newsize.EQ.PDM%maxParticleNumber) RETURN + +! IPWRITE(*,*) "Decrease",PDM%maxParticleNumber,nPart,NewSize +IF(.NOT.PDM%RearrangePartIDs) THEN + ! Search for highest occupied particle index and set Newsize to this Value + i=PDM%maxParticleNumber + DO WHILE(.NOT.PDM%ParticleInside(i).OR.i.EQ.NewSize) + i=i-1 + END DO + NewSize=i +ELSE + ! Rearrange particles with IDs>NewSize to lower IDs + DO i=NewSize+1,PDM%maxParticleNumber + IF(PDM%ParticleInside(i)) THEN + PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + 1 + ii = PDM%nextFreePosition(PDM%CurrentNextFreePosition) + IF(ii.EQ.0.OR.ii.GT.NewSize) THEN + CALL UpdateNextFreePosition() + PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + 1 + ii = PDM%nextFreePosition(PDM%CurrentNextFreePosition) + IF(ii.EQ.0.OR.ii.GT.NewSize) CALL ABORT(& + __STAMP__& + ,'This should not happen') + END IF + IF(PDM%ParticleVecLength.LT.ii) PDM%ParticleVecLength = ii + CALL ChangePartID(i,ii) + END IF + END DO +END IF + + + +IF(ALLOCATED(PEM%GlobalElemID)) CALL ChangeSizeArray(PEM%GlobalElemID,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PEM%pNext)) CALL ChangeSizeArray(PEM%pNext,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PEM%LastGlobalElemID)) CALL ChangeSizeArray(PEM%LastGlobalElemID,PDM%maxParticleNumber,NewSize) + +IF(ALLOCATED(PDM%ParticleInside)) CALL ChangeSizeArray(PDM%ParticleInside,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PDM%IsNewPart)) CALL ChangeSizeArray(PDM%IsNewPart,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PDM%dtFracPush)) CALL ChangeSizeArray(PDM%dtFracPush,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PDM%InRotRefFrame)) CALL ChangeSizeArray(PDM%InRotRefFrame,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PDM%PartInit)) CALL ChangeSizeArray(PDM%PartInit,PDM%maxParticleNumber,NewSize) + +IF(ALLOCATED(PartState)) CALL ChangeSizeArray(PartState,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(LastPartPos)) CALL ChangeSizeArray(LastPartPos,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartPosRef)) CALL ChangeSizeArray(PartPosRef,PDM%maxParticleNumber,NewSize,-888.) +IF(ALLOCATED(PartSpecies)) CALL ChangeSizeArray(PartSpecies,PDM%maxParticleNumber,NewSize,0) +IF(ALLOCATED(PartTimeStep)) CALL ChangeSizeArray(PartTimeStep,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartMPF)) CALL ChangeSizeArray(PartMPF,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartVeloRotRef)) CALL ChangeSizeArray(PartVeloRotRef,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(PartStateIntEn)) CALL ChangeSizeArray(PartStateIntEn,PDM%maxParticleNumber,NewSize) + +IF(ALLOCATED(Pt_temp)) CALL ChangeSizeArray(Pt_temp,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(Pt)) CALL ChangeSizeArray(Pt,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(FieldAtParticle)) CALL ChangeSizeArray(FieldAtParticle,PDM%maxParticleNumber,NewSize) + +IF(ALLOCATED(InterPlanePartIndx)) CALL ChangeSizeArray(InterPlanePartIndx,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(BGGas%PairingPartner)) CALL ChangeSizeArray(BGGas%PairingPartner,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(CollInf%OldCollPartner)) CALL ChangeSizeArray(CollInf%OldCollPartner,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(ElecRelaxPart)) CALL ChangeSizeArray(ElecRelaxPart,PDM%maxParticleNumber,NewSize) + +#if (PP_TimeDiscMethod==508) || (PP_TimeDiscMethod==509) +IF(ALLOCATED(velocityAtTime)) CALL ChangeSizeArray(velocityAtTime,PDM%maxParticleNumber,NewSize) +#endif + +#if USE_MPI +IF(ALLOCATED(PartTargetProc)) CALL ChangeSizeArray(PartTargetProc,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartShiftVector)) CALL ChangeSizeArray(PartShiftVector,PDM%maxParticleNumber,NewSize) +#endif + +#if defined(IMPA) || defined(ROS) +IF(ALLOCATED(PartXK)) CALL ChangeSizeArray(PartXK,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(R_PartXK)) CALL ChangeSizeArray(R_PartXK,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartStage)) CALL ChangeSizeArray(PartStage,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartStateN)) CALL ChangeSizeArray(PartStateN,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartQ)) CALL ChangeSizeArray(PartQ,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PEM%NormVec)) CALL ChangeSizeArray(PEM%NormVec,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(PEM%PeriodicMoved)) CALL ChangeSizeArray(PEM%PeriodicMoved,PDM%maxParticleNumber,NewSize,.FALSE.) +IF(ALLOCATED(PartDtFrac)) CALL ChangeSizeArray(PartDtFrac,PDM%maxParticleNumber,NewSize,1.) +#endif + +#ifdef IMPA +IF(ALLOCATED(F_PartX0)) CALL ChangeSizeArray(F_PartX0,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(F_PartXk)) CALL ChangeSizeArray(F_PartXk,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(Norm_F_PartX0)) CALL ChangeSizeArray(Norm_F_PartX0,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(Norm_F_PartXk)) CALL ChangeSizeArray(Norm_F_PartXk,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(Norm_F_PartXk_old)) CALL ChangeSizeArray(Norm_F_PartXk_old,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartDeltaX)) CALL ChangeSizeArray(PartDeltaX,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartLambdaAccept)) CALL ChangeSizeArray(PartLambdaAccept,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(DoPartInNewton)) CALL ChangeSizeArray(DoPartInNewton,PDM%maxParticleNumber,NewSize) +IF(ALLOCATED(PartIsImplicit)) CALL ChangeSizeArray(PartIsImplicit,PDM%maxParticleNumber,NewSize,.FALSE.) +#endif /* IMPA */ + +! __ __ __ __ ________ ______ ___________ ___ +! / / / /___ ____/ /___ _/ /____ /_ __/\ \/ / __ \/ ____/ ___/ / | ______________ ___ _______ +! / / / / __ \/ __ / __ `/ __/ _ \ / / \ / /_/ / __/ \__ \ / /| | / ___/ ___/ __ `/ / / / ___/ +! / /_/ / /_/ / /_/ / /_/ / /_/ __/ / / / / ____/ /___ ___/ / / ___ |/ / / / / /_/ / /_/ (__ ) +! \____/ .___/\__,_/\__,_/\__/\___/ /_/ /_/_/ /_____//____/ /_/ |_/_/ /_/ \__,_/\__, /____/ +! /_/ /____/ + +IF(ALLOCATED(AmbipolElecVelo)) THEN + ALLOCATE(AmbipolElecVelo_New(NewSize),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate increased Array in ReduceMaxParticleNumber') + DO i=1,NewSize + CALL MOVE_ALLOC(AmbipolElecVelo(i)%ElecVelo,AmbipolElecVelo_New(i)%ElecVelo) + END DO + DO i=NewSize+1,PDM%maxParticleNumber + SDEALLOCATE(AmbipolElecVelo(i)%ElecVelo) + END DO + DEALLOCATE(AmbipolElecVelo) + CALL MOVE_ALLOC(AmbipolElecVelo_New,AmbipolElecVelo) +END IF + +IF(ALLOCATED(ElectronicDistriPart)) THEN + ALLOCATE(ElectronicDistriPart_New(NewSize),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate increased Array in ReduceMaxParticleNumber') + DO i=1,NewSize + CALL MOVE_ALLOC(ElectronicDistriPart(i)%DistriFunc,ElectronicDistriPart_New(i)%DistriFunc) + END DO + DO i=NewSize+1,PDM%maxParticleNumber + SDEALLOCATE(ElectronicDistriPart(i)%DistriFunc) + END DO + DEALLOCATE(ElectronicDistriPart) + CALL MOVE_ALLOC(ElectronicDistriPart_New,ElectronicDistriPart) +END IF + +IF(ALLOCATED(VibQuantsPar)) THEN + ALLOCATE(VibQuantsPar_New(NewSize),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& +__STAMP__& +,'Cannot allocate increased Array in ReduceMaxParticleNumber') + DO i=1,NewSize + CALL MOVE_ALLOC(VibQuantsPar(i)%Quants,VibQuantsPar_New(i)%Quants) + END DO + DO i=NewSize+1,PDM%maxParticleNumber + SDEALLOCATE(VibQuantsPar(i)%Quants) + END DO + DEALLOCATE(VibQuantsPar) + CALL MOVE_ALLOC(VibQuantsPar_New,VibQuantsPar) +END IF + +IF(ALLOCATED(ClonedParticles)) THEN + SELECT CASE(RadialWeighting%CloneMode) + CASE(1) + ALLOCATE(ClonedParticles_new(1:INT(NewSize/RadialWeighting%CloneInputDelay),0:(RadialWeighting%CloneInputDelay-1)),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& + __STAMP__& + ,'Cannot allocate increased Array in ReduceMaxParticleNumber') + DO ii=0,RadialWeighting%CloneInputDelay-1 + DO i=1,RadialWeighting%ClonePartNum(ii) + IF(i.LE.INT(NewSize/RadialWeighting%CloneInputDelay)) THEN + ClonedParticles_new(i,ii)%Species=ClonedParticles(i,ii)%Species + ClonedParticles_new(i,ii)%PartState(1:6)=ClonedParticles(i,ii)%PartState(1:6) + ClonedParticles_new(i,ii)%PartStateIntEn(1:3)=ClonedParticles(i,ii)%PartStateIntEn(1:3) + ClonedParticles_new(i,ii)%Element=ClonedParticles(i,ii)%Element + ClonedParticles_new(i,ii)%LastPartPos(1:3)=ClonedParticles(i,ii)%LastPartPos(1:3) + ClonedParticles_new(i,ii)%WeightingFactor=ClonedParticles(i,ii)%WeightingFactor + CALL MOVE_ALLOC(ClonedParticles(i,ii)%VibQuants,ClonedParticles_new(i,ii)%VibQuants) + CALL MOVE_ALLOC(ClonedParticles(i,ii)%DistriFunc,ClonedParticles_new(i,ii)%DistriFunc) + CALL MOVE_ALLOC(ClonedParticles(i,ii)%AmbiPolVelo,ClonedParticles_new(i,ii)%AmbiPolVelo) + ELSE + SDEALLOCATE(ClonedParticles(i,ii)%VibQuants) + SDEALLOCATE(ClonedParticles(i,ii)%DistriFunc) + SDEALLOCATE(ClonedParticles(i,ii)%AmbiPolVelo) + END IF + END DO + IF(RadialWeighting%ClonePartNum(ii).GT.INT(NewSize/RadialWeighting%CloneInputDelay)) & + RadialWeighting%ClonePartNum(ii)=INT(NewSize/RadialWeighting%CloneInputDelay) + END DO + DEALLOCATE(ClonedParticles) + CALL MOVE_ALLOC(ClonedParticles_New,ClonedParticles) + CASE(2) + ALLOCATE(ClonedParticles_new(1:INT(NewSize/RadialWeighting%CloneInputDelay),0:RadialWeighting%CloneInputDelay),STAT=ALLOCSTAT) + IF (ALLOCSTAT.NE.0) CALL ABORT(& + __STAMP__& + ,'Cannot allocate increased Array in ReduceMaxParticleNumber') + DO ii=0,RadialWeighting%CloneInputDelay + DO i=1,RadialWeighting%ClonePartNum(ii) + IF(i.LE.INT(NewSize/RadialWeighting%CloneInputDelay)) THEN + ClonedParticles_new(i,ii)%Species=ClonedParticles(i,ii)%Species + ClonedParticles_new(i,ii)%PartState(1:6)=ClonedParticles(i,ii)%PartState(1:6) + ClonedParticles_new(i,ii)%PartStateIntEn(1:3)=ClonedParticles(i,ii)%PartStateIntEn(1:3) + ClonedParticles_new(i,ii)%Element=ClonedParticles(i,ii)%Element + ClonedParticles_new(i,ii)%LastPartPos(1:3)=ClonedParticles(i,ii)%LastPartPos(1:3) + ClonedParticles_new(i,ii)%WeightingFactor=ClonedParticles(i,ii)%WeightingFactor + CALL MOVE_ALLOC(ClonedParticles(i,ii)%VibQuants,ClonedParticles_new(i,ii)%VibQuants) + CALL MOVE_ALLOC(ClonedParticles(i,ii)%DistriFunc,ClonedParticles_new(i,ii)%DistriFunc) + CALL MOVE_ALLOC(ClonedParticles(i,ii)%AmbiPolVelo,ClonedParticles_new(i,ii)%AmbiPolVelo) + ELSE + SDEALLOCATE(ClonedParticles(i,ii)%VibQuants) + SDEALLOCATE(ClonedParticles(i,ii)%DistriFunc) + SDEALLOCATE(ClonedParticles(i,ii)%AmbiPolVelo) + END IF + END DO + IF(RadialWeighting%ClonePartNum(ii).GT.INT(NewSize/RadialWeighting%CloneInputDelay)) & + RadialWeighting%ClonePartNum(ii)=INT(NewSize/RadialWeighting%CloneInputDelay) + END DO + DEALLOCATE(ClonedParticles) + CALL MOVE_ALLOC(ClonedParticles_New,ClonedParticles) + CASE DEFAULT + CALL Abort(& + __STAMP__,& + 'ERROR in Radial Weighting of 2D/Axisymmetric: The selected cloning mode is not available! Choose between 1 and 2.'//& + ' CloneMode=1: Delayed insertion of clones; CloneMode=2: Delayed randomized insertion of clones') +END SELECT +END IF + +IF(ALLOCATED(PDM%nextFreePosition)) THEN + CALL ChangeSizeArray(PDM%nextFreePosition,PDM%maxParticleNumber,NewSize,0) + + !Set all NextFreePositions to zero which points to a partID>NewSize + DO i=1,NewSize + IF(PDM%nextFreePosition(i).GT.NewSize) PDM%nextFreePosition(i)=0 + END DO +END IF + +IF(PDM%ParticleVecLength.GT.NewSize) PDM%ParticleVecLength = NewSize +PDM%MaxParticleNumber=NewSize + +CALL UpdateNextFreePosition() +! read(*,*) + +END SUBROUTINE ReduceMaxParticleNumber + + +SUBROUTINE ChangePartID(OldID,NewID) +!=================================================================================================================================== +! Change PartID from OldID to NewID +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Particle_Vars +USE MOD_DSMC_Vars +#if USE_MPI +USE MOD_Particle_MPI_Vars ,ONLY: PartShiftVector, PartTargetProc +#endif +USE MOD_PICInterpolation_Vars ,ONLY: FieldAtParticle +#if defined(IMPA) || defined(ROS) +USE MOD_LinearSolver_Vars ,ONLY: PartXK, R_PartXK +#endif +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,INTENT(IN) :: OldID +INTEGER,INTENT(IN) :: NewID +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: i,TempPartID +!=================================================================================================================================== + +IF(ALLOCATED(PEM%GlobalElemID)) PEM%GlobalElemID(NewID)=PEM%GlobalElemID(OldID) +IF(ALLOCATED(PEM%pNext)) THEN + PEM%pNext(NewID)=PEM%pNext(OldID) + ! Update pNext onto this particle + TempPartID = PEM%pStart(PEM%LocalElemID(OldID)) + IF (TempPartID.EQ.OldID) THEN + PEM%pStart(PEM%LocalElemID(OldID)) = NewID + ELSE + DO i=1,PEM%pNumber(PEM%LocalElemID(OldID)) + IF(PEM%pNext(TempPartID).EQ.OldID) THEN + PEM%pNext(TempPartID) = NewID + EXIT + END IF + TempPartID = PEM%pNext(TempPartID) + END DO + END IF +END IF +IF(ALLOCATED(PEM%LastGlobalElemID)) PEM%LastGlobalElemID(NewID)=PEM%LastGlobalElemID(OldID) + +IF(ALLOCATED(PDM%ParticleInside)) THEN + PDM%ParticleInside(NewID)=PDM%ParticleInside(OldID) + PDM%ParticleInside(OldID)=.FALSE. +END IF +IF(ALLOCATED(PDM%IsNewPart)) THEN + PDM%IsNewPart(NewID)=PDM%IsNewPart(OldID) + PDM%IsNewPart(OldID)=.FALSE. +END IF +IF(ALLOCATED(PDM%dtFracPush)) THEN + PDM%dtFracPush(NewID)=PDM%dtFracPush(OldID) + PDM%dtFracPush(OldID)=.FALSE. +END IF +IF(ALLOCATED(PDM%InRotRefFrame)) THEN + PDM%InRotRefFrame(NewID)=PDM%InRotRefFrame(OldID) + PDM%InRotRefFrame(OldID)=.FALSE. +END IF +IF(ALLOCATED(PDM%PartInit)) PDM%PartInit(NewID)=PDM%PartInit(OldID) + +IF(ALLOCATED(PartState)) THEN + PartState(:,NewID)=PartState(:,OldID) + PartState(:,OldID) = 0.0 +END IF +IF(ALLOCATED(LastPartPos)) LastPartPos(:,NewID)=LastPartPos(:,OldID) +IF(ALLOCATED(PartPosRef)) THEN + PartPosRef(:,NewID)=PartPosRef(:,OldID) + PartPosRef(:,OldID) = -888. +END IF +IF(ALLOCATED(PartSpecies)) THEN + PartSpecies(NewID)=PartSpecies(OldID) + PartSpecies(OldID) = 0 +END IF +IF(ALLOCATED(PartTimeStep)) PartTimeStep(NewID)=PartTimeStep(OldID) +IF(ALLOCATED(PartMPF)) PartMPF(NewID)=PartMPF(OldID) +IF(ALLOCATED(PartVeloRotRef)) THEN + PartVeloRotRef(:,NewID)=PartVeloRotRef(:,OldID) + PartVeloRotRef(:,OldID) = 0.0 +END IF +IF(ALLOCATED(PartStateIntEn)) PartStateIntEn(:,NewID)=PartStateIntEn(:,OldID) + +IF(ALLOCATED(Pt_temp)) THEN + Pt_temp(:,NewID)=Pt_temp(:,OldID) + Pt_temp(:,OldID) = 0.0 +END IF +IF(ALLOCATED(Pt)) THEN + Pt(:,NewID)=Pt(:,OldID) + Pt(:,OldID) = 0 +END IF +IF(ALLOCATED(FieldAtParticle)) FieldAtParticle(:,NewID)=FieldAtParticle(:,OldID) + +IF(ALLOCATED(InterPlanePartIndx)) InterPlanePartIndx(NewID)=InterPlanePartIndx(OldID) +IF(ALLOCATED(BGGas%PairingPartner)) BGGas%PairingPartner(NewID)=BGGas%PairingPartner(OldID) +IF(ALLOCATED(CollInf%OldCollPartner)) THEN + CollInf%OldCollPartner(NewID)=CollInf%OldCollPartner(OldID) + IF(CollInf%OldCollPartner(NewID).GT.0.AND.CollInf%OldCollPartner(NewID).LE.PDM%maxParticleNumber) THEN + IF(CollInf%OldCollPartner(CollInf%OldCollPartner(NewID)).EQ.OldID) CollInf%OldCollPartner(CollInf%OldCollPartner(NewID))=NewID + END IF +END IF +IF(ALLOCATED(ElecRelaxPart)) ElecRelaxPart(NewID)=ElecRelaxPart(OldID) + +#if (PP_TimeDiscMethod==508) || (PP_TimeDiscMethod==509) +IF(ALLOCATED(velocityAtTime)) velocityAtTime(:,NewID)=velocityAtTime(:,OldID) +#endif + +#if USE_MPI +IF(ALLOCATED(PartTargetProc)) PartTargetProc(NewID)=PartTargetProc(OldID) +IF(ALLOCATED(PartShiftVector)) PartShiftVector(:,NewID)=PartShiftVector(:,OldID) +#endif + +#if defined(IMPA) || defined(ROS) +IF(ALLOCATED(PartXK)) PartXK(:,NewID)=PartXK(:,OldID) +IF(ALLOCATED(R_PartXK)) R_PartXK(:,NewID)=R_PartXK(:,OldID) +IF(ALLOCATED(PartStage)) PartStage(:,:,NewID)=PartStage(:,:,OldID) +IF(ALLOCATED(PartStateN)) PartStateN(:,NewID)=PartStateN(:,OldID) +IF(ALLOCATED(PartQ)) PartQ(:,NewID)=PartQ(:,OldID) +IF(ALLOCATED(PEM%NormVec)) THEN + PEM%NormVec(:,NewID)=PEM%NormVec(:,OldID) + PEM%NormVec(:,OldID) = 0. +END IF +IF(ALLOCATED(PEM%PeriodicMoved)) THEN + PEM%PeriodicMoved(NewID)=PEM%PeriodicMoved(OldID) + PEM%PeriodicMoved(OldID) = .FALSE. +END IF +IF(ALLOCATED(PartDtFrac)) THEN + PartDtFrac(NewID)=PartDtFrac(OldID) + PartDtFrac(OldID) = 1. +END IF +#endif + +#ifdef IMPA +IF(ALLOCATED(F_PartX0)) F_PartX0(:,NewID)=F_PartX0(:,OldID) +IF(ALLOCATED(F_PartXk)) F_PartXk(:,NewID)=F_PartXk(:,OldID) +IF(ALLOCATED(Norm_F_PartX0)) Norm_F_PartX0(NewID)=Norm_F_PartX0(OldID) +IF(ALLOCATED(Norm_F_PartXk)) Norm_F_PartXk(NewID)=Norm_F_PartXk(OldID) +IF(ALLOCATED(Norm_F_PartXk_old)) Norm_F_PartXk_old(NewID)=Norm_F_PartXk_old(OldID) +IF(ALLOCATED(PartDeltaX)) PartDeltaX(:,NewID)=PartDeltaX(:,OldID) +IF(ALLOCATED(PartLambdaAccept)) PartLambdaAccept(NewID)=PartLambdaAccept(OldID) +IF(ALLOCATED(DoPartInNewton)) DoPartInNewton(NewID)=DoPartInNewton(OldID) +IF(ALLOCATED(PartIsImplicit)) PartIsImplicit(NewID)=PartIsImplicit(OldID) +#endif /* IMPA */ + + +! __ __ __ __ ________ ______ ___________ ___ +! / / / /___ ____/ /___ _/ /____ /_ __/\ \/ / __ \/ ____/ ___/ / | ______________ ___ _______ +! / / / / __ \/ __ / __ `/ __/ _ \ / / \ / /_/ / __/ \__ \ / /| | / ___/ ___/ __ `/ / / / ___/ +! / /_/ / /_/ / /_/ / /_/ / /_/ __/ / / / / ____/ /___ ___/ / / ___ |/ / / / / /_/ / /_/ (__ ) +! \____/ .___/\__,_/\__,_/\__/\___/ /_/ /_/_/ /_____//____/ /_/ |_/_/ /_/ \__,_/\__, /____/ +! /_/ /____/ + +IF(ALLOCATED(AmbipolElecVelo)) THEN + IF(ALLOCATED(AmbipolElecVelo(OldID)%ElecVelo)) THEN + IF(ALLOCATED(AmbipolElecVelo(NewID)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(NewID)%ElecVelo) + CALL MOVE_ALLOC(AmbipolElecVelo(OldID)%ElecVelo,AmbipolElecVelo(NewID)%ElecVelo) + ELSE + IF(ALLOCATED(AmbipolElecVelo(NewID)%ElecVelo)) DEALLOCATE(AmbipolElecVelo(NewID)%ElecVelo) + END IF +END IF + +IF(ALLOCATED(ElectronicDistriPart)) THEN + IF(ALLOCATED(ElectronicDistriPart(OldID)%DistriFunc)) THEN + IF(ALLOCATED(ElectronicDistriPart(NewID)%DistriFunc)) DEALLOCATE(ElectronicDistriPart(NewID)%DistriFunc) + CALL MOVE_ALLOC(ElectronicDistriPart(OldID)%DistriFunc,ElectronicDistriPart(NewID)%DistriFunc) + ELSE + IF(ALLOCATED(ElectronicDistriPart(NewID)%DistriFunc)) DEALLOCATE(ElectronicDistriPart(NewID)%DistriFunc) + END IF +END IF + +IF(ALLOCATED(VibQuantsPar)) THEN + IF(ALLOCATED(VibQuantsPar(OldID)%Quants)) THEN + IF(ALLOCATED(VibQuantsPar(NewID)%Quants)) DEALLOCATE(VibQuantsPar(NewID)%Quants) + CALL MOVE_ALLOC(VibQuantsPar(OldID)%Quants,VibQuantsPar(NewID)%Quants) + ELSE + IF(ALLOCATED(VibQuantsPar(NewID)%Quants)) DEALLOCATE(VibQuantsPar(NewID)%Quants) + END IF +END IF + +END SUBROUTINE ChangePartID + + + END MODULE MOD_part_tools diff --git a/src/particles/particle_vMPF.f90 b/src/particles/particle_vMPF.f90 index 269ff71f6..5cbd54301 100644 --- a/src/particles/particle_vMPF.f90 +++ b/src/particles/particle_vMPF.f90 @@ -274,7 +274,7 @@ SUBROUTINE MergeParticles(iPartIndx_Node_in, nPart, nPartNew, iElem) END IF END IF DOF_rot = SpecDSMC(iSpec)%Xi_Rot - T_rot = 2.*E_rot/(DOF_rot*totalWeight*BoltzmannConst) + T_rot = 2.*E_rot/(DOF_rot*totalWeight*BoltzmannConst) END IF IF(DSMC%ElectronicModel.GT.0.AND.SpecDSMC(iSpec)%InterID.NE.4) THEN T_elec = CalcTelec(E_elec/totalWeight, iSpec) @@ -332,7 +332,7 @@ SUBROUTINE MergeParticles(iPartIndx_Node_in, nPart, nPartNew, iElem) IF(DSMC%ElectronicModel.GT.0.AND.SpecDSMC(iSpec)%InterID.NE.4) THEN Energy_Sum = E_elec IF (E_elec.GT.0.0) THEN - IF (E_elec_new.EQ.0.0) THEN + IF (E_elec_new.EQ.0.0) THEN ! E_elec_new = 0.0 DO iLoop = 1, nPartNew ! temporal continuous energy distribution iPart = iPartIndx_Node(iLoop) @@ -342,7 +342,7 @@ SUBROUTINE MergeParticles(iPartIndx_Node_in, nPart, nPartNew, iElem) partWeight = GetParticleWeight(iPart) E_elec_new = E_elec_new + partWeight * PartStateIntEn(3,iPart) END DO - END IF + END IF alpha = E_elec/E_elec_new DO iLoop = 1, nPartNew ! alpha = E_elec/E_elec_new @@ -380,7 +380,7 @@ SUBROUTINE MergeParticles(iPartIndx_Node_in, nPart, nPartNew, iElem) Energy_Sum = Energy_Sum + E_vib IF (E_vib.GT.0.0) THEN IF (E_vib_new.EQ.0.0) THEN - IF(SpecDSMC(iSpec)%PolyatomicMol) THEN + IF(SpecDSMC(iSpec)%PolyatomicMol) THEN DO iLoop = 1, nPartNew ! temporal continuous energy distribution iPart = iPartIndx_Node(iLoop) PartStateIntEn(1,iPart) = 0.0 @@ -411,7 +411,7 @@ SUBROUTINE MergeParticles(iPartIndx_Node_in, nPart, nPartNew, iElem) END DO END DO END IF - END IF + END IF alpha = E_vib/E_vib_new DO iLoop = 1, nPartNew iPart = iPartIndx_Node(iLoop) @@ -454,7 +454,7 @@ SUBROUTINE MergeParticles(iPartIndx_Node_in, nPart, nPartNew, iElem) Energy_Sum = 0.0 ! 6.3) ensuring rotational excitation - alpha = E_rot/E_rot_new + alpha = E_rot/E_rot_new DO iLoop = 1, nPartNew iPart = iPartIndx_Node(iLoop) partWeight = GetParticleWeight(iPart) @@ -567,6 +567,7 @@ SUBROUTINE SplitParticles(iPartIndx_Node, nPartIn, nPartNew) USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, PartTimeStep USE MOD_DSMC_Vars ,ONLY: PartStateIntEn, CollisMode, SpecDSMC, DSMC, PolyatomMolDSMC, VibQuantsPar USE MOD_Particle_Tracking_Vars,ONLY: TrackingMethod +USE MOD_Part_Tools ,ONLY: GetNextFreePosition !#ifdef CODE_ANALYZE !USE MOD_Globals ,ONLY: unit_stdout,myrank,abort !USE MOD_Particle_Vars ,ONLY: Symmetry @@ -588,13 +589,13 @@ SUBROUTINE SplitParticles(iPartIndx_Node, nPartIn, nPartNew) iNewPart = 0 nPart = nPartIn nSplit = nPartNew - nPart -DO WHILE(iNewPart.LT.nSplit) +DO iNewPart=1,nSplit ! Get a random particle (only from the initially available) CALL RANDOM_NUMBER(iRan) iPart = INT(iRan*nPart) + 1 PartIndx = iPartIndx_Node(iPart) ! Check whether the weighting factor would drop below vMPFSplitLimit (can be below one) - ! If the resulting MPF is less than 1 you go sub-atomic where the concept of time and space become irrelevant... + ! If the resulting MPF is less than 1 you go sub-atomic where the concept of time and space become irrelevant... IF((PartMPF(PartIndx) / 2.).LT.vMPFSplitLimit) THEN ! Skip this particle iPartIndx_Node(iPart) = iPartIndx_Node(nPart) @@ -607,9 +608,7 @@ SUBROUTINE SplitParticles(iPartIndx_Node, nPartIn, nPartNew) IF((PartMPF(PartIndx) / 2.).LT.vMPFSplitLimit) CALL abort(__STAMP__,'Particle split below limit: PartMPF(PartIndx) / 2.=',& RealInfoOpt=PartMPF(PartIndx) / 2.) PartMPF(PartIndx) = PartMPF(PartIndx) / 2. ! split particle - iNewPart = iNewPart + 1 - PositionNbr = PDM%nextFreePosition(iNewPart+PDM%CurrentNextFreePosition) - IF (PositionNbr.EQ.0) CALL Abort(__STAMP__,'ERROR in particle split: MaxParticleNumber reached!') + PositionNbr = GetNextFreePosition() PartState(1:6,PositionNbr) = PartState(1:6,PartIndx) IF(TrackingMethod.EQ.REFMAPPING) PartPosRef(1:3,PositionNbr)=PartPosRef(1:3,PartIndx) PartSpecies(PositionNbr) = PartSpecies(PartIndx) @@ -635,9 +634,6 @@ SUBROUTINE SplitParticles(iPartIndx_Node, nPartIn, nPartNew) PEM%pNumber(LocalElemID) = PEM%pNumber(LocalElemID) + 1 END DO -! Advance particle vector length and the current next free position with newly created particles -PDM%ParticleVecLength = PDM%ParticleVecLength + iNewPart -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + iNewPart END SUBROUTINE SplitParticles diff --git a/src/particles/particle_vars.f90 b/src/particles/particle_vars.f90 index 313c01b56..5ac4b4f13 100644 --- a/src/particles/particle_vars.f90 +++ b/src/particles/particle_vars.f90 @@ -195,8 +195,14 @@ PPURE INTEGER FUNCTION ElemID_INTERFACE(iPart) TYPE tParticleDataManagement INTEGER :: CurrentNextFreePosition ! Index of nextfree index in nextFreePosition-Array INTEGER :: maxParticleNumber ! Maximum Number of all Particles + INTEGER :: maxAllowedParticleNumber ! Maximum allowed number of PDM%maxParticleNumber + LOGICAL :: RearrangePartIDs ! Rearrange PartIDs during shrinking maxPartNum +#if USE_MPI + LOGICAL :: UNFPafterMPIPartSend ! UpdateNextFreePosition after MPI Part Send +#endif INTEGER :: ParticleVecLength ! Vector Length for Particle Push Calculation INTEGER :: ParticleVecLengthOld ! Vector Length for Particle Push Calculation + REAL :: MaxPartNumIncrease ! How much shall the PDM%MaxParticleNumber be incresed if it is full INTEGER , ALLOCATABLE :: PartInit(:) ! (1:NParts), initial emission condition number ! the calculation area INTEGER ,ALLOCATABLE :: nextFreePosition(:) ! =>NULL() ! next_free_Position(1:maxParticleNumber) @@ -219,9 +225,9 @@ PPURE INTEGER FUNCTION ElemID_INTERFACE(iPart) INTEGER :: MacroValSamplIterNum ! Number of iterations for sampling ! macroscopic values LOGICAL :: SampleElecExcitation ! Sampling the electronic excitation rate per species -INTEGER :: ExcitationLevelCounter ! -REAL, ALLOCATABLE :: ExcitationSampleData(:,:) ! -INTEGER, ALLOCATABLE :: ExcitationLevelMapping(:,:) ! +INTEGER :: ExcitationLevelCounter ! +REAL, ALLOCATABLE :: ExcitationSampleData(:,:) ! +INTEGER, ALLOCATABLE :: ExcitationLevelMapping(:,:) ! INTEGER, ALLOCATABLE :: vMPFMergeThreshold(:) ! Max particle number per cell and (iSpec) INTEGER, ALLOCATABLE :: vMPFSplitThreshold(:) ! Min particle number per cell and (iSpec) @@ -284,7 +290,7 @@ PPURE INTEGER FUNCTION ElemID_INTERFACE(iPart) REAL :: RotRefFrameFreq ! frequency of rotational frame of reference REAL :: RotRefFrameOmega(3) ! angular velocity of rotational frame of reference INTEGER :: nRefFrameRegions ! number of rotational frame of reference regions -REAL, ALLOCATABLE :: RotRefFramRegion(:,:) ! MIN/MAX defintion for multiple rotational frame of reference region +REAL, ALLOCATABLE :: RotRefFramRegion(:,:) ! MIN/MAX defintion for multiple rotational frame of reference region ! (i,RegionNumber), MIN:i=1, MAX:i=2 !=================================================================================================================================== END MODULE MOD_Particle_Vars diff --git a/src/particles/pic/deposition/pic_depo_tools.f90 b/src/particles/pic/deposition/pic_depo_tools.f90 index 035664778..72051ca13 100644 --- a/src/particles/pic/deposition/pic_depo_tools.f90 +++ b/src/particles/pic/deposition/pic_depo_tools.f90 @@ -53,7 +53,8 @@ SUBROUTINE DepositPhotonSEEHoles(iBC,NbrOfParticle) USE MOD_Particle_Boundary_Vars ,ONLY: PartBound USE MOD_PICDepo_Vars ,ONLY: DoDeposition USE MOD_Dielectric_Vars ,ONLY: DoDielectricSurfaceCharge -USE MOD_Particle_Vars ,ONLY: PEM, PDM, PartSpecies, PartState, Species, usevMPF, PartMPF +USE MOD_Particle_Vars ,ONLY: PEM, PartSpecies, PartState, Species, usevMPF, PartMPF +USE MOD_Part_Tools ,ONLY: GetNextFreePosition IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------! ! INPUT / OUTPUT VARIABLES @@ -72,7 +73,7 @@ SUBROUTINE DepositPhotonSEEHoles(iBC,NbrOfParticle) IF(DoDeposition.AND.DoDielectricSurfaceCharge.AND.PartBound%Dielectric(iBC))THEN DO iPart = 1, NbrOfParticle ! Get index from next free position array - ParticleIndex = PDM%nextFreePosition(iPart+PDM%CurrentNextFreePosition) + ParticleIndex = GetNextFreePosition(iPart) ! Get charge IF(usevMPF)THEN diff --git a/src/particles/pic/models/pic_models.f90 b/src/particles/pic/models/pic_models.f90 index 69954f8f2..9b86ca374 100644 --- a/src/particles/pic/models/pic_models.f90 +++ b/src/particles/pic/models/pic_models.f90 @@ -70,6 +70,7 @@ SUBROUTINE ADK_Bruhwiler2003() USE MOD_Particle_Vars ,ONLY: PDM, Species, PartSpecies, usevMPF, PartState, PEM, PartMPF USE MOD_DSMC_Vars ,ONLY: DSMC, SpecDSMC USE MOD_PICInterpolation_Vars ,ONLY: FieldAtParticle +USE MOD_Part_Tools ,ONLY: GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -78,7 +79,7 @@ SUBROUTINE ADK_Bruhwiler2003() ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: iPart, MaxElecQua, ChargedNum, SumOfFormedParticles, ElectronIndex +INTEGER :: iPart, MaxElecQua, ChargedNum, ElectronIndex REAL :: IonizationEnergy_eV, iRan, QuantumTunnelProb, EffQuantNum REAL :: CriticalValue_GV REAL :: E_GV @@ -90,7 +91,6 @@ SUBROUTINE ADK_Bruhwiler2003() REAL(KIND=8) :: b(KK) = (/(ii, ii=1,KK, 1)/) #endif /* CODE_ANALYZE */ !=================================================================================================================================== -SumOfFormedParticles = 0 DO iPart = 1, PDM%ParticleVecLength IF(PDM%ParticleInside(iPart)) THEN @@ -139,12 +139,7 @@ SUBROUTINE ADK_Bruhwiler2003() CALL RANDOM_NUMBER(iRan) IF(QuantumTunnelProb.GT.iRan) THEN !.... Get free particle index for the 3rd particle produced - SumOfFormedParticles = SumOfFormedParticles + 1 - ElectronIndex = PDM%nextFreePosition(SumOfFormedParticles+PDM%CurrentNextFreePosition) - IF (ElectronIndex.EQ.0) THEN - CALL abort(__STAMP__,& - 'New Particle Number greater max Part Num in Field Ionization.') - END IF + ElectronIndex = GetNextFreePosition() !Set new Species of new particle PDM%ParticleInside(ElectronIndex) = .TRUE. PartSpecies(ElectronIndex) = DSMC%ElectronSpecies @@ -162,9 +157,6 @@ SUBROUTINE ADK_Bruhwiler2003() END IF END DO -PDM%ParticleVecLength = PDM%ParticleVecLength + SumOfFormedParticles -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + SumOfFormedParticles - END SUBROUTINE ADK_Bruhwiler2003 @@ -182,6 +174,7 @@ SUBROUTINE ADK_Yu2018() USE MOD_Particle_Vars ,ONLY: PDM, Species, PartSpecies, usevMPF, PartState, PEM, PartMPF USE MOD_DSMC_Vars ,ONLY: DSMC, SpecDSMC USE MOD_PICInterpolation_Vars ,ONLY: FieldAtParticle +USE MOD_Part_Tools ,ONLY: GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -190,7 +183,7 @@ SUBROUTINE ADK_Yu2018() ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: iPart, MaxElecQua, SumOfFormedParticles, ElectronIndex +INTEGER :: iPart, MaxElecQua, ElectronIndex REAL :: IonizationEnergy_eV, iRan, QuantumTunnelProb REAL :: n #ifdef CODE_ANALYZE @@ -202,7 +195,6 @@ SUBROUTINE ADK_Yu2018() REAL :: E #endif /* CODE_ANALYZE */ !=================================================================================================================================== -SumOfFormedParticles = 0 DO iPart = 1, PDM%ParticleVecLength IF(PDM%ParticleInside(iPart)) THEN @@ -240,8 +232,7 @@ SUBROUTINE ADK_Yu2018() CALL RANDOM_NUMBER(iRan) IF(QuantumTunnelProb.GT.iRan) THEN !.... Get free particle index for the 3rd particle produced - SumOfFormedParticles = SumOfFormedParticles + 1 - ElectronIndex = PDM%nextFreePosition(SumOfFormedParticles+PDM%CurrentNextFreePosition) + ElectronIndex = GetNextFreePosition() IF (ElectronIndex.EQ.0) THEN CALL abort(__STAMP__,& 'New Particle Number greater max Part Num in Field Ionization.') @@ -263,9 +254,6 @@ SUBROUTINE ADK_Yu2018() END IF END DO -PDM%ParticleVecLength = PDM%ParticleVecLength + SumOfFormedParticles -PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + SumOfFormedParticles - END SUBROUTINE ADK_Yu2018 diff --git a/src/particles/restart/particle_restart.f90 b/src/particles/restart/particle_restart.f90 index cb5922bd3..843aae76e 100644 --- a/src/particles/restart/particle_restart.f90 +++ b/src/particles/restart/particle_restart.f90 @@ -83,7 +83,7 @@ SUBROUTINE ParticleRestart() #endif /*USE_LOADBALANCE*/ ! Rotational frame of reference USE MOD_Particle_Vars ,ONLY: UseRotRefFrame, PartVeloRotRef, RotRefFrameOmega -USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck +USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck, IncreaseMaxParticleNumber ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -136,6 +136,7 @@ SUBROUTINE ParticleRestart() LastElemInd = offsetElem+PP_nElems locnPart = PartInt(ELEM_LastPartInd,LastElemInd)-PartInt(ELEM_FirstPartInd,FirstElemInd) offsetnPart = PartInt(ELEM_FirstPartInd,FirstElemInd) + CALL IncreaseMaxParticleNumber(INT(locnPart)) DO iLoop = 1_IK,locnPart ! Sanity check: SpecID > 0 @@ -278,17 +279,18 @@ SUBROUTINE ParticleRestart() SDEALLOCATE(MapPartDataToReadin) PDM%ParticleVecLength = PDM%ParticleVecLength + iPart +#ifdef CODE_ANALYZE + IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL Abort(__STAMP__,'PDM%ParticleVeclength exceeds PDM%maxParticleNumber, Difference:',IntInfoOpt=PDM%ParticleVeclength-PDM%maxParticleNumber) + DO iPart=PDM%ParticleVecLength+1,PDM%maxParticleNumber + IF (PDM%ParticleInside(iPart)) THEN + IPWRITE(*,*) iPart,PDM%ParticleVecLength,PDM%maxParticleNumber + CALL Abort(__STAMP__,'Particle outside PDM%ParticleVeclength',IntInfoOpt=iPart) + END IF + END DO +#endif CALL UpdateNextFreePosition() LBWRITE(UNIT_stdOut,*)' DONE!' - ! if ParticleVecLength GT maxParticleNumber: Stop - IF (PDM%ParticleVecLength.GT.PDM%maxParticleNumber) THEN - SWRITE (UNIT_stdOut,*) "PDM%ParticleVecLength =", PDM%ParticleVecLength - SWRITE (UNIT_stdOut,*) "PDM%maxParticleNumber =", PDM%maxParticleNumber - CALL abort(__STAMP__& - ,' Number of Particles in Restart file is higher than MaxParticleNumber! Increase MaxParticleNumber!') - END IF ! PDM%ParticleVecLength.GT.PDM%maxParticleNumber - ! Since the elementside-local node number are NOT persistant and dependent on the location ! of the MPI borders, all particle-element mappings need to be checked after a restart ! Step 1: Identify particles that are not in the element in which they were before the restart @@ -717,6 +719,16 @@ SUBROUTINE ParticleRestart() END DO ! iPart = 1, TotalNbrOfMissingParticlesSum PDM%ParticleVecLength = PDM%ParticleVecLength + NbrOfFoundParts +#ifdef CODE_ANALYZE + IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL Abort(__STAMP__,'PDM%ParticleVeclength exceeds PDM%maxParticleNumber, Difference:',IntInfoOpt=PDM%ParticleVeclength-PDM%maxParticleNumber) + DO iPart=PDM%ParticleVecLength+1,PDM%maxParticleNumber + IF (PDM%ParticleInside(iPart)) THEN + IPWRITE(*,*) iPart,PDM%ParticleVecLength,PDM%maxParticleNumber + CALL Abort(__STAMP__,'Particle outside PDM%ParticleVeclength',IntInfoOpt=iPart) + END IF + END DO +#endif + ! IF(PDM%ParticleVecLength.GT.PDM%maxParticleNumber) CALL IncreaseMaxParticleNumber(PDM%ParticleVecLength*CEILING(1+0.5*PDM%MaxPartNumIncrease)-PDM%maxParticleNumber) ! Combine number of found particles to make sure none are lost completely or found twice IF(MPIroot)THEN diff --git a/src/particles/ttm/ttm_init.f90 b/src/particles/ttm/ttm_init.f90 index 9d7c0f9dd..7d0a19518 100644 --- a/src/particles/ttm/ttm_init.f90 +++ b/src/particles/ttm/ttm_init.f90 @@ -771,6 +771,7 @@ SUBROUTINE InitIMD_TTM_Coupling() USE MOD_part_emission_tools ,ONLY: CalcVelocity_maxwell_lpn USE MOD_DSMC_Vars ,ONLY: useDSMC USE MOD_Eval_xyz ,ONLY: TensorProductInterpolation +USE MOD_Part_Tools ,ONLY: GetNextFreePosition !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT VARIABLES @@ -864,10 +865,7 @@ SUBROUTINE InitIMD_TTM_Coupling() DO iElem=1,PP_nElems DO iPart=1,ElemCharge(iElem) ! 1 electron for each charge of each element - ! Set the next free position in the particle vector list - PDM%CurrentNextFreePosition = PDM%CurrentNextFreePosition + 1 - ParticleIndexNbr = PDM%nextFreePosition(PDM%CurrentNextFreePosition) - PDM%ParticleVecLength = PDM%ParticleVecLength + 1 + ParticleIndexNbr = GetNextFreePosition() !Set new SpeciesID of new particle (electron) PDM%ParticleInside(ParticleIndexNbr) = .true. diff --git a/src/timedisc/timedisc.f90 b/src/timedisc/timedisc.f90 index 98d2238f3..baa572080 100644 --- a/src/timedisc/timedisc.f90 +++ b/src/timedisc/timedisc.f90 @@ -123,6 +123,7 @@ SUBROUTINE TimeDisc() #endif /*defined(MEASURE_MPI_WAIT)*/ #if defined(PARTICLES) USE MOD_Particle_Analyze_Vars ,ONLY: CalcPointsPerDebyeLength,CalcPICTimeStep +USE MOD_Part_Tools ,ONLY: ReduceMaxParticleNumber #endif ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -491,6 +492,7 @@ SUBROUTINE TimeDisc() IF(time.GE.tEnd)EXIT ! done, worst case: one additional time step #ifdef PARTICLES + CALL ReduceMaxParticleNumber() ! Switch flag to false after the number of particles has been written to std out and before the time next step is started GlobalNbrOfParticlesUpdated = .FALSE. #endif /*PARTICLES*/ diff --git a/src/timedisc/timedisc_TimeStep_BGK.f90 b/src/timedisc/timedisc_TimeStep_BGK.f90 index 5db0421aa..c1d10e934 100644 --- a/src/timedisc/timedisc_TimeStep_BGK.f90 +++ b/src/timedisc/timedisc_TimeStep_BGK.f90 @@ -151,10 +151,6 @@ SUBROUTINE TimeStep_BGK() (Time.ge.(1-DSMC%TimeFracSamp)*TEnd) .OR. & WriteMacroVolumeValues ) THEN CALL UpdateNextFreePosition(.TRUE.) !postpone UNFP for CollisMode=0 to next IterDisplayStep or when needed for DSMC-Sampling -ELSE IF (PDM%nextFreePosition(PDM%CurrentNextFreePosition+1).GT.PDM%maxParticleNumber .OR. & - PDM%nextFreePosition(PDM%CurrentNextFreePosition+1).EQ.0) THEN - ! gaps in PartState are not filled until next UNFP and array might overflow more easily! - CALL abort(__STAMP__,'maximum nbr of particles reached!') END IF #if USE_MPI diff --git a/src/timedisc/timedisc_TimeStep_DSMC.f90 b/src/timedisc/timedisc_TimeStep_DSMC.f90 index 48eb6ea79..0be6be9d0 100644 --- a/src/timedisc/timedisc_TimeStep_DSMC.f90 +++ b/src/timedisc/timedisc_TimeStep_DSMC.f90 @@ -207,10 +207,6 @@ SUBROUTINE TimeStep_DSMC() (Time.ge.(1-DSMC%TimeFracSamp)*TEnd) .OR. & WriteMacroVolumeValues.OR.WriteMacroSurfaceValues ) THEN CALL UpdateNextFreePosition() !postpone UNFP for CollisMode=0 to next IterDisplayStep or when needed for DSMC-Sampling -ELSE IF (PDM%nextFreePosition(PDM%CurrentNextFreePosition+1).GT.PDM%maxParticleNumber .OR. & - PDM%nextFreePosition(PDM%CurrentNextFreePosition+1).EQ.0) THEN - ! gaps in PartState are not filled until next UNFP and array might overflow more easily! - CALL abort(__STAMP__,'maximum nbr of particles reached!') END IF IF(DSMC%UseOctree)THEN diff --git a/src/timedisc/timedisc_TimeStep_FPFlow.f90 b/src/timedisc/timedisc_TimeStep_FPFlow.f90 index 315969555..9bc410564 100644 --- a/src/timedisc/timedisc_TimeStep_FPFlow.f90 +++ b/src/timedisc/timedisc_TimeStep_FPFlow.f90 @@ -124,11 +124,6 @@ SUBROUTINE TimeStep_FPFlow() (Time.ge.(1-DSMC%TimeFracSamp)*TEnd) .OR. & WriteMacroVolumeValues ) THEN CALL UpdateNextFreePosition() !postpone UNFP for CollisMode=0 to next IterDisplayStep or when needed for DSMC-Sampling -ELSE IF (PDM%nextFreePosition(PDM%CurrentNextFreePosition+1).GT.PDM%maxParticleNumber .OR. & - PDM%nextFreePosition(PDM%CurrentNextFreePosition+1).EQ.0) THEN - CALL abort(& -__STAMP__,& -'maximum nbr of particles reached!') !gaps in PartState are not filled until next UNFP and array might overflow more easily! END IF IF (CoupledFPDSMC) THEN