Merge branch '222-bugfix-vmpf-for-pure-pic' into 'master.dev'

Resolve "Bugfix vMPF for pure PIC"

Closes #222

See merge request piclas/piclas!869

docs/documentation/userguide/features-and-models/features-particle-solver.md

@@ -224,7 +224,7 @@ only at the stagnation point, the time step defined during the initialization is
 (sec:variable-particle-weighting)=
 #### Variable Particle Weighting
 
-Variable particle weighting is currently supported with PIC and/or with a background gas (an additional trace species feature is described in Section {ref}`sec:background-gas`). The general functionality can be enabled with the following flag:
+Variable particle weighting is currently supported for PIC (with and without background gas) or a background gas (an additional trace species feature is described in Section {ref}`sec:background-gas`). The general functionality can be enabled with the following flag:
 
     Part-vMPF                           = T
 

regressioncheck/NIG_LoadBalance/sphere_soft_RK4_without_DSMC/analyze.ini

@@ -2,7 +2,7 @@
 analyze_L2=1000
 
 ! check if particles are outside of domain at tEnd
-check_hdf5_file        = implicit_one_State_000.00000010000000000.h5
+check_hdf5_file        = implicit_one_State_000.00000002000000000.h5
 check_hdf5_data_set    = PartData
 check_hdf5_span        = 1 ! check all rows
 check_hdf5_dimension   = 0:2

regressioncheck/NIG_LoadBalance/sphere_soft_RK4_without_DSMC/parameter.ini

@@ -8,11 +8,9 @@ IniExactFunc  = 0
 ! =============================================================================== !
 N             = 3  ! Polynomial degree
 NAnalyze      = 2 ! Number of analyze points
-nSample = 4
 ! =============================================================================== !
 ! MESH
 ! =============================================================================== !
-!MeshFile      = Sphere_Ngeo4_004_005_mesh.h5
 MeshFile      = Sphere_Ngeo4_001_001_mesh.h5
 useCurveds    = T
 ! if boundaries have to be changed (else they are used from Mesh directly):
@@ -30,7 +28,7 @@ DoCalcErrorNorms = T
 ! Load Balance
 ! =============================================================================== !
 DoLoadBalance           = T
-Load-DeviationThreshold = 0.1
+Load-DeviationThreshold = 1E-9
 Particles-MPIWeight     = 0.01
 
 ! =============================================================================== !
@@ -39,15 +37,14 @@ Particles-MPIWeight     = 0.01
 CFLscale   = 0.9  ! Scaling of theoretical CFL number
 c_corr     = 1
 BezierClipTolerance = 1e-12
-!BezierNewtonTolerance = 1e-4
 
 ! =============================================================================== !
 ! IMPLICIT
 ! =============================================================================== !
-tend       = 1E-7  ! End time
+tend       = 2E-8  ! End time
 Analyze_dt = 1E-8 ! Timestep of analyze outputs
 CalcPotentialEnergy      = TRUE
-
+IterDisplayStep = 50
 ! =============================================================================== !
 ! PARTICLES
 ! =============================================================================== !
@@ -57,52 +54,33 @@ Part-Boundary1-Condition=reflective
 
 Part-FIBGMdeltas=(/.5,.5,.5/)
 
-Part-vMPF=F
 Part-maxParticleNumber=20000
 Part-nSpecies=1
-PIC-externalField=(/0.,0.,0.,0.,0.,0./)
 
 Part-Species1-ChargeIC=-1.6022E-19
 Part-Species1-MassIC=9.1093826E-31
 Part-Species1-MacroParticleFactor=1000
-Part-Species1-nInits=1
 
+Part-Species1-nInits=1
 Part-Species1-Init1-SpaceIC=cuboid
-Part-Species1-Init1-velocityDistribution=maxwell                                                                                           
+Part-Species1-Init1-velocityDistribution=maxwell
 Part-Species1-Init1-MWTemperatureIC=1e8
 Part-Species1-Init1-ParticleNumber=500
-
 Part-Species1-Init1-BasePointIC=(/.25,.25,-0.25/)
 Part-Species1-Init1-BaseVector1IC=(/-.5,0.0,0.0/)
 Part-Species1-Init1-BaseVector2IC=(/0.0,-.5,0.0/)
 Part-Species1-Init1-CuboidHeightIC=0.5
-
 Part-Species1-Init1-NormalIC=(/0.,0.,1./)
-
 Part-Species1-Init1-VeloIC=0.
 Part-Species1-Init1-VeloVecIC=(/1.,0.,0./)
 ! =============================================================================== !
 ! tracking
 ! =============================================================================== !
-RefMappingGuess=1 !,3
-!BezierClipTolerance=1e-8  
-!BezierClipMaxIter  =105
-!BezierClipHit      =2e-4
-!BezierSplitLimit   =0.6
-!!epsilontol         =1e-12
-!BezierElevation=20
-!RefMappingEps      =1e-8
-PIC-DoInterpolation=F
-
-
 TrackingMethod = refmapping
+RefMappingGuess=1
+PIC-DoInterpolation=F
 ! =============================================================================== !
 ! DSMC
 ! =============================================================================== !
 UseDSMC=F
-Particles-DSMCReservoirSim=false
-Particles-DSMC-CollisMode=0 ! Collisionless flow
-Part-NumberOfRandomSeeds =2
-Particles-RandomSeed1= 1
-Particles-RandomSeed2= 2
 Particles-HaloEpsVelo=50000

regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/2D_Landmark/parameter.ini

@@ -80,7 +80,6 @@ DisplayLostParticles=T
 
 Part-Species1-MacroParticleFactor  = 1.67e4 ! 1.67e2 originally used for z=1e-4 m (case2: 75 parts per cell with dx=dy=5e-5 m)
 Part-Species2-MacroParticleFactor  = 1.67e4 ! 1.67e2 originally used for z=1e-4 m (case2: 75 parts per cell with dx=dy=5e-5 m)
-Part-Species3-MacroParticleFactor  = 1.67e4 ! 1.67e2 originally used for z=1e-4 m (case2: 75 parts per cell with dx=dy=5e-5 m)
 
 ! =============================================================================== !
 ! Particle Boundary Conditions
@@ -123,6 +122,10 @@ Particles-RandomSeed2                   = 2
 
 Particles-HaloEpsVelo                   = 300e6
 
+Part-vMPF= F,T
+Part-Species1-vMPFMergeThreshold = 200
+Part-Species2-vMPFMergeThreshold = 200
+
 ! =============================================================================== !
 ! Species1 | e
 ! =============================================================================== !

src/particles/particle_init.f90

@@ -595,34 +595,21 @@ SUBROUTINE AllocateParticleArrays()
 PDM%nextFreePosition(1:PDM%maxParticleNumber)=0
 
 ALLOCATE(PEM%GlobalElemID(1:PDM%maxParticleNumber), STAT=ALLOCSTAT)
-IF (ALLOCSTAT.NE.0) CALL abort(&
-__STAMP__&
-  ,' Cannot allocate PEM%GlobalElemID(1:PDM%maxParticleNumber) array!')
+IF (ALLOCSTAT.NE.0) CALL abort(__STAMP__,' Cannot allocate PEM%GlobalElemID(1:PDM%maxParticleNumber) array!')
 
 ALLOCATE(PEM%LastGlobalElemID(1:PDM%maxParticleNumber), STAT=ALLOCSTAT)
-IF (ALLOCSTAT.NE.0) CALL abort(&
-__STAMP__&
-  ,' Cannot allocate PEM%LastGlobalElemID(1:PDM%maxParticleNumber) array!')
+IF (ALLOCSTAT.NE.0) CALL abort(__STAMP__,' Cannot allocate PEM%LastGlobalElemID(1:PDM%maxParticleNumber) array!')
 
-IF (useDSMC) THEN
+IF (useDSMC.OR.usevMPF) THEN
   ALLOCATE(PEM%pStart(1:nElems)                         , &
            PEM%pNumber(1:nElems)                        , &
            PEM%pEnd(1:nElems)                           , &
            PEM%pNext(1:PDM%maxParticleNumber)           , STAT=ALLOCSTAT)
-
-  IF (ALLOCSTAT.NE.0) THEN
-    CALL abort(&
-__STAMP__&
-    , ' Cannot allocate DSMC PEM arrays!')
-  END IF
+  IF (ALLOCSTAT.NE.0) CALL abort(__STAMP__, ' Cannot allocate DSMC PEM arrays!')
 END IF
 IF (useDSMC) THEN
   ALLOCATE(PDM%PartInit(1:PDM%maxParticleNumber), STAT=ALLOCSTAT)
-  IF (ALLOCSTAT.NE.0) THEN
-    CALL abort(&
-__STAMP__&
-    ,' Cannot allocate DSMC PEM arrays!')
-  END IF
+  IF (ALLOCSTAT.NE.0) CALL abort(__STAMP__,' Cannot allocate PDM%PartInit array!')
 END IF
 
 END SUBROUTINE AllocateParticleArrays

src/particles/particle_vMPF.f90

@@ -40,9 +40,12 @@ MODULE MOD_vMPF
 !===================================================================================================================================
 SUBROUTINE SplitAndMerge()
 ! MODULES
-USE MOD_PARTICLE_Vars ,ONLY: vMPFMergeThreshold, vMPFSplitThreshold, PEM, nSpecies, PartSpecies,PDM
-USE MOD_Mesh_Vars     ,ONLY: nElems
-USE MOD_part_tools    ,ONLY: UpdateNextFreePosition
+USE MOD_PARTICLE_Vars         ,ONLY: vMPFMergeThreshold, vMPFSplitThreshold, PEM, nSpecies, PartSpecies,PDM
+USE MOD_Mesh_Vars             ,ONLY: nElems
+USE MOD_part_tools            ,ONLY: UpdateNextFreePosition
+#if USE_LOADBALANCE
+USE MOD_LoadBalance_Timers    ,ONLY: LBStartTime, LBElemSplitTime
+#endif /*USE_LOADBALANCE*/
 ! IMPLICIT VARIABLE HANDLING
 IMPLICIT NONE
 !-----------------------------------------------------------------------------------------------------------------------------------
@@ -52,15 +55,21 @@ SUBROUTINE SplitAndMerge()
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! LOCAL VARIABLES
 INTEGER               :: iElem, iLoop, iPart, nPartCell, iSpec
-INTEGER, ALLOCATABLE  :: iPartIndx_Node(:), nPart(:),iPartIndx_Node_Temp(:,:)
+INTEGER, ALLOCATABLE  :: iPartIndx_Node(:,:), nPart(:)
+#if USE_LOADBALANCE
+REAL                  :: tLBStart
+#endif /*USE_LOADBALANCE*/
 !===================================================================================================================================
+#if USE_LOADBALANCE
+CALL LBStartTime(tLBStart)
+#endif /*USE_LOADBALANCE*/
 ALLOCATE(nPart(nSpecies))
 DO iElem = 1, nElems
   nPart(:) = 0
   nPartCell = PEM%pNumber(iElem)
-  ALLOCATE(iPartIndx_Node_Temp(nSpecies,nPartCell))
+  ALLOCATE(iPartIndx_Node(nSpecies,nPartCell))
   DO iSpec = 1, nSpecies
-    iPartIndx_Node_Temp(iSpec,1:nPartCell) = 0
+    iPartIndx_Node(iSpec,1:nPartCell) = 0
   END DO
   iPart = PEM%pStart(iElem)
 
@@ -72,29 +81,29 @@ SUBROUTINE SplitAndMerge()
     END IF
     iSpec = PartSpecies(iPart)
     nPart(iSpec) = nPart(iSpec) + 1
-    iPartIndx_Node_Temp(iSpec,nPart(iSpec)) = iPart
+    iPartIndx_Node(iSpec,nPart(iSpec)) = iPart
     iPart = PEM%pNext(iPart)
   END DO
 
   DO iSpec = 1, nSpecies
-    IF((vMPFMergeThreshold(iSpec).EQ.0).AND.(vMPFSplitThreshold(iSpec).EQ.0)) CYCLE            ! Skip default values
-    IF(nPart(iSpec).EQ.0) CYCLE                     ! Skip when no particles are present
+    ! Skip default values
+    IF((vMPFMergeThreshold(iSpec).EQ.0).AND.(vMPFSplitThreshold(iSpec).EQ.0)) CYCLE
+    ! Skip when no particles are present
+    IF(nPart(iSpec).EQ.0) CYCLE
 
-    ! 2.) build partindx list for species
-    ALLOCATE(iPartIndx_Node(nPart(iSpec)))
-    iPartIndx_Node(1:nPart(iSpec)) = iPartIndx_Node_Temp(iSpec,1:nPart(iSpec))
-
-    ! 3.) Call split or merge routine
+    ! 2.) Call split or merge routine
     IF(nPart(iSpec).GT.vMPFMergeThreshold(iSpec).AND.(vMPFMergeThreshold(iSpec).NE.0)) THEN   ! Merge
-      CALL MergeParticles(iPartIndx_Node, nPart(iSpec), vMPFMergeThreshold(iSpec),iElem)
+      CALL MergeParticles(iPartIndx_Node(iSpec,1:nPart(iSpec)), nPart(iSpec), vMPFMergeThreshold(iSpec),iElem)
     ELSE IF(nPart(iSpec).LT.vMPFSplitThreshold(iSpec)) THEN                                   ! Split
-      CALL SplitParticles(iPartIndx_Node, nPart(iSpec), vMPFSplitThreshold(iSpec))
+      CALL SplitParticles(iPartIndx_Node(iSpec,1:nPart(iSpec)), nPart(iSpec), vMPFSplitThreshold(iSpec))
     END IF
-    DEALLOCATE(iPartIndx_Node)
 
   END DO
-  DEALLOCATE(iPartIndx_Node_Temp)
+  DEALLOCATE(iPartIndx_Node)
 
+#if USE_LOADBALANCE
+  CALL LBElemSplitTime(iElem,tLBStart)
+#endif /*USE_LOADBALANCE*/
 END DO
 CALL UpdateNextFreePosition()
 DEALLOCATE(nPart)

src/timedisc/timedisc_TimeStepByLSERK.f90

@@ -290,18 +290,15 @@ SUBROUTINE TimeStepByLSERK()
 
 IF ((time.GE.DelayTime).OR.(time.EQ.0)) CALL UpdateNextFreePosition()
 
-IF (useDSMC) THEN
-  IF (time.GE.DelayTime) THEN
+IF (time.GE.DelayTime) THEN
+  ! Direct Simulation Monte Carlo
+  IF (useDSMC) THEN
     CALL DSMC_main()
-#if USE_LOADBALANCE
-    CALL LBStartTime(tLBStart)
-#endif /*USE_LOADBALANCE*/
-    IF(UseSplitAndMerge) CALL SplitAndMerge()
-#if USE_LOADBALANCE
-    CALL LBPauseTime(LB_DSMC,tLBStart)
-#endif /*USE_LOADBALANCE*/
   END IF
+  ! Split & Merge: Variable particle weighting
+  IF(UseSplitAndMerge) CALL SplitAndMerge()
 END IF
+
 #ifdef EXTRAE
 CALL extrae_eventandcounters(int(9000001), int8(0))
 #endif /*EXTRAE*/

src/timedisc/timedisc_TimeStepPoisson.f90

@@ -251,17 +251,13 @@ SUBROUTINE TimeStepPoisson()
 
 IF ((time.GE.DelayTime).OR.(iter.EQ.0)) CALL UpdateNextFreePosition()
 
-IF (useDSMC) THEN
-  IF (time.GE.DelayTime) THEN
+IF (time.GE.DelayTime) THEN
+  ! Direct Simulation Monte Carlo
+  IF (useDSMC) THEN
     CALL DSMC_main()
-#if USE_LOADBALANCE
-    CALL LBStartTime(tLBStart)
-#endif /*USE_LOADBALANCE*/
-    IF(UseSplitAndMerge) CALL SplitAndMerge()
-#if USE_LOADBALANCE
-    CALL LBPauseTime(LB_DSMC,tLBStart)
-#endif /*USE_LOADBALANCE*/
   END IF
+  ! Split & Merge: Variable particle weighting
+  IF(UseSplitAndMerge) CALL SplitAndMerge()
 END IF
 #endif /*PARTICLES*/
 

src/timedisc/timedisc_TimeStepPoissonByBorisLeapfrog.f90

@@ -263,18 +263,15 @@ SUBROUTINE TimeStepPoissonByBorisLeapfrog()
 
 IF ((time.GE.DelayTime).OR.(iter.EQ.0)) CALL UpdateNextFreePosition()
 
-IF (useDSMC) THEN
-  IF (time.GE.DelayTime) THEN
+IF (time.GE.DelayTime) THEN
+  ! Direct Simulation Monte Carlo
+  IF (useDSMC) THEN
     CALL DSMC_main()
-#if USE_LOADBALANCE
-    CALL LBStartTime(tLBStart)
-#endif /*USE_LOADBALANCE*/
-    IF(UseSplitAndMerge) CALL SplitAndMerge()
-#if USE_LOADBALANCE
-    CALL LBPauseTime(LB_DSMC,tLBStart)
-#endif /*USE_LOADBALANCE*/
   END IF
+  ! Split & Merge: Variable particle weighting
+  IF(UseSplitAndMerge) CALL SplitAndMerge()
 END IF
+
 #ifdef EXTRAE
 CALL extrae_eventandcounters(int(9000001), int8(0))
 #endif /*EXTRAE*/

src/timedisc/timedisc_TimeStepPoissonByHigueraCary.f90

@@ -251,18 +251,15 @@ SUBROUTINE TimeStepPoissonByHigueraCary()
 
 IF ((time.GE.DelayTime).OR.(iter.EQ.0)) CALL UpdateNextFreePosition()
 
-IF (useDSMC) THEN
-  IF (time.GE.DelayTime) THEN
+IF (time.GE.DelayTime) THEN
+  ! Direct Simulation Monte Carlo
+  IF (useDSMC) THEN
     CALL DSMC_main()
-#if USE_LOADBALANCE
-    CALL LBStartTime(tLBStart)
-#endif /*USE_LOADBALANCE*/
-    IF(UseSplitAndMerge) CALL SplitAndMerge()
-#if USE_LOADBALANCE
-    CALL LBPauseTime(LB_DSMC,tLBStart)
-#endif /*USE_LOADBALANCE*/
   END IF
+  ! Split & Merge: Variable particle weighting
+  IF(UseSplitAndMerge) CALL SplitAndMerge()
 END IF
+
 #ifdef EXTRAE
 CALL extrae_eventandcounters(int(9000001), int8(0))
 #endif /*EXTRAE*/

src/timedisc/timedisc_TimeStepPoissonByLSERK.f90

@@ -369,17 +369,13 @@ SUBROUTINE TimeStepPoissonByLSERK()
 #ifdef PARTICLES
 IF ((time.GE.DelayTime).OR.(iter.EQ.0)) CALL UpdateNextFreePosition()
 
-IF (useDSMC) THEN
-  IF (time.GE.DelayTime) THEN
+IF (time.GE.DelayTime) THEN
+  ! Direct Simulation Monte Carlo
+  IF (useDSMC) THEN
     CALL DSMC_main()
-#if USE_LOADBALANCE
-    CALL LBStartTime(tLBStart)
-#endif /*USE_LOADBALANCE*/
-    IF(UseSplitAndMerge) CALL SplitAndMerge()
-#if USE_LOADBALANCE
-    CALL LBPauseTime(LB_DSMC,tLBStart)
-#endif /*USE_LOADBALANCE*/
   END IF
+  ! Split & Merge: Variable particle weighting
+  IF(UseSplitAndMerge) CALL SplitAndMerge()
 END IF
 #endif /*PARTICLES*/
 

src/timedisc/timedisc_TimeStep_DSMC.f90

@@ -237,16 +237,9 @@ SUBROUTINE TimeStep_DSMC()
 
 CALL DSMC_main()
 
-#if USE_LOADBALANCE
-CALL LBStartTime(tLBStart)
-#endif /*USE_LOADBALANCE*/
-
+! Split & Merge: Variable particle weighting
 IF(UseSplitAndMerge) CALL SplitAndMerge()
 
-#if USE_LOADBALANCE
-CALL LBPauseTime(LB_DSMC,tLBStart)
-#endif /*USE_LOADBALANCE*/
-
 END SUBROUTINE TimeStep_DSMC