More logical setup of the grid clustering (#85)

* added more logical tweaking for the grid clustering * changed `gtype` enumerator * minor documentation fix
CaNS-World · Sep 1, 2023 · 859e1d2 · 859e1d2
1 parent ff9920f
commit 859e1d2
Show file tree

Hide file tree

Showing 25 changed files with 77 additions and 48 deletions.
diff --git a/docs/INFO_INPUT.md b/docs/INFO_INPUT.md
@@ -5,7 +5,7 @@ Consider the following input file as example (corresponds to a turbulent plane c
 ~~~
 512 256 144              ! itot, jtot, ktot
 6. 3. 1.                 ! lx, ly, lz
-0.                       ! gr
+0 1.                     ! gtype, gr
 .95 1.0e5                ! cfl, dtmin
 5640.                    ! visci
 poi                      ! inivel
@@ -34,14 +34,19 @@ T F F                    ! is_forced(1:3)
 ~~~
 512 256 144              ! itot, jtot, ktot
 6. 3. 1.                 ! lx, ly, lz
-0.                       ! gr
+0 0.                     ! gtype, gr
 ~~~
 
 These lines set the computational grid.
 
 `itot, jtot, ktot ` and `lx, ly, lz` are the **number of points**  and **domain length** in each direction.
 
-`gr` is a **grid stretching parameter** that tweaks the non-uniform grid in the third direction; zero `gr` implies no stretching. See `initgrid.f90` for more details.
+`gtype` and `gr` are the **grid stretching type** and **grid stretching parameter** that tweak the non-uniform grid in the third direction; zero `gr` implies no stretching. See `initgrid.f90` for more details. The following options are available for `gtype`:
+
+* `1`: grid clustered towards both ends (default)
+* `2`: grid clustered towards the lower end
+* `3`: grid clustered towards the upper end
+* `4`: grid clustered towards the middle
 
 ---
 
@@ -106,9 +111,9 @@ These lines set the simulation termination criteria and whether the simulation s
 
 `stop_type` sets which criteria for terminating the simulation are to be used (more than one can be selected, and at least one of them must be `T`)
 
-* `stop_type(1)`, if true (`T`), the simulation will terminate after `nstep` time steps have been simulated;
-* `stop_type(2)`, if true (`T`), the simulation will terminate after `time_max` physical time units have been reached;
-* `stop_type(3)`, if true (`T`), the simulation will terminate after `tw_max` simulation wall-clock time (in hours) has been reached;
+* `stop_type(1)`, if true (`T`), the simulation will terminate after `nstep` time steps have been simulated
+* `stop_type(2)`, if true (`T`), the simulation will terminate after `time_max` physical time units have been reached
+* `stop_type(3)`, if true (`T`), the simulation will terminate after `tw_max` simulation wall-clock time (in hours) has been reached
 
 a checkpoint file `fld.bin` will be saved before the simulation is terminated.
 
@@ -126,12 +131,12 @@ a checkpoint file `fld.bin` will be saved before the simulation is terminated.
 
 These lines set the frequency of time step checking and output:
 
-* every `icheck` time steps **the new time step size** is computed according to the new stability criterion and cfl (above);
-* every `iout0d` time steps **history files with global scalar variables** are appended; currently the forcing pressure gradient and time step history are reported;
-* every `iout1d` time steps **1d profiles** are written (e.g. velocity and its moments) to a file;
-* every `iout2d` time steps **2d slices of a 3d scalar field** are written to a file;
-* every `iout3d` time steps **3d scalar fields** are written to a file;
-* every `isave`  time steps a **checkpoint file** is written (`fld_???????.bin`), and a symbolic link for the restart file, `fld.bin`, will point to this last save so that, by default, the last saved checkpoint file is used to restart the simulation.
+* every `icheck` time steps **the new time step size** is computed according to the new stability criterion and cfl (above)
+* every `iout0d` time steps **history files with global scalar variables** are appended; currently the forcing pressure gradient and time step history are reported
+* every `iout1d` time steps **1d profiles** are written (e.g. velocity and its moments) to a file
+* every `iout2d` time steps **2d slices of a 3d scalar field** are written to a file
+* every `iout3d` time steps **3d scalar fields** are written to a file
+* every `isave`  time steps a **checkpoint file** is written (`fld_???????.bin`), and a symbolic link for the restart file, `fld.bin`, will point to this last save so that, by default, the last saved checkpoint file is used to restart the simulation
 
 1d, 2d and 3d outputs can be tweaked modifying files `out?d.h90`, and re-compiling the source. See also `output.f90` for more details.
 
@@ -152,17 +157,17 @@ These lines set the boundary conditions (BC).
 
 The **type** (BC) for each field variable are set by a row of six characters, `X0 X1  Y0 Y1  Z0 Z1` where,
 
-* `X0` `X1` set the type of BC the field variable for the **lower** and **upper** boundaries in `x`;
-* `Y0` `Y1` set the type of BC the field variable for the **lower** and **upper** boundaries in `y`;
-* `Z0` `Z1` set the type of BC the field variable for the **lower** and **upper** boundaries in `z`.
+* `X0` `X1` set the type of BC the field variable for the **lower** and **upper** boundaries in `x`
+* `Y0` `Y1` set the type of BC the field variable for the **lower** and **upper** boundaries in `y`
+* `Z0` `Z1` set the type of BC the field variable for the **lower** and **upper** boundaries in `z`
 
 The four rows correspond to the three velocity components, and pressure, i.e. `u`, `v`, `w`, and `p`.
 
 The following options are available:
 
-* `P` periodic;
-* `D` Dirichlet;
-* `N` Neumann.
+* `P` periodic
+* `D` Dirichlet
+* `N` Neumann
 
 The **last four rows** follow the same logic, but now for the BC **values** (dummy for a periodic direction).
 

diff --git a/examples/_manuscript_lid_driven_cavity/dns.in b/examples/_manuscript_lid_driven_cavity/dns.in
@@ -1,6 +1,6 @@
 128 128 128              ! itot, jtot, ktot
 1. 1. 1.                 ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 zer                      ! inivel

diff --git a/examples/_manuscript_taylor_green_vortex/dns.in b/examples/_manuscript_taylor_green_vortex/dns.in
@@ -1,6 +1,6 @@
 512 512 512              ! itot, jtot, ktot
 6.2831853071795 6.283185307179586 6.283185307179586 ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1600.                    ! visci
 tgv                      ! inivel

diff --git a/examples/_manuscript_turbulent_channel/dns.in b/examples/_manuscript_turbulent_channel/dns.in
@@ -1,6 +1,6 @@
 512 256 144              ! itot, jtot, ktot
 6. 3. 1.                 ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 5640.                    ! visci
 poi                      ! inivel

diff --git a/examples/_manuscript_turbulent_duct/dns.in b/examples/_manuscript_turbulent_duct/dns.in
@@ -1,6 +1,6 @@
 512 128 128              ! itot, jtot, ktot
 10. 1. 1.                ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 4410.                    ! visci
 poi                      ! inivel

diff --git a/examples/closed_box/dns.in b/examples/closed_box/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 1. 1. 1.                 ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 zer                      ! inivel

diff --git a/examples/couette/dns.in b/examples/couette/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 1. 1.5 1.                ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 cou                      ! inivel

diff --git a/examples/developing_channel/dns.in b/examples/developing_channel/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 1. 1.5 1.                ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 zer                      ! inivel

diff --git a/examples/developing_duct/dns.in b/examples/developing_duct/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 1. 1.5 1.                ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 zer                      ! inivel

diff --git a/examples/half_channel/dns.in b/examples/half_channel/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 1. 1.5 1.                ! lx, ly, lz
-0.                       ! gr
+2 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 hcl                      ! inivel

diff --git a/examples/lid_driven_cavity/dns.in b/examples/lid_driven_cavity/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 1. 1. 1.                 ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 zer                      ! inivel

diff --git a/examples/periodic_channel/dns.in b/examples/periodic_channel/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 3. 1.5 1.                ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 log                      ! inivel

diff --git a/examples/periodic_duct/dns.in b/examples/periodic_duct/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 3. 1.5 1.                ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 log                      ! inivel

diff --git a/examples/taylor_green_vortex_2d/dns.in b/examples/taylor_green_vortex_2d/dns.in
@@ -1,6 +1,6 @@
 32 32 4                 ! itot, jtot, ktot
 6.2831853071795 6.283185307179586 0.125 ! lx, ly, lz
-0.                      ! gr
+1 0.                    ! gtype, gr
 .95 0.1                 ! cfl, dtmin
 100.                    ! visci
 tgw                     ! inivel

diff --git a/examples/temporal_boundary_layer/dns.in b/examples/temporal_boundary_layer/dns.in
@@ -1,6 +1,6 @@
 128 128 256                ! itot, jtot, ktot
 40. 20. 72.                ! lx, ly, lz
-4.                         ! gr
+2 4.                       ! gtype, gr
 .95 1.e5                   ! cfl, dtmin
 500.                       ! visci
 tbl                        ! inivel

diff --git a/examples/triperiodic/dns.in b/examples/triperiodic/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 1. 1. 1.                 ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 1000.                    ! visci
 zer                      ! inivel

diff --git a/examples/turbulent_channel_constant_pressure_gradient/dns.in b/examples/turbulent_channel_constant_pressure_gradient/dns.in
@@ -1,6 +1,6 @@
 512 256 144              ! itot, jtot, ktot
 12. 6. 2.                ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 180.                     ! visci
 pdc                      ! inivel

diff --git a/examples/turbulent_channel_convective_reference_frame/dns.in b/examples/turbulent_channel_convective_reference_frame/dns.in
@@ -1,6 +1,6 @@
 512 256 144              ! itot, jtot, ktot
 6. 3. 1.                 ! lx, ly, lz
-0.                       ! gr
+1 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 5640.                    ! visci
 iop                      ! inivel

diff --git a/examples/turbulent_half_channel_constant_pressure_gradient/dns.in b/examples/turbulent_half_channel_constant_pressure_gradient/dns.in
@@ -1,6 +1,6 @@
 512 256 72               ! itot, jtot, ktot
 12. 6. 1.                ! lx, ly, lz
-0.                       ! gr
+2 0.                     ! gtype, gr
 .95 1.e5                 ! cfl, dtmin
 180.                     ! visci
 hdc                      ! inivel

diff --git a/src/dns.in b/src/dns.in
@@ -1,6 +1,6 @@
 64 64 64                 ! itot, jtot, ktot
 3. 1.5 1.                ! lx, ly, lz
-1.                       ! gr
+1 1.                     ! gtype, gr
 .95 1.0e5                ! cfl, dtmin
 1000.                    ! visci
 poi                      ! inivel

diff --git a/src/initgrid.f90 b/src/initgrid.f90
@@ -11,23 +11,30 @@ module mod_initgrid
   private
   public initgrid
   contains
-  subroutine initgrid(inivel,n,gr,lz,dzc,dzf,zc,zf)
+  subroutine initgrid(gtype,n,gr,lz,dzc,dzf,zc,zf)
     !
     ! initializes the non-uniform grid along z
     !
     implicit none
-    character(len=3), intent(in) :: inivel
-    integer , intent(in ) :: n
+    integer, parameter :: CLUSTER_TWO_END              = 1, &
+                          CLUSTER_ONE_END              = 2, &
+                          CLUSTER_ONE_END_R            = 3, &
+                          CLUSTER_MIDDLE               = 4
+    integer , intent(in ) :: gtype,n
     real(rp), intent(in ) :: gr,lz
     real(rp), intent(out), dimension(0:n+1) :: dzc,dzf,zc,zf
     real(rp) :: z0
     integer :: k
     procedure (), pointer :: gridpoint => null()
-    select case(inivel)
-    case('zer','log','poi','cou','iop','pdc')
+    select case(gtype)
+    case(CLUSTER_TWO_END)
       gridpoint => gridpoint_cluster_two_end
-    case('hcl','hcp','hdc','tbl')
+    case(CLUSTER_ONE_END)
       gridpoint => gridpoint_cluster_one_end
+    case(CLUSTER_ONE_END_R)
+      gridpoint => gridpoint_cluster_one_end_r
+    case(CLUSTER_MIDDLE)
+      gridpoint => gridpoint_cluster_middle
     case default
       gridpoint => gridpoint_cluster_two_end
     end select
@@ -102,6 +109,20 @@ subroutine gridpoint_cluster_one_end(alpha,z0,z)
       z = z0
     end if
   end subroutine gridpoint_cluster_one_end
+  subroutine gridpoint_cluster_one_end_r(alpha,r0,r)
+    !
+    ! clustered at the upper side
+    !
+    implicit none
+    real(rp), intent(in ) :: alpha,r0
+    real(rp), intent(out) :: r
+    if(alpha /= 0._rp) then
+      r = 1._rp-1.0_rp*(1._rp+tanh((1._rp-r0-1.0_rp)*alpha)/tanh(alpha/1._rp))
+      !r = 1._rp-1.0_rp*(1._rp+erf( (1._rp-r0-1.0_rp)*alpha)/erf( alpha/1._rp))
+    else
+      r = r0
+    end if
+  end subroutine gridpoint_cluster_one_end_r
   subroutine gridpoint_cluster_middle(alpha,z0,z)
     !
     ! clustered in the middle

diff --git a/src/main.f90 b/src/main.f90
@@ -56,7 +56,7 @@ program cans
                                  cfl,dtmin, &
                                  inivel,    &
                                  dims, &
-                                 gr, &
+                                 gtype,gr, &
                                  is_forced,velf,bforce, &
                                  ng,l,dl,dli, &
                                  read_input
@@ -210,7 +210,7 @@ program cans
   if(myid == 0) print*, '*** Beginning of simulation ***'
   if(myid == 0) print*, '*******************************'
   if(myid == 0) print*, ''
-  call initgrid(inivel,ng(3),gr,lz,dzc_g,dzf_g,zc_g,zf_g)
+  call initgrid(gtype,ng(3),gr,lz,dzc_g,dzf_g,zc_g,zf_g)
   if(myid == 0) then
     inquire(iolength=rlen) 1._rp
     open(99,file=trim(datadir)//'grid.bin',access='direct',recl=4*ng(3)*rlen)

diff --git a/src/param.f90 b/src/param.f90
@@ -28,7 +28,9 @@ module mod_param
 ! variables to be determined from the input file 'dns.in'
 !
 integer , protected :: itot,jtot,ktot
-real(rp), protected :: lx,ly,lz,dx,dy,dz,dxi,dyi,dzi,gr
+real(rp), protected :: lx,ly,lz,dx,dy,dz,dxi,dyi,dzi
+integer , protected :: gtype
+real(rp), protected :: gr
 real(rp), protected :: cfl,dtmin
 real(rp), protected :: visci
 !
@@ -81,7 +83,7 @@ subroutine read_input(myid)
       if( ierr == 0 ) then
         read(iunit,*,iostat=ierr) itot,jtot,ktot
         read(iunit,*,iostat=ierr) lx,ly,lz
-        read(iunit,*,iostat=ierr) gr
+        read(iunit,*,iostat=ierr) gtype,gr
         read(iunit,*,iostat=ierr) cfl,dtmin
         read(iunit,*,iostat=ierr) visci
         read(iunit,*,iostat=ierr) inivel

diff --git a/src/solver.f90 b/src/solver.f90
@@ -56,6 +56,7 @@ subroutine solver(n,ng,arrplan,normfft,lambdaxy,a,b,c,bc,c_or_f,p)
     call fft(arrplan(1,2),py) ! fwd transform in y
     !
     call transpose_y_to_z(py,pz)
+    !
     q = 0
     if(c_or_f(3) == 'f'.and.bc(1,3) == 'D') q = 1
     if(bc(0,3)//bc(1,3) == 'PP') then

diff --git a/tests/lid_driven_cavity/dns.in b/tests/lid_driven_cavity/dns.in
@@ -1,6 +1,6 @@
 2 64 64                ! itot, jtot, ktot
 0.03125 1. 1.          ! lx, ly, lz
-0.                     ! gr
+1 0.                   ! gtype,gr
 .95 1.e5               ! cfl
 1000.                  ! visci
 zer                    ! inivel