diff --git a/REGGIE.md b/REGGIE.md index eb08bcf72..9d6f54280 100644 --- a/REGGIE.md +++ b/REGGIE.md @@ -58,6 +58,7 @@ Small test cases to check features with DSMC timedisc: [Link to build](regressio | | Rotational_Reference_Frame_Regions | | Rotational reference frame with several regions, switching between stationary and rotating frame | nProcs=1,2,3,4 | Particle trajectory | [Link](regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Regions/readme.md) | | | Rotational_Reference_Frame_RotBC | | Rotational reference frame in combination with the rotationally periodic BC | nProcs=1,2,3,4 | Particle trajectory | [Link](regressioncheck/CHE_DSMC/Rotational_Reference_Frame_RotBC/readme.md) | | | Rotational_Reference_Frame_Temperature | | Rotational reference frame: Many particles, multiple revolutions | nProcs=1,2,4 | Temperature | [Link](regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Temperature/readme.md) | +| | Rotational_Reference_Frame_Subcycling | | Time subcycling method within the rotational reference frame | nProcs=1 | Particle trajectory | [Link](regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/readme.md | | | SurfaceOutput | | Test of CalcSurfaceImpact and CalcBoundaryParticleOutput through defined electron flux | nProcs=1,4 | PartAnalyze, SurfaceAnalyze, DSMCSurfState | [Link](regressioncheck/CHE_DSMC/SurfaceOutput/readme.md) | | | DSMC_QualityFactors | | Quality factors: mean/max collision probability, MCS over MFP, mean free path, ResolvedCellPercentage | nProcs=1 | PartAnalyze | [Link](regressioncheck/CHE_DSMC/DSMC_QualityFactors/readme.md) | | | DSMC_QualityFactors_MPI | | Quality factors: ResolvedTimestep, max collision probability, MCS over MFP, ResolvedCellPercentage | nProcs=4 | PartAnalyze | [Link](regressioncheck/CHE_DSMC/DSMC_QualityFactors_MPI/readme.md) | @@ -77,7 +78,7 @@ Small test cases to check features with DSMC timedisc: [Link to build](regressio #### CHE_BGK/FPFlow -Both methods share the same regression tests in the different folders (CHE_BGK: [BGK build](regressioncheck/CHE_BGK/builds.ini), CHE_FPFlow: [FPFlow build](regressioncheck/CHE_FPFlow/builds.ini) +Both methods share the same regression tests in the different folders, CHE_BGK: [BGK build](regressioncheck/CHE_BGK/builds.ini), CHE_FPFlow: [FPFlow build](regressioncheck/CHE_FPFlow/builds.ini) | **No.** | **Case** | **CMAKE-CONFIG** | **Feature** | **Execution** | **Comparing** | **Readme** | | :-----: | :-----------------------: | :--------------: | :---------------------------------------------------------------------------------------------------------------------: | :-----------: | :------------------------: | :-----------------------------------------------------------------: | @@ -451,38 +452,38 @@ Test all features of radiation timedisc (cell-local emission using the radiation Overview of the test cases performed every week. -| **No.** | **Case** | **CMAKE-CONFIG** | **Feature** | **Execution** | **Comparing** | **Readme** | -| :-----: | :-----------------------------------------------------: | :-------------------------------------------------------: | :------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | :---------------------------: | :-------------------------------------------------------------------------------- | :------------------------------------------------------------------------------------------------: | -| 1 | plasma_wave | [PIC-Maxwell](regressioncheck/WEK_PIC_maxwell/builds.ini) | Maxwell-PIC,SF1D, FastPeriodic | nProcs=6, IMEX for ImplicitO4 | W_el LineIntegration (FieldAnalyze.csv) | [Link](regressioncheck/WEK_PIC_maxwell/plasma_wave/readme.md) | -| ** | 3D_periodic_shape_function | ** | Maxwell-PIC,shape function deposition over periodic sides 3D | nProcs= 1,2,6,10,20 | L2 error and PartAnalyze.csv | [Link](regressioncheck/WEK_PIC_maxwell/3D_periodic_shape_function/readme.md) | -| ** | 3D_periodic_CVWM | ** | Maxwell-PIC,CVWM over periodic sides 3D with 1000 elements | nProcs= 1,2,6,10,15,20,30 | L2 error and PartAnalyze.csv | [Link](regressioncheck/WEK_PIC_maxwell/3D_periodic_CVWM/readme.md) | -| ** | 3D_periodic_CVWM_split2hex | ** | Maxwell-PIC,CVWM over periodic sides 3D and split2hex grid with 768 hex elements | nProcs= 1,2,6,10,15,20,30 | L2 error and PartAnalyze.csv | [Link](regressioncheck/WEK_PIC_maxwell/3D_periodic_CVWM_split2hex/readme.md) | -| 2 | HEMPT-90deg-symmetry | [PIC-HDG](regressioncheck/WEK_PIC_poisson/builds.ini) | create mesh (hopr) and external magnetic field (superB) and use both in simulation | nProcs=1,10,20 | | [Link](regressioncheck/WEK_PIC_poisson/HEMPT-90deg-symmetry/readme.md) | -| 3 | CHEM_EQUI_diss_CH4 | [Reservoir](regressioncheck/WEK_Reservoir/builds.ini) | Relaxation into equilibrium with dissociation and recombination of CH4 | nProcs=2 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_diss_CH4/readme.md) | -| ** | CHEM_EQUI_exch_CH3-H | ** | Relaxation into equilibrium with exchange/radical reaction of CH3+H <-> CH2+H2 | nProcs=2 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_exch_CH3-H/readme.md) | -| ** | CHEM_EQUI_ionization_H | ** | Relaxation into equilibrium with ionization and recombination of H | nProcs=1 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_ionization_H/readme.md) | -| ** | CHEM_EQUI_diss_CH4_2DAxi_RadWeight | ** | Analogous to CHEM_EQUI_diss_CH4 with 2D axisymmetric mesh with radial weighting | nProcs=2 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_diss_CH4_2DAxi_RadWeight/readme.md) | -| ** | CHEM_EQUI_Titan_Chemistry | ** | Reservoir simulation with Titan's atmosphere (18 species, 28 reactions) | nProcs=6 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_Titan_Chemistry/readme.md) | -| ** | CHEM_EQUI_Titan_Chemistry_Database | ** | Reservoir simulation with Titan's atmosphere (18 species, 28 reactions) using species/reaction data from the species database | nProcs=6 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_Titan_Chemistry_Database/readme.md) | -| ** | MCC_MultiSpec_XSec | ** | Multi-species reservoir: Collision rates for neutral-electrons through cross-section data | nProcs=1 | | [Link](regressioncheck/WEK_Reservoir/MCC_MultiSpec_XSec/readme.md) | -| ** | MCC_MultiSpec_XSec_TCE_QK_Chem | ** | Multi-species reservoir: QK ionization and TCE dissociation | nProcs=2 | | [Link](regressioncheck/WEK_Reservoir/MCC_MultiSpec_XSec_Chem/readme.md) | -| ** | BGG_MultiSpec_XSec_Elec | ** | Background gas reservoir with VHS: Electronic excitation rates for neutral-electrons through cross-section data | nProcs=1 | | [Link](regressioncheck/WEK_Reservoir/BGG_MultiSpec_XSec_Elec/readme.md) | -| ** | MCC_N2_XSec_Elec | ** | Regular reservoir with MCC/VHS: Electronic excitation rates for neutral-electrons through cross-section data | nProcs=1 | | [Link](regressioncheck/WEK_Reservoir/MCC_N2_XSec_Elec/readme.md) | -| ** | 1D_Sod_Shocktube | [DSMC](regressioncheck/WEK_DSMC/builds.ini) | 1D test case shock tube | nProcs=6 | DSMCState | [Link](regressioncheck/WEK_Reservoir/1D_Sod_Shocktube/readme.md) | -| 4 | 2DAxi_ChannelFlow_ConstPressure_TruncAverage | ** | 2D axisymmetric: Pressure gradient driven pipe flow with adaptive surface flux, using a truncated running average | nProcs=6 | PartAnalyze: Average pressure and mass flow rate at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/2DAxi_ChannelFlow_ConstPressure_TruncAverage/readme.md) | -| 4 | ChannelFlow_AdaptiveBoundary_ConstMassflow | ** | Constant massflow driven channel flow with adaptive surface flux | nProcs=6 | PartAnalyze: Average pressure and mass flow rate at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/ChannelFlow_AdaptiveBoundary_ConstMassflow/readme.md) | -| 4 | ChannelFlow_AdaptiveBoundary_ConstPressure_FixedAverage | ** | Pressure gradient driven channel flow with adaptive surface flux, using a fixed average for the sampling | nProcs=6 | PartAnalyze: Average pressure at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/ChannelFlow_AdaptiveBoundary_ConstPressure_FixedAverage/readme.md) | -| 4 | ChannelFlow_AdaptiveBoundary_ConstPressure_Relaxation | ** | Pressure gradient driven channel flow with adaptive surface flux, using a relaxation factor for the sampling | nProcs=6 | PartAnalyze: Average pressure at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/ChannelFlow_AdaptiveBoundary_ConstPressure_Relaxation/readme.md) | -| 4 | ChannelFlow_AdaptiveBoundary_ConstPressure_TruncAverage | ** | Pressure gradient driven channel flow with adaptive surface flux, using a truncated running average for the sampling | nProcs=6 | PartAnalyze: Average pressure at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/ChannelFlow_AdaptiveBoundary_ConstPressure_TruncAverage/readme.md) | +| **No.** | **Case** | **CMAKE-CONFIG** | **Feature** | **Execution** | **Comparing** | **Readme** | +| :-----: | :-----------------------------------------------------: | :-------------------------------------------------------: | :------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | :---------------------------: | :--------------------------------------------------------------------------------- | :------------------------------------------------------------------------------------------------: | +| 1 | plasma_wave | [PIC-Maxwell](regressioncheck/WEK_PIC_maxwell/builds.ini) | Maxwell-PIC,SF1D, FastPeriodic | nProcs=6, IMEX for ImplicitO4 | W_el LineIntegration (FieldAnalyze.csv) | [Link](regressioncheck/WEK_PIC_maxwell/plasma_wave/readme.md) | +| ** | 3D_periodic_shape_function | ** | Maxwell-PIC,shape function deposition over periodic sides 3D | nProcs= 1,2,6,10,20 | L2 error and PartAnalyze.csv | [Link](regressioncheck/WEK_PIC_maxwell/3D_periodic_shape_function/readme.md) | +| ** | 3D_periodic_CVWM | ** | Maxwell-PIC,CVWM over periodic sides 3D with 1000 elements | nProcs= 1,2,6,10,15,20,30 | L2 error and PartAnalyze.csv | [Link](regressioncheck/WEK_PIC_maxwell/3D_periodic_CVWM/readme.md) | +| ** | 3D_periodic_CVWM_split2hex | ** | Maxwell-PIC,CVWM over periodic sides 3D and split2hex grid with 768 hex elements | nProcs= 1,2,6,10,15,20,30 | L2 error and PartAnalyze.csv | [Link](regressioncheck/WEK_PIC_maxwell/3D_periodic_CVWM_split2hex/readme.md) | +| 2 | HEMPT-90deg-symmetry | [PIC-HDG](regressioncheck/WEK_PIC_poisson/builds.ini) | create mesh (hopr) and external magnetic field (superB) and use both in simulation | nProcs=1,10,20 | | [Link](regressioncheck/WEK_PIC_poisson/HEMPT-90deg-symmetry/readme.md) | +| 3 | CHEM_EQUI_diss_CH4 | [Reservoir](regressioncheck/WEK_Reservoir/builds.ini) | Relaxation into equilibrium with dissociation and recombination of CH4 | nProcs=2 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_diss_CH4/readme.md) | +| ** | CHEM_EQUI_exch_CH3-H | ** | Relaxation into equilibrium with exchange/radical reaction of CH3+H <-> CH2+H2 | nProcs=2 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_exch_CH3-H/readme.md) | +| ** | CHEM_EQUI_ionization_H | ** | Relaxation into equilibrium with ionization and recombination of H | nProcs=1 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_ionization_H/readme.md) | +| ** | CHEM_EQUI_diss_CH4_2DAxi_RadWeight | ** | Analogous to CHEM_EQUI_diss_CH4 with 2D axisymmetric mesh with radial weighting | nProcs=2 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_diss_CH4_2DAxi_RadWeight/readme.md) | +| ** | CHEM_EQUI_Titan_Chemistry | ** | Reservoir simulation with Titan's atmosphere (18 species, 28 reactions) | nProcs=6 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_Titan_Chemistry/readme.md) | +| ** | CHEM_EQUI_Titan_Chemistry_Database | ** | Reservoir simulation with Titan's atmosphere (18 species, 28 reactions) using species/reaction data from the species database | nProcs=6 | PartAnalyze_ref.csv | [Link](regressioncheck/WEK_Reservoir/CHEM_EQUI_Titan_Chemistry_Database/readme.md) | +| ** | MCC_MultiSpec_XSec | ** | Multi-species reservoir: Collision rates for neutral-electrons through cross-section data | nProcs=1 | | [Link](regressioncheck/WEK_Reservoir/MCC_MultiSpec_XSec/readme.md) | +| ** | MCC_MultiSpec_XSec_TCE_QK_Chem | ** | Multi-species reservoir: QK ionization and TCE dissociation | nProcs=2 | | [Link](regressioncheck/WEK_Reservoir/MCC_MultiSpec_XSec_Chem/readme.md) | +| ** | BGG_MultiSpec_XSec_Elec | ** | Background gas reservoir with VHS: Electronic excitation rates for neutral-electrons through cross-section data | nProcs=1 | | [Link](regressioncheck/WEK_Reservoir/BGG_MultiSpec_XSec_Elec/readme.md) | +| ** | MCC_N2_XSec_Elec | ** | Regular reservoir with MCC/VHS: Electronic excitation rates for neutral-electrons through cross-section data | nProcs=1 | | [Link](regressioncheck/WEK_Reservoir/MCC_N2_XSec_Elec/readme.md) | +| ** | 1D_Sod_Shocktube | [DSMC](regressioncheck/WEK_DSMC/builds.ini) | 1D test case shock tube | nProcs=6 | DSMCState | [Link](regressioncheck/WEK_Reservoir/1D_Sod_Shocktube/readme.md) | +| 4 | 2DAxi_ChannelFlow_ConstPressure_TruncAverage | ** | 2D axisymmetric: Pressure gradient driven pipe flow with adaptive surface flux, using a truncated running average | nProcs=6 | PartAnalyze: Average pressure and mass flow rate at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/2DAxi_ChannelFlow_ConstPressure_TruncAverage/readme.md) | +| 4 | ChannelFlow_AdaptiveBoundary_ConstMassflow | ** | Constant massflow driven channel flow with adaptive surface flux | nProcs=6 | PartAnalyze: Average pressure and mass flow rate at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/ChannelFlow_AdaptiveBoundary_ConstMassflow/readme.md) | +| 4 | ChannelFlow_AdaptiveBoundary_ConstPressure_FixedAverage | ** | Pressure gradient driven channel flow with adaptive surface flux, using a fixed average for the sampling | nProcs=6 | PartAnalyze: Average pressure at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/ChannelFlow_AdaptiveBoundary_ConstPressure_FixedAverage/readme.md) | +| 4 | ChannelFlow_AdaptiveBoundary_ConstPressure_Relaxation | ** | Pressure gradient driven channel flow with adaptive surface flux, using a relaxation factor for the sampling | nProcs=6 | PartAnalyze: Average pressure at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/ChannelFlow_AdaptiveBoundary_ConstPressure_Relaxation/readme.md) | +| 4 | ChannelFlow_AdaptiveBoundary_ConstPressure_TruncAverage | ** | Pressure gradient driven channel flow with adaptive surface flux, using a truncated running average for the sampling | nProcs=6 | PartAnalyze: Average pressure at the adaptive surface flux BCs | [Link](regressioncheck/WEK_DSMC/ChannelFlow_AdaptiveBoundary_ConstPressure_TruncAverage/readme.md) | | | ChannelFlow_SurfChem_AdsorpDesorp_CO_O2 | ** | Channel flow with surface chemistry, testing adsorption/desorption of CO and O2 | nProcs=6 | Coverage (DSMCSurfChemStatePartAnalyze), Number density, temperature (PartAnalyze) | [Link](regressioncheck/WEK_DSMC/ChannelFlow_SurfChem_AdsorpDesorp_CO_O2/readme.md) | -| ** | Flow_Argon_Cylinder_Curved | ** | Hypersonic Argon flow around a cylinder (pseudo 2D) with DSMC on a curved mesh | nProcs=2 | | [Link](regressioncheck/WEK_DSMC/Flow_Argon_Cylinder_Curved/readme.md) | -| ** | Flow_Argon_Cylinder_LinearMesh | ** | Hypersonic Argon flow around a cylinder (2D) with DSMC on a linear mesh | nProcs=4 | | [Link](regressioncheck/WEK_DSMC/Flow_Argon_Cylinder_LinearMesh/readme.md) | -| ** | Flow_N2_70degCone | ** | 2D axisymmetric 70 degree cone | nProcs=6 | Surface Sampling, includes CalcSurfaceImpact and adaptive wall temperature | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | -| ** | fully_periodic_3D | ** | Periodic boundary conditions in all three directions | nProcs=10,20,30 | Check whether particles end up outside of the domain | [Link](regressioncheck/WEK_DSMC/fully_periodic_3D/readme.md) | -| ** | Surface_Sticking_Coefficient | ** | Channel flow with a sticking coefficient model | nProcs=5 | Surface sampling | [Link](regressioncheck/WEK_DSMC/Surface_Sticking_Coefficient/readme.md) | -| 5 | Flow_N2_70degCone | [BGK](regressioncheck/WEK_BGKFlow/builds.ini) | 2D axisymmetric 70 degree cone | nProcs=6 | | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | -| ** | MultiSpec_Supersonic_Couette_Ar-He | ** | Supersonic Couette flow with an Ar-He mixture | nProcs=5 | Temperature | [Link](regressioncheck/WEK_DSMC/MultiSpec_Supersonic_Couette_Ar-He/readme.md) | -| ** | MultiSpec_Supersonic_Couette_CO2-N2 | ** | Supersonic Couette flow with a CO2-N2 mixture | nProcs=5 | Temperature | [Link](regressioncheck/WEK_DSMC/MultiSpec_Supersonic_Couette_CO2-N2/readme.md) | -| 6 | Flow_N2_70degCone | [FP](regressioncheck/WEK_FPFlow/builds.ini) | 2D axisymmetric 70 degree cone | nProcs=6 | Surface Sampling, includes CalcSurfaceImpact | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | -| 7 | Flow_N2-N_70degConeHot | [DSMC](regressioncheck/WEK_DSMC/builds.ini) | 2D axisymmetric 70 degree cone (hotter and with N to get some radiation in the next step) | nProcs=6 | Surface Sampling | [Link](regressioncheck/WEK_DSMC/Flow_N2-N_70degConeHot/readme.md) | -| ** | Flow_N2-N_70degConeHot | [Radiation](regressioncheck/WEK_Radiation/builds.ini) | using previously simulated WEK_DSMC/Flow_N2_70degCone results to check radiation tool chain (write out DSMC results, readin those results, radiation solver, radiative transfer, piclas2vtk) | nProcs=6 | Surface heat flux | [Link](regressioncheck/WEK_Radiation/Flow_N2-N_70degConeHot/readme.md) | +| ** | Flow_Argon_Cylinder_Curved | ** | Hypersonic Argon flow around a cylinder (pseudo 2D) with DSMC on a curved mesh | nProcs=2 | | [Link](regressioncheck/WEK_DSMC/Flow_Argon_Cylinder_Curved/readme.md) | +| ** | Flow_Argon_Cylinder_LinearMesh | ** | Hypersonic Argon flow around a cylinder (2D) with DSMC on a linear mesh | nProcs=4 | | [Link](regressioncheck/WEK_DSMC/Flow_Argon_Cylinder_LinearMesh/readme.md) | +| ** | Flow_N2_70degCone | ** | 2D axisymmetric 70 degree cone | nProcs=6 | Surface Sampling, includes CalcSurfaceImpact and adaptive wall temperature | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | +| ** | fully_periodic_3D | ** | Periodic boundary conditions in all three directions | nProcs=10,20,30 | Check whether particles end up outside of the domain | [Link](regressioncheck/WEK_DSMC/fully_periodic_3D/readme.md) | +| ** | Surface_Sticking_Coefficient | ** | Channel flow with a sticking coefficient model | nProcs=5 | Surface sampling | [Link](regressioncheck/WEK_DSMC/Surface_Sticking_Coefficient/readme.md) | +| 5 | Flow_N2_70degCone | [BGK](regressioncheck/WEK_BGKFlow/builds.ini) | 2D axisymmetric 70 degree cone | nProcs=6 | | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | +| ** | MultiSpec_Supersonic_Couette_Ar-He | ** | Supersonic Couette flow with an Ar-He mixture | nProcs=5 | Temperature | [Link](regressioncheck/WEK_DSMC/MultiSpec_Supersonic_Couette_Ar-He/readme.md) | +| ** | MultiSpec_Supersonic_Couette_CO2-N2 | ** | Supersonic Couette flow with a CO2-N2 mixture | nProcs=5 | Temperature | [Link](regressioncheck/WEK_DSMC/MultiSpec_Supersonic_Couette_CO2-N2/readme.md) | +| 6 | Flow_N2_70degCone | [FP](regressioncheck/WEK_FPFlow/builds.ini) | 2D axisymmetric 70 degree cone | nProcs=6 | Surface Sampling, includes CalcSurfaceImpact | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | +| 7 | Flow_N2-N_70degConeHot | [DSMC](regressioncheck/WEK_DSMC/builds.ini) | 2D axisymmetric 70 degree cone (hotter and with N to get some radiation in the next step) | nProcs=6 | Surface Sampling | [Link](regressioncheck/WEK_DSMC/Flow_N2-N_70degConeHot/readme.md) | +| ** | Flow_N2-N_70degConeHot | [Radiation](regressioncheck/WEK_Radiation/builds.ini) | using previously simulated WEK_DSMC/Flow_N2_70degCone results to check radiation tool chain (write out DSMC results, readin those results, radiation solver, radiative transfer, piclas2vtk) | nProcs=6 | Surface heat flux | [Link](regressioncheck/WEK_Radiation/Flow_N2-N_70degConeHot/readme.md) | diff --git a/docs/documentation/userguide/features-and-models/particle-tracking.md b/docs/documentation/userguide/features-and-models/particle-tracking.md index dbe6576bd..f74b808a3 100644 --- a/docs/documentation/userguide/features-and-models/particle-tracking.md +++ b/docs/documentation/userguide/features-and-models/particle-tracking.md @@ -100,14 +100,14 @@ Per default the resting frame of reference is used and no further settings are r The rotating reference frame can be used to represent rotating geometries like e.g. turbine blades, since rotating/changing meshes are currently not supported. The corresponding rotational wall velocity has to be defined for the boundary as well, as shown in Section {ref}`sec:particle-boundary-conditions-reflective-wallvelo`. The distinction between a resting and rotating frame of reference is only important for the particle movement step. Here particles -are not moving on a straight line due to the pseudo forces, i.e. the centripetal force and the Coriolis force. +are not moving on a straight line due to the pseudo forces, i.e. the centrifugal and the Coriolis force. This means that particles follow a circular path towards a stationary boundary that represents a rotating geometry. The usage of the rotating reference frame is enabled by Part-UseRotationalReferenceFrame = T Additionally, the rotational axis (x-, y- or z-axis) and frequency have to be defiend by - Part-RotRefFrame-Axis = 1 ! x=1, y=2, z=3 + Part-RotRefFrame-Axis = 1 ! x=1, y=2, z=3 Part-RotRefFrame-Frequency = -100 ! [Hz, 1/s] The sign of the frequency (+/-) defines the direction of rotation according to the right-hand rule. @@ -125,4 +125,14 @@ boundary surfaces of these regions can only be defined perpendicular to the rota Part-RefFrameRegion2-MAX = 110 This allows to model systems of rotating and stationary geometries (e.g. pumps with stator and rotor blades) within a single simulation. For rotationally symmetric cases, the simulation domain can be reduced using the rotationally perodic boundary condition (as shown in Section {ref}`sec:particle-boundary-conditions-rotBC`). Examples are provided in the regression test directory, e.g. -`regressioncheck/CHE_DSMC/Rotational_Reference_Frame` and `regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Regions`. \ No newline at end of file +`regressioncheck/CHE_DSMC/Rotational_Reference_Frame` and `regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Regions`. + +### Time step subcycling for rotating frame of reference + +In PICLas, an explicit time stepping scheme is used for the DSMC method, with both collision and motion operators altering the particle distribution function within each iteration. This leads to changes in the particle positions, momentum, and energy due to motion and collisions. Operators can be sequentially executed through operator splitting, adjusting the particle positions based on velocities first, followed by collisions within a time step. It is crucial for the time step to resolve collision frequency adequately. External forces (i.e. the centrifugal and the Coriolis force in the case of a rotating reference frame) may require additional consideration for the time step determination, especially when particle acceleration needs to be modeled. To ensure that the existing time step requirement in DSMC, dictated by collisions, remains unaffected, a subcycling algorithm only for the particle motion can be used. This algorithm divides the motion and thus the modeling of external forces into smaller subtimesteps. Consequently, the time step can be chosen based on collision frequency, while the motion can be more finely resolved through subcycling. The usage of the subcycling algorithm is enabled by + + Part-RotRefFrame-UseSubCycling = T + +Additionally, the number of the subcycling steps can be defined by + + Part-RotRefFrame-SubCyclingSteps = 10 ! Default = 10 steps \ No newline at end of file diff --git a/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/DSMC.ini b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/DSMC.ini new file mode 100644 index 000000000..fa26a45d0 --- /dev/null +++ b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/DSMC.ini @@ -0,0 +1,10 @@ +! =============================================================================== ! +! Species1, O2 +! =============================================================================== ! +Part-Species1-InteractionID = 2 +Part-Species1-Tref = 273 +Part-Species1-dref = 4.07E-10 +Part-Species1-omega=0.27 +Part-Species1-CharaTempRot=2.1 +Part-Species1-CharaTempVib=2272.18 +Part-Species1-Ediss_eV=5.17 \ No newline at end of file diff --git a/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/ParticlePosition_Ref.csv b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/ParticlePosition_Ref.csv new file mode 100644 index 000000000..a8c7c3363 --- /dev/null +++ b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/ParticlePosition_Ref.csv @@ -0,0 +1,4 @@ +"001-time","PartNum","PartPosX","PartPosY","PartPosZ","PartVelX","PartVelY","PartVelZ","gamma","Element" +0.0000000000000000E+000,0.1000000000000000E+001,-.2500000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.7071067811865476E+001,0.7071067811865476E+001,0.0000000000000000E+000,0.1000000000000001E+001,0.3000000000000000E+001 +0.2000000000000000E-001,0.1000000000000000E+001,-.4716580473546355E-002,0.1782332060813231E+000,0.0000000000000000E+000,0.7071067811865476E+001,0.7071067811865476E+001,0.0000000000000000E+000,0.1000000000000001E+001,0.1000000000000000E+001 +0.4000000000000000E-001,0.1000000000000000E+001,-.1909047836164587E-001,0.4039505911426143E+000,0.0000000000000000E+000,-.1268564898553164E+002,0.7071067811865476E+001,0.0000000000000000E+000,0.1000000000000001E+001,0.1000000000000000E+001 diff --git a/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/analyze.ini b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/analyze.ini new file mode 100644 index 000000000..4ddd7ae19 --- /dev/null +++ b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/analyze.ini @@ -0,0 +1,6 @@ +! compare columns +compare_column_file = ParticlePosition.csv +compare_column_reference_file = ParticlePosition_Ref.csv +compare_column_index = 0,2,3 ! column index for comparison (starts at 0) +compare_column_tolerance_value = 0.0001 ! tolerance (depends on machine accuracy and MPI) +compare_column_tolerance_type = relative ! absolute or relative tolerance diff --git a/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/command_line.ini b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/command_line.ini new file mode 100644 index 000000000..2df6de114 --- /dev/null +++ b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/command_line.ini @@ -0,0 +1,2 @@ +MPI=1 +cmd_suffix=DSMC.ini diff --git a/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/externals.ini b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/externals.ini new file mode 100644 index 000000000..563c06535 --- /dev/null +++ b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/externals.ini @@ -0,0 +1,7 @@ +! --- Externals Tool Reggie +MPI = 1 +externalbinary = ./hopr/build/bin/hopr +externaldirectory = hopr.ini +externalruntime = pre + +nocrosscombination:MPI,externalbinary,externaldirectory,externalruntime diff --git a/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/hopr.ini b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/hopr.ini new file mode 100644 index 000000000..4a27dd61b --- /dev/null +++ b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/hopr.ini @@ -0,0 +1,55 @@ +ProjectName = cube +Debugvisu = T +DebugVisuLevel=2 +NVisu =1 +Mode =1 + +DEFVAR = (REAL): minus_x = -0.5 +DEFVAR = (REAL): plus_x = 0.5 +DEFVAR = (REAL): top_x_incline = 0.0 + +DEFVAR = (REAL): minus_y = -0.5 +DEFVAR = (REAL): plus_y = 0.5 + +DEFVAR = (REAL): minus_z = -0.1 +DEFVAR = (REAL): plus_z = 0.1 + +nZones = 3 + +Corner =(/minus_x,minus_y,minus_z ,, 0.0,minus_y,minus_z ,, 0.0,0.0,minus_z ,, minus_x,0.0,minus_z ,, minus_x,minus_y,plus_z ,, 0,minus_y,plus_z ,, 0.0,0.0,plus_z ,, minus_x,0.0,plus_z /) +nElems =(/1,1,1/) +BCIndex =(/6 ,4 ,0 ,0 ,2 ,5/) +elemtype =108 + +Corner =(/minus_x,0.0,minus_z ,, 0.0,0.0,minus_z ,,top_x_incline,plus_y,minus_z ,, minus_x,plus_y,minus_z ,, minus_x,0.0,plus_z ,,0.0,0.0,plus_z ,, top_x_incline,plus_y,plus_z ,, minus_x,plus_y,plus_z /) +nElems =(/1,1,1/) +BCIndex =(/6 ,0 ,7 ,3 ,2 ,5/) +elemtype =108 + +Corner =(/0.0,minus_y,minus_z ,, plus_x,minus_y,minus_z ,, plus_x,0.0,minus_z ,, 0.0,0.0,minus_z ,, 0.0,minus_y,plus_z ,, plus_x,minus_y,plus_z ,, plus_x,0.0,plus_z ,, 0.0,0.0,plus_z /) +nElems =(/1,1,1/) +BCIndex =(/6 ,4 ,1 ,7 ,0 ,5/) +elemtype =108 + + + + +nUserDefinedBoundaries=7 +BoundaryName=BC_Xplus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Xminus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Yplus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Yminus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Zplus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Zminus +BoundaryType=(/4,0,0,0/) +BoundaryName=BC_Wall +BoundaryType=(/4,0,0,0/) + +postscalemesh=true +!meshscale=1e-5 +jacobiantolerance=1e-20 diff --git a/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/parameter.ini b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/parameter.ini new file mode 100644 index 000000000..608334883 --- /dev/null +++ b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/parameter.ini @@ -0,0 +1,99 @@ +! =============================================================================== ! +! EQUATION (linearscalaradvection) +! =============================================================================== ! +IniExactFunc = 0 +! =============================================================================== ! +! DISCRETIZATION +! =============================================================================== ! +N = 1 ! Polynomial degree +NAnalyze = 1 ! Number of analyze points +! =============================================================================== ! +! MESH +! =============================================================================== ! +MeshFile = cube_mesh.h5 +useCurveds = T +! =============================================================================== ! +! OUTPUT / VISUALIZATION +! =============================================================================== ! +ProjectName = RotRefFrame +Logging = F +WriteErrorFiles = F +IterDisplayStep = 1 +DoCalcErrorNorms = T +Part-AnalyzeStep = 10 +! =============================================================================== ! +! CALCULATION +! =============================================================================== ! +ManualTimeStep = 2E-3 +tend = 4E-2 ! End time +Analyze_dt = 1.0 ! Timestep of analyze outputs +CFLscale = 0.5 + +Part-RotRefFrame-UseSubCycling = T +Part-RotRefFrame-SubCyclingSteps = 200 +! =============================================================================== ! +! BOUNDARIES +! =============================================================================== ! +Part-nBounds=7 +Part-Boundary1-SourceName=BC_Xplus +Part-Boundary1-Condition=open +Part-Boundary2-SourceName=BC_Xminus +Part-Boundary2-Condition=open +Part-Boundary3-SourceName=BC_Yplus +Part-Boundary3-Condition=open +Part-Boundary4-SourceName=BC_Yminus +Part-Boundary4-Condition=open +Part-Boundary5-SourceName=BC_Zplus +Part-Boundary5-Condition=open +Part-Boundary6-SourceName=BC_Zminus +Part-Boundary6-Condition=open +Part-Boundary7-SourceName=BC_Wall +Part-Boundary7-Condition=reflective +Part-Boundary7-WallTemp = 1.0 +Part-Boundary7-MomentumACC = 0. +Part-Boundary7-TransACC = 0. +Part-Boundary7-VibACC = 0. +Part-Boundary7-RotACC = 0. +Part-Boundary7-ElecACC = 0. +Part-Boundary7-RotVelo = T +Part-Boundary7-RotAxis = 3 +Part-Boundary7-RotFreq = 5 +Part-FIBGMdeltas=(/0.5,0.5,0.2/) +! =============================================================================== ! +! Tracking +! =============================================================================== ! +TrackingMethod = triatracking +! =============================================================================== ! +! PARTICLES +! =============================================================================== ! +Part-maxParticleNumber=10 +Part-nSpecies=1 +Part-Species1-ChargeIC=0. +Part-Species1-MassIC=1. +Part-Species1-MacroParticleFactor=1E0 + +Part-Species1-nInits = 1 +Part-Species1-Init1-SpaceIC=point +Part-Species1-Init1-ParticleNumber=1 +Part-Species1-Init1-BasePointIC=(/-0.25,0.0,0./) +Part-Species1-Init1-NormalIC=(/0.,0.,0.0/) +Part-Species1-Init1-velocityDistribution=constant +Part-Species1-Init1-VeloIC=10. +Part-Species1-Init1-VeloVecIC=(/1.0,1.0,0.0/) +! =============================================================================== ! +! Analysis +! =============================================================================== ! +Part-TrackPosition = T +! =============================================================================== ! +! DSMC +! =============================================================================== ! +UseDSMC=true +Particles-DSMC-CollisMode=0 ! Collisionless flow +Part-NumberOfRandomSeeds =2 +Particles-RandomSeed1= 1 +Particles-RandomSeed2= 2 +Particles-HaloEpsVelo=2 + +Part-UseRotationalReferenceFrame = T +Part-RotRefFrame-Axis = 3 +Part-RotRefFrame-Frequency = 5 diff --git a/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/readme.md b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/readme.md new file mode 100644 index 000000000..483e9da09 --- /dev/null +++ b/regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/readme.md @@ -0,0 +1,5 @@ +# Rotational frame of reference: Subcycling +* Frame of reference is rotating with 5 revolutions per second +* A single particle is initially placed at coordinates (x,y) = (-0.25,-0.0) with a velocity vector of (10,10,0). It is anticipated to traverse a circular-like trajectory within the rotating frame of reference due to fictitious forces. +* A wall, rotating synchronously with the frame of reference, induces a specular reflection upon collision. +* A relatively large time step is employed (Delta_t=2E-3). Reference positions are calculated using a smaller time step (1E-5). Utilization of a subcycling with 200 substeps ensures identical positions in the regression test as in the reference case. Disabling subcycling results in failure of the regression test. diff --git a/regressioncheck/NIG_code_analyze/Rotational_Reference_Frame_Wall_Specular/analyze.ini b/regressioncheck/NIG_code_analyze/Rotational_Reference_Frame_Wall_Specular/analyze.ini index 62b18d068..57d60d773 100644 --- a/regressioncheck/NIG_code_analyze/Rotational_Reference_Frame_Wall_Specular/analyze.ini +++ b/regressioncheck/NIG_code_analyze/Rotational_Reference_Frame_Wall_Specular/analyze.ini @@ -4,7 +4,7 @@ L2_error_name = L2_Part ! t-convergence test analyze_Convtest_t_rate = 0.4 ! allow 4 out of 7 unsuccesful tests -analyze_convtest_t_tolerance = 0.1 +analyze_convtest_t_tolerance = 0.76 analyze_convtest_t_order = 2 analyze_convtest_t_error_name = L2_Part diff --git a/regressioncheck/NIG_code_analyze/Rotational_Reference_Frame_Wall_Specular/parameter.ini b/regressioncheck/NIG_code_analyze/Rotational_Reference_Frame_Wall_Specular/parameter.ini index 0ec8a3b69..739b3ad56 100644 --- a/regressioncheck/NIG_code_analyze/Rotational_Reference_Frame_Wall_Specular/parameter.ini +++ b/regressioncheck/NIG_code_analyze/Rotational_Reference_Frame_Wall_Specular/parameter.ini @@ -25,7 +25,7 @@ Part-AnalyzeStep = 10 ! CALCULATION ! =============================================================================== ! ManualTimeStep = 1E-3,5E-4,1E-4,5E-5,1E-5,5E-6 -tend = 0.38 ! End time +tend = 0.25 ! End time Analyze_dt = 1.0 ! Timestep of analyze outputs CFLscale = 0.5 ! =============================================================================== ! @@ -75,9 +75,9 @@ Part-Species1-Init1-ParticleNumber=1 Part-Species1-Init1-BasePointIC=(/-0.25,0.0,0./) Part-Species1-Init1-NormalIC=(/1.,0.,0.0/) Part-Species1-Init1-velocityDistribution=constant -Part-Species1-Init1-VeloIC=0.0 +Part-Species1-Init1-VeloIC=1.0 -Part-Species1-Init1-VeloVecIC=(/0.0,1.0,0.0/) +Part-Species1-Init1-VeloVecIC=(/1.0,1.0,0.0/) ! =============================================================================== ! ! Analysis ! =============================================================================== ! diff --git a/src/particles/analyze/particle_analyze_code.f90 b/src/particles/analyze/particle_analyze_code.f90 index 761abe2e4..c1c3d1873 100644 --- a/src/particles/analyze/particle_analyze_code.f90 +++ b/src/particles/analyze/particle_analyze_code.f90 @@ -347,12 +347,13 @@ SUBROUTINE CalcAnalyticalParticleState(t,PartStateAnalytic,alpha_out,theta_out) IF(TimeReset.GT.0.0) THEN r_0Vec = r_WallVec v_0Vec = v_WallVec + v_0Vec = v_0Vec - CROSS(RotRefFrameOmega(1:3),r_0Vec) ! Transform velocity into RotRefFrame New_t = t - TimeReset ELSE r_0Vec = Species(iSpec)%Init(1)%BasePointIC(1:3) v_0Vec = Species(iSpec)%Init(1)%VeloVecIC(1:3) * Species(iSpec)%Init(1)%VeloIC New_t = t - v_0Vec = v_0Vec - CROSS(RotRefFrameOmega(1:3),r_0Vec) + v_0Vec = v_0Vec - CROSS(RotRefFrameOmega(1:3),r_0Vec) ! Transform velocity into RotRefFrame END IF r_tVec = r_0Vec + v_0Vec * New_t omega = 2.*PI*RotRefFrameFreq @@ -407,13 +408,14 @@ SUBROUTINE CalcAnalyticalParticleState(t,PartStateAnalytic,alpha_out,theta_out) - ( omega * (COS(omega * New_t) - omega * New_t * SIN(omega * New_t) ) ) & * ( TempArrayCross3(3) - 1/omega * TempArrayCross6(3) ) -! IF((PartStateAnalytic(1).GE.0.0).AND.(TimeReset.LE.0.0)) THEN - IF(ABS(PartStateAnalytic(1)).LT.1E-8) THEN + IF(PartStateAnalytic(1).GT.0.0) THEN TimeReset = t + PartStateAnalytic(1) = 0.0 ! set particle on wall r_WallVec = PartStateAnalytic(1:3) - v_WallVec(1) = -PartStateAnalytic(4) - v_WallVec(2) = PartStateAnalytic(5) - v_WallVec(3) = PartStateAnalytic(6) + v_0Vec = Species(iSpec)%Init(1)%VeloVecIC(1:3) * Species(iSpec)%Init(1)%VeloIC + v_WallVec = v_0Vec + v_WallVec(1) = - v_WallVec(1) ! mirror velocity (not equal push velocity or PartStateAnalytic(4:6)) + v_WallVec = v_WallVec + CROSS(RotRefFrameOmega(1:3),r_WallVec) ! add wall velocity END IF ! PartStateAnalytic(4:6) = 0.0 END SELECT diff --git a/src/particles/boundary/particle_boundary_condition.f90 b/src/particles/boundary/particle_boundary_condition.f90 index 45f4aaacb..6e20886cb 100644 --- a/src/particles/boundary/particle_boundary_condition.f90 +++ b/src/particles/boundary/particle_boundary_condition.f90 @@ -466,6 +466,8 @@ SUBROUTINE RotPeriodicInterPlaneBoundary(PartID,SideID,ElemID,IsInterPlanePart) USE MOD_Globals_Vars ,ONLY: PI USE MOD_Mesh_Tools ,ONLY: GetCNElemID USE MOD_Particle_Vars ,ONLY: PartState,LastPartPos,Species,PartSpecies +USE MOD_Particle_Vars ,ONLY: RotRefSubTimeStep, NewPosSubCycling, GlobalElemIDSubCycling, LastPartPosSubCycling,nSubCyclingSteps +USE MOD_Particle_Vars ,ONLY: InRotRefFrameSubCycling, PartVeloRotRefSubCycling, LastVeloRotRefSubCycling USE MOD_Particle_Mesh_Vars ,ONLY: ElemInfo_Shared, SideInfo_Shared, ElemSideNodeID_Shared, NodeCoords_Shared USE MOD_Particle_Boundary_Vars ,ONLY: PartBound, InterPlaneSideMapping USE MOD_TImeDisc_Vars ,ONLY: dt,RKdtFrac @@ -477,11 +479,12 @@ SUBROUTINE RotPeriodicInterPlaneBoundary(PartID,SideID,ElemID,IsInterPlanePart) USE MOD_DSMC_Vars ,ONLY: DSMC, AmbipolElecVelo USE MOD_part_operations ,ONLY: CreateParticle, RemoveParticle USE MOD_DSMC_Vars ,ONLY: CollisMode, useDSMC, PartStateIntEn -USE MOD_Particle_Vars ,ONLY: usevMPF,PartMPF,PDM,InterPlanePartNumber, InterPlanePartIndx -USE MOD_Particle_Vars ,ONLY: UseRotRefFrame, RotRefFrameOmega, PartVeloRotRef +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 USE MOD_Particle_Tracking_Vars ,ONLY: PartOut,MPIRankOut #endif /*CODE_ANALYZE*/ @@ -496,11 +499,11 @@ SUBROUTINE RotPeriodicInterPlaneBoundary(PartID,SideID,ElemID,IsInterPlanePart) LOGICAL,INTENT(IN),OPTIONAL :: IsInterPlanePart !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: iSide, InterSideID, NumInterPlaneSides,NewElemID, neighPartBound +INTEGER :: iSide, InterSideID, NumInterPlaneSides,NewElemID, neighPartBound, SpecID INTEGER :: BCType, nLocSides, newSideID, iLocSide, locSideID, TriNum, NodeNum REAL :: dtVar LOGICAL :: ParticleFound, DoCreateParticles, ThroughSide -REAL :: Velo_old(1:3), Velo_oldAmbi(1:3) +REAL :: Velo_old(1:3), Velo_oldAmbi(1:3), NewPos(1:3), NewVelo(1:3) REAL :: POI(1:3), POI_rotated(1:3), LastPartPos_rotated(1:3), PartState_rotated(1:3) REAL :: RanNum, RadiusPOI, RanAlpha, RotAlpha, RotDir REAL :: RadiusInterPlane(1:2) @@ -549,6 +552,17 @@ SUBROUTINE RotPeriodicInterPlaneBoundary(PartID,SideID,ElemID,IsInterPlanePart) ! (1.c.II) Create inter particles DO iNewPart = 1,NewPartNumber InterPlanePartNumber = InterPlanePartNumber + 1 + ! In case of sub cycling step particle information before sub cycling must be used => interplane particle can act like origin particle + IF(RotRefSubTimeStep) THEN + NewPos(1:3) = NewPosSubCycling(1:3) + NewElemID = GlobalElemIDSubCycling + ELSE + NewPos(1:3) = PartState(1:3,PartID) + NewElemID = ElemID + END IF + ! Store the values in a separate variable to avoid memory leaks, as these arrays might be resized during CreateParticle + SpecID = PartSpecies(PartID) + NewVelo(1:3) = PartState(4:6,PartID) IF (useDSMC.AND.(CollisMode.GT.1)) THEN VibEnergy = PartStateIntEn(1,PartID) RotEnergy = PartStateIntEn(2,PartID) @@ -562,33 +576,40 @@ SUBROUTINE RotPeriodicInterPlaneBoundary(PartID,SideID,ElemID,IsInterPlanePart) RotEnergy = 0.0 ElecEnergy = 0.0 END IF - ! For creating inter particles: - ! - LastPartPos(1:3,NewPartID) must be redefined as long as LastPartPos is set to PartState in CreateParticle routine - ! - ParticleInside for InterParticles must be .FALSE. in order to avoid error looping over the original PDM%ParticleVecLength - ! in ParticleTriaTracking() routine. The inside flag is set to .TRUE. - ! when we loop over all inter particle in SingleParticleTriaTracking routine - IF (usevMPF) THEN - CALL CreateParticle( PartSpecies(PartID), PartState(1:3,PartID), ElemID, PartState(4:6,PartID) & - , RotEnergy=RotEnergy,VibEnergy=VibEnergy,ElecEnergy=ElecEnergy & - , NewPartID=NewPartID, NewMPF=PartMPF(PartID) ) - LastPartPos(1:3,NewPartID) = LastPartPos(1:3,PartID) - PDM%ParticleInside(NewPartID) = .FALSE. - InterPlanePartIndx(InterPlanePartNumber) = NewPartID + ! === Creating inter particles: + CALL CreateParticle(SpecID,NewPos(1:3),NewElemID,NewVelo(1:3),RotEnergy=RotEnergy,VibEnergy=VibEnergy,ElecEnergy=ElecEnergy, & + OldPartID=PartID,NewPartID=NewPartID ) + ! LastPartPos(1:3,NewPartID) must be redefined as long as LastPartPos is set to PartState in CreateParticle routine + IF(RotRefSubTimeStep) THEN + LastPartPos(1:3,NewPartID) = LastPartPosSubCycling(1:3) ELSE - CALL CreateParticle( PartSpecies(PartID), PartState(1:3,PartID), ElemID,PartState(4:6,PartID) & - , RotEnergy=RotEnergy,VibEnergy=VibEnergy,ElecEnergy=ElecEnergy & - , NewPartID=NewPartID ) LastPartPos(1:3,NewPartID) = LastPartPos(1:3,PartID) - PDM%ParticleInside(NewPartID) = .FALSE. - InterPlanePartIndx(InterPlanePartNumber) = NewPartID END IF + ! ParticleInside for InterParticles must be .FALSE. in order to avoid error looping over the original PDM%ParticleVecLength + ! in ParticleTriaTracking() routine. The inside flag is set to .TRUE. when we loop over all inter particle in SingleParticleTriaTracking routine + PDM%ParticleInside(NewPartID) = .FALSE. + ! Storing the new particle index + InterPlanePartIndx(InterPlanePartNumber) = NewPartID ! Treatment for the rotational frame of reference: stored here, will be rotated together with the regular velocity later IF(UseRotRefFrame) THEN - PDM%InRotRefFrame(NewPartID) = PDM%InRotRefFrame(PartID) - IF(PDM%InRotRefFrame(NewPartID)) THEN - PartVeloRotRef(1:3,NewPartID) = PartVeloRotRef(1:3,PartID) + IF(RotRefSubTimeStep) THEN + PDM%InRotRefFrame(NewPartID) = InRotRefFrameSubCycling + IF(PDM%InRotRefFrame(NewPartID)) THEN + PartVeloRotRef(1:3,NewPartID) = PartVeloRotRefSubCycling(1:3) + LastPartVeloRotRef(1:3,NewPartID) = LastVeloRotRefSubCycling(1:3) + ELSE + PartVeloRotRef(1:3,NewPartID) = 0. + LastPartVeloRotRef(1:3,NewPartID) = 0. + END IF ELSE - PartVeloRotRef(1:3,NewPartID) = 0. + PDM%InRotRefFrame(NewPartID) = PDM%InRotRefFrame(PartID) + IF(PDM%InRotRefFrame(NewPartID)) THEN + PartVeloRotRef(1:3,NewPartID) = PartVeloRotRef(1:3,PartID) + LastPartVeloRotRef(1:3,NewPartID) = LastPartVeloRotRef(1:3,PartID) + ELSE + PartVeloRotRef(1:3,NewPartID) = 0. + LastPartVeloRotRef(1:3,NewPartID) = 0. + END IF END IF END IF END DO @@ -606,6 +627,8 @@ SUBROUTINE RotPeriodicInterPlaneBoundary(PartID,SideID,ElemID,IsInterPlanePart) ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific) dtVar = dtVar * Species(PartSpecies(PartID))%TimeStepFactor +IF(RotRefSubTimeStep) dtVar = dtVar / REAL(nSubCyclingSteps) + ! (2) Calculate the POI and a new random POI on corresponding inter plane using a random angle within the periodic segment POI(1:3) = LastPartPos(1:3,PartID) + TrackInfo%PartTrajectory(1:3)*TrackInfo%alpha IF (DSMC%DoAmbipolarDiff) THEN @@ -681,29 +704,19 @@ SUBROUTINE RotPeriodicInterPlaneBoundary(PartID,SideID,ElemID,IsInterPlanePart) ! (5) Treatment of velocity in rotational frame of reference IF(UseRotRefFrame) THEN - ! Setting the PartState to the POI to determine whether the particle moved into a RotRefFrame (is later overwritten anyway) - PartState(1:3,PartID) = POI_rotated(1:3) - ! Check is repeated in the FUNCTION InRotRefFrameCheck at the current PartState(1:3) - IF(InRotRefFrameCheck(PartID)) THEN - ! Particle moved into a RotRefFrame - IF(PDM%InRotRefFrame(PartID)) THEN - ! Particle comes from a RotRefFrame: rotate the old PartVeloRotRef - Velo_old(1:3) = PartVeloRotRef(1:3,PartID) - PartVeloRotRef(1:3,PartID) = RotateVectorAroundAxis(Velo_old(1:3),PartBound%RotPeriodicAxis,RotAlpha) - ELSE - ! Particle comes from an inertial frame: initialize the new PartVeloRotRef - PartVeloRotRef(1:3,PartID) = PartState(4:6,PartID) - CROSS(RotRefFrameOmega(1:3),PartState(1:3,PartID)) - END IF - ! Calculate the acceleration - PartVeloRotRef(1:3,PartID) = PartVeloRotRef(1:3,PartID) + CalcPartRHSRotRefFrame(PartState(1:3,PartID),PartVeloRotRef(1:3,PartID)) & - * dtVar * (1.0 - TrackInfo%alpha/TrackInfo%lengthPartTrajectory) - ELSE + ! Check from which frame are we coming from and adapt logical accordingly (assuming that an interplane is ALWAYS between two different reference frames) + IF(PDM%InRotRefFrame(PartID)) THEN + ! Particle comes from the rotational frame and CONSEQUENTLY moves to the inertial frame PartVeloRotRef(1:3,PartID) = 0. + PDM%InRotRefFrame(PartID) = .FALSE. + ELSE + ! Particle comes from an inertial frame: initialize the new PartVeloRotRef at the target location + PartVeloRotRef(1:3,PartID) = PartState(4:6,PartID) - CROSS(RotRefFrameOmega(1:3),PartState_rotated(1:3)) + PDM%InRotRefFrame(PartID) = .TRUE. END IF - PDM%InRotRefFrame(PartID) = InRotRefFrameCheck(PartID) END IF -! (7) Track the particle, moving inside the domain through the interplane BC +! (6) Track the particle, moving inside the domain through the interplane BC ParticleFound = .FALSE. PartState(1:3,PartID) = PartState_rotated(1:3) LastPartPos(1:3,PartID) = LastPartPos_rotated(1:3) diff --git a/src/particles/boundary/particle_boundary_init.f90 b/src/particles/boundary/particle_boundary_init.f90 index afce0706b..d7e1f5c0f 100644 --- a/src/particles/boundary/particle_boundary_init.f90 +++ b/src/particles/boundary/particle_boundary_init.f90 @@ -2243,17 +2243,13 @@ SUBROUTINE FinalizeParticleBoundary() IF(nRotPeriodicSides .GT.0) CALL UNLOCK_AND_FREE(NumRotPeriodicNeigh_Shared_Win) IF(MaxNumRotPeriodicNeigh.GT.0) CALL UNLOCK_AND_FREE(RotPeriodicSideMapping_Shared_Win) ADEALLOCATE(SurfSide2RotPeriodicSide_Shared) - ADEALLOCATE(SurfSide2RotPeriodicSide) ADEALLOCATE(NumRotPeriodicNeigh_Shared) - ADEALLOCATE(NumRotPeriodicNeigh) ADEALLOCATE(RotPeriodicSideMapping_Shared) - ADEALLOCATE(RotPeriodicSideMapping) - ADEALLOCATE(InterPlaneSideMapping) -#else - SDEALLOCATE(SurfSide2RotPeriodicSide) - SDEALLOCATE(NumRotPeriodicNeigh) - SDEALLOCATE(RotPeriodicSideMapping) #endif + ADEALLOCATE(SurfSide2RotPeriodicSide) + ADEALLOCATE(NumRotPeriodicNeigh) + ADEALLOCATE(RotPeriodicSideMapping) + SDEALLOCATE(InterPlaneSideMapping) END IF ! PartBound%UseRotPeriodicBC ! Adaptive wall temperature (e.g. calculate from sampled heat flux) diff --git a/src/particles/boundary/particle_boundary_vars.f90 b/src/particles/boundary/particle_boundary_vars.f90 index 1ae21844e..96323a0f1 100644 --- a/src/particles/boundary/particle_boundary_vars.f90 +++ b/src/particles/boundary/particle_boundary_vars.f90 @@ -132,7 +132,7 @@ MODULE MOD_Particle_Boundary_Vars INTEGER,ALLOCPOINT,DIMENSION(:) :: SurfSide2RotPeriodicSide ! Mapping between surf side and periodic sides. ! ==================================================================== ! Intermediate plane for multi rotational periodic sides -INTEGER,ALLOCPOINT,DIMENSION(:,:) :: InterPlaneSideMapping ! Mapping between inter plane BC_ID and SideID. +INTEGER,ALLOCATABLE :: InterPlaneSideMapping(:,:)! Mapping between inter plane BC_ID and SideID. ! ==================================================================== #if USE_MPI INTEGER,POINTER,DIMENSION(:) :: SurfSide2RotPeriodicSide_Shared diff --git a/src/particles/particle_init.f90 b/src/particles/particle_init.f90 index 5d5d8e086..59b42a512 100644 --- a/src/particles/particle_init.f90 +++ b/src/particles/particle_init.f90 @@ -212,7 +212,9 @@ SUBROUTINE DefineParametersParticles() , 'Set [T] to activate sampling of electronic energy excitation (currently only available for ElectronicModel = 3)', '.FALSE.') ! === Rotational frame of reference -CALL prms%CreateLogicalOption( 'Part-UseRotationalReferenceFrame', 'Activate rotational frame of reference', '.FALSE.') +CALL prms%CreateLogicalOption( 'Part-UseRotationalReferenceFrame', 'Activate the rotational frame of reference', '.FALSE.') +CALL prms%CreateLogicalOption( 'Part-RotRefFrame-UseSubCycling', 'Activate subcycling in the rotational frame of reference', '.FALSE.') +CALL prms%CreateIntOption( 'Part-RotRefFrame-SubCyclingSteps','Number of subcyling steps)','10') CALL prms%CreateIntOption( 'Part-RotRefFrame-Axis','Axis of rotational frame of reference (x=1, y=2, z=3)') CALL prms%CreateRealOption( 'Part-RotRefFrame-Frequency','Frequency of rotational frame of reference [1/s], sign according '//& 'to right-hand rule, e.g. positive: counter-clockwise, negative: clockwise') @@ -1472,6 +1474,7 @@ SUBROUTINE FinalizeParticles() SDEALLOCATE(PartState) SDEALLOCATE(LastPartPos) SDEALLOCATE(PartVeloRotRef) +SDEALLOCATE(LastPartVeloRotRef) SDEALLOCATE(PartSpecies) SDEALLOCATE(Pt) SDEALLOCATE(PDM%ParticleInside) @@ -1500,7 +1503,7 @@ SUBROUTINE FinalizeParticles() SDEALLOCATE(PEM%pNext) SDEALLOCATE(seeds) SDEALLOCATE(PartPosLandmark) -SDEALLOCATE(RotRefFramRegion) +SDEALLOCATE(RotRefFrameRegion) SDEALLOCATE(VirtMergedCells) SDEALLOCATE(PartDataVarNames) #if USE_MPI @@ -1727,6 +1730,8 @@ SUBROUTINE InitializeVariablesRotationalRefFrame() !=================================================================================================================================== UseRotRefFrame = GETLOGICAL('Part-UseRotationalReferenceFrame') +UseRotSubCycling = GETLOGICAL('Part-RotRefFrame-UseSubCycling') +nSubCyclingSteps = GETINT('Part-RotRefFrame-SubCyclingSteps') IF(UseRotRefFrame) THEN ! Abort for other timedisc except DSMC/BGK @@ -1735,6 +1740,8 @@ SUBROUTINE InitializeVariablesRotationalRefFrame() #endif ALLOCATE(PartVeloRotRef(1:3,1:PDM%maxParticleNumber)) PartVeloRotRef = 0.0 + ALLOCATE(LastPartVeloRotRef(1:3,1:PDM%maxParticleNumber)) + LastPartVeloRotRef = 0.0 RotRefFrameAxis = GETINT('Part-RotRefFrame-Axis') RotRefFrameFreq = GETREAL('Part-RotRefFrame-Frequency') omegaTemp = 2.*PI*RotRefFrameFreq @@ -1749,13 +1756,13 @@ SUBROUTINE InitializeVariablesRotationalRefFrame() CALL abort(__STAMP__,'ERROR Rotational Reference Frame: Axis must be between 1 and 3. Selected axis: ',IntInfoOpt=RotRefFrameAxis) END SELECT nRefFrameRegions = GETINT('Part-nRefFrameRegions') - ALLOCATE(RotRefFramRegion(1:2,1:nRefFrameRegions)) + ALLOCATE(RotRefFrameRegion(1:2,1:nRefFrameRegions)) IF(nRefFrameRegions.GT.0)THEN DO iRegion=1, nRefFrameRegions WRITE(UNIT=hilf,FMT='(I0)') iRegion - RotRefFramRegion(1,iRegion)= GETREAL('Part-RefFrameRegion'//TRIM(hilf)//'-MIN') - RotRefFramRegion(2,iRegion)= GETREAL('Part-RefFrameRegion'//TRIM(hilf)//'-MAX') - IF(RotRefFramRegion(1,iRegion).GE.RotRefFramRegion(2,iRegion)) THEN + RotRefFrameRegion(1,iRegion)= GETREAL('Part-RefFrameRegion'//TRIM(hilf)//'-MIN') + RotRefFrameRegion(2,iRegion)= GETREAL('Part-RefFrameRegion'//TRIM(hilf)//'-MAX') + IF(RotRefFrameRegion(1,iRegion).GE.RotRefFrameRegion(2,iRegion)) THEN CALL abort(__STAMP__,'ERROR Rotational Reference Frame: MIN > MAX in definition of region ',IntInfoOpt=iRegion) END IF END DO diff --git a/src/particles/particle_operations.f90 b/src/particles/particle_operations.f90 index eeb51a805..38a66cbce 100644 --- a/src/particles/particle_operations.f90 +++ b/src/particles/particle_operations.f90 @@ -31,9 +31,10 @@ MODULE MOD_part_operations CONTAINS -SUBROUTINE CreateParticle(SpecID,Pos,GlobElemID,Velocity,RotEnergy,VibEnergy,ElecEnergy,NewPartID,NewMPF) +SUBROUTINE CreateParticle(SpecID,Pos,GlobElemID,Velocity,RotEnergy,VibEnergy,ElecEnergy,OldPartID,NewPartID,NewMPF,NewTimestep) !=================================================================================================================================== !> creates a single particle at correct array position and assign properties +!> OldPartID can be supplied to get the MPF and Timestep, however, NewMPF and NewTimestep have priority !=================================================================================================================================== ! MODULES USE MOD_Globals @@ -57,8 +58,10 @@ SUBROUTINE CreateParticle(SpecID,Pos,GlobElemID,Velocity,RotEnergy,VibEnergy,Ele REAL, INTENT(IN) :: RotEnergy !< Rotational energy REAL, INTENT(IN) :: VibEnergy !< Vibrational energy REAL, INTENT(IN) :: ElecEnergy !< Electronic energy +INTEGER, INTENT(IN),OPTIONAL :: OldPartID !< ID of old particle (if available) INTEGER, INTENT(OUT),OPTIONAL :: NewPartID !< ID of newly created particle REAL, INTENT(IN),OPTIONAL :: NewMPF !< MPF of newly created particle +REAL, INTENT(IN),OPTIONAL :: NewTimestep !< Timestep of the newly created particle !----------------------------------------------------------------------------------------------------------------------------------! ! LOCAL VARIABLES INTEGER :: newParticleID @@ -96,7 +99,13 @@ SUBROUTINE CreateParticle(SpecID,Pos,GlobElemID,Velocity,RotEnergy,VibEnergy,Ele ! Set particle time step and weight (if required) IF (UseVarTimeStep) THEN - PartTimeStep(newParticleID) = GetParticleTimeStep(PartState(1,newParticleID),PartState(2,newParticleID),PEM%LocalElemID(newParticleID)) + IF(PRESENT(NewTimestep)) THEN + PartTimeStep(newParticleID) = NewTimestep + ELSE IF(PRESENT(OldPartID)) THEN + PartTimeStep(newParticleID) = PartTimeStep(OldPartID) + ELSE + PartTimeStep(newParticleID) = GetParticleTimeStep(PartState(1,newParticleID),PartState(2,newParticleID),PEM%LocalElemID(newParticleID)) + END IF END IF ! Set new particle MPF @@ -104,6 +113,9 @@ SUBROUTINE CreateParticle(SpecID,Pos,GlobElemID,Velocity,RotEnergy,VibEnergy,Ele IF(PRESENT(NewMPF))THEN ! MPF is already defined via input PartMPF(newParticleID) = NewMPF + ELSE IF(PRESENT(OldPartID)) THEN + ! MPF is available from the "original" particle + PartMPF(newParticleID) = PartMPF(OldPartID) ELSE ! Check if vMPF (and radial weighting is used) to determine the MPF of the new particle IF (RadialWeighting%DoRadialWeighting) THEN diff --git a/src/particles/particle_rhs.f90 b/src/particles/particle_rhs.f90 index 42edd3df0..5af5ab2b2 100644 --- a/src/particles/particle_rhs.f90 +++ b/src/particles/particle_rhs.f90 @@ -42,6 +42,7 @@ MODULE MOD_part_RHS PUBLIC :: PartRHS PUBLIC :: CalcPartRHSSingleParticle PUBLIC :: CalcPartRHSRotRefFrame +PUBLIC :: CalcPartPosInRotRef !---------------------------------------------------------------------------------------------------------------------------------- ABSTRACT INTERFACE @@ -794,4 +795,44 @@ SUBROUTINE PartVeloToImp(VeloToImp,doParticle_In) END SUBROUTINE PartVeloToImp + +SUBROUTINE CalcPartPosInRotRef(iPart, RotTimestep) +!=================================================================================================================================== +!> Particle push in rotational frame of reference using the midpoint method +!=================================================================================================================================== +! MODULES +USE MOD_Particle_Vars ,ONLY: PartState, PDM, PartVeloRotRef +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER, INTENT(IN) :: iPart +REAL, INTENT(IN) :: RotTimestep +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +REAL :: Pt_local(1:3), Pt_local_old(1:3), VeloRotRef_half(1:3), PartState_half(1:3) +!=================================================================================================================================== +IF(PDM%InRotRefFrame(iPart)) THEN + ! Midpoint method + ! calculate the acceleration (force / mass) at the current time step + Pt_local_old(1:3) = CalcPartRHSRotRefFrame(PartState(1:3,iPart), PartVeloRotRef(1:3,iPart)) + ! estimate the midpoint velocity in the rotational frame + VeloRotRef_half(1:3) = PartVeloRotRef(1:3,iPart) + 0.5*Pt_local_old(1:3)*RotTimestep + ! estimate the midpoint position + PartState_half(1:3) = PartState(1:3,iPart) + 0.5*PartVeloRotRef(1:3,iPart)*RotTimestep + ! calculate the acceleration (force / mass) at the midpoint + Pt_local(1:3) = CalcPartRHSRotRefFrame(PartState_half(1:3), VeloRotRef_half(1:3)) + ! update the position using the midpoint velocity in the rotational frame + PartState(1:3,iPart) = PartState(1:3,iPart) + VeloRotRef_half(1:3)*RotTimestep + ! update the velocity in the rotational frame using the midpoint acceleration + PartVeloRotRef(1:3,iPart) = PartVeloRotRef(1:3,iPart) + Pt_local(1:3)*RotTimestep +ELSE + PartState(1:3,iPart) = PartState(1:3,iPart) + PartState(4:6,iPart) * RotTimestep +END IF + +END SUBROUTINE CalcPartPosInRotRef + + END MODULE MOD_part_RHS diff --git a/src/particles/particle_tools.f90 b/src/particles/particle_tools.f90 index 049df9ae5..63c1f524b 100644 --- a/src/particles/particle_tools.f90 +++ b/src/particles/particle_tools.f90 @@ -718,7 +718,7 @@ PPURE FUNCTION InRotRefFrameCheck(iPart) !----------------------------------------------------------------------------------------------------------------------------------! ! MODULES ! !----------------------------------------------------------------------------------------------------------------------------------! -USE MOD_Particle_Vars ,ONLY: PartState,RotRefFramRegion,RotRefFrameAxis,nRefFrameRegions +USE MOD_Particle_Vars ,ONLY: PartState,RotRefFrameRegion,RotRefFrameAxis,nRefFrameRegions !----------------------------------------------------------------------------------------------------------------------------------! IMPLICIT NONE ! INPUT / OUTPUT VARIABLES @@ -732,8 +732,8 @@ PPURE FUNCTION InRotRefFrameCheck(iPart) IF(nRefFrameRegions.GT.0) THEN InRotRefFrameCheck = .FALSE. DO iRegion = 1, nRefFrameRegions - IF((PartState(RotRefFrameAxis,iPart).GT.RotRefFramRegion(1,iRegion)).AND. & - (PartState(RotRefFrameAxis,iPart).LT.RotRefFramRegion(2,iRegion))) THEN + IF((PartState(RotRefFrameAxis,iPart).GT.RotRefFrameRegion(1,iRegion)).AND. & + (PartState(RotRefFrameAxis,iPart).LT.RotRefFrameRegion(2,iRegion))) THEN InRotRefFrameCheck = .TRUE. EXIT END IF @@ -1868,6 +1868,7 @@ SUBROUTINE IncreaseMaxParticleNumber(Amount) IF(ALLOCATED(PartTimeStep)) CALL ChangeSizeArray(PartTimeStep,PDM%maxParticleNumber,NewSize) IF(ALLOCATED(PartMPF)) CALL ChangeSizeArray(PartMPF,PDM%maxParticleNumber,NewSize) IF(ALLOCATED(PartVeloRotRef)) CALL ChangeSizeArray(PartVeloRotRef,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(LastPartVeloRotRef)) CALL ChangeSizeArray(LastPartVeloRotRef,PDM%maxParticleNumber,NewSize,0.) IF(ALLOCATED(PartStateIntEn)) CALL ChangeSizeArray(PartStateIntEn,PDM%maxParticleNumber,NewSize) IF(ALLOCATED(Pt_temp)) CALL ChangeSizeArray(Pt_temp,PDM%maxParticleNumber,NewSize,0.) @@ -2073,6 +2074,7 @@ SUBROUTINE ReduceMaxParticleNumber() IF(ALLOCATED(PartTimeStep)) CALL ChangeSizeArray(PartTimeStep,PDM%maxParticleNumber,NewSize) IF(ALLOCATED(PartMPF)) CALL ChangeSizeArray(PartMPF,PDM%maxParticleNumber,NewSize) IF(ALLOCATED(PartVeloRotRef)) CALL ChangeSizeArray(PartVeloRotRef,PDM%maxParticleNumber,NewSize,0.) +IF(ALLOCATED(LastPartVeloRotRef)) CALL ChangeSizeArray(LastPartVeloRotRef,PDM%maxParticleNumber,NewSize,0.) IF(ALLOCATED(PartStateIntEn)) CALL ChangeSizeArray(PartStateIntEn,PDM%maxParticleNumber,NewSize) IF(ALLOCATED(Pt_temp)) CALL ChangeSizeArray(Pt_temp,PDM%maxParticleNumber,NewSize,0.) @@ -2266,6 +2268,10 @@ SUBROUTINE ChangePartID(OldID,NewID) PartVeloRotRef(:,NewID)=PartVeloRotRef(:,OldID) PartVeloRotRef(:,OldID) = 0.0 END IF +IF(ALLOCATED(LastPartVeloRotRef)) THEN + LastPartVeloRotRef(:,NewID)=LastPartVeloRotRef(:,OldID) + LastPartVeloRotRef(:,OldID) = 0.0 +END IF IF(ALLOCATED(PartStateIntEn)) PartStateIntEn(:,NewID)=PartStateIntEn(:,OldID) IF(ALLOCATED(Pt_temp)) THEN diff --git a/src/particles/particle_vars.f90 b/src/particles/particle_vars.f90 index 828bb99d8..c88d73707 100644 --- a/src/particles/particle_vars.f90 +++ b/src/particles/particle_vars.f90 @@ -72,6 +72,7 @@ MODULE MOD_Particle_Vars ! ! 2nd index: 1:NParts REAL , ALLOCATABLE :: PartPosRef(:,:) ! (1:3,1:NParts) particles pos mapped to -1|1 space REAL , ALLOCATABLE :: PartVeloRotRef(:,:) ! (1:3,1:NParts) Velocity in the rotational reference frame +REAL , ALLOCATABLE :: LastPartVeloRotRef(:,:) ! (1:3,1:NParts) Last Velocity in the rotational reference frame REAL , ALLOCATABLE :: Pt(:,:) ! Derivative of PartState (vx,xy,vz) only ! since temporal derivative of position ! is the velocity. Thus we can take @@ -212,7 +213,7 @@ PPURE INTEGER FUNCTION ElemID_INTERFACE(iPart) INTEGER :: ParticleVecLengthOld ! Vector Length for Particle Push Calculation REAL :: MaxPartNumIncrease ! How much shall the PDM%MaxParticleNumber be increased if it is full INTEGER ,ALLOCATABLE :: nextFreePosition(:) ! =>NULL() ! next_free_Position(1:maxParticleNumber) - ! List of free Positon + ! List of free Position LOGICAL ,ALLOCATABLE :: ParticleInside(:) ! Particle_inside (1:maxParticleNumber) LOGICAL ,ALLOCATABLE :: InRotRefFrame(:) ! Check for RotRefFrame (1:maxParticleNumber) LOGICAL ,ALLOCATABLE :: dtFracPush(:) ! Push random fraction only @@ -299,7 +300,17 @@ PPURE INTEGER FUNCTION ElemID_INTERFACE(iPart) REAL :: RotRefFrameFreq ! frequency of rotational frame of reference REAL :: RotRefFrameOmega(3) ! angular velocity of rotational frame of reference INTEGER :: nRefFrameRegions ! number of rotational frame of reference regions -REAL, ALLOCATABLE :: RotRefFramRegion(:,:) ! MIN/MAX defintion for multiple rotational frame of reference region +REAL, ALLOCATABLE :: RotRefFrameRegion(:,:) ! MIN/MAX defintion for multiple rotational frame of reference region ! (i,RegionNumber), MIN:i=1, MAX:i=2 +! Rotational frame of reference: Subcycling +LOGICAL :: UseRotSubCycling ! Flag if subcycling is active +INTEGER :: nSubCyclingSteps ! Number of subcycling steps +REAL :: LastPartPosSubCycling(3) ! Last position before subcycling +REAL :: NewPosSubCycling(3) ! New particle position before subcycling +REAL :: PartVeloRotRefSubCycling(3) ! Velocity in the rotational reference frame before subcycling +REAL :: LastVeloRotRefSubCycling(3) ! Last Velocity in the rotational reference frame before subcycling +INTEGER :: GlobalElemIDSubCycling ! Element ID before subcycling +LOGICAL :: RotRefSubTimeStep ! Flag for loop that defines that the current time step is a subcycling step +LOGICAL :: InRotRefFrameSubCycling ! Check for RotRefFrame before subcycling !=================================================================================================================================== END MODULE MOD_Particle_Vars diff --git a/src/particles/surfacemodel/surfacemodel_tools.f90 b/src/particles/surfacemodel/surfacemodel_tools.f90 index 589ed3b05..8240e5a71 100644 --- a/src/particles/surfacemodel/surfacemodel_tools.f90 +++ b/src/particles/surfacemodel/surfacemodel_tools.f90 @@ -86,7 +86,7 @@ SUBROUTINE PerfectReflection(PartID,SideID,n_Loc,opt_Symmetry) !----------------------------------------------------------------------------------------------------------------------------------! USE MOD_Globals USE MOD_Particle_Boundary_Vars ,ONLY: PartBound -USE MOD_Particle_Vars ,ONLY: PartState,LastPartPos,PartSpecies,Species,PartLorentzType +USE MOD_Particle_Vars ,ONLY: PartState,LastPartPos,PartSpecies,Species,PartLorentzType,RotRefSubTimeStep,nSubCyclingSteps USE MOD_DSMC_Vars ,ONLY: DSMC, AmbipolElecVelo USE MOD_Globals_Vars ,ONLY: c2_inv #if defined(LSERK) @@ -101,7 +101,7 @@ SUBROUTINE PerfectReflection(PartID,SideID,n_Loc,opt_Symmetry) USE MOD_Particle_Tracking_Vars ,ONLY: TrackInfo USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, PartTimeStep, VarTimeStep USE MOD_TimeDisc_Vars ,ONLY: dt,RKdtFrac -USE MOD_Particle_Vars ,ONLY: PDM, PartVeloRotRef +USE MOD_Particle_Vars ,ONLY: PDM, PartVeloRotRef, RotRefFrameOmega USE MOD_part_RHS ,ONLY: CalcPartRHSRotRefFrame ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -119,6 +119,7 @@ SUBROUTINE PerfectReflection(PartID,SideID,n_Loc,opt_Symmetry) INTEGER :: locBCID, SpecID LOGICAL :: Symmetry REAL :: POI_vec(1:3) +REAL :: NormNewVeloPush(1:3) REAL :: dtVar !=================================================================================================================================== ! Initialize @@ -144,6 +145,7 @@ SUBROUTINE PerfectReflection(PartID,SideID,n_Loc,opt_Symmetry) END IF ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific) dtVar = dtVar * Species(SpecID)%TimeStepFactor +IF(RotRefSubTimeStep) dtVar = dtVar / REAL(nSubCyclingSteps) IF(PDM%InRotRefFrame(PartID)) THEN ! In case of RotRefFrame utilize the respective velocity @@ -185,31 +187,31 @@ SUBROUTINE PerfectReflection(PartID,SideID,n_Loc,opt_Symmetry) END IF END IF -! set particle position on face +! Set particle position on face LastPartPos(1:3,PartID) = POI_vec(1:3) - -! Determine the correct velocity in case of a rotational frame of reference -NewVeloPush(1:3) = PartState(4:6,PartID) -! In case of RotRefFrame, the velocity in the rotational reference frame is mirrored as well -IF(PDM%InRotRefFrame(PartID)) THEN - ! Mirror the velocity in the rotational frame - PartVeloRotRef(1:3,PartID) = PartVeloRotRef(1:3,PartID) - 2.*DOT_PRODUCT(PartVeloRotRef(1:3,PartID),n_loc)*n_loc - ! ALTERNATIVE: Transform the inertial velocity (which was mirrored) - ! PartVeloRotRef(1:3,PartID) = PartState(4:6,PartID) - CROSS(RotRefFrameOmega(1:3),LastPartPos(1:3,PartID)) -END IF - TrackInfo%PartTrajectory(1:3) = TrackInfo%PartTrajectory(1:3)-2.*DOT_PRODUCT(TrackInfo%PartTrajectory(1:3),n_loc)*n_loc - -! Check if rotational frame of reference is used, otherwise mirror the LastPartPos +! Check if rotational frame of reference is used, otherwise mirror the LastPartPos for new particle position IF(PDM%InRotRefFrame(PartID)) THEN POI_fak = 1.- (TrackInfo%lengthPartTrajectory-TrackInfo%alpha) / VECNORM(OldVelo*dtVar) - ! Add the acceleration due to new velocity vector at the POI - PartVeloRotRef(1:3,PartID) = PartVeloRotRef(1:3,PartID) + CalcPartRHSRotRefFrame(LastPartPos(1:3,PartID),PartVeloRotRef(1:3,PartID)) * dtVar * (1.0 - POI_fak) - ! IF(DOT_PRODUCT(PartVeloRotRef(1:3,PartID),n_loc).GT.0.) THEN - ! PartVeloRotRef(1:3,PartID) = PartVeloRotRef(1:3,PartID) - 2.*DOT_PRODUCT(PartVeloRotRef(1:3,PartID),n_loc)*n_loc - ! END IF - PartState(1:3,PartID) = LastPartPos(1:3,PartID) + (1.0 - POI_fak) * dtVar * PartVeloRotRef(1:3,PartID) + ! Determine the correct velocity for the subsequent push in case of a rotational frame of reference at POI + NewVeloPush(1:3) = PartState(4:6,PartID) + NewVeloPush(1:3) = NewVeloPush(1:3) - CROSS(RotRefFrameOmega(1:3),LastPartPos(1:3,PartID)) + NewVeloPush(1:3) = NewVeloPush(1:3) + CalcPartRHSRotRefFrame(LastPartPos(1:3,PartID),NewVeloPush(1:3)) * (1.0 - POI_fak) * dtVar + ! Make sure the NewVeloPush is pointing away from the wall + IF(DOT_PRODUCT(n_loc,NewVeloPush(1:3)).GT.0.) THEN + ! Normal component of new velo push v = (v dot n / |n|^2) * n, |n| = 1 + NormNewVeloPush(1:3) = DOT_PRODUCT(n_loc,NewVeloPush(1:3)) * n_loc + ! Nullyfy normal component and keeping rest of NewVeloPush + NewVeloPush(1:3) = NewVeloPush(1:3) - NormNewVeloPush(1:3) + ! Move particle a little bit into the domain to avoid losing particles + NewVeloPush(1:3) = NewVeloPush(1:3) - 1E-6 * n_loc + END IF + ! Store the new rotational reference frame velocity + PartVeloRotRef(1:3,PartID) = NewVeloPush(1:3) + ! Calc new particle position with NewVeloPush + PartState(1:3,PartID) = LastPartPos(1:3,PartID) + (1.0 - POI_fak) * dtVar * NewVeloPush(1:3) ELSE + ! Mirror the LastPartPos for new particle position PartState(1:3,PartID) = LastPartPos(1:3,PartID) + TrackInfo%PartTrajectory(1:3)*(TrackInfo%lengthPartTrajectory - TrackInfo%alpha) END IF @@ -291,7 +293,7 @@ SUBROUTINE DiffuseReflection(PartID,SideID,n_loc) 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 +USE MOD_Particle_Vars ,ONLY: PDM, RotRefFrameOmega,RotRefSubTimeStep,nSubCyclingSteps USE MOD_Particle_Tracking_Vars ,ONLY: TrackInfo USE MOD_part_RHS ,ONLY: CalcPartRHSRotRefFrame ! IMPLICIT VARIABLE HANDLING @@ -340,6 +342,7 @@ SUBROUTINE DiffuseReflection(PartID,SideID,n_loc) ! Species-specific time step IF(VarTimeStep%UseSpeciesSpecific) dtVar = dtVar * Species(SpecID)%TimeStepFactor +IF(RotRefSubTimeStep) dtVar = dtVar / REAL(nSubCyclingSteps) IF(PDM%InRotRefFrame(PartID)) THEN ! In case of RotRefFrame utilize the respective velocity @@ -511,8 +514,6 @@ SUBROUTINE SurfaceModelParticleEmission(n_loc, PartID, SideID, ProductSpec, Prod USE MOD_Particle_Boundary_Vars ,ONLY: Partbound, GlobalSide2SurfSide USE MOD_Particle_Mesh_Vars ,ONLY: SideInfo_Shared USE MOD_DSMC_Vars ,ONLY: DSMC, SamplingActive -USE MOD_Particle_Vars ,ONLY: usevMPF,PartMPF -USE MOD_part_tools ,ONLY: CalcRadWeightMPF USE MOD_Particle_Mesh_Vars ,ONLY: BoundsOfElem_Shared ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -531,7 +532,7 @@ SUBROUTINE SurfaceModelParticleEmission(n_loc, PartID, SideID, ProductSpec, Prod !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: iNewPart, NewPartID, locBCID, SurfSideID -REAL :: tang1(1:3), tang2(1:3), WallVelo(1:3), WallTemp, NewVelo(3), OldMPF, BoundsOfElemCenter(1:3),NewPos(1:3) +REAL :: tang1(1:3), tang2(1:3), WallVelo(1:3), WallTemp, NewVelo(3), BoundsOfElemCenter(1:3),NewPos(1:3) REAL,PARAMETER :: eps=1e-6 REAL,PARAMETER :: eps2=1.0-eps !=================================================================================================================================== @@ -560,15 +561,8 @@ SUBROUTINE SurfaceModelParticleEmission(n_loc, PartID, SideID, ProductSpec, Prod NewVelo(1:3) = tang1(1:3)*NewVelo(1) + tang2(1:3)*NewVelo(2) - n_Loc(1:3)*NewVelo(3) + WallVelo(1:3) ! Create new position by using POI and moving the particle by eps in the direction of the element center NewPos(1:3) = eps*BoundsOfElemCenter(1:3) + eps2*POI_vec(1:3) - IF(usevMPF)THEN - ! Get MPF of old particle - OldMPF = PartMPF(PartID) - ! New particle acquires the MPF of the impacting particle (not necessarily the MPF of the newly created particle species) - CALL CreateParticle(ProductSpec,NewPos(1:3),GlobElemID,NewVelo(1:3),0.,0.,0.,NewPartID=NewPartID, NewMPF=OldMPF) - ELSE - ! New particle acquires the MPF of the new particle species - CALL CreateParticle(ProductSpec,NewPos(1:3),GlobElemID,NewVelo(1:3),0.,0.,0.,NewPartID=NewPartID) - END IF ! usevMPF + ! Create new particle: in case of vMPF oder VarTimeStep, new particle inherits the values of the old particle + CALL CreateParticle(ProductSpec,NewPos(1:3),GlobElemID,NewVelo(1:3),0.,0.,0.,OldPartID=PartID,NewPartID=NewPartID) ! Adding the energy that is transferred from the surface onto the internal energies of the particle CALL SurfaceModelEnergyAccommodation(NewPartID,locBCID,WallTemp) ! Sampling of newly created particles diff --git a/src/particles/tracking/particle_triatracking.f90 b/src/particles/tracking/particle_triatracking.f90 index bf6040499..acd8ea48b 100644 --- a/src/particles/tracking/particle_triatracking.f90 +++ b/src/particles/tracking/particle_triatracking.f90 @@ -40,10 +40,18 @@ SUBROUTINE ParticleTriaTracking() ! MODULES USE MOD_Preproc USE MOD_Globals -USE MOD_Particle_Vars ,ONLY: PEM,PDM,InterPlanePartNumber, InterPlanePartIndx +USE MOD_Particle_Vars ,ONLY: PEM,PDM,InterPlanePartNumber, InterPlanePartIndx, UseRotSubCycling,nSubCyclingSteps +USE MOD_Particle_Vars ,ONLY: RotRefSubTimeStep, NewPosSubCycling, GlobalElemIDSubCycling, LastPartPosSubCycling +USE MOD_Particle_Vars ,ONLY: InRotRefFrameSubCycling, PartVeloRotRefSubCycling, LastVeloRotRefSubCycling USE MOD_DSMC_Vars ,ONLY: RadialWeighting USE MOD_DSMC_Symmetry ,ONLY: DSMC_2D_RadialWeighting, DSMC_2D_SetInClones USE MOD_part_tools ,ONLY: ParticleOnProc +!----- Used for RotRef Subcycling +USE MOD_part_RHS ,ONLY: CalcPartPosInRotRef +USE MOD_TimeDisc_Vars ,ONLY: dt +USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, PartTimeStep, VarTimeStep, Species, PartState, LastPartPos, PartSpecies +USE MOD_Particle_Vars ,ONLY: PartVeloRotRef, LastPartVeloRotRef +!----- ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -59,7 +67,8 @@ SUBROUTINE ParticleTriaTracking() LOGICAL :: doParticle LOGICAL :: doPartInExists #endif -INTEGER :: i, InterPartID +INTEGER :: i, InterPartID, iStep +REAL :: dtVar !----------------------------------------------------------------------------------------------------------------------------------- !=================================================================================================================================== #ifdef IMPA @@ -81,7 +90,60 @@ SUBROUTINE ParticleTriaTracking() #else IF (PDM%ParticleInside(i)) THEN #endif /*IMPA*/ - CALL SingleParticleTriaTracking(i=i) + IF(UseRotSubCycling) THEN +!--- Store Particle informations befor tracking for sub-cycling. +!--- it must be stored before first call of "SingleParticleTriaTracking" +!--- because particle informations like LastPartPos & PartState can be changed within "SingleParticleTriaTracking" +!--- e.g. in RotPeriodicInterPlaneBoundary + RotRefSubTimeStep=.TRUE. + LastPartPosSubCycling(1:3) = LastPartPos(1:3,i) + NewPosSubCycling(1:3) = PartState(1:3,i) + PartVeloRotRefSubCycling(1:3) = PartVeloRotRef(1:3,i) + LastVeloRotRefSubCycling(1:3) = LastPartVeloRotRef(1:3,i) + GlobalElemIDSubCycling = PEM%LastGlobalElemID(i) + InRotRefFrameSubCycling = PDM%InRotRefFrame(i) + IF(RotRefSubTimeStep) THEN +!--- split time step in 10 sub-steps + IF (UseVarTimeStep) THEN + dtVar = dt * PartTimeStep(i) + ELSE + dtVar = dt + END IF + IF(VarTimeStep%UseSpeciesSpecific) dtVar = dtVar * Species(PartSpecies(i))%TimeStepFactor + dtVar = dtVar / REAL(nSubCyclingSteps) +!--- Reset Particle Push from timedisc + PartState(1:3,i) = LastPartPosSubCycling(1:3) + LastPartPos(1:3,i) = LastPartPosSubCycling(1:3) + PartVeloRotRef(1:3,i)=LastVeloRotRefSubCycling(1:3) + PEM%GlobalElemID(i) = GlobalElemIDSubCycling + PEM%LastGlobalElemID(i)= GlobalElemIDSubCycling + LastPartVeloRotRef(1:3,i) = LastVeloRotRefSubCycling(1:3) + PDM%ParticleInside(i) = .TRUE. + PDM%InRotRefFrame(i) = InRotRefFrameSubCycling +!--- Loop over sub time steps + DO iStep = 1, nSubCyclingSteps + IF (PDM%ParticleInside(i)) THEN + IF(iStep.GT.1) THEN + LastPartPos(1:3,i)=PartState(1:3,i) + PEM%LastGlobalElemID(i)=PEM%GlobalElemID(i) + LastPartVeloRotRef(1:3,i)=PartVeloRotRef(1:3,i) + END IF + CALL CalcPartPosInRotRef(i, dtVar) + CALL SingleParticleTriaTracking(i=i) + END IF + END DO + END IF +!--- Reset stored particle informations + RotRefSubTimeStep = .FALSE. + LastPartPosSubCycling = 0.0 + NewPosSubCycling = 0.0 + PartVeloRotRefSubCycling = 0.0 + LastVeloRotRefSubCycling = 0.0 + GlobalElemIDSubCycling = 0 + InRotRefFrameSubCycling = .FALSE. + ELSE + CALL SingleParticleTriaTracking(i=i) + END IF END IF ! Particle treatment for an axisymmetric simulation (cloning/deleting particles) IF(RadialWeighting%PerformCloning) THEN @@ -95,7 +157,53 @@ SUBROUTINE ParticleTriaTracking() DO i = 1,InterPlanePartNumber InterPartID = InterPlanePartIndx(i) PDM%ParticleInside(InterPartID) = .TRUE. - CALL SingleParticleTriaTracking(i=InterPartID,IsInterPlanePart=.TRUE.) + IF(UseRotSubCycling) THEN + RotRefSubTimeStep=.TRUE. + LastPartPosSubCycling(1:3) = LastPartPos(1:3,InterPartID) + NewPosSubCycling(1:3) = PartState(1:3,InterPartID) + PartVeloRotRefSubCycling(1:3) = PartVeloRotRef(1:3,InterPartID) + LastVeloRotRefSubCycling(1:3) = LastPartVeloRotRef(1:3,InterPartID) + GlobalElemIDSubCycling = PEM%LastGlobalElemID(InterPartID) + InRotRefFrameSubCycling = PDM%InRotRefFrame(InterPartID) + !--- split time step in 10 sub-steps + IF (UseVarTimeStep) THEN + dtVar = dt * PartTimeStep(InterPartID) + ELSE + dtVar = dt + END IF + IF(VarTimeStep%UseSpeciesSpecific) dtVar = dtVar * Species(PartSpecies(InterPartID))%TimeStepFactor + dtVar = dtVar / REAL(nSubCyclingSteps) + !--- reset Particle Push + PartState(1:3,InterPartID) = LastPartPosSubCycling(1:3) + LastPartPos(1:3,InterPartID) = LastPartPosSubCycling(1:3) + PartVeloRotRef(1:3,InterPartID)=LastVeloRotRefSubCycling(1:3) + PEM%GlobalElemID(InterPartID) = GlobalElemIDSubCycling + PEM%LastGlobalElemID(InterPartID)= GlobalElemIDSubCycling + LastPartVeloRotRef(1:3,InterPartID) = LastVeloRotRefSubCycling(1:3) + PDM%InRotRefFrame(InterPartID) = InRotRefFrameSubCycling + !--- Loop over sub time steps + DO iStep = 1, nSubCyclingSteps + IF (PDM%ParticleInside(InterPartID)) THEN + IF(iStep.GT.1) THEN + LastPartPos(1:3,InterPartID)=PartState(1:3,InterPartID) + PEM%LastGlobalElemID(InterPartID)=PEM%GlobalElemID(InterPartID) + LastPartVeloRotRef(1:3,InterPartID)=PartVeloRotRef(1:3,InterPartID) + END IF + CALL CalcPartPosInRotRef(InterPartID, dtVar) + CALL SingleParticleTriaTracking(i=InterPartID,IsInterPlanePart=.TRUE.) + END IF + END DO +!--- Reset stored particle informations + RotRefSubTimeStep = .FALSE. + LastPartPosSubCycling = 0.0 + NewPosSubCycling = 0.0 + PartVeloRotRefSubCycling = 0.0 + LastVeloRotRefSubCycling = 0.0 + GlobalElemIDSubCycling = 0 + InRotRefFrameSubCycling = .FALSE. + ELSE + CALL SingleParticleTriaTracking(i=InterPartID,IsInterPlanePart=.TRUE.) + END IF ! Particle treatment for an axisymmetric simulation (cloning/deleting particles) IF(RadialWeighting%PerformCloning) THEN IF(PDM%ParticleInside(InterPartID).AND.(ParticleOnProc(InterPartID))) THEN @@ -392,7 +500,6 @@ SUBROUTINE SingleParticleTriaTracking(i,IsInterPlanePart) ELSE CALL GetBoundaryInteraction(i,SideID,flip,ElemID,crossedBC,TriNum=TriNum) END IF - IF(.NOT.PDM%ParticleInside(i)) PartisDone = .TRUE. #if USE_LOADBALANCE IF (OldElemID.GE.offsetElem+1.AND.OldElemID.LE.offsetElem+PP_nElems) CALL LBElemSplitTime(OldElemID-offsetElem,tLBStart) @@ -423,7 +530,6 @@ SUBROUTINE SingleParticleTriaTracking(i,IsInterPlanePart) END IF END IF ! SideInfo_Shared(SIDE_BCID,SideID).GT./.LE. 0 CNElemID = GetCNElemID(ElemID) - IF (CNElemID.LT.1) THEN IPWRITE(UNIT_StdOut,*) "VECNORM(PartState(1:3,i)-LastPartPos(1:3,i)): ", VECNORM(PartState(1:3,i)-LastPartPos(1:3,i)) IPWRITE(UNIT_StdOut,*) " PartState(1:3,i) : ", PartState(1:3,i) diff --git a/src/timedisc/timedisc_TimeStep_DSMC.f90 b/src/timedisc/timedisc_TimeStep_DSMC.f90 index 73b551c80..68f0d2bbd 100644 --- a/src/timedisc/timedisc_TimeStep_DSMC.f90 +++ b/src/timedisc/timedisc_TimeStep_DSMC.f90 @@ -37,7 +37,7 @@ SUBROUTINE TimeStep_DSMC() #ifdef PARTICLES 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 +USE MOD_Particle_Vars ,ONLY: UseRotRefFrame, RotRefFrameOmega, PartVeloRotRef, LastPartVeloRotRef USE MOD_Particle_Vars ,ONLY: WriteMacroSurfaceValues, Symmetry, Species, PartSpecies USE MOD_Particle_Vars ,ONLY: UseVarTimeStep, PartTimeStep, VarTimeStep USE MOD_Particle_Vars ,ONLY: UseSplitAndMerge @@ -52,7 +52,7 @@ 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 +USE MOD_part_RHS ,ONLY: CalcPartRHSRotRefFrame, CalcPartPosInRotRef USE MOD_part_pos_and_velo ,ONLY: SetParticleVelocity USE MOD_Part_Tools ,ONLY: InRotRefFrameCheck USE MOD_Part_Tools ,ONLY: CalcPartSymmetryPos @@ -74,8 +74,6 @@ SUBROUTINE TimeStep_DSMC() ! LOCAL VARIABLES INTEGER :: iPart REAL :: timeEnd, timeStart, dtVar, RandVal -! Rotational frame of reference -REAL :: Pt_local(1:3), Pt_local_old(1:3), VeloRotRef_half(1:3), PartState_half(1:3) #if USE_LOADBALANCE REAL :: tLBStart #endif /*USE_LOADBALANCE*/ @@ -141,23 +139,8 @@ SUBROUTINE TimeStep_DSMC() PEM%LastGlobalElemID(iPart)=PEM%GlobalElemID(iPart) END IF IF(UseRotRefFrame) THEN - IF(PDM%InRotRefFrame(iPart)) THEN - ! Midpoint method - ! calculate the acceleration (force / mass) at the current time step - Pt_local_old(1:3) = CalcPartRHSRotRefFrame(PartState(1:3,iPart), PartVeloRotRef(1:3,iPart)) - ! estimate the mid-point velocity in the rotational frame - VeloRotRef_half(1:3) = PartVeloRotRef(1:3,iPart) + 0.5*Pt_local_old(1:3)*dtVar - ! estimate the mid-point position - PartState_half(1:3) = PartState(1:3,iPart) + 0.5*PartVeloRotRef(1:3,iPart)*dtVar - ! calculate the acceleration (force / mass) at the mid-point - Pt_local(1:3) = CalcPartRHSRotRefFrame(PartState_half(1:3), VeloRotRef_half(1:3)) - ! update the position using the mid-point velocity in the rotational frame - PartState(1:3,iPart) = PartState(1:3,iPart) + VeloRotRef_half(1:3)*dtVar - ! update the velocity in the rotational frame using the mid-point acceleration - PartVeloRotRef(1:3,iPart) = PartVeloRotRef(1:3,iPart) + Pt_local(1:3)*dtVar - ELSE - PartState(1:3,iPart) = PartState(1:3,iPart) + PartState(4:6,iPart) * dtVar - END IF + LastPartVeloRotRef(1:3,iPart)=PartVeloRotRef(1:3,iPart) + CALL CalcPartPosInRotRef(iPart, dtVar) ELSE PartState(1:3,iPart) = PartState(1:3,iPart) + PartState(4:6,iPart) * dtVar END IF @@ -266,6 +249,5 @@ SUBROUTINE TimeStep_DSMC() END SUBROUTINE TimeStep_DSMC - END MODULE MOD_TimeStep #endif