diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 44f348e98..8dd2d2dc6 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -752,6 +752,12 @@ NIG_PIC_poisson_Boris-Leapfrog: stage: reggie_nightly script: - cd build ; python ../reggie/reggie.py ../regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/ + +NIG_PIC_poisson_Boris-Leapfrog_PETSC: + <<: *defaults_nightly + stage: reggie_nightly + script: + - cd build ; python ../reggie/reggie.py ../regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/ NIG_PIC_poisson_Leapfrog: <<: *defaults_nightly diff --git a/REGGIE.md b/REGGIE.md index 9d6f54280..11c5e4033 100644 --- a/REGGIE.md +++ b/REGGIE.md @@ -339,12 +339,21 @@ Testing PIC compiled with Leapfrog integration (poisson,Leapfrog), solving Poiss Testing PIC compiled with Boris-Leapfrog integration (poisson,Boris-Leapfrog), solving Poisson's equation: [Link to build](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/builds.ini). -| **No.** | **Case** | **CMAKE-CONFIG** | **Feature** | **Execution** | **Comparing** | **Readme** | -| :-----: | :-------------------------------: | :------------------------------: | :------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------: | :-----------: | :--------------------------------------------------------------------: | :------------------------------------------------------------------------------------------------: | -| 1 | 2D_HET_Liu2010 | CMAKE_BUILD_TYPE = Release,Debug | 2D Poisson-PIC, BGGas distribution, null collision on/off, pre-defined external magnetic field, neutralization BC, SEE model with variable electron bulk temperature, particle flux, total electric current and emitted SEE over time into SurfaceAnalyze.csv | nProcs=3,6,12 | integrate number of electrons impinging the anode (SurfaceAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/2D_HET_Liu2010/readme.md) | -| 2 | 2D_Landmark | CMAKE_BUILD_TYPE = Release,Debug | 2D Poisson-PIC, emission models for Landmark (volumetric ionization and neutralizer) | nProcs=4 | integrate number of electrons impinging the anode (SurfaceAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/2D_Landmark/readme.md) | -| 3 | 3D_HET_Liu2010 | CMAKE_BUILD_TYPE = Release,Debug | 3D Poisson-PIC, BGGas distribution, null collision on/off, pre-defined external magnetic field, neutralization BC, SEE model with variable electron bulk temperature, dielectric surface charging (hollow cylinder), vMPF=T restart from vMPF=F restart file | nProcs=6 | integrate number of electrons impinging the anode (SurfaceAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/3D_HET_Liu2010/readme.md) | -| 4 | MCC_EBeam_SpeciesSpecificTimestep | CMAKE_BUILD_TYPE = Release,Debug | 1D-PIC-MCC electron beam, emission current surface flux and species-specific time step for electrons, using ManualTimeStep for MCC | nProcs=4 | Number density (PartAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_SpeciesSpecificTimestep/readme.md) | +| **No.** | **Case** | **CMAKE-CONFIG** | **Feature** | **Execution** | **Comparing** | **Readme** | +| :-----: | :-------------------------------: | :--------------: | :------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------: | :-----------: | :--------------------------------------------------------------------: | :------------------------------------------------------------------------------------------------: | +| 1 | 2D_HET_Liu2010 | | 2D Poisson-PIC, BGGas distribution, null collision on/off, pre-defined external magnetic field, neutralization BC, SEE model with variable electron bulk temperature, particle flux, total electric current and emitted SEE over time into SurfaceAnalyze.csv | nProcs=3,6,12 | integrate number of electrons impinging the anode (SurfaceAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/2D_HET_Liu2010/readme.md) | +| 2 | 2D_Landmark | | 2D Poisson-PIC, emission models for Landmark (volumetric ionization and neutralizer) | nProcs=4 | integrate number of electrons impinging the anode (SurfaceAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/2D_Landmark/readme.md) | +| 3 | 3D_HET_Liu2010 | | 3D Poisson-PIC, BGGas distribution, null collision on/off, pre-defined external magnetic field, neutralization BC, SEE model with variable electron bulk temperature, dielectric surface charging (hollow cylinder), vMPF=T restart from vMPF=F restart file | nProcs=6 | integrate number of electrons impinging the anode (SurfaceAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/3D_HET_Liu2010/readme.md) | +| 4 | MCC_EBeam_3D_and_2D-axisym | | 1D-PIC-MCC electron beam, emission current surface flux, using ManualTimeStep for MCC, 3D cylinder geometry or 2D axisymmetric HDG | nProcs=4 | Number density (PartAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/readme.md) | +| 5 | MCC_EBeam_SpeciesSpecificTimestep | | 1D-PIC-MCC electron beam, emission current surface flux and species-specific time step for electrons, using ManualTimeStep for MCC | nProcs=4 | Number density (PartAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_SpeciesSpecificTimestep/readme.md) | + +### NIG_PIC_poisson_Boris-Leapfrog_PETSC + +Testing PIC compiled with Boris-Leapfrog integration (poisson,Boris-Leapfrog) and PETSC, solving Poisson's equation: [Link to build](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/builds.ini). + +| **No.** | **Case** | **CMAKE-CONFIG** | **Feature** | **Execution** | **Comparing** | **Readme** | +| :-----: | :-------------------------------: | :--------------: | :------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------: | :-----------: | :--------------------------------------------------------------------: | :------------------------------------------------------------------------------------------------: | +| 1 | 2D_axisymmetricHDG_OConner | | 2D Poisson-PIC with axisymmetric HDG, electron beam expanding over time, testing of different depositions (cell_mean, cell_volweight, cell_volweight_mean), with PETSC (Precond-T) | nProcs=1,2 | kinetic energy (PartAnalyze.csv) and potential energy (FieldAnalyze.csv) of electrons over time | [Link](regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/readme.md) | ### NIG_PIC_poisson_RK3 diff --git a/docs/documentation/userguide/features-and-models/poisson-field-solver.md b/docs/documentation/userguide/features-and-models/poisson-field-solver.md index c3b58d322..d4ab9e531 100644 --- a/docs/documentation/userguide/features-and-models/poisson-field-solver.md +++ b/docs/documentation/userguide/features-and-models/poisson-field-solver.md @@ -28,19 +28,6 @@ To see the available parameter input file options, simply run ./bin/piclas --help HDG -## CG Solver -The default numerical method for solving the resulting system of linear equations, is the Conjugate Gradient Method. The following -input parameter can be set to control the simulation - - epsCG = 1e-6 ! default value for the abort residual - MaxIterCG = 500 ! default value for the number of CG solver iterations - -where `epsCG` is the residual of the CG solver and `MaxIterCG` are the maximum number of iteration performed in one time step to -solve the system. -Furthermore, the residual can be either set absolute (default) or relative via - - useRelativeAbortCrit = F ! default - ## PETSc Solver A multitude of different numerical methods to solve the resulting system of linear equations is given by the implemented PETSc library {cite}`petsc-web-page`, {cite}`petsc-user-ref`, {cite}`petsc-efficient`. For detailed installation steps of PETSc within PICLas, see Section {ref}`sec:petsc-installation`. @@ -61,7 +48,33 @@ where the following options are possible | `3` | iterative | Krylov subspace | | `10` | direct | | -Note that the same parameter setting for `epsCG` will result in a smaller residual with PETSc as compared with the default CG solver -without using the PETSc library. +## CG Solver +The formerly used numerical method for solving the resulting system of linear equations, is the Conjugate Gradient Method. The following +input parameter can be set to control the simulation + + epsCG = 1e-6 ! default value for the abort residual + MaxIterCG = 500 ! default value for the number of CG solver iterations + +where `epsCG` is the residual of the CG solver and `MaxIterCG` are the maximum number of iteration performed in one time step to +solve the system. +Furthermore, the residual can be either set absolute (default) or relative via + + useRelativeAbortCrit = F ! default + +Note that the same parameter setting for `epsCG` will result in a higher residual without the PETSC library as compared with PETSC solver. + +## Symmetric Simulations + +The Poisson solver is suitable for 2D/1D and axially symmetric 2D simulations. Use the following parameter for these types of simulations: + + ! 2D + Particles-Symmetry-Order = 2 + + ! 2D axisymmetrc + Particles-Symmetry-Order = 2 + Particles-Symmetry2DAxisymmetric = true + + ! 1D + Particles-Symmetry-Order = 1 diff --git a/regressioncheck/NIG_PIC_Deposition/Plasma_Ball_cell_volweight_mean_save_CVWM/PartAnalyze_ref_save.csv b/regressioncheck/NIG_PIC_Deposition/Plasma_Ball_cell_volweight_mean_save_CVWM/PartAnalyze_ref_save.csv index c38d96c39..7a98639c9 100644 --- a/regressioncheck/NIG_PIC_Deposition/Plasma_Ball_cell_volweight_mean_save_CVWM/PartAnalyze_ref_save.csv +++ b/regressioncheck/NIG_PIC_Deposition/Plasma_Ball_cell_volweight_mean_save_CVWM/PartAnalyze_ref_save.csv @@ -1,2 +1,3 @@ 001-TIME,002-Charge,003-Charge-absError,004-Charge-relError -0.1000000000000000E-011,0.1069379900794236E+002,0.1369025896277520E-001,0.1281846401056845E-002 +0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000 +0.1000000000000000E-011,0.1068010874898000E+002,0.4121147867408581E-012,0.3858713393533871E-013 diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/PartAnalyze_ref.csv b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/PartAnalyze_ref.csv new file mode 100644 index 000000000..ea25ea444 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/PartAnalyze_ref.csv @@ -0,0 +1,21 @@ +001-TIME,002-NumDens-Spec-001,003-NumDens-Spec-002,004-NumDens-Spec-003,005-NumDens-Spec-004,006-Current-Spec-001-SF-001 +0,0,2E+019,0,2E+019,0 +1E-09,3028100000000000,2E+019,6600000000000,2.00030347E+019,0.000201073154515 +2.00000000000007E-09,3029000000000000,2E+019,13500000000000,2.00030425E+019,0.00020027206625 +2.99999999999997E-09,3029500000000000,2E+019,19300000000000,2.00030488E+019,0.000199470977985 +3.99999999999997E-09,3028700000000000,2E+019,25900000000000,2.00030546E+019,0.00020027206625 +5.00000000000028E-09,3028000000000000,2E+019,33700000000000,2.00030617E+019,0.00020027206625 +6.00000000000059E-09,3028700000000000,2E+019,41500000000000,2.00030702E+019,0.000199470977985 +7.0000000000009E-09,3029100000000000,2E+019,47900000000000,2.0003077E+019,0.00020027206625 +8.00000000000075E-09,3027400000000000,2E+019,54800000000000,2.00030822E+019,0.00020027206625 +9.00000000000023E-09,3027900000000000,2E+019,62700000000000,2.00030906E+019,0.000199470977985 +0.00000001,3027100000000000,2E+019,70000000000000,2.00030971E+019,0.00020027206619146 +0.00000002,3027600000000000,2E+019,133300000000000,2.00031609E+019,0.000200272066896514 +0.00000003,3027300000000000,2E+019,197700000000000,2.0003225E+019,0.000199470980245914 +0.00000004,3025900000000000,2E+019,265700000000000,2.00032916E+019,0.000200272068521319 +0.00000005,3024000000000000,2E+019,326300000000000,2.00033503E+019,0.000200272068521319 +0.00000006,3024200000000000,2E+019,392200000000000,2.00034164E+019,0.000200272068518668 +0.00000007,3023700000000000,2E+019,457800000000000,2.00034815E+019,0.000200272068521319 +0.00000008,3022200000000000,2E+019,523400000000000,2.00035456E+019,0.000200272068521319 +0.00000009,3022000000000000,2E+019,585700000000000,2.00036077E+019,0.000199470980247234 +0.0000001,3021400000000000,2E+019,655700000000000,2.00036771E+019,0.000200272068521319 diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/analyze.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/analyze.ini new file mode 100644 index 000000000..fd247575c --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/analyze.ini @@ -0,0 +1,5 @@ +! compare the last line of PartAnalyze.csv with a reference file +compare_data_file_name = PartAnalyze.csv +compare_data_file_reference = PartAnalyze_ref.csv +compare_data_file_tolerance = 15e-2 +compare_data_file_tolerance_type = relative diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/command_line.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/command_line.ini new file mode 100644 index 000000000..f0056b1d6 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/command_line.ini @@ -0,0 +1,2 @@ +MPI=4 +database = ../../../SpeciesDatabase.h5 diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/externals.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/externals.ini new file mode 100644 index 000000000..b83d17c0d --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/externals.ini @@ -0,0 +1,6 @@ +! --- Externals Tool Reggie +MPI = 1 , 1 ! Single execution +externalbinary = ./hopr/build/bin/hopr , ./hopr/build/bin/hopr ! Relative binary path in build directory +externaldirectory = pre-hopr-3D-cylinder , pre-hopr-2D-axisym ! Directory name, where the files are located for the external tool reggie +externalruntime = pre , pre ! Run after piclas is completed (post: after, pre: before) +nocrosscombination:MPI,externalbinary,externaldirectory,externalruntime diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/parameter.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/parameter.ini new file mode 100644 index 000000000..faf698d99 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/parameter.ini @@ -0,0 +1,149 @@ +! =============================================================================== ! +! EQUATION (linearscalaradvection) +! =============================================================================== ! +IniExactFunc = 0 +! =============================================================================== ! +! DISCRETIZATION +! =============================================================================== ! +N = 1 ! Polynomial degree +NAnalyze = 1 ! Number of analyze points +NVisu = 1 +VisuParticles = T +TimeStampLength = 15 +! =============================================================================== ! +! MESH +! =============================================================================== ! +MeshFile = ./pre-hopr-3D-cylinder/Cylinder_mesh.h5,./pre-hopr-2D-axisym/2D_axisym_mortonZ_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 +CalcSurfFluxInfo = T +CalcNumDens = T +! =============================================================================== ! +! Load Balance +! =============================================================================== ! +DoLoadBalance = T +Load-DeviationThreshold = 1e-2 +LoadBalanceMaxSteps = 10 +! =============================================================================== ! +! CALCULATION +! =============================================================================== ! +tend = 1.0E-7 +Analyze_dt = 1.0E-7 +! =============================================================================== ! +! FIELD SOLVER +! =============================================================================== ! +epsCG = 1e-4 ! 1e-6 +MaxIterCG = 2000 +PrecondType = 2 + +PIC-Deposition-Type = cell_volweight_mean +! =============================================================================== ! +! FIELD BOUNDARY +! =============================================================================== ! +BoundaryName = BC_Xminus +BoundaryType = (/5,1/) +RefState = (/-20000.0 , 0.0 , 0.0/) + +BoundaryName = BC_Xplus +BoundaryType = (/4,0/) + +! =============================================================================== ! +! PARTICLE BOUNDARY +! =============================================================================== ! +Part-nBounds=6 +Part-Boundary1-SourceName = BC_Xplus ! Cathode: -20 kV +Part-Boundary1-Condition = open + +Part-Boundary2-SourceName = BC_Xminus ! Anode: 0V +Part-Boundary2-Condition = reflective + +Part-Boundary3-SourceName = BC_Yplus +Part-Boundary3-Condition = symmetric + +Part-Boundary4-SourceName = BC_Yminus +Part-Boundary4-Condition = symmetric,symmetric_axis + +Part-Boundary5-SourceName = BC_Zplus +Part-Boundary5-Condition = symmetric,symmetric_dim + +Part-Boundary6-SourceName = BC_Zminus +Part-Boundary6-Condition = symmetric,symmetric_dim + +Part-FIBGMdeltas=(/1e-4,1e-4,1e-4/) +! =============================================================================== ! +! PARTICLES +! =============================================================================== ! +Particles-Species-Database = SpeciesDatabase.h5 +Particles-CollXSec-NullCollision = T + +Part-Species1-UseCollXSec = T + +Part-maxParticleNumber=50000 +Part-nSpecies = 3 +Part-Species$-MacroParticleFactor = 200.0, 300.0 + +Part-AnalyzeStep = 5 +IterDisplayStep = 100 +ManualTimeStep = 10e-11 + +Particles-Symmetry-Order = 3,2 +Particles-Symmetry2DAxisymmetric = F,T +! =============================================================================== ! +! Species1 - electron +! =============================================================================== ! +Part-Species1-SpeciesName = electron +Part-Species1-nSurfaceFluxBCs=1 + +Part-Species1-Surfaceflux1-BC=2 +Part-Species1-Surfaceflux1-VeloIC = 0. +Part-Species1-Surfaceflux1-VeloVecIC = (/0,0,1/),(/1,0,0/) +Part-Species1-Surfaceflux1-velocityDistribution = maxwell_lpn +Part-Species1-Surfaceflux1-MWTemperatureIC = 5. +Part-Species1-Surfaceflux1-EmissionCurrent = 2E-4 + +nocrosscombination:MeshFile,Part-Species$-MacroParticleFactor,Particles-Symmetry-Order,Particles-Symmetry2DAxisymmetric,Part-Species1-Surfaceflux1-VeloVecIC,Part-Boundary4-Condition,Part-Boundary5-Condition,Part-Boundary6-Condition +! =============================================================================== ! +! 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 = 2E19 +Part-Species2-Init1-VeloIC = 0 +Part-Species2-Init1-VeloVecIC = (/0.,0.,1./) +! =============================================================================== ! +! Species3 - ArgonIon +! =============================================================================== ! +Part-Species3-SpeciesName = ArIon1 +! =============================================================================== ! +! DSMC +! =============================================================================== ! +Particles-HaloEpsVelo=5.0E+08 +Particles-DSMC-CalcSurfaceVal=F +UseDSMC=true +Particles-DSMC-CollisMode = 3 +Part-NumberOfRandomSeeds=2 +Particles-RandomSeed1=1 +Particles-RandomSeed2=2 +Particles-DSMC-UseOctree=F +Particles-DSMC-UseNearestNeighbour = F +Particles-DSMC-CalcQualityFactors = F +! =============================================================================== ! +! Reactions +! =============================================================================== ! +DSMC-NumOfReactions = 1 +! ---------------------------------------------------- +! Electron impact +! ---------------------------------------------------- +! Ionization: Ar + e --> ArIon1 + e + e +DSMC-Reaction1-ReactionModel = XSec +DSMC-Reaction1-Reactants = (/2,1,0/) +DSMC-Reaction1-Products = (/3,1,1,0/) diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/pre-hopr-2D-axisym/hopr.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/pre-hopr-2D-axisym/hopr.ini new file mode 100644 index 000000000..2467ab7c0 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/pre-hopr-2D-axisym/hopr.ini @@ -0,0 +1,67 @@ +!==================================================================================================================================================== +! Variables +!==================================================================================================================================================== +DEFVAR = (REAL): f1 = 0.0 +DEFVAR = (REAL): f2 = 0.0 +DEFVAR = (REAL): f3 = 0.0 + +DEFVAR=(INT): ix = 50 +DEFVAR=(INT): iy1 = 4 +DEFVAR=(INT): iz = 1 + +DEFVAR = (REAL): x0 = 0.0 +DEFVAR = (REAL): x1 = 5.0 + +DEFVAR = (REAL): y0 = 0.0 +DEFVAR = (REAL): y1 = 0.05 + +DEFVAR = (REAL): z0 = -0.00005 ! changed from 0.0 due to axisym requirement +DEFVAR = (REAL): z1 = 0.00005 ! changed from 0.01 due to axisym requirement +!==================================================================================================================================================== +! General +!==================================================================================================================================================== +ProjectName = 2D_axisym_mortonZ +Debugvisu = T +DebugVisuLevel = 1 +NVisu = 1 +Mode = 1 +postscalemesh = T +meshscale = 1e-3 +!jacobiantolerance = 1e-20 +useCurveds = F +SpaceQuandt = 100.0 + +sfc_type = mortonZ + +!==================================================================================================================================================== +! Boundary Conditions +!==================================================================================================================================================== +BoundaryName = BC_Xminus +BoundaryType = (/3,0,0,0/) + +BoundaryName = BC_Xplus +BoundaryType = (/3,0,0,0/) + +BoundaryName = BC_Yminus +BoundaryType = (/10,0,0,0/) + +BoundaryName = BC_Yplus +BoundaryType = (/10,0,0,0/) + +BoundaryName = BC_Zminus +BoundaryType = (/10,0,0,0/) + +BoundaryName = BC_Zplus +BoundaryType = (/10,0,0,0/) + +nZones = 1 + +!==================================================================================================================================================== +! Section 1 of 3: y0 to y1 +!==================================================================================================================================================== +Corner =(/x0,y0,z0,, x1,y0,z0,, x1,y1,z0,, x0,y1,z0,, x0,y0,z1,, x1,y0,z1,, x1,y1,z1,, x0,y1,z1 /) +nElems =(/ix,iy1,iz/) +factor =(/f1,f2,f3/) +BCIndex =(/5 ,3 ,2 ,4 ,1 ,6/) ! +! =(/z-,y-,x+,y+,x-,z+/) ! Indices of Boundary Conditions +elemtype =108 diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/pre-hopr-3D-cylinder/hopr.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/pre-hopr-3D-cylinder/hopr.ini new file mode 100644 index 000000000..758d96367 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/pre-hopr-3D-cylinder/hopr.ini @@ -0,0 +1,141 @@ +DEFVAR=(INT): icenter = 2 ! no. elems in inner square i0xi0 +DEFVAR=(INT): i01 = 2 ! Number of elements in azimuthal direction i.e., the number of elements per 45° of the cylinder.\n ! The total number will result in 2*i01 (quarter cylinder), 4*i01 (half cylinder) or 8*i01 (full cylinder) for the total number of elements in azimuthal direction +DEFVAR=(INT): ir1 = 1 ! Number of elements in radial direction +DEFVAR=(INT): iz0 = 50 ! Number of elements in z-direction + +DEFVAR=(REAL): r01 = 2.0 ! middle square dim +DEFVAR=(REAL): r02 = 4.0 ! middle square dim + +DEFVAR=(REAL): z0 = 0 ! half length of domain in z +DEFVAR=(REAL): z1 = 5e-3 ! half length of domain in z + +DEFVAR=(REAL): f1 = 1.0 ! stretching factor in radial direction (a larger value than 1.0 will create small elements at the inner cylinder) + +!================================================================================================================================= ! +! OUTPUT +!================================================================================================================================= ! +ProjectName = Cylinder +Debugvisu = T ! Visualiz1e mesh and boundary conditions (tecplot ascii) +checkElemJacobians = T + +!================================================================================================================================= ! +! MESH +!================================================================================================================================= ! +Mode = 1 ! Mode for Cartesian boxes +nZones = 9 ! number of boxes + +useCurveds = F +BoundaryOrder = 3 ! NGeo+1 + +MeshPostDeform = 1 ! deforms [-1,1]^2 to a cylinder with radius Postdeform_R0 +PostDeform_R0 = 0.01250e-3 ! here domain is [-4,4]^2 mapped to a cylinder with radius 0.25*4 = 1 + +!================================================================================================================================= ! +! BOUNDARY CONDITIONS +!================================================================================================================================= ! +BoundaryName = BC_Xminus ! BC index 1 +BoundaryType = (/3,0,0,0/) ! (/ Type, curveIndex, State, alpha /) + +BoundaryName = BC_Xplus ! BC index 2 +BoundaryType = (/3,0,0,0/) ! (/ Type, curveIndex, State, alpha /) + +BoundaryName = BC_Yminus ! BC index 3 +BoundaryType = (/10,0,0,0/) ! (/ Type, curveIndex, State, alpha /) + +BoundaryName = BC_Yplus ! BC index 4 +BoundaryType = (/10,0,0,0/) ! (/ Type, curveIndex, State, alpha /) + +BoundaryName = BC_Zminus ! BC index 5 +BoundaryType = (/10,0,0,0/) ! (/ Type, curveIndex, State, alpha /) + +BoundaryName = BC_Zplus ! BC index 6 +BoundaryType = (/10,0,0,0/) ! (/ Type, curveIndex, State, alpha /) + +! -------------------------------------------------------------------------------------------------------------------------------- ! +! MESH +! -------------------------------------------------------------------------------------------------------------------------------- ! +Mode = 1 ! Mode for Cartesian boxes +nZones = 13 ! number of boxes +!center +Corner = (/-r01 , -r01 , z0 , , r01 , -r01 , z0 , , r01 , r01 , z0 , , -r01 , r01 , z0 , , -r01 , -r01 , z1 , , r01 , -r01 , z1 , , r01 , r01 , z1 , , -r01 , r01 , z1 /) +nElems = (/icenter,icenter,iz0/) ! number of elements in each direction +!nElems = (/i0,i0,i0/) ! number of elements in each direction +BCIndex = (/1 , 0 , 0 , 0 , 0 , 2/) ! Indices of Boundary Conditions +! = (/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype = 108 ! element type (108: Hexahedral) + + + +! --------------------------------------------------------------- +! Upper cylinder half +! --------------------------------------------------------------- + +!left-lower (x-) +Corner =(/-r01 , 0. ,z0 ,, -r02 , 0. ,z0 ,, -r02 , r02 , z0 ,, -r01 , r01 , z0 ,, -r01 , 0. , z1 ,, -r02 , 0. , z1 ,, -r02 , r02 , z1 ,, -r01 , r01 , z1 /) +nElems =(/ir1,i01,iz0/) ! number of elements in each direction +BCIndex =(/1 , 0 , 5 , 0 , 0 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) +! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype =108 ! element type (108: Hexahedral) +factor =(/f1,1.,1./) ! stretching + +!left-upper (y+) +Corner =(/0. , r01 , z0,, -r01 , r01 , z0 ,, -r02 , r02 , z0 ,, 0. , r02 , z0 ,, 0. , r01 , z1 ,, -r01 , r01 , z1 ,, -r02 , r02 , z1 ,, 0. , r02 , z1 /) +nElems =(/i01,ir1,iz0/) ! number of elements in each direction +BCIndex =(/1 , 0 , 0 , 4 , 0 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) +! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype =108 ! element type (108: Hexahedral) +factor =(/1.,f1,1./) ! stretching + +!right-lower (x+) +Corner =(/r01 , 0. , z0 ,, r02 , 0. , z0 ,, r02 , r02 , z0 ,, r01 , r01 , z0 ,, r01 , 0. , z1 ,, r02 , 0. , z1 ,, r02 , r02 , z1 ,, r01 , r01 , z1 /) +nElems =(/ir1,i01,iz0/) ! number of elements in each direction +BCIndex =(/1 , 0 , 6 , 0 , 0 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) + +! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype =108 ! element type (108: Hexahedral) +factor =(/f1,1.,1./) ! stretching + +!right-upper (y+) +Corner =(/0. , r01 , z0 ,, r01 , r01 , z0 ,, r02 , r02 , z0 ,, 0. , r02 , z0 ,, 0. , r01 , z1 ,, r01 , r01 , z1 ,, r02 , r02 , z1 ,, 0. , r02 , z1 /) +nElems =(/i01,ir1,iz0/) ! number of elements in each direction +BCIndex =(/1 , 0 , 0 , 4 , 0 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) +! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype =108 ! element type (108: Hexahedral) +factor =(/1.,f1,1./) ! stretching + +! --------------------------------------------------------------- +! Bottom cylinder half +! --------------------------------------------------------------- +!left-lower (x-) +Corner =(/-r01 , 0. , z0 ,, -r02 , 0. , z0 ,, -r02 , -r02 , z0 ,, -r01 , -r01 , z0 ,, -r01 , 0. , z1 ,, -r02 , 0. , z1 ,, -r02 , -r02 , z1 ,, -r01 , -r01 , z1 /) +nElems =(/ir1,i01,iz0/) ! number of elements in each direction +BCIndex =(/1 , 0 , 5 , 0 , 0 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) +! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype =108 ! element type (108: Hexahedral) +factor =(/f1,1.,1./) ! stretching + +!left-upper (y+) +Corner =(/0. , -r01 , z0 ,, -r01 , -r01 , z0 ,, -r02 , -r02 , z0 ,, 0. , -r02 , z0 ,, 0. , -r01 , z1 ,, -r01 , -r01 , z1 ,, -r02 , -r02 , z1 ,, 0. , -r02 , z1 /) +nElems =(/i01,ir1,iz0/) ! number of elements in each direction +BCIndex =(/1 , 0 , 0 , 3 , 0 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) +! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype =108 ! element type (108: Hexahedral) +factor =(/1.,f1,1./) ! stretching + +!right-lower (x+) +Corner =(/r01 , 0. , z0 ,, r02 , 0. , z0 ,, r02 , -r02 , z0 ,, r01 , -r01 , z0 ,, r01 , 0. , z1 ,, r02 , 0. , z1 ,, r02 , -r02 , z1 ,, r01 , -r01 , z1 /) +nElems =(/ir1,i01,iz0/) ! number of elements in each direction +BCIndex =(/1 , 0 , 6 , 0 , 0 , 2/) ! Indices of Boundary Conditions for six Boundary Faces (z- , y- , x+ , y+ , x- , z+) +! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype =108 ! element type (108: Hexahedral) +factor =(/f1,1.,1./) ! stretching + +!right-upper (y+) +Corner =(/0. , -r01 , z0 ,, r01 , -r01 ,z0 ,, r02 , -r02 , z0 ,, 0. , -r02 , z0 ,, 0. , -r01 , z1 ,, r01 , -r01 , z1 ,, r02 , -r02 , z1 ,, 0. , -r02 , z1 /) +nElems =(/i01,ir1,iz0/) ! number of elements in each direction +BCIndex =(/1 , 0 , 0 , 3 , 0 , 2/) ! Indices of Boundary Conditions +! =(/z- , y- , x+ , y+ , x- , z+/) ! Indices of Boundary Conditions +elemtype =108 ! element type (108: Hexahedral) +factor =(/1.,f1,1./) ! stretching + + diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/readme.md b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/readme.md new file mode 100644 index 000000000..ac35556ce --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog/MCC_EBeam_3D_and_2D-axisym/readme.md @@ -0,0 +1,8 @@ +# PIC-MCC - Electron beam acceleration and ionization of background gas +* Inserting electrons with a zero velocity and at a low temperature +* Acceleration with a potential difference of 20 kV +* Ionization of background gas through MCC +* ATTENTION: two different meshes are tested + * 3D cylinder ... with === axial direction in z === + * 2D axisymmetric with === axial direction in x === +* Comparing the number density of the electron and ions in the channel diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/FieldAnalyze_ref.csv b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/FieldAnalyze_ref.csv new file mode 100644 index 000000000..e0ebcdf03 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/FieldAnalyze_ref.csv @@ -0,0 +1,9 @@ +001-time,002-E-El +0,0 +5E-10,3.58933445054784E-09 +0.000000001,9.75543419771557E-09 +0.0000000015,1.54323429703684E-08 +0.000000002,1.96271838047281E-08 +0.0000000025,1.98073944003939E-08 +0.000000003,1.98238640067753E-08 +0.0000000035,1.99446790647449E-08 diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/PartAnalyze_ref.csv b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/PartAnalyze_ref.csv new file mode 100644 index 000000000..5fb2ae597 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/PartAnalyze_ref.csv @@ -0,0 +1,9 @@ +001-TIME,002-Ekin-001 +0.0000000000000000E+000,0.0000000000000000E+000 +0.5000000000000000E-009,0.9085364027978898E-006 +0.1000000000000000E-008,0.1816794672013311E-005 +0.1500000000000000E-008,0.2710430761980031E-005 +0.2000000000000000E-008,0.3559750823868184E-005 +0.2500000000000000E-008,0.3598174446363344E-005 +0.3000000000000000E-008,0.3599385605377411E-005 +0.3500000000000000E-008,0.3610420487270387E-005 diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/analyze.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/analyze.ini new file mode 100644 index 000000000..4ad9116f8 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/analyze.ini @@ -0,0 +1,9 @@ +! compare the last line of PartAnalyze.csv with a reference file +compare_column_file = PartAnalyze.csv , FieldAnalyze.csv +compare_column_reference_file = PartAnalyze_ref.csv , FieldAnalyze_ref.csv +compare_column_tolerance_value = 7.2E-8 , 4E-10 +compare_column_tolerance_type = absolute , absolute +compare_column_index = 1 + +compare_column_one_diff_per_run=F +!nocrosscombination:compare_column_file,compare_column_reference_file,compare_column_tolerance_value,compare_column_tolerance_type,compare_column_index diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/command_line.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/command_line.ini new file mode 100644 index 000000000..4b20f5225 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/command_line.ini @@ -0,0 +1 @@ +MPI=1,2 diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/externals.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/externals.ini new file mode 100644 index 000000000..a3ff6aa61 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/externals.ini @@ -0,0 +1,6 @@ +! --- Externals Tool Reggie +MPI = 1 ! Single execution +externalbinary = ./hopr/build/bin/hopr ! Relative binary path in build directory +externaldirectory = hopr.ini ! Directory name, where the files are located for the external tool reggie +externalruntime = pre ! Run after piclas is completed (post: after, pre: before) +nocrosscombination:MPI,externalbinary,externaldirectory,externalruntime diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/hopr.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/hopr.ini new file mode 100644 index 000000000..50d9ba11a --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/hopr.ini @@ -0,0 +1,67 @@ +!==================================================================================================================================================== +! Variables +!==================================================================================================================================================== +DEFVAR = (REAL): f1 = 0.0 +DEFVAR = (REAL): f2 = 0.0 +DEFVAR = (REAL): f3 = 0.0 + +DEFVAR=(INT): ix = 80 +DEFVAR=(INT): iy1 = 16 +DEFVAR=(INT): iz = 1 + +DEFVAR = (REAL): x0 = 0.0 +DEFVAR = (REAL): x1 = 0.1 + +DEFVAR = (REAL): y0 = 0.0 +DEFVAR = (REAL): y1 = 0.02 + +DEFVAR = (REAL): z0 = -0.0005 +DEFVAR = (REAL): z1 = 0.0005 +!==================================================================================================================================================== +! General +!==================================================================================================================================================== +ProjectName = Cylinder_Axisym_mortonZ +Debugvisu = T +DebugVisuLevel = 1 +NVisu = 1 +Mode = 1 +postscalemesh = T +!meshscale = 1e-3 +!jacobiantolerance = 1e-20 +useCurveds = F +SpaceQuandt = 100.0 + +sfc_type = mortonZ + +!==================================================================================================================================================== +! Boundary Conditions +!==================================================================================================================================================== +BoundaryName = BC_Inflow +BoundaryType = (/3,0,0,0/) + +BoundaryName = BC_Outflow +BoundaryType = (/3,0,0,0/) + +BoundaryName = BC_Yminus +BoundaryType = (/10,0,0,0/) + +BoundaryName = BC_Wall +BoundaryType = (/4,0,0,0/) + +BoundaryName = BC_Zminus +BoundaryType = (/10,0,0,0/) + +BoundaryName = BC_Zplus +BoundaryType = (/10,0,0,0/) + +nZones = 1 + +!==================================================================================================================================================== +! Section 1 of 3: y0 to y1 +!==================================================================================================================================================== +Corner =(/x0,y0,z0,, x1,y0,z0,, x1,y1,z0,, x0,y1,z0,, x0,y0,z1,, x1,y0,z1,, x1,y1,z1,, x0,y1,z1 /) +nElems =(/ix,iy1,iz/) +factor =(/f1,f2,f3/) +BCIndex =(/5 ,3 ,2 ,4 ,1 ,6/) ! +! =(/z-,y-,x+,y+,x-,z+/) ! Indices of Boundary Conditions +elemtype =108 diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/parameter.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/parameter.ini new file mode 100644 index 000000000..bbf626dd7 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/parameter.ini @@ -0,0 +1,137 @@ +! =============================================================================== ! +! EQUATION +! =============================================================================== ! +IniExactFunc = 0 +! =============================================================================== ! +! DISCRETIZATION +! =============================================================================== ! +N = 2 +PrecondType = 10 +GeometricNGeo = 6 ! Degree of mesh representation +NAnalyze = 8 ! Number of analyze points +! =============================================================================== ! +! MESH +! =============================================================================== ! +MeshFile = Cylinder_Axisym_mortonZ_mesh.h5 +useCurveds = F +! =============================================================================== ! +! OUTPUT / VISUALIZATION +! =============================================================================== ! +TrackingMethod = TriaTracking +ProjectName = ExpandingParticleBeam +NVisu = 2 +Logging = F +WriteErrorFiles = F +printRandomSeeds = F +DoCalcErrorNorms = F +VisuParticles = T +! =============================================================================== ! +! CALCULATION +! =============================================================================== ! +tend = 3.5e-9 +Analyze_dt = 0.5e-9 +ManualTimeStep = 5e-12 +IterDisplayStep = 50 + +!CFLscale = 0.9 ! Scaling of theoretical CFL number +!c_corr = 1 +! =============================================================================== ! +! Load Balance +! =============================================================================== ! +DoLoadBalance = F +PartWeightLoadBalance = T +!Load-DeviationThreshold = 1.e-2 +Particles-MPIWeight = 0.01 +LoadBalanceMaxSteps = 5 + +! =============================================================================== ! +! Field Boundaries +! =============================================================================== ! +BoundaryName = BC_Inflow +BoundaryType = (/4,0/) ! 4: Dirichlet with zero potential + +BoundaryName = BC_Outflow +BoundaryType = (/4,0/) ! 4: Dirichlet with zero potential + +BoundaryName = BC_Yminus +BoundaryType = (/10,0/) + +BoundaryName = BC_Wall +BoundaryType = (/5,1/) ! 5: Dirichlet, 1: Nbr of RefState +RefState = (/0.0, 0.0, 0.0/) ! RefState Nbr 1: Voltage and Frequency + +BoundaryName = BC_Zminus +BoundaryType = (/10,0/) + +BoundaryName = BC_Zplus +BoundaryType = (/10,0/) +! =============================================================================== ! +! PARTICLES +! =============================================================================== ! +Part-nSpecies = 1 + +Part-nBounds = 6 +Part-Boundary1-SourceName = BC_Inflow +Part-Boundary1-Condition = open +Part-Boundary2-SourceName = BC_Outflow +Part-Boundary2-Condition = open +Part-Boundary3-SourceName = BC_Yminus +Part-Boundary3-Condition = symmetric_axis + +Part-Boundary4-SourceName = BC_Wall +Part-Boundary4-Condition = open + +Part-Boundary5-SourceName = BC_Zminus +Part-Boundary5-Condition = symmetric_dim + +Part-Boundary6-SourceName = BC_Zplus +Part-Boundary6-Condition = symmetric_dim + +PIC-DoInterpolation = T +PIC-Interpolation-Type = particle_position +Part-LorentzType = 3 + +PIC-DoDeposition = T +PIC-Deposition-Type = cell_mean, cell_volweight, cell_volweight_mean +!PIC-Deposition-Type = shape_function_cc +!PIC-shapefunction-radius = 5e-3 + +Part-FIBGMdeltas = (/0.1 , 0.02 , 0.001/) +Part-FactorFIBGM = (/40 , 8 , 1/) + +Particles-Symmetry-Order = 2 +Particles-Symmetry2DAxisymmetric = T + +! =============================================================================== ! +! Species1 - electrons +! =============================================================================== ! +Part-Species1-MassIC = 9.1093826E-31 +Part-Species1-ChargeIC = -1.60217653E-19 +Part-Species1-MacroParticleFactor = 5E5 + +Part-Species1-nSurfaceFluxBCs = 1 +Part-Species1-Surfaceflux1-BC = 1 +Part-Species1-Surfaceflux1-VeloIC = 5e7 +Part-Species1-Surfaceflux1-VeloVecIC = (/1,0,0/) +Part-Species1-Surfaceflux1-velocityDistribution = maxwell_lpn +Part-Species1-Surfaceflux1-MWTemperatureIC = 1. +Part-Species1-Surfaceflux1-EmissionCurrent = 0.25 +Part-Species1-Surfaceflux1-CircularInflow = TRUE +Part-Species1-Surfaceflux1-axialDir = 1 +Part-Species1-Surfaceflux1-origin = (/0.,0./) +Part-Species1-Surfaceflux1-rmax = 8e-3 +! =============================================================================== ! +! Analysis +! =============================================================================== ! +Part-AnalyzeStep = 701!1 +Field-AnalyzeStep = 701!1 +CalcKineticEnergy = T +CalcPartBalance = F +CalcCharge = F +CalcNumSpec = F +CalcPotentialEnergy = T +Part-NumberOfRandomSeeds = 2 +Particles-RandomSeed1 = 1180520427 +Particles-RandomSeed2 = 1708457652 +PIC-OutputSource = T ! HDF5 output of maxwell source terms + diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/readme.md b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/readme.md new file mode 100644 index 000000000..dcaa58591 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/2D_axisymmetricHDG_OConner/readme.md @@ -0,0 +1,10 @@ +# PIC-axisymmetric HDG - O'Connor Benchmark Expanding electron beam +* A beam of electrons is injected into a cylindrical drift tube and evolved over time +* Taken from S.O'Connor et al., A Set of Benchmark Tests for Validation of 3D Particle In Cell Methods (https://doi.org/10.48550/arXiv.2101.09299) +* Beam is expanding over time due to the mutual repulsion (analytical values for beam radius over time are given +* Axisymmetric depositions of HDG are tested using +* Three different depositions are tested (all with PETSc and a combination of PrecondType = 2, 10) + * cell_mean + * cell_volweight + * cell_volweight_mean +* Comparing the kinetic (PartAnalyze.csv) and potential energy (FieldAnalyze.csv) of the electrons over time diff --git a/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/builds.ini b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/builds.ini new file mode 100644 index 000000000..14ab052d1 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Boris-Leapfrog_PETSC/builds.ini @@ -0,0 +1,12 @@ +! relative binary path in build directory +binary=./bin/piclas + +! compiler flags +CMAKE_BUILD_TYPE = Release,Debug +LIBS_BUILD_HDF5 = OFF +PICLAS_POLYNOMIAL_DEGREE = N +PICLAS_EQNSYSNAME = poisson +PICLAS_TIMEDISCMETHOD = Boris-Leapfrog +LIBS_USE_MPI = ON +PICLAS_NODETYPE = GAUSS +PICLAS_PETSC = ON diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8867e202b..1ece73da8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -433,6 +433,7 @@ FILE(GLOB_RECURSE piclasF90 ./src/init/*.f90 ./src/recordpoints/*.f90 ./src/restart/*.f90 ./src/utils/*.f90 + ./src/symmetry/*.f90 ./unitTests/unittest.f90 ./unitTests/unittest_vars.f90) diff --git a/src/equations/poisson/equation.f90 b/src/equations/poisson/equation.f90 index 8a47c2434..e6d073bea 100644 --- a/src/equations/poisson/equation.f90 +++ b/src/equations/poisson/equation.f90 @@ -995,12 +995,13 @@ PPURE SUBROUTINE CalcSourceHDG(i,j,k,iElem,resu, Phi, warning_linear, warning_li !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES REAL :: xvec(3) -REAL :: r1,r2,dx +REAL :: r1,r2 REAL,DIMENSION(3) :: dx1,dx2,dr1dx,dr2dx,dr1dx2,dr2dx2 #ifdef PARTICLES REAL :: source_e ! Electron charge density for Boltzmann relation (electrons as isothermal fluid!) INTEGER :: RegionID #if defined(CODE_ANALYZE) +REAL :: dx REAL :: ElemCharLengthX #endif /*defined(CODE_ANALYZE)*/ #endif /*PARTICLES*/ diff --git a/src/init/define_parameters_init.f90 b/src/init/define_parameters_init.f90 index 0809bf00b..a8f81040e 100644 --- a/src/init/define_parameters_init.f90 +++ b/src/init/define_parameters_init.f90 @@ -32,6 +32,7 @@ SUBROUTINE InitDefineParameters() USE MOD_Globals ,ONLY: MPIRoot USE MOD_MPI_Shared ,ONLY: DefineParametersMPIShared #endif /*USE_MPI*/ +USE MOD_Symmetry ,ONLY: DefineParametersSymmetry USE MOD_Globals_Init ,ONLY: DefineParametersGlobals USE MOD_ReadInTools ,ONLY: prms USE MOD_MPI ,ONLY: DefineParametersMPI @@ -109,6 +110,7 @@ SUBROUTINE InitDefineParameters() #endif /*USE_MPI*/ CALL DefineParametersIO() CALL DefineParametersGlobals() +CALL DefineParametersSymmetry() CALL DefineParametersLoadBalance() CALL DefineParametersInterpolation() CALL DefineParametersRestart() diff --git a/src/init/piclas_init.f90 b/src/init/piclas_init.f90 index cbf95e347..89319922c 100644 --- a/src/init/piclas_init.f90 +++ b/src/init/piclas_init.f90 @@ -85,9 +85,10 @@ SUBROUTINE InitPiclas(IsLoadBalance) #if USE_MPI USE MOD_MPI ,ONLY: InitMPIvars #endif /*USE_MPI*/ +USE MOD_Symmetry ,ONLY: InitSymmetry #ifdef PARTICLES USE MOD_DSMC_Vars ,ONLY: UseDSMC -USE MOD_ParticleInit ,ONLY: InitParticleGlobals,InitParticles,InitSymmetry +USE MOD_ParticleInit ,ONLY: InitParticleGlobals,InitParticles USE MOD_TTMInit ,ONLY: InitTTM,InitIMD_TTM_Coupling USE MOD_TTM_Vars ,ONLY: DoImportTTMFile USE MOD_Particle_Analyze ,ONLY: InitParticleAnalyze @@ -135,11 +136,10 @@ SUBROUTINE InitPiclas(IsLoadBalance) #ifdef PARTICLES ! DSMC handling: useDSMC=GETLOGICAL('UseDSMC') +#endif /*PARTICLES*/ CALL InitSymmetry() -#endif /*PARTICLES*/ - ! Initialization IF(IsLoadBalance)THEN DoRestart=.TRUE. diff --git a/src/interpolation/eval_xyz.f90 b/src/interpolation/eval_xyz.f90 index e7c91b51d..30316a951 100644 --- a/src/interpolation/eval_xyz.f90 +++ b/src/interpolation/eval_xyz.f90 @@ -464,7 +464,7 @@ SUBROUTINE RefElemNewton(Xi,X_In,wBaryCL_N_In,XiCL_N_In,XCL_N_In,dXCL_N_In,N_In, USE MOD_Particle_Vars ,ONLY: PartIsImplicit #endif USE MOD_Particle_Vars ,ONLY: LastPartPos -USE MOD_TimeDisc_Vars ,ONLY: iter,time +USE MOD_TimeDisc_Vars ,ONLY: iter !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT VARIABLES @@ -538,25 +538,7 @@ SUBROUTINE RefElemNewton(Xi,X_In,wBaryCL_N_In,XiCL_N_In,XCL_N_In,dXCL_N_In,N_In, ! Compute inverse of Jacobian sdetJac=getDet(Jac) IF(sdetJac.GT.0.) THEN - sdetJac=1./sdetJac - ELSE - ! Newton has not converged !?!? - IF(Mode.EQ.1)THEN - IPWRITE(UNIT_stdOut,*) ' Particle not inside of element!' - IPWRITE(UNIT_stdOut,*) ' time ', time - IPWRITE(UNIT_stdOut,*) ' Timestep-Iter', iter - IPWRITE(UNIT_stdOut,*) ' sdetJac ', sdetJac - IPWRITE(UNIT_stdOut,*) ' Newton-Iter ', NewtonIter - IPWRITE(UNIT_stdOut,*) ' xi ', xi(1:3) - IPWRITE(UNIT_stdOut,*) ' PartPos ', X_in - IPWRITE(UNIT_stdOut,*) ' GlobalElemID ', ElemID - CALL abort(__STAMP__,'Newton in FindXiForPartPos singular. iter,sdetJac',NewtonIter,sdetJac) - ELSE - Xi(1)=HUGE(1.0) - Xi(2)=Xi(1) - Xi(3)=Xi(1) - RETURN - END IF + sdetJac=1./sdetJac ENDIF sJac=getInv(Jac,sdetJac) @@ -694,6 +676,7 @@ PPURE SUBROUTINE GetRefNewtonStartValue(X_in,Xi,ElemID) #else USE MOD_Mesh_Vars, ONLY:ElemBaryNGeo #endif +USE MOD_Symmetry_Vars ,ONLY: Symmetry !----------------------------------------------------------------------------------------------------------------------------------! ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -719,6 +702,9 @@ PPURE SUBROUTINE GetRefNewtonStartValue(X_in,Xi,ElemID) epsOne = 1.0 + RefMappingEps RefMappingGuessLoc = RefMappingGuess +! AXISYMMETRIC HDG +IF (Symmetry%Axisymmetric) RefMappingGuessLoc = 4 + SELECT CASE(RefMappingGuessLoc) CASE(1) diff --git a/src/io_hdf5/hdf5_input_field.f90 b/src/io_hdf5/hdf5_input_field.f90 index 2dc5df6aa..322882ade 100644 --- a/src/io_hdf5/hdf5_input_field.f90 +++ b/src/io_hdf5/hdf5_input_field.f90 @@ -75,7 +75,7 @@ SUBROUTINE ReadExternalFieldFromHDF5( DataSet, ExternalField, DeltaExternalField INTEGER(HID_T) :: file_id_loc ! File identifier INTEGER(HID_T) :: dset_id_loc ! Dataset identifier INTEGER(HID_T) :: filespace ! filespace identifier -LOGICAL :: DatasetFound,AttribtueFound,NaNDetected +LOGICAL :: DatasetFound,AttributeFound,NaNDetected REAL :: delta,deltaOld INTEGER :: iDirMax !=================================================================================================================================== @@ -116,8 +116,8 @@ SUBROUTINE ReadExternalFieldFromHDF5( DataSet, ExternalField, DeltaExternalField ExternalFieldRadInd = -1 ExternalFieldAxisSym = .FALSE. AttributeName = 'r' -CALL H5AEXISTS_F(file_id_loc, TRIM(AttributeName), AttribtueFound, iError) -IF(AttribtueFound) CALL ReadAttribute(file_id_loc,AttributeName,1,IntScalar=ExternalFieldRadInd) +CALL H5AEXISTS_F(file_id_loc, TRIM(AttributeName), AttributeFound, iError) +IF(AttributeFound) CALL ReadAttribute(file_id_loc,AttributeName,1,IntScalar=ExternalFieldRadInd) IF(ExternalFieldRadInd.GT.0) ExternalFieldAxisSym=.TRUE. ! Set spatial dimension @@ -129,7 +129,7 @@ SUBROUTINE ReadExternalFieldFromHDF5( DataSet, ExternalField, DeltaExternalField ! Check if axial symmetric and 2D IF(ExternalFieldDim.EQ.2)THEN - IF(.NOT.ExternalFieldAxisSym) CALL abort(__STAMP__,'Only 2D axis symmetric external field imeplemented.') + IF(.NOT.ExternalFieldAxisSym) CALL abort(__STAMP__,'Only 2D axis symmetric external field implemented.') ELSE ! 3D and symmetric is not possible ExternalFieldAxisSym = .FALSE. @@ -140,11 +140,11 @@ SUBROUTINE ReadExternalFieldFromHDF5( DataSet, ExternalField, DeltaExternalField ExternalFieldAxisDir=-1 ! Check z-dir AttributeName = 'z' - CALL H5AEXISTS_F(file_id_loc, TRIM(AttributeName), AttribtueFound, iError) - IF(AttribtueFound) CALL ReadAttribute(file_id_loc,AttributeName,1,IntScalar=ExternalFieldAxisDir) + CALL H5AEXISTS_F(file_id_loc, TRIM(AttributeName), AttributeFound, iError) + IF(AttributeFound) CALL ReadAttribute(file_id_loc,AttributeName,1,IntScalar=ExternalFieldAxisDir) IF(ExternalFieldAxisDir.GT.0) ExternalFieldAxisDir=3 ! Check if not axial symmetric with z-direction - IF(ExternalFieldAxisDir.NE.3) CALL abort(__STAMP__,'Only z-axis symmetric external field imeplemented currently.') + IF(ExternalFieldAxisDir.NE.3) CALL abort(__STAMP__,'Only z-axis symmetric external field implemented currently.') END IF ! ExternalFieldAxisSym ! Calculate the deltas and make sure that they are equidistant diff --git a/src/mesh/mesh.f90 b/src/mesh/mesh.f90 index 173534b68..122a23ca7 100644 --- a/src/mesh/mesh.f90 +++ b/src/mesh/mesh.f90 @@ -105,6 +105,7 @@ SUBROUTINE InitMesh(meshMode,MeshFile_IN) USE MOD_Prepare_Mesh ,ONLY: setLocalSideIDs,fillMeshInfo USE MOD_ReadInTools ,ONLY: PrintOption USE MOD_ReadInTools ,ONLY: GETLOGICAL,GETSTR,GETREAL,GETINT,GETREALARRAY +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_Prepare_Mesh ,ONLY: exchangeFlip #endif @@ -391,42 +392,47 @@ SUBROUTINE InitMesh(meshMode,MeshFile_IN) #endif /*ROS or IMPA*/ #endif /*maxwell*/ - ! assign all metrics Metrics_fTilde,Metrics_gTilde,Metrics_hTilde - ! assign 1/detJ (sJ) - ! assign normal and tangential vectors and surfElems on faces - crossProductMetrics=GETLOGICAL('crossProductMetrics','.FALSE.') - LBWRITE(UNIT_stdOut,'(A)') "NOW CALLING calcMetrics..." - CALL InitMeshBasis(NGeo,PP_N,xGP) - - ! get XCL_NGeo - ALLOCATE(XCL_NGeo(1:3,0:NGeo,0:NGeo,0:NGeo,1:nElems)) - XCL_NGeo = 0. -#ifdef PARTICLES - ALLOCATE(dXCL_NGeo(1:3,1:3,0:NGeo,0:NGeo,0:NGeo,1:nElems)) - dXCL_NGeo = 0. - CALL CalcMetrics(XCL_NGeo_Out=XCL_NGeo,dXCL_NGeo_Out=dXCL_NGeo) -#else - CALL CalcMetrics(XCL_NGeo_Out=XCL_NGeo) -#endif /*PARTICLES*/ -#if USE_LOADBALANCE - END IF -#endif /*USE_LOADBALANCE*/ - - ! surface data - ALLOCATE(Face_xGP (3,0:PP_N,0:PP_N,1:nSides)) - ALLOCATE(NormVec (3,0:PP_N,0:PP_N,1:nSides)) - ALLOCATE(TangVec1 (3,0:PP_N,0:PP_N,1:nSides)) - ALLOCATE(TangVec2 (3,0:PP_N,0:PP_N,1:nSides)) - ALLOCATE(SurfElem ( 0:PP_N,0:PP_N,1:nSides)) - ALLOCATE( Ja_Face(3,3,0:PP_N,0:PP_N,1:nSides)) ! temp - Face_xGP = 0. - NormVec = 0. - TangVec1 = 0. - TangVec2 = 0. - SurfElem = 0. - - DO iElem=1,nElems - CALL CalcSurfMetrics(JaCL_N(:,:,:,:,:,iElem),iElem) + ! assign all metrics Metrics_fTilde,Metrics_gTilde,Metrics_hTilde + ! assign 1/detJ (sJ) + ! assign normal and tangential vectors and surfElems on faces + + crossProductMetrics=GETLOGICAL('crossProductMetrics','.FALSE.') + IF(Symmetry%axisymmetric.AND..NOT.crossProductMetrics) THEN + crossProductMetrics = .TRUE. + LBWRITE(UNIT_stdOut,'(A)') 'WARNING: setting ""crossProductMetrics" to true for axisymmetric simulations!' + END IF + LBWRITE(UNIT_stdOut,'(A)') "NOW CALLING calcMetrics..." + CALL InitMeshBasis(NGeo,PP_N,xGP) + + ! get XCL_NGeo + ALLOCATE(XCL_NGeo(1:3,0:NGeo,0:NGeo,0:NGeo,1:nElems)) + XCL_NGeo = 0. +#ifdef PARTICLES + ALLOCATE(dXCL_NGeo(1:3,1:3,0:NGeo,0:NGeo,0:NGeo,1:nElems)) + dXCL_NGeo = 0. + CALL CalcMetrics(XCL_NGeo_Out=XCL_NGeo,dXCL_NGeo_Out=dXCL_NGeo) +#else + CALL CalcMetrics(XCL_NGeo_Out=XCL_NGeo) +#endif /*PARTICLES*/ +#if USE_LOADBALANCE + END IF +#endif /*USE_LOADBALANCE*/ + + ! surface data + ALLOCATE(Face_xGP (3,0:PP_N,0:PP_N,1:nSides)) + ALLOCATE(NormVec (3,0:PP_N,0:PP_N,1:nSides)) + ALLOCATE(TangVec1 (3,0:PP_N,0:PP_N,1:nSides)) + ALLOCATE(TangVec2 (3,0:PP_N,0:PP_N,1:nSides)) + ALLOCATE(SurfElem ( 0:PP_N,0:PP_N,1:nSides)) + ALLOCATE( Ja_Face(3,3,0:PP_N,0:PP_N,1:nSides)) ! temp + Face_xGP = 0. + NormVec = 0. + TangVec1 = 0. + TangVec2 = 0. + SurfElem = 0. + + DO iElem=1,nElems + CALL CalcSurfMetrics(JaCL_N(:,:,:,:,:,iElem),iElem) END DO ! Compute element bary and element radius for processor-local elements (without halo region) diff --git a/src/mesh/mesh_readin.f90 b/src/mesh/mesh_readin.f90 index 0654b7e7d..9bc6b4459 100644 --- a/src/mesh/mesh_readin.f90 +++ b/src/mesh/mesh_readin.f90 @@ -239,6 +239,10 @@ SUBROUTINE ReadMesh(FileString,ReadNodes) #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance,UseH5IOLoadBalance,offsetElemMPIOld #endif /*USE_LOADBALANCE*/ +#if USE_HDG +! Axisymmetric HDG +USE MOD_Symmetry_Vars ,ONLY: Symmetry +#endif /*USE_HDG*/ ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -495,6 +499,22 @@ SUBROUTINE ReadMesh(FileString,ReadNodes) ! ALLOCATE MORTAR ElemID=SideInfo(SIDE_NBELEMID,iSide) !IF nbElemID <0, this marks a mortar master side. ! The number (-1,-2,-3) is the Type of mortar + +#if USE_HDG + ! AXISYMMETRIC HDG + IF(Symmetry%Axisymmetric) THEN + ! In 2D check that there is only one layer of elements in z-direction + IF ((iLocSide.EQ.1).OR.(iLocSide.EQ.6)) THEN + BCindex = SideInfo(SIDE_BCID,iSide) + IF ((iElem.NE.ElemID).AND.(BCindex.LE.0)) THEN + CALL Abort(__STAMP__, & + "Mesh not oriented in z-direction or more than one layer of elements in z-direction! " // & + "Please set 'orientZ = T' or change number of element in z-direction in HOPR parameter file.") + END IF + END IF + END IF +#endif /*USE_HDG*/ + IF(ElemID.LT.0)THEN ! mortar Sides attached! aSide%MortarType=ABS(ElemID) SELECT CASE(aSide%MortarType) @@ -1345,7 +1365,7 @@ SUBROUTINE FinalizeMeshReadin(meshMode) #if USE_MPI USE MOD_MPI_Shared USE MOD_MPI_Shared_Vars -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #endif #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance diff --git a/src/mesh/metrics.f90 b/src/mesh/metrics.f90 index fbf74a6b1..59167cf37 100644 --- a/src/mesh/metrics.f90 +++ b/src/mesh/metrics.f90 @@ -143,6 +143,8 @@ SUBROUTINE CalcMetrics(XCL_NGeo_Out,dXCL_NGeo_out) #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_LOADBALANCE*/ +USE MOD_Symmetry_Vars ,ONLY: Symmetry +USE MOD_Globals_Vars ,ONLY: PI !----------------------------------------------------------------------------------------------------------------------------------- ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -252,6 +254,8 @@ SUBROUTINE CalcMetrics(XCL_NGeo_Out,dXCL_NGeo_out) dXCL_Ngeo(2,:,i,j,k)=dXCL_Ngeo(2,:,i,j,k) + DCL_Ngeo(j,ll)*XCL_Ngeo(:,i,ll,k) dXCL_Ngeo(3,:,i,j,k)=dXCL_Ngeo(3,:,i,j,k) + DCL_Ngeo(k,ll)*XCL_Ngeo(:,i,j,ll) END DO !l=0,N + ! AXISYMMETRIC HDG + IF(Symmetry%Axisymmetric) dXCL_Ngeo(3,3,i,j,k)=PI*XCL_Ngeo(2,i,j,k) END DO; END DO; END DO !i,j,k=0,Ngeo ! 1.c)Jacobians! grad(X_1) (grad(X_2) x grad(X_3)) @@ -321,6 +325,8 @@ SUBROUTINE CalcMetrics(XCL_NGeo_Out,dXCL_NGeo_out) dXCL(2,:)=dXCL(2,:) + DCL_N(j,ll)*XCL_N(:,i,ll,k,iElem) dXCL(3,:)=dXCL(3,:) + DCL_N(k,ll)*XCL_N(:,i,j,ll,iElem) END DO !l=0,N + ! AXISYMMETRIC HDG + IF(Symmetry%Axisymmetric) dXCL(:,3)=PI * XCL_N(2,i,j,k,iElem) END ASSOCIATE END DO; END DO; END DO !i,j,k=0,N END IF !N>=Ngeo @@ -345,6 +351,7 @@ SUBROUTINE CalcMetrics(XCL_NGeo_Out,dXCL_NGeo_out) END ASSOCIATE END DO; END DO; END DO !i,j,k=0,N ELSE ! curl metrics + ! AXISYMMETRIC HDG: Does not work without some work :( ! invariant curl form, as cross product: R^dd = 1/2( XCL_N(:) x (d/dxi_dd XCL_N(:))) ! !R_CL_N(dd,nn)=1/2*( XCL_N(nn+2)* d/dxi_dd XCL_N(nn+1) - XCL_N(nn+1)* d/dxi_dd XCL_N(nn+2)) @@ -546,7 +553,8 @@ END SUBROUTINE CalcSurfMetrics !================================================================================================================================== SUBROUTINE SurfMetricsFromJa(Nloc,NormalDir,TangDir,NormalSign,Ja_Face,NormVec,TangVec1,TangVec2,SurfElem) ! MODULES -USE MOD_Globals, ONLY: CROSS +USE MOD_Globals, ONLY: CROSS +USE MOD_Symmetry_Vars ,ONLY: Symmetry !---------------------------------------------------------------------------------------------------------------------------------- IMPLICIT NONE !---------------------------------------------------------------------------------------------------------------------------------- @@ -566,11 +574,19 @@ SUBROUTINE SurfMetricsFromJa(Nloc,NormalDir,TangDir,NormalSign,Ja_Face,NormVec,T !================================================================================================================================== DO q=0,Nloc; DO p=0,Nloc SurfElem( p,q) = SQRT(SUM(Ja_Face(NormalDir,:,p,q)**2)) - NormVec( :,p,q) = NormalSign*Ja_Face(NormalDir,:,p,q)/SurfElem(p,q) - TangVec1(:,p,q) = Ja_Face(TangDir,:,p,q) - SUM(Ja_Face(TangDir,:,p,q)*NormVec(:,p,q)) & - *NormVec(:,p,q) - TangVec1(:,p,q) = TangVec1(:,p,q)/SQRT(SUM(TangVec1(:,p,q)**2)) - TangVec2(:,p,q) = CROSS(NormVec(:,p,q),TangVec1(:,p,q)) + ! AXISYMMETRIC HDG: Fixed SurfElem=0 at r=0... + IF(Symmetry%Axisymmetric .AND. SurfElem(p,q).EQ.0.) THEN + NormVec(:,p,q) = (/0., -1., 0./) + TangVec1(:,p,q) = (/1., 0., 0./) + TangVec2(:,p,q) = (/0., 0., 1./) + ELSE + NormVec( :,p,q) = NormalSign*Ja_Face(NormalDir,:,p,q)/SurfElem(p,q) + TangVec1(:,p,q) = Ja_Face(TangDir,:,p,q) - SUM(Ja_Face(TangDir,:,p,q)*NormVec(:,p,q)) & + *NormVec(:,p,q) + TangVec1(:,p,q) = TangVec1(:,p,q)/SQRT(SUM(TangVec1(:,p,q)**2)) + TangVec2(:,p,q) = CROSS(NormVec(:,p,q),TangVec1(:,p,q)) + END IF + END DO; END DO ! p,q END SUBROUTINE SurfMetricsFromJa diff --git a/src/particles/bgk/bgk_adaptation.f90 b/src/particles/bgk/bgk_adaptation.f90 index f9d39cb2c..cc58ff765 100644 --- a/src/particles/bgk/bgk_adaptation.f90 +++ b/src/particles/bgk/bgk_adaptation.f90 @@ -396,7 +396,7 @@ SUBROUTINE BGK_AllocateAveragingNode(Averaging) !=================================================================================================================================== ! MODULES USE MOD_BGK_Vars, ONLY: tNodeAverage,BGKCollModel -USE MOD_Particle_Vars, ONLY: Symmetry +USE MOD_Symmetry_Vars, ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- diff --git a/src/particles/bgk/bgk_init.f90 b/src/particles/bgk/bgk_init.f90 index 9b6ff5b1c..54050a849 100644 --- a/src/particles/bgk/bgk_init.f90 +++ b/src/particles/bgk/bgk_init.f90 @@ -92,7 +92,8 @@ SUBROUTINE InitBGK() USE MOD_BGK_Vars USE MOD_Preproc USE MOD_Mesh_Vars ,ONLY: nElems, NGeo -USE MOD_Particle_Vars ,ONLY: nSpecies, Species, DoVirtualCellMerge, Symmetry +USE MOD_Particle_Vars ,ONLY: nSpecies, Species, DoVirtualCellMerge +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_DSMC_Vars ,ONLY: DSMC, RadialWeighting, CollInf USE MOD_DSMC_ParticlePairing ,ONLY: DSMC_init_octree USE MOD_Globals_Vars ,ONLY: Pi, BoltzmannConst @@ -287,7 +288,7 @@ RECURSIVE SUBROUTINE DeleteNodeAverage(NodeAverage) !----------------------------------------------------------------------------------------------------------------------------------! USE MOD_Globals USE MOD_BGK_Vars -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT VARIABLES diff --git a/src/particles/bgk/bgk_main.f90 b/src/particles/bgk/bgk_main.f90 index d856dd6dc..5a44cd380 100644 --- a/src/particles/bgk/bgk_main.f90 +++ b/src/particles/bgk/bgk_main.f90 @@ -43,7 +43,7 @@ SUBROUTINE BGK_DSMC_main(stage_opt) ! MODULES USE MOD_Globals USE MOD_BGK_Adaptation ,ONLY: BGK_octree_adapt, BGK_quadtree_adapt -USE MOD_Particle_Vars ,ONLY: PEM, Species, WriteMacroVolumeValues, Symmetry, usevMPF +USE MOD_Particle_Vars ,ONLY: PEM, Species, WriteMacroVolumeValues, usevMPF USE MOD_BGK_Vars ,ONLY: DoBGKCellAdaptation,BGKDSMCSwitchDens USE MOD_BGK_Vars ,ONLY: BGKMovingAverage,ElemNodeAveraging USE MOD_BGK_Vars ,ONLY: BGK_MeanRelaxFactor,BGK_MeanRelaxFactorCounter,BGK_MaxRelaxFactor,BGK_QualityFacSamp @@ -57,6 +57,7 @@ SUBROUTINE BGK_DSMC_main(stage_opt) USE MOD_TimeDisc_Vars ,ONLY: TEnd, Time USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared USE MOD_Mesh_Tools ,ONLY: GetCNElemID +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_Particle_Mesh_Vars ,ONLY: IsExchangeElem #endif @@ -169,7 +170,7 @@ SUBROUTINE BGK_main(stage_opt) USE MOD_TimeDisc_Vars ,ONLY: TEnd, Time USE MOD_Mesh_Vars ,ONLY: nElems, offsetElem USE MOD_BGK_Adaptation ,ONLY: BGK_octree_adapt, BGK_quadtree_adapt -USE MOD_Particle_Vars ,ONLY: PEM, WriteMacroVolumeValues, WriteMacroSurfaceValues, Symmetry, DoVirtualCellMerge, VirtMergedCells +USE MOD_Particle_Vars ,ONLY: PEM, WriteMacroVolumeValues, WriteMacroSurfaceValues, DoVirtualCellMerge, VirtMergedCells USE MOD_BGK_Vars ,ONLY: DoBGKCellAdaptation, BGKMovingAverage, ElemNodeAveraging USE MOD_BGK_Vars ,ONLY: BGK_MeanRelaxFactor,BGK_MeanRelaxFactorCounter,BGK_MaxRelaxFactor,BGK_QualityFacSamp USE MOD_BGK_Vars ,ONLY: BGK_MaxRotRelaxFactor, BGK_PrandtlNumber, BGK_ExpectedPrandtlNumber @@ -179,6 +180,7 @@ SUBROUTINE BGK_main(stage_opt) USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared USE MOD_DSMC_Vars ,ONLY: DSMC USE MOD_Mesh_Tools ,ONLY: GetCNElemID +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_Particle_Mesh_Vars ,ONLY: IsExchangeElem #endif diff --git a/src/particles/boundary/particle_boundary_condition.f90 b/src/particles/boundary/particle_boundary_condition.f90 index 88f3f0bd2..dafd81764 100644 --- a/src/particles/boundary/particle_boundary_condition.f90 +++ b/src/particles/boundary/particle_boundary_condition.f90 @@ -481,8 +481,6 @@ SUBROUTINE RotPeriodicInterPlaneBoundary(PartID,SideID,ElemID,IsInterPlanePart) USE MOD_DSMC_Vars ,ONLY: CollisMode, useDSMC, PartStateIntEn USE MOD_Particle_Vars ,ONLY: PDM,InterPlanePartNumber, InterPlanePartIndx USE MOD_Particle_Vars ,ONLY: UseRotRefFrame, RotRefFrameOmega, PartVeloRotRef, LastPartVeloRotRef -USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck -USE MOD_part_RHS ,ONLY: CalcPartRHSRotRefFrame USE MOD_part_tools ,ONLY: RotateVectorAroundAxis USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, PartTimeStep #ifdef CODE_ANALYZE diff --git a/src/particles/boundary/particle_boundary_sampling.f90 b/src/particles/boundary/particle_boundary_sampling.f90 index 40ca4784b..f8354bb5f 100644 --- a/src/particles/boundary/particle_boundary_sampling.f90 +++ b/src/particles/boundary/particle_boundary_sampling.f90 @@ -74,7 +74,6 @@ SUBROUTINE InitParticleBoundarySampling() !----------------------------------------------------------------------------------------------------------------------------------! USE MOD_Globals USE MOD_Basis ,ONLY: LegendreGaussNodesAndWeights -USE MOD_DSMC_Symmetry ,ONLY: DSMC_2D_CalcSymmetryArea, DSMC_1D_CalcSymmetryArea USE MOD_Mesh_Vars ,ONLY: NGeo,nBCs,BoundaryName USE MOD_Particle_Boundary_Vars ,ONLY: SurfTotalSideOnNode USE MOD_Particle_Boundary_Vars ,ONLY: nSurfSample,dXiEQ_SurfSample,PartBound,XiEQ_SurfSample @@ -99,8 +98,9 @@ SUBROUTINE InitParticleBoundarySampling() USE MOD_Particle_Surfaces_Vars ,ONLY: BezierControlPoints3D USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod USE MOD_Particle_Vars ,ONLY: nSpecies,UseVarTimeStep,VarTimeStep -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_ReadInTools ,ONLY: GETINT,GETLOGICAL,GETINTARRAY +USE MOD_Particle_Mesh_Tools ,ONLY: DSMC_2D_CalcSymmetryArea, DSMC_1D_CalcSymmetryArea #if USE_MPI USE MOD_MPI_Shared USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED @@ -470,7 +470,8 @@ SUBROUTINE CalcSurfaceValues(during_dt_opt) USE MOD_Particle_Boundary_vars ,ONLY: SurfOutputSize, SWIVarTimeStep, SWIStickingCoefficient USE MOD_Particle_Boundary_Vars ,ONLY: MacroSurfaceVal, MacroSurfaceSpecVal USE MOD_Particle_Mesh_Vars ,ONLY: SideInfo_Shared -USE MOD_Particle_Vars ,ONLY: WriteMacroSurfaceValues,nSpecies,MacroValSampTime,UseVarTimeStep,Symmetry,VarTimeStep +USE MOD_Particle_Vars ,ONLY: WriteMacroSurfaceValues,nSpecies,MacroValSampTime,UseVarTimeStep,VarTimeStep +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Particle_Vars ,ONLY: Species USE MOD_Restart_Vars ,ONLY: RestartTime USE MOD_TimeDisc_Vars ,ONLY: TEnd diff --git a/src/particles/boundary/particle_boundary_tools.f90 b/src/particles/boundary/particle_boundary_tools.f90 index 39f61f53a..6a91876f8 100644 --- a/src/particles/boundary/particle_boundary_tools.f90 +++ b/src/particles/boundary/particle_boundary_tools.f90 @@ -323,7 +323,7 @@ SUBROUTINE GetRadialDistance2D(GlobalSideID,dir,origin,rmin,rmax) USE MOD_Particle_Surfaces ,ONLY: GetSideBoundingBox USE MOD_Particle_Mesh_Tools ,ONLY: GetSideBoundingBoxTria USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry !----------------------------------------------------------------------------------------------------------------------------------! ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE diff --git a/src/particles/dsmc/dsmc_ambipolardiff.f90 b/src/particles/dsmc/dsmc_ambipolardiff.f90 index e66075cec..d8ad50a23 100644 --- a/src/particles/dsmc/dsmc_ambipolardiff.f90 +++ b/src/particles/dsmc/dsmc_ambipolardiff.f90 @@ -156,9 +156,10 @@ SUBROUTINE AD_InsertParticles(iPartIndx_Node, nPart, iPartIndx_NodeTotalAmbi, To 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: PDM, PartSpecies, PartState, PEM, Species, PartMPF, usevMPF USE MOD_PARTICLE_Vars ,ONLY: UseVarTimeStep, PartTimeStep USE MOD_Particle_Tracking ,ONLY: ParticleInsideCheck +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Part_Tools ,ONLY: GetNextFreePosition ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE diff --git a/src/particles/dsmc/dsmc_analyze.f90 b/src/particles/dsmc/dsmc_analyze.f90 index 9f0537b8e..691659e80 100644 --- a/src/particles/dsmc/dsmc_analyze.f90 +++ b/src/particles/dsmc/dsmc_analyze.f90 @@ -524,7 +524,7 @@ SUBROUTINE DSMC_output_calc(nVar,nVar_quality,nVarloc,nVar_HeatPress,DSMC_MacroV USE MOD_DSMC_Vars ,ONLY: CollisMode, SpecDSMC, DSMC, useDSMC, RadialWeighting, BGGas USE MOD_FPFlow_Vars ,ONLY: FPInitDone, FP_QualityFacSamp USE MOD_Mesh_Vars ,ONLY: nElems -USE MOD_Particle_Vars ,ONLY: Species, nSpecies, WriteMacroVolumeValues, usevMPF, Symmetry +USE MOD_Particle_Vars ,ONLY: Species, nSpecies, WriteMacroVolumeValues, usevMPF USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, VarTimeStep USE MOD_Particle_Vars ,ONLY: DoVirtualCellMerge, VirtMergedCells, SamplePressTensHeatflux USE MOD_Particle_TimeStep ,ONLY: GetParticleTimeStep @@ -534,6 +534,7 @@ SUBROUTINE DSMC_output_calc(nVar,nVar_quality,nVarloc,nVar_HeatPress,DSMC_MacroV USE MOD_Mesh_Vars ,ONLY: offSetElem USE MOD_Mesh_Tools ,ONLY: GetCNElemID USE MOD_Particle_Analyze_Tools ,ONLY: CalcTelec,CalcTVibPoly +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- diff --git a/src/particles/dsmc/dsmc_bg_gas.f90 b/src/particles/dsmc/dsmc_bg_gas.f90 index 8b18739cd..a1d1f79ac 100644 --- a/src/particles/dsmc/dsmc_bg_gas.f90 +++ b/src/particles/dsmc/dsmc_bg_gas.f90 @@ -87,7 +87,7 @@ SUBROUTINE BGGas_Initialize() ! 1.) Check compatibility with other features and whether required parameters have been read-in IF(UseVarTimeStep) THEN IF(.NOT.VarTimeStep%UseSpeciesSpecific) CALL abort(__STAMP__, & - 'ERROR: 2D/Axisymmetric and variable timestep (except species-specific) are not implemented with a background gas yet!') + 'ERROR: Variable timestep (except species-specific) are not implemented with a background gas yet!') END IF DO iSpec = 1, nSpecies diff --git a/src/particles/dsmc/dsmc_chemical_reactions.f90 b/src/particles/dsmc/dsmc_chemical_reactions.f90 index 3f36e5584..579b85884 100644 --- a/src/particles/dsmc/dsmc_chemical_reactions.f90 +++ b/src/particles/dsmc/dsmc_chemical_reactions.f90 @@ -363,7 +363,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac) USE MOD_part_operations ,ONLY: RemoveParticle #ifdef CODE_ANALYZE USE MOD_Globals ,ONLY: unit_stdout,myrank -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #endif /* CODE_ANALYZE */ USE MOD_Particle_Analyze_Vars ,ONLY: CalcPartBalance,nPartIn,nPartOut,PartEkinIn,PartEkinOut USE MOD_Particle_Analyze_Tools ,ONLY: CalcEkinPart diff --git a/src/particles/dsmc/dsmc_collis_mode.f90 b/src/particles/dsmc/dsmc_collis_mode.f90 index 9a9f0787a..84c9c7f8c 100644 --- a/src/particles/dsmc/dsmc_collis_mode.f90 +++ b/src/particles/dsmc/dsmc_collis_mode.f90 @@ -47,7 +47,7 @@ SUBROUTINE DSMC_Elastic_Col(iPair) #ifdef CODE_ANALYZE USE MOD_Globals ,ONLY: Abort USE MOD_Globals ,ONLY: unit_stdout,myrank -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #endif /* CODE_ANALYZE */ USE MOD_DSMC_Vars ,ONLY: DSMC ! IMPLICIT VARIABLE HANDLING @@ -922,8 +922,9 @@ SUBROUTINE DSMC_perform_collision(iPair, iElem, NodeVolume, NodePartNum) USE MOD_Globals ,ONLY: Abort, CROSS USE MOD_DSMC_Vars ,ONLY: CollisMode, Coll_pData, SelectionProc USE MOD_DSMC_Vars ,ONLY: DSMC -USE MOD_Particle_Vars ,ONLY: PartState, Symmetry +USE MOD_Particle_Vars ,ONLY: PartState USE MOD_Particle_Vars ,ONLY: UseRotRefFrame, PDM, PartVeloRotRef, RotRefFrameOmega +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_DSMC_Vars ,ONLY: RadialWeighting USE MOD_Particle_Vars ,ONLY: usevMPF, Species, PartSpecies USE MOD_Particle_Analyze_Vars ,ONLY: CalcCollRates diff --git a/src/particles/dsmc/dsmc_init.f90 b/src/particles/dsmc/dsmc_init.f90 index c27b185d0..73bb7dd51 100644 --- a/src/particles/dsmc/dsmc_init.f90 +++ b/src/particles/dsmc/dsmc_init.f90 @@ -279,8 +279,9 @@ SUBROUTINE InitDSMC() USE MOD_DSMC_Vars USE MOD_Mesh_Vars ,ONLY: nElems, NGEo USE MOD_Globals_Vars ,ONLY: Pi, BoltzmannConst, ElementaryCharge -USE MOD_Particle_Vars ,ONLY: nSpecies, Species, PDM, Symmetry, UseVarTimeStep, usevMPF -USE MOD_Particle_Vars ,ONLY: DoFieldIonization, SpeciesDatabase, SampleElecExcitation +USE MOD_Particle_Vars ,ONLY: nSpecies, Species, PDM, UseVarTimeStep, usevMPF +USE MOD_Symmetry_Vars ,ONLY: Symmetry +USE MOD_Particle_Vars ,ONLY: DoFieldIonization, SpeciesDatabase, SampleElecExcitation USE MOD_DSMC_ParticlePairing ,ONLY: DSMC_init_octree USE MOD_DSMC_ChemInit ,ONLY: DSMC_chemical_init USE MOD_DSMC_PolyAtomicModel ,ONLY: InitPolyAtomicMolecs @@ -1720,7 +1721,6 @@ SUBROUTINE FinalizeDSMC() SDEALLOCATE(RadialWeighting%ClonePartNum) SDEALLOCATE(ClonedParticles) -SDEALLOCATE(SymmetrySide) SDEALLOCATE(AmbiPolarSFMapping) END SUBROUTINE FinalizeDSMC @@ -1760,7 +1760,7 @@ RECURSIVE SUBROUTINE DeleteNodeVolume(Node) !----------------------------------------------------------------------------------------------------------------------------------! USE MOD_Globals USE MOD_DSMC_Vars -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT VARIABLES diff --git a/src/particles/dsmc/dsmc_main.f90 b/src/particles/dsmc/dsmc_main.f90 index 057c3db91..26426026f 100644 --- a/src/particles/dsmc/dsmc_main.f90 +++ b/src/particles/dsmc/dsmc_main.f90 @@ -46,9 +46,10 @@ SUBROUTINE DSMC_main(DoElement) USE MOD_DSMC_Vars ,ONLY: DSMC, CollInf, DSMCSumOfFormedParticles, BGGas, CollisMode, ElecRelaxPart USE MOD_DSMC_Analyze ,ONLY: SummarizeQualityFactors, DSMCMacroSampling USE MOD_DSMC_Relaxation ,ONLY: FinalizeCalcVibRelaxProb, InitCalcVibRelaxProb -USE MOD_Particle_Vars ,ONLY: PEM, PDM, WriteMacroVolumeValues, Symmetry +USE MOD_Particle_Vars ,ONLY: PEM, PDM, WriteMacroVolumeValues USE MOD_DSMC_ParticlePairing ,ONLY: DSMC_pairing_standard, DSMC_pairing_octree, DSMC_pairing_quadtree, DSMC_pairing_dotree USE MOD_Particle_Vars ,ONLY: WriteMacroSurfaceValues +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_LOADBALANCE USE MOD_LoadBalance_Timers ,ONLY: LBStartTime, LBElemSplitTime #endif /*USE_LOADBALANCE*/ diff --git a/src/particles/dsmc/dsmc_particle_pairing.f90 b/src/particles/dsmc/dsmc_particle_pairing.f90 index 0a9392e38..ad8308198 100644 --- a/src/particles/dsmc/dsmc_particle_pairing.f90 +++ b/src/particles/dsmc/dsmc_particle_pairing.f90 @@ -442,7 +442,7 @@ SUBROUTINE PerformPairingAndCollision(iPartIndx_Node, PartNum, iElem, NodeVolume USE MOD_DSMC_Vars ,ONLY: Coll_pData,CollInf,CollisMode,PartStateIntEn,ChemReac,DSMC,RadialWeighting USE MOD_DSMC_Vars ,ONLY: SelectionProc, useRelaxProbCorrFactor, iPartIndx_NodeNewElecRelax, newElecRelaxParts USE MOD_DSMC_Vars ,ONLY: iPartIndx_NodeElecRelaxChem,nElecRelaxChemParts -USE MOD_Particle_Vars ,ONLY: PartSpecies, nSpecies, PartState, UseVarTimeStep, Symmetry, usevMPF +USE MOD_Particle_Vars ,ONLY: PartSpecies, nSpecies, PartState, UseVarTimeStep, usevMPF USE MOD_DSMC_Analyze ,ONLY: CalcGammaVib, CalcInstantTransTemp, CalcMeanFreePath, CalcInstantElecTempXi USE MOD_part_tools ,ONLY: GetParticleWeight USE MOD_DSMC_Relaxation ,ONLY: CalcMeanVibQuaDiatomic,SumVibRelaxProb @@ -450,6 +450,7 @@ SUBROUTINE PerformPairingAndCollision(iPartIndx_Node, PartNum, iElem, NodeVolume USE MOD_DSMC_AmbipolarDiffusion ,ONLY: AD_InsertParticles, AD_DeleteParticles USE MOD_DSMC_ElectronicModel ,ONLY: LT_ElectronicEnergyExchange, LT_ElectronicExc_ConstructPartList USE MOD_DSMC_ElectronicModel ,ONLY: CalcProbCorrFactorElec, LT_ElectronicEnergyExchangeChem +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -1242,7 +1243,7 @@ SUBROUTINE DSMC_init_octree() USE MOD_ReadInTools USE MOD_DSMC_Vars ,ONLY: DSMC, ElemNodeVol USE MOD_Mesh_Vars ,ONLY: nElems -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -1298,7 +1299,8 @@ SUBROUTINE DSMC_CalcSubNodeVolumes2D(iElem, NodeDepth, Node) ! Pairing subroutine for octree and nearest neighbour, decides whether to create a new octree node or start nearest neighbour search !=================================================================================================================================== ! MODULES -USE MOD_DSMC_Vars, ONLY : OctreeVdm, tNodeVolume, SymmetrySide +USE MOD_DSMC_Vars, ONLY : OctreeVdm, tNodeVolume +USE MOD_Particle_Mesh_Vars ,ONLY: SymmetrySide USE MOD_Mesh_Vars, ONLY : SurfElem, Face_xGP USE MOD_Preproc USE MOD_ChangeBasis, ONLY : ChangeBasis2D @@ -1515,7 +1517,7 @@ RECURSIVE SUBROUTINE AddNodeVolumes2D(NodeDepth, Node, DetJac, VdmLocal, iElem, ! MODULES USE MOD_Globals_Vars ,ONLY: Pi USE MOD_DSMC_Vars ,ONLY: tOctreeVdm, tNodeVolume -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -1902,8 +1904,7 @@ SUBROUTINE GeoCoordToMap2D(x_in,xi_Out,iElem) !> xi is defined in the 1DrefElem xi=[-1,1] !=================================================================================================================================== ! MODULES -USE MOD_DSMC_Vars ,ONLY: SymmetrySide -USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared +USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared,SymmetrySide USE MOD_Mesh_Vars ,ONLY: offsetElem USE MOD_Mesh_Tools ,ONLY: GetCNElemID ! IMPLICIT VARIABLE HANDLING @@ -1988,8 +1989,7 @@ FUNCTION MapToGeo2D(xi,iElem) !> !=================================================================================================================================== ! MODULES -USE MOD_DSMC_Vars ,ONLY: SymmetrySide -USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared +USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared, SymmetrySide USE MOD_Mesh_Vars ,ONLY: offsetElem USE MOD_Mesh_Tools ,ONLY: GetCNElemID ! IMPLICIT VARIABLE HANDLING @@ -2107,7 +2107,7 @@ RECURSIVE SUBROUTINE AddNodeVolumes1D(NodeDepth, Node, iElem, SubNodesIn) USE MOD_Globals USE MOD_Globals_Vars ,ONLY: Pi USE MOD_DSMC_Vars ,ONLY: tOctreeVdm, tNodeVolume -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Particle_Mesh_Vars ,ONLY: GEO ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE diff --git a/src/particles/dsmc/dsmc_symmetry.f90 b/src/particles/dsmc/dsmc_symmetry.f90 index 0a3279531..7b179399a 100644 --- a/src/particles/dsmc/dsmc_symmetry.f90 +++ b/src/particles/dsmc/dsmc_symmetry.f90 @@ -26,9 +26,8 @@ MODULE MOD_DSMC_Symmetry !----------------------------------------------------------------------------------------------------------------------------------- ! Private Part --------------------------------------------------------------------------------------------------------------------- ! Public Part ---------------------------------------------------------------------------------------------------------------------- -PUBLIC :: DSMC_1D_InitVolumes, DSMC_2D_InitVolumes, DSMC_2D_InitRadialWeighting, DSMC_2D_RadialWeighting, DSMC_2D_SetInClones -PUBLIC :: DSMC_2D_CalcSymmetryArea, DSMC_1D_CalcSymmetryArea, DSMC_2D_CalcSymmetryAreaSubSides -PUBLIC :: DefineParametersParticleSymmetry, DSMC_2D_TreatIdenticalParticles +PUBLIC :: DefineParametersParticleSymmetry +PUBLIC :: DSMC_2D_InitRadialWeighting, DSMC_2D_RadialWeighting, DSMC_2D_SetInClones, DSMC_2D_TreatIdenticalParticles !=================================================================================================================================== CONTAINS @@ -41,11 +40,7 @@ SUBROUTINE DefineParametersParticleSymmetry() USE MOD_ReadInTools ,ONLY: prms,addStrListEntry IMPLICIT NONE -CALL prms%SetSection("Particle Symmetry") -CALL prms%CreateIntOption( 'Particles-Symmetry-Order', & - 'Order of the Simulation 1, 2 or 3 D', '3') -CALL prms%CreateLogicalOption('Particles-Symmetry2DAxisymmetric', 'Activating an axisymmetric simulation with the same mesh '//& - 'requirements as for the 2D case (y is then the radial direction)', '.FALSE.') +CALL prms%SetSection("Particle Radial Weighting") CALL prms%CreateLogicalOption('Particles-RadialWeighting', 'Activates a radial weighting in y for the axisymmetric '//& 'simulation based on the particle position.', '.FALSE.') CALL prms%CreateRealOption( 'Particles-RadialWeighting-PartScaleFactor', 'Axisymmetric radial weighting factor, defining '//& @@ -66,361 +61,6 @@ SUBROUTINE DefineParametersParticleSymmetry() END SUBROUTINE DefineParametersParticleSymmetry -SUBROUTINE DSMC_2D_InitVolumes() -!=================================================================================================================================== -!> Routine determines the symmetry sides and calculates the 2D (area faces in symmetry plane) and axisymmetric volumes (cells are -!> revolved around the symmetry axis). The symmetry side (SymmetrySide array) will be used later on to determine in which two -!> directions the quadtree shall refine the mesh, skipping the z-dimension to avoid an unnecessary refinement. Additionally, -!> symmetry sides will be skipped during tracking (SideIsSymSide_Shared array). -!=================================================================================================================================== -! MODULES -USE MOD_Globals -USE MOD_Globals_Vars ,ONLY: Pi -USE MOD_PreProc -USE MOD_Mesh_Vars ,ONLY: nElems,offsetElem,nBCSides,SideToElem -USE MOD_Particle_Vars ,ONLY: Symmetry -USE MOD_Particle_Boundary_Vars ,ONLY: PartBound -USE MOD_Particle_Mesh_Vars ,ONLY: GEO,LocalVolume,MeshVolume -USE MOD_DSMC_Vars ,ONLY: SymmetrySide -USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared,ElemCharLength_Shared -USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared, SideInfo_Shared, SideIsSymSide -USE MOD_Mesh_Tools ,ONLY: GetCNElemID -USE MOD_Particle_Mesh_Tools ,ONLY: GetGlobalNonUniqueSideID -USE MOD_Particle_Surfaces ,ONLY: CalcNormAndTangTriangle -#if USE_MPI -USE MOD_MPI_Shared -USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED -USE MOD_Particle_Mesh_Vars ,ONLY: nNonUniqueGlobalSides, offsetComputeNodeElem, ElemInfo_Shared -USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared_Win,ElemCharLength_Shared_Win,SideIsSymSide_Shared,SideIsSymSide_Shared_Win -USE MOD_MPI_Shared_Vars ,ONLY: myComputeNodeRank,nComputeNodeProcessors,MPI_COMM_LEADERS_SHARED -#else -USE MOD_Particle_Mesh_Vars ,ONLY: nComputeNodeSides -#endif /*USE_MPI*/ -! IMPLICIT VARIABLE HANDLING - IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -INTEGER :: SideID, iLocSide, iNode, BCSideID, locElemID, CNElemID, iSide -REAL :: radius, triarea(2) -#if USE_MPI -REAL :: CNVolume -INTEGER :: offsetElemCNProc -#endif /*USE_MPI*/ -LOGICAL :: SymmetryBCExists -INTEGER :: firstElem, lastElem, firstSide, lastSide -!=================================================================================================================================== - -#if USE_MPI -CALL Allocate_Shared((/nNonUniqueGlobalSides/),SideIsSymSide_Shared_Win,SideIsSymSide_Shared) -CALL MPI_WIN_LOCK_ALL(0,SideIsSymSide_Shared_Win,IERROR) -SideIsSymSide => SideIsSymSide_Shared -! only CN root nullifies -IF(myComputeNodeRank.EQ.0) SideIsSymSide = .FALSE. -! This sync/barrier is required as it cannot be guaranteed that the zeros have been written to memory by the time the MPI_REDUCE -! is executed (see MPI specification). Until the Sync is complete, the status is undefined, i.e., old or new value or utter nonsense. -CALL BARRIER_AND_SYNC(SideIsSymSide_Shared_Win,MPI_COMM_SHARED) -#else -ALLOCATE(SideIsSymSide(nComputeNodeSides)) -SideIsSymSide = .FALSE. -#endif /*USE_MPI*/ - -! Flag of symmetry sides to be skipped during tracking -#if USE_MPI - firstSide = INT(REAL( myComputeNodeRank *nNonUniqueGlobalSides)/REAL(nComputeNodeProcessors))+1 - lastSide = INT(REAL((myComputeNodeRank+1)*nNonUniqueGlobalSides)/REAL(nComputeNodeProcessors)) -#else - firstSide = 1 - lastSide = nComputeNodeSides -#endif - -DO iSide = firstSide, lastSide - SideIsSymSide(iSide) = .FALSE. - ! ignore non-BC sides - IF (SideInfo_Shared(SIDE_BCID,iSide).LE.0) CYCLE -#if USE_MPI - ! ignore sides outside of halo region - IF (ElemInfo_Shared(ELEM_HALOFLAG,SideInfo_Shared(SIDE_ELEMID,iSide)).EQ.0) CYCLE -#endif /*USE_MPI*/ - IF (PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,iSide))).EQ.PartBound%SymmetryDim) & - SideIsSymSide(iSide) = .TRUE. -END DO - -SymmetryBCExists = .FALSE. -ALLOCATE(SymmetrySide(1:nElems,1:2)) ! 1: GlobalSide, 2: LocalSide -SymmetrySide = -1 - -! Sanity check: mesh has to be centered at z = 0 -IF(.NOT.ALMOSTEQUALRELATIVE(GEO%zmaxglob,ABS(GEO%zminglob),1e-5)) THEN - SWRITE(*,*) 'Maximum dimension in z:', GEO%zmaxglob - SWRITE(*,*) 'Minimum dimension in z:', GEO%zminglob - SWRITE(*,*) 'Deviation', (ABS(GEO%zmaxglob)-ABS(GEO%zminglob))/ABS(GEO%zminglob), ' > 1e-5' - CALL abort(__STAMP__& - ,'ERROR: Please orient your mesh with one cell in z-direction around 0, |z_min| = z_max !') -END IF - -! Calculation of the correct volume and characteristic length -DO BCSideID=1,nBCSides - locElemID = SideToElem(S2E_ELEM_ID,BCSideID) - iLocSide = SideToElem(S2E_LOC_SIDE_ID,BCSideID) - SideID=GetGlobalNonUniqueSideID(offsetElem+locElemID,iLocSide) - IF (PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).EQ.PartBound%SymmetryDim) THEN - CNElemID = GetCNElemID(SideInfo_Shared(SIDE_ELEMID,SideID)) - iLocSide = SideInfo_Shared(SIDE_LOCALID,SideID) - ! Exclude the symmetry axis (y=0) - IF(Symmetry%Axisymmetric) THEN - IF(MAXVAL(NodeCoords_Shared(2,ElemSideNodeID_Shared(:,iLocSide,CNElemID)+1)).LE.0.0) CYCLE - END IF - ! The z-plane with the positive z component is chosen - IF(MINVAL(NodeCoords_Shared(3,ElemSideNodeID_Shared(:,iLocSide,CNElemID)+1)).GT.(GEO%zmaxglob+GEO%zminglob)/2.) THEN - IF(SymmetrySide(locElemID,1).GT.0) THEN - CALL abort(__STAMP__& - ,'ERROR: PICLas could not determine a unique symmetry surface for 2D/axisymmetric calculation!'//& - ' Please orient your mesh with x as the symmetry axis and positive y as the second/radial direction!') - END IF - SymmetrySide(locElemID,1) = BCSideID - SymmetrySide(locElemID,2) = iLocSide - ! The volume calculated at this point (final volume for the 2D case) corresponds to the cell face area (z-dimension=1) in - ! the xy-plane. - ElemVolume_Shared(CNElemID) = 0.0 - CALL CalcNormAndTangTriangle(area=triarea(1),TriNum=1, SideID=SideID) - CALL CalcNormAndTangTriangle(area=triarea(2),TriNum=2, SideID=SideID) - ElemVolume_Shared(CNElemID) = triarea(1) + triarea(2) - ! Characteristic length is compared to the mean free path as the condition to refine the mesh. For the 2D/axisymmetric case - ! the third dimension is not considered as particle interaction occurs in the xy-plane, effectively reducing the refinement - ! requirement. - ElemCharLength_Shared(CNElemID) = SQRT(ElemVolume_Shared(CNElemID)) - ! Axisymmetric case: The volume is multiplied by the circumference to get the volume of the ring. The cell face in the - ! xy-plane is rotated around the x-axis. The radius is the middle point of the cell face. - IF (Symmetry%Axisymmetric) THEN - radius = 0. - DO iNode = 1, 4 - radius = radius + NodeCoords_Shared(2,ElemSideNodeID_Shared(iNode,iLocSide,CNElemID)+1) - END DO - radius = radius / 4. - ElemVolume_Shared(CNElemID) = ElemVolume_Shared(CNElemID) * 2. * Pi * radius - END IF - SymmetryBCExists = .TRUE. - END IF ! Greater z-coord - END IF -END DO - -IF(.NOT.SymmetryBCExists) THEN - CALL abort(__STAMP__& - ,'At least one symmetric BC (in the xy-plane) has to be defined for 2D simulations') -END IF - -! LocalVolume & MeshVolume: Recalculate the volume of the mesh of a single process and the total mesh volume -#if USE_MPI -FirstElem = offsetElem - offsetComputeNodeElem + 1 -LastElem = offsetElem - offsetComputeNodeElem + nElems -#else -firstElem = 1 -lastElem = nElems -#endif /*USE_MPI*/ - -LocalVolume = SUM(ElemVolume_Shared(FirstElem:LastElem)) - -#if USE_MPI -CALL BARRIER_AND_SYNC(SideIsSymSide_Shared_Win ,MPI_COMM_SHARED) -CALL BARRIER_AND_SYNC(ElemVolume_Shared_Win ,MPI_COMM_SHARED) -CALL BARRIER_AND_SYNC(ElemCharLength_Shared_Win,MPI_COMM_SHARED) -! Compute-node mesh volume -offsetElemCNProc = offsetElem - offsetComputeNodeElem -CNVolume = SUM(ElemVolume_Shared(offsetElemCNProc+1:offsetElemCNProc+nElems)) -CALL MPI_ALLREDUCE(MPI_IN_PLACE,CNVolume,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_SHARED,iError) -IF (myComputeNodeRank.EQ.0) THEN - ! All-reduce between node leaders - CALL MPI_ALLREDUCE(CNVolume,MeshVolume,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_LEADERS_SHARED,IERROR) -END IF -! Broadcast from node leaders to other processors on the same node -CALL MPI_BCAST(MeshVolume,1, MPI_DOUBLE_PRECISION,0,MPI_COMM_SHARED,iERROR) -#else -MeshVolume = LocalVolume -#endif /*USE_MPI*/ - -END SUBROUTINE DSMC_2D_InitVolumes - - -SUBROUTINE DSMC_1D_InitVolumes() -!=================================================================================================================================== -!> Routine determines a symmetry side and calculates the 1D (area faces at x axis) and axisymmetric volumes (cells are -!> revolved around the symmetry axis). The symmetry side will be used later on to determine in which two directions the quadtree -!> shall refine the mesh, skipping the z-dimension to avoid an unnecessary refinement. -!=================================================================================================================================== -! MODULES -USE MOD_Globals -USE MOD_Mesh_Vars ,ONLY: nElems, offsetElem -USE MOD_Particle_Boundary_Vars ,ONLY: PartBound -USE MOD_Particle_Mesh_Vars ,ONLY: GEO,LocalVolume,MeshVolume, SideIsSymSide -USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared,ElemCharLength_Shared -USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared, SideInfo_Shared -USE MOD_Mesh_Tools ,ONLY: GetCNElemID -USE MOD_Particle_Mesh_Tools ,ONLY: GetGlobalNonUniqueSideID -#if USE_MPI -USE MOD_MPI_Shared -USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED -USE MOD_Particle_Mesh_Vars ,ONLY: offsetComputeNodeElem,SideIsSymSide_Shared,SideIsSymSide_Shared_Win, nNonUniqueGlobalSides -USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared_Win,ElemCharLength_Shared_Win, ElemInfo_Shared -USE MOD_MPI_Shared_Vars ,ONLY: myComputeNodeRank,nComputeNodeProcessors,MPI_COMM_LEADERS_SHARED -#else -USE MOD_Particle_Mesh_Vars ,ONLY: nComputeNodeSides -#endif /*USE_MPI*/ -! IMPLICIT VARIABLE HANDLING - IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -INTEGER :: iLocSide, iElem, SideID, iOrder, CNElemID, iSide -REAL :: X(2), MaxCoord, MinCoord -LOGICAL :: SideInPlane, X1Occupied -INTEGER :: firstElem,lastElem, firstSide, lastSide -#if USE_MPI -REAL :: CNVolume -INTEGER :: offsetElemCNProc -#endif /*USE_MPI*/ -!=================================================================================================================================== - -#if USE_MPI -firstElem = offsetElem - offsetComputeNodeElem + 1 -lastElem = offsetElem - offsetComputeNodeElem + nElems -#else -firstElem = 1 -lastElem = nElems -#endif /*USE_MPI*/ - -ALLOCATE(GEO%XMinMax(2,nElems)) - -IF(.NOT.ALMOSTEQUALRELATIVE(GEO%ymaxglob,ABS(GEO%yminglob),1e-5)) THEN - SWRITE(*,*) 'Maximum dimension in y:', GEO%ymaxglob - SWRITE(*,*) 'Minimum dimension in y:', GEO%yminglob - SWRITE(*,*) 'Deviation', (ABS(GEO%ymaxglob)-ABS(GEO%yminglob))/ABS(GEO%yminglob), ' > 1e-5' - CALL abort(__STAMP__& - ,'ERROR: Please orient your mesh with one cell in y-direction around 0, |y_min| = y_max !') -END IF - -IF(.NOT.ALMOSTEQUALRELATIVE(GEO%zmaxglob,ABS(GEO%zminglob),1e-5)) THEN - SWRITE(*,*) 'Maximum dimension in z:', GEO%zmaxglob - SWRITE(*,*) 'Minimum dimension in z:', GEO%zminglob - SWRITE(*,*) 'Deviation', (ABS(GEO%zmaxglob)-ABS(GEO%zminglob))/ABS(GEO%zminglob), ' > 1e-5' - CALL abort(__STAMP__& - ,'ERROR: Please orient your mesh with one cell in z-direction around 0, |z_min| = z_max !') -END IF - -DO iElem = 1,nElems - ! Check if all sides of the element are parallel to xy-, xz-, or yz-plane and Sides parallel to xy-,and xz-plane are symmetric - ! And determine xmin and xmax of the element - X = 0 - X1Occupied = .FALSE. - DO iLocSide = 1,6 - SideInPlane=.FALSE. - SideID = GetGlobalNonUniqueSideID(offsetElem+iElem,iLocSide) - CNElemID = GetCNElemID(SideInfo_Shared(SIDE_ELEMID,SideID)) - DO iOrder=1,3 - MaxCoord = MAXVAL(NodeCoords_Shared(iOrder,ElemSideNodeID_Shared(:,iLocSide,CNElemID)+1)) - MinCoord = MINVAL(NodeCoords_Shared(iOrder,ElemSideNodeID_Shared(:,iLocSide,CNElemID)+1)) - IF(ALMOSTALMOSTEQUAL(MaxCoord,MinCoord).OR.(ALMOSTZERO(MaxCoord).AND.ALMOSTZERO(MinCoord))) THEN - IF(SideInPlane) CALL abort(__STAMP__& - ,'ERROR: Please orient your mesh with all element sides parallel to xy-,xz-,or yz-plane') - SideInPlane=.TRUE. - IF(iOrder.GE.2)THEN - IF(PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).NE.PartBound%SymmetryDim) & - CALL abort(__STAMP__,& - 'ERROR: Sides parallel to xy-,and xz-plane has to be the symmetric boundary condition') - END IF - IF(iOrder.EQ.1) THEN - IF(X1Occupied) THEN - X(2) = NodeCoords_Shared(1,ElemSideNodeID_Shared(1,iLocSide,CNElemID)+1) - ELSE - X(1) = NodeCoords_Shared(1,ElemSideNodeID_Shared(1,iLocSide,CNElemID)+1) - X1Occupied = .TRUE. - END IF - END IF !iOrder.EQ.1 - END IF - END DO !iOrder=1,3 - IF(.NOT.SideInPlane) THEN - IPWRITE(*,*) 'ElemID:',iElem,'SideID:',iLocSide - DO iOrder=1,4 - IPWRITE(*,*) 'Node',iOrder,'x:',NodeCoords_Shared(1,ElemSideNodeID_Shared(iOrder,iLocSide,CNElemID)+1), & - 'y:',NodeCoords_Shared(2,ElemSideNodeID_Shared(iOrder,iLocSide,CNElemID)+1), & - 'z:',NodeCoords_Shared(3,ElemSideNodeID_Shared(iOrder,iLocSide,CNElemID)+1) - END DO - CALL abort(__STAMP__& - ,'ERROR: Please orient your mesh with all element sides parallel to xy-,xz-,or yz-plane') - END IF - END DO ! iLocSide = 1,6 - GEO%XMinMax(1,iElem) = MINVAL(X) - GEO%XMinMax(2,iElem) = MAXVAL(X) - ElemVolume_Shared(CNElemID) = ABS(X(1)-X(2)) - ElemCharLength_Shared(CNElemID) = ElemVolume_Shared(CNElemID) -END DO -! LocalVolume & MeshVolume: Recalculate the volume of the mesh of a single process and the total mesh volume -LocalVolume = SUM(ElemVolume_Shared(firstElem:lastElem)) -#if USE_MPI -CALL BARRIER_AND_SYNC(ElemVolume_Shared_Win ,MPI_COMM_SHARED) -CALL BARRIER_AND_SYNC(ElemCharLength_Shared_Win,MPI_COMM_SHARED) -! Compute-node mesh volume -offsetElemCNProc = offsetElem - offsetComputeNodeElem -CNVolume = SUM(ElemVolume_Shared(offsetElemCNProc+1:offsetElemCNProc+nElems)) -CALL MPI_ALLREDUCE(MPI_IN_PLACE,CNVolume,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_SHARED,iError) -IF (myComputeNodeRank.EQ.0) THEN - ! All-reduce between node leaders - CALL MPI_ALLREDUCE(CNVolume,MeshVolume,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_LEADERS_SHARED,IERROR) -END IF -! Broadcast from node leaders to other processors on the same node -CALL MPI_BCAST(MeshVolume,1, MPI_DOUBLE_PRECISION,0,MPI_COMM_SHARED,iERROR) -#else -MeshVolume = SUM(ElemVolume_Shared) -#endif /*USE_MPI*/ - -#if USE_MPI -CALL Allocate_Shared((/nNonUniqueGlobalSides/),SideIsSymSide_Shared_Win,SideIsSymSide_Shared) -CALL MPI_WIN_LOCK_ALL(0,SideIsSymSide_Shared_Win,IERROR) -SideIsSymSide => SideIsSymSide_Shared -! only CN root nullifies -IF(myComputeNodeRank.EQ.0) SideIsSymSide = .FALSE. -! This sync/barrier is required as it cannot be guaranteed that the zeros have been written to memory by the time the MPI_REDUCE -! is executed (see MPI specification). Until the Sync is complete, the status is undefined, i.e., old or new value or utter nonsense. -CALL BARRIER_AND_SYNC(SideIsSymSide_Shared_Win,MPI_COMM_SHARED) -#else -ALLOCATE(SideIsSymSide(nComputeNodeSides)) -SideIsSymSide = .FALSE. -#endif /*USE_MPI*/ - -! Flag of symmetry sides to be skipped during tracking -#if USE_MPI - firstSide = INT(REAL( myComputeNodeRank *nNonUniqueGlobalSides)/REAL(nComputeNodeProcessors))+1 - lastSide = INT(REAL((myComputeNodeRank+1)*nNonUniqueGlobalSides)/REAL(nComputeNodeProcessors)) -#else - firstSide = 1 - lastSide = nComputeNodeSides -#endif - -DO iSide = firstSide, lastSide - SideIsSymSide(iSide) = .FALSE. - ! ignore non-BC sides - IF (SideInfo_Shared(SIDE_BCID,iSide).LE.0) CYCLE -#if USE_MPI - ! ignore sides outside of halo region - IF (ElemInfo_Shared(ELEM_HALOFLAG,SideInfo_Shared(SIDE_ELEMID,iSide)).EQ.0) CYCLE -#endif /*USE_MPI*/ - IF (PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,iSide))).EQ.PartBound%SymmetryDim) & - SideIsSymSide(iSide) = .TRUE. -END DO - -#if USE_MPI -CALL BARRIER_AND_SYNC(SideIsSymSide_Shared_Win ,MPI_COMM_SHARED) -#endif - -END SUBROUTINE DSMC_1D_InitVolumes - - SUBROUTINE DSMC_2D_InitRadialWeighting() !=================================================================================================================================== !> Read-in and initialize the variables required for the cloning procedures. Two modes with a delayed clone insertion are available: @@ -743,147 +383,6 @@ SUBROUTINE DSMC_2D_SetInClones() END SUBROUTINE DSMC_2D_SetInClones -REAL FUNCTION DSMC_2D_CalcSymmetryArea(iLocSide,iElem, ymin, ymax) -!=================================================================================================================================== -!> Calculates the actual area of an element for 2D simulations (plane/axisymmetric) regardless of the mesh dimension in z -!> Utilized in the particle emission (surface flux) and boundary sampling -!=================================================================================================================================== -! MODULES -USE MOD_Globals -USE MOD_Globals_Vars ,ONLY: Pi -USE MOD_Particle_Vars ,ONLY: Symmetry -USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared, ElemSideNodeID_Shared -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -INTEGER,INTENT(IN) :: iElem,iLocSide !> iElem is the compute-node element ID -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -REAL, OPTIONAL, INTENT(OUT) :: ymax,ymin -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -INTEGER :: iNode -REAL :: P(1:2,1:4), Pmin(2), Pmax(2), Length, MidPoint -!=================================================================================================================================== - -Pmin = HUGE(Pmin) -Pmax = -HUGE(Pmax) - -DO iNode = 1,4 - P(1:2,iNode) = NodeCoords_Shared(1:2,ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1) -END DO - -Pmax(1) = MAXVAL(P(1,:)) -Pmax(2) = MAXVAL(P(2,:)) -Pmin(1) = MINVAL(P(1,:)) -Pmin(2) = MINVAL(P(2,:)) - -IF (PRESENT(ymax).AND.PRESENT(ymin)) THEN - ymin = Pmin(2) - ymax = Pmax(2) -END IF - -Length = SQRT((Pmax(1)-Pmin(1))**2 + (Pmax(2)-Pmin(2))**2) - -MidPoint = (Pmax(2)+Pmin(2)) / 2. -IF(Symmetry%Axisymmetric) THEN - DSMC_2D_CalcSymmetryArea = Length * MidPoint * Pi * 2. - ! Area of the cells on the rotational symmetry axis is set to one - IF(.NOT.(DSMC_2D_CalcSymmetryArea.GT.0.0)) DSMC_2D_CalcSymmetryArea = 1. -ELSE - DSMC_2D_CalcSymmetryArea = Length -END IF -RETURN - -END FUNCTION DSMC_2D_CalcSymmetryArea - - -REAL FUNCTION DSMC_1D_CalcSymmetryArea(iLocSide,iElem) -!=================================================================================================================================== -!> Calculates the actual area of an element for 1D simulations regardless of the mesh dimension in z and y -!> Utilized in the particle emission (surface flux) and boundary sampling -!=================================================================================================================================== -! MODULES -USE MOD_Globals -USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared, ElemSideNodeID_Shared -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -INTEGER,INTENT(IN) :: iElem,iLocSide -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -REAL :: Pmin, Pmax, Length -!=================================================================================================================================== - -Pmax = MAXVAL(NodeCoords_Shared(1,ElemSideNodeID_Shared(:,iLocSide,iElem)+1)) -Pmin = MINVAL(NodeCoords_Shared(1,ElemSideNodeID_Shared(:,iLocSide,iElem)+1)) - -Length = ABS(Pmax-Pmin) - -! IF(Symmetry%Axisymmetric) THEN -! MidPoint = (Pmax(2)+Pmin(2)) / 2. -! DSMC_1D_CalcSymmetryArea = Length * MidPoint * Pi * 2. -! ! Area of the cells on the rotational symmetry axis is set to one -! IF(.NOT.(DSMC_1D_CalcSymmetryArea.GT.0.0)) DSMC_1D_CalcSymmetryArea = 1. -! ELSE - DSMC_1D_CalcSymmetryArea = Length - IF (DSMC_1D_CalcSymmetryArea.EQ.0.) DSMC_1D_CalcSymmetryArea = 1. -! END IF -RETURN - -END FUNCTION DSMC_1D_CalcSymmetryArea - - -FUNCTION DSMC_2D_CalcSymmetryAreaSubSides(iLocSide,iElem) -!=================================================================================================================================== -!> Calculates the area of the subsides for the insertion with the surface flux -!=================================================================================================================================== -! MODULES -USE MOD_Globals -USE MOD_Globals_Vars ,ONLY: Pi -USE MOD_DSMC_Vars ,ONLY: RadialWeighting -USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -INTEGER,INTENT(IN) :: iLocSide,iElem !> iElem is the compute-node element ID -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -REAL :: DSMC_2D_CalcSymmetryAreaSubSides(RadialWeighting%nSubSides) -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -INTEGER :: iNode -REAL :: P(1:2,1:4), Pmin(2), Pmax(2), MidPoint, PminTemp, PmaxTemp, Length -!=================================================================================================================================== - -Pmin = HUGE(Pmin) -Pmax = -HUGE(Pmax) - -DO iNode = 1,4 - P(1:2,iNode) = NodeCoords_Shared(1:2,ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1) -END DO - -Pmax(1) = MAXVAL(P(1,:)) -Pmax(2) = MAXVAL(P(2,:)) -Pmin(1) = MINVAL(P(1,:)) -Pmin(2) = MINVAL(P(2,:)) -Length = SQRT((Pmax(1)-Pmin(1))**2 + (Pmax(2)-Pmin(2))**2) - -DO iNode = 1, RadialWeighting%nSubSides - PminTemp = Pmin(2) + (Pmax(2) - Pmin(2))/RadialWeighting%nSubSides*(iNode-1.) - PmaxTemp = Pmin(2) + (Pmax(2) - Pmin(2))/RadialWeighting%nSubSides*iNode - MidPoint = (PmaxTemp+PminTemp) / 2. - DSMC_2D_CalcSymmetryAreaSubSides(iNode) = Length/RadialWeighting%nSubSides * MidPoint * Pi * 2. -END DO - -RETURN - -END FUNCTION DSMC_2D_CalcSymmetryAreaSubSides - - SUBROUTINE DSMC_2D_TreatIdenticalParticles(iPair, nPair, nPart, iElem, iPartIndx_Node) !=================================================================================================================================== !> Check if particle pairs have a zero relative velocity (and thus a collision probability of zero), if they do, break up the pair diff --git a/src/particles/dsmc/dsmc_vars.f90 b/src/particles/dsmc/dsmc_vars.f90 index 0328b71fa..1b1803457 100644 --- a/src/particles/dsmc/dsmc_vars.f90 +++ b/src/particles/dsmc/dsmc_vars.f90 @@ -505,7 +505,6 @@ MODULE MOD_DSMC_Vars REAL :: EElec ! Energy of Cell (Electronic State) END TYPE -INTEGER, ALLOCATABLE :: SymmetrySide(:,:) REAL,ALLOCATABLE :: DSMC_Solution(:,:,:) ! 1:3 v, 4:6 v^2, 7 dens, 8 Evib, 9 erot, 10 eelec REAL,ALLOCATABLE :: DSMC_SolutionPressTens(:,:) ! 1:3 pressure tensor, 4:6 heatflux diff --git a/src/particles/emission/particle_emission_init.f90 b/src/particles/emission/particle_emission_init.f90 index a56d47702..81797279f 100644 --- a/src/particles/emission/particle_emission_init.f90 +++ b/src/particles/emission/particle_emission_init.f90 @@ -149,6 +149,7 @@ SUBROUTINE InitializeVariablesSpeciesInits() USE MOD_Particle_Vars USE MOD_DSMC_Vars ,ONLY: useDSMC, BGGas USE MOD_DSMC_BGGas ,ONLY: BGGas_Initialize +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_MPI*/ @@ -821,9 +822,10 @@ 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: Species,nSpecies,SpecReset,Symmetry +USE MOD_Particle_Vars ,ONLY: Species,nSpecies,SpecReset USE MOD_ReadInTools USE MOD_Restart_Vars ,ONLY: DoRestart +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !---------------------------------------------------------------------------------------------------------------------------------- diff --git a/src/particles/emission/particle_emission_tools.f90 b/src/particles/emission/particle_emission_tools.f90 index f90c0d062..9210cea7c 100644 --- a/src/particles/emission/particle_emission_tools.f90 +++ b/src/particles/emission/particle_emission_tools.f90 @@ -1074,7 +1074,8 @@ SUBROUTINE SetParticlePositionCuboidCylinder(FractNbr,iInit,chunkSize,particle_p !=================================================================================================================================== ! modules USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Species, Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry +USE MOD_Particle_Vars ,ONLY: Species USE MOD_Part_Tools ,ONLY: CalcPartSymmetryPos, CalcRadWeightMPF USE MOD_DSMC_Vars ,ONLY: RadialWeighting !USE MOD_Particle_Mesh_Vars ,ONLY: GEO @@ -1166,9 +1167,10 @@ SUBROUTINE SetParticlePositionSphere(FractNbr,iInit,chunkSize,particle_positions !=================================================================================================================================== ! modules USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Species, Symmetry -USE MOD_Part_tools ,ONLY: DICEUNITVECTOR, CalcPartSymmetryPos, CalcRadWeightMPF -USE MOD_DSMC_Vars ,ONLY: RadialWeighting +USE MOD_Particle_Vars ,ONLY: Species +USE MOD_Part_tools ,ONLY: DICEUNITVECTOR, CalcPartSymmetryPos, CalcRadWeightMPF +USE MOD_DSMC_Vars ,ONLY: RadialWeighting +USE MOD_Symmetry_Vars ,ONLY: Symmetry !---------------------------------------------------------------------------------------------------------------------------------- ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE diff --git a/src/particles/emission/particle_macroscopic_restart.f90 b/src/particles/emission/particle_macroscopic_restart.f90 index b2e05118e..eb0a32f2a 100644 --- a/src/particles/emission/particle_macroscopic_restart.f90 +++ b/src/particles/emission/particle_macroscopic_restart.f90 @@ -42,10 +42,11 @@ SUBROUTINE MacroRestart_InsertParticles() 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_Particle_Vars ,ONLY: Species, PDM, nSpecies, PartState, UseVarTimeStep USE MOD_Restart_Vars ,ONLY: MacroRestartValues USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared,BoundsOfElem_Shared, ElemMidPoint_Shared USE MOD_Particle_Tracking ,ONLY: ParticleInsideCheck +USE MOD_Symmetry_Vars ,ONLY: Symmetry !----------------------------------------------------------------------------------------------------------------------------------- ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE diff --git a/src/particles/emission/particle_position_and_velocity.f90 b/src/particles/emission/particle_position_and_velocity.f90 index a2c958b15..35a5a5211 100644 --- a/src/particles/emission/particle_position_and_velocity.f90 +++ b/src/particles/emission/particle_position_and_velocity.f90 @@ -44,13 +44,15 @@ SUBROUTINE ParticleEmissionCellLocal(iSpec,iInit,NbrOfParticle) USE MOD_Particle_Mesh_Vars ,ONLY: BoundsOfElem_Shared,ElemVolume_Shared,ElemMidPoint_Shared USE MOD_Mesh_Tools ,ONLY: GetCNElemID USE MOD_Particle_Tracking ,ONLY: ParticleInsideCheck -USE MOD_Particle_Vars ,ONLY: Species, PDM, PartState, PEM, Symmetry, UseVarTimeStep, PartTimeStep, PartMPF, PartSpecies +USE MOD_Particle_Vars ,ONLY: Species, PDM, PartState, PEM, UseVarTimeStep, PartTimeStep, PartMPF, PartSpecies USE MOD_Particle_Vars ,ONLY: usevMPF, UseSplitAndMerge, vMPFSplitThreshold USE MOD_Particle_TimeStep ,ONLY: GetParticleTimeStep USE MOD_Part_Tools ,ONLY: IncreaseMaxParticleNumber, GetNextFreePosition +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_Particle_MPI_Vars ,ONLY: PartMPIInitGroup #endif /*USE_MPI*/ +!---------------------------------------------------------------------------------------------------------------------------------- ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- diff --git a/src/particles/emission/particle_surface_chemflux.f90 b/src/particles/emission/particle_surface_chemflux.f90 index 5508f794c..1a5df7aa1 100644 --- a/src/particles/emission/particle_surface_chemflux.f90 +++ b/src/particles/emission/particle_surface_chemflux.f90 @@ -47,7 +47,6 @@ SUBROUTINE ParticleSurfChemFlux() USE MOD_Mesh_Vars ,ONLY: SideToElem, offsetElem USE MOD_Particle_Mesh_Vars ,ONLY: ElemMidPoint_Shared USE MOD_Mesh_Tools ,ONLY: GetCNElemID -USE MOD_Part_Emission_Tools ,ONLY: SetParticleMPF USE MOD_Particle_Analyze_Vars ,ONLY: CalcPartBalance, nPartIn, PartEkinIn USE MOD_Particle_Analyze_Tools ,ONLY: CalcEkinPart USE MOD_Particle_Mesh_Tools ,ONLY: GetGlobalNonUniqueSideID @@ -57,6 +56,7 @@ SUBROUTINE ParticleSurfChemFlux() USE MOD_SurfaceModel_Vars ,ONLY: SurfChem, SurfChemReac, ChemWallProp, ChemDesorpWall USE MOD_Particle_SurfFlux ,ONLY: CalcPartPosTriaSurface, DefineSideDirectVec2D USE MOD_DSMC_PolyAtomicModel ,ONLY: DSMC_SetInternalEnr +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_MPI_Shared_vars ,ONLY: MPI_COMM_SHARED USE MOD_MPI_Shared ,ONLY: BARRIER_AND_SYNC diff --git a/src/particles/emission/particle_surface_flux.f90 b/src/particles/emission/particle_surface_flux.f90 index d4cad3300..2302c493d 100644 --- a/src/particles/emission/particle_surface_flux.f90 +++ b/src/particles/emission/particle_surface_flux.f90 @@ -48,6 +48,7 @@ SUBROUTINE ParticleSurfaceflux() USE MOD_Particle_Surfaces_Vars ,ONLY: SurfFluxSideSize, TriaSurfaceFlux, BCdata_auxSF USE MOD_Particle_TimeStep ,ONLY: GetParticleTimeStep USE MOD_Timedisc_Vars ,ONLY: RKdtFrac, dt +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_DSMC_PolyAtomicModel ,ONLY: DSMC_SetInternalEnr #if defined(IMPA) || defined(ROS) USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod diff --git a/src/particles/emission/particle_surface_flux_init.f90 b/src/particles/emission/particle_surface_flux_init.f90 index 68329a5dc..806db7e1f 100644 --- a/src/particles/emission/particle_surface_flux_init.f90 +++ b/src/particles/emission/particle_surface_flux_init.f90 @@ -139,7 +139,7 @@ SUBROUTINE InitializeParticleSurfaceflux() !USE MOD_Particle_SurfChemFlux_Init USE MOD_Restart_Vars ,ONLY: DoRestart, RestartTime USE MOD_DSMC_Vars ,ONLY: AmbiPolarSFMapping, DSMC, useDSMC -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_Particle_Vars ,ONLY: DoPoissonRounding, DoTimeDepInflow #endif /*USE_MPI*/ @@ -702,6 +702,7 @@ SUBROUTINE BCSurfMeshSideAreasandNormals() #endif /*CODE_ANALYZE*/ END SUBROUTINE BCSurfMeshSideAreasandNormals + SUBROUTINE CreateSideListAndFinalizeAreasSurfFlux(nDataBC, BCdata_auxSFTemp) !=================================================================================================================================== ! SideList for SurfaceFlux in BCdata_auxSF is created. Furthermore, the side areas are corrected for the 1D/2D case and finally @@ -719,10 +720,11 @@ SUBROUTINE CreateSideListAndFinalizeAreasSurfFlux(nDataBC, BCdata_auxSFTemp) USE MOD_Mesh_Vars ,ONLY: nBCSides, offsetElem, BC, SideToElem USE MOD_Particle_Mesh_Vars ,ONLY: GEO, ElemMidPoint_Shared, SideInfo_Shared USE MOD_Mesh_Tools ,ONLY: GetCNElemID -USE MOD_Particle_Vars ,ONLY: UseCircularInflow, Species, DoSurfaceFlux, nSpecies, Symmetry -USE MOD_DSMC_Symmetry ,ONLY: DSMC_1D_CalcSymmetryArea, DSMC_2D_CalcSymmetryArea, DSMC_2D_CalcSymmetryAreaSubSides +USE MOD_Particle_Vars ,ONLY: UseCircularInflow, Species, DoSurfaceFlux, nSpecies +USE MOD_Particle_Mesh_Tools ,ONLY: DSMC_1D_CalcSymmetryArea, DSMC_2D_CalcSymmetryArea, DSMC_2D_CalcSymmetryAreaSubSides USE MOD_DSMC_Vars ,ONLY: RadialWeighting USE MOD_Particle_Surfaces ,ONLY: CalcNormAndTangTriangle +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- diff --git a/src/particles/fp_flow/fpflow_init.f90 b/src/particles/fp_flow/fpflow_init.f90 index 086244f10..a821bc0c2 100644 --- a/src/particles/fp_flow/fpflow_init.f90 +++ b/src/particles/fp_flow/fpflow_init.f90 @@ -72,7 +72,8 @@ SUBROUTINE InitFPFlow() USE MOD_ReadInTools USE MOD_DSMC_Vars ,ONLY: DSMC, CollInf USE MOD_DSMC_ParticlePairing ,ONLY: DSMC_init_octree -USE MOD_PARTICLE_Vars ,ONLY: nSpecies, Species, DoVirtualCellMerge,Symmetry +USE MOD_PARTICLE_Vars ,ONLY: nSpecies, Species, DoVirtualCellMerge +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_FPFlow_Vars USE MOD_BGK_Vars ,ONLY: DoBGKCellAdaptation, BGKMinPartPerCell #if USE_LOADBALANCE @@ -109,9 +110,7 @@ SUBROUTINE InitFPFlow() BGKMinPartPerCell = GETINT('Particles-FP-MinPartsPerCell') IF(.NOT.DSMC%UseOctree) THEN DSMC%UseOctree = .TRUE. - IF(NGeo.GT.PP_N) CALL abort(& -__STAMP__& -,' Set PP_N to NGeo, otherwise the volume is not computed correctly.') + IF(NGeo.GT.PP_N) CALL abort(__STAMP__,' Set PP_N to NGeo, otherwise the volume is not computed correctly.') CALL DSMC_init_octree() END IF END IF diff --git a/src/particles/fp_flow/fpflow_main.f90 b/src/particles/fp_flow/fpflow_main.f90 index 1fc481e17..0798eb468 100644 --- a/src/particles/fp_flow/fpflow_main.f90 +++ b/src/particles/fp_flow/fpflow_main.f90 @@ -47,7 +47,7 @@ SUBROUTINE FP_DSMC_main() USE MOD_Globals USE MOD_TimeDisc_Vars ,ONLY: TEnd, Time USE MOD_Mesh_Vars ,ONLY: nElems, offsetElem -USE MOD_Particle_Vars ,ONLY: PEM, Species, WriteMacroVolumeValues, Symmetry, usevMPF +USE MOD_Particle_Vars ,ONLY: PEM, Species, WriteMacroVolumeValues, usevMPF USE MOD_FP_CollOperator ,ONLY: FP_CollisionOperator USE MOD_FPFlow_Vars ,ONLY: FPDSMCSwitchDens, FP_QualityFacSamp, FP_PrandtlNumber USE MOD_FPFlow_Vars ,ONLY: FP_MaxRelaxFactor, FP_MaxRotRelaxFactor, FP_MeanRelaxFactor, FP_MeanRelaxFactorCounter @@ -58,6 +58,7 @@ SUBROUTINE FP_DSMC_main() USE MOD_Part_Tools ,ONLY: GetParticleWeight USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared USE MOD_Mesh_Tools ,ONLY: GetCNElemID +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -143,7 +144,7 @@ SUBROUTINE FPFlow_main() USE MOD_Globals USE MOD_TimeDisc_Vars ,ONLY: TEnd, Time USE MOD_Mesh_Vars ,ONLY: nElems, offsetElem -USE MOD_Particle_Vars ,ONLY: PEM, WriteMacroVolumeValues, WriteMacroSurfaceValues, Symmetry, DoVirtualCellMerge, VirtMergedCells +USE MOD_Particle_Vars ,ONLY: PEM, WriteMacroVolumeValues, WriteMacroSurfaceValues, DoVirtualCellMerge, VirtMergedCells USE MOD_FP_CollOperator ,ONLY: FP_CollisionOperator USE MOD_DSMC_Vars ,ONLY: DSMC USE MOD_BGK_Vars ,ONLY: DoBGKCellAdaptation @@ -153,6 +154,7 @@ SUBROUTINE FPFlow_main() USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared USE MOD_Mesh_Tools ,ONLY: GetCNElemID USE MOD_DSMC_Analyze ,ONLY: DSMCMacroSampling +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- diff --git a/src/particles/particle_init.f90 b/src/particles/particle_init.f90 index 6c8e1373c..1c82c05e4 100644 --- a/src/particles/particle_init.f90 +++ b/src/particles/particle_init.f90 @@ -60,7 +60,6 @@ END FUNCTION GetPID_C PUBLIC::FinalizeParticles PUBLIC::DefineParametersParticles PUBLIC::InitialIonization -PUBLIC::InitSymmetry !=================================================================================================================================== CONTAINS @@ -237,7 +236,7 @@ SUBROUTINE InitParticleGlobals(IsLoadBalance) USE MOD_Globals_Vars ,ONLY: ProjectName USE MOD_ReadInTools USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_LOADBALANCE*/ @@ -350,11 +349,11 @@ SUBROUTINE InitParticles() #endif #if (PP_TimeDiscMethod==300) USE MOD_FPFlow_Init ,ONLY: InitFPFlow -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #endif #if (PP_TimeDiscMethod==400) USE MOD_BGK_Init ,ONLY: InitBGK -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #endif USE MOD_Particle_Vars ,ONLY: BulkElectronTemp #if USE_LOADBALANCE @@ -483,7 +482,7 @@ SUBROUTINE InitializeVariables() USE MOD_Globals USE MOD_ReadInTools USE MOD_Particle_Vars -USE MOD_DSMC_Symmetry ,ONLY: DSMC_1D_InitVolumes, DSMC_2D_InitVolumes, DSMC_2D_InitRadialWeighting +USE MOD_DSMC_Symmetry ,ONLY: DSMC_2D_InitRadialWeighting USE MOD_DSMC_Vars ,ONLY: RadialWeighting USE MOD_Part_RHS ,ONLY: InitPartRHS USE MOD_Particle_Mesh ,ONLY: InitParticleMesh @@ -495,6 +494,7 @@ SUBROUTINE InitializeVariables() USE MOD_PICInit ,ONLY: InitPIC USE MOD_PICDepo_Vars ,ONLY: DoDeposition USE MOD_PICInterpolation_Vars ,ONLY: DoInterpolation +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_SurfaceModel_Vars ,ONLY: SurfChem #if USE_MPI USE MOD_Particle_MPI_Emission ,ONLY: InitEmissionComm @@ -586,9 +586,6 @@ SUBROUTINE InitializeVariables() CALL InitPIC() ! === 2D/1D/Axisymmetric initialization -! Calculate the volumes for 2D simulation (requires the GEO%zminglob/GEO%zmaxglob from InitFIBGM) -IF(Symmetry%Order.EQ.2) CALL DSMC_2D_InitVolumes() -IF(Symmetry%Order.EQ.1) CALL DSMC_1D_InitVolumes() IF(Symmetry%Axisymmetric) THEN IF(RadialWeighting%DoRadialWeighting) THEN ! Initialization of RadialWeighting in 2D axisymmetric simulations @@ -802,6 +799,7 @@ SUBROUTINE InitializeVariablesVarTimeStep() USE MOD_Particle_Vars USE MOD_Mesh_Vars ,ONLY: nElems USE MOD_Particle_Tracking_Vars ,ONLY: TrackingMethod +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Particle_TimeStep ,ONLY: VarTimeStep_CalcElemFacs ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -1873,51 +1871,5 @@ SUBROUTINE InitPartDataSize() END SUBROUTINE InitPartDataSize -SUBROUTINE InitSymmetry() -!=================================================================================================================================== -!> Initialize if a 2D/1D Simulation is performed and which type -!=================================================================================================================================== -! MODULES -USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Symmetry -USE MOD_DSMC_Vars ,ONLY: RadialWeighting -USE MOD_ReadInTools ,ONLY: GETLOGICAL,GETINT -USE MOD_Particle_Mesh_Tools ,ONLY: InitParticleInsideQuad -USE MOD_Particle_TriaTracking ,ONLY: InitSingleParticleTriaTracking -USE MOD_Particle_InterSection ,ONLY: InitParticleThroughSideCheck1D2D -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -!=================================================================================================================================== -Symmetry%Order = GETINT('Particles-Symmetry-Order') - -IF((Symmetry%Order.LE.0).OR.(Symmetry%Order.GE.4)) CALL ABORT(__STAMP__& - ,'Particles-Symmetry-Order (space dimension) has to be in the range of 1 to 3') - -! Initialize the function pointers for triatracking -CALL InitParticleInsideQuad() -CALL InitSingleParticleTriaTracking() -IF(Symmetry%Order.LE.2) CALL InitParticleThroughSideCheck1D2D() - -Symmetry%Axisymmetric = GETLOGICAL('Particles-Symmetry2DAxisymmetric') -IF(Symmetry%Axisymmetric.AND.(Symmetry%Order.EQ.3)) CALL ABORT(__STAMP__& - ,'ERROR: Axisymmetric simulations only for 1D or 2D') -IF(Symmetry%Axisymmetric.AND.(Symmetry%Order.EQ.1))CALL ABORT(__STAMP__& - ,'ERROR: Axisymmetric simulations are only implemented for Particles-Symmetry-Order=2 !') -IF(Symmetry%Axisymmetric) THEN - RadialWeighting%DoRadialWeighting = GETLOGICAL('Particles-RadialWeighting') -ELSE - RadialWeighting%DoRadialWeighting = .FALSE. - RadialWeighting%PerformCloning = .FALSE. -END IF - -END SUBROUTINE InitSymmetry - END MODULE MOD_ParticleInit diff --git a/src/particles/particle_mesh/particle_mesh.f90 b/src/particles/particle_mesh/particle_mesh.f90 index 6e1b40425..15dfe9fe1 100644 --- a/src/particles/particle_mesh/particle_mesh.f90 +++ b/src/particles/particle_mesh/particle_mesh.f90 @@ -171,7 +171,8 @@ SUBROUTINE InitParticleMesh() USE MOD_PICInterpolation_Vars ,ONLY: DoInterpolation USE MOD_PICDepo_Vars ,ONLY: DoDeposition,DepositionType USE MOD_ReadInTools ,ONLY: GETREAL,GETINT,GETLOGICAL,GetRealArray, GETINTFROMSTR -USE MOD_Particle_Vars ,ONLY: Symmetry, DoVirtualCellMerge +USE MOD_Symmetry_Vars ,ONLY: Symmetry +USE MOD_Particle_Vars ,ONLY: DoVirtualCellMerge USE MOD_Particle_Boundary_Vars ,ONLY: PartBound #ifdef CODE_ANALYZE !USE MOD_Particle_Surfaces_Vars ,ONLY: SideBoundingBoxVolume @@ -199,6 +200,7 @@ SUBROUTINE InitParticleMesh() USE MOD_Photon_TrackingVars ,ONLY: PhotonModeBPO,UsePhotonTriaTracking USE MOD_RayTracing_Vars ,ONLY: UseRayTracing USE MOD_RayTracing_Vars ,ONLY: PerformRayTracing +USE MOD_Particle_Mesh_Tools ,ONLY: InitVolumes_1D,InitVolumes_2D !USE MOD_DSMC_Vars ,ONLY: DSMC ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -348,7 +350,7 @@ SUBROUTINE InitParticleMesh() END IF WRITE(tmpStr,'(I2.2)') RefMappingGuessProposal RefMappingGuess = GETINT('RefMappingGuess',tmpStr) -IF((RefMappingGuess.LT.1).AND.(UseCurveds)) THEN ! Linear intial guess on curved meshes might cause problems. +IF((RefMappingGuess.LT.1).AND.(UseCurveds)) THEN ! Linear initial guess on curved meshes might cause problems. LBWRITE(UNIT_stdOut,'(A)')' WARNING: read-in [RefMappingGuess=1] when using [UseCurveds=T] may create problems!' END IF RefMappingEps = GETREAL('RefMappingEps','1e-4') @@ -456,6 +458,10 @@ SUBROUTINE InitParticleMesh() ! Build mappings UniqueNodeID->CN Element IDs and CN Element ID -> CN Element IDs IF(FindNeighbourElems) CALL BuildNodeNeighbourhood() +! Calculate the volumes for 2D/1D simulation (requires the GEO%zminglob/GEO%zmaxglob from InitFIBGM) +IF(Symmetry%Order.EQ.2) CALL InitVolumes_2D() +IF(Symmetry%Order.EQ.1) CALL InitVolumes_1D() + ! BezierAreaSample stuff: IF (TriaSurfaceFlux) THEN BezierSampleN = 1 @@ -489,7 +495,7 @@ SUBROUTINE FinalizeParticleMesh() ! MODULES USE MOD_Globals USE MOD_Particle_Mesh_Vars -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_RayTracing_Vars ,ONLY: UseRayTracing USE MOD_Particle_Surfaces_Vars ,ONLY: BezierElevation @@ -823,7 +829,7 @@ SUBROUTINE FinalizeParticleMesh() ADEALLOCATE(ElemSideNodeID2D_Shared) ADEALLOCATE(SideNormalEdge2D_Shared) ELSE IF (Symmetry%Order.EQ.1) THEN - ADEALLOCATE(ElemSideNodeID1D_Shared) + ADEALLOCATE(ElemSideNodeID1D_Shared) END IF ! BuildEpsOneCell @@ -899,6 +905,7 @@ SUBROUTINE FinalizeParticleMesh() ADEALLOCATE(ConcaveElemSide_Shared) ADEALLOCATE(ElemSideNodeID_Shared) ADEALLOCATE(ElemMidPoint_Shared) +SDEALLOCATE(SymmetrySide) ! Load Balance !#if !USE_LOADBALANCE diff --git a/src/particles/particle_mesh/particle_mesh_build.f90 b/src/particles/particle_mesh/particle_mesh_build.f90 index d6ed221ed..9bc2e0765 100644 --- a/src/particles/particle_mesh/particle_mesh_build.f90 +++ b/src/particles/particle_mesh/particle_mesh_build.f90 @@ -35,9 +35,7 @@ MODULE MOD_Particle_Mesh_Build SUBROUTINE BuildMesh2DInfo() !=================================================================================================================================== -!> Routine determines a symmetry side and calculates the 2D (area faces in symmetry plane) and axisymmetric volumes (cells are -!> revolved around the symmetry axis). The symmetry side will be used later on to determine in which two directions the quadtree -!> shall refine the mesh, skipping the z-dimension to avoid an unnecessary refinement. +!> Build the ElemSideNodeID and SideNormalEdge array for 2D mesh !=================================================================================================================================== ! MODULES USE MOD_Globals @@ -49,7 +47,7 @@ SUBROUTINE BuildMesh2DInfo() USE MOD_Mesh_Tools ,ONLY: GetGlobalElemID USE MOD_Particle_Mesh_Tools ,ONLY: GetGlobalNonUniqueSideID #if USE_MPI -USE MOD_MPI_Shared +USE MOD_MPI_Shared USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED, nComputeNodeTotalElems USE MOD_Particle_Mesh_Vars ,ONLY: ElemSideNodeID2D_Shared_Win, SideNormalEdge2D_Shared_Win USE MOD_MPI_Shared_Vars ,ONLY: myComputeNodeRank, nComputeNodeProcessors @@ -87,9 +85,10 @@ SUBROUTINE BuildMesh2DInfo() DO iElem = firstElem, lastElem GlobalElemID = GetGlobalElemID(iElem) - DO iLocSide = 1, 6 + DO iLocSide = 1, 6 DefineSide = .TRUE. SideID=GetGlobalNonUniqueSideID(GlobalElemID,iLocSide) + ! Check if side is a symmetry_dim boundary IF (SideInfo_Shared(SIDE_BCID,SideID).GT.0) THEN IF (PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).EQ.PartBound%SymmetryDim) THEN ElemSideNodeID2D_Shared(:,iLocSide, iElem) = -1 @@ -97,18 +96,18 @@ SUBROUTINE BuildMesh2DInfo() DefineSide = .FALSE. END IF END IF - + IF (DefineSide) THEN - tmpNode = 0 + tmpNode = 0 DO iNode = 1, 4 IF(NodeCoords_Shared(3,ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1).GT.(GEO%zmaxglob+GEO%zminglob)/2.) THEN tmpNode = tmpNode + 1 ElemSideNodeID2D_Shared(tmpNode,iLocSide, iElem) = ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1 END IF - END DO - EdgeVec(1:2) = NodeCoords_Shared(1:2,ElemSideNodeID2D_Shared(2,ilocSide,iElem))-NodeCoords_Shared(1:2,ElemSideNodeID2D_Shared(1,ilocSide,iElem)) + END DO + EdgeVec(1:2) = NodeCoords_Shared(1:2,ElemSideNodeID2D_Shared(2,iLocSide,iElem))-NodeCoords_Shared(1:2,ElemSideNodeID2D_Shared(1,iLocSide,iElem)) NormVec(1) = -EdgeVec(2) - NormVec(2) = EdgeVec(1) + NormVec(2) = EdgeVec(1) FaceMidPoint(1:2) = (NodeCoords_Shared(1:2,ElemSideNodeID_Shared(1,iLocSide,iElem)+1) & + NodeCoords_Shared(1:2,ElemSideNodeID_Shared(2,iLocSide,iElem)+1)& + NodeCoords_Shared(1:2,ElemSideNodeID_Shared(3,iLocSide,iElem)+1) & @@ -139,9 +138,7 @@ END SUBROUTINE BuildMesh2DInfo SUBROUTINE BuildMesh1DInfo() !=================================================================================================================================== -!> Routine determines a symmetry side and calculates the 1D (area faces in symmetry plane) -!> The symmetry sides will be used later on to determine in which direction the quadtree -!> shall refine the mesh, skipping the z- and y-dimension to avoid an unnecessary refinement. +!> Build the ElemSideNodeID array for 1D mesh !=================================================================================================================================== ! MODULES USE MOD_Globals @@ -152,7 +149,7 @@ SUBROUTINE BuildMesh1DInfo() USE MOD_Mesh_Tools ,ONLY: GetGlobalElemID USE MOD_Particle_Mesh_Tools ,ONLY: GetGlobalNonUniqueSideID #if USE_MPI -USE MOD_MPI_Shared +USE MOD_MPI_Shared USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED, nComputeNodeTotalElems USE MOD_Particle_Mesh_Vars ,ONLY: ElemSideNodeID1D_Shared_Win USE MOD_MPI_Shared_Vars ,ONLY: myComputeNodeRank, nComputeNodeProcessors @@ -186,7 +183,7 @@ SUBROUTINE BuildMesh1DInfo() DO iElem = firstElem, lastElem GlobalElemID = GetGlobalElemID(iElem) - DO iLocSide = 1, 6 + DO iLocSide = 1, 6 DefineSide = .TRUE. SideID=GetGlobalNonUniqueSideID(GlobalElemID,iLocSide) IF (SideInfo_Shared(SIDE_BCID,SideID).GT.0) THEN @@ -195,14 +192,14 @@ SUBROUTINE BuildMesh1DInfo() DefineSide = .FALSE. END IF END IF - + IF (DefineSide) THEN DO iNode = 1, 4 IF((NodeCoords_Shared(3,ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1).GT.(GEO%zmaxglob+GEO%zminglob)/2.).AND.& (NodeCoords_Shared(2,ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1).GT.(GEO%ymaxglob+GEO%yminglob)/2.)) THEN ElemSideNodeID1D_Shared(iLocSide, iElem) = ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1 END IF - END DO + END DO END IF END DO END DO diff --git a/src/particles/particle_mesh/particle_mesh_tools.f90 b/src/particles/particle_mesh/particle_mesh_tools.f90 index 99705792c..7ce3f9cd9 100644 --- a/src/particles/particle_mesh/particle_mesh_tools.f90 +++ b/src/particles/particle_mesh/particle_mesh_tools.f90 @@ -45,8 +45,9 @@ SUBROUTINE ParticleInsideQuadInterface(PartStateLoc,ElemID,InElementCheck,Det_Ou PUBLIC :: ParticleInsideQuad3D, InitPEM_LocalElemID, InitPEM_CNElemID, GetGlobalNonUniqueSideID, GetSideBoundingBoxTria PUBLIC :: GetMeshMinMax, IdentifyElemAndSideType, WeirdElementCheck, CalcParticleMeshMetrics, CalcXCL_NGeo -PUBLIC :: CalcBezierControlPoints, InitParticleGeometry, ComputePeriodicVec, ParticleInsideQuad2D, ParticleInsideQuad -PUBLIC :: InitParticleInsideQuad +PUBLIC :: CalcBezierControlPoints, InitParticleGeometry, ComputePeriodicVec +PUBLIC :: InitParticleInsideQuad, ParticleInsideQuad +PUBLIC :: InitVolumes_2D, InitVolumes_1D, DSMC_2D_CalcSymmetryArea, DSMC_1D_CalcSymmetryArea, DSMC_2D_CalcSymmetryAreaSubSides !=================================================================================================================================== CONTAINS @@ -56,7 +57,7 @@ SUBROUTINE ParticleInsideQuadInterface(PartStateLoc,ElemID,InElementCheck,Det_Ou SUBROUTINE InitParticleInsideQuad() ! MODULES USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !================================================================================================================================== @@ -2052,4 +2053,494 @@ SUBROUTINE ComputePeriodicVec() END SUBROUTINE ComputePeriodicVec +SUBROUTINE InitVolumes_2D() +!=================================================================================================================================== +!> Routine determines the symmetry sides and calculates the 2D (area faces in symmetry plane) and axisymmetric volumes (cells are +!> revolved around the symmetry axis). The symmetry side (SymmetrySide array) will be used later on to determine in which two +!> directions the quadtree shall refine the mesh, skipping the z-dimension to avoid an unnecessary refinement. Additionally, +!> symmetry sides will be skipped during tracking (SideIsSymSide_Shared array). +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Globals_Vars ,ONLY: Pi +USE MOD_PreProc +USE MOD_Mesh_Vars ,ONLY: nElems,offsetElem,nBCSides,SideToElem +USE MOD_Particle_Boundary_Vars ,ONLY: PartBound +USE MOD_Particle_Mesh_Vars ,ONLY: GEO,LocalVolume,MeshVolume +USE MOD_Particle_Mesh_Vars ,ONLY: SymmetrySide +USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared,ElemCharLength_Shared +USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared, SideInfo_Shared, SideIsSymSide +USE MOD_Mesh_Tools ,ONLY: GetCNElemID +USE MOD_Particle_Surfaces ,ONLY: CalcNormAndTangTriangle +#if USE_MPI +USE MOD_MPI_Shared +USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED +USE MOD_Particle_Mesh_Vars ,ONLY: nNonUniqueGlobalSides, offsetComputeNodeElem, ElemInfo_Shared +USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared_Win,ElemCharLength_Shared_Win,SideIsSymSide_Shared,SideIsSymSide_Shared_Win +USE MOD_MPI_Shared_Vars ,ONLY: myComputeNodeRank,nComputeNodeProcessors,MPI_COMM_LEADERS_SHARED +#else +USE MOD_Particle_Mesh_Vars ,ONLY: nComputeNodeSides +#endif /*USE_MPI*/ +USE MOD_Symmetry_Vars ,ONLY: Symmetry +! IMPLICIT VARIABLE HANDLING + IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: SideID, iLocSide, iNode, BCSideID, locElemID, CNElemID, iSide +REAL :: radius, triarea(2) +#if USE_MPI +REAL :: CNVolume +INTEGER :: offsetElemCNProc +#endif /*USE_MPI*/ +LOGICAL :: SymmetryBCExists +INTEGER :: firstElem, lastElem, firstSide, lastSide +!=================================================================================================================================== + +#if USE_MPI +CALL Allocate_Shared((/nNonUniqueGlobalSides/),SideIsSymSide_Shared_Win,SideIsSymSide_Shared) +CALL MPI_WIN_LOCK_ALL(0,SideIsSymSide_Shared_Win,IERROR) +SideIsSymSide => SideIsSymSide_Shared +! only CN root nullifies +IF(myComputeNodeRank.EQ.0) SideIsSymSide = .FALSE. +! This sync/barrier is required as it cannot be guaranteed that the zeros have been written to memory by the time the MPI_REDUCE +! is executed (see MPI specification). Until the Sync is complete, the status is undefined, i.e., old or new value or utter nonsense. +CALL BARRIER_AND_SYNC(SideIsSymSide_Shared_Win,MPI_COMM_SHARED) +#else +ALLOCATE(SideIsSymSide(nComputeNodeSides)) +SideIsSymSide = .FALSE. +#endif /*USE_MPI*/ + +! Flag of symmetry sides to be skipped during tracking +#if USE_MPI + firstSide = INT(REAL( myComputeNodeRank *nNonUniqueGlobalSides)/REAL(nComputeNodeProcessors))+1 + lastSide = INT(REAL((myComputeNodeRank+1)*nNonUniqueGlobalSides)/REAL(nComputeNodeProcessors)) +#else + firstSide = 1 + lastSide = nComputeNodeSides +#endif + +DO iSide = firstSide, lastSide + SideIsSymSide(iSide) = .FALSE. + ! ignore non-BC sides + IF (SideInfo_Shared(SIDE_BCID,iSide).LE.0) CYCLE +#if USE_MPI + ! ignore sides outside of halo region + IF (ElemInfo_Shared(ELEM_HALOFLAG,SideInfo_Shared(SIDE_ELEMID,iSide)).EQ.0) CYCLE +#endif /*USE_MPI*/ + IF (PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,iSide))).EQ.PartBound%SymmetryDim) & + SideIsSymSide(iSide) = .TRUE. +END DO + +SymmetryBCExists = .FALSE. +ALLOCATE(SymmetrySide(1:nElems,1:2)) ! 1: GlobalSide, 2: LocalSide +SymmetrySide = -1 + +! Sanity check: mesh has to be centered at z = 0 +IF(.NOT.ALMOSTEQUALRELATIVE(GEO%zmaxglob,ABS(GEO%zminglob),1e-5)) THEN + SWRITE(*,*) 'Maximum dimension in z:', GEO%zmaxglob + SWRITE(*,*) 'Minimum dimension in z:', GEO%zminglob + SWRITE(*,*) 'Deviation', (ABS(GEO%zmaxglob)-ABS(GEO%zminglob))/ABS(GEO%zminglob), ' > 1e-5' + CALL abort(__STAMP__,'ERROR: Please orient your mesh with one cell in z-direction around 0, |z_min| = z_max !') +END IF + +! Calculation of the correct volume and characteristic length +DO BCSideID=1,nBCSides + locElemID = SideToElem(S2E_ELEM_ID,BCSideID) + iLocSide = SideToElem(S2E_LOC_SIDE_ID,BCSideID) + SideID=GetGlobalNonUniqueSideID(offsetElem+locElemID,iLocSide) + IF (PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).EQ.PartBound%SymmetryDim) THEN + CNElemID = GetCNElemID(SideInfo_Shared(SIDE_ELEMID,SideID)) + iLocSide = SideInfo_Shared(SIDE_LOCALID,SideID) + ! Exclude the symmetry axis (y=0) + IF(Symmetry%Axisymmetric) THEN + IF(MAXVAL(NodeCoords_Shared(2,ElemSideNodeID_Shared(:,iLocSide,CNElemID)+1)).LE.0.0) CYCLE + END IF + ! The z-plane with the positive z component is chosen + IF(MINVAL(NodeCoords_Shared(3,ElemSideNodeID_Shared(:,iLocSide,CNElemID)+1)).GT.(GEO%zmaxglob+GEO%zminglob)/2.) THEN + IF(SymmetrySide(locElemID,1).GT.0) THEN + CALL abort(__STAMP__& + ,'ERROR: PICLas could not determine a unique symmetry surface for 2D/axisymmetric calculation!'//& + ' Please orient your mesh with x as the symmetry axis and positive y as the second/radial direction!') + END IF + SymmetrySide(locElemID,1) = BCSideID + SymmetrySide(locElemID,2) = iLocSide + ! The volume calculated at this point (final volume for the 2D case) corresponds to the cell face area (z-dimension=1) in + ! the xy-plane. + ElemVolume_Shared(CNElemID) = 0.0 + CALL CalcNormAndTangTriangle(area=triarea(1),TriNum=1, SideID=SideID) + CALL CalcNormAndTangTriangle(area=triarea(2),TriNum=2, SideID=SideID) + ElemVolume_Shared(CNElemID) = triarea(1) + triarea(2) + ! Characteristic length is compared to the mean free path as the condition to refine the mesh. For the 2D/axisymmetric case + ! the third dimension is not considered as particle interaction occurs in the xy-plane, effectively reducing the refinement + ! requirement. + ElemCharLength_Shared(CNElemID) = SQRT(ElemVolume_Shared(CNElemID)) + ! Axisymmetric case: The volume is multiplied by the circumference to get the volume of the ring. The cell face in the + ! xy-plane is rotated around the x-axis. The radius is the middle point of the cell face. + IF (Symmetry%Axisymmetric) THEN + radius = 0. + DO iNode = 1, 4 + radius = radius + NodeCoords_Shared(2,ElemSideNodeID_Shared(iNode,iLocSide,CNElemID)+1) + END DO + radius = radius / 4. + ElemVolume_Shared(CNElemID) = ElemVolume_Shared(CNElemID) * 2. * Pi * radius + END IF + SymmetryBCExists = .TRUE. + END IF ! Greater z-coord + END IF +END DO + +IF(.NOT.SymmetryBCExists) CALL abort(__STAMP__,'At least one symmetric BC (in the xy-plane) has to be defined for 2D simulations') + +! LocalVolume & MeshVolume: Recalculate the volume of the mesh of a single process and the total mesh volume +#if USE_MPI +FirstElem = offsetElem - offsetComputeNodeElem + 1 +LastElem = offsetElem - offsetComputeNodeElem + nElems +#else +firstElem = 1 +lastElem = nElems +#endif /*USE_MPI*/ + +LocalVolume = SUM(ElemVolume_Shared(FirstElem:LastElem)) + +#if USE_MPI +CALL BARRIER_AND_SYNC(SideIsSymSide_Shared_Win ,MPI_COMM_SHARED) +CALL BARRIER_AND_SYNC(ElemVolume_Shared_Win ,MPI_COMM_SHARED) +CALL BARRIER_AND_SYNC(ElemCharLength_Shared_Win,MPI_COMM_SHARED) +! Compute-node mesh volume +offsetElemCNProc = offsetElem - offsetComputeNodeElem +CNVolume = SUM(ElemVolume_Shared(offsetElemCNProc+1:offsetElemCNProc+nElems)) +CALL MPI_ALLREDUCE(MPI_IN_PLACE,CNVolume,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_SHARED,iError) +IF (myComputeNodeRank.EQ.0) THEN + ! All-reduce between node leaders + CALL MPI_ALLREDUCE(CNVolume,MeshVolume,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_LEADERS_SHARED,IERROR) +END IF +! Broadcast from node leaders to other processors on the same node +CALL MPI_BCAST(MeshVolume,1, MPI_DOUBLE_PRECISION,0,MPI_COMM_SHARED,iERROR) +#else +MeshVolume = LocalVolume +#endif /*USE_MPI*/ + +END SUBROUTINE InitVolumes_2D + + +SUBROUTINE InitVolumes_1D() +!=================================================================================================================================== +!> Routine determines a symmetry side and calculates the 1D (area faces at x axis) and axisymmetric volumes (cells are +!> revolved around the symmetry axis). The symmetry side will be used later on to determine in which two directions the quadtree +!> shall refine the mesh, skipping the z-dimension to avoid an unnecessary refinement. +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Mesh_Vars ,ONLY: nElems, offsetElem +USE MOD_Particle_Boundary_Vars ,ONLY: PartBound +USE MOD_Particle_Mesh_Vars ,ONLY: GEO,LocalVolume,MeshVolume, SideIsSymSide +USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared,ElemCharLength_Shared +USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared, SideInfo_Shared +USE MOD_Mesh_Tools ,ONLY: GetCNElemID +#if USE_MPI +USE MOD_MPI_Shared +USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED +USE MOD_Particle_Mesh_Vars ,ONLY: offsetComputeNodeElem,SideIsSymSide_Shared,SideIsSymSide_Shared_Win,nNonUniqueGlobalSides +USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared_Win,ElemCharLength_Shared_Win, ElemInfo_Shared +USE MOD_MPI_Shared_Vars ,ONLY: myComputeNodeRank,nComputeNodeProcessors,MPI_COMM_LEADERS_SHARED +#else +USE MOD_Particle_Mesh_Vars ,ONLY: nComputeNodeSides +#endif /*USE_MPI*/ +! IMPLICIT VARIABLE HANDLING + IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: iLocSide, iElem, SideID, iOrder, CNElemID, iSide +REAL :: X(2), MaxCoord, MinCoord +LOGICAL :: SideInPlane, X1Occupied +INTEGER :: firstElem,lastElem, firstSide, lastSide +#if USE_MPI +REAL :: CNVolume +INTEGER :: offsetElemCNProc +#endif /*USE_MPI*/ +!=================================================================================================================================== + +#if USE_MPI +firstElem = offsetElem - offsetComputeNodeElem + 1 +lastElem = offsetElem - offsetComputeNodeElem + nElems +#else +firstElem = 1 +lastElem = nElems +#endif /*USE_MPI*/ + +ALLOCATE(GEO%XMinMax(2,nElems)) + +IF(.NOT.ALMOSTEQUALRELATIVE(GEO%ymaxglob,ABS(GEO%yminglob),1e-5)) THEN + SWRITE(*,*) 'Maximum dimension in y:', GEO%ymaxglob + SWRITE(*,*) 'Minimum dimension in y:', GEO%yminglob + SWRITE(*,*) 'Deviation', (ABS(GEO%ymaxglob)-ABS(GEO%yminglob))/ABS(GEO%yminglob), ' > 1e-5' + CALL abort(__STAMP__& + ,'ERROR: Please orient your mesh with one cell in y-direction around 0, |y_min| = y_max !') +END IF + +IF(.NOT.ALMOSTEQUALRELATIVE(GEO%zmaxglob,ABS(GEO%zminglob),1e-5)) THEN + SWRITE(*,*) 'Maximum dimension in z:', GEO%zmaxglob + SWRITE(*,*) 'Minimum dimension in z:', GEO%zminglob + SWRITE(*,*) 'Deviation', (ABS(GEO%zmaxglob)-ABS(GEO%zminglob))/ABS(GEO%zminglob), ' > 1e-5' + CALL abort(__STAMP__& + ,'ERROR: Please orient your mesh with one cell in z-direction around 0, |z_min| = z_max !') +END IF + +DO iElem = 1,nElems + ! Check if all sides of the element are parallel to xy-, xz-, or yz-plane and Sides parallel to xy-,and xz-plane are symmetric + ! And determine xmin and xmax of the element + X = 0 + X1Occupied = .FALSE. + DO iLocSide = 1,6 + SideInPlane=.FALSE. + SideID = GetGlobalNonUniqueSideID(offsetElem+iElem,iLocSide) + CNElemID = GetCNElemID(SideInfo_Shared(SIDE_ELEMID,SideID)) + DO iOrder=1,3 + MaxCoord = MAXVAL(NodeCoords_Shared(iOrder,ElemSideNodeID_Shared(:,iLocSide,CNElemID)+1)) + MinCoord = MINVAL(NodeCoords_Shared(iOrder,ElemSideNodeID_Shared(:,iLocSide,CNElemID)+1)) + IF(ALMOSTALMOSTEQUAL(MaxCoord,MinCoord).OR.(ALMOSTZERO(MaxCoord).AND.ALMOSTZERO(MinCoord))) THEN + IF(SideInPlane) CALL abort(__STAMP__& + ,'ERROR: Please orient your mesh with all element sides parallel to xy-,xz-,or yz-plane') + SideInPlane=.TRUE. + IF(iOrder.GE.2)THEN + IF(PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,SideID))).NE.PartBound%SymmetryDim) & + CALL abort(__STAMP__,& + 'ERROR: Sides parallel to xy-,and xz-plane has to be the symmetric boundary condition') + END IF + IF(iOrder.EQ.1) THEN + IF(X1Occupied) THEN + X(2) = NodeCoords_Shared(1,ElemSideNodeID_Shared(1,iLocSide,CNElemID)+1) + ELSE + X(1) = NodeCoords_Shared(1,ElemSideNodeID_Shared(1,iLocSide,CNElemID)+1) + X1Occupied = .TRUE. + END IF + END IF !iOrder.EQ.1 + END IF + END DO !iOrder=1,3 + IF(.NOT.SideInPlane) THEN + IPWRITE(*,*) 'ElemID:',iElem,'SideID:',iLocSide + DO iOrder=1,4 + IPWRITE(*,*) 'Node',iOrder,'x:',NodeCoords_Shared(1,ElemSideNodeID_Shared(iOrder,iLocSide,CNElemID)+1), & + 'y:',NodeCoords_Shared(2,ElemSideNodeID_Shared(iOrder,iLocSide,CNElemID)+1), & + 'z:',NodeCoords_Shared(3,ElemSideNodeID_Shared(iOrder,iLocSide,CNElemID)+1) + END DO + CALL abort(__STAMP__& + ,'ERROR: Please orient your mesh with all element sides parallel to xy-,xz-,or yz-plane') + END IF + END DO ! iLocSide = 1,6 + GEO%XMinMax(1,iElem) = MINVAL(X) + GEO%XMinMax(2,iElem) = MAXVAL(X) + ElemVolume_Shared(CNElemID) = ABS(X(1)-X(2)) + ElemCharLength_Shared(CNElemID) = ElemVolume_Shared(CNElemID) +END DO +! LocalVolume & MeshVolume: Recalculate the volume of the mesh of a single process and the total mesh volume +LocalVolume = SUM(ElemVolume_Shared(firstElem:lastElem)) +#if USE_MPI +CALL BARRIER_AND_SYNC(ElemVolume_Shared_Win ,MPI_COMM_SHARED) +CALL BARRIER_AND_SYNC(ElemCharLength_Shared_Win,MPI_COMM_SHARED) +! Compute-node mesh volume +offsetElemCNProc = offsetElem - offsetComputeNodeElem +CNVolume = SUM(ElemVolume_Shared(offsetElemCNProc+1:offsetElemCNProc+nElems)) +CALL MPI_ALLREDUCE(MPI_IN_PLACE,CNVolume,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_SHARED,iError) +IF (myComputeNodeRank.EQ.0) THEN + ! All-reduce between node leaders + CALL MPI_ALLREDUCE(CNVolume,MeshVolume,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_LEADERS_SHARED,IERROR) +END IF +! Broadcast from node leaders to other processors on the same node +CALL MPI_BCAST(MeshVolume,1, MPI_DOUBLE_PRECISION,0,MPI_COMM_SHARED,iERROR) +#else +MeshVolume = SUM(ElemVolume_Shared) +#endif /*USE_MPI*/ + +#if USE_MPI +CALL Allocate_Shared((/nNonUniqueGlobalSides/),SideIsSymSide_Shared_Win,SideIsSymSide_Shared) +CALL MPI_WIN_LOCK_ALL(0,SideIsSymSide_Shared_Win,IERROR) +SideIsSymSide => SideIsSymSide_Shared +! only CN root nullifies +IF(myComputeNodeRank.EQ.0) SideIsSymSide = .FALSE. +! This sync/barrier is required as it cannot be guaranteed that the zeros have been written to memory by the time the MPI_REDUCE +! is executed (see MPI specification). Until the Sync is complete, the status is undefined, i.e., old or new value or utter nonsense. +CALL BARRIER_AND_SYNC(SideIsSymSide_Shared_Win,MPI_COMM_SHARED) +#else +ALLOCATE(SideIsSymSide(nComputeNodeSides)) +SideIsSymSide = .FALSE. +#endif /*USE_MPI*/ + +! Flag of symmetry sides to be skipped during tracking +#if USE_MPI + firstSide = INT(REAL( myComputeNodeRank *nNonUniqueGlobalSides)/REAL(nComputeNodeProcessors))+1 + lastSide = INT(REAL((myComputeNodeRank+1)*nNonUniqueGlobalSides)/REAL(nComputeNodeProcessors)) +#else + firstSide = 1 + lastSide = nComputeNodeSides +#endif + +DO iSide = firstSide, lastSide + SideIsSymSide(iSide) = .FALSE. + ! ignore non-BC sides + IF (SideInfo_Shared(SIDE_BCID,iSide).LE.0) CYCLE +#if USE_MPI + ! ignore sides outside of halo region + IF (ElemInfo_Shared(ELEM_HALOFLAG,SideInfo_Shared(SIDE_ELEMID,iSide)).EQ.0) CYCLE +#endif /*USE_MPI*/ + IF (PartBound%TargetBoundCond(PartBound%MapToPartBC(SideInfo_Shared(SIDE_BCID,iSide))).EQ.PartBound%SymmetryDim) & + SideIsSymSide(iSide) = .TRUE. +END DO + +#if USE_MPI +CALL BARRIER_AND_SYNC(SideIsSymSide_Shared_Win ,MPI_COMM_SHARED) +#endif + +END SUBROUTINE InitVolumes_1D + + +REAL FUNCTION DSMC_2D_CalcSymmetryArea(iLocSide,iElem, ymin, ymax) +!=================================================================================================================================== +!> Calculates the actual area of an element for 2D simulations (plane/axisymmetric) regardless of the mesh dimension in z +!> Utilized in the particle emission (surface flux) and boundary sampling +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Globals_Vars ,ONLY: Pi +USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared, ElemSideNodeID_Shared +USE MOD_Symmetry_Vars ,ONLY: Symmetry +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,INTENT(IN) :: iElem,iLocSide !> iElem is the compute-node element ID +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +REAL, OPTIONAL, INTENT(OUT) :: ymax,ymin +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: iNode +REAL :: P(1:2,1:4), Pmin(2), Pmax(2), Length, MidPoint +!=================================================================================================================================== + +Pmin = HUGE(Pmin) +Pmax = -HUGE(Pmax) + +DO iNode = 1,4 + P(1:2,iNode) = NodeCoords_Shared(1:2,ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1) +END DO + +Pmax(1) = MAXVAL(P(1,:)) +Pmax(2) = MAXVAL(P(2,:)) +Pmin(1) = MINVAL(P(1,:)) +Pmin(2) = MINVAL(P(2,:)) + +IF (PRESENT(ymax).AND.PRESENT(ymin)) THEN + ymin = Pmin(2) + ymax = Pmax(2) +END IF + +Length = SQRT((Pmax(1)-Pmin(1))**2 + (Pmax(2)-Pmin(2))**2) + +MidPoint = (Pmax(2)+Pmin(2)) / 2. +IF(Symmetry%Axisymmetric) THEN + DSMC_2D_CalcSymmetryArea = Length * MidPoint * Pi * 2. + ! Area of the cells on the rotational symmetry axis is set to one + IF(.NOT.(DSMC_2D_CalcSymmetryArea.GT.0.0)) DSMC_2D_CalcSymmetryArea = 1. +ELSE + DSMC_2D_CalcSymmetryArea = Length +END IF +RETURN + +END FUNCTION DSMC_2D_CalcSymmetryArea + + +REAL FUNCTION DSMC_1D_CalcSymmetryArea(iLocSide,iElem) +!=================================================================================================================================== +!> Calculates the actual area of an element for 1D simulations regardless of the mesh dimension in z and y +!> Utilized in the particle emission (surface flux) and boundary sampling +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared, ElemSideNodeID_Shared +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,INTENT(IN) :: iElem,iLocSide +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +REAL :: Pmin, Pmax, Length +!=================================================================================================================================== + +Pmax = MAXVAL(NodeCoords_Shared(1,ElemSideNodeID_Shared(:,iLocSide,iElem)+1)) +Pmin = MINVAL(NodeCoords_Shared(1,ElemSideNodeID_Shared(:,iLocSide,iElem)+1)) + +Length = ABS(Pmax-Pmin) + +! IF(Symmetry%Axisymmetric) THEN +! MidPoint = (Pmax(2)+Pmin(2)) / 2. +! DSMC_1D_CalcSymmetryArea = Length * MidPoint * Pi * 2. +! ! Area of the cells on the rotational symmetry axis is set to one +! IF(.NOT.(DSMC_1D_CalcSymmetryArea.GT.0.0)) DSMC_1D_CalcSymmetryArea = 1. +! ELSE + DSMC_1D_CalcSymmetryArea = Length + IF (DSMC_1D_CalcSymmetryArea.EQ.0.) DSMC_1D_CalcSymmetryArea = 1. +! END IF +RETURN + +END FUNCTION DSMC_1D_CalcSymmetryArea + + +FUNCTION DSMC_2D_CalcSymmetryAreaSubSides(iLocSide,iElem) +!=================================================================================================================================== +!> Calculates the area of the subsides for the insertion with the surface flux +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_Globals_Vars ,ONLY: Pi +USE MOD_DSMC_Vars ,ONLY: RadialWeighting +USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemSideNodeID_Shared +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,INTENT(IN) :: iLocSide,iElem !> iElem is the compute-node element ID +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +REAL :: DSMC_2D_CalcSymmetryAreaSubSides(RadialWeighting%nSubSides) +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: iNode +REAL :: P(1:2,1:4), Pmin(2), Pmax(2), MidPoint, PminTemp, PmaxTemp, Length +!=================================================================================================================================== + +Pmin = HUGE(Pmin) +Pmax = -HUGE(Pmax) + +DO iNode = 1,4 + P(1:2,iNode) = NodeCoords_Shared(1:2,ElemSideNodeID_Shared(iNode,iLocSide,iElem)+1) +END DO + +Pmax(1) = MAXVAL(P(1,:)) +Pmax(2) = MAXVAL(P(2,:)) +Pmin(1) = MINVAL(P(1,:)) +Pmin(2) = MINVAL(P(2,:)) +Length = SQRT((Pmax(1)-Pmin(1))**2 + (Pmax(2)-Pmin(2))**2) + +DO iNode = 1, RadialWeighting%nSubSides + PminTemp = Pmin(2) + (Pmax(2) - Pmin(2))/RadialWeighting%nSubSides*(iNode-1.) + PmaxTemp = Pmin(2) + (Pmax(2) - Pmin(2))/RadialWeighting%nSubSides*iNode + MidPoint = (PmaxTemp+PminTemp) / 2. + DSMC_2D_CalcSymmetryAreaSubSides(iNode) = Length/RadialWeighting%nSubSides * MidPoint * Pi * 2. +END DO + +RETURN + +END FUNCTION DSMC_2D_CalcSymmetryAreaSubSides + + END MODULE MOD_Particle_Mesh_Tools diff --git a/src/particles/particle_mesh/particle_mesh_vars.f90 b/src/particles/particle_mesh/particle_mesh_vars.f90 index e63723f45..bc69ecfec 100644 --- a/src/particles/particle_mesh/particle_mesh_vars.f90 +++ b/src/particles/particle_mesh/particle_mesh_vars.f90 @@ -293,6 +293,7 @@ MODULE MOD_Particle_Mesh_Vars ! periodic case INTEGER,ALLOCATABLE :: PeriodicSFCaseMatrix(:,:) ! matrix to compute periodic cases INTEGER :: NbrOfPeriodicSFCases ! Number of periodic cases +LOGICAL :: AxisymmetricSF ! ==================================================================== INTEGER :: RefMappingGuess ! select guess for mapping into reference ! element @@ -394,6 +395,8 @@ MODULE MOD_Particle_Mesh_Vars REAL,ALLOCATABLE :: ElemTolerance(:) INTEGER, ALLOCATABLE :: ElemToGlobalElemID(:) ! mapping form local-elemid to global-id is built via nodes +INTEGER, ALLOCATABLE :: SymmetrySide(:,:) + !=================================================================================================================================== END MODULE MOD_Particle_Mesh_Vars diff --git a/src/particles/particle_timestep.f90 b/src/particles/particle_timestep.f90 index 062c7fcd2..5185b8246 100644 --- a/src/particles/particle_timestep.f90 +++ b/src/particles/particle_timestep.f90 @@ -108,7 +108,8 @@ SUBROUTINE InitPartTimeStep() !----------------------------------------------------------------------------------------------------------------------------------! USE MOD_Globals USE MOD_ReadInTools ,ONLY: GETLOGICAL, GETINT, GETREAL, GETREALARRAY -USE MOD_Particle_Vars ,ONLY: Symmetry, VarTimeStep +USE MOD_Particle_Vars ,ONLY: VarTimeStep +USE MOD_Symmetry_Vars ,ONLY: Symmetry !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT / OUTPUT VARIABLES @@ -181,7 +182,8 @@ SUBROUTINE VarTimeStep_InitDistribution() USE MOD_io_hdf5 USE MOD_Mesh_Vars ,ONLY: nGlobalElems USE MOD_HDF5_Input, ONLY: OpenDataFile,CloseDataFile,DatasetExists,ReadArray,ReadAttribute,GetDataProps -USE MOD_PARTICLE_Vars, ONLY: VarTimeStep, Symmetry +USE MOD_PARTICLE_Vars, ONLY: VarTimeStep +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Restart_Vars, ONLY: DoRestart, RestartFile, DoMacroscopicRestart, MacroRestartFileName USE MOD_StringTools ,ONLY:STRICMP !----------------------------------------------------------------------------------------------------------------------------------! @@ -410,7 +412,8 @@ REAL FUNCTION GetParticleTimeStep(xPos, yPos, iElem) !=================================================================================================================================== ! MODULES USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: VarTimeStep, Symmetry +USE MOD_Particle_Vars ,ONLY: VarTimeStep +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Particle_Mesh_Vars ,ONLY: GEO ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE diff --git a/src/particles/particle_tools.f90 b/src/particles/particle_tools.f90 index 9674e6c45..5369ae212 100644 --- a/src/particles/particle_tools.f90 +++ b/src/particles/particle_tools.f90 @@ -620,7 +620,7 @@ PPURE REAL FUNCTION CalcRadWeightMPF(yPos, iSpec, iPart) USE MOD_Particle_Vars ,ONLY: Species, PEM USE MOD_Particle_Mesh_Vars ,ONLY: GEO USE MOD_Particle_Mesh_Vars ,ONLY: ElemMidPoint_Shared -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -1318,7 +1318,6 @@ SUBROUTINE MergeCells() END SUBROUTINE MergeCells - SUBROUTINE InitializeParticleMaxwell(iPart,iSpec,iElem,Mode,iInit) !=================================================================================================================================== !> Initialize a particle from a given macroscopic result, requires the macroscopic velocity, translational and internal temperatures @@ -1334,10 +1333,6 @@ SUBROUTINE InitializeParticleMaxwell(iPart,iSpec,iElem,Mode,iInit) !USE MOD_part_tools ,ONLY: CalcRadWeightMPF, CalcEElec_particle, CalcEVib_particle, CalcERot_particle !USE MOD_part_tools ,ONLY: CalcVelocity_maxwell_particle !----------------------------------------------------------------------------------------------------------------------------------- -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES INTEGER, INTENT(IN) :: iPart, iSpec, iElem INTEGER, INTENT(IN) :: Mode !< 1: Macroscopic restart (data for each element) !< 2: Emission distribution (equidistant data from .h5 file) @@ -1555,7 +1550,7 @@ SUBROUTINE CalcPartSymmetryPos(Pos,Velo,ElectronVelo) !=================================================================================================================================== ! MODULES USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- diff --git a/src/particles/particle_vMPF.f90 b/src/particles/particle_vMPF.f90 index bc47e1053..0546e4637 100644 --- a/src/particles/particle_vMPF.f90 +++ b/src/particles/particle_vMPF.f90 @@ -141,7 +141,7 @@ SUBROUTINE MergeParticles(iPartIndx_Node_in, nPart, nPartNew, iElem) USE MOD_Globals ,ONLY: LOG_RAN #ifdef CODE_ANALYZE USE MOD_Globals ,ONLY: unit_stdout,myrank,abort -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #endif /* CODE_ANALYZE */ ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -573,7 +573,7 @@ SUBROUTINE SplitParticles(iPartIndx_Node, nPartIn, nPartNew) USE MOD_Part_Tools ,ONLY: GetNextFreePosition !#ifdef CODE_ANALYZE !USE MOD_Globals ,ONLY: unit_stdout,myrank,abort -!USE MOD_Particle_Vars ,ONLY: Symmetry +!USE MOD_Symmetry_Vars ,ONLY: Symmetry !#endif /* CODE_ANALYZE */ ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE diff --git a/src/particles/particle_vars.f90 b/src/particles/particle_vars.f90 index 59d4faf23..4ae30977d 100644 --- a/src/particles/particle_vars.f90 +++ b/src/particles/particle_vars.f90 @@ -32,12 +32,6 @@ MODULE MOD_Particle_Vars #endif !----------------------------------------------------------------------------------------------------------------------------------- -TYPE tSymmetry - INTEGER :: Order ! 1-3 D - LOGICAL :: Axisymmetric -END TYPE tSymmetry - -TYPE(tSymmetry) :: Symmetry LOGICAL :: DoFieldIonization ! Do Field Ionization by quantum tunneling INTEGER :: FieldIonizationModel !'Field Ionization models. Implemented models are: diff --git a/src/particles/pic/deposition/pic_depo.f90 b/src/particles/pic/deposition/pic_depo.f90 index e15b5ce31..340186f51 100644 --- a/src/particles/pic/deposition/pic_depo.f90 +++ b/src/particles/pic/deposition/pic_depo.f90 @@ -70,7 +70,7 @@ SUBROUTINE DefineParametersPICDeposition() '1D: 2*(N+1)\n'//& '2D: Pi*(N+1)^2\n'//& '3D: (4/3)*Pi*(N+1)^3\n') -CALL prms%CreateLogicalOption('PIC-shapefunction-adaptive-smoothing', 'Enable smooth transition of element-dependent radius when'//& +CALL prms%CreateLogicalOption( 'PIC-shapefunction-adaptive-smoothing', 'Enable smooth transition of element-dependent radius when'//& ' using shape_function_adaptive.', '.FALSE.') END SUBROUTINE DefineParametersPICDeposition @@ -86,8 +86,8 @@ SUBROUTINE InitializeDeposition() USE MOD_Basis ,ONLY: LegendreGaussNodesAndWeights,LegGaussLobNodesAndWeights USE MOD_ChangeBasis ,ONLY: ChangeBasis3D USE MOD_Dielectric_Vars ,ONLY: DoDielectricSurfaceCharge -USE MOD_Interpolation_Vars ,ONLY: xGP,wBary -USE MOD_Mesh_Vars ,ONLY: nElems,sJ +USE MOD_Interpolation_Vars ,ONLY: xGP, wGP, NodeType +USE MOD_Mesh_Vars ,ONLY: nElems, sJ USE MOD_Particle_Vars USE MOD_Particle_Mesh_Vars ,ONLY: nUniqueGlobalNodes, GEO USE MOD_Particle_Mesh_Tools ,ONLY: GetGlobalNonUniqueSideID @@ -97,6 +97,8 @@ SUBROUTINE InitializeDeposition() USE MOD_Preproc USE MOD_ReadInTools ,ONLY: GETREAL,GETINT,GETLOGICAL,GETSTR,GETREALARRAY,GETINTARRAY USE MOD_Mesh_Tools ,ONLY: GetGlobalElemID, GetCNElemID +USE MOD_Interpolation ,ONLY: GetVandermonde +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_Mesh_Vars ,ONLY: offsetElem USE MOD_Particle_Mesh_Vars ,ONLY: NodeToElemInfo,NodeToElemMapping,ElemNodeID_Shared,NodeInfo_Shared @@ -121,8 +123,6 @@ SUBROUTINE InitializeDeposition() ! LOCAL VARIABLES REAL,ALLOCATABLE :: xGP_tmp(:),wGP_tmp(:) INTEGER :: ALLOCSTAT, iElem, i, j, k, kk, ll, mm, iNode -REAL :: DetLocal(1,0:PP_N,0:PP_N,0:PP_N), DetJac(1,0:1,0:1,0:1) -REAL, ALLOCATABLE :: Vdm_tmp(:,:) CHARACTER(255) :: TimeAverageFile #if USE_MPI INTEGER :: UniqueNodeID @@ -210,28 +210,19 @@ SUBROUTINE InitializeDeposition() ALLOCATE(CellVolWeight_Volumes(0:1,0:1,0:1,nElems)) CellVolWeightFac(0:PP_N) = xGP(0:PP_N) CellVolWeightFac(0:PP_N) = (CellVolWeightFac(0:PP_N)+1.0)/2.0 - CALL LegendreGaussNodesAndWeights(1,xGP_tmp,wGP_tmp) - ALLOCATE( Vdm_tmp(0:1,0:PP_N)) - CALL InitializeVandermonde(PP_N,1,wBary,xGP,xGP_tmp,Vdm_tmp) + CellVolWeight_Volumes=0.0 DO iElem=1, nElems - DO k=0,PP_N - DO j=0,PP_N - DO i=0,PP_N - DetLocal(1,i,j,k)=1./sJ(i,j,k,iElem) - END DO ! i=0,PP_N - END DO ! j=0,PP_N - END DO ! k=0,PP_N - CALL ChangeBasis3D(1,PP_N, 1,Vdm_tmp, DetLocal(:,:,:,:),DetJac(:,:,:,:)) - DO k=0,1 - DO j=0,1 - DO i=0,1 - CellVolWeight_Volumes(i,j,k,iElem) = DetJac(1,i,j,k)*wGP_tmp(i)*wGP_tmp(j)*wGP_tmp(k) - END DO ! i=0,PP_N - END DO ! j=0,PP_N - END DO ! k=0,PP_N + DO i=0,PP_N;DO j=0,PP_N;DO k=0,PP_N + CellVolWeight_Volumes(0,0,0,iElem) = CellVolWeight_Volumes(0,0,0,iElem) + 1/sJ(i,j,k,iElem)*((1.-xGP(i))*(1.-xGP(j))*(1.-xGP(k))*wGP(i)*wGP(j)*wGP(k)/8.) + CellVolWeight_Volumes(0,0,1,iElem) = CellVolWeight_Volumes(0,0,1,iElem) + 1/sJ(i,j,k,iElem)*((1.-xGP(i))*(1.-xGP(j))*(1.+xGP(k))*wGP(i)*wGP(j)*wGP(k)/8.) + CellVolWeight_Volumes(0,1,0,iElem) = CellVolWeight_Volumes(0,1,0,iElem) + 1/sJ(i,j,k,iElem)*((1.-xGP(i))*(1.+xGP(j))*(1.-xGP(k))*wGP(i)*wGP(j)*wGP(k)/8.) + CellVolWeight_Volumes(0,1,1,iElem) = CellVolWeight_Volumes(0,1,1,iElem) + 1/sJ(i,j,k,iElem)*((1.-xGP(i))*(1.+xGP(j))*(1.+xGP(k))*wGP(i)*wGP(j)*wGP(k)/8.) + CellVolWeight_Volumes(1,0,0,iElem) = CellVolWeight_Volumes(1,0,0,iElem) + 1/sJ(i,j,k,iElem)*((1.+xGP(i))*(1.-xGP(j))*(1.-xGP(k))*wGP(i)*wGP(j)*wGP(k)/8.) + CellVolWeight_Volumes(1,0,1,iElem) = CellVolWeight_Volumes(1,0,1,iElem) + 1/sJ(i,j,k,iElem)*((1.+xGP(i))*(1.-xGP(j))*(1.+xGP(k))*wGP(i)*wGP(j)*wGP(k)/8.) + CellVolWeight_Volumes(1,1,0,iElem) = CellVolWeight_Volumes(1,1,0,iElem) + 1/sJ(i,j,k,iElem)*((1.+xGP(i))*(1.+xGP(j))*(1.-xGP(k))*wGP(i)*wGP(j)*wGP(k)/8.) + CellVolWeight_Volumes(1,1,1,iElem) = CellVolWeight_Volumes(1,1,1,iElem) + 1/sJ(i,j,k,iElem)*((1.+xGP(i))*(1.+xGP(j))*(1.+xGP(k))*wGP(i)*wGP(j)*wGP(k)/8.) + END DO; END DO; END DO END DO - DEALLOCATE(Vdm_tmp) - DEALLOCATE(wGP_tmp, xGP_tmp) CASE('cell_volweight_mean') #if USE_MPI ALLOCATE(DoNodeMapping(0:nProcessors_Global-1),SendNode(1:nUniqueGlobalNodes)) @@ -503,9 +494,12 @@ SUBROUTINE InitializeDeposition() SELECT CASE(TRIM(DepositionType)) CASE('shape_function_cc', 'shape_function_adaptive') w_sf = 1.0 ! set dummy value + CASE('shape_function') + IF(Symmetry%axisymmetric) CALL abort(__STAMP__,'Axisymmetric simulations only with shape_function_cc or shape_function_adaptive!') END SELECT ! --- Set periodic case matrix for shape function deposition (virtual displacement of particles in the periodic directions) + CALL InitAxisymmetrySF() CALL InitPeriodicSFCaseMatrix() ! --- Set element flag for cycling already completed elements @@ -526,7 +520,7 @@ SUBROUTINE InitializeDeposition() #else ALLOCATE(ChargeSFDone(1:nElems)) #endif /*USE_MPI*/ - +CASE('cell_mean') CASE DEFAULT CALL abort(__STAMP__,'Unknown DepositionType in pic_depo.f90') END SELECT @@ -764,7 +758,7 @@ SUBROUTINE InitPeriodicSFCaseMatrix() IF (GEO%nPeriodicVectors.LE.0) THEN ! Set defaults and return in non-periodic case - NbrOfPeriodicSFCases = 1 + NbrOfPeriodicSFCases = 0 ALLOCATE(PeriodicSFCaseMatrix(1:1,1:3)) PeriodicSFCaseMatrix(:,:) = 0 @@ -825,6 +819,38 @@ SUBROUTINE InitPeriodicSFCaseMatrix() END SUBROUTINE InitPeriodicSFCaseMatrix +!=================================================================================================================================== +!> Fill PeriodicSFCaseMatrix when using shape function deposition in combination with periodic boundaries +!=================================================================================================================================== +SUBROUTINE InitAxisymmetrySF() +! MODULES +USE MOD_Particle_Boundary_Vars ,ONLY: PartBound,nPartBound +USE MOD_Symmetry_Vars ,ONLY: Symmetry +USE MOD_Particle_Mesh_Vars ,ONLY: AxisymmetricSF +USE MOD_ReadInTools ,ONLY: PrintOption +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------! +! INPUT / OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: iPartBound +!=================================================================================================================================== +AxisymmetricSF = .FALSE. +IF (Symmetry%Axisymmetric) THEN + ! Check boundaries for symmetry axis + DO iPartBound=1,nPartBound + IF(PartBound%TargetBoundCond(iPartBound).EQ.PartBound%SymmetryAxis) THEN + AxisymmetricSF = .TRUE. + CALL PrintOption('Found symmetry axis for shape function deposition','INFO',LogOpt=AxisymmetricSF) + RETURN + END IF + END DO +END IF + +END SUBROUTINE InitAxisymmetrySF + + SUBROUTINE Deposition(doParticle_In, stage_opt) !============================================================================================================================ ! This subroutine performs the deposition of the particle charge and current density to the grid diff --git a/src/particles/pic/deposition/pic_depo_method.f90 b/src/particles/pic/deposition/pic_depo_method.f90 index 5052a9f3f..cef3e5d47 100644 --- a/src/particles/pic/deposition/pic_depo_method.f90 +++ b/src/particles/pic/deposition/pic_depo_method.f90 @@ -42,6 +42,7 @@ SUBROUTINE DepositionMethodInterface(doParticle_In,stage_opt) INTEGER,PARAMETER :: PRM_DEPO_SF_ADAPTIVE= 2 ! shape_function_adaptive INTEGER,PARAMETER :: PRM_DEPO_CVW = 6 ! cell_volweight INTEGER,PARAMETER :: PRM_DEPO_CVWM = 12 ! cell_volweight_mean +INTEGER,PARAMETER :: PRM_DEPO_CM = 27 ! cell_mean INTERFACE InitDepositionMethod MODULE PROCEDURE InitDepositionMethod @@ -79,7 +80,9 @@ SUBROUTINE DefineParametersDepositionMethod() '1.2) shape_function_cc ('//TRIM(int2strf(PRM_DEPO_SF_CC))//')\n' //& '1.3) shape_function_adaptive ('//TRIM(int2strf(PRM_DEPO_SF_ADAPTIVE))//')\n' //& '2.) cell_volweight ('//TRIM(int2strf(PRM_DEPO_CVW))//')\n' //& - '3.) cell_volweight_mean ('//TRIM(int2strf(PRM_DEPO_CVWM))//')' & + '3.) cell_volweight_mean ('//TRIM(int2strf(PRM_DEPO_CVWM))//')\n' //& + '4.) cell_mean ('//TRIM(int2strf(PRM_DEPO_CM))//')' & + ,'cell_volweight') CALL addStrListEntry('PIC-Deposition-Type' , 'shape_function' , PRM_DEPO_SF) @@ -87,6 +90,7 @@ SUBROUTINE DefineParametersDepositionMethod() CALL addStrListEntry('PIC-Deposition-Type' , 'shape_function_adaptive' , PRM_DEPO_SF_ADAPTIVE) CALL addStrListEntry('PIC-Deposition-Type' , 'cell_volweight' , PRM_DEPO_CVW) CALL addStrListEntry('PIC-Deposition-Type' , 'cell_volweight_mean' , PRM_DEPO_CVWM) +CALL addStrListEntry('PIC-Deposition-Type' , 'cell_mean' , PRM_DEPO_CM) END SUBROUTINE DefineParametersDepositionMethod @@ -154,6 +158,9 @@ SUBROUTINE InitDepositionMethod() Case(PRM_DEPO_CVWM) ! cell_volweight_mean DepositionType = 'cell_volweight_mean' DepositionMethod => DepositionMethod_CVWM +Case(PRM_DEPO_CM) ! cell_mean + DepositionType = 'cell_mean' + DepositionMethod => DepositionMethod_CM CASE DEFAULT CALL CollectiveStop(__STAMP__,'Unknown DepositionMethod!' ,IntInfo=DepositionType_loc) END SELECT @@ -208,6 +215,7 @@ SUBROUTINE InitDepositionMethod() CALL DepositionMethod_SF() CALL DepositionMethod_CVW() CALL DepositionMethod_CVWM() +CALL DepositionMethod_CM() END SUBROUTINE InitDepositionMethod @@ -742,6 +750,113 @@ SUBROUTINE DepositionMethod_CVWM(doParticle_In, stage_opt) END SUBROUTINE DepositionMethod_CVWM +SUBROUTINE DepositionMethod_CM(doParticle_In, stage_opt) +!=================================================================================================================================== +! 'cell_volweight' +! Linear charge density distribution within a cell (discontinuous across cell interfaces) +!=================================================================================================================================== +! MODULES +USE MOD_Preproc +USE MOD_Particle_Vars ,ONLY: Species, PartSpecies,PDM,PEM,usevMPF,PartMPF +USE MOD_Particle_Vars ,ONLY: PartState +USE MOD_PICDepo_Vars ,ONLY: PartSource +USE MOD_Part_Tools ,ONLY: isDepositParticle +USE MOD_Mesh_Tools ,ONLY: GetCNElemID +USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared +USE MOD_Mesh_Vars ,ONLY: offSetElem +#if USE_LOADBALANCE +USE MOD_LoadBalance_Timers ,ONLY: LBStartTime,LBPauseTime,LBElemSplitTime,LBElemPauseTime_avg +USE MOD_LoadBalance_Timers ,ONLY: LBElemSplitTime_avg +#endif /*USE_LOADBALANCE*/ +USE MOD_Mesh_Vars ,ONLY: nElems +#if ((USE_HDG) && (PP_nVar==1)) +USE MOD_TimeDisc_Vars ,ONLY: dt,dt_Min +#endif +#if USE_MPI +USE MOD_MPI_Shared ,ONLY: BARRIER_AND_SYNC +#endif /*USE_MPI*/ +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +LOGICAL,INTENT(IN),OPTIONAL :: doParticle_In(1:PDM%ParticleVecLength) ! TODO: definition of this variable +INTEGER,INTENT(IN),OPTIONAL :: stage_opt +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +REAL :: Charge, TSource(1:4) +#if USE_LOADBALANCE +REAL :: tLBStart +#endif /*USE_LOADBALANCE*/ +#if !((USE_HDG) && (PP_nVar==1)) +INTEGER, PARAMETER :: SourceDim=1 +LOGICAL, PARAMETER :: doCalculateCurrentDensity=.TRUE. +#else +LOGICAL :: doCalculateCurrentDensity +INTEGER :: SourceDim +#endif +INTEGER :: iPart,iElem, iDim +!=================================================================================================================================== + +#if USE_LOADBALANCE + CALL LBStartTime(tLBStart) ! Start time measurement +#endif /*USE_LOADBALANCE*/ + +! Check whether charge and current density have to be computed or just the charge density +#if ((USE_HDG) && (PP_nVar==1)) +IF(ALMOSTEQUAL(dt,dt_Min(DT_ANALYZE)).OR.ALMOSTEQUAL(dt,dt_Min(DT_END)))THEN + doCalculateCurrentDensity=.TRUE. + SourceDim=1 +ELSE ! do not calculate current density + doCalculateCurrentDensity=.FALSE. + SourceDim=4 +END IF +#endif + +DO iPart = 1,PDM%ParticleVecLength + ! TODO: The cycle with .AND.doParticle_In is used for analysis or IMPA + IF(PRESENT(doParticle_In))THEN + IF (.NOT.(PDM%ParticleInside(iPart).AND.doParticle_In(iPart))) CYCLE + ELSE + IF (.NOT.PDM%ParticleInside(iPart)) CYCLE + END IF + ! Don't deposit neutral particles! + IF(.NOT.isDepositParticle(iPart)) CYCLE + IF (usevMPF) THEN + Charge= Species(PartSpecies(iPart))%ChargeIC * PartMPF(iPart) + ELSE + Charge= Species(PartSpecies(iPart))%ChargeIC * Species(PartSpecies(iPart))%MacroParticleFactor + END IF ! usevMPF + IF(doCalculateCurrentDensity)THEN + TSource(1:3) = PartState(4:6,iPart)*Charge + ELSE + TSource(1:3) = 0.0 + END IF + iElem = PEM%LocalElemID(iPart) + TSource(4) = Charge + DO iDim=SourceDim,4 + PartSource(iDim,:,:,:,iElem) = PartSource(iDim,:,:,:,iElem) + TSource(iDim) + END DO +#if USE_LOADBALANCE + CALL LBElemSplitTime(iElem,tLBStart) ! Split time measurement (Pause/Stop and Start again) and add time to iElem +#endif /*USE_LOADBALANCE*/ +END DO + +DO iElem=1, nElems + PartSource(SourceDim:4,:,:,:,iElem) = PartSource(SourceDim:4,:,:,:,iElem) / ElemVolume_Shared(GetCNElemID(iElem+offSetElem)) +END DO + +#if USE_LOADBALANCE +CALL LBElemSplitTime_avg(tLBStart) ! Average over the number of elems (and Start again) +#endif /*USE_LOADBALANCE*/ + +! Suppress compiler warnings +RETURN +iPart=stage_opt +END SUBROUTINE DepositionMethod_CM + + SUBROUTINE DepositionMethod_SF(doParticle_In, stage_opt) !=================================================================================================================================== ! 'shape_function' @@ -877,7 +992,7 @@ END SUBROUTINE DepositionMethod_SF !=================================================================================================================================== !> Set NodeSource to zero on Dirichlet sides to account for mirror charges !=================================================================================================================================== -SUBROUTINE NullifyNodeSourceDirichletSides(SourceDim ) +SUBROUTINE NullifyNodeSourceDirichletSides(SourceDim) ! MODULES USE MOD_HDG_Vars ,ONLY: nDirichletBCSides,DirichletBC USE MOD_Particle_Mesh_Vars ,ONLY: ElemSideNodeID_Shared,NodeInfo_Shared @@ -918,5 +1033,4 @@ SUBROUTINE NullifyNodeSourceDirichletSides(SourceDim ) END SUBROUTINE NullifyNodeSourceDirichletSides #endif /*USE_HDG*/ - END MODULE MOD_PICDepo_Method diff --git a/src/particles/pic/deposition/pic_depo_shapefunction_tools.f90 b/src/particles/pic/deposition/pic_depo_shapefunction_tools.f90 index a0edfc31b..1dabc68b7 100644 --- a/src/particles/pic/deposition/pic_depo_shapefunction_tools.f90 +++ b/src/particles/pic/deposition/pic_depo_shapefunction_tools.f90 @@ -42,7 +42,7 @@ SUBROUTINE calcSfSource(SourceSize_in,ChargeMPF,PartPos,PartID,PartVelo) USE MOD_Globals USE MOD_PICDepo_Vars ,ONLY: DepositionType,dimFactorSF,sfDepo3D,r_sf,r2_sf,r2_sf_inv,SFElemr2_Shared,w_sf USE MOD_PICDepo_Vars ,ONLY: totalChargePeriodicSF,SFAdaptiveSmoothing -USE MOD_Particle_Mesh_Vars ,ONLY: NbrOfPeriodicSFCases +USE MOD_Particle_Mesh_Vars ,ONLY: NbrOfPeriodicSFCases,AxisymmetricSF USE MOD_Particle_Vars ,ONLY: PEM USE MOD_Mesh_Tools ,ONLY: GetCNElemID !----------------------------------------------------------------------------------------------------------------------------------- @@ -86,9 +86,9 @@ SUBROUTINE calcSfSource(SourceSize_in,ChargeMPF,PartPos,PartID,PartVelo) END IF ! Check if periodic sides are present, otherwise simply use PartPos for deposition -IF(NbrOfPeriodicSFCases.GT.1)THEN +IF((NbrOfPeriodicSFCases.GT.1) .OR. AxisymmetricSF) THEN - IF(TRIM(DepositionType).EQ.'shape_function_adaptive')THEN + IF(TRIM(DepositionType).EQ.'shape_function_adaptive') THEN ! Get radius/radius squared/inverse radius squared for each CN element when using adaptive SF OrigElem = PEM%GlobalElemID(PartID) OrigCNElemID = GetCNElemID(OrigElem) @@ -101,7 +101,7 @@ SUBROUTINE calcSfSource(SourceSize_in,ChargeMPF,PartPos,PartID,PartVelo) r2_sf_inv_tmp= r2_sf_inv END IF ! TRIM(DepositionType).EQ.'shape_function_adaptive' - IF(TRIM(DepositionType).EQ.'shape_function')THEN + IF(TRIM(DepositionType).EQ.'shape_function') THEN ! Incorporate the coefficient that considers the volume integral of the kernel Fac = Fac*w_sf ELSE @@ -113,7 +113,13 @@ SUBROUTINE calcSfSource(SourceSize_in,ChargeMPF,PartPos,PartID,PartVelo) PartPosShifted(1:3) = GetPartPosShifted(iCase,PartPos(1:3)) CALL calcTotalChargePeriodic_cc(PartPosShifted , Fac(4) , totalChargePeriodicSF , r_sf_tmp , r2_sf_tmp , r2_sf_inv_tmp) END DO - + IF(AxisymmetricSF) THEN + IF(NbrOfPeriodicSFCases.EQ.0) CALL calcTotalChargePeriodic_cc(PartPos , Fac(4) , totalChargePeriodicSF , r_sf_tmp , r2_sf_tmp , r2_sf_inv_tmp) + PartPosShifted(1) = PartPos(1) + PartPosShifted(2) = -PartPos(2) + PartPosShifted(3) = PartPos(3) + CALL calcTotalChargePeriodic_cc(PartPosShifted , Fac(4) , totalChargePeriodicSF , r_sf_tmp , r2_sf_tmp , r2_sf_inv_tmp) + END IF ! Check if the charge is to be distributed over a line (1D) or area (2D) IF(.NOT.sfDepo3D)THEN totalChargePeriodicSF = totalChargePeriodicSF / dimFactorSF @@ -128,6 +134,13 @@ SUBROUTINE calcSfSource(SourceSize_in,ChargeMPF,PartPos,PartID,PartVelo) PartPosShifted(1:3) = GetPartPosShifted(iCase,PartPos(1:3)) CALL depoChargeOnDOFsSF(PartPosShifted , SourceSize , Fac , r_sf_tmp , r2_sf_tmp , r2_sf_inv_tmp) END DO + IF(AxisymmetricSF) THEN + IF(NbrOfPeriodicSFCases.EQ.0) CALL depoChargeOnDOFsSF(PartPos , SourceSize , Fac , r_sf_tmp , r2_sf_tmp , r2_sf_inv_tmp) + PartPosShifted(1) = PartPos(1) + PartPosShifted(2) = -PartPos(2) + PartPosShifted(3) = PartPos(3) + CALL depoChargeOnDOFsSF(PartPosShifted , SourceSize , Fac , r_sf_tmp , r2_sf_tmp , r2_sf_inv_tmp) + END IF ELSE ! Select deposition type @@ -179,6 +192,7 @@ SUBROUTINE depoChargeOnDOFsSF_RGetAccumulate(Position,SourceSize,Fac) #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: nDeposPerElem #endif /*USE_LOADBALANCE*/ +USE MOD_Symmetry_Vars ,ONLY: Symmetry !----------------------------------------------------------------------------------------------------------------------------------- IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -824,7 +838,7 @@ SUBROUTINE UpdatePartSource(dim1,k,l,m,globElemID,Source) !#endif #if USE_MPI ELSE -! IF (myComputeNodeRank.EQ.0) THEN +! IF (myComputeNodeRank.EQ.0) THEN ! PartSourceGlob(dim1:4,k,l,m,CNElemID) = PartSourceGlob(dim1:4,k,l,m,CNElemID) + Source(dim1:4) ! ELSE ASSOCIATE( ShapeID => SendElemShapeID(CNElemID)) @@ -835,7 +849,7 @@ SUBROUTINE UpdatePartSource(dim1,k,l,m,globElemID,Source) END IF ! IF (CNElemID.GT.nComputeNodeElems) THEN ! ExRankID = CNRankToSendRank(0) -! ELSE +! ELSE ! ExRankID = CNRankToSendRank(ElemInfo_Shared(ELEM_RANK,globElemID)-ComputeNodeRootRank) ExRankID = CNRankToSendRank(ElemInfo_Shared(ELEM_RANK,globElemID)) ! END IF @@ -1023,6 +1037,7 @@ SUBROUTINE InitShapeFunctionDimensionalty() USE MOD_Globals_Vars ,ONLY: PI USE MOD_Mesh_Vars ,ONLY: nGlobalElems USE MOD_PICDepo_Tools ,ONLY: beta +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_LOADBALANCE*/ @@ -1059,9 +1074,17 @@ SUBROUTINE InitShapeFunctionDimensionalty() CASE (1) ! --- 1D shape function ------------------------------------------------------------------------------------------------- ! Set perpendicular directions IF(dim_sf_dir.EQ.1)THEN ! Shape function deposits charge in x-direction - dimFactorSF = (GEO%ymaxglob-GEO%yminglob)*(GEO%zmaxglob-GEO%zminglob) + IF(Symmetry%axisymmetric) THEN + dimFactorSF = (GEO%ymaxglob-GEO%yminglob)*2*PI !(GEO%zmaxglob-GEO%zminglob) ! AXISYMMETRIC HDG + ELSE + dimFactorSF = (GEO%ymaxglob-GEO%yminglob)*(GEO%zmaxglob-GEO%zminglob) + END IF ELSE IF (dim_sf_dir.EQ.2)THEN ! Shape function deposits charge in y-direction - dimFactorSF = (GEO%xmaxglob-GEO%xminglob)*(GEO%zmaxglob-GEO%zminglob) + IF(Symmetry%axisymmetric) THEN + dimFactorSF = (GEO%xmaxglob-GEO%xminglob)*2*PI !(GEO%zmaxglob-GEO%zminglob) ! AXISYMMETRIC HDG + ELSE + dimFactorSF = (GEO%xmaxglob-GEO%xminglob)*(GEO%zmaxglob-GEO%zminglob) ! AXISYMMETRIC HDG + END IF ELSE IF (dim_sf_dir.EQ.3)THEN ! Shape function deposits charge in z-direction dimFactorSF = (GEO%xmaxglob-GEO%xminglob)*(GEO%ymaxglob-GEO%yminglob) END IF @@ -1086,7 +1109,11 @@ SUBROUTINE InitShapeFunctionDimensionalty() ELSE IF (dim_sf_dir.EQ.2)THEN ! Shape function deposits charge in x-z-direction (const. in y) dimFactorSF = (GEO%ymaxglob-GEO%yminglob) ELSE IF (dim_sf_dir.EQ.3)THEN! Shape function deposits charge in x-y-direction (const. in z) - dimFactorSF = (GEO%zmaxglob-GEO%zminglob) + IF(Symmetry%axisymmetric) THEN + dimFactorSF = 1 + ELSE + dimFactorSF = (GEO%zmaxglob-GEO%zminglob) ! AXISYMMETRIC HDG + END IF END IF ! Set prefix factor @@ -1105,6 +1132,9 @@ SUBROUTINE InitShapeFunctionDimensionalty() ! Set shape function length (3D volume) VolumeShapeFunction=PI*(r_sf_loc**2)*dimFactorSF + IF(Symmetry%axisymmetric) THEN + VolumeShapeFunction=VolumeShapeFunction*2*PI*GEO%ymaxglob + END IF ! Calculate number of 2D DOF (assume third direction with 1 cell layer and width dimFactorSF) nTotalDOF=nGlobalElems*(PP_N+1)**2 diff --git a/src/particles/pic/deposition/pic_depo_tools.f90 b/src/particles/pic/deposition/pic_depo_tools.f90 index 72051ca13..d01441eb5 100644 --- a/src/particles/pic/deposition/pic_depo_tools.f90 +++ b/src/particles/pic/deposition/pic_depo_tools.f90 @@ -197,9 +197,7 @@ SUBROUTINE CalcCellLocNodeVolumes() ! MODULES USE MOD_Globals USE MOD_PreProc -USE MOD_Basis ,ONLY: InitializeVandermonde -USE MOD_ChangeBasis ,ONLY: ChangeBasis3D -USE MOD_Interpolation_Vars ,ONLY: wGP, xGP, wBary +USE MOD_Interpolation_Vars ,ONLY: wGP, xGP #if USE_MPI USE MOD_MPI_Shared USE MOD_MPI_Shared_Vars ,ONLY: nComputeNodeTotalElems, nComputeNodeProcessors, myComputeNodeRank, MPI_COMM_SHARED @@ -217,8 +215,6 @@ SUBROUTINE CalcCellLocNodeVolumes() ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -REAL :: Vdm_loc(0:1,0:PP_N),wGP_loc,xGP_loc(0:1),DetJac(1,0:1,0:1,0:1) -REAL :: DetLocal(1,0:PP_N,0:PP_N,0:PP_N) INTEGER :: j,k,l,iElem, firstElem, lastElem, iNode, jNode REAL :: NodeVolumeLoc(1:nUniqueGlobalNodes) #if USE_MPI @@ -248,38 +244,23 @@ SUBROUTINE CalcCellLocNodeVolumes() lastElem = nElems #endif /*USE_MPI*/ -IF (PP_N.NE.1) THEN - xGP_loc(0) = -0.5 - xGP_loc(1) = 0.5 - wGP_loc = 1. - CALL InitializeVandermonde(PP_N,1,wBary,xGP,xGP_loc, Vdm_loc) -END IF ! ElemNodeID and ElemsJ use compute node elems DO iElem = firstElem, lastElem - IF (PP_N.EQ.1) THEN - wGP_loc = wGP(0) - DO j=0, PP_N; DO k=0, PP_N; DO l=0, PP_N - DetJac(1,j,k,l)=1./ElemsJ(j,k,l,iElem) - END DO; END DO; END DO - ELSE - DO j=0, PP_N; DO k=0, PP_N; DO l=0, PP_N - DetLocal(1,j,k,l)=1./ElemsJ(j,k,l,iElem) - END DO; END DO; END DO - CALL ChangeBasis3D(1,PP_N, 1, Vdm_loc, DetLocal(:,:,:,:),DetJac(:,:,:,:)) - END IF #if USE_MPI ASSOCIATE( NodeVolume => NodeVolumeLoc ) #endif /*USE_MPI*/ ! Get UniqueNodeIDs NodeID = NodeInfo_Shared(ElemNodeID_Shared(1:8,iElem)) - NodeVolume(NodeID(1)) = NodeVolume(NodeID(1)) + DetJac(1,0,0,0) - NodeVolume(NodeID(2)) = NodeVolume(NodeID(2)) + DetJac(1,1,0,0) - NodeVolume(NodeID(3)) = NodeVolume(NodeID(3)) + DetJac(1,1,1,0) - NodeVolume(NodeID(4)) = NodeVolume(NodeID(4)) + DetJac(1,0,1,0) - NodeVolume(NodeID(5)) = NodeVolume(NodeID(5)) + DetJac(1,0,0,1) - NodeVolume(NodeID(6)) = NodeVolume(NodeID(6)) + DetJac(1,1,0,1) - NodeVolume(NodeID(7)) = NodeVolume(NodeID(7)) + DetJac(1,1,1,1) - NodeVolume(NodeID(8)) = NodeVolume(NodeID(8)) + DetJac(1,0,1,1) + DO j=0,PP_N;DO k=0,PP_N;DO l=0,PP_N + NodeVolume(NodeID(1)) = NodeVolume(NodeID(1)) + 1/ElemsJ(j,k,l,iElem)*((1.-xGP(j))*(1.-xGP(k))*(1.-xGP(l))*wGP(j)*wGP(k)*wGP(l)/8.) + NodeVolume(NodeID(2)) = NodeVolume(NodeID(2)) + 1/ElemsJ(j,k,l,iElem)*((1.+xGP(j))*(1.-xGP(k))*(1.-xGP(l))*wGP(j)*wGP(k)*wGP(l)/8.) + NodeVolume(NodeID(3)) = NodeVolume(NodeID(3)) + 1/ElemsJ(j,k,l,iElem)*((1.+xGP(j))*(1.+xGP(k))*(1.-xGP(l))*wGP(j)*wGP(k)*wGP(l)/8.) + NodeVolume(NodeID(4)) = NodeVolume(NodeID(4)) + 1/ElemsJ(j,k,l,iElem)*((1.-xGP(j))*(1.+xGP(k))*(1.-xGP(l))*wGP(j)*wGP(k)*wGP(l)/8.) + NodeVolume(NodeID(5)) = NodeVolume(NodeID(5)) + 1/ElemsJ(j,k,l,iElem)*((1.-xGP(j))*(1.-xGP(k))*(1.+xGP(l))*wGP(j)*wGP(k)*wGP(l)/8.) + NodeVolume(NodeID(6)) = NodeVolume(NodeID(6)) + 1/ElemsJ(j,k,l,iElem)*((1.+xGP(j))*(1.-xGP(k))*(1.+xGP(l))*wGP(j)*wGP(k)*wGP(l)/8.) + NodeVolume(NodeID(7)) = NodeVolume(NodeID(7)) + 1/ElemsJ(j,k,l,iElem)*((1.+xGP(j))*(1.+xGP(k))*(1.+xGP(l))*wGP(j)*wGP(k)*wGP(l)/8.) + NodeVolume(NodeID(8)) = NodeVolume(NodeID(8)) + 1/ElemsJ(j,k,l,iElem)*((1.-xGP(j))*(1.+xGP(k))*(1.+xGP(l))*wGP(j)*wGP(k)*wGP(l)/8.) + END DO; END DO; END DO #if USE_MPI END ASSOCIATE #endif /*USE_MPI*/ diff --git a/src/particles/pic/interpolation/pic_interpolation_tools.f90 b/src/particles/pic/interpolation/pic_interpolation_tools.f90 index eb4625175..877ef4db9 100644 --- a/src/particles/pic/interpolation/pic_interpolation_tools.f90 +++ b/src/particles/pic/interpolation/pic_interpolation_tools.f90 @@ -498,6 +498,7 @@ PPURE FUNCTION InterpolateVariableExternalField2D(Pos) USE MOD_PICInterpolation_Vars ,ONLY: DeltaExternalField,VariableExternalField USE MOD_PICInterpolation_Vars ,ONLY: VariableExternalFieldN USE MOD_PICInterpolation_Vars ,ONLY: VariableExternalFieldMin,VariableExternalFieldMax +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -509,76 +510,88 @@ PPURE FUNCTION InterpolateVariableExternalField2D(Pos) !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: iPos,jPos !< index in array (equidistant subdivision assumed) -REAL :: r,delta,f(1:2),mat(2,2),dx(1:2),dy(1:2),vec(1:2) +REAL :: r,delta,f(1:2),mat(2,2),dx(1:2),dy(1:2),vec(1:2),x,y,z INTEGER :: idx1,idx2,idx3,idx4,i !=================================================================================================================================== -ASSOCIATE(& - x => Pos(1) ,& - y => Pos(2) ,& - z => Pos(3) & - ) + +IF(Symmetry%Axisymmetric) THEN + ! 2D axisymmetry requires x as the rotational axis and a positive y + ! Particle position has been rotated onto the symmetry plane -> z = 0 and y-coordinate is equal to the radius + y = Pos(2) + r = Pos(2) + z = Pos(1) +ELSE + x = Pos(1) + y = Pos(2) + z = Pos(3) r = SQRT(x*x + y*y) +END IF - IF(r.GT.VariableExternalFieldMax(1))THEN - InterpolateVariableExternalField2D = 0. - ELSEIF(r.LT.VariableExternalFieldMin(1))THEN - InterpolateVariableExternalField2D = 0. - ELSEIF(z.GT.VariableExternalFieldMax(2))THEN - InterpolateVariableExternalField2D = 0. - ELSEIF(z.LT.VariableExternalFieldMin(2))THEN - InterpolateVariableExternalField2D = 0. - ELSE +IF(r.GT.VariableExternalFieldMax(1))THEN + InterpolateVariableExternalField2D = 0. +ELSEIF(r.LT.VariableExternalFieldMin(1))THEN + InterpolateVariableExternalField2D = 0. +ELSEIF(z.GT.VariableExternalFieldMax(2))THEN + InterpolateVariableExternalField2D = 0. +ELSEIF(z.LT.VariableExternalFieldMin(2))THEN + InterpolateVariableExternalField2D = 0. +ELSE - ! Get index in r and z - iPos = INT((r-VariableExternalField(1,1))/DeltaExternalField(1)) + 1 ! dr = DeltaExternalField(1) - jPos = INT((z-VariableExternalField(2,1))/DeltaExternalField(2)) + 1 ! dz = DeltaExternalField(2) - - ! Catch problem when r or z are exactly at the upper boundary and INT() does not round to the lower integer (do not add +1 in - ! this case) - iPos = MIN(iPos, VariableExternalFieldN(2) - 1 ) - jPos = MIN(jPos, VariableExternalFieldN(1) - 1 ) - - - ! Shift all points by Nz = EmissionDistributionNum(1) - ASSOCIATE( Nz => VariableExternalFieldN(1) ) - ! 1.1 - idx1 = (iPos-1)*Nz + jPos - ! 2.1 - idx2 = (iPos-1)*Nz + jPos + 1 - ! 1.2 - idx3 = iPos*Nz + jPos - ! 2.2 - idx4 = iPos*Nz + jPos + 1 - END ASSOCIATE + ! Get index in r and z + iPos = INT((r-VariableExternalField(1,1))/DeltaExternalField(1)) + 1 ! dr = DeltaExternalField(1) + jPos = INT((z-VariableExternalField(2,1))/DeltaExternalField(2)) + 1 ! dz = DeltaExternalField(2) + + ! Catch problem when r or z are exactly at the upper boundary and INT() does not round to the lower integer (do not add +1 in + ! this case) + iPos = MIN(iPos, VariableExternalFieldN(2) - 1 ) + jPos = MIN(jPos, VariableExternalFieldN(1) - 1 ) + + + ! Shift all points by Nz = EmissionDistributionNum(1) + ASSOCIATE( Nz => VariableExternalFieldN(1) ) + ! 1.1 + idx1 = (iPos-1)*Nz + jPos + ! 2.1 + idx2 = (iPos-1)*Nz + jPos + 1 + ! 1.2 + idx3 = iPos*Nz + jPos + ! 2.2 + idx4 = iPos*Nz + jPos + 1 + END ASSOCIATE - ! Interpolate - delta = DeltaExternalField(1)*DeltaExternalField(2) - delta = 1./delta + ! Interpolate + delta = DeltaExternalField(1)*DeltaExternalField(2) + delta = 1./delta - dx(1) = VariableExternalField(1,idx4)-r - dx(2) = DeltaExternalField(1) - dx(1) + dx(1) = VariableExternalField(1,idx4)-r + dx(2) = DeltaExternalField(1) - dx(1) - dy(1) = VariableExternalField(2,idx4)-z - dy(2) = DeltaExternalField(2) - dy(1) + dy(1) = VariableExternalField(2,idx4)-z + dy(2) = DeltaExternalField(2) - dy(1) - DO i = 1, 2 - mat(1,1) = VariableExternalField(2+i,idx1) - mat(2,1) = VariableExternalField(2+i,idx2) - mat(1,2) = VariableExternalField(2+i,idx3) - mat(2,2) = VariableExternalField(2+i,idx4) + DO i = 1, 2 + mat(1,1) = VariableExternalField(2+i,idx1) + mat(2,1) = VariableExternalField(2+i,idx2) + mat(1,2) = VariableExternalField(2+i,idx3) + mat(2,2) = VariableExternalField(2+i,idx4) - vec(1) = dx(1) - vec(2) = dx(2) - f(i) = delta * DOT_PRODUCT( vec, MATMUL(mat, (/dy(1),dy(2)/) ) ) - END DO ! i = 1, 2 + vec(1) = dx(1) + vec(2) = dx(2) + f(i) = delta * DOT_PRODUCT( vec, MATMUL(mat, (/dy(1),dy(2)/) ) ) + END DO ! i = 1, 2 - ! Transform from Br, Bz to Bx, By, Bz - r=1./r + ! Transform from Br, Bz to Bx, By, Bz + r=1./r + IF(Symmetry%Axisymmetric) THEN + InterpolateVariableExternalField2D(1) = f(2) + InterpolateVariableExternalField2D(2) = f(1)*y*r + InterpolateVariableExternalField2D(3) = 0. + ELSE InterpolateVariableExternalField2D(1) = f(1)*x*r InterpolateVariableExternalField2D(2) = f(1)*y*r InterpolateVariableExternalField2D(3) = f(2) - END IF ! r.GT.VariableExternalFieldMax(1) -END ASSOCIATE + END IF +END IF ! r.GT.VariableExternalFieldMax(1) END FUNCTION InterpolateVariableExternalField2D diff --git a/src/particles/surfacemodel/surfacemodel_porous.f90 b/src/particles/surfacemodel/surfacemodel_porous.f90 index ff5b73c1f..ff0f9da86 100644 --- a/src/particles/surfacemodel/surfacemodel_porous.f90 +++ b/src/particles/surfacemodel/surfacemodel_porous.f90 @@ -98,7 +98,8 @@ SUBROUTINE InitPorousBoundaryCondition() USE MOD_ReadInTools USE MOD_Particle_Boundary_Tools ,ONLY: GetRadialDistance2D USE MOD_Particle_Mesh_Vars ,ONLY: SideInfo_Shared -USE MOD_Particle_Vars ,ONLY: Symmetry, VarTimeStep +USE MOD_Symmetry_Vars ,ONLY: Symmetry +USE MOD_Particle_Vars ,ONLY: VarTimeStep USE MOD_SurfaceModel_Vars ,ONLY: nPorousBC, PorousBC USE MOD_Particle_Boundary_Vars ,ONLY: PartBound, nPorousSides, MapSurfSideToPorousSide_Shared, PorousBCSampWall USE MOD_Particle_Boundary_Vars ,ONLY: PorousBCInfo_Shared,PorousBCProperties_Shared,PorousBCSampWall_Shared diff --git a/src/particles/surfacemodel/surfacemodel_tools.f90 b/src/particles/surfacemodel/surfacemodel_tools.f90 index a99812ed2..d9b5789e9 100644 --- a/src/particles/surfacemodel/surfacemodel_tools.f90 +++ b/src/particles/surfacemodel/surfacemodel_tools.f90 @@ -285,19 +285,18 @@ SUBROUTINE DiffuseReflection(PartID,SideID,n_loc) ! MODULES ! !----------------------------------------------------------------------------------------------------------------------------------! USE MOD_Globals +USE MOD_Globals_Vars ,ONLY: TwoepsMach USE MOD_Particle_Mesh_Vars -USE MOD_Globals ,ONLY: ABORT, OrthoNormVec, VECNORM, DOTPRODUCT USE MOD_DSMC_Vars ,ONLY: DSMC, AmbipolElecVelo USE MOD_Particle_Boundary_Vars ,ONLY: PartBound -USE MOD_Particle_Vars ,ONLY: PartState,LastPartPos,Species,PartSpecies,Symmetry,PartVeloRotRef +USE MOD_Particle_Vars ,ONLY: PartState,LastPartPos,Species,PartSpecies,PartVeloRotRef USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, PartTimeStep, VarTimeStep USE MOD_TimeDisc_Vars ,ONLY: dt,RKdtFrac USE MOD_Mesh_Tools ,ONLY: GetCNElemID USE MOD_Particle_Vars ,ONLY: PDM, RotRefFrameOmega,RotRefSubTimeStep,nSubCyclingSteps USE MOD_Particle_Tracking_Vars ,ONLY: TrackInfo USE MOD_part_RHS ,ONLY: CalcPartRHSRotRefFrame -USE MOD_Globals_Vars ,ONLY: TwoepsMach -USE MOD_Part_Tools ,ONLY: CalcPartSymmetryPos +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------! diff --git a/src/particles/tracking/particle_intersection.f90 b/src/particles/tracking/particle_intersection.f90 index 34e13353b..498cf8d8a 100644 --- a/src/particles/tracking/particle_intersection.f90 +++ b/src/particles/tracking/particle_intersection.f90 @@ -274,7 +274,7 @@ END SUBROUTINE ParticleThroughSideCheck3DFast SUBROUTINE InitParticleThroughSideCheck1D2D() ! MODULES USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !================================================================================================================================== diff --git a/src/particles/tracking/particle_triatracking.f90 b/src/particles/tracking/particle_triatracking.f90 index 014164dd4..13d6b2feb 100644 --- a/src/particles/tracking/particle_triatracking.f90 +++ b/src/particles/tracking/particle_triatracking.f90 @@ -44,7 +44,7 @@ SUBROUTINE SingleParticleTriaTrackingInterface(i,IsInterPlanePart) SUBROUTINE InitSingleParticleTriaTracking() ! MODULES USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !================================================================================================================================== @@ -267,7 +267,7 @@ SUBROUTINE SingleParticleTriaTracking3D(i,IsInterPlanePart) USE MOD_Mesh_Tools ,ONLY: GetCNElemID USE MOD_Particle_Vars ,ONLY: PEM,PDM,PartSpecies USE MOD_Particle_Vars ,ONLY: PartState,LastPartPos -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Particle_Vars ,ONLY: UseRotRefFrame, RotRefFrameOmega, PartVeloRotRef USE MOD_Particle_Mesh_Tools ,ONLY: ParticleInsideQuad3D USE MOD_Particle_Mesh_Vars diff --git a/src/posti/piclas2vtk/piclas2vtk.f90 b/src/posti/piclas2vtk/piclas2vtk.f90 index 11c562fd1..133786c97 100644 --- a/src/posti/piclas2vtk/piclas2vtk.f90 +++ b/src/posti/piclas2vtk/piclas2vtk.f90 @@ -45,7 +45,7 @@ PROGRAM piclas2vtk #ifdef PARTICLES USE MOD_Particle_Mesh_Tools ,ONLY: InitParticleGeometry USE MOD_Particle_Mesh_Vars ,ONLY: ConcaveElemSide_Shared,ElemSideNodeID_Shared,ElemMidPoint_Shared -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #endif /*PARTICLES*/ USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared, ElemNodeID_Shared, NodeInfo_Shared USE MOD_Mesh_Tools ,ONLY: InitGetCNElemID, InitGetGlobalElemID diff --git a/src/posti/superB/superB.f90 b/src/posti/superB/superB.f90 index 971b86f1e..3a71febf5 100644 --- a/src/posti/superB/superB.f90 +++ b/src/posti/superB/superB.f90 @@ -39,7 +39,7 @@ PROGRAM SuperB_standalone USE MOD_Interpolation_Vars ,ONLY: BGField,BGFieldAnalytic USE MOD_Mesh ,ONLY: InitMesh #ifdef PARTICLES -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #endif /*PARTICLES*/ #if USE_MPI USE MOD_MPI_Shared diff --git a/src/radiation/radiative_transfer/radtrans_init.f90 b/src/radiation/radiative_transfer/radtrans_init.f90 index bd008f564..4c1fcdb3b 100644 --- a/src/radiation/radiative_transfer/radtrans_init.f90 +++ b/src/radiation/radiative_transfer/radtrans_init.f90 @@ -531,7 +531,7 @@ SUBROUTINE ElemInObsCone(ElemID, ElemInCone) USE MOD_Mesh_Tools ,ONLY: GetGlobalElemID USE MOD_Particle_Mesh_Vars ,ONLY: NodeCoords_Shared,ElemInfo_Shared, BoundsOfElem_Shared, SideIsSymSide, SideInfo_Shared USE MOD_RadiationTrans_Vars ,ONLY: RadObservationPoint, RadTransObsVolumeFrac -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Particle_Mesh_Tools ,ONLY: ParticleInsideQuad USE MOD_Photon_TrackingTools ,ONLY: PhotonIntersectionWithSide2DDir, PhotonThroughSideCheck3DDir ! IMPLICIT VARIABLE HANDLING @@ -655,7 +655,7 @@ SUBROUTINE ElemOnLineOfSight(ElemID, ElemInCone) USE MOD_Mesh_Tools ,ONLY: GetGlobalElemID USE MOD_Particle_Mesh_Vars ,ONLY: ElemInfo_Shared, SideInfo_Shared, SideIsSymSide USE MOD_RadiationTrans_Vars ,ONLY: RadObservationPoint, RadObservationPOI -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Photon_TrackingTools ,ONLY: PhotonIntersectionWithSide2DDir, PhotonThroughSideCheck3DDir USE MOD_Particle_Boundary_Vars ,ONLY: PartBound ! IMPLICIT VARIABLE HANDLING @@ -795,10 +795,7 @@ SUBROUTINE FinalizeRadiationTransport() #if USE_MPI USE MOD_MPI_Shared_Vars ,ONLY: MPI_COMM_SHARED USE MOD_MPI_Shared -USE MOD_Particle_Mesh_Vars ,ONLY: ElemSideNodeID2D_Shared_Win,SideNormalEdge2D_Shared_Win #endif -USE MOD_Particle_Vars ,ONLY: Symmetry -USE MOD_Particle_Mesh_Vars ,ONLY: ElemSideNodeID2D_Shared,SideNormalEdge2D_Shared !----------------------------------------------------------------------------------------------------------------------------------- IMPLICIT NONE ! INPUT VARIABLES @@ -825,8 +822,6 @@ SUBROUTINE FinalizeRadiationTransport() ADEALLOCATE(Radiation_Emission_Spec_Total_Shared) ADEALLOCATE(RadTransObsVolumeFrac_Shared) ADEALLOCATE(Radiation_Emission_Spec_Max_Shared) -ADEALLOCATE(ElemSideNodeID2D_Shared) -ADEALLOCATE(SideNormalEdge2D_Shared) END SUBROUTINE FinalizeRadiationTransport diff --git a/src/radiation/radiative_transfer/radtrans_main.f90 b/src/radiation/radiative_transfer/radtrans_main.f90 index 3ad3672a1..0d60d7b00 100644 --- a/src/radiation/radiative_transfer/radtrans_main.f90 +++ b/src/radiation/radiative_transfer/radtrans_main.f90 @@ -54,7 +54,7 @@ SUBROUTINE RadTrans_main() USE MOD_Output ,ONLY: PrintStatusLineRadiation USE MOD_MPI_Shared_Vars USE MOD_MPI_Shared -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry #if USE_MPI USE MOD_RadiationTrans_Vars ,ONLY: RadTransPhotPerCell_Shared_Win #else diff --git a/src/radiation/radiative_transfer/tracking/radtrans_tracking.f90 b/src/radiation/radiative_transfer/tracking/radtrans_tracking.f90 index af4e42171..78277d658 100644 --- a/src/radiation/radiative_transfer/tracking/radtrans_tracking.f90 +++ b/src/radiation/radiative_transfer/tracking/radtrans_tracking.f90 @@ -52,9 +52,9 @@ SUBROUTINE InitPhotonSurfSample() #else USE MOD_Particle_Boundary_Vars ,ONLY: nGlobalSurfSides #endif /*USE_MPI*/ -USE MOD_Particle_Vars ,ONLY: Symmetry +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Basis ,ONLY: LegendreGaussNodesAndWeights -USE MOD_DSMC_Symmetry ,ONLY: DSMC_2D_CalcSymmetryArea, DSMC_1D_CalcSymmetryArea +USE MOD_Particle_Mesh_Tools ,ONLY: DSMC_2D_CalcSymmetryArea, DSMC_1D_CalcSymmetryArea USE MOD_Mesh_Vars ,ONLY: NGeo USE MOD_Mesh_Tools ,ONLY: GetCNElemID USE MOD_Particle_Mesh_Vars ,ONLY: SideInfo_Shared,NodeCoords_Shared diff --git a/src/symmetry/symmetry.f90 b/src/symmetry/symmetry.f90 new file mode 100644 index 000000000..11b3df6a5 --- /dev/null +++ b/src/symmetry/symmetry.f90 @@ -0,0 +1,105 @@ +!================================================================================================================================== +! 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_Symmetry +!=================================================================================================================================== +!> Routines for 2D (planar/axisymmetric) simulations +!=================================================================================================================================== +! MODULES +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +PRIVATE + +!----------------------------------------------------------------------------------------------------------------------------------- +! GLOBAL VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! Private Part --------------------------------------------------------------------------------------------------------------------- +! Public Part ---------------------------------------------------------------------------------------------------------------------- +PUBLIC :: DefineParametersSymmetry, InitSymmetry +!=================================================================================================================================== + +CONTAINS + +!================================================================================================================================== +!> Define parameters for particles +!================================================================================================================================== +SUBROUTINE DefineParametersSymmetry() +! MODULES +USE MOD_ReadInTools ,ONLY: prms +IMPLICIT NONE + +CALL prms%SetSection("Particle Symmetry") +CALL prms%CreateIntOption( 'Particles-Symmetry-Order', & + 'Order of the Simulation 1, 2 or 3 D', '3') +CALL prms%CreateLogicalOption('Particles-Symmetry2DAxisymmetric', 'Activating an axisymmetric simulation with the same mesh '//& + 'requirements as for the 2D case (y is then the radial direction)', '.FALSE.') + +END SUBROUTINE DefineParametersSymmetry + + +SUBROUTINE InitSymmetry() +!=================================================================================================================================== +!> Initialize if a 2D/1D Simulation is performed and which type +!=================================================================================================================================== +! MODULES +USE MOD_Globals +USE MOD_ReadInTools ,ONLY: GETLOGICAL,GETINT +USE MOD_Symmetry_Vars ,ONLY: Symmetry +#if defined(PARTICLES) +USE MOD_Particle_Mesh_Tools ,ONLY: InitParticleInsideQuad +USE MOD_Particle_TriaTracking ,ONLY: InitSingleParticleTriaTracking +USE MOD_Particle_InterSection ,ONLY: InitParticleThroughSideCheck1D2D +USE MOD_DSMC_Vars ,ONLY: RadialWeighting +#endif /*defined(PARTICLES)*/ +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- + +!=================================================================================================================================== +Symmetry%Order = GETINT('Particles-Symmetry-Order') + +IF((Symmetry%Order.LE.0).OR.(Symmetry%Order.GE.4)) CALL ABORT(__STAMP__& + ,'Particles-Symmetry-Order (space dimension) has to be in the range of 1 to 3') + +#if defined(PARTICLES) +! Initialize the function pointers for triatracking +CALL InitParticleInsideQuad() +CALL InitSingleParticleTriaTracking() +IF(Symmetry%Order.LE.2) CALL InitParticleThroughSideCheck1D2D() +#endif /*defined(PARTICLES)*/ + +Symmetry%Axisymmetric = GETLOGICAL('Particles-Symmetry2DAxisymmetric') +IF(Symmetry%Axisymmetric.AND.(Symmetry%Order.EQ.3)) CALL ABORT(__STAMP__& + ,'ERROR: Axisymmetric simulations only for 1D or 2D') +IF(Symmetry%Axisymmetric.AND.(Symmetry%Order.EQ.1))CALL ABORT(__STAMP__& + ,'ERROR: Axisymmetric simulations are only implemented for Particles-Symmetry-Order=2 !') + +#if defined(PARTICLES) +IF(Symmetry%Axisymmetric) THEN + RadialWeighting%DoRadialWeighting = GETLOGICAL('Particles-RadialWeighting') +ELSE + RadialWeighting%DoRadialWeighting = .FALSE. + RadialWeighting%PerformCloning = .FALSE. +END IF +#endif /*defined(PARTICLES)*/ + +END SUBROUTINE InitSymmetry + +END MODULE MOD_Symmetry \ No newline at end of file diff --git a/src/symmetry/symmetry_vars.f90 b/src/symmetry/symmetry_vars.f90 new file mode 100644 index 000000000..de6c6c7af --- /dev/null +++ b/src/symmetry/symmetry_vars.f90 @@ -0,0 +1,35 @@ +!================================================================================================================================== +! 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_Symmetry_Vars +!=================================================================================================================================== +! Contains the Particles' variables (general for all modules: PIC, DSMC, FP) +!=================================================================================================================================== +! MODULES +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +PUBLIC +SAVE +!----------------------------------------------------------------------------------------------------------------------------------- +! GLOBAL VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- + +TYPE tSymmetry +INTEGER :: Order ! 1-3 D +LOGICAL :: Axisymmetric +END TYPE tSymmetry + +TYPE(tSymmetry) :: Symmetry + +END MODULE MOD_Symmetry_Vars \ No newline at end of file diff --git a/src/timedisc/timedisc_TimeStepPoisson.f90 b/src/timedisc/timedisc_TimeStepPoisson.f90 index a3b72ae15..5bf8c6def 100644 --- a/src/timedisc/timedisc_TimeStepPoisson.f90 +++ b/src/timedisc/timedisc_TimeStepPoisson.f90 @@ -62,7 +62,7 @@ SUBROUTINE TimeStepPoisson() #if USE_MPI USE MOD_Particle_MPI ,ONLY: IRecvNbOfParticles, MPIParticleSend,MPIParticleRecv,SendNbOfparticles #endif -USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition,isPushParticle +USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition,isPushParticle,CalcPartSymmetryPos USE MOD_Particle_Tracking ,ONLY: PerformTracking USE MOD_vMPF ,ONLY: SplitAndMerge USE MOD_Particle_Vars ,ONLY: UseSplitAndMerge diff --git a/src/timedisc/timedisc_TimeStepPoissonByBorisLeapfrog.f90 b/src/timedisc/timedisc_TimeStepPoissonByBorisLeapfrog.f90 index 3202d8148..6deb0bf1c 100644 --- a/src/timedisc/timedisc_TimeStepPoissonByBorisLeapfrog.f90 +++ b/src/timedisc/timedisc_TimeStepPoissonByBorisLeapfrog.f90 @@ -60,7 +60,7 @@ SUBROUTINE TimeStepPoissonByBorisLeapfrog() #if USE_MPI USE MOD_Particle_MPI ,ONLY: IRecvNbOfParticles, MPIParticleSend,MPIParticleRecv,SendNbOfparticles #endif -USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition,isPushParticle +USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition,isPushParticle,CalcPartSymmetryPos USE MOD_Particle_Tracking ,ONLY: PerformTracking USE MOD_vMPF ,ONLY: SplitAndMerge USE MOD_Particle_Vars ,ONLY: UseSplitAndMerge @@ -190,7 +190,6 @@ SUBROUTINE TimeStepPoissonByBorisLeapfrog() !-- x(n) => x(n+1) by v(n+0.5): PartState(1:3,iPart) = PartState(1:3,iPart) + PartState(4:6,iPart) * dtFrac - CALL CalcPartSymmetryPos(PartState(1:3,iPart),PartState(4:6,iPart)) ! If coupled power output is active and particle carries charge, calculate energy difference and add to output variable IF (CalcCoupledPower) CALL CalcCoupledPowerPart(iPart,'after') diff --git a/src/timedisc/timedisc_TimeStepPoissonByHigueraCary.f90 b/src/timedisc/timedisc_TimeStepPoissonByHigueraCary.f90 index 24c742a6a..8a6e19296 100644 --- a/src/timedisc/timedisc_TimeStepPoissonByHigueraCary.f90 +++ b/src/timedisc/timedisc_TimeStepPoissonByHigueraCary.f90 @@ -62,7 +62,7 @@ SUBROUTINE TimeStepPoissonByHigueraCary() USE MOD_Particle_MPI ,ONLY: IRecvNbOfParticles, MPIParticleSend,MPIParticleRecv,SendNbOfparticles USE MOD_Particle_MPI_Vars ,ONLY: PartMPIExchange #endif -USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition,isPushParticle +USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition,isPushParticle,CalcPartSymmetryPos USE MOD_Particle_Tracking ,ONLY: PerformTracking USE MOD_vMPF ,ONLY: SplitAndMerge USE MOD_Particle_Vars ,ONLY: UseSplitAndMerge @@ -114,6 +114,7 @@ SUBROUTINE TimeStepPoissonByHigueraCary() ! 1st part of the position update (also neutral particles) !-- x(n) => x(n+1/2) by v(n) = u(n)/gamma(n): PartState(1:3,iPart) = PartState(1:3,iPart) + 0.5 * PartState(4:6,iPart) * dt + CALL CalcPartSymmetryPos(PartState(1:3,iPart),PartState(4:6,iPart)) END IF ! PDM%ParticleInside(iPart) END DO ! iPart=1,PDM%ParticleVecLength #if USE_LOADBALANCE diff --git a/src/timedisc/timedisc_TimeStepPoissonByLSERK.f90 b/src/timedisc/timedisc_TimeStepPoissonByLSERK.f90 index eaf2c8af7..bb6bd71fd 100644 --- a/src/timedisc/timedisc_TimeStepPoissonByLSERK.f90 +++ b/src/timedisc/timedisc_TimeStepPoissonByLSERK.f90 @@ -63,7 +63,7 @@ SUBROUTINE TimeStepPoissonByLSERK() USE MOD_Particle_MPI ,ONLY: IRecvNbOfParticles, MPIParticleSend,MPIParticleRecv,SendNbOfparticles #endif USE MOD_Particle_Localization ,ONLY: CountPartsPerElem -USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition,isPushParticle +USE MOD_Part_Tools ,ONLY: UpdateNextFreePosition,isPushParticle,CalcPartSymmetryPos USE MOD_Particle_Tracking ,ONLY: PerformTracking USE MOD_vMPF ,ONLY: SplitAndMerge USE MOD_Particle_Vars ,ONLY: UseSplitAndMerge @@ -204,9 +204,7 @@ SUBROUTINE TimeStepPoissonByLSERK() PDM%dtFracPush(iPart) = .FALSE. IF (.NOT.DoForceFreeSurfaceFlux) PDM%IsNewPart(iPart) = .FALSE. !change to false: Pt_temp is now rebuilt... END IF !IsNewPart - CALL CalcPartSymmetryPos(PartState(1:3,iPart),PartState(4:6,iPart)) - ! If coupled power output is active and particle carries charge, calculate energy difference and add to output variable IF (CalcCoupledPower) CALL CalcCoupledPowerPart(iPart,'after') END IF @@ -334,9 +332,7 @@ SUBROUTINE TimeStepPoissonByLSERK() END IF IF (.NOT.DoForceFreeSurfaceFlux .OR. iStage.EQ.nRKStages) PDM%IsNewPart(iPart) = .FALSE. !change to false: Pt_temp is now rebuilt... END IF !IsNewPart - CALL CalcPartSymmetryPos(PartState(1:3,iPart),PartState(4:6,iPart)) - ! If coupled power output is active and particle carries charge, calculate energy difference and add to output variable IF (CalcCoupledPower) CALL CalcCoupledPowerPart(iPart,'after') END IF diff --git a/src/timedisc/timedisc_TimeStep_BGK.f90 b/src/timedisc/timedisc_TimeStep_BGK.f90 index 2abd4efb8..befa46804 100644 --- a/src/timedisc/timedisc_TimeStep_BGK.f90 +++ b/src/timedisc/timedisc_TimeStep_BGK.f90 @@ -36,7 +36,7 @@ SUBROUTINE TimeStep_BGK() USE MOD_PreProc USE MOD_TimeDisc_Vars ,ONLY: dt, IterDisplayStep, iter, TEnd, Time USE MOD_Globals ,ONLY: abort, CROSS -USE MOD_Particle_Vars ,ONLY: PartState, LastPartPos, PDM, PEM, DoSurfaceFlux, WriteMacroVolumeValues, Symmetry +USE MOD_Particle_Vars ,ONLY: PartState, LastPartPos, PDM, PEM, DoSurfaceFlux, WriteMacroVolumeValues USE MOD_Particle_Vars ,ONLY: UseRotRefFrame, RotRefFrameOmega, PartVeloRotRef USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, PartTimeStep USE MOD_DSMC_Vars ,ONLY: DSMC, CollisMode @@ -61,6 +61,7 @@ SUBROUTINE TimeStep_BGK() USE MOD_SurfaceModel_Vars ,ONLY: nPorousBC, DoChemSurface USE MOD_Particle_SurfChemFlux USE MOD_DSMC_ParticlePairing ,ONLY: GeoCoordToMap2D +USE MOD_Symmetry_Vars ,ONLY: Symmetry ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- diff --git a/src/timedisc/timedisc_TimeStep_DSMC.f90 b/src/timedisc/timedisc_TimeStep_DSMC.f90 index 321a10a91..dc4afed01 100644 --- a/src/timedisc/timedisc_TimeStep_DSMC.f90 +++ b/src/timedisc/timedisc_TimeStep_DSMC.f90 @@ -39,7 +39,7 @@ SUBROUTINE TimeStep_DSMC() USE MOD_Globals ,ONLY: abort, CROSS USE MOD_Particle_Vars ,ONLY: PartState, LastPartPos, PDM, PEM, DoSurfaceFlux, WriteMacroVolumeValues USE MOD_Particle_Vars ,ONLY: UseRotRefFrame, RotRefFrameOmega, PartVeloRotRef, LastPartVeloRotRef -USE MOD_Particle_Vars ,ONLY: WriteMacroSurfaceValues, Symmetry, Species, PartSpecies +USE MOD_Particle_Vars ,ONLY: WriteMacroSurfaceValues, Species, PartSpecies USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, PartTimeStep, VarTimeStep USE MOD_Particle_Vars ,ONLY: UseSplitAndMerge USE MOD_DSMC_Vars ,ONLY: DSMC, CollisMode, AmbipolElecVelo @@ -54,7 +54,8 @@ SUBROUTINE TimeStep_DSMC() USE MOD_SurfaceModel_Porous ,ONLY: PorousBoundaryRemovalProb_Pressure USE MOD_SurfaceModel_Vars ,ONLY: nPorousBC, DoChemSurface USE MOD_vMPF ,ONLY: SplitAndMerge -USE MOD_part_RHS ,ONLY: CalcPartRHSRotRefFrame, CalcPartPosInRotRef +USE MOD_Symmetry_Vars ,ONLY: Symmetry +USE MOD_part_RHS ,ONLY: CalcPartPosInRotRef USE MOD_part_pos_and_velo ,ONLY: SetParticleVelocity USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck USE MOD_Part_Tools ,ONLY: CalcPartSymmetryPos diff --git a/src/timedisc/timedisc_TimeStep_FPFlow.f90 b/src/timedisc/timedisc_TimeStep_FPFlow.f90 index c4bc38671..b09b705b7 100644 --- a/src/timedisc/timedisc_TimeStep_FPFlow.f90 +++ b/src/timedisc/timedisc_TimeStep_FPFlow.f90 @@ -44,6 +44,7 @@ SUBROUTINE TimeStep_FPFlow() USE MOD_Particle_SurfFlux ,ONLY: ParticleSurfaceflux USE MOD_Particle_Tracking ,ONLY: PerformTracking USE MOD_Particle_Tracking_vars ,ONLY: tTracking,MeasureTrackTime +USE MOD_Symmetry_Vars ,ONLY: Symmetry USE MOD_Restart_Vars ,ONLY: DoRestart USE MOD_Part_Tools ,ONLY: CalcPartSymmetryPos #if USE_MPI