Merge branch 'feature.RotRefFrame.Subcycling' into 'master.dev'

[feature.RotRefFrame.Subcycling] Feature.rot ref frame.subcycling

Closes #235

See merge request piclas/piclas!862

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`. 
+`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

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

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

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

regressioncheck/CHE_DSMC/Rotational_Reference_Frame_Subcycling/command_line.ini

@@ -0,0 +1,2 @@
+MPI=1
+cmd_suffix=DSMC.ini

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

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

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

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.

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
 

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
 ! =============================================================================== !

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