Skip to content
Snippets Groups Projects
Commit a3a94db9 authored by Evgeny Mortikov's avatar Evgeny Mortikov
Browse files

turbulence closure setup update

parent 23144121
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
obl_k_epsilon.f90
+ 128
− 10
View file @ a3a94db9
......@@ -4,6 +4,11 @@ module obl_k_epsilon
! modules used
! --------------------------------------------------------------------------------
#ifdef USE_CONFIG_PARSER
use iso_c_binding, only: C_NULL_CHAR
use config_parser
#endif
use obl_math
use obl_grid
use obl_bc
......@@ -25,20 +30,24 @@ module obl_k_epsilon
public :: TKE_init, EPS_init
public :: TKE_bc, EPS_bc
public :: set_config_param
! --------------------------------------------------------------------------------
!> @brief k-epsilon parameters
type, public :: kepsilonParamType
real :: C_mu = 0.09 !< momentum stability function constant
real :: PrT0 = 1.25 !< turbulent Prandtl number
real :: TKE_min = 0.000001 !< minimal valur for TKE, [J/kg]
real :: eps_min = (0.000001)**(3.0/2.0) !< minimal value for epsilon (TKE dissipation rate), [W/kg]
real :: kappa = 0.4 !< von Karman constant
real :: C1_eps = 1.44 !< eps-eq.: production constant
real :: C2_eps = 1.92 !< eps-eq.: dissipation constant
real :: C3_eps = -0.40 !< eps-eq.: buoyancy constant (-0.4 stable stratification or 1.14 unstable stratification)
real :: Sch_TKE = 1.0 !< Schmidt number for TKE
real :: Sch_eps = 1.11 !< Scmidt number for eps
real :: ScT_TKE = 1.0 !< Schmidt number for TKE
real :: ScT_EPS = 1.11 !< Scmidt number for eps
real :: TKE_min = 0.000001 !< minimal valur for TKE, [J/kg]
real :: EPS_min = (0.000001)**(3.0/2.0) !< minimal value for epsilon (TKE dissipation rate), [W/kg]
real :: kappa = 0.4 !< von Karman constant
real :: PrT0 = 1.25 !< turbulent Prandtl number
real :: c_S2_min = 1e-5 !< min shear frequency, [(1/s)**2]
end type
......@@ -270,7 +279,7 @@ module obl_k_epsilon
real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
real :: eps(nz) !< TKE dissipation rate, [W/kg]
Kvisc(:) = Km(:) / param%Sch_TKE
Kvisc(:) = Km(:) / param%ScT_TKE
gamma = 0.0
......@@ -344,7 +353,7 @@ module obl_k_epsilon
real :: eps(nz) !< TKE dissipation rate, [W/kg]
real :: TKE(nz) !< TKE, [J/kg]
Kvisc(:) = Km(:) / param%Sch_eps
Kvisc(:) = Km(:) / param%ScT_EPS
gamma = 0.0
......@@ -416,7 +425,7 @@ module obl_k_epsilon
real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
real :: eps(nz) !< TKE dissipation rate, [W/kg]
Kvisc(:) = Km(:) / param%Sch_TKE
Kvisc(:) = Km(:) / param%ScT_TKE
gamma = 0.0
......@@ -491,7 +500,7 @@ module obl_k_epsilon
real :: TKE(nz) !< TKE, [J/kg]
Kvisc(:) = Km(:) / param%Sch_eps
Kvisc(:) = Km(:) / param%ScT_EPS
gamma = 0.0
......@@ -540,4 +549,113 @@ module obl_k_epsilon
end subroutine solve_eps_eq_semiimplicit
! --------------------------------------------------------------------------------
subroutine set_config_param(param, tag, ierr)
!< @brief set turbulence closure parameters
! ----------------------------------------------------------------------------
type(kepsilonParamType), intent(out) :: param
integer, intent(out) :: ierr
character(len = *), intent(in) :: tag
integer :: status
! --------------------------------------------------------------------------------
ierr = 0 ! = OK
#ifdef USE_CONFIG_PARSER
call c_config_is_varname(trim(tag)//".C_mu"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".C_mu"//C_NULL_CHAR, param%C_mu, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".C1_eps"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".C1_eps"//C_NULL_CHAR, param%C1_eps, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".C2_eps"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".C2_eps"//C_NULL_CHAR, param%C2_eps, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".C3_eps"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".C3_eps"//C_NULL_CHAR, param%C3_eps, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".ScT_TKE"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".ScT_TKE"//C_NULL_CHAR, param%ScT_TKE, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".ScT_EPS"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".ScT_EPS"//C_NULL_CHAR, param%ScT_EPS, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".TKE_min"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".TKE_min"//C_NULL_CHAR, param%TKE_min, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".EPS_min"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".EPS_min"//C_NULL_CHAR, param%EPS_min, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".kappa"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".kappa"//C_NULL_CHAR, param%kappa, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".PrT0"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".PrT0"//C_NULL_CHAR, param%PrT0, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".c_S2_min"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".c_S2_min"//C_NULL_CHAR, param%c_S2_min, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
#else
!> *: just skipping setup without configuration file
#endif
end subroutine
end module
\ No newline at end of file
obl_main.f90
+ 22
− 0
View file @ a3a94db9
......@@ -26,14 +26,17 @@ program obl_main
use obl_diag
use obl_bc
use obl_pph, only: pphParamType, &
set_config_param_pph => set_config_param, &
define_stability_functions_pph => define_stability_functions, &
init_turbulence_closure_pph => init_turbulence_closure, &
advance_turbulence_closure_pph => advance_turbulence_closure
use obl_pph_dyn, only: pphDynParamType, &
set_config_param_pph_dyn => set_config_param, &
define_stability_functions_pph_dyn => define_stability_functions, &
init_turbulence_closure_pph_dyn => init_turbulence_closure, &
advance_turbulence_closure_pph_dyn => advance_turbulence_closure
use obl_k_epsilon, only: kepsilonParamType, &
set_config_param_k_epsilon => set_config_param, &
define_stability_functions_k_epsilon => define_stability_functions, &
init_turbulence_closure_k_epsilon => init_turbulence_closure, &
advance_turbulence_closure_k_epsilon => advance_turbulence_closure, &
......@@ -238,6 +241,25 @@ program obl_main
endif
! ----------------------------------------------------------------------------
!< setting turbulence closure parameters
! ----------------------------------------------------------------------------
call set_config_param_pph(param_pph, "turbulence_model.pph", ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set ''pph'' parameters '
return
endif
call set_config_param_pph_dyn(param_pph_dyn, "turbulence_model.pph_dyn", ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set ''pph-dyn'' parameters '
return
endif
call set_config_param_k_epsilon(param_k_epsilon, "turbulence_model.k_epsilon", ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set ''k-epsilon'' parameters '
return
endif
! ----------------------------------------------------------------------------
!< initialization of main fields
! ----------------------------------------------------------------------------
call set_initial_conditions(grid, obl_config_id, ierr)
......
obl_pph.f90
+ 102
− 0
View file @ a3a94db9
#include "obl_def.fi"
module obl_pph
!< @brief standard pacanowski-philander scheme
! --------------------------------------------------------------------------------
! modules used
! --------------------------------------------------------------------------------
#ifdef USE_CONFIG_PARSER
use iso_c_binding, only: C_NULL_CHAR
use config_parser
#endif
use obl_grid
use obl_bc
use obl_turb_common
......@@ -21,6 +28,8 @@ module obl_pph
public :: get_eddy_viscosity, get_eddy_diffusivity
public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec
public :: set_config_param
! --------------------------------------------------------------------------------
!> @brief Pacanowski-Philander parameters
......@@ -155,4 +164,97 @@ module obl_pph
end do
end subroutine
! --------------------------------------------------------------------------------
subroutine set_config_param(param, tag, ierr)
!< @brief set turbulence closure parameters
! ----------------------------------------------------------------------------
type(pphParamType), intent(out) :: param
integer, intent(out) :: ierr
character(len = *), intent(in) :: tag
integer :: status
! --------------------------------------------------------------------------------
ierr = 0 ! = OK
#ifdef USE_CONFIG_PARSER
call c_config_is_varname(trim(tag)//".Km_0"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Km_0"//C_NULL_CHAR, param%Km_0, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Kh_0"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Kh_0"//C_NULL_CHAR, param%Kh_0, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".alpha"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".alpha"//C_NULL_CHAR, param%alpha, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".n"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".n"//C_NULL_CHAR, param%n, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Km_b"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Km_b"//C_NULL_CHAR, param%Km_b, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Kh_b"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Kh_b"//C_NULL_CHAR, param%Kh_b, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".kappa"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".kappa"//C_NULL_CHAR, param%kappa, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".PrT0"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".PrT0"//C_NULL_CHAR, param%PrT0, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".c_S2_min"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".c_S2_min"//C_NULL_CHAR, param%c_S2_min, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
#else
!> *: just skipping setup without configuration file
#endif
end subroutine
end module
obl_pph_dyn.f90
+ 94
− 2
View file @ a3a94db9
......@@ -5,6 +5,11 @@ module obl_pph_dyn
! modules used
! --------------------------------------------------------------------------------
#ifdef USE_CONFIG_PARSER
use iso_c_binding, only: C_NULL_CHAR
use config_parser
#endif
use obl_grid
use obl_bc
use obl_turb_common
......@@ -23,6 +28,8 @@ module obl_pph_dyn
public :: get_eddy_viscosity, get_eddy_diffusivity
public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec
public :: set_config_param
! --------------------------------------------------------------------------------
!> @brief Pacanowski-Philander dynamic parameters
......@@ -32,11 +39,11 @@ module obl_pph_dyn
real :: Km_b = 0.0001 !< background eddy viscosity, [m**2 / s], *: INMCM: 5 * 10**(-6)
real :: Kh_b = 0.00001 !< background eddy diffusivity, [m**2 / s]
real :: C_l = 4.0 / 9.0 !< length scale tuning constant, [n/d], *: = (4.0 / pi^2) in Obukhov length scale
real :: kappa = 0.4 !< von Karman constant
real :: PrT0 = 1.0 !< turbulent Prandtl number, *: probably redundant
real :: C_l = 4.0 / 9.0 !< length scale tuning constant, [n/d], *: = (4.0 / pi^2) in Obukhov length scale
real :: c_S2_min = 1e-5 !< min shear frequency, [(1/s)**2]
end type
......@@ -187,4 +194,89 @@ module obl_pph_dyn
end do
end subroutine
! --------------------------------------------------------------------------------
subroutine set_config_param(param, tag, ierr)
!< @brief set turbulence closure parameters
! ----------------------------------------------------------------------------
type(pphDynParamType), intent(out) :: param
integer, intent(out) :: ierr
character(len = *), intent(in) :: tag
integer :: status
! --------------------------------------------------------------------------------
ierr = 0 ! = OK
#ifdef USE_CONFIG_PARSER
call c_config_is_varname(trim(tag)//".alpha"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".alpha"//C_NULL_CHAR, param%alpha, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".n"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".n"//C_NULL_CHAR, param%n, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Km_b"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Km_b"//C_NULL_CHAR, param%Km_b, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Kh_b"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Kh_b"//C_NULL_CHAR, param%Kh_b, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".C_l"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".C_l"//C_NULL_CHAR, param%C_l, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".kappa"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".kappa"//C_NULL_CHAR, param%kappa, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".PrT0"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".PrT0"//C_NULL_CHAR, param%PrT0, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".c_S2_min"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".c_S2_min"//C_NULL_CHAR, param%c_S2_min, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
#else
!> *: just skipping setup without configuration file
#endif
end subroutine
end module
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment