xcompact3d · fangjian19 · Sep 27, 2024 · Sep 27, 2024 · Sep 27, 2024 · Oct 8, 2024
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
@@ -17,6 +17,7 @@ if (${BUILD_TESTING_FULL})
   add_subdirectory(Gravity-current)
   add_subdirectory(Mixing-layer)
   add_subdirectory(MHD)
+  add_subdirectory(Particle-Tracking)
   add_subdirectory(TBL-Turbulent-Boundary-Layer)
   add_subdirectory(Wind-Turbine/NREL-5MW_ALM)
 endif (${BUILD_TESTING_FULL})
diff --git a/examples/MHD/input_otv.i3d b/examples/MHD/input_otv.i3d
@@ -42,7 +42,7 @@ inflow_noise = 0.0    ! Turbulence intensity (1=100%) !! Inflow condition
 ! Time stepping
 dt = 0.0001           ! Time step
 ifirst = 1            ! First iteration
-ilast  = 20           ! Last iteration
+ilast  = 20000        ! Last iteration
 
 ! Enable modelling tools
 ilesmod=0             ! if 0 then DNS

diff --git a/examples/Particle-Tracking/CMakeLists.txt b/examples/Particle-Tracking/CMakeLists.txt
@@ -0,0 +1,33 @@
+# TGV flow with passive particle trackers
+set(case "TGV-Particles")
+set(case_dir "${test_dir}/${case}")
+file(MAKE_DIRECTORY ${case_dir})
+set(input_file "input_tgv_particles.i3d")
+install(DIRECTORY DESTINATION ${CMAKE_INSTALL_PREFIX}/examples/${case})
+install(FILES ${input_file} DESTINATION ${CMAKE_INSTALL_PREFIX}/examples/${case})
+if(ADIOS2_FOUND)
+  install(FILES adios2_config.xml DESTINATION ${CMAKE_INSTALL_PREFIX}/examples/${case})
+endif()
+# If testing active add test
+file(COPY ${input_file} DESTINATION ${case_dir})
+if(ADIOS2_FOUND)
+  file(COPY adios2_config.xml DESTINATION ${case_dir})
+endif()
+add_test(NAME ${case} COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} $<TARGET_FILE:xcompact3d> ${input_file} WORKING_DIRECTORY ${case_dir})
+
+# Channel flow with constantly injected passive particle trackers
+set(case "Channel-Particles")
+set(case_dir "${test_dir}/${case}")
+file(MAKE_DIRECTORY ${case_dir})
+set(input_file "input_channel_particles.i3d")
+install(DIRECTORY DESTINATION ${CMAKE_INSTALL_PREFIX}/examples/${case})
+install(FILES ${input_file} DESTINATION ${CMAKE_INSTALL_PREFIX}/examples/${case})
+if(ADIOS2_FOUND)
+  install(FILES adios2_config.xml DESTINATION ${CMAKE_INSTALL_PREFIX}/examples/${case})
+endif()
+# If testing active add test
+file(COPY ${input_file} DESTINATION ${case_dir})
+if(ADIOS2_FOUND)
+  file(COPY adios2_config.xml DESTINATION ${case_dir})
+endif()
+add_test(NAME ${case} COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} $<TARGET_FILE:xcompact3d> ${input_file} WORKING_DIRECTORY ${case_dir})
diff --git a/examples/Particle-Tracking/input_channel_particles.i3d b/examples/Particle-Tracking/input_channel_particles.i3d
@@ -0,0 +1,163 @@
+! -*- mode: f90 -*-
+
+!===================
+&BasicParam
+!===================
+
+! Flow type (1=Lock-exchange, 2=TGV, 3=Channel, 4=Periodic hill, 5=Cylinder, 6=dbg-schemes)
+itype = 3
+
+! Domain decomposition
+p_row=0               ! Row partition
+p_col=0                ! Column partition
+
+! Mesh
+nx=32                 ! X-direction nodes
+ny=65                 ! Y-direction nodes
+nz=32                 ! Z-direction nodes
+istret = 2            ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
+beta = 0.259065151    ! Refinement parameter (beta)
+
+! Domain
+xlx = 6.0         ! Lx (Size of the box in x-direction)
+yly = 2.0        ! Ly (Size of the box in y-direction)
+zlz = 3.0        ! Lz (Size of the box in z-direction)
+
+! Flow parameters
+re  = 4200.0           ! nu=1/re (Kinematic Viscosity)
+cpg = F                ! if cpg=T, then re is friction Reynolds number
+
+! Time stepping
+dt = 0.005            ! Time step
+ifirst = 1            ! First iteration
+ilast  = 200        ! Last iteration
+
+! Enable modelling tools
+ilesmod=0             ! if 0 then DNS
+numscalar = 0         ! How many scalars? (Set to zero to disable scalars)
+iibm=0                ! Flag for immersed boundary method
+
+! Boundary and initial conditions
+iin = 4               ! Inflow conditions (1: classic, 2: turbinit, 4: SEM)
+u1  = 1.              ! u1 (max velocity) (for inflow condition)
+u2  = 1.              ! u2 (min velocity) (for inflow condition)
+init_noise = 0.03     ! Turbulence intensity (1=100%) !! Initial condition
+inflow_noise = 0.0    ! Turbulence intensity (1=100%) !! Inflow condition
+idir_stream  = 1      ! Index of the streamwise direction (1=X, 3=Z, 2 not available)
+
+nclx1 = 0
+nclxn = 0
+ncly1 = 2
+nclyn = 2
+nclz1 = 0
+nclzn = 0
+
+mhd_active = .false.
+particle_active =.true.
+/End
+
+!====================
+&NumOptions
+!====================
+
+! Spatial derivatives
+ifirstder = 4         ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact)
+isecondder = 5        ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact, 5->hyperviscous 6th)
+ipinter = 3           ! interpolation scheme (1: classic, 2: optimized, 3: optimized agressive)
+
+! Time scheme
+iimplicit = 0
+itimescheme = 5       ! Time integration scheme (1->Euler,2->AB2, 3->AB3, 4->AB4,5->RK3,6->RK4, 7-->CN2+AB3)
+
+! Dissipation control
+nu0nu = 4.0           ! Ratio between hyperviscosity/viscosity at nu
+cnu = 0.44            ! Ratio between hypervisvosity at k_m=2/3pi and k_c= pi
+
+/End
+
+!=================
+&InOutParam
+!=================
+
+! Basic I/O
+irestart = 0           ! Read initial flow field 
+icheckpoint = 10000    ! Frequency for writing backup file
+ioutput = 20         ! Frequency for visualization
+ilist    = 10         ! Frequency for the output to screen
+nvisu = 1              ! Size for visualisation collection
+
+/End
+
+!=================
+&Statistics
+!=================
+
+wrotation = 0.d0      ! rotation speed to trigger turbulence
+spinup_time = 0    ! number of time steps with a rotation to trigger turbulence
+nstat = 1             ! Size arrays for statistic collection
+initstat = 1000000      ! Time steps after which statistics are collected
+
+/End
+
+!########################
+! OPTIONAL PARAMETERS
+!#######################
+
+!================
+&ScalarParam
+!================
+
+Sc(1) = 1.0           ! Schmidt number
+
+nclxS1 = 0
+nclxSn = 0
+nclyS1 = 2
+nclySn = 2
+nclzS1 = 0
+nclzSn = 0
+
+/End
+
+!====================
+&MHDParam
+!====================
+
+mhd_equation = .false.
+hartmann = 4.d0
+
+/End
+
+!====================
+&ParTrack
+!====================
+
+bc_particle = 'peri','peri','refl','refl','peri','peri'
+initype_particle = 'uniform'
+n_particles = 64
+particle_inject_period = 0.1
+
+/End
+
+!================
+&LESModel
+!================
+
+iles = 0              ! LES Model (1: Phys Smag, 2: Phys WALE, 3: Phys dyn. Smag, 4: iSVV, 5: dyn SEV)
+smagcst = 0.14        ! Smagorinsky constant
+walecst = 0.5         ! WALES Model Coefficient
+iwall = 0             ! Enable wall modeling
+
+/End
+
+!================
+&WallModel
+!================
+
+SmagWallDamp = 0      ! Smagorinsky damping function for ABL if 1
+
+/End
+
+&CASE
+/End
+
+
diff --git a/examples/Particle-Tracking/input_tgv_particles.i3d b/examples/Particle-Tracking/input_tgv_particles.i3d
@@ -0,0 +1,134 @@
+! -*- mode: f90 -*-
+
+!===================
+&BasicParam
+!===================
+
+! Flow type (1=Lock-exchange, 2=TGV, 3=Channel, 4=Periodic hill, 5=Cylinder, 6=dbg-schemes)
+itype = 2
+
+! Domain decomposition
+p_row=0               ! Row partition
+p_col=0               ! Column partition
+
+! Mesh
+nx=32                ! X-direction nodes
+ny=32                ! Y-direction nodes
+nz=32                ! Z-direction nodes
+istret = 0            ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
+beta = 0.259065151    ! Refinement parameter (beta)
+
+! Domain
+xlx = 6.28318530717958 ! Lx (Size of the box in x-direction)
+yly = 6.28318530717958 ! Ly (Size of the boy in y-direction)
+zlz = 6.28318530717958 ! Lz (Size of the boz in z-direction)
+
+! Boundary conditions
+nclx1 = 0 
+nclxn = 0 
+ncly1 = 0 
+nclyn = 0 
+nclz1 = 0 
+nclzn = 0 
+
+! Flow parameters
+iin = 1               ! Inflow conditions (1: classic, 2: turbinit)
+re = 1600.             ! nu=1/re (Kinematic Viscosity)
+u1 = 8.               ! u1 (max velocity) (for inflow condition)
+u2 = 8.               ! u2 (min velocity) (for inflow condition)
+init_noise  = 0.0     ! Turbulence intensity (1=100%) !! Initial condition
+inflow_noise = 0.0    ! Turbulence intensity (1=100%) !! Inflow condition
+
+! Time stepping
+dt = 0.001            ! Time step
+ifirst= 1            ! First iteration
+ilast = 200         ! Last iteration
+
+! Enable modelling tools
+ilesmod=0             ! if 0 then DNS
+iscalar=0             ! If iscalar=0 (no scalar), if iscalar=1 (scalar)
+iibm=0                ! Flag for immersed boundary method
+
+! Enable io
+ivisu=1               ! Store snapshots
+ipost=1               ! Do online postprocessing
+
+ifilter=0
+C_filter=0.49
+
+particle_active=.true.
+/End
+
+!====================
+&NumOptions
+!====================
+
+! Spatial derivatives
+ifirstder = 4         ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact)
+isecondder = 5        ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact, 5->hyperviscous 6th)
+
+! Time scheme
+itimescheme = 5       ! Time integration scheme (1->Euler,2->AB2, 3->AB3, 4->AB4,5->RK3,6->RK4)
+
+/End
+
+!=================
+&InOutParam
+!=================
+
+! Basic I/O
+irestart = 0          ! Read initial flow field ?
+icheckpoint = 1000     ! Frequency for writing backup file
+ioutput = 100         ! Frequency for visualization
+ilist = 5             ! Frequency for the output
+nvisu = 1             ! Size for visualisation collection
+
+/End
+
+!=================
+&Statistics
+!=================
+
+nstat = 1             ! Size arrays for statistic collection
+
+/End
+
+!====================
+&MHDParam
+!====================
+
+mhd_active = .false.
+mhd_equation=.true.
+Stuart = 100.d0
+Rem    = 100.d0
+
+/End
+
+!====================
+&ParTrack
+!====================
+
+bc_particle = 'peri','peri','peri','peri','peri','peri' 
+initype_particle = 'uniform'
+n_particles = 512
+particle_inject_period = 0.0
+
+/End
+
+!#######################
+! OPTIONAL PARAMETERS
+!#######################
+
+!=================
+&LESModel
+!=================
+jles=0
+smagcst=0.1
+walecst=0.5 
+iwall=0
+/End
+
+!=================
+&CASE
+!=================
+/End
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -44,6 +44,7 @@ add_executable(xcompact3d
                     mptool.f90
                     navier.f90
                     parameters.f90
+                    particle.f90
                     poisson.f90
                     probes.f90
                     schemes.f90

diff --git a/src/Case-TGV.f90 b/src/Case-TGV.f90
@@ -30,12 +30,17 @@ subroutine init_tgv (ux1,uy1,uz1,ep1,phi1)
     use MPI
     use mhd,  only: Bm,Bmean
     use visu, only: write_snapshot, end_snapshot
+    use decomp_2d, only : ph1
+    use var, only : nzmsize, numscalar, nrhotime, npress
 
     implicit none
 
     real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ux1,uy1,uz1,ep1
     real(mytype),dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi1
 
+    real(mytype),dimension(xsize(1),xsize(2),xsize(3),nrhotime) :: rho1
+    real(mytype),dimension(ph1%zst(1):ph1%zen(1), ph1%zst(2):ph1%zen(2), nzmsize, npress) :: pp3
+
     real(mytype) :: y,r,um,r3,x,z,h,ct
     real(mytype) :: cx0,cy0,cz0,hg,lg
     integer :: k,j,i,ierror,is,code
@@ -104,6 +109,7 @@ subroutine init_tgv (ux1,uy1,uz1,ep1,phi1)
 #endif
 
     return
+
   end subroutine init_tgv
   !********************************************************************
 
@@ -594,6 +600,7 @@ subroutine visu_tgv(ux1, uy1, uz1, num)
     endif
 
   end subroutine visu_tgv
+
   !############################################################################
   subroutine suspended(phi1,vol1,mp1)