From b0dc06aa7cd3c0de1e287def43cb05622339d120 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Fri, 15 Mar 2019 07:01:39 +0100 Subject: [PATCH 01/32] Pointer for grid was missing and grid % vol causing an error. --- Sources/Process/Turbulence/Source_Zeta_K_Eps_Zeta_F.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Sources/Process/Turbulence/Source_Zeta_K_Eps_Zeta_F.f90 b/Sources/Process/Turbulence/Source_Zeta_K_Eps_Zeta_F.f90 index 5f99fbcd4..c16e56f03 100644 --- a/Sources/Process/Turbulence/Source_Zeta_K_Eps_Zeta_F.f90 +++ b/Sources/Process/Turbulence/Source_Zeta_K_Eps_Zeta_F.f90 @@ -15,10 +15,10 @@ subroutine Source_Zeta_K_Eps_Zeta_F(grid, sol, n_step) !------------------------------------------------------------------------------! implicit none !--------------------------------[Arguments]-----------------------------------! - type(Grid_Type) :: grid type(Solver_Type), target :: sol integer :: n_step !----------------------------------[Locals]------------------------------------! + type(Grid_Type), pointer :: grid type(Matrix_Type), pointer :: a real, pointer :: b(:) integer :: c @@ -50,7 +50,7 @@ subroutine Source_Zeta_K_Eps_Zeta_F(grid, sol, n_step) if(n_step > 500) then b(c) = b(c) + f22 % n(c) * grid % vol(c) * density else - b(c) = b(c) + max(0.0, f22 % n(c)*grid % vol(c)) * density + b(c) = b(c) + max(0.0, f22 % n(c) * grid % vol(c)) * density a % val(a % dia(c)) = a % val(a % dia(c)) & + max(0.0, -f22 % n(c) * grid % vol(c) & / (zeta % n(c) + TINY)) * density From a64dc2023600f26e241b5548c675a3393b231faf Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Fri, 15 Mar 2019 14:42:27 +0100 Subject: [PATCH 02/32] I added heat_flux calculation needed by subroutine. Previously it was done in Update_Boundary_Condition.f90. --- Tests/Les/Channel_Re_Tau_180/User_Mod/Save_Results.f90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Tests/Les/Channel_Re_Tau_180/User_Mod/Save_Results.f90 b/Tests/Les/Channel_Re_Tau_180/User_Mod/Save_Results.f90 index e39593765..8788e19cd 100644 --- a/Tests/Les/Channel_Re_Tau_180/User_Mod/Save_Results.f90 +++ b/Tests/Les/Channel_Re_Tau_180/User_Mod/Save_Results.f90 @@ -83,6 +83,13 @@ subroutine User_Mod_Save_Results(flow, save_name) nu_max = 0.0 n_points = 0 + if(heat_transfer) then + call Comm_Mod_Global_Sum_Real(heat_flux) + call Comm_Mod_Global_Sum_Real(heated_area) + heat_flux = heat_flux / (heated_area + TINY) + heat = heat_flux * heated_area + end if + open(9, file=coord_name) ! Write the number of searching intervals From a6554531c5d85db1d16626efb9ec3a720dabe48d Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Fri, 15 Mar 2019 14:45:33 +0100 Subject: [PATCH 03/32] I added heat_flux calculation, needed by the subroutine. Previously it was calculated in Update_Boundary_Condition. --- Tests/Les/Pipe_Re_Tau_180/User_Mod/Save_Results.f90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Tests/Les/Pipe_Re_Tau_180/User_Mod/Save_Results.f90 b/Tests/Les/Pipe_Re_Tau_180/User_Mod/Save_Results.f90 index a495ff3a1..62e654ffe 100644 --- a/Tests/Les/Pipe_Re_Tau_180/User_Mod/Save_Results.f90 +++ b/Tests/Les/Pipe_Re_Tau_180/User_Mod/Save_Results.f90 @@ -84,6 +84,13 @@ subroutine User_Mod_Save_Results(flow, save_name) nu_max = 0.0 n_points = 0 + if(heat_transfer) then + call Comm_Mod_Global_Sum_Real(heat_flux) + call Comm_Mod_Global_Sum_Real(heated_area) + heat_flux = heat_flux / (heated_area + TINY) + heat = heat_flux * heated_area + end if + open(9, file=coord_name) ! Write the number of searching intervals From d56807055ffa62ec49596dc9c32e8d0c29e3bb19 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sun, 17 Mar 2019 08:02:04 +0100 Subject: [PATCH 04/32] t2 transport equation for zeta-f model is moved after zeta and F22 equations. This solved the problem that occurs otherwise. It is still important find out what is really happening. --- Sources/Process/Main_Pro.f90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/Sources/Process/Main_Pro.f90 b/Sources/Process/Main_Pro.f90 index a2d7194c9..834cc393e 100644 --- a/Sources/Process/Main_Pro.f90 +++ b/Sources/Process/Main_Pro.f90 @@ -270,11 +270,6 @@ program Processor call Compute_Turbulent(flow, sol, dt, ini, kin, n) call Compute_Turbulent(flow, sol, dt, ini, eps, n) - if(heat_transfer) then - call Calculate_Heat_Flux(flow) - call Compute_Turbulent(flow, sol, dt, ini, t2, n) - end if - call Update_Boundary_Values(flow) call Compute_F22(grid, sol, ini, f22) @@ -282,6 +277,10 @@ program Processor call Calculate_Vis_T_K_Eps_Zeta_F(flow) + if(heat_transfer) then + call Calculate_Heat_Flux(flow) + call Compute_Turbulent(flow, sol, dt, ini, t2, n) + end if end if if(turbulence_model .eq. RSM_MANCEAU_HANJALIC .or. & From 94eff2e6ed237f18922270cb3b288b079f65a7f9 Mon Sep 17 00:00:00 2001 From: Egor Palkin Date: Mon, 18 Mar 2019 16:59:13 +0700 Subject: [PATCH 05/32] This commit is designed to fix "Problem with Source_T2 subroutine #128" reported by Muhamed. It was caused by missing t2 and p_t2 in heat_transfer expressions, when turbulence_model .eq. HYBRID_LES_RANS. --- Sources/Process/Backup_Mod/Load.f90 | 5 ++--- Sources/Process/Backup_Mod/Save.f90 | 5 ++--- Sources/Process/Initialize_Variables.f90 | 7 ++++--- Sources/Process/Load_Boundary_Conditions.f90 | 15 ++++++++++++--- Sources/Process/Main_Pro.f90 | 9 +++++---- Sources/Process/Save_Cgns_Results.f90 | 3 ++- Sources/Process/Save_Vtu_Results.f90 | 6 +++--- Sources/Process/Turbulence/Allocate.f90 | 3 +++ Sources/Process/Update_Boundary_Values.f90 | 7 ++++--- 9 files changed, 37 insertions(+), 23 deletions(-) diff --git a/Sources/Process/Backup_Mod/Load.f90 b/Sources/Process/Backup_Mod/Load.f90 index 097dee0a7..aebf3d1b7 100644 --- a/Sources/Process/Backup_Mod/Load.f90 +++ b/Sources/Process/Backup_Mod/Load.f90 @@ -184,12 +184,11 @@ subroutine Backup_Mod_Load(fld, time_step, time_step_stat, backup) ! Turbulence quantities connected with heat transfer end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then call Backup_Mod_Read_Variable(fh, d, vc, 't2', t2) call Backup_Mod_Read_Cell_Bnd(fh, d, vc, 'p_t2', p_t2 (-nb_s:nc_s)) call Backup_Mod_Read_Cell_Bnd(fh, d, vc, 'con_wall', con_wall(-nb_s:nc_s)) - else if (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) then - call Backup_Mod_Read_Cell_Bnd(fh, d, vc, 'con_wall', con_wall(-nb_s:nc_s)) end if !----------------------------! diff --git a/Sources/Process/Backup_Mod/Save.f90 b/Sources/Process/Backup_Mod/Save.f90 index 0b6cb9351..0df439042 100644 --- a/Sources/Process/Backup_Mod/Save.f90 +++ b/Sources/Process/Backup_Mod/Save.f90 @@ -163,12 +163,11 @@ subroutine Backup_Mod_Save(fld, time_step, time_step_stat, name_save) call Backup_Mod_Write_Cell_Bnd(fh, d, vc, 'l_scale', l_scale(-nb_s:nc_s)) end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then call Backup_Mod_Write_Variable(fh, d, vc, 't2', t2) call Backup_Mod_Write_Cell_Bnd(fh, d, vc, 'p_t2', p_t2 (-nb_s:nc_s)) call Backup_Mod_Write_Cell_Bnd(fh, d, vc, 'con_wall', con_wall(-nb_s:nc_s)) - else if (heat_transfer .and. turbulence_model .eq. HYBRID_LES_RANS) then - call Backup_Mod_Write_Cell_Bnd(fh, d, vc, 'con_wall', con_wall(-nb_s:nc_s)) end if diff --git a/Sources/Process/Initialize_Variables.f90 b/Sources/Process/Initialize_Variables.f90 index d9fa9ded9..45a2d5e6d 100644 --- a/Sources/Process/Initialize_Variables.f90 +++ b/Sources/Process/Initialize_Variables.f90 @@ -158,7 +158,8 @@ subroutine Initialize_Variables(flow) i=Key_Ind('F22', keys,nks);prof(k,0)=f22_def; f22 %n(c)=prof(k,i) end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then i=Key_Ind('T2', keys,nks);prof(k,0)=t2_def; t2 %n(c)=prof(k,i) end if @@ -309,8 +310,8 @@ subroutine Initialize_Variables(flow) y_plus(c) = 0.001 end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. & - heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then vals(0) = t2_def; t2 % n(c) = vals(Key_Ind('T2', keys, nks)) t2 % o(c) = t2 % n(c) t2 % oo(c) = t2 % n(c) diff --git a/Sources/Process/Load_Boundary_Conditions.f90 b/Sources/Process/Load_Boundary_Conditions.f90 index c4711f56c..2315eb6b3 100644 --- a/Sources/Process/Load_Boundary_Conditions.f90 +++ b/Sources/Process/Load_Boundary_Conditions.f90 @@ -281,7 +281,10 @@ subroutine Load_Boundary_Conditions(flow, backup) i = Key_Ind('F22', keys, nks); if(i > 0) f22 % n(c) = vals(i) end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. & + heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. & + heat_transfer) ) then i = Key_Ind('T2', keys, nks); if(i > 0) t2 % n(c) = vals(i) end if @@ -420,7 +423,10 @@ subroutine Load_Boundary_Conditions(flow, backup) i = Key_Ind('F22', keys, nks); if(i>0) f22 % n(c) = prof(k,i) end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. & + heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. & + heat_transfer) ) then i = Key_Ind('T2', keys, nks); if(i>0) t2 % n(c) = prof(k,i) end if @@ -640,7 +646,10 @@ subroutine Load_Boundary_Conditions(flow, backup) if(i > 0) f22 % n(c) = wi*prof(m,i) + (1.-wi)*prof(m+1,i) end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. & + heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. & + heat_transfer) ) then i = Key_Ind('T2',keys,nks) if(i > 0) t2 % n(c) = wi*prof(m,i) + (1.-wi)*prof(m+1,i) end if diff --git a/Sources/Process/Main_Pro.f90 b/Sources/Process/Main_Pro.f90 index 834cc393e..a2d7194c9 100644 --- a/Sources/Process/Main_Pro.f90 +++ b/Sources/Process/Main_Pro.f90 @@ -270,6 +270,11 @@ program Processor call Compute_Turbulent(flow, sol, dt, ini, kin, n) call Compute_Turbulent(flow, sol, dt, ini, eps, n) + if(heat_transfer) then + call Calculate_Heat_Flux(flow) + call Compute_Turbulent(flow, sol, dt, ini, t2, n) + end if + call Update_Boundary_Values(flow) call Compute_F22(grid, sol, ini, f22) @@ -277,10 +282,6 @@ program Processor call Calculate_Vis_T_K_Eps_Zeta_F(flow) - if(heat_transfer) then - call Calculate_Heat_Flux(flow) - call Compute_Turbulent(flow, sol, dt, ini, t2, n) - end if end if if(turbulence_model .eq. RSM_MANCEAU_HANJALIC .or. & diff --git a/Sources/Process/Save_Cgns_Results.f90 b/Sources/Process/Save_Cgns_Results.f90 index 10f9df08d..bae578a4c 100644 --- a/Sources/Process/Save_Cgns_Results.f90 +++ b/Sources/Process/Save_Cgns_Results.f90 @@ -186,7 +186,8 @@ subroutine Save_Results(flow, name_save) f22 % n(1), "TurbulentQuantityF22") end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then call Cgns_Mod_Write_Field(base, block, solution, field, grid, & t2 % n(1), "TurbulentQuantityT2") call Cgns_Mod_Write_Field(base, block, solution, field, grid, & diff --git a/Sources/Process/Save_Vtu_Results.f90 b/Sources/Process/Save_Vtu_Results.f90 index c7793303b..04591851e 100644 --- a/Sources/Process/Save_Vtu_Results.f90 +++ b/Sources/Process/Save_Vtu_Results.f90 @@ -263,10 +263,10 @@ subroutine Save_Results(flow, name_save) call Save_Vtu_Scalar(grid, IN_4, IN_5, "TurbulentQuantityF22", f22 % n(1)) end if - if (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then call Save_Vtu_Scalar(grid, IN_4, IN_5, "TurbulentQuantityT2", t2 % n(1)) - call Save_Vtu_Scalar(grid, IN_4, IN_5, "TurbulentT2Production", & - p_t2(1)) + call Save_Vtu_Scalar(grid, IN_4, IN_5, "TurbulentT2Production", p_t2(1)) end if ! Save vis and vis_t if(turbulence_model .eq. DES_SPALART .or. & diff --git a/Sources/Process/Turbulence/Allocate.f90 b/Sources/Process/Turbulence/Allocate.f90 index d8f6b96a5..169b4683e 100644 --- a/Sources/Process/Turbulence/Allocate.f90 +++ b/Sources/Process/Turbulence/Allocate.f90 @@ -565,10 +565,12 @@ subroutine Turbulence_Allocate(flow) if(heat_transfer) then + call Var_Mod_Allocate_Solution('T2', '', t2, grid) call Var_Mod_Allocate_New_Only('UT', ut, grid) call Var_Mod_Allocate_New_Only('VT', vt, grid) call Var_Mod_Allocate_New_Only('WT', wt, grid) allocate(con_wall(-grid % n_bnd_cells:grid % n_cells)); con_wall = 0. + allocate(p_t2 (-grid % n_bnd_cells:grid % n_cells)); p_t2 = 0. end if ! heat_transfer @@ -603,6 +605,7 @@ subroutine Turbulence_Allocate(flow) call Var_Mod_Allocate_Statistics(wt) ! new value allocated above call Var_Mod_Allocate_New_Only('TT', tt, grid) call Var_Mod_Allocate_Statistics(tt) + call Var_Mod_Allocate_Statistics(t2) end if ! heat_transfer diff --git a/Sources/Process/Update_Boundary_Values.f90 b/Sources/Process/Update_Boundary_Values.f90 index b7a6e9d3d..a30b4c113 100644 --- a/Sources/Process/Update_Boundary_Values.f90 +++ b/Sources/Process/Update_Boundary_Values.f90 @@ -113,7 +113,8 @@ subroutine Update_Boundary_Values(flow) end if end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. OUTFLOW .or. & Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. CONVECT .or. & Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. PRESSURE .or. & @@ -236,8 +237,8 @@ subroutine Update_Boundary_Values(flow) f22 % n(c2) = f22 % n(grid % bnd_cond % copy_c(c2)) end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. & - heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then t2 % n(c2) = t2 % n(grid % bnd_cond % copy_c(c2)) end if From fdb61a93a64766fd7bee3e0bfaea33d019b49b85 Mon Sep 17 00:00:00 2001 From: Egor Palkin Date: Mon, 18 Mar 2019 17:11:54 +0700 Subject: [PATCH 06/32] This commit is designed to fix "Problem with Source_T2 subroutine #128" reported by Muhamed. It was caused by missing t2 and p_t2 in heat_transfer expressions, when turbulence_model .eq. HYBRID_LES_RANS. (a minor fix to the previous commit) --- Sources/Process/Turbulence/Compute_Turbulent.f90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Sources/Process/Turbulence/Compute_Turbulent.f90 b/Sources/Process/Turbulence/Compute_Turbulent.f90 index dac481f42..228c06069 100644 --- a/Sources/Process/Turbulence/Compute_Turbulent.f90 +++ b/Sources/Process/Turbulence/Compute_Turbulent.f90 @@ -212,7 +212,8 @@ subroutine Compute_Turbulent(flow, sol, dt, ini, phi, n_step) if(phi % name .eq. 'ZETA') call Source_Zeta_K_Eps_Zeta_F(flow, sol, n_step) end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then if(phi % name .eq. 'T2') call Source_T2(flow, sol) end if @@ -242,7 +243,7 @@ subroutine Compute_Turbulent(flow, sol, dt, ini, phi, n_step) do c = 1, grid % n_cells if( phi % n(c) < 0.0 ) phi % n(c) = phi % o(c) - if(phi % name .eq. 'ZETA') phi % n(c) = min(phi % n(c), 1.8) + if(phi % name .eq. 'ZETA') phi % n(c) = min(phi % n(c), 1.8) end do ! Print info on the screen @@ -257,7 +258,8 @@ subroutine Compute_Turbulent(flow, sol, dt, ini, phi, n_step) call Info_Mod_Iter_Fill_At(3, 3, phi % name, exec_iter, phi % res) end if - if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then + if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. & + (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then if(phi % name .eq. 'T2') & call Info_Mod_Iter_Fill_At(3, 5, phi % name, exec_iter, phi % res) end if From beca7235e424e4f659092d0a24e5ddaaac65f1f4 Mon Sep 17 00:00:00 2001 From: Egor Palkin Date: Mon, 18 Mar 2019 17:47:00 +0700 Subject: [PATCH 07/32] Improved test_build.sh and fixed a couple of issues. --- Tests/test_build.sh | 106 +++++++++++++++++++++++++++++++------------- 1 file changed, 75 insertions(+), 31 deletions(-) diff --git a/Tests/test_build.sh b/Tests/test_build.sh index 3439138e4..cd3ae9bcd 100755 --- a/Tests/test_build.sh +++ b/Tests/test_build.sh @@ -1,16 +1,25 @@ #!/bin/bash -# it is an useful script in debug perposus to automatically build and run most -# cases in T-Flows +# Description: This script is made for debug purposes. +# It automatically builds and runs most cases in T-Flows -# compilation flags +# Requires: gnuplot, texlive-base, mpi, gfortran + +# Compilation flags used in makefiles FCOMP="gnu" -# conduct tests with DEBUG=yes +# Conduct tests with DEBUG=yes DEBUG="no" -# repeat tests with CGNS_HDF5=yes +# Repeat tests with CGNS_HDF5=yes CGNS="yes" -# folder structure +# A small reminder how to set up alternatives if you have mpich and openmpi: +#update-alternatives --install /usr/bin/mpif90 mpif90 /usr/bin/mpif90.openmpi 20 +#update-alternatives --install /usr/bin/mpif90 mpif90 /usr/bin/mpif90.mpich 80 + +#update-alternatives --install /usr/bin/mpirun mpirun /usr/bin/mpirun.openmpi 20 +#update-alternatives --install /usr/bin/mpirun mpirun /usr/bin/mpirun.mpich 80 + +# Folder structure TEST_DIR=$PWD # dir with tests GENE_DIR=$PWD/../Sources/Generate # Generate src folder CONV_DIR=$PWD/../Sources/Convert # Convert src folder @@ -18,14 +27,14 @@ DIVI_DIR=$PWD/../Sources/Divide # Divide src folder PROC_DIR=$PWD/../Sources/Process # Process src folder BINA_DIR=$PWD/../Binaries/ # binaries folder -# executables +# Executables GENE_EXE=$BINA_DIR/Generate # Generate ex CONV_EXE=$BINA_DIR/Convert # Convert ex DIVI_EXE=$BINA_DIR/Divide # Divide ex PROC_EXE=$BINA_DIR/Process # Process ex -# folders with geometry -# generator file (.dom) +# Folders with geometry +# Generator file (.dom) LAMINAR_BACKSTEP_ORTH_DIR=$TEST_DIR/Laminar/Backstep/Orthogonal LAMINAR_BACKSTEP_NON_ORTH_DIR=$TEST_DIR/Laminar/Backstep/Nonorthogonal @@ -33,7 +42,7 @@ LES_CAVITY_LID_DRIVEN_DIR=$TEST_DIR/Laminar/Cavity/Lid_Driven/Re_1000 LES_CAVITY_THERM_DRIVEN_DIR_106=$TEST_DIR/Laminar/Cavity/Thermally_Driven/Ra_10e6 LES_CAVITY_THERM_DRIVEN_DIR_108=$TEST_DIR/Laminar/Cavity/Thermally_Driven/Ra_10e8 -RANS_BACKSTEP_5100_DIR=$TEST_DIR/Rans/Backstep_Re_5100 +RANS_BACKSTEP_5100_DIR=$TEST_DIR/Rans/Backstep_Re_05100 RANS_BACKSTEP_28000_DIR=$TEST_DIR/Rans/Backstep_Re_28000 RANS_CHANNEL_LR_LONG_DIR=$TEST_DIR/Rans/Channel_Re_Tau_590/Long_Domain @@ -46,7 +55,7 @@ $TEST_DIR/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh HYB_CHANNEL_HR_STRETCHED_DIR=\ $TEST_DIR/Hybrid_Les_Rans/Channel_Re_Tau_2000/Stretched_Mesh -# mesh file (.cgns/.neu) +# Mesh file (.cgns/.neu) LES_PIPE_DIR=$TEST_DIR/Les/Pipe_Re_Tau_180 RANS_IMPINGING_JET_DIR=$TEST_DIR/Rans/Impinging_Jet_2d_Distant_Re_23000 @@ -54,10 +63,10 @@ RANS_FUEL_BUNDLE_DIR=$TEST_DIR/Rans/Fuel_Bundle #RANS_PIPE_DIR=$TEST_DIR/Rans/Pipe_Re_Tau_2000 #------------------------------------------------------------------------------# -# start time measurements from this moment +# Start time measurements from this moment current_time=$(date +%s) -# script logs +# Script logs FULL_LOG=$TEST_DIR/test_build.log # logs of current script if [ -f $FULL_LOG ]; then cp /dev/null $FULL_LOG; fi echo "Full log is being written in file" "$FULL_LOG" @@ -65,10 +74,10 @@ echo "Full log is being written in file" "$FULL_LOG" # exit when any command fails set -e -# keep track of the last executed command -trap 'last_command=$current_command; current_command=$BASH_COMMAND' DEBUG +# Keep track of the last executed command +# trap 'last_command=$current_command; current_command=$BASH_COMMAND' DEBUG # echo an error message before exiting -trap 'echo "\"${last_command}\" command filed with exit code $?."' EXIT +# trap 'echo "\"${last_command}\" command filed with exit code $?."' EXIT #------------------------------------------------------------------------------# # time in seconds @@ -464,6 +473,7 @@ function processor_backup_tests { # Grasp/embrace as many different model combinations as you can + echo "================================= TEST 1 ==============================" #-- Channel_Re_Tau_590 [k_eps model + T] replace_line_with_first_occurence_in_file "TURBULENCE_MODEL" \ "TURBULENCE_MODEL k_eps" $RANS_CHANNEL_LR_UNIFORM_DIR/control @@ -472,6 +482,7 @@ function processor_backup_tests { process_backup_test yes $RANS_CHANNEL_LR_UNIFORM_DIR fi + echo "================================= TEST 2 ==============================" #-- Channel_Re_Tau_590 [k_eps_zeta_f model + T] replace_line_with_first_occurence_in_file "TURBULENCE_MODEL" \ "TURBULENCE_MODEL k_eps_zeta_f" $RANS_CHANNEL_LR_UNIFORM_DIR/control @@ -480,6 +491,7 @@ function processor_backup_tests { process_backup_test yes $RANS_CHANNEL_LR_UNIFORM_DIR fi + echo "================================= TEST 3 ==============================" #-- Channel_Re_Tau_590_Rsm [rsm_hanjalic_jakirlic model + T] replace_line_with_first_occurence_in_file "TURBULENCE_MODEL" \ "TURBULENCE_MODEL rsm_hanjalic_jakirlic" $RANS_CHANNEL_LR_RSM_DIR/control @@ -488,6 +500,7 @@ function processor_backup_tests { process_backup_test yes $RANS_CHANNEL_LR_RSM_DIR fi + echo "================================= TEST 4 ==============================" #-- Channel_Re_Tau_590_Rsm [rsm_manceau_hanjalic model + T] replace_line_with_first_occurence_in_file "TURBULENCE_MODEL" \ "TURBULENCE_MODEL rsm_manceau_hanjalic" $RANS_CHANNEL_LR_RSM_DIR/control @@ -496,12 +509,14 @@ function processor_backup_tests { process_backup_test yes $RANS_CHANNEL_LR_RSM_DIR fi + echo "================================= TEST 5 ==============================" #-- Pipe_Re_Tau_180 [les_dynamic] process_backup_test no $LES_PIPE_DIR if [ "$CGNS" = "yes" ]; then process_backup_test yes $LES_PIPE_DIR fi + echo "================================= TEST 6 ==============================" #-- Cavity_Lid_Driven_Re_1000 [none] process_backup_test no $LES_CAVITY_LID_DRIVEN_DIR if [ "$CGNS" = "yes" ]; then @@ -526,10 +541,19 @@ function process_save_exit_now_test { echo "Test: save_now & exit now on " $2 - name_in_div=$(head -n1 "$2"/divide.scr) - nproc_in_div=$(head -n2 "$2"/divide.scr | tail -n1) + cd "$2" + name_in_div=$(head -n1 divide.scr) + nproc_in_div=$(head -n2 divide.scr | tail -n1) - #----------------------------------------# + # change number of timesteps to 3 + replace_line_with_first_occurence_in_file "NUMBER_OF_TIME_STEPS" \ + "NUMBER_OF_TIME_STEPS 3" control + + # change backup interval to 1 ts + replace_line_with_first_occurence_in_file "BACKUP_SAVE_INTERVAL" \ + "BACKUP_SAVE_INTERVAL 1" control + + echo "================================= TEST 1 ==============================" echo "np=1, MPI=no" clean_compile $PROC_DIR $1 no # dir CGNS_HDF5 MPI cd $2 @@ -544,7 +568,7 @@ function process_save_exit_now_test { echo "save_exit_now_test was successfull" fi - #----------------------------------------# + echo "================================= TEST 2 ==============================" echo "np=1, MPI=yes" clean_compile $PROC_DIR $1 yes # dir CGNS_HDF5 MPI cd $2 @@ -559,7 +583,7 @@ function process_save_exit_now_test { echo "save_exit_now_test was successfull" fi - #----------------------------------------# + echo "================================= TEST 3 ==============================" echo "np=2, MPI=yes" echo "save_now" @@ -583,9 +607,9 @@ function process_save_exit_now_tests { echo " !!" echo " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - process_save_exit_now_test no $RANS_CHANNEL_DIR + process_save_exit_now_test no $LAMINAR_BACKSTEP_ORTH_DIR if [ "$CGNS" = "yes" ]; then - process_save_exit_now_test yes $RANS_CHANNEL_DIR + process_save_exit_now_test yes $LAMINAR_BACKSTEP_ORTH_DIR fi } @@ -677,62 +701,71 @@ function processor_compilation_tests { echo " !!" echo " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" + echo "================================= TEST 1 ==============================" # no User_Mod/ dir !!! processor_compilation_test \ "$LES_CAVITY_LID_DRIVEN_DIR" \ "none" \ "$LES_CAVITY_LID_DRIVEN_DIR/Xmgrace" + echo "================================= TEST 2 ==============================" # no User_Mod/ dir !!! processor_compilation_test \ "$LES_CAVITY_THERM_DRIVEN_DIR_106" \ "none" \ "$LES_CAVITY_THERM_DRIVEN_DIR_106/Xmgrace" + echo "================================= TEST 3 ==============================" # no User_Mod/ dir !!! processor_compilation_test \ "$LES_CAVITY_THERM_DRIVEN_DIR_108" \ "none" \ "$LES_CAVITY_THERM_DRIVEN_DIR_108/Xmgrace" - # [~2 min test] + echo "================================= TEST 4 ==============================" + # User_Mod/ dir exists processor_compilation_test \ "$RANS_CHANNEL_LR_UNIFORM_DIR" \ "k_eps" \ "$RANS_CHANNEL_LR_UNIFORM_DIR/Xmgrace" - # [~2 min test] + echo "================================= TEST 5 ==============================" + # User_Mod/ dir exists processor_compilation_test \ "$RANS_CHANNEL_LR_STRETCHED_DIR" \ "k_eps_zeta_f" \ "$RANS_CHANNEL_LR_STRETCHED_DIR/Xmgrace" - # [~5 min test] + echo "================================= TEST 6 ==============================" + # User_Mod/ dir exists processor_compilation_test \ "$RANS_CHANNEL_LR_RSM_DIR" \ "rsm_manceau_hanjalic" \ "$RANS_CHANNEL_LR_RSM_DIR/Xmgrace" - # [~5 min test] + echo "================================= TEST 7 ==============================" + # User_Mod/ dir exists processor_compilation_test \ "$RANS_CHANNEL_LR_RSM_DIR" \ "rsm_hanjalic_jakirlic" \ "$RANS_CHANNEL_LR_RSM_DIR/Xmgrace" - # [~12 HOURS test] + echo "================================= TEST 8 ==============================" + # User_Mod/ dir exists processor_compilation_test \ "$HYB_CHANNEL_HR_UNIFORM_DIR" \ "hybrid_les_rans" \ "$HYB_CHANNEL_HR_UNIFORM_DIR/Xmgrace" - # [~12 HOURS test] + echo "================================= TEST 9 ==============================" + # User_Mod/ dir exists processor_compilation_test \ "$HYB_CHANNEL_HR_STRETCHED_DIR" \ "hybrid_les_rans" \ "$HYB_CHANNEL_HR_STRETCHED_DIR/Xmgrace" # # Issue: pipe does not pass processor_backup_tests -# processor_full_length_test \ +# processor_compilation_test \ # "$LES_PIPE_DIR" \ # "les_dynamic" \ # "$LES_PIPE_DIR/Xmgrace" @@ -794,7 +827,7 @@ function processor_full_length_test { launch_process par $nproc_in_div - if [ -f *-res-plus.dat ];then + if ls *-res-plus.dat 1> /dev/null 2>&1; then # case-ts??????-res-plus.dat # extract essential data from produced .dat files last_results_plus_dat_file=$(realpath --relative-to="$3" \ @@ -805,6 +838,8 @@ function processor_full_length_test { launch_gnuplot "$3" gnuplot_script_template.sh \ "$last_results_plus_dat_file" "result_plus_"$2"" + else + echo "Warning: file *-res-plus.dat does not exist" fi } #------------------------------------------------------------------------------# @@ -822,54 +857,63 @@ function processor_full_length_tests { echo " !!" echo " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" + echo "================================= TEST 1 ==============================" # no User_Mod/ dir !!! processor_full_length_test \ "$LES_CAVITY_LID_DRIVEN_DIR" \ "none" \ "$LES_CAVITY_LID_DRIVEN_DIR/Xmgrace" + echo "================================= TEST 2 ==============================" # no User_Mod/ dir !!! processor_full_length_test \ "$LES_CAVITY_THERM_DRIVEN_DIR_106" \ "none" \ "$LES_CAVITY_THERM_DRIVEN_DIR_106/Xmgrace" + echo "================================= TEST 3 ==============================" # no User_Mod/ dir !!! processor_full_length_test \ "$LES_CAVITY_THERM_DRIVEN_DIR_108" \ "none" \ "$LES_CAVITY_THERM_DRIVEN_DIR_108/Xmgrace" + echo "================================= TEST 4 ==============================" # [~2 min test] processor_full_length_test \ "$RANS_CHANNEL_LR_UNIFORM_DIR" \ "k_eps" \ "$RANS_CHANNEL_LR_UNIFORM_DIR/Xmgrace" + echo "================================= TEST 5 ==============================" # [~2 min test] processor_full_length_test \ "$RANS_CHANNEL_LR_STRETCHED_DIR" \ "k_eps_zeta_f" \ "$RANS_CHANNEL_LR_STRETCHED_DIR/Xmgrace" + echo "================================= TEST 6 ==============================" # [~5 min test] processor_full_length_test \ "$RANS_CHANNEL_LR_RSM_DIR" \ "rsm_manceau_hanjalic" \ "$RANS_CHANNEL_LR_RSM_DIR/Xmgrace" + echo "================================= TEST 7 ==============================" # [~5 min test] processor_full_length_test \ "$RANS_CHANNEL_LR_RSM_DIR" \ "rsm_hanjalic_jakirlic" \ "$RANS_CHANNEL_LR_RSM_DIR/Xmgrace" + echo "================================= TEST 8 ==============================" # [~12 HOURS test] processor_full_length_test \ "$HYB_CHANNEL_HR_UNIFORM_DIR" \ "hybrid_les_rans" \ "$HYB_CHANNEL_HR_UNIFORM_DIR/Xmgrace" + echo "================================= TEST 9 ==============================" # [~12 HOURS test] processor_full_length_test \ "$HYB_CHANNEL_HR_STRETCHED_DIR" \ From 494500fc3c1819be91145a7e6268234916799918 Mon Sep 17 00:00:00 2001 From: Egor Palkin Date: Tue, 19 Mar 2019 17:35:13 +0700 Subject: [PATCH 08/32] Updated Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh --- .../Xmgrace/gnuplot_script_template.sh | 8 +-- .../Xmgrace/k_plus_Re_tau_2000.dat | 65 +++++++++++++++++++ Tests/test_build.sh | 2 +- 3 files changed, 70 insertions(+), 5 deletions(-) create mode 100644 Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/k_plus_Re_tau_2000.dat diff --git a/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh index 05f80417b..578856866 100644 --- a/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh +++ b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh @@ -103,10 +103,10 @@ set ytics offset 0,0.0 border 0,1,7 scale 5 set mytics 2 plot \ -1/0 with lp ls 1 ps 5 t '$k_{res}$', \ -'DAT_FILE_WITH_RESULTS_MACRO' usi ($1):($3) with lp ls 1 not, \ -1/0 with lp ls 2 ps 5 t '$k_{tot}$', \ -'DAT_FILE_WITH_RESULTS_MACRO' usi ($1):($5) with lp ls 2 not +1/0 with lp ls 1 ps 5 t '$k_{tot} \, Current$', \ +'DAT_FILE_WITH_RESULTS_MACRO' usi ($1):($5) with lp ls 1 not, \ +1/0 with lp ls 2 ps 5 t '$k_{tot} \, ref. LES$', \ +'k_plus_Re_tau_2000.dat' usi ($1):($2) with lp ls 2 not # necessary line unset multiplot \ No newline at end of file diff --git a/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/k_plus_Re_tau_2000.dat b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/k_plus_Re_tau_2000.dat new file mode 100644 index 000000000..f1569f946 --- /dev/null +++ b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/k_plus_Re_tau_2000.dat @@ -0,0 +1,65 @@ + 0.76 0.05341118800000001 + 2.32 0.46102494 + 4.02 1.2829949 + 5.88 2.4515291 + 7.9 3.7440052 + 10.1 4.8625406 + 12.48 5.6116377 + 15.08 5.9721373 + 17.92 6.0368282 + 20.98 5.9186594 + 24.34 5.7072191 + 27.98 5.4597938 + 31.92 5.2073395 + 36.22 4.966163 +40.90000000000001 4.744191499999999 + 45.98 4.5423094 + 51.48 4.359308599999999 + 57.48 4.1928366 + 63.98 4.042262099999999 +71.04000000000001 3.9066771 + 78.7 3.7850351 + 87.02 3.6758656 +96.04000000000001 3.5790707 + 105.8 3.4925206 + 116.4 3.4157166 + 127.86 3.348173 + 140.28 3.2893865 + 153.72 3.2381728 + 168.24 3.1929787 + 183.94 3.1504534 + 200.9 3.1074085 + 219.2 3.0652902 + 238.94 3.0227965 + 260.2 2.9778608 + 283.08 2.9364654 + 307.7 2.8894625 + 334.16 2.8333964 + 362.54 2.7733809 + 392.96 2.7078265 + 425.52 2.6379633 + 460.34 2.5613083 + 497.5 2.4824873 + 537.1 2.4085961 + 579.22 2.3303677 + 623.98 2.2573193 +671.4000000000001 2.183273 + 721.58 2.0958828 + 774.54 2.0071174 + 830.34 1.9231549 + 888.96 1.848519 +950.4000000000001 1.7529777 + 1014.64 1.65451 + 1081.62 1.5366693 + 1151.26 1.4209837 + 1223.44 1.290933 + 1298.04 1.1705889 + 1374.88 1.0512143 + 1453.78 0.92785728 + 1534.52 0.82042875 + 1616.88 0.72361343 + 1700.58 0.64354871 + 1785.36 0.58830353 + 1870.902 0.5446869999999999 + 1956.914 0.52843626 + diff --git a/Tests/test_build.sh b/Tests/test_build.sh index cd3ae9bcd..75b2db698 100755 --- a/Tests/test_build.sh +++ b/Tests/test_build.sh @@ -907,7 +907,7 @@ function processor_full_length_tests { "$RANS_CHANNEL_LR_RSM_DIR/Xmgrace" echo "================================= TEST 8 ==============================" - # [~12 HOURS test] + # [~4.5 HOURS test] processor_full_length_test \ "$HYB_CHANNEL_HR_UNIFORM_DIR" \ "hybrid_les_rans" \ From 390bb73965596e73bd3850257cd27f9e9e7de858 Mon Sep 17 00:00:00 2001 From: Egor Palkin Date: Tue, 19 Mar 2019 17:55:27 +0700 Subject: [PATCH 09/32] Updated Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh [minor update] --- .../Uniform_Mesh/Xmgrace/gnuplot_script_template.sh | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh index 578856866..9e7e19448 100644 --- a/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh +++ b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh @@ -91,10 +91,11 @@ plot \ # ----- k_plus unset logscale x -set format x "$%2.1t\\\cdot10^{%L}$" +set format x "$%2.1t\\!\\cdot\\!10^{%L}$" set xlabel '$y^+$' offset 0., 0.5 set xrange [1e-1:2e3] -set xtics offset 0,0.0 border 0,500,1500 scale 5 +set xtics offset 0,0.0 border 0,500,1900 scale 5 +set xtics add ("0" 0) set mxtics 5 set ylabel "$k^+$" offset 0.5, 0.0 From 1b998935c42e216d75956c9cf62365aa9e881149 Mon Sep 17 00:00:00 2001 From: Egor Palkin Date: Wed, 20 Mar 2019 10:29:47 +0700 Subject: [PATCH 10/32] Updated test_build.sh for Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Stretched_Mesh/ --- .../Stretched_Mesh/Xmgrace/k_plus_Re_tau_2000.dat | 1 + Tests/test_build.sh | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) create mode 120000 Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Stretched_Mesh/Xmgrace/k_plus_Re_tau_2000.dat diff --git a/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Stretched_Mesh/Xmgrace/k_plus_Re_tau_2000.dat b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Stretched_Mesh/Xmgrace/k_plus_Re_tau_2000.dat new file mode 120000 index 000000000..5e3960d55 --- /dev/null +++ b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Stretched_Mesh/Xmgrace/k_plus_Re_tau_2000.dat @@ -0,0 +1 @@ +../../Uniform_Mesh/Xmgrace/k_plus_Re_tau_2000.dat \ No newline at end of file diff --git a/Tests/test_build.sh b/Tests/test_build.sh index 75b2db698..c813acf99 100755 --- a/Tests/test_build.sh +++ b/Tests/test_build.sh @@ -10,7 +10,7 @@ FCOMP="gnu" # Conduct tests with DEBUG=yes DEBUG="no" # Repeat tests with CGNS_HDF5=yes -CGNS="yes" +CGNS="no" # A small reminder how to set up alternatives if you have mpich and openmpi: #update-alternatives --install /usr/bin/mpif90 mpif90 /usr/bin/mpif90.openmpi 20 @@ -914,7 +914,7 @@ function processor_full_length_tests { "$HYB_CHANNEL_HR_UNIFORM_DIR/Xmgrace" echo "================================= TEST 9 ==============================" - # [~12 HOURS test] + # [~4 HOURS test] processor_full_length_test \ "$HYB_CHANNEL_HR_STRETCHED_DIR" \ "hybrid_les_rans" \ From da31ea231653257c4413be632a4b92075ca404bc Mon Sep 17 00:00:00 2001 From: Egor Palkin Date: Wed, 20 Mar 2019 11:04:56 +0700 Subject: [PATCH 11/32] Visually improved gnuplot_script_template.sh --- .../Uniform_Mesh/Xmgrace/gnuplot_script_template.sh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh index 9e7e19448..fd36fdc35 100644 --- a/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh +++ b/Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Uniform_Mesh/Xmgrace/gnuplot_script_template.sh @@ -104,10 +104,10 @@ set ytics offset 0,0.0 border 0,1,7 scale 5 set mytics 2 plot \ -1/0 with lp ls 1 ps 5 t '$k_{tot} \, Current$', \ -'DAT_FILE_WITH_RESULTS_MACRO' usi ($1):($5) with lp ls 1 not, \ -1/0 with lp ls 2 ps 5 t '$k_{tot} \, ref. LES$', \ -'k_plus_Re_tau_2000.dat' usi ($1):($2) with lp ls 2 not +1/0 with lp ls 1 ps 5 t '$k_{tot} \: ref. \: LES$', \ +'k_plus_Re_tau_2000.dat' usi ($1):($2) with lp ls 1 not, \ +1/0 with lp ls 2 ps 5 t '$k_{tot} \: Current$', \ +'DAT_FILE_WITH_RESULTS_MACRO' usi ($1):($5) with lp ls 2 not # necessary line unset multiplot \ No newline at end of file From b4f21c4562d798929f21c3e737c16894d9ca23a0 Mon Sep 17 00:00:00 2001 From: palkinev <31729290+palkinev@users.noreply.github.com> Date: Thu, 21 Mar 2019 14:41:52 +0700 Subject: [PATCH 12/32] statistics in RSM models caused a crash Process/Turbulence/Allocate.f90 function did not have allocation block for the first moments statistics (u % mean, v % mean,w % mean, p % mean). --- Sources/Process/Turbulence/Allocate.f90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Sources/Process/Turbulence/Allocate.f90 b/Sources/Process/Turbulence/Allocate.f90 index 169b4683e..dbf100206 100644 --- a/Sources/Process/Turbulence/Allocate.f90 +++ b/Sources/Process/Turbulence/Allocate.f90 @@ -228,6 +228,13 @@ subroutine Turbulence_Allocate(flow) if(turbulence_statistics) then + ! First moments + call Var_Mod_Allocate_Statistics(u) + call Var_Mod_Allocate_Statistics(v) + call Var_Mod_Allocate_Statistics(w) + call Var_Mod_Allocate_Statistics(p) + + ! Second moments call Var_Mod_Allocate_Statistics(uu) call Var_Mod_Allocate_Statistics(vv) call Var_Mod_Allocate_Statistics(ww) From e175c3db20012efbcd55edc3c7f16bb50ecc9256 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Wed, 27 Mar 2019 14:13:06 +0100 Subject: [PATCH 13/32] Initial commit. --- .../User_Mod/Cut_lines_ht_cp.f90 | 335 ++++++++++++++++++ .../Askevein_hill/User_Mod/Save_Results.f90 | 20 ++ .../Askevein_hill/User_Mod/Tecplot_plane.f90 | 106 ++++++ 3 files changed, 461 insertions(+) create mode 100755 Tests/Rans/Askevein_hill/User_Mod/Cut_lines_ht_cp.f90 create mode 100755 Tests/Rans/Askevein_hill/User_Mod/Save_Results.f90 create mode 100755 Tests/Rans/Askevein_hill/User_Mod/Tecplot_plane.f90 diff --git a/Tests/Rans/Askevein_hill/User_Mod/Cut_lines_ht_cp.f90 b/Tests/Rans/Askevein_hill/User_Mod/Cut_lines_ht_cp.f90 new file mode 100755 index 000000000..d90685b8c --- /dev/null +++ b/Tests/Rans/Askevein_hill/User_Mod/Cut_lines_ht_cp.f90 @@ -0,0 +1,335 @@ +!==============================================================================! + subroutine User_Mod_Cut_lines_ht_cp(flow, save_name) +!------------------------------------------------------------------------------! +! This is a prototype of a function for customized source for scalar. ! +! It is called from "Compute_Scalar" function, just before calling the ! +! linear solver. Both system matrix ("a_matrix") and right hand side ! +! vector ("b_vector") are sent should the user want to stabilize the ! +! system for always positive variables, for example. ! +!------------------------------------------------------------------------------! +!----------------------------------[Modules]-----------------------------------! + use Field_Mod, only: Field_Type, heat_transfer, heat_flux, & + density, viscosity, capacity, conductivity + use Grid_Mod, only: Grid_Type + use Var_Mod, only: Var_Type + use Matrix_Mod, only: Matrix_Type + use Bulk_Mod, only: Bulk_Type + use Rans_Mod + use Comm_Mod ! parallel stuff + use Name_Mod, only: problem_name + use Const_Mod ! constants +!------------------------------------------------------------------------------! + implicit none +!---------------------------------[Arguments]----------------------------------! + type(Field_Type), target :: flow + character(len=*) :: save_name +! type(Var_Type), target :: phi +! type(Matrix_Type) :: a_matrix +! real, dimension(:) :: b_vector +!----------------------------------[Locals]------------------------------------! + type(Grid_Type), pointer :: grid + type(Bulk_Type), pointer :: bulk + type(Var_Type), pointer :: u, v, w, t + real, pointer :: flux(:) + !------------------------------------------------------------------------------! + INTERFACE + LOGICAL FUNCTION Approx_Real(A,B,tol) + REAL :: A,B + REAL, OPTIONAL :: tol + END FUNCTION Approx_Real + END INTERFACE + +!-----------------------------[Parameters]-----------------------------! + + REAL :: Rad_2, Ufric, Rad_1 +!------------------------------[Calling]-------------------------------! + INTEGER :: Nprob, pl, c, i, count, k, c1, c2, s, N_hor, n + INTEGER :: Nfound + CHARACTER :: name_ht*12, name_cp*12, name_inflow*12 + character(len=80) :: res_name_cp, res_name_ht + character(len=80) :: store_name + REAL,ALLOCATABLE :: y_p(:), x_p(:), z_p(:), Ump(:), Vmp(:), Wmp(:), & + uup(:), vvp(:), wwp(:), & + uvp(:), uwp(:), vwp(:), & + Tmp(:), TTp(:), & + uTp(:), vTp(:), wTp(:), & + Ksgsp(:), & + var_1(:), var_2(:), var_3(:), Rad_mp(:), & + var_4(:), var_5(:), R_p(:) + INTEGER,ALLOCATABLE :: Np(:), Ncount(:) + REAL :: R, Urad_n, Utan_n, R1, R2, Urad, Utan, dummy + REAL :: xc_max, yc_max, zc_max, range + logical :: there +!==============================================================================! + ! Take aliases + grid => flow % pnt_grid + flux => flow % flux + bulk => flow % bulk + u => flow % u + v => flow % v + w => flow % w + t => flow % t + +!--------------------------------------------------------------------! + + ! Store the name + store_name = problem_name + problem_name = save_name + + call Name_File(0, res_name_ht, "-ht.dat") + call Name_File(0, res_name_cp, "-cp.dat") + + N_hor = 2 + Nfound = 0 + zc_max = 0.0 + range = HUGE + + allocate(x_p(N_hor)) + allocate(y_p(N_hor)) + allocate(R_p(N_hor)) + + do c = -1, -grid % n_bnd_cells, -1 + if( Grid_Mod_Bnd_Cond_Type(grid, c) .eq. WALL .or. & + Grid_Mod_Bnd_Cond_Type(grid, c) .eq. WALLFL) then + if(grid % zc(c) > zc_max) then + zc_max = grid % zc(c) + end if + end if + end do + + call Comm_Mod_Global_Max_Real(zc_max) + + if(this_proc< 2) write(*,*) 'Hill top height is ', zc_max + + call Comm_Mod_Wait + + xc_max = 0.0 + yc_max = 0.0 + + do c = 1, grid % n_cells + range = min(range, grid % delta(c)) + end do + + call Comm_Mod_Global_Min_Real(range) + + do c = -1, -grid % n_bnd_cells, -1 + if( Grid_Mod_Bnd_Cond_Type(grid, c) .eq. WALL .or. & + Grid_Mod_Bnd_Cond_Type(grid, c) .eq. WALLFL) then + if(Approx_Real( grid % zc(c), zc_max)) then + xc_max = grid % xc(c) + yc_max = grid % yc(c) + end if + end if + end do + + call Comm_Mod_Global_Sum_Real(xc_max) + call Comm_Mod_Global_Sum_Real(yc_max) + + if(this_proc< 2) write(*,*)'Coordinate of hill top (HT)', xc_max, yc_max , zc_max + + call Comm_Mod_Wait + + ! HT + x_p(1) = xc_max + y_p(1) = yc_max + + !CP ! + !distance between HT and CP is 400 m ! + +! x_p(2) = xc_max - 410.0*cos(78.0*3.14159/180.0) +! y_p(2) = yc_max - 410.0*sin(73.0*3.14159/180.0) + + x_p(2) = xc_max - .5*cos(78.0*3.14159/180.0) + y_p(2) = yc_max - .5*sin(73.0*3.14159/180.0) +!-Radius for averaging --adjust radius according to number of counts + R_p(1) = range + R_p(2) = range + + !------------------! + ! Read 1d file ! + !------------------! + inquire(file='Z_coordinate.dat', exist=there) + if(.not. there) then + if(this_proc < 2) then + print *, '===============================================================' + print *, 'In order to extract profiles and write them in ascii files' + print *, 'the code has to read cell-faces coordinates ' + print *, 'in wall-normal direction in the ascii file ''case_name.1d.''' + print *, 'The file format should be as follows:' + print *, '10 ! number of cells + 1' + print *, '1 0.0' + print *, '2 0.1' + print *, '3 0.2' + print *, '... ' + print *, '===============================================================' + end if + + ! Restore the name and return + problem_name = store_name + return + end if + + if(this_proc< 2) write(6, *) '# Now reading Z_coordinate.dat ' + + open(9, FILE='Z_coordinate.dat') +!---- write the number of probes + read(9,*) Nprob + allocate(z_p(Nprob)) + +!---- write the probe coordinates out + do pl=1,Nprob + read(9,*) dummy, z_p(pl) + end do + close(9) + + allocate(Np(Nprob)); Np = 0 + allocate(Ump(Nprob)); Ump = 0.0 + allocate(uwp(Nprob)); uwp = 0.0 + allocate(Rad_mp(Nprob)); Rad_mp = 0.0 + allocate(var_1(Nprob)); var_1 = 0.0 + allocate(var_2(Nprob)); var_2 = 0.0 + allocate(var_3(Nprob)); var_3 = 0.0 + allocate(var_4(Nprob)); var_4 = 0.0 + + allocate(Ncount(Nprob)); Ncount=0 + count = 0 + +!+++++++++++++++++++++++++++++! +! average the results ! +!+++++++++++++++++++++++++++++! + + do k = 1, 2 + do i = 1, Nprob-1 + 101 continue + Ncount(i) = 0 + do c = 1, grid % n_cells - grid % comm % n_buff_cells + Rad_1 = sqrt((grid % xc(c)-x_p(k))**2 + (grid % yc(c)-y_p(k))**2) + Rad_2 = grid % wall_dist(c) !grid % zc(c) + if(Rad_1 < R_p(k)) then + if(Rad_2 > z_p(i) .and. Rad_2 < z_p(i+1)) then + Rad_mp(i) = Rad_mp(i) + Rad_2 + Ncount(i) = Ncount(i) + 1 + end if + end if + end do + + call Comm_Mod_Wait + call Comm_Mod_Global_Sum_Int(Ncount(i)) + + if(Ncount(i) /= 0) then + Nfound = Ncount(i) !max(Nfound, Ncount(i)) + end if + + call Comm_Mod_Wait + + if(Nfound > 6) then + R_p(k) = R_p(k) * 0.9 + write(*,*) 'vise od 6', Nfound, R_p(k), k + goto 101 + else + continue + end if +! write(*,*) 'It found ', Nfound, ' points!' + + end do + end do + + do i = 1, Nprob + Ncount(i) = 0 + Rad_mp(i) = 0.0 + end do + +! write(*,*) x_p(2), x_p(1) +! write(*,*) y_p(2), y_p(1) + + do k = 1, 2 + do i = 1, Nprob-1 + do c = 1, grid % n_cells + Rad_1 = sqrt((grid % xc(c)-x_p(k))**2 + (grid % yc(c)-y_p(k))**2) + Rad_2 = grid % wall_dist(c) + if(Rad_1 < R_p(k)) then + if(Rad_2 > z_p(i) .and. Rad_2 < z_p(i+1)) then + Ncount(i) = Ncount(i) + 1 + Ump(i) = Ump(i) + u % n(c) + uwp(i) = uwp(i) + vis_t(c)*(u % z(c) + w % x(c)) + var_1(i) = var_1(i) + kin % n(c) + var_2(i) = var_2(i) + eps % n(c) + var_3(i) = var_3(i) + 5.64 + 1.41*log(grid % wall_dist(c)) + var_4(i) = var_4(i) + grid % wall_dist(c) + Rad_mp(i) = Rad_mp(i) + Rad_2 + end if + end if + end do + end do + +!---- average over all processors + do pl=1, Nprob + call Comm_Mod_Global_Sum_Int(Ncount(pl)) + call Comm_Mod_Global_Sum_Real(Ump(pl)) + call Comm_Mod_Global_Sum_Real(Rad_mp(pl)) + call Comm_Mod_Global_Sum_Real(uwp(pl)) + call Comm_Mod_Global_Sum_Real(var_1(pl)) + call Comm_Mod_Global_Sum_Real(var_2(pl)) + call Comm_Mod_Global_Sum_Real(var_3(pl)) + call Comm_Mod_Global_Sum_Real(var_4(pl)) + + count = count + Ncount(pl) + end do + + if(k == 1) then + open(k,FILE=res_name_ht) + else if(k == 2) then + open(k,FILE=res_name_cp) + end if + + do i = 1, Nprob + if(Ncount(i) /= 0) then + Ump(i) = Ump(i)/Ncount(i) + uwp(i) = uwp(i)/Ncount(i) + var_1(i) = var_1(i)/Ncount(i) + var_2(i) = var_2(i)/Ncount(i) + var_3(i) = var_3(i)/Ncount(i) + var_4(i) = var_4(i)/Ncount(i) + Rad_mp(i) = Rad_mp(i)/Ncount(i) + +!--------------Write variables ------------------------------------------------------------! + + write(k,'(3E15.7, I7)') var_4(i), (Ump(i) - var_3(i))/var_3(i), & + (var_1(i)/(8.9**2)), Ncount(i) + +!------------------------------------------------------------------------------------------! + + Ump(i) = 0.0 + uwp(i) = 0.0 + var_1(i) = 0.0 + var_2(i) = 0.0 + var_3(i) = 0.0 + var_4(i) = 0.0 + Rad_mp(i) = 0.0 + + Ncount(i) = 0 + end if + end do + + close(k) + end do + + deallocate(Np) + deallocate(z_p) + deallocate(Ump) + deallocate(uwp) + deallocate(var_1) + deallocate(var_2) + deallocate(var_3) + deallocate(var_4) + deallocate(Rad_mp) + deallocate(Ncount) + + call Comm_Mod_Wait + + if(this_proc< 2) write(*,*) 'Finished with user function for Askervein hill ! ' + + ! Restore the name + problem_name = store_name + + END SUBROUTINE User_Mod_Cut_lines_ht_cp diff --git a/Tests/Rans/Askevein_hill/User_Mod/Save_Results.f90 b/Tests/Rans/Askevein_hill/User_Mod/Save_Results.f90 new file mode 100755 index 000000000..3516d71c0 --- /dev/null +++ b/Tests/Rans/Askevein_hill/User_Mod/Save_Results.f90 @@ -0,0 +1,20 @@ +include '../User_Mod/Tecplot_plane.f90' +include '../User_Mod/Cut_lines_ht_cp.f90' + +!==============================================================================! + subroutine User_Mod_Save_Results(flow, save_name) +!------------------------------------------------------------------------------! +! Calls User_Impinging_Jet_Nu and User_Impinging_Jet_Profile functions. ! +!------------------------------------------------------------------------------! + use Field_Mod +!------------------------------------------------------------------------------! + implicit none +!---------------------------------[Arguments]----------------------------------! + type(Field_Type) :: flow + character(len=*) :: save_name +!==============================================================================! + + call User_Mod_Tecplot_plane (flow, save_name) + call User_Mod_Cut_lines_ht_cp(flow, save_name) + + end subroutine diff --git a/Tests/Rans/Askevein_hill/User_Mod/Tecplot_plane.f90 b/Tests/Rans/Askevein_hill/User_Mod/Tecplot_plane.f90 new file mode 100755 index 000000000..ca8a032af --- /dev/null +++ b/Tests/Rans/Askevein_hill/User_Mod/Tecplot_plane.f90 @@ -0,0 +1,106 @@ +!==============================================================================! + subroutine User_Mod_Tecplot_plane(flow, save_name) +!------------------------------------------------------------------------------! +! This is a prototype of a function for customized source for scalar. ! +! It is called from "Compute_Scalar" function, just before calling the ! +! linear solver. Both system matrix ("a_matrix") and right hand side ! +! vector ("b_vector") are sent should the user want to stabilize the ! +! system for always positive variables, for example. ! +!------------------------------------------------------------------------------! +!----------------------------------[Modules]-----------------------------------! + use Field_Mod, only: Field_Type, heat_flux + use Grid_Mod, only: Grid_Type + use Var_Mod, only: Var_Type + use Matrix_Mod, only: Matrix_Type + use Bulk_Mod, only: Bulk_Type + use Rans_Mod + use Comm_Mod ! parallel stuff + use Name_Mod, only: problem_name + use Const_Mod ! constants + +!------------------------------------------------------------------------------! + implicit none +!---------------------------------[Arguments]----------------------------------! + type(Field_Type), target :: flow + character(len=*) :: save_name +!-----------------------------------[Locals]-----------------------------------! + type(Grid_Type), pointer :: grid + type(Bulk_Type), pointer :: bulk + type(Var_Type), pointer :: u, v, w, t + real, pointer :: flux(:) + + INTEGER :: n + REAL :: x1_dis, x0_dis, y1_dis, y0_dis, z1_dis, z0_dis + REAL :: TauWup, TauWdown + REAL :: Utan, UnorSq, Unor, UtotSq, dely, Stot, R, UtanX +!-----------------------------------[Locals]-----------------------------------! + INTEGER :: c, s, c1, c2 + CHARACTER*80 :: namout, res_name, name_out_9 + CHARACTER*39 :: path + character(len=80) :: store_name +!------------------------------------------------------------------------------! + ! Take aliases + grid => flow % pnt_grid + flux => flow % flux + bulk => flow % bulk + u => flow % u + v => flow % v + w => flow % w + t => flow % t +!------------------------------------------------------------------------------! + + ! Store the name + store_name = problem_name + problem_name = save_name + + call Name_File(this_proc, res_name, "-tec360.dat") +! call Name_File(this_proc, name_out_9, '-tec360.dat') +! open(9, FILE='slice.dat') + +! if(this_proc < 2) write(*,*) '# Now reading the file: slice.dat ' +! read(9,*) z1_dis, z0_dis +! if(this_proc < 2) write(*,*) '# X:[ ', x0_dis, " ;", x1_dis, "]", & +! ' Y:[ ', y0_dis, " ;", y1_dis, "]", ' Z:[ ', z0_dis, " ;", z1_dis, "]" + +! inflowfile = 'tecplt-10m-xxxxxx-xx.dat' + +! write(inflowfile(12:17),'(I6.6)') n +! write(inflowfile(19:20),'(I2.2)') this_proc + open(500+this_proc, FILE=res_name) + + if( this_proc < 2 ) then + write(*,*) 'Capturing field.. ', res_name + end if + + + do s = 1, grid % n_faces + c1 = grid % faces_c(1,s) + c2 = grid % faces_c(2,s) + + if(c2 < 0) then + if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. WALL .or. & + Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. WALLFL) then + + write (500+this_proc,'(8E17.7E3)') grid % xc(c1), grid % yc(c1),& + grid % zc(c1), U % n(c1), V % n(c1), W % n(c1), Kin % n(c1), & + sqrt( (U % n(c1))**2+ (V % n(c1))**2+ (W % n(c1))**2) + + endif + endif + end do + + call Comm_Mod_Wait + + close(9) + close(500+this_proc) + + call Comm_Mod_Wait + + if ( this_proc < 2 ) then + write(*,*) 'It is done' + end if + + ! Restore the name + problem_name = store_name + + END SUBROUTINE From 1800f8281c5b8f6e1d0bb0913618509a9d46329b Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sat, 13 Apr 2019 16:40:59 +0200 Subject: [PATCH 14/32] Limiter for hydraulic roughness z_o is added in order to prevent unphysical values of u_plus. --- Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 b/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 index 884098df0..c77158a3b 100644 --- a/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 +++ b/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 @@ -126,12 +126,12 @@ subroutine Calculate_Vis_T_K_Eps(flow) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) + z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) y_plus(c1) = Y_Plus_Rough_Walls(u_tau(c1), & grid % wall_dist(c1), & kin_vis) u_plus = U_Plus_Rough_Walls(grid % wall_dist(c1)) - vis_wall(c1) = y_plus(c1) * viscosity * kappa & - / log((grid % wall_dist(c1)+z_o)/z_o) ! is this U+? + vis_wall(c1) = y_plus(c1) * viscosity / u_plus end if if(heat_transfer) then From 7db6c1926f29cf58279d3197512ce8c42384c4d4 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sat, 13 Apr 2019 16:42:05 +0200 Subject: [PATCH 15/32] Limiter for hydraulic roughness z_o is added in order to prevent unphysical values of u_plus. --- Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 b/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 index 782716188..b02d3904f 100644 --- a/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 +++ b/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 @@ -119,12 +119,12 @@ subroutine Calculate_Vis_T_K_Eps_Zeta_F(flow) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) + z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) y_plus(c1) = Y_Plus_Rough_Walls(u_tau(c1), & grid % wall_dist(c1), & kin_vis) u_plus = U_Plus_Rough_Walls(grid % wall_dist(c1)) - vis_wall(c1) = y_plus(c1) * viscosity * kappa & - / log((grid % wall_dist(c1)+z_o)/z_o) ! is this U+? + vis_wall(c1) = y_plus(c1) * viscosity / u_plus end if if(heat_transfer) then From 8485a8ee2dddf13bb6b9cddf61e0c6d8861c433b Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sat, 13 Apr 2019 16:42:20 +0200 Subject: [PATCH 16/32] Limiter for hydraulic roughness z_o is added in order to prevent unphysical values of u_plus. --- Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 b/Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 index 9b5772362..9d2fdaebb 100644 --- a/Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 +++ b/Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 @@ -132,7 +132,8 @@ subroutine Source_Kin_K_Eps_Zeta_F(flow, sol) end if if(rough_walls) then - z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) + z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) + z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) u_tau(c1) = c_mu25 * sqrt(kin % n(c1)) y_plus(c1) = Y_Plus_Rough_Walls(u_tau(c1), & grid % wall_dist(c1), kin_vis) @@ -153,7 +154,7 @@ subroutine Source_Kin_K_Eps_Zeta_F(flow, sol) ebf = max(0.01 * y_plus(c1)**4 / (1.0 + 5.0*y_plus(c1)),tiny) - p_kin_wf = tau_wall(c1) * 0.07**0.25 * sqrt(kin % n(c1)) & + p_kin_wf = tau_wall(c1) * c_mu25 * sqrt(kin % n(c1)) & / (grid % wall_dist(c1) * kappa) p_kin_int = vis_t(c1) * shear(c1)**2 From 8b75ae3fc294b19ad2e85ab11b6ec65317b7c866 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sat, 13 Apr 2019 16:42:48 +0200 Subject: [PATCH 17/32] Limiter for hydraulic roughness z_o is added in order to prevent unphysical values of u_plus. Bug in rough_wall section is corrected. --- Sources/Process/Turbulence/Source_Kin_K_Eps.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 b/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 index ea3c507c5..cfbf737bf 100644 --- a/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 +++ b/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 @@ -107,6 +107,7 @@ subroutine Source_Kin_K_Eps(flow, sol) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) + z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) u_tau(c1) = c_mu25 * sqrt(kin % n(c1)) y_plus(c1) = u_tau(c1) * (grid % wall_dist(c1) + z_o) & / kin_vis @@ -116,7 +117,7 @@ subroutine Source_Kin_K_Eps(flow, sol) p_kin(c1) = density * tau_wall(c1) * c_mu25 * sqrt(kin % n(c1)) & / (kappa*(grid % wall_dist(c1) + z_o)) - b(c1) = b(c1) + (p_kin(c1) - p_kin(c1)) & + b(c1) = b(c1) + (p_kin(c1) - vis_t(c1) * shear(c1)**2) & * grid % vol(c1) else u_tau(c1) = c_mu25 * sqrt(kin % n(c1)) @@ -127,7 +128,7 @@ subroutine Source_Kin_K_Eps(flow, sol) ebf = 0.01 * y_plus(c1)**4 / (1.0 + 5.0*y_plus(c1)) - p_kin_wf = tau_wall(c1) * 0.07**0.25 * sqrt(kin % n(c1)) & + p_kin_wf = tau_wall(c1) * c_mu25 * sqrt(kin % n(c1)) & / (grid % wall_dist(c1) * kappa) p_kin_int = vis_t(c1) * shear(c1)**2 From 459bcd362ac4acabde19f68251f8775ebd09bc0a Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sat, 13 Apr 2019 16:43:44 +0200 Subject: [PATCH 18/32] Limiter for hydraulic roughness z_o is added in order to prevent unphysical values of u_plus. EBF used in old definition of blending function is deleted. --- Sources/Process/Turbulence/Source_Eps_K_Eps.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Sources/Process/Turbulence/Source_Eps_K_Eps.f90 b/Sources/Process/Turbulence/Source_Eps_K_Eps.f90 index d0fb4a320..2db535060 100644 --- a/Sources/Process/Turbulence/Source_Eps_K_Eps.f90 +++ b/Sources/Process/Turbulence/Source_Eps_K_Eps.f90 @@ -35,7 +35,7 @@ subroutine Source_Eps_K_Eps(flow, sol) integer :: s, c, c1, c2, j real :: u_tan, u_nor_sq, u_nor, u_tot_sq real :: re_t, f_mu, u_tau_new, fa, kin_vis - real :: eps_wf, eps_int, ebf, y_star + real :: eps_wf, eps_int, y_star !==============================================================================! ! Dimensions: ! ! ! @@ -116,7 +116,8 @@ subroutine Source_Eps_K_Eps(flow, sol) end if if(rough_walls) then - z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) + z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) + z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) eps % n(c1) = c_mu75 * kin % n(c1)**1.5 & / ((grid % wall_dist(c1) + z_o) * kappa) @@ -136,7 +137,6 @@ subroutine Source_Eps_K_Eps(flow, sol) u_tau_new = sqrt(tau_wall(c1)/density) y_plus(c1) = Y_Plus_Low_Re(u_tau_new, grid % wall_dist(c1), kin_vis) - ebf = 0.01 * y_plus(c1)**4 / (1.0 + 5.0*y_plus(c1)) eps_int = 2.0*viscosity/density * kin % n(c1) & / grid % wall_dist(c1)**2 From 535157c64b2276237a6006f520568f074bb6aee2 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sat, 13 Apr 2019 16:44:36 +0200 Subject: [PATCH 19/32] Limiter for hydraulic roughness z_o is added in order to prevent unphysical values of u_plus. EBF used in old definition of blending function is deleted. --- Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 b/Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 index 1591ec1e4..bca16b2e4 100644 --- a/Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 +++ b/Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 @@ -27,7 +27,7 @@ subroutine Source_Eps_K_Eps_Zeta_F(flow, sol) real, pointer :: b(:) integer :: c, s, c1, c2, j real :: u_tan, u_nor_sq, u_nor, u_tot_sq - real :: e_sor, c_11e, ebf + real :: e_sor, c_11e real :: eps_wf, eps_int real :: fa, u_tau_new, kin_vis !==============================================================================! @@ -115,6 +115,7 @@ subroutine Source_Eps_K_Eps_Zeta_F(flow, sol) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) + z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) eps % n(c1) = c_mu75 * kin % n(c1)**1.5 / & ((grid % wall_dist(c1) + z_o) * kappa) @@ -135,8 +136,6 @@ subroutine Source_Eps_K_Eps_Zeta_F(flow, sol) u_tau_new = sqrt(tau_wall(c1)/density) y_plus(c1) = Y_Plus_Low_Re(u_tau_new, grid % wall_dist(c1), kin_vis) - ebf = 0.001 * y_plus(c1)**4 / (1.0 + y_plus(c1)) - eps_int = 2.0* kin_vis * kin % n(c1) & / grid % wall_dist(c1)**2 eps_wf = c_mu75 * kin % n(c1)**1.5 & @@ -148,7 +147,6 @@ subroutine Source_Eps_K_Eps_Zeta_F(flow, sol) / (kappa*grid % wall_dist(c1) * p_kin(c1)), & 1.0) eps % n(c1) = (1.0 - fa) * eps_int + fa * eps_wf - ! Adjusting coefficient to fix eps value in near wall calls do j = a % row(c1), a % row(c1 + 1) - 1 a % val(j) = 0.0 From a018305d877295ca2dad08e509b0e48a5b79b5ca Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sat, 13 Apr 2019 16:45:06 +0200 Subject: [PATCH 20/32] Coefficient 0.6 is changed to 0.8 as 0.6 is too low and is changing time scale far from the wall which is not supposed to happened. 0.6 value was influencing rough wall results as well as channel flow with very high Reynolds number. The motivation for adopting 0.6 was in improving Nusselt number in impinging jets. However, theoretical number of this coefficient is around 1. --- Sources/Process/Turbulence/Time_And_Length_Scale.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sources/Process/Turbulence/Time_And_Length_Scale.f90 b/Sources/Process/Turbulence/Time_And_Length_Scale.f90 index 1e546174b..bb3a634cc 100644 --- a/Sources/Process/Turbulence/Time_And_Length_Scale.f90 +++ b/Sources/Process/Turbulence/Time_And_Length_Scale.f90 @@ -45,7 +45,7 @@ subroutine Time_And_Length_Scale(grid) t_1(c) = kin % n(c)/eps_l(c) t_2(c) = c_t*sqrt(kin_vis/eps_l(c)) - t_3(c) = 0.6/(sqrt(3.0)*c_mu_d * zeta % n(c) * shear(c) + TINY) + t_3(c) = 0.8/(sqrt(3.0)*c_mu_d * zeta % n(c) * shear(c) + TINY) l_1(c) = kin % n(c)**1.5/eps_l(c) l_2(c) = c_nu * (kin_vis**3 / eps_l(c))**0.25 From 1eee314d55c148835555d92b5d3eac067a889e31 Mon Sep 17 00:00:00 2001 From: palkinev <31729290+palkinev@users.noreply.github.com> Date: Mon, 15 Apr 2019 15:22:34 +0700 Subject: [PATCH 21/32] Updated test_build.sh: process_save_exit_now_test In the issue #133 Bojan fairly noticed that save_exit_now test case behaved somehow incorrectly. In particular, string "save_exit_now_test was successful" has never appeared. It turned out that test relied on the output in stdout rather than in $LOG_FILE It is fixed now. Moreover, scenario which this test case reproduce is now more intricate and accurate than before. --- Tests/test_build.sh | 138 ++++++++++++++++++++++++++++++++++++++------ 1 file changed, 121 insertions(+), 17 deletions(-) diff --git a/Tests/test_build.sh b/Tests/test_build.sh index c813acf99..178a0dc46 100755 --- a/Tests/test_build.sh +++ b/Tests/test_build.sh @@ -545,55 +545,159 @@ function process_save_exit_now_test { name_in_div=$(head -n1 divide.scr) nproc_in_div=$(head -n2 divide.scr | tail -n1) + if [ -f "save_now" ]; then rm save_now; fi + if [ -f "exit_now" ]; then rm exit_now; fi + # change number of timesteps to 3 replace_line_with_first_occurence_in_file "NUMBER_OF_TIME_STEPS" \ "NUMBER_OF_TIME_STEPS 3" control - # change backup interval to 1 ts + # change backup interval to 10 ts replace_line_with_first_occurence_in_file "BACKUP_SAVE_INTERVAL" \ - "BACKUP_SAVE_INTERVAL 1" control + "BACKUP_SAVE_INTERVAL 10" control + + # comment line with LOAD_BACKUP_NAME + n1=$(printf "%06d" 1) + replace_line_with_first_occurence_in_file "LOAD_BACKUP_NAME" \ + "#LOAD_BACKUP_NAME "$name_in_div"-ts"$n1".backup" control echo "================================= TEST 1 ==============================" echo "np=1, MPI=no" clean_compile $PROC_DIR $1 no # dir CGNS_HDF5 MPI cd $2 - echo "save_now" + echo "Forcing to save: save_now" touch save_now - if launch_process seq 1 | grep -q ""$name_in_div"-ts000001"; then - echo "exit_now" + # get current line count where search starts + n_start="$(echo "$(wc -l $FULL_LOG | cut -d" " -f1) + 1" | bc -l)" + + # start from scratch + launch_process seq 1 + + # find if save was made in the range [n_start: n_finish] + if tail -n+$n_start $FULL_LOG | \ + grep -q "# Creating file: "$name_in_div"-ts"$n1""; then + + echo "save_now was successfull" + + echo "Forcing to exit: exit_now" touch exit_now - launch_process seq 1 | grep -q "# Exiting !" - echo "save_exit_now_test was successfull" + + # uncomment line with LOAD_BACKUP_NAME + replace_line_with_first_occurence_in_file "LOAD_BACKUP_NAME" \ + "LOAD_BACKUP_NAME "$name_in_div"-ts"$n1".backup" control + + # start from ts=1 + n_start="$(echo "$(wc -l $FULL_LOG | cut -d" " -f1) + 1" | bc -l)" + + launch_process seq 1 + + if tail -n+$n_start $FULL_LOG | \ + tr -s " " | \ + grep -q "Time step : 3"; then + + echo "exit_now was NOT successfull" + else + echo "exit_now was successfull" + fi + else + echo "save_now was NOT successfull" fi echo "================================= TEST 2 ==============================" echo "np=1, MPI=yes" + + # comment line with LOAD_BACKUP_NAME + replace_line_with_first_occurence_in_file "LOAD_BACKUP_NAME" \ + "#LOAD_BACKUP_NAME "$name_in_div"-ts"$n1".backup" control + clean_compile $PROC_DIR $1 yes # dir CGNS_HDF5 MPI cd $2 - echo "save_now" + echo "Forcing to save: save_now" touch save_now - if launch_process par 1 | grep -q ""$name_in_div"-ts000001"; then - echo "exit_now" + # get current line count where search starts + n_start="$(echo "$(wc -l $FULL_LOG | cut -d" " -f1) + 1" | bc -l)" + + # start from scratch + launch_process par 1 + + # find if save was made in the range [n_start: n_finish] + if tail -n+$n_start $FULL_LOG | \ + grep -q "# Creating file: "$name_in_div"-ts"$n1""; then + + echo "save_now was successfull" + + echo "Forcing to exit: exit_now" touch exit_now - launch_process par 1 | grep -q "# Exiting !" - echo "save_exit_now_test was successfull" + + # uncomment line with LOAD_BACKUP_NAME + replace_line_with_first_occurence_in_file "LOAD_BACKUP_NAME" \ + "LOAD_BACKUP_NAME "$name_in_div"-ts"$n1".backup" control + + # start from ts=1 + n_start="$(echo "$(wc -l $FULL_LOG | cut -d" " -f1) + 1" | bc -l)" + + launch_process par 1 + + if tail -n+$n_start $FULL_LOG | \ + tr -s " " | \ + grep -q "Time step : 3"; then + + echo "exit_now was NOT successfull" + else + echo "exit_now was successfull" + fi + else + echo "save_now was NOT successfull" fi echo "================================= TEST 3 ==============================" echo "np=2, MPI=yes" - echo "save_now" + # comment line with LOAD_BACKUP_NAME + replace_line_with_first_occurence_in_file "LOAD_BACKUP_NAME" \ + "#LOAD_BACKUP_NAME "$name_in_div"-ts"$n1".backup" control + + echo "Forcing to save: save_now" touch save_now - if launch_process par $nproc_in_div | grep -q ""$name_in_div"-ts000001"; then - echo "exit_now" + # get current line count where search starts + n_start="$(echo "$(wc -l $FULL_LOG | cut -d" " -f1) + 1" | bc -l)" + + # start from scratch + launch_process par $nproc_in_div + + # find if save was made in the range [n_start: n_finish] + if tail -n+$n_start $FULL_LOG | \ + grep -q "# Creating file: "$name_in_div"-ts"$n1""; then + + echo "save_now was successfull" + + echo "Forcing to exit: exit_now" touch exit_now - launch_process par $nproc_in_div | grep -q "# Exiting !" - echo "save_exit_now_test was successfull" + + # uncomment line with LOAD_BACKUP_NAME + replace_line_with_first_occurence_in_file "LOAD_BACKUP_NAME" \ + "LOAD_BACKUP_NAME "$name_in_div"-ts"$n1".backup" control + + # start from ts=1 + n_start="$(echo "$(wc -l $FULL_LOG | cut -d" " -f1) + 1" | bc -l)" + + launch_process par $nproc_in_div + + if tail -n+$n_start $FULL_LOG | \ + tr -s " " | \ + grep -q "Time step : 3"; then + + echo "exit_now was NOT successfull" + else + echo "exit_now was successfull" + fi + else + echo "save_now was NOT successfull" fi } #------------------------------------------------------------------------------# From d3b0d2e97211fe8041cf7fca2babaad8057484e2 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Wed, 17 Apr 2019 15:45:31 +0200 Subject: [PATCH 22/32] Inconsistency in dimension is corrected and max limiter is removed in rough wall part as newly introduced limiter in z_o is sufficient. --- Sources/Process/Turbulence/Source_Kin_K_Eps.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 b/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 index cfbf737bf..a40542a6b 100644 --- a/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 +++ b/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 @@ -113,9 +113,9 @@ subroutine Source_Kin_K_Eps(flow, sol) / kin_vis tau_wall(c1) = density*kappa*u_tau(c1)*u_tan & - / log(max(((grid % wall_dist(c1)+z_o)/z_o), 1.05)) + / log(((grid % wall_dist(c1)+z_o) / z_o)) - p_kin(c1) = density * tau_wall(c1) * c_mu25 * sqrt(kin % n(c1)) & + p_kin(c1) = tau_wall(c1) * c_mu25 * sqrt(kin % n(c1)) & / (kappa*(grid % wall_dist(c1) + z_o)) b(c1) = b(c1) + (p_kin(c1) - vis_t(c1) * shear(c1)**2) & * grid % vol(c1) From ec6d0e7c1dd10897379d01ab285ebd6eb9486147 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Wed, 17 Apr 2019 15:54:34 +0200 Subject: [PATCH 23/32] z_o_f is iniciated as a negative number in order to allowed z_o to be zero in Roughness_Coefficient.f90. --- Sources/Process/Turbulence/Allocate.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sources/Process/Turbulence/Allocate.f90 b/Sources/Process/Turbulence/Allocate.f90 index dbf100206..991a5e582 100644 --- a/Sources/Process/Turbulence/Allocate.f90 +++ b/Sources/Process/Turbulence/Allocate.f90 @@ -123,7 +123,7 @@ subroutine Turbulence_Allocate(flow) ! Hydraulic roughness given by formula if(rough_walls) then - allocate(z_o_f(-grid % n_bnd_cells:grid % n_cells)); z_o_f = 0. + allocate(z_o_f(-grid % n_bnd_cells:grid % n_cells)); z_o_f = -1. end if if(heat_transfer) then From efadd25c7d3337909ce53ccd45ff530f0187c26d Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Wed, 17 Apr 2019 15:56:10 +0200 Subject: [PATCH 24/32] It is now possible to have z_o equals zero and recover smooth wall. --- Sources/Process/Turbulence/Roughness_Coefficient.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sources/Process/Turbulence/Roughness_Coefficient.f90 b/Sources/Process/Turbulence/Roughness_Coefficient.f90 index 154db4126..b02e26009 100644 --- a/Sources/Process/Turbulence/Roughness_Coefficient.f90 +++ b/Sources/Process/Turbulence/Roughness_Coefficient.f90 @@ -17,7 +17,7 @@ real function Roughness_Coefficient(grid, z_o_function, c) Roughness_Coefficient = z_o - if(z_o_function > tiny) then + if(z_o_function > -tiny) then Roughness_Coefficient = z_o_function end if From ebbf93c058d00d42905e5f0cbc920594f5ca3fe3 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Sun, 21 Apr 2019 20:46:16 +0200 Subject: [PATCH 25/32] Askervein hill case is removed. --- .../User_Mod/Cut_lines_ht_cp.f90 | 335 ------------------ .../Askevein_hill/User_Mod/Save_Results.f90 | 20 -- .../Askevein_hill/User_Mod/Tecplot_plane.f90 | 106 ------ 3 files changed, 461 deletions(-) delete mode 100755 Tests/Rans/Askevein_hill/User_Mod/Cut_lines_ht_cp.f90 delete mode 100755 Tests/Rans/Askevein_hill/User_Mod/Save_Results.f90 delete mode 100755 Tests/Rans/Askevein_hill/User_Mod/Tecplot_plane.f90 diff --git a/Tests/Rans/Askevein_hill/User_Mod/Cut_lines_ht_cp.f90 b/Tests/Rans/Askevein_hill/User_Mod/Cut_lines_ht_cp.f90 deleted file mode 100755 index d90685b8c..000000000 --- a/Tests/Rans/Askevein_hill/User_Mod/Cut_lines_ht_cp.f90 +++ /dev/null @@ -1,335 +0,0 @@ -!==============================================================================! - subroutine User_Mod_Cut_lines_ht_cp(flow, save_name) -!------------------------------------------------------------------------------! -! This is a prototype of a function for customized source for scalar. ! -! It is called from "Compute_Scalar" function, just before calling the ! -! linear solver. Both system matrix ("a_matrix") and right hand side ! -! vector ("b_vector") are sent should the user want to stabilize the ! -! system for always positive variables, for example. ! -!------------------------------------------------------------------------------! -!----------------------------------[Modules]-----------------------------------! - use Field_Mod, only: Field_Type, heat_transfer, heat_flux, & - density, viscosity, capacity, conductivity - use Grid_Mod, only: Grid_Type - use Var_Mod, only: Var_Type - use Matrix_Mod, only: Matrix_Type - use Bulk_Mod, only: Bulk_Type - use Rans_Mod - use Comm_Mod ! parallel stuff - use Name_Mod, only: problem_name - use Const_Mod ! constants -!------------------------------------------------------------------------------! - implicit none -!---------------------------------[Arguments]----------------------------------! - type(Field_Type), target :: flow - character(len=*) :: save_name -! type(Var_Type), target :: phi -! type(Matrix_Type) :: a_matrix -! real, dimension(:) :: b_vector -!----------------------------------[Locals]------------------------------------! - type(Grid_Type), pointer :: grid - type(Bulk_Type), pointer :: bulk - type(Var_Type), pointer :: u, v, w, t - real, pointer :: flux(:) - !------------------------------------------------------------------------------! - INTERFACE - LOGICAL FUNCTION Approx_Real(A,B,tol) - REAL :: A,B - REAL, OPTIONAL :: tol - END FUNCTION Approx_Real - END INTERFACE - -!-----------------------------[Parameters]-----------------------------! - - REAL :: Rad_2, Ufric, Rad_1 -!------------------------------[Calling]-------------------------------! - INTEGER :: Nprob, pl, c, i, count, k, c1, c2, s, N_hor, n - INTEGER :: Nfound - CHARACTER :: name_ht*12, name_cp*12, name_inflow*12 - character(len=80) :: res_name_cp, res_name_ht - character(len=80) :: store_name - REAL,ALLOCATABLE :: y_p(:), x_p(:), z_p(:), Ump(:), Vmp(:), Wmp(:), & - uup(:), vvp(:), wwp(:), & - uvp(:), uwp(:), vwp(:), & - Tmp(:), TTp(:), & - uTp(:), vTp(:), wTp(:), & - Ksgsp(:), & - var_1(:), var_2(:), var_3(:), Rad_mp(:), & - var_4(:), var_5(:), R_p(:) - INTEGER,ALLOCATABLE :: Np(:), Ncount(:) - REAL :: R, Urad_n, Utan_n, R1, R2, Urad, Utan, dummy - REAL :: xc_max, yc_max, zc_max, range - logical :: there -!==============================================================================! - ! Take aliases - grid => flow % pnt_grid - flux => flow % flux - bulk => flow % bulk - u => flow % u - v => flow % v - w => flow % w - t => flow % t - -!--------------------------------------------------------------------! - - ! Store the name - store_name = problem_name - problem_name = save_name - - call Name_File(0, res_name_ht, "-ht.dat") - call Name_File(0, res_name_cp, "-cp.dat") - - N_hor = 2 - Nfound = 0 - zc_max = 0.0 - range = HUGE - - allocate(x_p(N_hor)) - allocate(y_p(N_hor)) - allocate(R_p(N_hor)) - - do c = -1, -grid % n_bnd_cells, -1 - if( Grid_Mod_Bnd_Cond_Type(grid, c) .eq. WALL .or. & - Grid_Mod_Bnd_Cond_Type(grid, c) .eq. WALLFL) then - if(grid % zc(c) > zc_max) then - zc_max = grid % zc(c) - end if - end if - end do - - call Comm_Mod_Global_Max_Real(zc_max) - - if(this_proc< 2) write(*,*) 'Hill top height is ', zc_max - - call Comm_Mod_Wait - - xc_max = 0.0 - yc_max = 0.0 - - do c = 1, grid % n_cells - range = min(range, grid % delta(c)) - end do - - call Comm_Mod_Global_Min_Real(range) - - do c = -1, -grid % n_bnd_cells, -1 - if( Grid_Mod_Bnd_Cond_Type(grid, c) .eq. WALL .or. & - Grid_Mod_Bnd_Cond_Type(grid, c) .eq. WALLFL) then - if(Approx_Real( grid % zc(c), zc_max)) then - xc_max = grid % xc(c) - yc_max = grid % yc(c) - end if - end if - end do - - call Comm_Mod_Global_Sum_Real(xc_max) - call Comm_Mod_Global_Sum_Real(yc_max) - - if(this_proc< 2) write(*,*)'Coordinate of hill top (HT)', xc_max, yc_max , zc_max - - call Comm_Mod_Wait - - ! HT - x_p(1) = xc_max - y_p(1) = yc_max - - !CP ! - !distance between HT and CP is 400 m ! - -! x_p(2) = xc_max - 410.0*cos(78.0*3.14159/180.0) -! y_p(2) = yc_max - 410.0*sin(73.0*3.14159/180.0) - - x_p(2) = xc_max - .5*cos(78.0*3.14159/180.0) - y_p(2) = yc_max - .5*sin(73.0*3.14159/180.0) -!-Radius for averaging --adjust radius according to number of counts - R_p(1) = range - R_p(2) = range - - !------------------! - ! Read 1d file ! - !------------------! - inquire(file='Z_coordinate.dat', exist=there) - if(.not. there) then - if(this_proc < 2) then - print *, '===============================================================' - print *, 'In order to extract profiles and write them in ascii files' - print *, 'the code has to read cell-faces coordinates ' - print *, 'in wall-normal direction in the ascii file ''case_name.1d.''' - print *, 'The file format should be as follows:' - print *, '10 ! number of cells + 1' - print *, '1 0.0' - print *, '2 0.1' - print *, '3 0.2' - print *, '... ' - print *, '===============================================================' - end if - - ! Restore the name and return - problem_name = store_name - return - end if - - if(this_proc< 2) write(6, *) '# Now reading Z_coordinate.dat ' - - open(9, FILE='Z_coordinate.dat') -!---- write the number of probes - read(9,*) Nprob - allocate(z_p(Nprob)) - -!---- write the probe coordinates out - do pl=1,Nprob - read(9,*) dummy, z_p(pl) - end do - close(9) - - allocate(Np(Nprob)); Np = 0 - allocate(Ump(Nprob)); Ump = 0.0 - allocate(uwp(Nprob)); uwp = 0.0 - allocate(Rad_mp(Nprob)); Rad_mp = 0.0 - allocate(var_1(Nprob)); var_1 = 0.0 - allocate(var_2(Nprob)); var_2 = 0.0 - allocate(var_3(Nprob)); var_3 = 0.0 - allocate(var_4(Nprob)); var_4 = 0.0 - - allocate(Ncount(Nprob)); Ncount=0 - count = 0 - -!+++++++++++++++++++++++++++++! -! average the results ! -!+++++++++++++++++++++++++++++! - - do k = 1, 2 - do i = 1, Nprob-1 - 101 continue - Ncount(i) = 0 - do c = 1, grid % n_cells - grid % comm % n_buff_cells - Rad_1 = sqrt((grid % xc(c)-x_p(k))**2 + (grid % yc(c)-y_p(k))**2) - Rad_2 = grid % wall_dist(c) !grid % zc(c) - if(Rad_1 < R_p(k)) then - if(Rad_2 > z_p(i) .and. Rad_2 < z_p(i+1)) then - Rad_mp(i) = Rad_mp(i) + Rad_2 - Ncount(i) = Ncount(i) + 1 - end if - end if - end do - - call Comm_Mod_Wait - call Comm_Mod_Global_Sum_Int(Ncount(i)) - - if(Ncount(i) /= 0) then - Nfound = Ncount(i) !max(Nfound, Ncount(i)) - end if - - call Comm_Mod_Wait - - if(Nfound > 6) then - R_p(k) = R_p(k) * 0.9 - write(*,*) 'vise od 6', Nfound, R_p(k), k - goto 101 - else - continue - end if -! write(*,*) 'It found ', Nfound, ' points!' - - end do - end do - - do i = 1, Nprob - Ncount(i) = 0 - Rad_mp(i) = 0.0 - end do - -! write(*,*) x_p(2), x_p(1) -! write(*,*) y_p(2), y_p(1) - - do k = 1, 2 - do i = 1, Nprob-1 - do c = 1, grid % n_cells - Rad_1 = sqrt((grid % xc(c)-x_p(k))**2 + (grid % yc(c)-y_p(k))**2) - Rad_2 = grid % wall_dist(c) - if(Rad_1 < R_p(k)) then - if(Rad_2 > z_p(i) .and. Rad_2 < z_p(i+1)) then - Ncount(i) = Ncount(i) + 1 - Ump(i) = Ump(i) + u % n(c) - uwp(i) = uwp(i) + vis_t(c)*(u % z(c) + w % x(c)) - var_1(i) = var_1(i) + kin % n(c) - var_2(i) = var_2(i) + eps % n(c) - var_3(i) = var_3(i) + 5.64 + 1.41*log(grid % wall_dist(c)) - var_4(i) = var_4(i) + grid % wall_dist(c) - Rad_mp(i) = Rad_mp(i) + Rad_2 - end if - end if - end do - end do - -!---- average over all processors - do pl=1, Nprob - call Comm_Mod_Global_Sum_Int(Ncount(pl)) - call Comm_Mod_Global_Sum_Real(Ump(pl)) - call Comm_Mod_Global_Sum_Real(Rad_mp(pl)) - call Comm_Mod_Global_Sum_Real(uwp(pl)) - call Comm_Mod_Global_Sum_Real(var_1(pl)) - call Comm_Mod_Global_Sum_Real(var_2(pl)) - call Comm_Mod_Global_Sum_Real(var_3(pl)) - call Comm_Mod_Global_Sum_Real(var_4(pl)) - - count = count + Ncount(pl) - end do - - if(k == 1) then - open(k,FILE=res_name_ht) - else if(k == 2) then - open(k,FILE=res_name_cp) - end if - - do i = 1, Nprob - if(Ncount(i) /= 0) then - Ump(i) = Ump(i)/Ncount(i) - uwp(i) = uwp(i)/Ncount(i) - var_1(i) = var_1(i)/Ncount(i) - var_2(i) = var_2(i)/Ncount(i) - var_3(i) = var_3(i)/Ncount(i) - var_4(i) = var_4(i)/Ncount(i) - Rad_mp(i) = Rad_mp(i)/Ncount(i) - -!--------------Write variables ------------------------------------------------------------! - - write(k,'(3E15.7, I7)') var_4(i), (Ump(i) - var_3(i))/var_3(i), & - (var_1(i)/(8.9**2)), Ncount(i) - -!------------------------------------------------------------------------------------------! - - Ump(i) = 0.0 - uwp(i) = 0.0 - var_1(i) = 0.0 - var_2(i) = 0.0 - var_3(i) = 0.0 - var_4(i) = 0.0 - Rad_mp(i) = 0.0 - - Ncount(i) = 0 - end if - end do - - close(k) - end do - - deallocate(Np) - deallocate(z_p) - deallocate(Ump) - deallocate(uwp) - deallocate(var_1) - deallocate(var_2) - deallocate(var_3) - deallocate(var_4) - deallocate(Rad_mp) - deallocate(Ncount) - - call Comm_Mod_Wait - - if(this_proc< 2) write(*,*) 'Finished with user function for Askervein hill ! ' - - ! Restore the name - problem_name = store_name - - END SUBROUTINE User_Mod_Cut_lines_ht_cp diff --git a/Tests/Rans/Askevein_hill/User_Mod/Save_Results.f90 b/Tests/Rans/Askevein_hill/User_Mod/Save_Results.f90 deleted file mode 100755 index 3516d71c0..000000000 --- a/Tests/Rans/Askevein_hill/User_Mod/Save_Results.f90 +++ /dev/null @@ -1,20 +0,0 @@ -include '../User_Mod/Tecplot_plane.f90' -include '../User_Mod/Cut_lines_ht_cp.f90' - -!==============================================================================! - subroutine User_Mod_Save_Results(flow, save_name) -!------------------------------------------------------------------------------! -! Calls User_Impinging_Jet_Nu and User_Impinging_Jet_Profile functions. ! -!------------------------------------------------------------------------------! - use Field_Mod -!------------------------------------------------------------------------------! - implicit none -!---------------------------------[Arguments]----------------------------------! - type(Field_Type) :: flow - character(len=*) :: save_name -!==============================================================================! - - call User_Mod_Tecplot_plane (flow, save_name) - call User_Mod_Cut_lines_ht_cp(flow, save_name) - - end subroutine diff --git a/Tests/Rans/Askevein_hill/User_Mod/Tecplot_plane.f90 b/Tests/Rans/Askevein_hill/User_Mod/Tecplot_plane.f90 deleted file mode 100755 index ca8a032af..000000000 --- a/Tests/Rans/Askevein_hill/User_Mod/Tecplot_plane.f90 +++ /dev/null @@ -1,106 +0,0 @@ -!==============================================================================! - subroutine User_Mod_Tecplot_plane(flow, save_name) -!------------------------------------------------------------------------------! -! This is a prototype of a function for customized source for scalar. ! -! It is called from "Compute_Scalar" function, just before calling the ! -! linear solver. Both system matrix ("a_matrix") and right hand side ! -! vector ("b_vector") are sent should the user want to stabilize the ! -! system for always positive variables, for example. ! -!------------------------------------------------------------------------------! -!----------------------------------[Modules]-----------------------------------! - use Field_Mod, only: Field_Type, heat_flux - use Grid_Mod, only: Grid_Type - use Var_Mod, only: Var_Type - use Matrix_Mod, only: Matrix_Type - use Bulk_Mod, only: Bulk_Type - use Rans_Mod - use Comm_Mod ! parallel stuff - use Name_Mod, only: problem_name - use Const_Mod ! constants - -!------------------------------------------------------------------------------! - implicit none -!---------------------------------[Arguments]----------------------------------! - type(Field_Type), target :: flow - character(len=*) :: save_name -!-----------------------------------[Locals]-----------------------------------! - type(Grid_Type), pointer :: grid - type(Bulk_Type), pointer :: bulk - type(Var_Type), pointer :: u, v, w, t - real, pointer :: flux(:) - - INTEGER :: n - REAL :: x1_dis, x0_dis, y1_dis, y0_dis, z1_dis, z0_dis - REAL :: TauWup, TauWdown - REAL :: Utan, UnorSq, Unor, UtotSq, dely, Stot, R, UtanX -!-----------------------------------[Locals]-----------------------------------! - INTEGER :: c, s, c1, c2 - CHARACTER*80 :: namout, res_name, name_out_9 - CHARACTER*39 :: path - character(len=80) :: store_name -!------------------------------------------------------------------------------! - ! Take aliases - grid => flow % pnt_grid - flux => flow % flux - bulk => flow % bulk - u => flow % u - v => flow % v - w => flow % w - t => flow % t -!------------------------------------------------------------------------------! - - ! Store the name - store_name = problem_name - problem_name = save_name - - call Name_File(this_proc, res_name, "-tec360.dat") -! call Name_File(this_proc, name_out_9, '-tec360.dat') -! open(9, FILE='slice.dat') - -! if(this_proc < 2) write(*,*) '# Now reading the file: slice.dat ' -! read(9,*) z1_dis, z0_dis -! if(this_proc < 2) write(*,*) '# X:[ ', x0_dis, " ;", x1_dis, "]", & -! ' Y:[ ', y0_dis, " ;", y1_dis, "]", ' Z:[ ', z0_dis, " ;", z1_dis, "]" - -! inflowfile = 'tecplt-10m-xxxxxx-xx.dat' - -! write(inflowfile(12:17),'(I6.6)') n -! write(inflowfile(19:20),'(I2.2)') this_proc - open(500+this_proc, FILE=res_name) - - if( this_proc < 2 ) then - write(*,*) 'Capturing field.. ', res_name - end if - - - do s = 1, grid % n_faces - c1 = grid % faces_c(1,s) - c2 = grid % faces_c(2,s) - - if(c2 < 0) then - if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. WALL .or. & - Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. WALLFL) then - - write (500+this_proc,'(8E17.7E3)') grid % xc(c1), grid % yc(c1),& - grid % zc(c1), U % n(c1), V % n(c1), W % n(c1), Kin % n(c1), & - sqrt( (U % n(c1))**2+ (V % n(c1))**2+ (W % n(c1))**2) - - endif - endif - end do - - call Comm_Mod_Wait - - close(9) - close(500+this_proc) - - call Comm_Mod_Wait - - if ( this_proc < 2 ) then - write(*,*) 'It is done' - end if - - ! Restore the name - problem_name = store_name - - END SUBROUTINE From 677fce3c244cbb6a3af0038d96bb6fbc73d43a1f Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Mon, 6 May 2019 14:54:34 +0200 Subject: [PATCH 26/32] Limter for roughness coefficient is removed as it might lead to instabilities if massive separation occures and small y_plus is present in domain. --- Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 b/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 index c77158a3b..c38fc05c6 100644 --- a/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 +++ b/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90 @@ -126,7 +126,6 @@ subroutine Calculate_Vis_T_K_Eps(flow) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) - z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) y_plus(c1) = Y_Plus_Rough_Walls(u_tau(c1), & grid % wall_dist(c1), & kin_vis) From 8854ea8a09bc4fa63e45c990483a9bea87f74b35 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Mon, 6 May 2019 14:56:25 +0200 Subject: [PATCH 27/32] Limter for roughness coefficient is removed as it might lead to instabilities if massive separation occures and small y_plus is present in domain. --- Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 b/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 index b02d3904f..7d06fd9e0 100644 --- a/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 +++ b/Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90 @@ -119,7 +119,7 @@ subroutine Calculate_Vis_T_K_Eps_Zeta_F(flow) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) - z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) +! z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) y_plus(c1) = Y_Plus_Rough_Walls(u_tau(c1), & grid % wall_dist(c1), & kin_vis) From 521885566b0a3b8f607e51e4e68aacdf682aa834 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Mon, 6 May 2019 14:56:45 +0200 Subject: [PATCH 28/32] Limter for roughness coefficient is removed as it might lead to instabilities if massive separation occures and small y_plus is present in domain. --- Sources/Process/Turbulence/Source_Kin_K_Eps.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 b/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 index a40542a6b..ea612f024 100644 --- a/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 +++ b/Sources/Process/Turbulence/Source_Kin_K_Eps.f90 @@ -107,7 +107,6 @@ subroutine Source_Kin_K_Eps(flow, sol) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) - z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) u_tau(c1) = c_mu25 * sqrt(kin % n(c1)) y_plus(c1) = u_tau(c1) * (grid % wall_dist(c1) + z_o) & / kin_vis From 4121830d08d8bd746c60d29923240953b495624e Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Mon, 6 May 2019 14:56:56 +0200 Subject: [PATCH 29/32] Limter for roughness coefficient is removed as it might lead to instabilities if massive separation occures and small y_plus is present in domain. --- Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 b/Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 index 9d2fdaebb..1cb7c7952 100644 --- a/Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 +++ b/Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90 @@ -133,7 +133,6 @@ subroutine Source_Kin_K_Eps_Zeta_F(flow, sol) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) - z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) u_tau(c1) = c_mu25 * sqrt(kin % n(c1)) y_plus(c1) = Y_Plus_Rough_Walls(u_tau(c1), & grid % wall_dist(c1), kin_vis) From a2d8d4be19ee3cf703a5457ed72545a21dc2b14e Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Mon, 6 May 2019 14:57:08 +0200 Subject: [PATCH 30/32] Limter for roughness coefficient is removed as it might lead to instabilities if massive separation occures and small y_plus is present in domain. --- Sources/Process/Turbulence/Source_Eps_K_Eps.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/Sources/Process/Turbulence/Source_Eps_K_Eps.f90 b/Sources/Process/Turbulence/Source_Eps_K_Eps.f90 index 2db535060..35db95908 100644 --- a/Sources/Process/Turbulence/Source_Eps_K_Eps.f90 +++ b/Sources/Process/Turbulence/Source_Eps_K_Eps.f90 @@ -117,7 +117,6 @@ subroutine Source_Eps_K_Eps(flow, sol) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) - z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) eps % n(c1) = c_mu75 * kin % n(c1)**1.5 & / ((grid % wall_dist(c1) + z_o) * kappa) From 1cdcf1c0e7c67527e90dbb8b19e610029a8303a3 Mon Sep 17 00:00:00 2001 From: mhadziabdic Date: Mon, 6 May 2019 14:57:19 +0200 Subject: [PATCH 31/32] Limter for roughness coefficient is removed as it might lead to instabilities if massive separation occures and small y_plus is present in domain. --- Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 b/Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 index bca16b2e4..0f390196f 100644 --- a/Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 +++ b/Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90 @@ -115,7 +115,6 @@ subroutine Source_Eps_K_Eps_Zeta_F(flow, sol) if(rough_walls) then z_o = Roughness_Coefficient(grid, z_o_f(c1), c1) - z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) eps % n(c1) = c_mu75 * kin % n(c1)**1.5 / & ((grid % wall_dist(c1) + z_o) * kappa) From 581857fa42f8e66d1a8dab87d45497111a7d0fdc Mon Sep 17 00:00:00 2001 From: muhammmedaly Date: Tue, 12 Sep 2023 06:36:42 +0200 Subject: [PATCH 32/32] Addig a branch for Swarm_Mode --- Sources/Process/Swarm_Mod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sources/Process/Swarm_Mod.f90 b/Sources/Process/Swarm_Mod.f90 index 4abe695ad..56dbc9031 100644 --- a/Sources/Process/Swarm_Mod.f90 +++ b/Sources/Process/Swarm_Mod.f90 @@ -9,7 +9,7 @@ module Swarm_Mod !------------------------------------------------------------------------------! implicit none !==============================================================================! - +! Flag for branch "Mohamed_Particles" !-------------------! ! Particle type ! !-------------------!