diff --git a/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/DSMC.ini b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/DSMC.ini new file mode 100644 index 000000000..1b4af9b7f --- /dev/null +++ b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/DSMC.ini @@ -0,0 +1,17 @@ +! =============================================================================== ! +! Species1, N+ +! =============================================================================== ! +Part-Species1-SpeciesName=NIon1 +Part-Species1-InteractionID = 10 +Part-Species1-TRef =273 +Part-Species1-DRef = 3.0E-10 +Part-Species1-omega=0.24 +! =============================================================================== ! +! Species2, e +! =============================================================================== ! +Part-Species2-SpeciesName = electron +Part-Species2-InteractionID = 4 +Part-Species2-TRef =273 +Part-Species2-DRef = 2.817920E-15 +Part-Species2-omega=0.24 + diff --git a/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/analyze.ini b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/analyze.ini new file mode 100644 index 000000000..e69de29bb diff --git a/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/command_line.ini b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/command_line.ini new file mode 100644 index 000000000..416ce7a09 --- /dev/null +++ b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/command_line.ini @@ -0,0 +1,2 @@ +MPI=1,6 +cmd_suffix=DSMC.ini diff --git a/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/hopr.ini b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/hopr.ini new file mode 100644 index 000000000..4821f5bf3 --- /dev/null +++ b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/hopr.ini @@ -0,0 +1,47 @@ +!=============================================================================== ! +! MAKEFILE PARAMETER (put a "#" in front, NO blanks!) +!=============================================================================== ! +! This is only a dummy parameter needed for the regression check +#MPI= + +!=============================================================================== ! +! OUTPUT +!=============================================================================== ! + ProjectName =cube ! name of the project (used for filenames) + Debugvisu =F ! Write debug mesh to tecplot file + Logging =F ! Write log files + +!=============================================================================== ! +! MESH +!=============================================================================== ! + Mode =1 ! 1 Cartesian 2 gambit file 3 CGNS + nZones =1 ! number of zones + Corner =(/0.,0.,0.,,1E-6,0.,0.,,1E-6,1E-6,0.,,0.,1E-6,0. ,,0.,0.,1E-6,,1E-6,0.,1E-6,,1E-6,1E-6,1E-6,,0.,1E-6,1E-6/) ! [0,1]x[0,1]x[0,0.05] + nElems =(/10,10,10/) ! Anzahl der Elemente in jede Richtung (nfine 4:16 5:32 6:64 7:128) + BCIndex =(/1,1,2,2,3,3/) ! Indices of UserDefinedBoundaries + elemtype =108 ! Elementform (108: Hexaeder) + useCurveds =F ! T if curved boundaries defined + SpaceQuandt =1. ! characteristic length of the mesh + ConformConnect=T + jacobianTolerance = 1e-27 +!=============================================================================== ! +! BOUNDARY CONDITIONS +!=============================================================================== ! + nUserDefinedBoundaries=3 + BoundaryName=BC1 ! Outflow: open (absorbing) [for MAXWELL] + BoundaryType=(/4,0,0,0/) ! Outflow: open (absorbing) [for MAXWELL] + BoundaryName=BC2 ! Outflow: open (absorbing) [for MAXWELL] + BoundaryType=(/4,0,0,0/) ! Outflow: open (absorbing) [for MAXWELL] + BoundaryName=BC3 ! Outflow: open (absorbing) [for MAXWELL] + BoundaryType=(/4,0,0,0/) ! Outflow: open (absorbing) [for MAXWELL] + +!=============================================================================== ! +! BASIS +!=============================================================================== ! + NVisu = 7 + +!=============================================================================== ! +! SEARCH +!=============================================================================== ! +! nElemsNodeSearch=50 +! RefineSideSearch=50 diff --git a/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/parameter.ini b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/parameter.ini new file mode 100644 index 000000000..b12304aeb --- /dev/null +++ b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/parameter.ini @@ -0,0 +1,108 @@ +CFLscale = 0.2 +IniExactFunc = 0 +N = 1 +NAnalyze = 1 +! =============================================================================== ! +! MESH +! =============================================================================== ! +MeshFile = cube_mesh.h5 +Logging = F +useCurveds = F +! if boundaries have to be changed (else they are used from Mesh directly): +TrackingMethod = triatracking!, tracing, refmapping +! =============================================================================== ! +! OUTPUT / VISUALIZATION +! =============================================================================== ! +ProjectName = Reservoir +IterDisplayStep = 1 +Part-AnalyzeStep = 1 +CalcNumDens = T +CalcTemp = T +! =============================================================================== ! +! CALCULATION +! =============================================================================== ! +Particles-ManualTimeStep = 1E-12 +tend = 1E-10 +Analyze_dt = 1E-7 ! Timestep of analyze outputs +! =============================================================================== ! +! PARTICLES +! =============================================================================== ! +Part-maxParticleNumber=3500000 +Part-nSpecies=2 +Part-nBounds=3 +Part-Boundary1-SourceName=BC1 +Part-Boundary1-Condition=open +Part-Boundary2-SourceName=BC2 +Part-Boundary2-Condition=open +Part-Boundary3-SourceName=BC3 +Part-Boundary3-Condition=open +Part-FIBGMdeltas=(/1E-6,1E-6,1E-6/) +Particles-HaloEpsVelo=5000 +! =============================================================================== ! +! DSMC +! =============================================================================== ! +UseDSMC = T +Particles-DSMC-CollisMode = 1 +Particles-DSMC-CalcQualityFactors = F + +Particles-DSMC-AmbipolarDiffusion = T + +Particles-DSMC-ElectronicModel = 0 +! =============================================================================== ! +! Weighting Factor +! =============================================================================== ! +Part-Species1-MacroParticleFactor = 5 +Part-Species2-MacroParticleFactor = 5 +! =============================================================================== ! +! Species1 | N+ +! =============================================================================== ! +Part-Species1-ChargeIC = 1.60217653E-19 +Part-Species1-MassIC = 2.32600E-26 + +Part-Species1-nInits = 0 +Part-Species1-nSurfacefluxBCs=3 +Part-Species1-Surfaceflux1-BC=1 +Part-Species1-SurfaceFlux1-velocityDistribution = maxwell_lpn +Part-Species1-SurfaceFlux1-PartDensity = 2E23 +Part-Species1-SurfaceFlux1-VeloIC = 1000. +Part-Species1-SurfaceFlux1-VeloVecIC = (/0.,1.,0./) +Part-Species1-SurfaceFlux1-MWTemperatureIC = 1000. +Part-Species1-Surfaceflux2-BC=2 +Part-Species1-SurfaceFlux2-velocityDistribution = maxwell_lpn +Part-Species1-SurfaceFlux2-PartDensity = 2E23 +Part-Species1-SurfaceFlux2-VeloIC = 1000. +Part-Species1-SurfaceFlux2-VeloVecIC = (/0.,1.,0./) +Part-Species1-SurfaceFlux2-MWTemperatureIC = 1000. +Part-Species1-Surfaceflux3-BC=3 +Part-Species1-SurfaceFlux3-velocityDistribution = maxwell_lpn +Part-Species1-SurfaceFlux3-PartDensity = 2E23 +Part-Species1-SurfaceFlux3-VeloIC = 1000. +Part-Species1-SurfaceFlux3-VeloVecIC = (/0.,1.,0./) +Part-Species1-SurfaceFlux3-MWTemperatureIC = 1000. +! =============================================================================== ! +! Species2 | e +! =============================================================================== ! +Part-Species2-ChargeIC = -1.60217653E-19 +Part-Species2-MassIC = 9.1093826E-31 +Part-Species2-nInits = 0 +Part-Species2-nSurfacefluxBCs=3 +Part-Species2-Surfaceflux1-BC=1 +Part-Species2-SurfaceFlux1-velocityDistribution = maxwell_lpn +Part-Species2-SurfaceFlux1-PartDensity = 2E23 +Part-Species2-SurfaceFlux1-VeloIC = 1000. +Part-Species2-SurfaceFlux1-VeloVecIC = (/0.,1.,0./) +Part-Species2-SurfaceFlux1-MWTemperatureIC = 1000. +Part-Species2-Surfaceflux2-BC=2 +Part-Species2-SurfaceFlux2-velocityDistribution = maxwell_lpn +Part-Species2-SurfaceFlux2-PartDensity = 2E23 +Part-Species2-SurfaceFlux2-VeloIC = 1000. +Part-Species2-SurfaceFlux2-VeloVecIC = (/0.,1.,0./) +Part-Species2-SurfaceFlux2-MWTemperatureIC = 1000. +Part-Species2-Surfaceflux3-BC=3 +Part-Species2-SurfaceFlux3-velocityDistribution = maxwell_lpn +Part-Species2-SurfaceFlux3-PartDensity = 2E23 +Part-Species2-SurfaceFlux3-VeloIC = 1000. +Part-Species2-SurfaceFlux3-VeloVecIC = (/0.,1.,0./) +Part-Species2-SurfaceFlux3-MWTemperatureIC = 1000. + + diff --git a/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/readme.md b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/readme.md new file mode 100644 index 000000000..e36c45423 --- /dev/null +++ b/regressioncheck/NIG_DSMC/Ambipolar_Diffusion_SF/readme.md @@ -0,0 +1,4 @@ +# Surface Flux, Ambipolar Diffusion +* Empty computational Domains with surface flux of ions and electrons applied on all sides +* Electrons are treated with the ambipolar diffusion model (attached to the ion movement, but with their own velocity vector) +* Sanaty check will fail if there are a different number of electrons and ions in an element and leads to an abort diff --git a/src/particles/dsmc/dsmc_ambipolardiff.f90 b/src/particles/dsmc/dsmc_ambipolardiff.f90 index f6fb70ff4..106a82aa1 100644 --- a/src/particles/dsmc/dsmc_ambipolardiff.f90 +++ b/src/particles/dsmc/dsmc_ambipolardiff.f90 @@ -26,7 +26,7 @@ MODULE MOD_DSMC_AmbipolarDiffusion !----------------------------------------------------------------------------------------------------------------------------------- ! Private Part --------------------------------------------------------------------------------------------------------------------- ! Public Part ---------------------------------------------------------------------------------------------------------------------- -PUBLIC :: InitializeVariablesAmbipolarDiff, AD_SetInitElectronVelo, AD_InsertParticles, AD_DeleteParticles +PUBLIC :: InitializeVariablesAmbipolarDiff, AD_SetInitElectronVelo, AD_InsertParticles, AD_DeleteParticles, AD_SetSFElectronVelo !=================================================================================================================================== CONTAINS @@ -154,6 +154,301 @@ SUBROUTINE AD_SetInitElectronVelo(FractNbr,iInit,NbrOfParticle) END SUBROUTINE AD_SetInitElectronVelo +SUBROUTINE AD_SetSFElectronVelo(iSpec,iSFIon,iSample,jSample,iSide,BCSideID,SideID,ElemID,NbrOfParticle,PartIns,particle_xis) +!=================================================================================================================================== +!> Initialize the electron velocity vector during the surface flux particle insertion +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Globals_Vars, ONLY : PI, BoltzmannConst +USE MOD_Particle_Vars +USE MOD_Particle_Surfaces_Vars, ONLY : SurfMeshSubSideData, TriaSurfaceFlux +USE MOD_Particle_Surfaces, ONLY : CalcNormAndTangBezier +USE MOD_Particle_Sampling_Vars ,ONLY: AdaptBCMapElemToSample, AdaptBCMacroVal +USE MOD_DSMC_Vars ,ONLY: AmbiPolarSFMapping, AmbipolElecVelo +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,INTENT(IN) :: iSpec,iSFIon,iSample,jSample,iSide,BCSideID,SideID,ElemID,NbrOfParticle,PartIns +REAL,INTENT(IN) :: particle_xis(:) +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: i,PositionNbr,envelope,currentBC,SampleElemID, iSF,iPart +REAL :: Vec3D(3), vec_nIn(1:3), vec_t1(1:3), vec_t2(1:3) +REAL :: a,zstar,RandVal1,RandVal2(2),RandVal3(3),u,RandN,RandN_save,Velo1,Velo2,Velosq,T,beta,z +LOGICAL :: RandN_in_Mem +REAL :: projFak ! VeloVecIC projected to inwards normal of tria +REAL :: Velo_t1 ! Velo comp. of first orth. vector in tria +REAL :: Velo_t2 ! Velo comp. of second orth. vector in tria +REAL :: VeloIC +REAL :: VeloVec(1:3) +REAL :: VeloVecIC(1:3),v_thermal, pressure +!=================================================================================================================================== + +IF(PartIns.LT.1) RETURN + +iSF = AmbiPolarSFMapping(iSpec,iSFIon) + +RandN_in_Mem=.FALSE. +envelope=-1 +currentBC = Species(iSpec)%Surfaceflux(iSF)%BC + +IF (.NOT.Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal) 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) +END IF !.NOT.VeloIsNormal + +IF(.NOT.Species(iSpec)%Surfaceflux(iSF)%Adaptive) THEN + VeloIC = Species(iSpec)%Surfaceflux(iSF)%VeloIC + T = Species(iSpec)%Surfaceflux(iSF)%MWTemperatureIC + a = Species(iSpec)%Surfaceflux(iSF)%SurfFluxSubSideData(iSample,jSample,iSide)%a_nIn + projFak = Species(iSpec)%Surfaceflux(iSF)%SurfFluxSubSideData(iSample,jSample,iSide)%projFak + Velo_t1 = Species(iSpec)%Surfaceflux(iSF)%SurfFluxSubSideData(iSample,jSample,iSide)%Velo_t1 + Velo_t2 = Species(iSpec)%Surfaceflux(iSF)%SurfFluxSubSideData(iSample,jSample,iSide)%Velo_t2 +ELSE !Species(iSpec)%Surfaceflux(iSF)%Adaptive + SampleElemID = AdaptBCMapElemToSample(ElemID) + SELECT CASE(Species(iSpec)%Surfaceflux(iSF)%AdaptiveType) + CASE(1,3,4) ! Pressure and massflow inlet (pressure/massflow, temperature const) + T = Species(iSpec)%Surfaceflux(iSF)%MWTemperatureIC + CASE(2) ! adaptive Outlet/freestream + pressure = Species(iSpec)%Surfaceflux(iSF)%AdaptivePressure + T = pressure / (BoltzmannConst * AdaptBCMacroVal(4,SampleElemID,iSpec)) + CASE DEFAULT + CALL abort(__STAMP__,'ERROR in SurfaceFlux: Wrong adaptive type for Surfaceflux velocities!') + END SELECT + VeloVec(1) = AdaptBCMacroVal(DSMC_VELOX,SampleElemID,iSpec) + VeloVec(2) = AdaptBCMacroVal(DSMC_VELOY,SampleElemID,iSpec) + VeloVec(3) = AdaptBCMacroVal(DSMC_VELOZ,SampleElemID,iSpec) + vec_nIn(1:3) = SurfMeshSubSideData(iSample,jSample,BCSideID)%vec_nIn(1:3) + VeloVec(1:3) = DOT_PRODUCT(VeloVec,vec_nIn)*vec_nIn(1:3) + VeloIC = SQRT(DOT_PRODUCT(VeloVec,VeloVec)) + IF (ABS(VeloIC).GT.0.) THEN + VeloVecIC = VeloVec / VeloIC + ELSE + VeloVecIC = (/1.,0.,0./) + END IF + projFak = DOT_PRODUCT(vec_nIn,VeloVecIC) !VeloVecIC projected to inwards normal + v_thermal = SQRT(2.*BoltzmannConst*T/Species(iSpec)%MassIC) !thermal speed + IF ( ALMOSTEQUAL(v_thermal,0.)) THEN + v_thermal = 1. + END IF + a = VeloIC * projFak / v_thermal !speed ratio proj. to inwards n (can be negative!) + Velo_t1 = VeloIC * DOT_PRODUCT(vec_t1,VeloVecIC) !v in t1-dir + Velo_t2 = VeloIC * DOT_PRODUCT(vec_t2,VeloVecIC) !v in t2-dir +END IF !Adaptive SurfaceFlux + +! Set velocities +SELECT CASE(TRIM(Species(iSpec)%Surfaceflux(iSF)%velocityDistribution)) +CASE('constant') + IF (.NOT.Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal) THEN + VeloVecIC(1:3) = Species(iSpec)%Surfaceflux(iSF)%VeloVecIC(1:3) + VeloVecIC(1:3) = VeloVecIC(1:3) / VECNORM(VeloVecIC(1:3)) + END IF + 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(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 ) + 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(iSpec)%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 + END DO !i = ...NbrOfParticle +CASE('maxwell','maxwell_lpn') + !-- determine envelope for most efficient ARM [Garcia and Wagner 2006, JCP217-2] + IF (ALMOSTZERO(VeloIC*projFak)) THEN + ! Rayleigh distri + envelope = 0 + ELSE IF (-0.4.LT.a .AND. a.LT.1.3) THEN + ! low speed flow + IF (a.LE.0.) THEN + envelope = 1 + ELSE + envelope = 3 + END IF !choose envelope based on flow direction + ELSE + ! high speed / general flow + IF (a.LT.0.) THEN + envelope = 2 + ELSE + envelope = 4 + END IF !choose envelope based on flow direction + END IF !low speed / high speed / rayleigh flow + + 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(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 ) + 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 + 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 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 +! 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))) + 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 +! 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 +! 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) +! 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) +! 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) ) + 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 + END DO !i = ...NbrOfParticle +CASE DEFAULT + CALL abort(__STAMP__,'ERROR in SurfaceFlux: Wrong velocity distribution!') +END SELECT + +END SUBROUTINE AD_SetSFElectronVelo + + SUBROUTINE AD_InsertParticles(iPartIndx_Node, nPart, iPartIndx_NodeTotalAmbi, TotalPartNum) !=================================================================================================================================== !> Creating electrons for each actual ion simulation particle, using the stored velocity vector for the electron diff --git a/src/particles/dsmc/dsmc_vars.f90 b/src/particles/dsmc/dsmc_vars.f90 index cc83180a2..5bdd498a6 100644 --- a/src/particles/dsmc/dsmc_vars.f90 +++ b/src/particles/dsmc/dsmc_vars.f90 @@ -458,6 +458,7 @@ MODULE MOD_DSMC_Vars END TYPE TYPE (tAmbipolElecVelo), ALLOCATABLE :: AmbipolElecVelo(:) +INTEGER, ALLOCATABLE :: AmbiPolarSFMapping(:,:) INTEGER, ALLOCATABLE :: iPartIndx_NodeNewAmbi(:) INTEGER :: newAmbiParts diff --git a/src/particles/emission/particle_surface_flux.f90 b/src/particles/emission/particle_surface_flux.f90 index 730443788..4ab248cf7 100644 --- a/src/particles/emission/particle_surface_flux.f90 +++ b/src/particles/emission/particle_surface_flux.f90 @@ -36,7 +36,7 @@ SUBROUTINE ParticleSurfaceflux() USE MOD_Globals USE MOD_Particle_Vars USE MOD_part_tools ,ONLY: CalcRadWeightMPF -USE MOD_DSMC_Vars ,ONLY: useDSMC, CollisMode, RadialWeighting +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 @@ -48,6 +48,7 @@ SUBROUTINE ParticleSurfaceflux() USE MOD_Particle_Surfaces_Vars ,ONLY: SurfFluxSideSize, TriaSurfaceFlux, BCdata_auxSF USE MOD_Particle_VarTimeStep ,ONLY: CalcVarTimeStep USE MOD_Timedisc_Vars ,ONLY: RKdtFrac, dt +USE MOD_DSMC_AmbipolarDiffusion ,ONLY: AD_SetSFElectronVelo #if defined(IMPA) || defined(ROS) USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod #endif /*IMPA*/ @@ -81,6 +82,9 @@ SUBROUTINE ParticleSurfaceflux() #endif /*USE_LOADBALANCE*/ !=================================================================================================================================== DO iSpec=1,nSpecies + IF (DSMC%DoAmbipolarDiff) THEN + IF (iSpec.EQ.DSMC%AmbiDiffElecSpec) CYCLE + END IF DO iSF=1,Species(iSpec)%nSurfacefluxBCs PartsEmitted = 0 currentBC = Species(iSpec)%Surfaceflux(iSF)%BC @@ -246,12 +250,24 @@ SUBROUTINE ParticleSurfaceflux() CALL abort(__STAMP__,'ERROR in ParticleSurfaceflux: ParticleIndexNbr.EQ.0 - maximum nbr of particles reached?') END IF END DO - DEALLOCATE(particle_positions) - IF (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. .NOT.TriaSurfaceFlux) DEALLOCATE(particle_xis) !----- 2a.: set velocities if special for each subside CALL SetSurfacefluxVelocities(iSpec,iSF,iSample,jSample,iSide,BCSideID,SideID,ElemID,NbrOfParticle,PartInsSubSide) PartsEmitted = PartsEmitted + PartInsSubSide + + 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 (Species(iSpec)%Surfaceflux(iSF)%VeloIsNormal .AND. .NOT.TriaSurfaceFlux) DEALLOCATE(particle_xis) + DEALLOCATE(particle_positions) #if USE_LOADBALANCE !used for calculating LoadBalance of tCurrent(LB_SURFFLUX) ==> "2b.: set remaining properties" nSurfacefluxPerElem(ElemID)=nSurfacefluxPerElem(ElemID)+PartInsSubSide diff --git a/src/particles/emission/particle_surface_flux_init.f90 b/src/particles/emission/particle_surface_flux_init.f90 index c7efe6bfd..d680cfbae 100644 --- a/src/particles/emission/particle_surface_flux_init.f90 +++ b/src/particles/emission/particle_surface_flux_init.f90 @@ -122,6 +122,7 @@ SUBROUTINE InitializeParticleSurfaceflux() USE MOD_Particle_Vars ,ONLY: UseCircularInflow, DoForceFreeSurfaceFlux USE MOD_Particle_Sampling_Vars ,ONLY: UseAdaptive USE MOD_Restart_Vars ,ONLY: DoRestart, RestartTime +USE MOD_DSMC_Vars ,ONLY: AmbiPolarSFMapping, DSMC, useDSMC #if USE_MPI USE MOD_Particle_Vars ,ONLY: DoPoissonRounding, DoTimeDepInflow USE MOD_Particle_MPI_Vars ,ONLY: PartMPI @@ -138,7 +139,7 @@ SUBROUTINE InitializeParticleSurfaceflux() !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES ! Local variable declaration -INTEGER :: iSpec,iSF,SideID,BCSideID,iSide,ElemID,iLocSide,iSample,jSample,currentBC +INTEGER :: iSpec,iSF,SideID,BCSideID,iSide,ElemID,iLocSide,iSample,jSample,currentBC, MaxSF, iSFElec INTEGER :: iCopy1, iCopy2, iCopy3, MaxSurfacefluxBCs,nDataBC REAL :: tmp_SubSideDmax(SurfFluxSideSize(1),SurfFluxSideSize(2)) REAL :: tmp_SubSideAreas(SurfFluxSideSize(1),SurfFluxSideSize(2)) @@ -313,6 +314,30 @@ SUBROUTINE InitializeParticleSurfaceflux() END IF DoForceFreeSurfaceFlux = GETLOGICAL('DoForceFreeSurfaceFlux','.FALSE.') +IF (useDSMC) THEN + IF (DSMC%DoAmbipolarDiff) THEN + MaxSF = 0 + DO iSpec = 1,nSpecies + MaxSF = MAX(MaxSF,Species(iSpec)%nSurfacefluxBCs) + END DO + + ALLOCATE(AmbiPolarSFMapping(nSpecies,MaxSF)) + AmbiPolarSFMapping = -1 + DO iSpec = 1,nSpecies + IF (Species(iSpec)%ChargeIC.LE.0.0) CYCLE + DO iSF = 1,Species(iSpec)%nSurfacefluxBCs + DO iSFElec = 1,Species(DSMC%AmbiDiffElecSpec)%nSurfacefluxBCs + IF(Species(iSpec)%Surfaceflux(iSF)%BC.EQ.Species(DSMC%AmbiDiffElecSpec)%Surfaceflux(iSFElec)%BC) THEN + AmbiPolarSFMapping(iSpec,iSF) = iSFElec + END IF + END DO + IF(AmbiPolarSFMapping(iSpec,iSF).EQ.-1) CALL abort(__STAMP__,& + 'ERROR: No corresponding Electron Surface Flux found for Species: ',IntInfoOpt=iSpec) + END DO + END DO + END IF +END IF + END SUBROUTINE InitializeParticleSurfaceflux