From e41a022dd48d4734e8cd1dde2c1a863771ea33c2 Mon Sep 17 00:00:00 2001 From: Paul Nizenkov Date: Thu, 7 Mar 2024 00:58:25 +0100 Subject: [PATCH 1/2] Fix of BoundaryParticleOutput in combination with species-specific time step, additional regression test --- REGGIE.md | 1 + .../BPO_SpeciesTimeStep/PartAnalyze_ref.csv | 3 + .../SurfaceAnalyze_ref.csv | 3 + .../CHE_DSMC/BPO_SpeciesTimeStep/analyze.ini | 6 + .../BPO_SpeciesTimeStep/command_line.ini | 2 + .../BPO_SpeciesTimeStep/externals.ini | 6 + .../CHE_DSMC/BPO_SpeciesTimeStep/hopr.ini | 39 +++++++ .../BPO_SpeciesTimeStep/parameter.ini | 109 ++++++++++++++++++ .../CHE_DSMC/BPO_SpeciesTimeStep/readme.md | 4 + .../analyze/particle_analyze_tools.f90 | 14 +-- .../emission/particle_surface_flux.f90 | 8 +- src/particles/particle_operations.f90 | 31 ++--- src/particles/particle_timestep.f90 | 2 +- .../surfacemodel/surfacemodel_analyze.f90 | 5 +- 14 files changed, 199 insertions(+), 34 deletions(-) create mode 100644 regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/PartAnalyze_ref.csv create mode 100644 regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/SurfaceAnalyze_ref.csv create mode 100644 regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/analyze.ini create mode 100644 regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/command_line.ini create mode 100644 regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/externals.ini create mode 100644 regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/hopr.ini create mode 100644 regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/parameter.ini create mode 100644 regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/readme.md diff --git a/REGGIE.md b/REGGIE.md index 3a91e50d5..3872508b4 100644 --- a/REGGIE.md +++ b/REGGIE.md @@ -53,6 +53,7 @@ Small test cases to check features with DSMC timedisc: [Link to build](regressio | | BC_PorousBC | | PorousBC as a pump with 2 species, hard compiled N=1 | nProcs=3 | Total # of removed part through BC | | | | BC_PorousBC_2DAxi | | PorousBC as a pump with 2 species (axisymmetric, with/without radial weighting), hard compiled N=1 | nProcs=1,2 | Total number density | [Link](regressioncheck/CHE_DSMC/BC_PorousBC_2DAxi/readme.md) | | | BC_RotationalPeriodic | | Rotationally periodic BC with "worst-case" mesh based on tetrahedrons | nProcs=1,5 | Particle number | [Link](regressioncheck/CHE_DSMC/BC_RotationalPeriodic/readme.md) | +| | BPO_SpeciesTimeStep | | Species-specific time step with BoundaryParticleOutput | nProcs=4 | PartAnalyze, SurfaceAnalyze | [Link](regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/readme.md) | | | cube | | Collismode=2,3, hard compiled N=1 | nProcs=2 | | | | | Rotational_Reference_Frame_Regions | | Rotational reference frame with several regions, switching between stationary and rotating frame | nProcs=1,2,3,4 | Particle trajectory | [Link](regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Regions/readme.md) | | | Rotational_Reference_Frame_RotBC | | Rotational reference frame in combination with the rotationally periodic BC | nProcs=1,2,3,4 | Particle trajectory | [Link](regressioncheck/CHE_DSMC/Rotational_Reference_Frame_RotBC/readme.md) | diff --git a/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/PartAnalyze_ref.csv b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/PartAnalyze_ref.csv new file mode 100644 index 000000000..d2341d0a2 --- /dev/null +++ b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/PartAnalyze_ref.csv @@ -0,0 +1,3 @@ +001-TIME,002-Current-Spec-001-SF-001 +0.0000000000000000E+000,0.0000000000000000E+000 +0.1000000000000000E-008,0.2000000166752067E+001 diff --git a/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/SurfaceAnalyze_ref.csv b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/SurfaceAnalyze_ref.csv new file mode 100644 index 000000000..13f32dedf --- /dev/null +++ b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/SurfaceAnalyze_ref.csv @@ -0,0 +1,3 @@ +001-TIME,002-Flux-Spec-001-BC_Xplus,003-TotalElectricCurrent-BC_Xplus +0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000 +0.1000000000000000E-008,0.1248281052631579E+020,-.1999966605370011E+001 diff --git a/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/analyze.ini b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/analyze.ini new file mode 100644 index 000000000..c50cda33f --- /dev/null +++ b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/analyze.ini @@ -0,0 +1,6 @@ +! compare the last line of PartAnalyze.csv with a reference file +compare_data_file_name = PartAnalyze.csv, SurfaceAnalyze.csv +compare_data_file_reference = PartAnalyze_ref.csv, SurfaceAnalyze_ref.csv +compare_data_file_tolerance = 1e-3 +compare_data_file_tolerance_type = relative +compare_data_file_one_diff_per_run=F \ No newline at end of file diff --git a/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/command_line.ini b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/command_line.ini new file mode 100644 index 000000000..f0056b1d6 --- /dev/null +++ b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/command_line.ini @@ -0,0 +1,2 @@ +MPI=4 +database = ../../../SpeciesDatabase.h5 diff --git a/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/externals.ini b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/externals.ini new file mode 100644 index 000000000..0bdc5a5c4 --- /dev/null +++ b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/externals.ini @@ -0,0 +1,6 @@ +! --- Externals Tool Reggie +MPI = 1 +externalbinary = ./hopr/build/bin/hopr +externaldirectory = hopr.ini +externalruntime = pre +cmd_suffix = \ No newline at end of file diff --git a/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/hopr.ini b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/hopr.ini new file mode 100644 index 000000000..267b996a5 --- /dev/null +++ b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/hopr.ini @@ -0,0 +1,39 @@ +ProjectName = channel +Debugvisu = T +DebugVisuLevel=1 +NVisu =1 +Mode =1 + +DEFVAR = (REAL): minus_x = 0.0 +DEFVAR = (REAL): plus_x = 10.0 + +DEFVAR = (REAL): minus_y = 0.0 +DEFVAR = (REAL): plus_y = 1.0 + +DEFVAR = (REAL): minus_z = 0.0 +DEFVAR = (REAL): plus_z = 1.0 + +Corner =(/minus_x,minus_y,minus_z ,, plus_x,minus_y,minus_z ,, plus_x,plus_y,minus_z ,, minus_x,plus_y,minus_z ,, minus_x,minus_y,plus_z ,, plus_x,minus_y,plus_z ,, plus_x,plus_y,plus_z ,, minus_x,plus_y,plus_z /) +nElems =(/2,2,2/) +elemtype =108 + +BCIndex =(/6 ,4 ,1 ,3 ,2 ,5/) +! =(/z-,y-,x+,y+,x-,z+/) +nZones = 1 +nUserDefinedBoundaries=6 +BoundaryName=BC_Xplus +BoundaryType=(/3,0,0,0/) +BoundaryName=BC_Xminus +BoundaryType=(/3,0,0,0/) +BoundaryName=BC_Yplus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Yminus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Zplus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Zminus +BoundaryType=(/4,0,0,0/) + +postscalemesh=true +meshscale=1e-5 +jacobiantolerance=1e-20 diff --git a/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/parameter.ini b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/parameter.ini new file mode 100644 index 000000000..8fcb2817a --- /dev/null +++ b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/parameter.ini @@ -0,0 +1,109 @@ +IniExactFunc = 0 +N = 1 ! Polynomial degree +NAnalyze = 1 ! Number of analyze points +! =============================================================================== ! +! PICLAS2VTK +! =============================================================================== ! +NVisu = 1 +VisuParticles = T +TimeStampLength = 14 +! =============================================================================== ! +! MESH +! =============================================================================== ! +MeshFile = channel_mesh.h5 +useCurveds = F +! if boundaries have to be changed (else they are used from Mesh directly): +TrackingMethod = triatracking +! =============================================================================== ! +! OUTPUT / VISUALIZATION +! =============================================================================== ! +ProjectName = SurfFlux_Tria_EmissionCurrent +IterDisplayStep = 20 +Part-AnalyzeStep = 20 +Surface-AnalyzeStep = 20 +CalcSurfFluxInfo = T +! =============================================================================== ! +! CALCULATION +! =============================================================================== ! +tend = 1.0E-9 +Analyze_dt = 1.0E-9 +! =============================================================================== ! +! Load Balance +! =============================================================================== ! +DoLoadBalance = T +DoInitialAutoRestart = T +Load-DeviationThreshold = 1e-9 +! =============================================================================== ! +! PARTICLE BOUNDARY +! =============================================================================== ! +Part-nBounds=6 +Part-Boundary1-SourceName=BC_Xplus +Part-Boundary1-Condition=open +Part-Boundary2-SourceName=BC_Xminus +Part-Boundary2-Condition=reflective +Part-Boundary3-SourceName=BC_Yplus +Part-Boundary3-Condition=reflective +Part-Boundary4-SourceName=BC_Yminus +Part-Boundary4-Condition=reflective +Part-Boundary5-SourceName=BC_Zplus +Part-Boundary5-Condition=symmetric +Part-Boundary6-SourceName=BC_Zminus +Part-Boundary6-Condition=symmetric +Part-FIBGMdeltas=(/1e-5,5e-6,5e-6/) + +CalcBoundaryParticleOutput = T +BPO-NPartBoundaries = 1 +BPO-PartBoundaries = (/1/) +BPO-NSpecies = 1 +BPO-Species = (/1/) +! =============================================================================== ! +! DISCRETIZATION +! =============================================================================== ! +Part-maxParticleNumber=500000 +Part-Species$-MacroParticleFactor = 5E3,1E4 + +ManualTimeStep = 5.0000E-11 +Part-Species1-TimeStepFactor = 0.2,0.5 +! =============================================================================== ! +! SPECIES +! =============================================================================== ! +Part-nSpecies = 2 +Particles-Species-Database = SpeciesDatabase.h5 +! =============================================================================== ! +! Species1 - electron +! =============================================================================== ! +Part-Species1-SpeciesName = electron +Part-Species1-UseCollXSec = T + +Part-Species1-nSurfaceFluxBCs=1 +Part-Species1-Surfaceflux1-BC=2 +Part-Species1-Surfaceflux1-VeloIC = 1E7 +Part-Species1-Surfaceflux1-VeloVecIC = (/1,0,0/) +Part-Species1-Surfaceflux1-velocityDistribution = maxwell_lpn +Part-Species1-Surfaceflux1-MWTemperatureIC = 5. +Part-Species1-Surfaceflux1-EmissionCurrent = 2. +! =============================================================================== ! +! Species2 - Argon +! =============================================================================== ! +Part-Species2-SpeciesName = Ar + +Part-Species2-nInits = 1 +Part-Species2-Init1-SpaceIC = background +Part-Species2-Init1-velocityDistribution = maxwell_lpn +Part-Species2-Init1-MWTemperatureIC = 300 +Part-Species2-Init1-PartDensity = 1E10 +Part-Species2-Init1-VeloIC = 0 +Part-Species2-Init1-VeloVecIC = (/0.,0.,1./) +! =============================================================================== ! +! DSMC +! =============================================================================== ! +Particles-HaloEpsVelo=5.0E+07 +Particles-DSMC-CalcSurfaceVal=F +UseDSMC=true +Particles-DSMC-CollisMode = 2 +Part-NumberOfRandomSeeds=2 +Particles-RandomSeed1=1 +Particles-RandomSeed2=2 +Particles-DSMC-UseOctree=F +Particles-DSMC-UseNearestNeighbour = F +Particles-DSMC-CalcQualityFactors = F \ No newline at end of file diff --git a/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/readme.md b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/readme.md new file mode 100644 index 000000000..5c543d0f7 --- /dev/null +++ b/regressioncheck/CHE_DSMC/BPO_SpeciesTimeStep/readme.md @@ -0,0 +1,4 @@ +# Boundary Particle Output with Species Time Step +* Particle emission using surface flux by defining an emission current of 2A +* Species specific time step for the electrons (factor 0.2) +* Check whether the resulting current at the oppsite boundary is determined correctly \ No newline at end of file diff --git a/src/particles/analyze/particle_analyze_tools.f90 b/src/particles/analyze/particle_analyze_tools.f90 index 0e88f9702..d1556008a 100644 --- a/src/particles/analyze/particle_analyze_tools.f90 +++ b/src/particles/analyze/particle_analyze_tools.f90 @@ -1284,8 +1284,7 @@ SUBROUTINE CalcSurfaceFluxInfo() USE MOD_Globals USE MOD_TimeDisc_Vars ,ONLY: dt, iter USE MOD_Particle_Analyze_Vars ,ONLY: FlowRateSurfFlux,PressureAdaptiveBC,PartAnalyzeStep -USE MOD_DSMC_Vars ,ONLY: RadialWeighting -USE MOD_Particle_Vars ,ONLY: Species,nSpecies,usevMPF,VarTimeStep +USE MOD_Particle_Vars ,ONLY: Species,nSpecies,VarTimeStep USE MOD_Particle_Surfaces_Vars ,ONLY: BCdata_auxSF, SurfFluxSideSize, SurfMeshSubSideData USE MOD_Particle_Sampling_Vars ,ONLY: UseAdaptiveBC, AdaptBCMacroVal, AdaptBCMapElemToSample, AdaptBCAreaSurfaceFlux USE MOD_Particle_Sampling_Vars ,ONLY: AdaptBCAverageValBC, AdaptBCAverageMacroVal @@ -1302,7 +1301,7 @@ SUBROUTINE CalcSurfaceFluxInfo() !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: iSpec, iSF, ElemID, SampleElemID, SurfSideID, iSide, iSample, jSample, currentBC -REAL :: MacroParticleFactor, dtVar +REAL :: dtVar #if USE_MPI INTEGER :: MaxSurfaceFluxBCs #endif /*USE_MPI*/ @@ -1313,12 +1312,7 @@ SUBROUTINE CalcSurfaceFluxInfo() IF(UseAdaptiveBC) PressureAdaptiveBC = 0. ! 1) Calculate the processor-local mass flow rate and sum-up the area weighted pressure DO iSpec = 1, nSpecies - ! If usevMPF or DoRadialWeighting then the MacroParticleFactor is already included in the GetParticleWeight - IF(usevMPF.OR.RadialWeighting%DoRadialWeighting) THEN - MacroParticleFactor = 1. - ELSE - MacroParticleFactor = Species(iSpec)%MacroParticleFactor - END IF + ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific) THEN dtVar = dt * Species(iSpec)%TimeStepFactor ELSE @@ -1326,7 +1320,7 @@ SUBROUTINE CalcSurfaceFluxInfo() END IF DO iSF = 1, Species(iSpec)%nSurfacefluxBCs ! SampledMassFlow contains the weighted particle number balance (in - out) - FlowRateSurfFlux(iSpec,iSF) = Species(iSpec)%Surfaceflux(iSF)%SampledMassflow * MacroParticleFactor / dtVar + FlowRateSurfFlux(iSpec,iSF) = Species(iSpec)%Surfaceflux(iSF)%SampledMassflow / dtVar IF(Species(iSpec)%Surfaceflux(iSF)%UseEmissionCurrent) THEN FlowRateSurfFlux(iSpec,iSF) = FlowRateSurfFlux(iSpec,iSF) * ABS(Species(iSpec)%ChargeIC) ELSE diff --git a/src/particles/emission/particle_surface_flux.f90 b/src/particles/emission/particle_surface_flux.f90 index d9ca9af02..39b090fba 100644 --- a/src/particles/emission/particle_surface_flux.f90 +++ b/src/particles/emission/particle_surface_flux.f90 @@ -71,6 +71,7 @@ SUBROUTINE ParticleSurfaceflux() INTEGER :: PartInsSideRadWeight(1:RadialWeighting%nSubSides) REAL :: Particle_pos(3), RandVal1, xyzNod(3), RVec(2), minPos(2), xi(2), Vector1(3), Vector2(3) REAL :: ndist(3), midpoint(3) +REAL :: MPF LOGICAL :: AcceptPos REAL,ALLOCATABLE :: particle_positions(:), particle_xis(:) INTEGER,ALLOCATABLE :: PartInsSubSides(:,:,:) @@ -244,7 +245,12 @@ SUBROUTINE ParticleSurfaceflux() PartMPF(ParticleIndexNbr) = CalcRadWeightMPF(PartState(2,ParticleIndexNbr), iSpec,ParticleIndexNbr) END IF IF(CalcSurfFluxInfo) THEN - SF%SampledMassflow = SF%SampledMassflow + GetParticleWeight(ParticleIndexNbr) + IF(usevMPF.OR.RadialWeighting%DoRadialWeighting) THEN + MPF = GetParticleWeight(ParticleIndexNbr) + ELSE + MPF = GetParticleWeight(ParticleIndexNbr) * Species(iSpec)%MacroParticleFactor + END IF + SF%SampledMassflow = SF%SampledMassflow + MPF END IF #ifdef CODE_ANALYZE CALL AnalyzePartPos(ParticleIndexNbr) diff --git a/src/particles/particle_operations.f90 b/src/particles/particle_operations.f90 index 868f01312..eeb51a805 100644 --- a/src/particles/particle_operations.f90 +++ b/src/particles/particle_operations.f90 @@ -137,7 +137,7 @@ SUBROUTINE RemoveParticle(PartID,BCID,alpha,crossedBC) !=================================================================================================================================== ! MODULES USE MOD_Globals_Vars ,ONLY: ElementaryCharge -USE MOD_Particle_Vars ,ONLY: PDM, PartSpecies, Species, PartMPF, usevMPF, PartState, PartPosRef, Pt +USE MOD_Particle_Vars ,ONLY: PDM, PartSpecies, Species, usevMPF, PartState, PartPosRef, Pt USE MOD_Particle_Sampling_Vars ,ONLY: UseAdaptiveBC, AdaptBCPartNumOut USE MOD_Particle_Vars ,ONLY: UseNeutralization, NeutralizationSource, NeutralizationBalance,nNeutralizationElems USE MOD_Particle_Boundary_Vars ,ONLY: PartBound @@ -152,14 +152,14 @@ SUBROUTINE RemoveParticle(PartID,BCID,alpha,crossedBC) #endif USE MOD_Particle_Analyze_Tools ,ONLY: CalcEkinPart USE MOD_part_tools ,ONLY: GetParticleWeight -USE MOD_DSMC_Vars ,ONLY: CollInf, AmbipolElecVelo, ElectronicDistriPart, VibQuantsPar +USE MOD_DSMC_Vars ,ONLY: CollInf, AmbipolElecVelo, ElectronicDistriPart, VibQuantsPar, RadialWeighting USE MOD_Mesh_Vars ,ONLY: BoundaryName #if USE_HDG USE MOD_Globals ,ONLY: abort USE MOD_HDG_Vars ,ONLY: UseFPC,FPC,UseEPC,EPC USE MOD_Mesh_Vars ,ONLY: BoundaryType #endif /*USE_HDG*/ -USE MOD_Particle_Vars ,ONLY: PartState, LastPartPos +USE MOD_Particle_Vars ,ONLY: PartState!, LastPartPos !USE MOD_Particle_Mesh_Vars ,ONLY: GEO !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE @@ -229,13 +229,17 @@ SUBROUTINE RemoveParticle(PartID,BCID,alpha,crossedBC) ! - the mass flow through the boundary shall be calculated or ! - the charges impinging on the boundary are to be summed (thruster neutralization) IF(PRESENT(BCID)) THEN - + ! Determine the particle weight + IF(usevMPF.OR.RadialWeighting%DoRadialWeighting) THEN + MPF = GetParticleWeight(PartID) + ELSE + MPF = GetParticleWeight(PartID) * Species(iSpec)%MacroParticleFactor + END IF ! Check if adaptive BC or surface flux info IF(UseAdaptiveBC.OR.CalcSurfFluxInfo) THEN DO iSF=1,Species(iSpec)%nSurfacefluxBCs IF(Species(iSpec)%Surfaceflux(iSF)%BC.EQ.BCID) THEN - Species(iSpec)%Surfaceflux(iSF)%SampledMassflow = Species(iSpec)%Surfaceflux(iSF)%SampledMassflow & - - GetParticleWeight(PartID) + Species(iSpec)%Surfaceflux(iSF)%SampledMassflow = Species(iSpec)%Surfaceflux(iSF)%SampledMassflow - MPF IF(Species(iSpec)%Surfaceflux(iSF)%AdaptiveType.EQ.4) AdaptBCPartNumOut(iSpec,iSF) = AdaptBCPartNumOut(iSpec,iSF) + 1 END IF END DO @@ -254,11 +258,6 @@ SUBROUTINE RemoveParticle(PartID,BCID,alpha,crossedBC) ! Check if BPO boundary is encountered IF(CalcBoundaryParticleOutput)THEN - IF(usevMPF)THEN - MPF = PartMPF(PartID) - ELSE - MPF = Species(iSpec)%MacroParticleFactor - END IF ASSOCIATE( iBPOBC => BPO%BCIDToBPOBCID(BCID),& iBPOSpec => BPO%SpecIDToBPOSpecID(iSpec)) IF(iBPOBC.GT.0.AND.iBPOSpec.GT.0)THEN! count this species on this BC @@ -273,11 +272,6 @@ SUBROUTINE RemoveParticle(PartID,BCID,alpha,crossedBC) iBC = PartBound%MapToFieldBC(BCID) IF(iBC.LE.0) CALL abort(__STAMP__,'iBC = PartBound%MapToFieldBC(BCID) must be >0',IntInfoOpt=iBC) IF(BoundaryType(iBC,BC_TYPE).EQ.20)THEN ! BCType = BoundaryType(iBC,BC_TYPE) - IF(usevMPF)THEN - MPF = PartMPF(PartID) - ELSE - MPF = Species(iSpec)%MacroParticleFactor - END IF BCState = BoundaryType(iBC,BC_STATE) ! State is iFPC iUniqueFPCBC = FPC%Group(BCState,2) FPC%ChargeProc(iUniqueFPCBC) = FPC%ChargeProc(iUniqueFPCBC) + Species(iSpec)%ChargeIC * MPF @@ -289,11 +283,6 @@ SUBROUTINE RemoveParticle(PartID,BCID,alpha,crossedBC) iBC = PartBound%MapToFieldBC(BCID) IF(iBC.LE.0) CALL abort(__STAMP__,'iBC = PartBound%MapToFieldBC(BCID) must be >0',IntInfoOpt=iBC) IF(BoundaryType(iBC,BC_TYPE).EQ.8)THEN ! BCType = BoundaryType(iBC,BC_TYPE) - IF(usevMPF)THEN - MPF = PartMPF(PartID) - ELSE - MPF = Species(iSpec)%MacroParticleFactor - END IF BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC iUniqueEPCBC = EPC%Group(BCState,2) EPC%ChargeProc(iUniqueEPCBC) = EPC%ChargeProc(iUniqueEPCBC) + Species(iSpec)%ChargeIC * MPF diff --git a/src/particles/particle_timestep.f90 b/src/particles/particle_timestep.f90 index 6c95a5849..062c7fcd2 100644 --- a/src/particles/particle_timestep.f90 +++ b/src/particles/particle_timestep.f90 @@ -466,7 +466,7 @@ END FUNCTION GetParticleTimeStep PURE REAL FUNCTION GetSpeciesTimeStep(iCase) !=================================================================================================================================== -!> Determines the species-specific time step from the collision case +!> Determines the species-specific time step from the collision case (works only in combination with a background gas) !=================================================================================================================================== ! MODULES USE MOD_Globals diff --git a/src/particles/surfacemodel/surfacemodel_analyze.f90 b/src/particles/surfacemodel/surfacemodel_analyze.f90 index 50c9e6265..6c0833e76 100644 --- a/src/particles/surfacemodel/surfacemodel_analyze.f90 +++ b/src/particles/surfacemodel/surfacemodel_analyze.f90 @@ -151,7 +151,7 @@ SUBROUTINE AnalyzeSurface(Time) USE MOD_Restart_Vars ,ONLY: DoRestart USE MOD_Particle_Boundary_Vars ,ONLY: PartBound USE MOD_SurfaceModel_Vars ,ONLY: nPorousBC, PorousBC -USE MOD_Particle_Vars ,ONLY: nSpecies,UseNeutralization,NeutralizationBalanceGlobal,Species +USE MOD_Particle_Vars ,ONLY: nSpecies,UseNeutralization,NeutralizationBalanceGlobal,Species,VarTimeStep #if USE_MPI USE MOD_Particle_Boundary_Vars ,ONLY: SurfCOMM #endif /*USE_MPI*/ @@ -312,6 +312,9 @@ SUBROUTINE AnalyzeSurface(Time) IF(ABS(SurfModelAnalyzeSampleTime).LE.0.0)THEN CALL WriteDataInfo(unit_index,RealScalar=0.0) ELSE + ! Scaling the number of particles depending on the species-specific timestep factor + IF(VarTimeStep%UseSpeciesSpecific) BPO%RealPartOut(iBPO,iSpec) = BPO%RealPartOut(iBPO,iSpec) & + / Species(BPO%Species(iSpec))%TimeStepFactor CALL WriteDataInfo(unit_index,RealScalar=BPO%RealPartOut(iBPO,iSpec)/SurfModelAnalyzeSampleTime) END IF ! ABS(SurfModelAnalyzeSampleTime).LE.0.0 END DO From 104699d6b316770fb07b89403d23a3da7066a016 Mon Sep 17 00:00:00 2001 From: Paul Nizenkov Date: Mon, 11 Mar 2024 23:13:42 +0100 Subject: [PATCH 2/2] Utilize ParticleAnalyzeSampleTime to calculate per second variables (collision/reaction rates, mass flow etc.) instead of using the time step and PartAnalyzeStep --- .../analyze/particle_analyze_tools.f90 | 89 ++++--------------- 1 file changed, 17 insertions(+), 72 deletions(-) diff --git a/src/particles/analyze/particle_analyze_tools.f90 b/src/particles/analyze/particle_analyze_tools.f90 index d1556008a..dea073b50 100644 --- a/src/particles/analyze/particle_analyze_tools.f90 +++ b/src/particles/analyze/particle_analyze_tools.f90 @@ -1282,8 +1282,8 @@ SUBROUTINE CalcSurfaceFluxInfo() !=================================================================================================================================== ! MODULES ! USE MOD_Globals -USE MOD_TimeDisc_Vars ,ONLY: dt, iter -USE MOD_Particle_Analyze_Vars ,ONLY: FlowRateSurfFlux,PressureAdaptiveBC,PartAnalyzeStep +USE MOD_TimeDisc_Vars ,ONLY: iter +USE MOD_Particle_Analyze_Vars ,ONLY: FlowRateSurfFlux,PressureAdaptiveBC,ParticleAnalyzeSampleTime USE MOD_Particle_Vars ,ONLY: Species,nSpecies,VarTimeStep USE MOD_Particle_Surfaces_Vars ,ONLY: BCdata_auxSF, SurfFluxSideSize, SurfMeshSubSideData USE MOD_Particle_Sampling_Vars ,ONLY: UseAdaptiveBC, AdaptBCMacroVal, AdaptBCMapElemToSample, AdaptBCAreaSurfaceFlux @@ -1314,9 +1314,9 @@ SUBROUTINE CalcSurfaceFluxInfo() DO iSpec = 1, nSpecies ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific) THEN - dtVar = dt * Species(iSpec)%TimeStepFactor + dtVar = ParticleAnalyzeSampleTime * Species(iSpec)%TimeStepFactor ELSE - dtVar = dt + dtVar = ParticleAnalyzeSampleTime END IF DO iSF = 1, Species(iSpec)%nSurfacefluxBCs ! SampledMassFlow contains the weighted particle number balance (in - out) @@ -1371,15 +1371,8 @@ SUBROUTINE CalcSurfaceFluxInfo() END IF #endif /*USE_MPI*/ -! 3) Consider Part-AnalyzeStep for FlowRateSurfFlux and determine the average pressure (value does not depend on the Part-AnalyzeStep) +! 3) Determine the average pressure IF (MPIRoot) THEN - IF(PartAnalyzeStep.GT.1)THEN - IF(PartAnalyzeStep.EQ.HUGE(PartAnalyzeStep))THEN - FlowRateSurfFlux = FlowRateSurfFlux / iter - ELSE - FlowRateSurfFlux = FlowRateSurfFlux / MIN(PartAnalyzeStep,iter) - END IF - END IF IF(UseAdaptiveBC) THEN IF(AdaptBCAverageValBC) THEN PressureAdaptiveBC(:,:) = AdaptBCAverageMacroVal(3,:,:) @@ -2125,10 +2118,9 @@ SUBROUTINE CollRates(CRate) !=================================================================================================================================== ! MODULES USE MOD_Globals +USE MOD_Particle_Analyze_Vars ,ONLY: ParticleAnalyzeSampleTime USE MOD_DSMC_Vars ,ONLY: CollInf, DSMC -USE MOD_TimeDisc_Vars ,ONLY: dt, iter USE MOD_Particle_Vars ,ONLY: VarTimeStep -USE MOD_Particle_Analyze_Vars ,ONLY: PartAnalyzeStep USE MOD_Particle_TimeStep ,ONLY: GetSpeciesTimeStep ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -2156,26 +2148,14 @@ SUBROUTINE CollRates(CRate) DO iCase=1, CollInf%NumCase ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific.AND..NOT.VarTimeStep%DisableForMCC) THEN - dtVar = dt * GetSpeciesTimeStep(iCase) + dtVar = ParticleAnalyzeSampleTime * GetSpeciesTimeStep(iCase) ELSE - dtVar = dt + dtVar = ParticleAnalyzeSampleTime END IF CRate(iCase) = DSMC%NumColl(iCase) / dtVar END DO ! Total collision rate is the sum of the case-specific rates CRate(CollInf%NumCase + 1) = SUM(CRate(1:CollInf%NumCase)) - ! Consider Part-AnalyzeStep - IF(PartAnalyzeStep.GT.1)THEN - IF(PartAnalyzeStep.EQ.HUGE(PartAnalyzeStep))THEN - DO iCase=1, CollInf%NumCase + 1 - CRate(iCase) = CRate(iCase) / iter - END DO ! iCase=1, CollInf%NumCase + 1 - ELSE - DO iCase=1, CollInf%NumCase + 1 - CRate(iCase) = CRate(iCase) / MIN(PartAnalyzeStep,iter) - END DO ! iCase=1, CollInf%NumCase + 1 - END IF - END IF END IF DSMC%NumColl = 0. @@ -2194,8 +2174,7 @@ SUBROUTINE CalcRelaxRates(NumSpec,VibRelaxProbCase) USE MOD_DSMC_Vars ,ONLY: CollisMode, CollInf USE MOD_MCC_Vars ,ONLY: SpecXSec, XSec_Relaxation USE MOD_Particle_Mesh_Vars ,ONLY: MeshVolume -USE MOD_TimeDisc_Vars ,ONLY: dt, iter -USE MOD_Particle_Analyze_Vars ,ONLY: PartAnalyzeStep +USE MOD_Particle_Analyze_Vars ,ONLY: ParticleAnalyzeSampleTime USE MOD_Particle_TimeStep ,ONLY: GetSpeciesTimeStep ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -2231,27 +2210,15 @@ SUBROUTINE CalcRelaxRates(NumSpec,VibRelaxProbCase) iCase = CollInf%Coll_Case(iSpec,jSpec) ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific.AND..NOT.VarTimeStep%DisableForMCC) THEN - dtVar = dt * GetSpeciesTimeStep(iCase) + dtVar = ParticleAnalyzeSampleTime * GetSpeciesTimeStep(iCase) ELSE - dtVar = dt + dtVar = ParticleAnalyzeSampleTime END IF VibRelaxProbCase(iCase) = SpecXSec(iCase)%VibCount * MPF_1 * MeshVolume & / (dtVar * MPF_1*NumSpec(iSpec) * MPF_2*NumSpec(jSpec)) END DO END DO END IF - ! Consider Part-AnalyzeStep - IF(PartAnalyzeStep.GT.1) THEN - IF(PartAnalyzeStep.EQ.HUGE(PartAnalyzeStep))THEN - DO iCase = 1, CollInf%NumCase - VibRelaxProbCase(iCase) = VibRelaxProbCase(iCase) / iter - END DO - ELSE - DO iCase = 1, CollInf%NumCase - VibRelaxProbCase(iCase) = VibRelaxProbCase(iCase) / MIN(PartAnalyzeStep,iter) - END DO - END IF - END IF END IF SpecXSec(:)%VibCount = 0. @@ -2267,8 +2234,7 @@ SUBROUTINE CalcRelaxRatesElec(ElecRelaxRate) USE MOD_Particle_Vars ,ONLY: VarTimeStep USE MOD_DSMC_Vars ,ONLY: CollInf USE MOD_MCC_Vars ,ONLY: SpecXSec -USE MOD_TimeDisc_Vars ,ONLY: dt, iter -USE MOD_Particle_Analyze_Vars ,ONLY: PartAnalyzeStep +USE MOD_Particle_Analyze_Vars ,ONLY: ParticleAnalyzeSampleTime USE MOD_Particle_TimeStep ,ONLY: GetSpeciesTimeStep ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -2308,21 +2274,13 @@ SUBROUTINE CalcRelaxRatesElec(ElecRelaxRate) IF(SpecXSec(iCase)%UseElecXSec) THEN ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific.AND..NOT.VarTimeStep%DisableForMCC) THEN - dtVar = dt * GetSpeciesTimeStep(iCase) + dtVar = ParticleAnalyzeSampleTime * GetSpeciesTimeStep(iCase) ELSE - dtVar = dt + dtVar = ParticleAnalyzeSampleTime END IF ElecRelaxRate(iCase,:) = ElecRelaxRate(iCase,:) / dtVar END IF END DO - ! Consider Part-AnalyzeStep - IF(PartAnalyzeStep.GT.1)THEN - IF(PartAnalyzeStep.EQ.HUGE(PartAnalyzeStep))THEN - ElecRelaxRate = ElecRelaxRate / iter - ELSE - ElecRelaxRate = ElecRelaxRate / MIN(PartAnalyzeStep,iter) - END IF - END IF END IF DO iCase=1, CollInf%NumCase @@ -2342,11 +2300,10 @@ SUBROUTINE ReacRates(NumSpec, RRate) !=================================================================================================================================== ! MODULES USE MOD_Globals +USE MOD_Particle_Analyze_Vars ,ONLY: ParticleAnalyzeSampleTime USE MOD_DSMC_Vars ,ONLY: ChemReac, DSMC -USE MOD_TimeDisc_Vars ,ONLY: dt, iter USE MOD_Particle_Vars ,ONLY: Species, nSpecies, VarTimeStep USE MOD_Particle_Mesh_Vars ,ONLY: MeshVolume -USE MOD_Particle_Analyze_Vars ,ONLY: PartAnalyzeStep USE MOD_Particle_TimeStep ,ONLY: GetSpeciesTimeStep ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -2378,9 +2335,9 @@ SUBROUTINE ReacRates(NumSpec, RRate) iCase = ChemReac%ReactCase(iReac) ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific.AND..NOT.VarTimeStep%DisableForMCC) THEN - dtVar = dt * GetSpeciesTimeStep(iCase) + dtVar = ParticleAnalyzeSampleTime * GetSpeciesTimeStep(iCase) ELSE - dtVar = dt + dtVar = ParticleAnalyzeSampleTime END IF IF ((NumSpec(ChemReac%Reactants(iReac,1)).GT.0).AND.(NumSpec(ChemReac%Reactants(iReac,2)).GT.0)) THEN IF(ChemReac%Reactants(iReac,3).NE.0) THEN @@ -2420,18 +2377,6 @@ SUBROUTINE ReacRates(NumSpec, RRate) ChemReac%NumReac = 0. ChemReac%ReacCount = 0 ChemReac%ReacCollMean = 0.0 -! Consider Part-AnalyzeStep -IF(PartAnalyzeStep.GT.1)THEN - IF(PartAnalyzeStep.EQ.HUGE(PartAnalyzeStep))THEN - DO iReac=1, ChemReac%NumOfReact - RRate(iReac) = RRate(iReac) / iter - END DO ! iReac=1, ChemReac%NumOfReact - ELSE - DO iReac=1, ChemReac%NumOfReact - RRate(iReac) = RRate(iReac) / REAL(MIN(PartAnalyzeStep,iter)) - END DO ! iReac=1, ChemReac%NumOfReact - END IF -END IF END SUBROUTINE ReacRates #endif