From 859e1d2f6767dcd59aa36107794f62b5e8fa59e1 Mon Sep 17 00:00:00 2001 From: Pedro Costa Date: Fri, 1 Sep 2023 13:49:12 +0200 Subject: [PATCH] More logical setup of the grid clustering (#85) * added more logical tweaking for the grid clustering * changed `gtype` enumerator * minor documentation fix --- docs/INFO_INPUT.md | 41 +++++++++++-------- examples/_manuscript_lid_driven_cavity/dns.in | 2 +- .../_manuscript_taylor_green_vortex/dns.in | 2 +- examples/_manuscript_turbulent_channel/dns.in | 2 +- examples/_manuscript_turbulent_duct/dns.in | 2 +- examples/closed_box/dns.in | 2 +- examples/couette/dns.in | 2 +- examples/developing_channel/dns.in | 2 +- examples/developing_duct/dns.in | 2 +- examples/half_channel/dns.in | 2 +- examples/lid_driven_cavity/dns.in | 2 +- examples/periodic_channel/dns.in | 2 +- examples/periodic_duct/dns.in | 2 +- examples/taylor_green_vortex_2d/dns.in | 2 +- examples/temporal_boundary_layer/dns.in | 2 +- examples/triperiodic/dns.in | 2 +- .../dns.in | 2 +- .../dns.in | 2 +- .../dns.in | 2 +- src/dns.in | 2 +- src/initgrid.f90 | 33 ++++++++++++--- src/main.f90 | 4 +- src/param.f90 | 6 ++- src/solver.f90 | 1 + tests/lid_driven_cavity/dns.in | 2 +- 25 files changed, 77 insertions(+), 48 deletions(-) diff --git a/docs/INFO_INPUT.md b/docs/INFO_INPUT.md index 572e564b..b8a04fc1 100644 --- 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 index bcc145f1..9c469b71 100644 --- 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 index 1cc7f61c..d4eb4fb9 100644 --- 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 index ae08aab4..cd505b3d 100644 --- 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 index af7fe596..7ad5422e 100644 --- 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 index 7cf58514..9f9b2f37 100644 --- 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 index 0aa58ff0..bdc4c6da 100644 --- 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 index e1066301..761b3870 100644 --- 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 index 769490c5..7fe49453 100644 --- 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 index 05d5cb82..5a428cc2 100644 --- 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 index e2c2fbc7..a15fe10d 100644 --- 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 index 24c92cad..612017a7 100644 --- 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 index 3133c8b7..3babc506 100644 --- 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 index 9d82beb1..962fad83 100644 --- 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 index 57a43e52..476071db 100644 --- 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 index dc5b3d71..74363edd 100644 --- 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 index b1134be5..73605306 100644 --- 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 index 3cf3d992..8758e724 100644 --- 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 index ce107884..6b5d575a 100644 --- 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 index b714ca3e..e90d54ef 100644 --- 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 index 09553232..4774ac91 100644 --- 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 index c04e59fb..a8c33324 100644 --- 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 index 07dc2937..0251ec61 100644 --- 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 index a0274e51..ba37c1ce 100644 --- 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 index 99f0d80a..45a7a3bb 100644 --- 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