Merge branch 'fix.changeArraySize.default' into 'master.dev'

[fix.changeArraySize.default] Initialize arrays during ChangeSizeArray

Closes #265

See merge request piclas/piclas!939

CONTRIBUTORS.md

@@ -1,6 +1,6 @@
 # List of Contributors
 
-This is a (possibly incomplete) list of the people who contributed to piclas.
+This is a (possibly incomplete) list of the people who contributed to PICLas.
 
 * Christoph Altmann
 * Muhammed Atak
@@ -15,29 +15,36 @@ This is a (possibly incomplete) list of the people who contributed to piclas.
 * David Flad
 * Fabian Föll
 * Hannes Frank
+* Félix Garmirian
 * Gregor Gassner
 * Florian Hindenlang
 * Andrea Hinkel
 * Konstantin Hinsberger
+* Franziska Hild
 * Timon Hitz
 * Malte Hoffmann
 * Serena Keller
 * Johannes Kleinert
 * Marius Koch
 * Patrick Kopper
 * Nico Krais
+* Simone Lauterbach
 * Robert Mansk
+* Claudia Marianowski
 * Claus-Dieter Munz
 * Asim Mirza
 * Jonathan Neudorfer
 * Paul Nizenkov
 * Jonatan Núñez-de la Rosa
 * Philipp Offenhaeuser
 * Philip Ortwein
+* Tobias Ott
 * Marcel Pfeiffer
 * Wladimir Reschke
 * Matthias Sonntag
 * Marcel "Mitch" Scherrmann
+* Valentin Schunck
+* Miklas Schütte
 * Anna Schwarz
 * Marc Staudenmaier
 * Torsten Stindl

src/particles/boundary/particle_boundary_init.f90

@@ -1952,6 +1952,7 @@ SUBROUTINE InitRotPeriodicInterPlaneMapping()
 
 ALLOCATE(InterPlanePartIndx(1:PDM%maxParticleNumber), STAT=ALLOCSTAT)
 IF (ALLOCSTAT.NE.0) CALL abort(__STAMP__,'ERROR in particle_boundary_init.f90: Cannot allocate InterPlanePartIndx array!')
+InterPlanePartIndx = 0
 
 HasInterPlaneOnProc = .FALSE.
 

src/particles/dsmc/dsmc_chemical_init.f90

@@ -1205,10 +1205,10 @@ SUBROUTINE ReadReacFromDatabase(ReadInNumOfReact)
     DO iSpec = 1, 3
       IF(ChemReac%Reactants(iReac,iSpec).NE.0) SumReactMass = SumReactMass + Species(ChemReac%Reactants(iReac,iSpec))%MassIC
     END DO
-    ! Santiy mass check for reactions of selected chemistry model, real compare with RelMassTol
+    ! Sanity mass check for reactions of selected chemistry model, real compare with RelMassTol
     IF(.NOT.ALMOSTEQUALRELATIVE(SumProdMass,SumReactMass,RelMassTol)) THEN
-      CALL PrintOption('DSMC_Chemistry might not mass conserving for chemical reaction:','WARNING',StrOpt=TRIM(ChemReac%ReactionName(iReac)))
-      CALL abort(__STAMP__,'DSMC_Chemistry might not mass conserving for current chemical reaction')
+      CALL PrintOption('DSMC_Chemistry might not be mass conserving for chemical reaction:','WARNING',StrOpt=TRIM(ChemReac%ReactionName(iReac)))
+      CALL abort(__STAMP__,'DSMC_Chemistry might not be mass conserving for current chemical reaction!')
     END IF
       ! Read-in of the reaction parameters, depending on the model
     SELECT CASE (TRIM(ChemReac%ReactModel(iReac)))

src/particles/mcc/mcc.f90

@@ -612,12 +612,16 @@ END SUBROUTINE MonteCarloCollision
 
 
 !===================================================================================================================================
-!> Calculate the reaction probability if collision cross-section data is used (only with a background gas from the MCC routine)
+!> Calculate the reaction probability if cross-section data is used to be added to the DSMC-style collision probability
+!> (only with a background gas from the MCC routine)
+!> NOTE: Relativistic collision energy is only utilized to determine the cross-section, everything else is done using the classical
+!>       approach
 !===================================================================================================================================
 SUBROUTINE MCC_CalcReactionProb(iCase,bgSpec,CRela2,CollEnergy_in,PartIndex,bggPartIndex,iElem)
 ! MODULES
+USE MOD_Globals_Vars          ,ONLY: RelativisticLimit
 USE MOD_Particle_Vars         ,ONLY: Species, PartSpecies, VarTimeStep
-USE MOD_DSMC_Vars             ,ONLY: SpecDSMC, BGGas, ChemReac, DSMC, PartStateIntEn
+USE MOD_DSMC_Vars             ,ONLY: SpecDSMC, BGGas, ChemReac, DSMC, PartStateIntEn, CollInf
 USE MOD_MCC_Vars              ,ONLY: SpecXSec
 USE MOD_Particle_Vars         ,ONLY: Species
 USE MOD_TimeDisc_Vars         ,ONLY: dt
@@ -632,7 +636,7 @@ SUBROUTINE MCC_CalcReactionProb(iCase,bgSpec,CRela2,CollEnergy_in,PartIndex,bggP
 ! LOCAL VARIABLES
 INTEGER                       :: jSpec, iPath, ReacTest, EductReac(1:3), ProductReac(1:4), iProd
 INTEGER                       :: NumWeightProd
-REAL                          :: EZeroPoint_Educt, EZeroPoint_Prod, CollEnergy
+REAL                          :: EZeroPoint_Educt, EZeroPoint_Prod, CollEnergy, CollEnergyNonRela
 REAL                          :: CrossSection, dtVar
 REAL                          :: Temp_Rot, Temp_Vib, Temp_Elec, BGGasNumDens, BGGasFraction
 !===================================================================================================================================
@@ -673,8 +677,8 @@ SUBROUTINE MCC_CalcReactionProb(iCase,bgSpec,CRela2,CollEnergy_in,PartIndex,bggP
         EZeroPoint_Prod = EZeroPoint_Prod + SpecDSMC(ProductReac(iProd))%EZeroPoint
       END IF
     END DO
-    ! Adding the internal energy of particle species
-    CollEnergy = CollEnergy_in + PartStateIntEn(1,PartIndex) + PartStateIntEn(2,PartIndex)
+    ! Adding the internal energy of particle species (relative translational energy is added at the end)
+    CollEnergy = PartStateIntEn(1,PartIndex) + PartStateIntEn(2,PartIndex)
     ! Internal energy of background species
     IF((Species(jSpec)%InterID.EQ.2).OR.(Species(jSpec)%InterID.EQ.20)) THEN
       IF(BGGas%UseDistribution) THEN
@@ -697,8 +701,19 @@ SUBROUTINE MCC_CalcReactionProb(iCase,bgSpec,CRela2,CollEnergy_in,PartIndex,bggP
       PartStateIntEn(3,bggPartIndex) = CalcEElec_particle(jSpec,Temp_Elec,bggPartIndex)
       CollEnergy = CollEnergy + PartStateIntEn(3,PartIndex) + PartStateIntEn(3,bggPartIndex)
     END IF
+    ! Work-around for relativistic energies: since the energy distribution after the reaction is not done relativistically yet,
+    ! we have to check whether sufficient collision energy is available in the classical manner
+    IF(CRela2 .LT. RelativisticLimit) THEN
+      ! Classical
+      CollEnergyNonRela = CollEnergy_in + CollEnergy
+      CollEnergy = CollEnergyNonRela
+    ELSE
+      ! Relativistic
+      CollEnergyNonRela = 0.5 * CollInf%MassRed(iCase) * CRela2 + CollEnergy
+      CollEnergy = CollEnergy_in + CollEnergy
+    END IF
     ! Check first if sufficient energy is available for the products after the reaction
-    IF(((CollEnergy-EZeroPoint_Prod).GE.-ChemReac%EForm(ReacTest))) THEN
+    IF(((CollEnergyNonRela-EZeroPoint_Prod).GE.-ChemReac%EForm(ReacTest))) THEN
       CollEnergy = CollEnergy - EZeroPoint_Educt
       CrossSection = InterpolateCrossSection(SpecXSec(iCase)%ReactionPath(iPath)%XSecData,CollEnergy)
       ASSOCIATE( ReactionProb => ChemReac%CollCaseInfo(iCase)%ReactionProb(iPath) )

src/particles/particle_init.f90

@@ -1027,6 +1027,7 @@ SUBROUTINE InitializeVariablesvMPF()
 
   ALLOCATE(PartMPF(1:PDM%maxParticleNumber), STAT=ALLOCSTAT)
   IF (ALLOCSTAT.NE.0) CALL abort(__STAMP__,'ERROR in particle_init.f90: Cannot allocate Particle arrays!')
+  PartMPF = 1.
 END IF
 END SUBROUTINE InitializeVariablesvMPF
 

src/particles/particle_tools.f90

@@ -1840,13 +1840,11 @@ SUBROUTINE IncreaseMaxParticleNumber(Amount)
   IF(Amount.EQ.0) RETURN
   NewSize=PDM%MaxParticleNumber+Amount
   ! IPWRITE(*,*) "Increase by amount",PDM%MaxParticleNumber,NewSize
-  IF(NewSize.GT.PDM%maxAllowedParticleNumber)CALL ABORT(&
-  __STAMP__&
+  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__&
+  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
@@ -1862,26 +1860,26 @@ SUBROUTINE IncreaseMaxParticleNumber(Amount)
 IF(ALLOCATED(PDM%InRotRefFrame)) CALL ChangeSizeArray(PDM%InRotRefFrame,PDM%maxParticleNumber,NewSize,.FALSE.)
 
 IF(ALLOCATED(PartState)) CALL ChangeSizeArray(PartState,PDM%maxParticleNumber,NewSize,0.)
-IF(ALLOCATED(LastPartPos)) CALL ChangeSizeArray(LastPartPos,PDM%maxParticleNumber,NewSize)
+IF(ALLOCATED(LastPartPos)) CALL ChangeSizeArray(LastPartPos,PDM%maxParticleNumber,NewSize,0.)
 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(PartTimeStep)) CALL ChangeSizeArray(PartTimeStep,PDM%maxParticleNumber,NewSize,1.)
+IF(ALLOCATED(PartMPF)) CALL ChangeSizeArray(PartMPF,PDM%maxParticleNumber,NewSize,1.)
 IF(ALLOCATED(PartVeloRotRef)) CALL ChangeSizeArray(PartVeloRotRef,PDM%maxParticleNumber,NewSize,0.)
 IF(ALLOCATED(LastPartVeloRotRef)) CALL ChangeSizeArray(LastPartVeloRotRef,PDM%maxParticleNumber,NewSize,0.)
-IF(ALLOCATED(PartStateIntEn)) CALL ChangeSizeArray(PartStateIntEn,PDM%maxParticleNumber,NewSize)
+IF(ALLOCATED(PartStateIntEn)) CALL ChangeSizeArray(PartStateIntEn,PDM%maxParticleNumber,NewSize,0.)
 
 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(FieldAtParticle)) CALL ChangeSizeArray(FieldAtParticle,PDM%maxParticleNumber,NewSize,0.)
 
-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(InterPlanePartIndx)) CALL ChangeSizeArray(InterPlanePartIndx,PDM%maxParticleNumber,NewSize,0)
+IF(ALLOCATED(BGGas%PairingPartner)) CALL ChangeSizeArray(BGGas%PairingPartner,PDM%maxParticleNumber,NewSize,0)
+IF(ALLOCATED(CollInf%OldCollPartner)) CALL ChangeSizeArray(CollInf%OldCollPartner,PDM%maxParticleNumber,NewSize,0)
 IF(ALLOCATED(ElecRelaxPart)) CALL ChangeSizeArray(ElecRelaxPart,PDM%maxParticleNumber,NewSize,.TRUE.)
 
 #if (PP_TimeDiscMethod==508) || (PP_TimeDiscMethod==509)
-IF(ALLOCATED(velocityAtTime)) CALL ChangeSizeArray(velocityAtTime,PDM%maxParticleNumber,NewSize)
+IF(ALLOCATED(velocityAtTime)) CALL ChangeSizeArray(velocityAtTime,PDM%maxParticleNumber,NewSize,0.)
 #endif
 
 #if USE_MPI
@@ -1921,9 +1919,7 @@ SUBROUTINE IncreaseMaxParticleNumber(Amount)
 
 IF(ALLOCATED(AmbipolElecVelo)) THEN
   ALLOCATE(AmbipolElecVelo_New(NewSize),STAT=ALLOCSTAT)
-  IF (ALLOCSTAT.NE.0) CALL ABORT(&
-__STAMP__&
-,'Cannot allocate increased Array in IncreaseMaxParticleNumber')
+  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
@@ -1933,9 +1929,7 @@ SUBROUTINE IncreaseMaxParticleNumber(Amount)
 
 IF(ALLOCATED(ElectronicDistriPart)) THEN
   ALLOCATE(ElectronicDistriPart_New(NewSize),STAT=ALLOCSTAT)
-  IF (ALLOCSTAT.NE.0) CALL ABORT(&
-__STAMP__&
-,'Cannot allocate increased Array in IncreaseMaxParticleNumber')
+  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
@@ -1945,9 +1939,7 @@ SUBROUTINE IncreaseMaxParticleNumber(Amount)
 
 IF(ALLOCATED(VibQuantsPar)) THEN
   ALLOCATE(VibQuantsPar_New(NewSize),STAT=ALLOCSTAT)
-  IF (ALLOCSTAT.NE.0) CALL ABORT(&
-__STAMP__&
-,'Cannot allocate increased Array in IncreaseMaxParticleNumber')
+  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
@@ -2046,9 +2038,7 @@ SUBROUTINE ReduceMaxParticleNumber()
         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')
+        IF(ii.EQ.0.OR.ii.GT.NewSize) CALL ABORT(__STAMP__,'Error in ReduceMaxParticleNumber: New particle index outside of possible range')
       END IF
       IF(PDM%ParticleVecLength.LT.ii) PDM%ParticleVecLength = ii
       CALL ChangePartID(i,ii)
@@ -2068,26 +2058,26 @@ SUBROUTINE ReduceMaxParticleNumber()
 IF(ALLOCATED(PDM%InRotRefFrame)) CALL ChangeSizeArray(PDM%InRotRefFrame,PDM%maxParticleNumber,NewSize,.FALSE.)
 
 IF(ALLOCATED(PartState)) CALL ChangeSizeArray(PartState,PDM%maxParticleNumber,NewSize,0.)
-IF(ALLOCATED(LastPartPos)) CALL ChangeSizeArray(LastPartPos,PDM%maxParticleNumber,NewSize)
+IF(ALLOCATED(LastPartPos)) CALL ChangeSizeArray(LastPartPos,PDM%maxParticleNumber,NewSize,0.)
 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(PartTimeStep)) CALL ChangeSizeArray(PartTimeStep,PDM%maxParticleNumber,NewSize,1.)
+IF(ALLOCATED(PartMPF)) CALL ChangeSizeArray(PartMPF,PDM%maxParticleNumber,NewSize,1.)
 IF(ALLOCATED(PartVeloRotRef)) CALL ChangeSizeArray(PartVeloRotRef,PDM%maxParticleNumber,NewSize,0.)
 IF(ALLOCATED(LastPartVeloRotRef)) CALL ChangeSizeArray(LastPartVeloRotRef,PDM%maxParticleNumber,NewSize,0.)
-IF(ALLOCATED(PartStateIntEn)) CALL ChangeSizeArray(PartStateIntEn,PDM%maxParticleNumber,NewSize)
+IF(ALLOCATED(PartStateIntEn)) CALL ChangeSizeArray(PartStateIntEn,PDM%maxParticleNumber,NewSize,0.)
 
 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(FieldAtParticle)) CALL ChangeSizeArray(FieldAtParticle,PDM%maxParticleNumber,NewSize,0.)
 
-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(InterPlanePartIndx)) CALL ChangeSizeArray(InterPlanePartIndx,PDM%maxParticleNumber,NewSize,0)
+IF(ALLOCATED(BGGas%PairingPartner)) CALL ChangeSizeArray(BGGas%PairingPartner,PDM%maxParticleNumber,NewSize,0)
+IF(ALLOCATED(CollInf%OldCollPartner)) CALL ChangeSizeArray(CollInf%OldCollPartner,PDM%maxParticleNumber,NewSize,0)
 IF(ALLOCATED(ElecRelaxPart)) CALL ChangeSizeArray(ElecRelaxPart,PDM%maxParticleNumber,NewSize,.TRUE.)
 
 #if (PP_TimeDiscMethod==508) || (PP_TimeDiscMethod==509)
-IF(ALLOCATED(velocityAtTime)) CALL ChangeSizeArray(velocityAtTime,PDM%maxParticleNumber,NewSize)
+IF(ALLOCATED(velocityAtTime)) CALL ChangeSizeArray(velocityAtTime,PDM%maxParticleNumber,NewSize,0.)
 #endif
 
 #if USE_MPI
@@ -2127,9 +2117,7 @@ SUBROUTINE ReduceMaxParticleNumber()
 
 IF(ALLOCATED(AmbipolElecVelo)) THEN
   ALLOCATE(AmbipolElecVelo_New(NewSize),STAT=ALLOCSTAT)
-  IF (ALLOCSTAT.NE.0) CALL ABORT(&
-__STAMP__&
-,'Cannot allocate increased Array in ReduceMaxParticleNumber')
+  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
@@ -2142,9 +2130,7 @@ SUBROUTINE ReduceMaxParticleNumber()
 
 IF(ALLOCATED(ElectronicDistriPart)) THEN
   ALLOCATE(ElectronicDistriPart_New(NewSize),STAT=ALLOCSTAT)
-  IF (ALLOCSTAT.NE.0) CALL ABORT(&
-__STAMP__&
-,'Cannot allocate increased Array in ReduceMaxParticleNumber')
+  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
@@ -2157,9 +2143,7 @@ SUBROUTINE ReduceMaxParticleNumber()
 
 IF(ALLOCATED(VibQuantsPar)) THEN
   ALLOCATE(VibQuantsPar_New(NewSize),STAT=ALLOCSTAT)
-  IF (ALLOCSTAT.NE.0) CALL ABORT(&
-__STAMP__&
-,'Cannot allocate increased Array in ReduceMaxParticleNumber')
+  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

src/particles/pic/interpolation/pic_interpolation.f90

@@ -181,6 +181,7 @@ SUBROUTINE InitializeParticleInterpolation
 !--- Allocate arrays for interpolation of fields to particles
 ALLOCATE(FieldAtParticle(1:6,1:PDM%maxParticleNumber), STAT=ALLOCSTAT)
 IF (ALLOCSTAT.NE.0) CALL abort(__STAMP__ ,'ERROR in pic_interpolation.f90: Cannot allocate FieldAtParticle array!',ALLOCSTAT)
+FieldAtParticle = 0.
 
 SELECT CASE(TRIM(InterpolationType))
 CASE('particle_position')