45
45
use sparsity_patterns_meshes
52
54
! locally allocatad fields
53
type(scalar_field), save :: tke_old, tke_src_old
54
real, save :: c_mu, c_eps_1, c_eps_2, sigma_eps, sigma_k, lmax
55
real, save :: fields_min = 1.e-10
56
character(len=FIELD_NAME_LEN), save :: src_abs
55
real, save :: fields_min = 1.0e-11
56
logical, save :: low_Re = .false.
58
public :: keps_init, keps_cleanup, keps_tke, keps_eps, keps_eddyvisc, keps_bcs, keps_adapt_mesh, keps_check_options
58
public :: keps_advdif_diagnostics, keps_momentum_diagnostics, keps_bcs, &
59
& k_epsilon_check_options, tensor_inner_product
61
! - Init in populate_state.
62
! - call keps_tke (which calculates production and sets source/absorption/diffusivity for solve).
63
! - call keps_eps (which sets source/absorption/diffusivity for solve).
64
! - After keps_eps solve, keps_eddyvisc recalculates the eddy viscosity and adds it to the viscosity field.
65
! - Wall functions are added to selected boundaries in keps_bcs and wall_functions.
66
! - keps_adapt_options repopulates the fields after an adapt.
67
! - When done, clean-up.
62
! - call diagnostics to obtain source terms and calculate eddy viscosity
63
! - after each scalar field solve recalculates the eddy viscosity
64
! - wall functions are added to selected boundaries in keps_bcs and wall_functions
72
! initialise parameters based on options
75
subroutine keps_init(state)
77
type(state_type), intent(inout) :: state
78
type(scalar_field), pointer :: Field
79
character(len=OPTION_PATH_LEN) :: keps_path
82
ewrite(1,*)'Now in k_epsilon turbulence model - keps_init'
83
keps_path = trim(state%option_path)//"/subgridscale_parameterisations/k-epsilon"
84
ewrite(2,*)'keps_path: ', trim(keps_path)
86
! Allocate the temporary, module-level variables
87
call keps_allocate_fields(state)
89
! Are source and absorption terms implicit/explicit?
90
call get_option(trim(keps_path)//'/source_absorption', src_abs)
92
! What is the maximum physical lengthscale in the domain?
93
call get_option(trim(keps_path)//'/lengthscale_limit', lmax)
95
! Get the 5 model constants
96
call get_option(trim(keps_path)//'/C_mu', C_mu, default = 0.09)
97
call get_option(trim(keps_path)//'/C_eps_1', c_eps_1, default = 1.44)
98
call get_option(trim(keps_path)//'C_eps_2', c_eps_2, default = 1.92)
99
call get_option(trim(keps_path)//'/sigma_k', sigma_k, default = 1.0)
100
call get_option(trim(keps_path)//'sigma_eps', sigma_eps, default = 1.3)
102
! Get background viscosity
103
call get_option(trim(keps_path)//"/tensor_field::BackgroundViscosity/prescribed/&
104
&value::WholeMesh/isotropic/constant", visc)
106
! initialise eddy viscosity
107
if (have_option(trim(keps_path)//"/scalar_field::ScalarEddyViscosity/diagnostic")) then
108
Field => extract_scalar_field(state, "ScalarEddyViscosity")
109
call set(Field, visc)
112
! initialise lengthscale
113
Field => extract_scalar_field(state, "LengthScale")
116
ewrite(2,*) "k-epsilon parameters"
117
ewrite(2,*) "--------------------------------------------"
118
ewrite(2,*) "c_mu: ", c_mu
119
ewrite(2,*) "c_eps_1: ", c_eps_1
120
ewrite(2,*) "c_eps_2: ", c_eps_2
121
ewrite(2,*) "sigma_k: ", sigma_k
122
ewrite(2,*) "sigma_eps: ",sigma_eps
123
ewrite(2,*) "fields_min: ",fields_min
124
ewrite(2,*) "background visc: ", visc
125
ewrite(2,*) "implicit/explicit source/absorption terms: ", trim(src_abs)
126
ewrite(2,*) "max turbulence lengthscale: ", lmax
127
ewrite(2,*) "--------------------------------------------"
129
end subroutine keps_init
131
!------------------------------------------------------------------------------
133
subroutine keps_tke(state)
135
type(state_type), intent(inout) :: state
136
type(scalar_field), pointer :: src_kk, abs_kk, kk, eps, EV, lumped_mass
137
type(scalar_field) :: src_rhs, abs_rhs, prescribed_src_kk, prescribed_abs_kk
138
type(vector_field), pointer :: positions, nu, u
139
type(tensor_field), pointer :: kk_diff
140
type(element_type), pointer :: shape_kk
142
integer, pointer, dimension(:) :: nodes_kk
144
real, allocatable, dimension(:) :: detwei, strain_ngi, rhs_addto
145
real, allocatable, dimension(:,:,:):: dshape_kk
146
logical :: prescribed_src, prescribed_abs
148
ewrite(1,*) "In keps_tke"
150
positions => extract_vector_field(state, "Coordinate")
151
nu => extract_vector_field(state, "NonlinearVelocity")
152
u => extract_vector_field(state, "Velocity")
153
src_kk => extract_scalar_field(state, "TurbulentKineticEnergySource")
154
abs_kk => extract_scalar_field(state, "TurbulentKineticEnergyAbsorption")
155
kk_diff => extract_tensor_field(state, "TurbulentKineticEnergyDiffusivity")
156
kk => extract_scalar_field(state, "TurbulentKineticEnergy")
157
eps => extract_scalar_field(state, "TurbulentDissipation")
158
EV => extract_scalar_field(state, "ScalarEddyViscosity")
160
! Set copy of old kk for eps solve
161
call set(tke_old, kk)
163
prescribed_src=.false.
164
prescribed_abs=.false.
165
call allocate(src_rhs, kk%mesh, name="KKSRCRHS")
166
call allocate(abs_rhs, kk%mesh, name="KKABSRHS")
167
call zero(src_rhs); call zero(abs_rhs)
170
do ele = 1, ele_count(kk)
171
shape_kk => ele_shape(kk, ele)
172
nodes_kk => ele_nodes(kk, ele)
174
allocate(dshape_kk (size(nodes_kk), ele_ngi(kk, ele), positions%dim))
175
allocate(detwei (ele_ngi(kk, ele)))
176
allocate(strain_ngi (ele_ngi(kk, ele)))
177
allocate(rhs_addto(ele_loc(kk, ele)))
178
call transform_to_physical( positions, ele, shape_kk, dshape=dshape_kk, detwei=detwei )
180
! Calculate TKE production at ngi using strain rate (double_dot_product) function
181
strain_ngi = double_dot_product(dshape_kk, ele_val(u, ele) )
184
rhs_addto = shape_rhs(shape_kk, detwei*strain_ngi*ele_val_at_quad(EV, ele))
185
call addto(src_rhs, nodes_kk, rhs_addto)
188
rhs_addto = shape_rhs(shape_kk, detwei*ele_val_at_quad(eps,ele)/ele_val_at_quad(kk,ele))
189
call addto(abs_rhs, nodes_kk, rhs_addto)
191
deallocate(dshape_kk, detwei, strain_ngi, rhs_addto)
194
lumped_mass => get_lumped_mass(state, kk%mesh)
196
! This allows user-specified source and absorption terms, so that an MMS test can be set up.
197
prescribed_src = (have_option("/material_phase[0]/subgridscale_parameterisations/k-epsilon/&
198
&scalar_field::TurbulentKineticEnergy/prognostic/&
199
&scalar_field::Source/prescribed"))
201
if(prescribed_src) then
202
ewrite(2,*) "Prescribed k source"
203
call allocate(prescribed_src_kk, kk%mesh, name="PRSKKSRC")
204
call zero(prescribed_src_kk)
205
call set(prescribed_src_kk, src_kk)
206
ewrite_minmax(prescribed_src_kk)
210
prescribed_abs = (have_option("/material_phase[0]/subgridscale_parameterisations/k-epsilon/&
211
&scalar_field::TurbulentKineticEnergy/prognostic/&
212
&scalar_field::Absorption/prescribed"))
214
if(prescribed_abs) then
215
ewrite(2,*) "Prescribed k absorption"
216
call allocate(prescribed_abs_kk, kk%mesh, name="PRSKKABS")
217
call zero(prescribed_abs_kk)
218
call set(prescribed_abs_kk, abs_kk)
219
ewrite_minmax(prescribed_abs_kk)
223
ewrite(2,*) "Calculating k source and absorption"
224
do i = 1, node_count(kk)
225
select case (src_abs)
227
call set(src_kk, i, node_val(src_rhs,i)/node_val(lumped_mass,i))
228
call set(abs_kk, i, node_val(abs_rhs,i)/node_val(lumped_mass,i))
230
residual = (node_val(abs_rhs,i) - node_val(src_rhs,i))/node_val(lumped_mass,i)
231
call set(src_kk, i, -min(0.0, residual) )
232
call set(abs_kk, i, max(0.0, residual) )
234
FLAbort("Invalid implicitness option for k")
238
! Set copy of old kk for eps solve
239
call set(tke_src_old, src_kk)
241
if(prescribed_src) then
242
ewrite_minmax(src_kk)
243
! Set copy of old kk for eps solve
244
call set(tke_src_old, src_kk)
245
call addto(src_kk, prescribed_src_kk)
246
ewrite_minmax(src_kk)
247
call deallocate(prescribed_src_kk)
250
if(prescribed_abs) then
251
call set(abs_kk, prescribed_abs_kk)
252
call deallocate(prescribed_abs_kk)
255
call deallocate(src_rhs); call deallocate(abs_rhs)
257
! Set diffusivity for k equation.
259
do i = 1, kk_diff%dim(1)
260
call set(kk_diff, i, i, EV, scale=1. / sigma_k)
263
ewrite_minmax(kk_diff)
264
ewrite_minmax(tke_old)
266
ewrite_minmax(src_kk)
267
ewrite_minmax(tke_src_old)
268
ewrite_minmax(abs_kk)
270
end subroutine keps_tke
272
!----------------------------------------------------------------------------------
274
subroutine keps_eps(state)
276
type(state_type), intent(inout) :: state
277
type(scalar_field), pointer :: src_eps, src_kk, abs_eps, eps, EV, lumped_mass
278
type(scalar_field) :: src_rhs, abs_rhs, prescribed_src_eps, prescribed_abs_eps
279
type(vector_field), pointer :: positions
280
type(tensor_field), pointer :: eps_diff
281
type(element_type), pointer :: shape_eps
284
real, allocatable, dimension(:) :: detwei, rhs_addto
285
integer, pointer, dimension(:) :: nodes_eps
286
logical :: prescribed_src, prescribed_abs
288
ewrite(1,*) "In keps_eps"
289
eps => extract_scalar_field(state, "TurbulentDissipation")
290
src_eps => extract_scalar_field(state, "TurbulentDissipationSource")
291
src_kk => extract_scalar_field(state, "TurbulentKineticEnergySource")
292
abs_eps => extract_scalar_field(state, "TurbulentDissipationAbsorption")
293
eps_diff => extract_tensor_field(state, "TurbulentDissipationDiffusivity")
294
EV => extract_scalar_field(state, "ScalarEddyViscosity")
295
positions => extract_vector_field(state, "Coordinate")
297
call allocate(src_rhs, eps%mesh, name="EPSSRCRHS")
298
call allocate(abs_rhs, eps%mesh, name="EPSABSRHS")
299
call zero(src_rhs); call zero(abs_rhs)
302
do ele = 1, ele_count(eps)
303
shape_eps => ele_shape(eps, ele)
304
nodes_eps => ele_nodes(eps, ele)
306
allocate(detwei (ele_ngi(eps, ele)))
307
allocate(rhs_addto(ele_loc(eps, ele)))
308
call transform_to_physical(positions, ele, detwei=detwei)
311
rhs_addto = shape_rhs(shape_eps, detwei*c_eps_1*ele_val_at_quad(eps,ele)/ &
312
ele_val_at_quad(tke_old,ele)*ele_val_at_quad(tke_src_old,ele))
313
call addto(src_rhs, nodes_eps, rhs_addto)
316
rhs_addto = shape_rhs(shape_eps, detwei*c_eps_2*ele_val_at_quad(eps,ele)/ &
317
ele_val_at_quad(tke_old,ele))
318
call addto(abs_rhs, nodes_eps, rhs_addto)
320
deallocate(detwei, rhs_addto)
323
lumped_mass => get_lumped_mass(state, eps%mesh)
325
! This allows user-specified source and absorption terms, so that an MMS test can be set up.
326
prescribed_src = (have_option("/material_phase[0]/subgridscale_parameterisations/k-epsilon/&
327
&scalar_field::TurbulentKineticEnergy/prognostic/&
328
&scalar_field::Source/prescribed"))
330
if(prescribed_src) then
331
ewrite(2,*) "Prescribed eps source"
332
call allocate(prescribed_src_eps, eps%mesh, name="PRSEPSSRC")
333
call zero(prescribed_src_eps)
334
call set(prescribed_src_eps, src_eps)
335
ewrite_minmax(prescribed_src_eps)
339
prescribed_abs = (have_option("/material_phase[0]/subgridscale_parameterisations/k-epsilon/&
340
&scalar_field::TurbulentDissipation/prognostic/scalar_field::Absorption/prescribed"))
342
if(prescribed_abs) then
343
ewrite(2,*) "Prescribed eps absorption"
344
call allocate(prescribed_abs_eps, eps%mesh, name="PRSEPSABS")
345
call zero(prescribed_abs_eps)
346
call set(prescribed_abs_eps, abs_eps)
347
ewrite_minmax(prescribed_abs_eps)
351
ewrite(2,*) "Calculating epsilon source and absorption"
352
do i = 1, node_count(eps)
353
select case (src_abs)
355
call set(src_eps, i, node_val(src_rhs,i)/node_val(lumped_mass,i))
356
call set(abs_eps, i, node_val(abs_rhs,i)/node_val(lumped_mass,i))
358
residual = (node_val(abs_rhs,i) - node_val(src_rhs,i))/node_val(lumped_mass,i)
359
call set(src_eps, i, -min(0.0, residual) )
360
call set(abs_eps, i, max(0.0, residual) )
362
FLAbort("Invalid implicitness option for epsilon")
366
if(prescribed_src) then
367
ewrite_minmax(src_eps)
368
! Set copy of old eps for eps solve
369
call set(tke_src_old, src_eps)
370
call addto(src_eps, prescribed_src_eps)
371
ewrite_minmax(src_eps)
372
call deallocate(prescribed_src_eps)
375
if(prescribed_abs) then
376
call set(abs_eps, prescribed_abs_eps)
377
call deallocate(prescribed_abs_eps)
380
call deallocate(src_rhs); call deallocate(abs_rhs)
382
! Set diffusivity for Eps
384
do i = 1, eps_diff%dim(1)
385
call set(eps_diff, i,i, EV, scale=1./sigma_eps)
388
ewrite_minmax(eps_diff)
389
ewrite_minmax(tke_old)
390
ewrite_minmax(src_kk)
392
ewrite_minmax(src_eps)
393
ewrite_minmax(abs_eps)
395
end subroutine keps_eps
398
! eddyvisc calculates the lengthscale, and then the eddy viscosity.
68
subroutine keps_advdif_diagnostics(state)
70
type(state_type), intent(inout) :: state
72
call keps_damping_functions(state, advdif=.true.)
73
call keps_eddyvisc(state, advdif=.true.)
74
call keps_diffusion(state)
75
call keps_tracer_diffusion(state)
76
call keps_calculate_rhs(state)
78
end subroutine keps_advdif_diagnostics
80
subroutine keps_momentum_diagnostics(state)
82
type(state_type), intent(inout) :: state
84
call keps_damping_functions(state, advdif=.false.)
85
call keps_eddyvisc(state, advdif=.false.)
87
end subroutine keps_momentum_diagnostics
89
!--------------------------------------------------------------------------------!
91
subroutine keps_damping_functions(state, advdif)
93
type(state_type), intent(in) :: state
94
logical, intent(in) :: advdif
96
type(scalar_field), pointer :: f_1, f_2, f_mu, y, dummydensity, density
97
type(scalar_field) :: k, eps
98
type(tensor_field), pointer :: bg_visc
100
real :: f_mu_val, f_1_val, f_2_val, Re_T, R_y, fields_max
101
character(len=FIELD_NAME_LEN) :: equation_type
102
character(len=OPTION_PATH_LEN) :: option_path
104
ewrite(1,*) 'in keps_damping_functions'
106
option_path = trim(state%option_path)//'/subgridscale_parameterisations/k-epsilon/'
108
f_1 => extract_scalar_field(state, "f_1")
109
f_2 => extract_scalar_field(state, "f_2")
110
f_mu => extract_scalar_field(state, "f_mu")
112
! initialise low_Re damping functions
117
call get_option(trim(state%option_path)// &
118
& "/subgridscale_parameterisations/k-epsilon/max_damping_value", fields_max)
120
! Low Reynolds damping functions
121
! Check for low reynolds boundary condition and calculate damping functions
122
! Lam-Bremhorst model (Wilcox 1998 - Turbulence modelling for CFD)
124
have_option(trim(option_path)//"debugging_options/enable_lowRe_damping")) then
126
bg_visc => extract_tensor_field(state, "BackgroundViscosity")
127
y => extract_scalar_field(state, "DistanceToWall", stat = stat)
129
FLAbort("I need the distance to the wall - enable a DistanceToWall field")
132
call time_averaged_value(state, k, 'TurbulentKineticEnergy', advdif, option_path)
133
call time_averaged_value(state, eps, 'TurbulentDissipation', advdif, option_path)
135
allocate(dummydensity)
136
call allocate(dummydensity, f_1%mesh, "DummyDensity", field_type=FIELD_TYPE_CONSTANT)
137
call set(dummydensity, 1.0)
138
dummydensity%option_path = ""
140
! Depending on the equation type, extract the density or set it to some dummy field allocated above
141
call get_option(trim(state%option_path)//&
142
"/vector_field::Velocity/prognostic/equation[0]/name", equation_type)
143
select case(equation_type)
144
case("LinearMomentum")
145
density=>extract_scalar_field(state, "Density")
147
density=>dummydensity
149
density=>dummydensity
151
! developer error... out of sync options input and code
152
FLAbort("Unknown equation type for velocity")
155
node_loop: do node = 1, node_count(k)
157
! calc of damping values with error catching
158
if (node_val(bg_visc,1,1,node) <= fields_min) then
162
else if (node_val(eps,node) <= fields_min) then
163
R_y = (node_val(density,node) * node_val(k,node)**0.5 * node_val(y,node)) / &
164
node_val(bg_visc,1,1,node)
166
f_mu_val = (1.0 - exp(- 0.0165*R_y))**2.0
167
f_1_val = (0.05/node_val(f_mu,node))**3.0 + 1.0
170
Re_T = (node_val(density,node) * node_val(k,node)**2.0) / &
171
(node_val(eps,node) * node_val(bg_visc,1,1,node))
172
R_y = (node_val(density,node) * node_val(k,node)**0.5 * node_val(y,node)) / &
173
node_val(bg_visc,1,1,node)
175
f_mu_val = (1.0 - exp(- 0.0165*R_y))**2.0 * (20.5/Re_T + 1.0)
176
f_1_val = (0.05/f_mu_val)**3.0 + 1.0
177
f_2_val = 1.0 - exp(- Re_T**2.0)
180
! limit values of damping functions
181
call set(f_mu, node, min(f_mu_val, 1.0))
182
call set(f_1, node, min(f_1_val, fields_max))
183
call set(f_2, node, min(f_2_val, fields_max))
189
call deallocate(dummydensity)
190
deallocate(dummydensity)
194
end subroutine keps_damping_functions
196
!------------------------------------------------------------------------------!
198
subroutine keps_calculate_rhs(state)
200
type(state_type), intent(inout) :: state
202
type(scalar_field), dimension(3) :: src_abs_terms
203
type(scalar_field), dimension(2) :: fields
204
type(scalar_field), pointer :: src, abs, f_1, f_2, debug
205
type(scalar_field) :: src_to_abs
206
type(vector_field), pointer :: x, u, g
207
type(scalar_field), pointer :: dummydensity, density, buoyancy_density, scalar_eddy_visc
208
integer :: i, ele, term, stat
209
real :: g_magnitude, c_eps_1, c_eps_2, sigma_p
210
logical :: have_buoyancy_turbulence = .true., lump_mass
211
character(len=OPTION_PATH_LEN) :: option_path
212
character(len=FIELD_NAME_LEN), dimension(2) :: field_names
213
character(len=FIELD_NAME_LEN) :: equation_type, implementation
215
option_path = trim(state%option_path)//'/subgridscale_parameterisations/k-epsilon/'
217
if (.not. have_option(trim(option_path))) then
221
ewrite(1,*) 'In calculate k-epsilon rhs'
223
! get model constants
224
call get_option(trim(option_path)//'/C_eps_1', c_eps_1, default = 1.44)
225
call get_option(trim(option_path)//'/C_eps_2', c_eps_2, default = 1.92)
226
call get_option(trim(option_path)//'/sigma_p', sigma_p, default = 1.0)
229
x => extract_vector_field(state, "Coordinate")
230
u => extract_vector_field(state, "NonlinearVelocity")
231
scalar_eddy_visc => extract_scalar_field(state, "ScalarEddyViscosity")
232
f_1 => extract_scalar_field(state, "f_1")
233
f_2 => extract_scalar_field(state, "f_2")
234
g => extract_vector_field(state, "GravityDirection", stat)
235
call get_option('/physical_parameters/gravity/magnitude', g_magnitude, stat)
237
have_buoyancy_turbulence = .false.
240
allocate(dummydensity)
241
call allocate(dummydensity, X%mesh, "DummyDensity", field_type=FIELD_TYPE_CONSTANT)
242
call set(dummydensity, 1.0)
243
dummydensity%option_path = ""
245
! Depending on the equation type, extract the density or set it to some dummy field allocated above
246
call get_option(trim(state%option_path)//&
247
"/vector_field::Velocity/prognostic/equation[0]/name", equation_type)
248
select case(equation_type)
249
case("LinearMomentum")
250
density=>extract_scalar_field(state, "Density")
252
density=>dummydensity
254
density=>dummydensity
256
! developer error... out of sync options input and code
257
FLAbort("Unknown equation type for velocity")
260
if(have_buoyancy_turbulence) then
261
buoyancy_density => extract_scalar_field(state, 'VelocityBuoyancyDensity')
264
field_names(1) = 'TurbulentKineticEnergy'
265
field_names(2) = 'TurbulentDissipation'
267
field_loop: do i = 1, 2
268
if (have_option(trim(option_path)//'scalar_field::'// &
269
trim(field_names(i))//'/prescribed')) then
273
!-----------------------------------------------------------------------------------
276
src => extract_scalar_field(state, trim(field_names(i))//"Source")
277
abs => extract_scalar_field(state, trim(field_names(i))//"Absorption")
279
call time_averaged_value(state, fields(1), trim(field_names(i)), .true., option_path)
280
call time_averaged_value(state, fields(2), trim(field_names(3-i)), .true., option_path)
282
call allocate(src_abs_terms(1), fields(1)%mesh, name="production_term")
283
call allocate(src_abs_terms(2), fields(1)%mesh, name="destruction_term")
284
call allocate(src_abs_terms(3), fields(1)%mesh, name="buoyancy_term")
285
call zero(src_abs_terms(1)); call zero(src_abs_terms(2)); call zero(src_abs_terms(3))
286
call zero(src); call zero(abs)
288
!-----------------------------------------------------------------------------------
291
do ele = 1, ele_count(fields(1))
292
call assemble_rhs_ele(src_abs_terms, fields(i), fields(3-i), scalar_eddy_visc, u, &
293
density, buoyancy_density, have_buoyancy_turbulence, g, g_magnitude, x, &
294
c_eps_1, c_eps_2, sigma_p, f_1, f_2, ele, i)
297
! For non-DG we apply inverse mass globally
298
if(continuity(fields(1))>=0) then
299
lump_mass = have_option(trim(option_path)//'mass_terms/lump_mass')
301
call solve_cg_inv_mass(state, src_abs_terms(term), lump_mass, option_path)
304
!-----------------------------------------------------------------------------------
306
! Source disabling for debugging purposes
308
if(have_option(trim(option_path)//'debugging_options/disable_'//&
309
trim(src_abs_terms(term)%name))) then
310
call zero(src_abs_terms(term))
313
!-----------------------------------------------------------------------------------
315
! Produce debugging output
317
debug => extract_scalar_field(state, &
318
trim(field_names(i))//"_"//trim(src_abs_terms(term)%name), stat)
320
call set(debug, src_abs_terms(term))
323
!-----------------------------------------------------------------------------------
325
! Implement terms as source or absorbtion
327
call get_option(trim(option_path)//&
328
'time_discretisation/source_term_implementation/'//&
329
trim(src_abs_terms(term)%name), implementation)
330
select case(implementation)
332
call addto(src, src_abs_terms(term))
334
call allocate(src_to_abs, fields(1)%mesh, name='SourceToAbsorbtion')
335
call set(src_to_abs, fields(1))
336
where (src_to_abs%val >= fields_min)
337
src_to_abs%val=1./src_to_abs%val
339
src_to_abs%val=1./fields_min
341
call scale(src_abs_terms(term), src_to_abs)
342
call addto(abs, src_abs_terms(term), -1.0)
343
call deallocate(src_to_abs)
345
! developer error... out of sync options input and code
346
FLAbort("Unknown implementation type for k-epsilon source terms")
349
!-----------------------------------------------------------------------------------
351
! This allows user-specified source and absorption terms, so that an MMS test can be
353
debug => extract_scalar_field(state, &
354
trim(field_names(i))//"PrescribedSource", stat)
356
call addto(src, debug)
358
!-----------------------------------------------------------------------------------
362
call deallocate(src_abs_terms(term))
364
call deallocate(fields(1))
365
call deallocate(fields(2))
369
call deallocate(dummydensity)
370
deallocate(dummydensity)
372
end subroutine keps_calculate_rhs
374
!------------------------------------------------------------------------------!
376
subroutine assemble_rhs_ele(src_abs_terms, k, eps, scalar_eddy_visc, u, density, &
377
buoyancy_density, have_buoyancy_turbulence, g, g_magnitude, &
378
X, c_eps_1, c_eps_2, sigma_p, f_1, f_2, ele, field_id)
380
type(scalar_field), dimension(3), intent(inout) :: src_abs_terms
381
type(scalar_field), intent(in) :: k, eps, scalar_eddy_visc, f_1, f_2
382
type(vector_field), intent(in) :: X, u, g
383
type(scalar_field), intent(in) :: density, buoyancy_density
384
real, intent(in) :: g_magnitude, c_eps_1, c_eps_2, sigma_p
385
logical, intent(in) :: have_buoyancy_turbulence
386
integer, intent(in) :: ele, field_id
388
real, dimension(ele_loc(k, ele), ele_ngi(k, ele), x%dim) :: dshape
389
real, dimension(ele_ngi(k, ele)) :: detwei, rhs, scalar_eddy_visc_ele, k_ele, eps_ele
390
real, dimension(3, ele_loc(k, ele)) :: rhs_addto
391
integer, dimension(ele_loc(k, ele)) :: nodes
392
real, dimension(ele_loc(k, ele), ele_loc(k, ele)) :: invmass
393
real, dimension(u%dim, u%dim, ele_ngi(k, ele)) :: reynolds_stress, grad_u
394
type(element_type), pointer :: shape
395
integer :: term, ngi, dim, gi, i
397
! For buoyancy turbulence stuff
398
real, dimension(u%dim, ele_ngi(u, ele)) :: vector, u_quad, g_quad
400
real, dimension(ele_ngi(u, ele)) :: scalar, c_eps_3
401
type(element_type), pointer :: shape_density
402
real, dimension(:, :, :), allocatable :: dshape_density
404
shape => ele_shape(k, ele)
405
nodes = ele_nodes(k, ele)
407
call transform_to_physical( X, ele, shape, dshape=dshape, detwei=detwei )
409
! get bounded values of k and epsilon for source terms
410
! this doesn't change the field values of k and epsilon
411
k_ele = ele_val_at_quad(k,ele)
412
eps_ele = ele_val_at_quad(eps, ele)
413
ngi = ele_ngi(u, ele)
415
k_ele(gi) = max(k_ele(gi), fields_min)
416
eps_ele(gi) = max(eps_ele(gi), fields_min)
419
! Compute Reynolds stress
420
grad_u = ele_grad_at_quad(u, ele, dshape)
421
scalar_eddy_visc_ele = ele_val_at_quad(scalar_eddy_visc, ele)
424
reynolds_stress(:,:,gi) = scalar_eddy_visc_ele(gi)*(grad_u(:,:,gi) + transpose(grad_u(:,:,gi)))
427
reynolds_stress(i,i,:) = reynolds_stress(i,i,:) - (2./3.)*k_ele*ele_val_at_quad(density, ele)
431
rhs = tensor_inner_product(reynolds_stress, grad_u)
432
if (field_id==2) then
433
rhs = rhs*c_eps_1*ele_val_at_quad(f_1,ele)*eps_ele/k_ele
435
rhs_addto(1,:) = shape_rhs(shape, detwei*rhs)
438
rhs = -1.0*eps_ele*ele_val_at_quad(density, ele)
439
if (field_id==2) then
440
rhs = rhs*c_eps_2*ele_val_at_quad(f_2,ele)*eps_ele/k_ele
442
rhs_addto(2,:) = shape_rhs(shape, detwei*rhs)
445
! Calculate buoyancy turbulence term and add to addto array
446
if(have_buoyancy_turbulence) then
448
! calculate scalar and vector components of the source term
449
allocate(dshape_density(ele_loc(buoyancy_density, ele), ele_ngi(buoyancy_density, ele), X%dim))
450
if(.not.(buoyancy_density%mesh == k%mesh)) then
451
shape_density => ele_shape(buoyancy_density, ele)
452
call transform_to_physical( X, ele, shape_density, dshape=dshape_density )
454
dshape_density = dshape
457
scalar = -1.0*g_magnitude*ele_val_at_quad(scalar_eddy_visc, ele)/(sigma_p*ele_val_at_quad(density,ele))
458
vector = ele_val_at_quad(g, ele)*ele_grad_at_quad(buoyancy_density, ele, dshape_density)
460
! multiply vector component by scalar and sum across dimensions - note that the
461
! vector part has been multiplied by the gravitational direction so that it is
462
! zero everywhere apart from in this direction.
464
scalar(gi) = sum(scalar(gi) * vector(:, gi))
467
if (field_id == 2) then
468
! calculate c_eps_3 = tanh(v/u)
469
g_quad = ele_val_at_quad(g, ele)
470
u_quad = ele_val_at_quad(u, ele)
472
! get components of velocity in direction of gravity and in other directions
473
u_z = dot_product(g_quad(:, gi), u_quad(:, gi))
474
u_xy = (norm2(u_quad(:, gi))**2.0 - u_z**2.0)**0.5
475
if (u_xy > fields_min) then
476
c_eps_3(gi) = tanh(u_z/u_xy)
481
scalar = scalar*c_eps_1*ele_val_at_quad(f_1,ele)*c_eps_3*eps_ele/k_ele
484
! multiply by determinate weights, integrate and assign to rhs
485
rhs_addto(3,:) = shape_rhs(shape, scalar * detwei)
487
deallocate(dshape_density)
490
! No buoyancy term, so set this part of the array to zero.
494
! In the DG case we apply the inverse mass locally.
495
if(continuity(k)<0) then
496
invmass = inverse(shape_shape(shape, shape, detwei))
498
rhs_addto(term,:) = matmul(rhs_addto(term,:), invmass)
503
call addto(src_abs_terms(term), nodes, rhs_addto(term,:))
506
end subroutine assemble_rhs_ele
508
!------------------------------------------------------------------------------!
509
! calculate inner product for 2xN matrices dim,dim,N !
510
!------------------------------------------------------------------------------!
511
function tensor_inner_product(A, B)
513
real, dimension(:,:,:), intent(in) :: A, B
515
real, dimension(size(A,1), size(A,2), size(A,3)) :: C
516
real, dimension(size(A,3)) :: tensor_inner_product
521
tensor_inner_product(i) = sum(C(:,:,i))
524
end function tensor_inner_product
527
! eddyvisc calculates the lengthscale and the eddy viscosity
399
528
! Eddy viscosity is added to the background viscosity.
402
subroutine keps_eddyvisc(state)
404
type(state_type), intent(inout) :: state
405
type(tensor_field), pointer :: eddy_visc, viscosity, bg_visc
406
type(vector_field), pointer :: positions
407
type(scalar_field), pointer :: kk, eps, EV, ll, lumped_mass
408
type(scalar_field) :: ev_rhs
409
type(element_type), pointer :: shape_ev
411
integer, pointer, dimension(:) :: nodes_ev
412
real, allocatable, dimension(:) :: detwei, rhs_addto
414
ewrite(1,*) "In keps_eddyvisc"
415
kk => extract_scalar_field(state, "TurbulentKineticEnergy")
416
eps => extract_scalar_field(state, "TurbulentDissipation")
417
positions => extract_vector_field(state, "Coordinate")
418
eddy_visc => extract_tensor_field(state, "EddyViscosity")
419
viscosity => extract_tensor_field(state, "Viscosity")
420
bg_visc => extract_tensor_field(state, "BackgroundViscosity")
421
EV => extract_scalar_field(state, "ScalarEddyViscosity")
422
ll => extract_scalar_field(state, "LengthScale")
428
call allocate(ev_rhs, EV%mesh, name="EVRHS")
431
! Initialise viscosity to background value
432
call set(viscosity, bg_visc)
434
!Clip fields: can't allow negative/zero epsilon or k
435
do i = 1, node_count(EV)
436
call set(kk, i, max(node_val(kk,i), fields_min))
437
call set(eps, i, max(node_val(eps,i), fields_min))
438
! Limit lengthscale to prevent instablilities.
439
call set(ll, i, min(node_val(kk,i)**1.5 / node_val(eps,i), lmax))
442
! Calculate scalar eddy viscosity by integration over element
443
do ele = 1, ele_count(EV)
444
nodes_ev => ele_nodes(EV, ele)
445
shape_ev => ele_shape(EV, ele)
446
allocate(detwei (ele_ngi(EV, ele)))
447
allocate(rhs_addto (ele_loc(EV, ele)))
448
call transform_to_physical(positions, ele, detwei=detwei)
449
rhs_addto = shape_rhs(shape_ev, detwei*C_mu*ele_val_at_quad(kk,ele)**0.5*ele_val_at_quad(ll,ele))
450
call addto(ev_rhs, nodes_ev, rhs_addto)
451
deallocate(detwei, rhs_addto)
454
lumped_mass => get_lumped_mass(state, EV%mesh)
457
do i = 1, node_count(EV)
458
call set(EV, i, node_val(ev_rhs,i)/node_val(lumped_mass,i))
461
call deallocate(ev_rhs)
463
ewrite(2,*) "Setting k-epsilon eddy-viscosity tensor"
466
! eddy viscosity tensor is isotropic
467
do i = 1, eddy_visc%dim(1)
468
call set(eddy_visc, i, i, EV)
471
! Add turbulence model contribution to viscosity field
472
call addto(viscosity, eddy_visc)
474
ewrite_minmax(eddy_visc)
475
ewrite_minmax(viscosity)
530
subroutine keps_eddyvisc(state, advdif)
532
type(state_type), intent(inout) :: state
533
logical, intent(in) :: advdif
535
type(tensor_field), pointer :: eddy_visc, viscosity, bg_visc
536
type(vector_field), pointer :: x, u
537
type(scalar_field) :: kk, eps
538
type(scalar_field), pointer :: scalar_eddy_visc, ll, f_mu, density, dummydensity
539
type(scalar_field) :: ev_rhs
540
integer :: i, j, ele, stat
542
! Options grabbed from the options tree
544
character(len=OPTION_PATH_LEN) :: option_path
545
logical :: lump_mass, have_visc = .true.
546
character(len=FIELD_NAME_LEN) :: equation_type
548
option_path = trim(state%option_path)//'/subgridscale_parameterisations/k-epsilon/'
550
if (.not. have_option(trim(option_path))) then
554
ewrite(1,*) "In keps_eddyvisc"
557
call get_option(trim(option_path)//'/C_mu', c_mu, default = 0.09)
560
call time_averaged_value(state, kk, "TurbulentKineticEnergy", advdif, option_path)
561
call time_averaged_value(state, eps, "TurbulentDissipation", advdif, option_path)
562
x => extract_vector_field(state, "Coordinate")
563
u => extract_vector_field(state, "NonlinearVelocity")
564
eddy_visc => extract_tensor_field(state, "EddyViscosity")
565
f_mu => extract_scalar_field(state, "f_mu")
566
bg_visc => extract_tensor_field(state, "BackgroundViscosity")
567
scalar_eddy_visc => extract_scalar_field(state, "ScalarEddyViscosity")
568
ll => extract_scalar_field(state, "LengthScale")
569
viscosity => extract_tensor_field(state, "Viscosity", stat)
576
ewrite_minmax(scalar_eddy_visc)
578
allocate(dummydensity)
579
call allocate(dummydensity, X%mesh, "DummyDensity", field_type=FIELD_TYPE_CONSTANT)
580
call set(dummydensity, 1.0)
581
dummydensity%option_path = ""
583
! Depending on the equation type, extract the density or set it to some dummy field allocated above
584
call get_option(trim(state%option_path)//&
585
"/vector_field::Velocity/prognostic/equation[0]/name", equation_type)
586
select case(equation_type)
587
case("LinearMomentum")
588
density=>extract_scalar_field(state, "Density")
590
density=>dummydensity
592
density=>dummydensity
594
! developer error... out of sync options input and code
595
FLAbort("Unknown equation type for velocity")
598
call allocate(ev_rhs, scalar_eddy_visc%mesh, name="EVRHS")
601
! Initialise viscosity to background value
603
call set(viscosity, bg_visc)
606
! Compute the length scale diagnostic field here.
607
do i = 1, node_count(scalar_eddy_visc)
608
call set(ll, i, max(node_val(kk,i), fields_min)**1.5 / max(node_val(eps,i), fields_min))
611
! Calculate scalar eddy viscosity by integration over element
612
do ele = 1, ele_count(scalar_eddy_visc)
613
call keps_eddyvisc_ele(ele, X, kk, eps, scalar_eddy_visc, f_mu, density, ev_rhs)
616
! For non-DG we apply inverse mass globally
617
if(continuity(scalar_eddy_visc)>=0) then
618
lump_mass = have_option(trim(option_path)//'mass_terms/lump_mass')
619
call solve_cg_inv_mass(state, ev_rhs, lump_mass, option_path)
622
! Allow for prescribed eddy-viscosity
623
if (.not. have_option(trim(option_path)//'/scalar_field::ScalarEddyViscosity/prescribed')) then
624
call set(scalar_eddy_visc, ev_rhs)
627
call deallocate(ev_rhs)
631
call deallocate(dummydensity)
632
deallocate(dummydensity)
634
ewrite(2,*) "Setting k-epsilon eddy-viscosity tensor"
637
! this is skipped if zero_eddy_viscosity is set - this is the easiest way to
638
! disable feedback from the k-epsilon model back into the rest of the model
639
if (.not. have_option(trim(option_path)//'debugging_options/zero_reynolds_stress_tensor')) then
640
do i = 1, eddy_visc%dim(1)
641
do j = 1, eddy_visc%dim(1)
642
call set(eddy_visc, i, j, scalar_eddy_visc)
647
! Add turbulence model contribution to viscosity field
649
call addto(viscosity, eddy_visc)
652
ewrite_minmax(eddy_visc)
653
ewrite_minmax(viscosity)
654
ewrite_minmax(scalar_eddy_visc)
659
subroutine keps_eddyvisc_ele(ele, X, kk, eps, scalar_eddy_visc, f_mu, density, ev_rhs)
661
type(vector_field), intent(in) :: x
662
type(scalar_field), intent(in) :: kk, eps, scalar_eddy_visc, f_mu, density
663
type(scalar_field), intent(inout):: ev_rhs
664
integer, intent(in) :: ele
666
type(element_type), pointer :: shape_ev
667
integer, pointer, dimension(:) :: nodes_ev
668
real, dimension(ele_ngi(scalar_eddy_visc, ele)) :: detwei
669
real, dimension(ele_loc(scalar_eddy_visc, ele)) :: rhs_addto
670
real, dimension(ele_loc(scalar_eddy_visc, ele), ele_loc(scalar_eddy_visc, ele)) :: invmass
671
real, dimension(ele_ngi(kk, ele)) :: kk_at_quad, eps_at_quad
674
nodes_ev => ele_nodes(scalar_eddy_visc, ele)
675
shape_ev => ele_shape(scalar_eddy_visc, ele)
678
call transform_to_physical(X, ele, detwei=detwei)
680
! Get the k and epsilon values at the Gauss points
681
kk_at_quad = ele_val_at_quad(kk,ele)
682
eps_at_quad = ele_val_at_quad(eps,ele)
684
! Clip the field values at the Gauss points.
685
! Note 1: This isn't a permanent change directly to the field itself,
686
! only to the values used in the computation of the eddy viscosity.
687
! Note 2: Can't allow negative/zero epsilon or k.
688
! Note 3: Here we assume all fields have the same number of
689
! Gauss points per element.
690
where (kk_at_quad < fields_min)
691
kk_at_quad = fields_min
693
where (eps_at_quad < fields_min)
694
eps_at_quad = fields_min
697
! Compute the eddy viscosity
698
rhs_addto = shape_rhs(shape_ev, detwei*C_mu*ele_val_at_quad(density,ele)*&
699
ele_val_at_quad(f_mu,ele)*(kk_at_quad**2.0)/eps_at_quad)
701
! In the DG case we will apply the inverse mass locally.
702
if(continuity(scalar_eddy_visc)<0) then
703
invmass = inverse(shape_shape(shape_ev, shape_ev, detwei))
704
rhs_addto = matmul(rhs_addto, invmass)
707
! Add the element's contribution to the nodes of ev_rhs
708
call addto(ev_rhs, nodes_ev, rhs_addto)
710
end subroutine keps_eddyvisc_ele
481
712
end subroutine keps_eddyvisc
714
!---------------------------------------------------------------------------------
716
subroutine keps_diffusion(state)
718
! calculates k and epsilon field diffusivities
719
type(state_type), intent(inout) :: state
721
type(tensor_field), pointer :: diff, bg_visc, eddy_visc
722
real :: sigma_k, sigma_eps
724
character(len=OPTION_PATH_LEN) :: option_path
726
ewrite(1,*) 'in keps_diffusion'
728
option_path = trim(state%option_path)//'/subgridscale_parameterisations/k-epsilon/'
730
eddy_visc => extract_tensor_field(state, "EddyViscosity")
731
bg_visc => extract_tensor_field(state, "BackgroundViscosity")
732
call get_option(trim(option_path)//'/sigma_k', sigma_k, default = 1.0)
733
call get_option(trim(option_path)//'/sigma_eps', sigma_eps, default = 1.3)
736
diff => extract_tensor_field(state, "TurbulentKineticEnergyDiffusivity")
738
do i = 1, node_count(diff)
739
do j = 1, diff%dim(1)
740
call addto(diff, j, j, i, node_val(bg_visc, j, j, i))
741
call addto(diff, j, j, i, node_val(eddy_visc, j, j, i) / sigma_k)
744
diff => extract_tensor_field(state, "TurbulentDissipationDiffusivity")
746
do i = 1, node_count(diff)
747
do j = 1, diff%dim(1)
748
call addto(diff, j, j, i, node_val(bg_visc, j, j, i))
749
call addto(diff, j, j, i, node_val(eddy_visc, j, j, i) / sigma_eps)
753
end subroutine keps_diffusion
755
!---------------------------------------------------------------------------------
757
subroutine keps_tracer_diffusion(state)
759
! calculates scalar field diffusivity based upon eddy viscosity and background
761
type(state_type), intent(inout) :: state
763
type(tensor_field), pointer :: t_field
764
integer :: i_field, i, stat
765
real :: sigma_p, local_background_diffusivity
766
type(scalar_field) :: local_background_diffusivity_field
767
type(scalar_field), pointer :: scalar_eddy_viscosity, s_field
768
type(tensor_field), pointer :: global_background_diffusivity
769
type(tensor_field) :: background_diffusivity
771
ewrite(1,*) 'In keps_tracer_diffusion'
773
do i_field = 1, scalar_field_count(state)
774
s_field => extract_scalar_field(state, i_field)
776
if (have_option(trim(s_field%option_path)//&
777
'/prognostic/subgridscale_parameterisation::k-epsilon')) then
779
ewrite(1,*) 'Calculating turbulent diffusivity for field: ', s_field%name
782
if (.not.(have_option(trim(state%option_path)//'/subgridscale_parameterisations/k-epsilon')))&
784
FLExit('you must have /subgridscale_parameterisations/k-epsilon to be able to calculate diffusivity based upon the k-epsilon model')
787
t_field => extract_tensor_field(state, trim(s_field%name)//'Diffusivity', stat=stat)
789
FLExit('you must have a Diffusivity field to be able to calculate diffusivity based upon the k-epsilon model')
790
else if (.not. have_option(trim(t_field%option_path)//"/diagnostic/algorithm::Internal")) then
791
FLExit('you must have a diagnostic Diffusivity field with algorithm::Internal to be able to calculate diffusivity based upon the k-epsilon model')
795
call get_option(trim(state%option_path)//'/subgridscale_parameterisations/k-epsilon/sigma_p', sigma_p)
797
! allocate and zero required fields
798
call allocate(background_diffusivity, t_field%mesh, name="background_diffusivity")
799
call zero(background_diffusivity)
800
call allocate(local_background_diffusivity_field, t_field%mesh, &
801
name="local_background_diffusivity_field")
802
call zero(local_background_diffusivity_field)
804
! set background_diffusivity (local takes precendence over global)
805
call get_option(trim(s_field%option_path)//&
806
'/prognostic/subgridscale_parameterisation::k-epsilon/background_diffusivity', &
807
local_background_diffusivity, stat=stat)
809
! set local isotropic background diffusivity
810
call addto(local_background_diffusivity_field, local_background_diffusivity)
811
do i = 1, background_diffusivity%dim(1)
812
call set(background_diffusivity, i, i, local_background_diffusivity_field)
815
global_background_diffusivity => extract_tensor_field(state, 'BackgroundDiffusivity', stat=stat)
817
call set(background_diffusivity, global_background_diffusivity)
822
scalar_eddy_viscosity => extract_scalar_field(state, 'ScalarEddyViscosity', stat)
825
call addto(t_field, background_diffusivity)
826
do i = 1, t_field%dim(1)
827
call addto(t_field, i, i, scalar_eddy_viscosity, 1.0/sigma_p)
830
call deallocate(background_diffusivity)
831
call deallocate(local_background_diffusivity_field)
836
end subroutine keps_tracer_diffusion
483
838
!--------------------------------------------------------------------------------!
484
839
! This gets and applies locally defined boundary conditions (wall functions) !
485
840
!--------------------------------------------------------------------------------!
487
842
subroutine keps_bcs(state)
489
type(state_type), intent(in) :: state
490
type(scalar_field), pointer :: field1, field2 ! k or epsilon
491
type(scalar_field), pointer :: surface_field, EV
492
type(vector_field), pointer :: positions, u
493
type(tensor_field), pointer :: bg_visc
494
type(scalar_field) :: rhs_field, surface_values
495
type(mesh_type), pointer :: surface_mesh
496
integer :: i, j, ele, sele, index, nbcs
497
integer, dimension(:), pointer :: surface_elements, surface_node_list
498
character(len=FIELD_NAME_LEN) :: bc_type, bc_name, wall_fns
499
character(len=OPTION_PATH_LEN) :: bc_path, bc_path_i
500
integer, allocatable, dimension(:) :: nodes_bdy
501
real, dimension(:), allocatable :: rhs
504
ewrite(2,*) "In keps_bcs"
506
positions => extract_vector_field(state, "Coordinate")
507
u => extract_vector_field(state, "Velocity")
508
EV => extract_scalar_field(state, "ScalarEddyViscosity")
509
bg_visc => extract_tensor_field(state, "BackgroundViscosity")
511
! THIS IS NOT AVAILABLE BEFORE KEPS_INIT HAS BEEN CALLED! Grab it here for initialising BCs.
512
call get_option(trim(state%option_path)//"/subgridscale_parameterisations/k-epsilon/C_mu", cmu, default = 0.09)
513
ewrite(2,*) "cmu: ", cmu
515
field_loop: do index=1,2
844
type(state_type), intent(in) :: state
845
type(scalar_field), pointer :: field1, field2 ! k or epsilon
846
type(scalar_field), pointer :: f_1, f_2, f_mu
847
type(scalar_field), pointer :: surface_field, scalar_eddy_visc
848
type(scalar_field), pointer :: density, dummydensity
849
type(vector_field), pointer :: X, u
850
type(tensor_field), pointer :: bg_visc
851
type(scalar_field) :: rhs_field, surface_values, masslump
852
type(mesh_type), pointer :: surface_mesh
853
integer :: i, j, sele, index, nbcs, stat
854
integer, dimension(:), pointer :: surface_elements, surface_node_list
855
character(len=FIELD_NAME_LEN) :: bc_type, bc_name, wall_fns
856
character(len=OPTION_PATH_LEN) :: bc_path, bc_path_i, option_path
858
character(len=FIELD_NAME_LEN) :: equation_type
860
option_path = trim(state%option_path)//'/subgridscale_parameterisations/k-epsilon/'
862
ewrite(2,*) "In keps_bcs"
864
X => extract_vector_field(state, "Coordinate")
865
u => extract_vector_field(state, "Velocity")
866
scalar_eddy_visc => extract_scalar_field(state, "ScalarEddyViscosity")
867
bg_visc => extract_tensor_field(state, "BackgroundViscosity")
868
f_1 => extract_scalar_field(state, "f_1")
869
f_2 => extract_scalar_field(state, "f_2")
870
f_mu => extract_scalar_field(state, "f_mu")
872
allocate(dummydensity)
873
call allocate(dummydensity, X%mesh, "DummyDensity", field_type=FIELD_TYPE_CONSTANT)
874
call set(dummydensity, 1.0)
875
dummydensity%option_path = ""
877
! Depending on the equation type, extract the density or set it to some dummy field allocated above
878
call get_option(trim(u%option_path)//"/prognostic/equation[0]/name", equation_type)
879
select case(equation_type)
880
case("LinearMomentum")
881
density=>extract_scalar_field(state, "Density")
883
density=>dummydensity
885
density=>dummydensity
887
! developer error... out of sync options input and code
888
FLAbort("Unknown equation type for velocity")
891
call get_option(trim(option_path)//"C_mu", c_mu, default = 0.09)
893
field_loop: do index=1,2
518
896
field1 => extract_scalar_field(state, "TurbulentKineticEnergy")
521
899
field1 => extract_scalar_field(state, "TurbulentDissipation")
522
900
field2 => extract_scalar_field(state, "TurbulentKineticEnergy")
525
bc_path=trim(field1%option_path)//'/prognostic/boundary_conditions'
526
nbcs=option_count(trim(bc_path))
528
! Loop over boundary conditions for field1
529
boundary_conditions: do i=0, nbcs-1
903
bc_path=trim(field1%option_path)//'/prognostic/boundary_conditions'
904
nbcs=option_count(trim(bc_path))
906
! Loop over boundary conditions for field1
907
boundary_conditions: do i=0, nbcs-1
531
909
bc_path_i=trim(bc_path)//"["//int2str(i)//"]"
533
911
! Get name and type of boundary condition
534
912
call get_option(trim(bc_path_i)//"/name", bc_name)
535
913
call get_option(trim(bc_path_i)//"/type[0]/name", bc_type)
537
if (trim(bc_type) .eq. "k_epsilon") then
538
! Get bc by name. Get type just to make sure it's now dirichlet
539
call get_boundary_condition(field1, name=bc_name, type=bc_type, surface_node_list=surface_node_list, &
540
surface_element_list=surface_elements, surface_mesh=surface_mesh)
542
! Do we have high- or low-Reynolds-number wall functions?
543
call get_option(trim(bc_path_i)//"/type::k_epsilon/", wall_fns)
545
ewrite(2,*) "Calculating field BC: ",trim(field1%name),' ',trim(bc_name),' ',trim(bc_type),' ',trim(wall_fns)
547
! Get surface field already created in bcs_from_options
548
surface_field => extract_surface_field(field1, bc_name=bc_name, name="value")
549
call zero(surface_field)
551
call allocate(surface_values, surface_mesh, name="surfacevalues")
552
call allocate(rhs_field, field1%mesh, name="rhs")
553
call zero(surface_values); call zero(rhs_field)
555
do j = 1, ele_count(surface_mesh)
556
sele = surface_elements(j)
557
ele = face_ele(rhs_field, sele)
559
allocate(rhs(face_loc(rhs_field, sele)))
560
allocate(nodes_bdy(face_loc(rhs_field, sele)))
561
nodes_bdy = face_global_nodes(rhs_field, sele)
563
! Calculate wall function
564
call keps_wall_function(field1,field2,positions,u,bg_visc,EV,ele,sele,index,wall_fns,cmu,rhs)
566
! Add element contribution to rhs field
567
call addto(rhs_field, nodes_bdy, rhs)
569
deallocate(rhs); deallocate(nodes_bdy)
572
! Put values onto surface mesh
573
call remap_field_to_surface(rhs_field, surface_values, surface_elements)
574
ewrite_minmax(rhs_field)
575
do j = 1, size(surface_node_list)
576
call set(surface_field, j, node_val(surface_values, j))
579
call deallocate(surface_values); call deallocate(rhs_field)
914
! Do we have high- or low-Reynolds-number wall functions?
915
call get_option(trim(bc_path_i)//"/type::k_epsilon/", wall_fns, stat=stat)
917
if (trim(bc_type)=="k_epsilon" .and. wall_fns=="low_Re") then
919
! lowRe BC's are just zero Dirichlet or Neumann - damping functions get calculated in
923
else if (trim(bc_type)=="k_epsilon" .and. wall_fns=="high_Re") then
925
! Get bc by name. Get type just to make sure it's now dirichlet
926
call get_boundary_condition(field1, name=bc_name, type=bc_type, surface_node_list=surface_node_list, &
927
surface_element_list=surface_elements, surface_mesh=surface_mesh)
929
ewrite(2,*) "Calculating highRe k-epsilon BC: ",trim(field1%name),' ',trim(bc_name),' ',trim(bc_type),' ',trim(wall_fns)
931
! Get surface field already created in bcs_from_options
932
surface_field => extract_surface_field(field1, bc_name=bc_name, name="value")
933
call zero(surface_field)
935
call allocate(surface_values, surface_mesh, name="surfacevalues")
936
call allocate(rhs_field, field1%mesh, name="rhs")
937
call zero(surface_values); call zero(rhs_field)
939
if(continuity(field1)>=0) then
940
call allocate(masslump, field1%mesh, 'Masslump')
944
! Calculate high-Re wall function
945
do sele = 1, ele_count(surface_mesh)
946
call keps_wall_function(field1, X, u, masslump, scalar_eddy_visc, density, sele, &
947
index, c_mu, rhs_field)
950
! In the continuous case we globally apply inverse mass
951
if(continuity(field1)>=0) then
952
where (masslump%val/=0.0)
953
masslump%val=1./masslump%val
955
call scale(rhs_field, masslump)
956
call deallocate(masslump)
959
! Put values onto surface mesh
960
call remap_field_to_surface(rhs_field, surface_values, surface_elements)
961
ewrite_minmax(rhs_field)
962
do j = 1, size(surface_node_list)
963
call set(surface_field, j, node_val(surface_values, j))
966
call deallocate(surface_values); call deallocate(rhs_field)
582
end do boundary_conditions
969
end do boundary_conditions
972
call deallocate(dummydensity)
973
deallocate(dummydensity)
585
975
end subroutine keps_bcs
587
!------------------------------------------------------------------------------------
589
subroutine keps_cleanup()
591
ewrite(1,*) "In keps_cleanup"
593
call deallocate(tke_old)
594
call deallocate(tke_src_old)
596
end subroutine keps_cleanup
598
!------------------------------------------------------------------------------------!
599
! Called after an adapt to reset the fields and arrays within the module !
600
!------------------------------------------------------------------------------------!
602
subroutine keps_adapt_mesh(state)
604
type(state_type), intent(inout) :: state
606
ewrite(1,*) "In keps_adapt_mesh"
608
call keps_allocate_fields(state)
609
call keps_eddyvisc(state)
612
call keps_eddyvisc(state)
614
end subroutine keps_adapt_mesh
616
!---------------------------------------------------------------------------------
618
subroutine keps_check_options(state)
620
type(state_type), intent(in) :: state
621
type(vector_field), pointer :: u
622
character(len=OPTION_PATH_LEN) :: option_path
623
character(len=FIELD_NAME_LEN) :: kmsh, emsh, vmsh
626
ewrite(1,*) "In keps_check_options"
627
option_path = trim(state%option_path)//"/subgridscale_parameterisations/k-epsilon"
628
u => extract_vector_field(state, "Velocity")
630
! one dimensional problems not supported
631
call get_option("/geometry/dimension/", dimension)
632
if (dimension .eq. 1 .and. have_option(trim(option_path))) then
977
!--------------------------------------------------------------------------------!
978
! Only used if bc type == k_epsilon for field and high_Re. !
979
!--------------------------------------------------------------------------------!
980
subroutine keps_wall_function(field1,X,U,masslump,scalar_eddy_visc,density,sele,index,c_mu,rhs_field)
982
type(scalar_field), pointer, intent(in) :: field1, scalar_eddy_visc, density
983
type(vector_field), pointer, intent(in) :: X, U
984
type(scalar_field), intent(inout) :: masslump, rhs_field
985
integer, intent(in) :: sele, index
986
type(element_type), pointer :: shape, f_shape
987
type(element_type) :: augmented_shape
988
integer :: i, j, gi, ele, dim
989
integer, dimension(face_loc(rhs_field, sele)) :: nodes_bdy
990
real :: kappa, h, c_mu
991
real, dimension(1,1) :: hb
992
real, dimension(face_ngi(field1,sele)) :: detwei, tau_wall_at_quad, visc_at_quad, density_at_quad
993
real, dimension(X%dim,1) :: n
994
real, dimension(X%dim,X%dim) :: G
995
real, dimension(face_loc(rhs_field, sele)) :: rhs
996
real, dimension(X%dim,face_ngi(field1,sele)) :: normal, normal_shear_at_quad
997
real, dimension(X%dim, X%dim, face_ngi(field1, sele)):: f_invJ, grad_U_at_quad, shear_at_quad
998
real, dimension(face_loc(field1, sele), face_loc(field1, sele)) :: fmass
999
real, dimension(X%dim, X%dim, ele_ngi(field1, face_ele(field1, sele))) :: invJ
1000
real, dimension(ele_loc(X, face_ele(field1, sele)), face_ngi(field1, sele), X%dim) :: ele_dshape_at_face_quad
1002
nodes_bdy = face_global_nodes(rhs_field, sele) ! nodes in rhs_field
1003
ele = face_ele(field1, sele) ! ele number for volume mesh
1004
dim = mesh_dim(X) ! field dimension
1006
! get shape functions
1007
f_shape => face_shape(field1, sele)
1008
shape => ele_shape(field1, ele)
1010
! generate shape functions that include quadrature points on the face required
1011
! check that the shape does not already have these first
1012
assert(shape%degree == 1)
1013
if(associated(shape%dn_s)) then
1014
augmented_shape = shape
1015
call incref(augmented_shape)
1017
augmented_shape = make_element_shape(shape%loc, shape%dim, shape%degree, shape%quadrature, quad_s=f_shape%quadrature)
1020
! assumes that the jacobian is the same for all quadrature points
1021
! this is not valid for spheres!
1022
call compute_inverse_jacobian(ele_val(X, ele), ele_shape(X, ele), invj = invJ)
1023
assert(ele_numbering_family(shape) == FAMILY_SIMPLEX)
1024
f_invJ = spread(invJ(:, :, 1), 3, size(f_invJ, 3))
1026
call transform_facet_to_physical(X, sele, detwei_f = detwei, normal = normal)
1028
! Evaluate the volume element shape function derivatives at the surface
1029
! element quadrature points
1030
ele_dshape_at_face_quad = eval_volume_dshape_at_face_quad(augmented_shape, &
1031
& local_face_number(X, sele), f_invJ)
1033
! Calculate grad of U at the surface element quadrature points
1036
grad_U_at_quad(i, j, :) = matmul(ele_val(U, j, ele), ele_dshape_at_face_quad(:,:,i))
1040
density_at_quad = face_val_at_quad(density, sele)
1041
visc_at_quad = face_val_at_quad(scalar_eddy_visc, sele)
1042
do gi = 1, face_ngi(X, sele)
1043
! Multiply by visosity
1044
shear_at_quad(:,:,gi) = grad_U_at_quad(:,:,gi)*visc_at_quad(gi)
1045
! Multiply by surface normal (dim,sgi) to obtain shear in direction normal
1046
! to surface - transpose (because fluidity stores data in row-major order??)
1047
normal_shear_at_quad(:,gi) = matmul(shear_at_quad(:,:,gi), normal(:,gi))
1048
! Wall shear stress tau_wall/rho = nu_T*du/dn.
1049
! Get streamwise velocity gradient by taking sqrt(grad_n.grad_n).
1050
! Eddy viscosity is dynamic not kinematic so divide by density.
1051
tau_wall_at_quad(gi) = norm2(normal_shear_at_quad(:,gi))
1055
! k = tau_wall/rho/c_mu**0.5
1056
rhs = shape_rhs(f_shape, tau_wall_at_quad/density_at_quad/c_mu**0.5*detwei)
1057
else if (index==2) then
1058
! calculate wall-normal element size
1059
G = matmul(transpose(invJ(:,:,1)), invJ(:,:,1))
1060
n(:,1) = normal(:,1)
1061
hb = 1. / sqrt( matmul(matmul(transpose(n), G), n) )
1063
! Von Karman's constant
1065
! epsilon = (tau_wall/rho)**1.5/kappa/h
1066
rhs = shape_rhs(f_shape, (tau_wall_at_quad/density_at_quad)**1.5/kappa/h*detwei)
1069
fmass = shape_shape(f_shape, f_shape, detwei)
1070
! In the CG case we need to calculate a global lumped mass
1071
if(continuity(field1)>=0) then
1072
call addto(masslump, nodes_bdy, sum(fmass,1))
1073
else ! In the DG case we will apply the inverse mass locally.
1075
rhs = matmul(inverse(fmass), rhs)
1080
call addto(rhs_field, nodes_bdy, rhs)
1082
call deallocate(augmented_shape)
1084
end subroutine keps_wall_function
1086
!---------------------------------------------------------------------------------
1088
subroutine time_averaged_value(state, A, field_name, advdif, option_path)
1090
type(state_type), intent(in) :: state
1091
type(scalar_field), intent(inout) :: A
1092
character(len=*), intent(in) :: field_name
1093
logical, intent(in) :: advdif ! advdif or mom - whether to use old or iterated values
1094
character(len=OPTION_PATH_LEN), intent(in) :: option_path
1097
type(scalar_field), pointer :: old, iterated
1099
call get_option(trim(option_path)//'time_discretisation/theta', theta)
1101
old => extract_scalar_field(state, "Old"//trim(field_name))
1103
iterated => extract_scalar_field(state, "Iterated"//trim(field_name))
1105
iterated => extract_scalar_field(state, trim(field_name))
1108
call allocate(A, old%mesh, name="Nonlinear"//trim(field_name))
1110
call addto(A, old, 1.0-theta)
1111
call addto(A, iterated, theta)
1114
ewrite_minmax(iterated)
1117
end subroutine time_averaged_value
1119
!---------------------------------------------------------------------------------
1121
subroutine solve_cg_inv_mass(state, A, lump, option_path)
1123
type(state_type), intent(inout) :: state
1124
type(scalar_field), intent(inout) :: A
1125
logical, intent(in) :: lump
1126
character(len=OPTION_PATH_LEN), intent(in) :: option_path
1128
type(scalar_field), pointer :: lumped_mass
1129
type(csr_matrix), pointer :: mass_matrix
1130
type(scalar_field) :: inv_lumped_mass, x
1133
call allocate(inv_lumped_mass, A%mesh)
1134
lumped_mass => get_lumped_mass(state, A%mesh)
1135
call invert(lumped_mass, inv_lumped_mass)
1136
call scale(A, inv_lumped_mass)
1137
call deallocate(inv_lumped_mass)
1139
call allocate(x, A%mesh)
1140
mass_matrix => get_mass_matrix(state, A%mesh)
1141
call petsc_solve(x, mass_matrix, A, &
1142
trim(option_path)//&
1143
'mass_terms/use_consistent_mass_matrix/')
1148
end subroutine solve_cg_inv_mass
1150
!---------------------------------------------------------------------------------
1152
subroutine solve_cg_inv_mass_vector(state, A, lump, option_path)
1154
type(state_type), intent(inout) :: state
1155
type(vector_field), intent(inout) :: A
1156
logical, intent(in) :: lump
1157
character(len=OPTION_PATH_LEN), intent(in) :: option_path
1159
type(scalar_field), pointer :: lumped_mass
1160
type(csr_matrix), pointer :: mass_matrix
1161
type(scalar_field) :: inv_lumped_mass
1162
type(vector_field) :: x
1165
call allocate(inv_lumped_mass, A%mesh)
1166
lumped_mass => get_lumped_mass(state, A%mesh)
1167
call invert(lumped_mass, inv_lumped_mass)
1168
call scale(A, inv_lumped_mass)
1169
call deallocate(inv_lumped_mass)
1171
call allocate(x, A%dim, A%mesh)
1172
mass_matrix => get_mass_matrix(state, A%mesh)
1173
call petsc_solve(x, mass_matrix, A, &
1174
trim(option_path)//&
1175
'mass_terms/use_consistent_mass_matrix/')
1180
end subroutine solve_cg_inv_mass_vector
1182
!---------------------------------------------------------------------------------
1184
subroutine k_epsilon_check_options
1186
character(len=OPTION_PATH_LEN) :: option_path
1187
character(len=FIELD_NAME_LEN) :: kmsh, emsh, vmsh
1188
integer :: dimension, stat, n_phases, istate
1190
ewrite(1,*) "In keps_check_options"
1192
n_phases = option_count("/material_phase")
1194
if(option_count("/material_phase/subgridscale_parameterisations/k-epsilon") > 1) then
1195
FLExit("The k-epsilon model can only be applied to a single-phase.")
1198
do istate = 0, n_phases-1
1200
option_path = "/material_phase["//int2str(istate)//"]/subgridscale_parameterisations/k-epsilon"
1202
! one dimensional problems not supported
1203
call get_option("/geometry/dimension/", dimension)
1204
if (dimension .eq. 1 .and. have_option(trim(option_path))) then
633
1205
FLExit("k-epsilon model is only supported for dimension > 1")
635
! Don't do k-epsilon if it's not included in the model!
636
if (.not.have_option(trim(option_path))) return
638
! checking for required fields
639
if (.not.have_option(trim(option_path)//"/scalar_field::TurbulentKineticEnergy")) then
640
FLExit("You need TurbulentKineticEnergy field for k-epsilon")
642
if (.not.have_option(trim(option_path)//"/scalar_field::TurbulentDissipation")) then
643
FLExit("You need TurbulentDissipation field for k-epsilon")
645
! Check that TurbulentKineticEnergy and TurbulentDissipation fields are on the same mesh
646
call get_option(trim(option_path)//"/scalar_field::TurbulentKineticEnergy/prognostic/mesh", kmsh)
647
call get_option(trim(option_path)//"/scalar_field::TurbulentDissipation/prognostic/mesh", emsh)
648
call get_option(trim(state%option_path)//"/vector_field::Velocity/prognostic/mesh", vmsh)
649
if(.not. kmsh==emsh .or. .not. kmsh==vmsh .or. .not. emsh==vmsh) then
1207
! Don't do k-epsilon if it's not included in the model!
1208
if (.not.have_option(trim(option_path))) return
1210
! checking for required fields
1211
if (have_option(trim(option_path)//"/scalar_field::TurbulentKineticEnergy/prognostic")) then
1212
! diffusivity is on and diagnostic
1213
if (.not.have_option(trim(option_path)//"/scalar_field::TurbulentKineticEnergy"//&
1214
&"/prognostic/tensor_field::Diffusivity")) then
1215
FLExit("You need TurbulentKineticEnergy Diffusivity field for k-epsilon")
1217
if (.not.have_option(trim(option_path)//&
1218
&"/scalar_field::TurbulentKineticEnergy/prognostic/"//&
1219
&"/tensor_field::Diffusivity/diagnostic/algorithm::Internal")) then
1220
FLExit("You need TurbulentKineticEnergy Diffusivity field set to diagnostic/internal")
1223
if (.not.have_option(trim(option_path)//&
1224
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
1225
&"/scalar_field::Source")) then
1226
FLExit("You need TurbulentKineticEnergy Source field for k-epsilon")
1228
if (.not. have_option(trim(option_path)//&
1229
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
1230
&"/scalar_field::Source/diagnostic/algorithm::Internal")) then
1231
FLExit("You need TurbulentKineticEnergy Source field set to diagnostic/internal")
1234
if (.not.have_option(trim(option_path)//&
1235
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
1236
&"/scalar_field::Absorption")) then
1237
FLExit("You need TurbulentKineticEnergy Absorption field for k-epsilon")
1239
if (.not.have_option(trim(option_path)//&
1240
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
1241
&"/scalar_field::Absorption/diagnostic/algorithm::Internal")) then
1242
FLExit("You need TurbulentKineticEnergy Absorption field set to diagnostic/internal")
1244
else if (have_option(trim(option_path)// &
1245
"/scalar_field::TurbulentKineticEnergy/prescribed")) then
1246
ewrite(0,*) "WARNING: TurbulentKineticEnergy field is prescribed"
1248
FLExit("You need prognostic/prescribed TurbulentKineticEnergy field for k-epsilon")
1250
if (have_option(trim(option_path)//"/scalar_field::TurbulentDissipation/prognostic")) then
1251
! diffusivity is on and diagnostic
1252
if (.not.have_option(trim(option_path)//"/scalar_field::TurbulentDissipation"//&
1253
&"/prognostic/tensor_field::Diffusivity")) then
1254
FLExit("You need TurbulentDissipation Diffusivity field for k-epsilon")
1256
if (.not.have_option(trim(option_path)//&
1257
&"/scalar_field::TurbulentDissipation/prognostic/"//&
1258
&"/tensor_field::Diffusivity/diagnostic/algorithm::Internal")) then
1259
FLExit("You need TurbulentDissipation Diffusivity field set to diagnostic/internal")
1262
if (.not.have_option(trim(option_path)//&
1263
&"/scalar_field::TurbulentDissipation/prognostic"//&
1264
&"/scalar_field::Source")) then
1265
FLExit("You need TurbulentDissipation Source field for k-epsilon")
1267
if (.not. have_option(trim(option_path)//&
1268
&"/scalar_field::TurbulentDissipation/prognostic"//&
1269
&"/scalar_field::Source/diagnostic/algorithm::Internal")) then
1270
FLExit("You need TurbulentDissipation Source field set to diagnostic/internal")
1273
if (.not.have_option(trim(option_path)//&
1274
&"/scalar_field::TurbulentDissipation/prognostic"//&
1275
&"/scalar_field::Absorption")) then
1276
FLExit("You need TurbulentDissipation Absorption field for k-epsilon")
1278
if (.not.have_option(trim(option_path)//&
1279
&"/scalar_field::TurbulentDissipation/prognostic"//&
1280
&"/scalar_field::Absorption/diagnostic/algorithm::Internal")) then
1281
FLExit("You need TurbulentDissipation Absorption field set to diagnostic/internal")
1283
else if (have_option(trim(option_path)// &
1284
"/scalar_field::TurbulentDissipation/prescribed")) then
1285
ewrite(0,*) "WARNING: TurbulentDissipation field is prescribed"
1287
FLExit("You need prognostic/prescribed TurbulentDissipation field for k-epsilon")
1290
! Check that TurbulentKineticEnergy and TurbulentDissipation fields are on the same
1291
! mesh as the velocity
1292
call get_option(trim(option_path)//&
1293
&"/scalar_field::TurbulentKineticEnergy/prognostic/mesh/name", kmsh, stat)
1295
call get_option(trim(option_path)//&
1296
&"/scalar_field::TurbulentKineticEnergy/prescribed/mesh/name", kmsh,&
1299
call get_option(trim(option_path)//&
1300
&"/scalar_field::TurbulentDissipation/prognostic/mesh/name", emsh, stat)
1302
call get_option(trim(option_path)//&
1303
&"/scalar_field::TurbulentDissipation/prescribed/mesh/name", emsh,&
1306
call get_option("/material_phase["//int2str(istate)//"]/vector_field::Velocity/prognostic/mesh/name", vmsh,&
1309
call get_option("/material_phase["//int2str(istate)//"]/vector_field::Velocity/prescribed/mesh/name", vmsh,&
1312
FLExit("You must use a prognostic or prescribed Velocity field")
1315
if(.not. kmsh==emsh .or. .not. kmsh==vmsh .or. .not. emsh==vmsh) then
650
1316
FLExit("You must use the Velocity mesh for TurbulentKineticEnergy and TurbulentDissipation fields")
652
! check that diffusivity is on for the two turbulent fields, and diagnostic
653
if (.not.have_option(trim(option_path)//"/scalar_field::TurbulentKineticEnergy"//&
654
&"/prognostic/tensor_field::Diffusivity")) then
655
FLExit("You need TurbulentKineticEnergy Diffusivity field for k-epsilon")
657
if (.not.have_option(trim(option_path)//&
658
&"/scalar_field::TurbulentKineticEnergy/prognostic/"//&
659
&"/tensor_field::Diffusivity/diagnostic/algorithm::Internal")) then
660
FLExit("You need TurbulentKineticEnergy Diffusivity field set to diagnostic/internal")
662
if (.not.have_option(trim(option_path)//&
663
&"/scalar_field::TurbulentDissipation/prognostic"//&
664
&"/tensor_field::Diffusivity")) then
665
FLExit("You need TurbulentDissipation Diffusivity field for k-epsilon")
667
if (.not.have_option(trim(option_path)//&
668
&"/scalar_field::TurbulentDissipation/prognostic"//&
669
&"/tensor_field::Diffusivity/diagnostic/algorithm::Internal")) then
670
FLExit("You need TurbulentDissipation Diffusivity field set to diagnostic/internal")
673
if (.not.have_option(trim(option_path)//&
674
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
675
&"/scalar_field::Source")) then
676
FLExit("You need TurbulentKineticEnergy Source field for k-epsilon")
678
if (.not. have_option(trim(option_path)//&
679
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
680
&"/scalar_field::Source/diagnostic/algorithm::Internal")&
681
.or. .not. have_option(trim(option_path)//&
682
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
683
&"/scalar_field::Source/prescribed")) then
684
FLExit("You need TurbulentKineticEnergy Source field set to diagnostic/internal or prescribed")
686
if (.not.have_option(trim(option_path)//&
687
&"/scalar_field::TurbulentDissipation/prognostic"//&
688
&"/scalar_field::Source")) then
689
FLExit("You need TurbulentDissipation Source field for k-epsilon")
691
if (.not. have_option(trim(option_path)//&
692
&"/scalar_field::TurbulentDissipation/prognostic"//&
693
&"/scalar_field::Source/diagnostic/algorithm::Internal")&
694
.or. .not. have_option(trim(option_path)//&
695
&"/scalar_field::TurbulentDissipation/prognostic"//&
696
&"/scalar_field::Source/prescribed")) then
697
FLExit("You need TurbulentDissipation Source field set to diagnostic/internal or prescribed")
700
if (.not.have_option(trim(option_path)//&
701
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
702
&"/scalar_field::Absorption")) then
703
FLExit("You need TurbulentKineticEnergy Absorption field for k-epsilon")
705
if (.not.have_option(trim(option_path)//&
706
&"/scalar_field::TurbulentKineticEnergy/prognostic"//&
707
&"/scalar_field::Absorption/diagnostic/algorithm::Internal")) then
708
FLExit("You need TurbulentKineticEnergy Absorption field set to diagnostic/internal")
710
if (.not.have_option(trim(option_path)//&
711
&"/scalar_field::TurbulentDissipation/prognostic"//&
712
&"/scalar_field::Absorption")) then
713
FLExit("You need TurbulentDissipation Absorption field for k-epsilon")
715
if (.not.have_option(trim(option_path)//&
716
&"/scalar_field::TurbulentDissipation/prognostic"//&
717
&"/scalar_field::Absorption/diagnostic/algorithm::Internal")) then
718
FLExit("You need TurbulentDissipation Absorption field set to diagnostic/internal")
720
! check there's a viscosity somewhere
721
if (.not.have_option(trim(state%option_path)//"/vector_field::Velocity/prognostic"//&
722
&"/tensor_field::Viscosity/")) then
1319
! Velocity field options
1320
if (.not.have_option("/material_phase["//int2str(istate)//"]/vector_field::Velocity/prognostic"//&
1321
"/tensor_field::Viscosity/") .and. &
1322
.not.have_option("/material_phase["//int2str(istate)//"]/vector_field::Velocity/prescribed")) then
723
1323
FLExit("Need viscosity switched on under the Velocity field for k-epsilon.")
725
! check that the user has switched Velocity/viscosity to diagnostic
726
if (.not.have_option(trim(state%option_path)//"/vector_field::Velocity/prognostic"//&
727
&"/tensor_field::Viscosity/diagnostic/")) then
1325
! check that the user has switched Velocity/viscosity to diagnostic
1326
if (have_option("/material_phase["//int2str(istate)//"]/vector_field::Velocity/prognostic") .and. &
1327
.not.have_option("/material_phase["//int2str(istate)//"]/vector_field::Velocity/prognostic"//&
1328
"/tensor_field::Viscosity/diagnostic/")) then
728
1329
FLExit("You need to switch the viscosity field under Velocity to diagnostic/internal")
730
! check that the background viscosity is an isotropic constant
731
if (.not.have_option(trim(option_path)//"/tensor_field::BackgroundViscosity/prescribed"//&
732
&"/value::WholeMesh/isotropic/constant")) then
733
FLExit("You need to switch the BackgroundViscosity field to isotropic/constant")
736
end subroutine keps_check_options
738
!------------------------------------------------------------------!
739
! Private subroutines !
740
!------------------------------------------------------------------!
742
subroutine keps_allocate_fields(state)
744
type(state_type), intent(inout) :: state
745
type(vector_field), pointer :: vectorField
747
ewrite(1,*) "In keps_allocate_fields"
749
vectorField => extract_vector_field(state, "Velocity")
751
! allocate some space for the fields we need for calculations
752
call allocate(tke_old, vectorField%mesh, "Old_TKE")
753
call allocate(tke_src_old,vectorField%mesh, "Old_TKE_Source")
754
call zero(tke_old); call zero(tke_src_old)
756
end subroutine keps_allocate_fields
758
!--------------------------------------------------------------------------------!
759
! Only used if bc type == k_epsilon for field. !
760
!--------------------------------------------------------------------------------!
762
subroutine keps_wall_function(field1,field2,positions,u,bg_visc,EV,ele,sele,index,wall_fns,cmu,rhs)
764
type(scalar_field), pointer, intent(in) :: field1, field2, EV
765
type(vector_field), pointer, intent(in) :: positions, u
766
type(tensor_field), pointer, intent(in) :: bg_visc
767
integer, intent(in) :: ele, sele, index
768
character(len=FIELD_NAME_LEN), intent(in) :: wall_fns
769
real, dimension(face_loc(field1,sele)), intent(inout):: rhs
771
type(element_type), pointer :: shape, fshape
772
integer :: i, j, gi, sgi, sloc
773
real :: kappa, h, cmu
774
real, dimension(1,1) :: hb
775
real, dimension(ele_ngi(field1,ele)) :: detwei
776
real, dimension(face_ngi(field1,sele)) :: detwei_bdy, ustar, q_sgin, visc_sgi
777
real, dimension(face_loc(field1,sele)) :: lumpedfmass
778
real, dimension(ele_loc(field1,ele)) :: sqrt_k
779
real, dimension(positions%dim,1) :: n
780
real, dimension(positions%dim,positions%dim) :: G
781
real, dimension(positions%dim,ele_ngi(field1,ele)) :: grad_k
782
real, dimension(positions%dim,face_ngi(field1,sele)) :: normal_bdy, q_sgi, qq_sgin
783
real, dimension(positions%dim,ele_loc(field1,ele)) :: q
784
real, dimension(positions%dim,face_loc(field1,ele)) :: q_s
785
real, dimension(face_loc(field1,ele),face_loc(field1,ele)) :: fmass
786
real, dimension(ele_loc(field1,ele),ele_loc(field1,ele)) :: invmass
787
real, dimension(positions%dim,positions%dim,ele_loc(field1,sele)) :: qq
788
real, dimension(positions%dim,positions%dim,ele_ngi(field1,sele)) :: grad_u, invJ
789
real, dimension(positions%dim,positions%dim,face_loc(field1,sele)) :: qq_s
790
real, dimension(positions%dim,positions%dim,face_ngi(field1,sele)) :: bg_visc_sgi, qq_sgi
791
real, dimension(ele_loc(field1,ele),ele_ngi(field1,ele),positions%dim) :: dshape
793
! Get ids, lists and shape functions
794
sgi = face_ngi(field1, sele) ! no. of gauss points in surface element
795
sloc = face_loc(field1, sele) ! no. of nodes on surface element
796
shape => ele_shape(field1, ele) ! scalar field shape functions in volume element
797
fshape => face_shape(field1, sele) ! scalar field shape functions in surface element
799
! Get shape fn gradients, element/face quadrature weights, and surface normal
800
call transform_to_physical( positions, ele, shape, dshape=dshape, detwei=detwei, invJ=invJ )
801
call transform_facet_to_physical( positions, sele, detwei_f=detwei_bdy, normal=normal_bdy )
803
invmass = shape_shape(shape, shape, detwei)
806
fmass = shape_shape(fshape, fshape, detwei_bdy) !*density_gi*vfrac_gi to be safe?
807
lumpedfmass = sum(fmass, 2)
809
! low Re wall functions for k and epsilon: see e.g. Wilcox (1994)
810
if(wall_fns=="low_Re") then
813
else if (index==2) then
814
bg_visc_sgi = face_val_at_quad(bg_visc,sele)
815
visc_sgi = bg_visc_sgi(1,1,:)
817
! grad(k**0.5) (dim, ngi)
818
sqrt_k = ele_val(field2, ele)
819
sqrt_k = sqrt(abs(sqrt_k))
820
! Lifted from ele_grad_at_quad function:
821
do i=1, positions%dim
822
grad_k(i,:) = matmul(sqrt_k, dshape(:,:,i))
825
! grad(k**0.5) at ele_nodes (dim,loc)
826
q = shape_vector_rhs(shape, grad_k, detwei)
828
q = matmul(q,invmass)
829
! Pick surface nodes (dim,sloc)
830
q_s = q(:,face_local_nodes(field1,sele))
831
q_sgi = matmul(q_s, fshape%n)
833
!dot with surface normal
835
q_sgin(gi) = dot_product(q_sgi(:,gi),normal_bdy(:,gi))
838
! integral of 2*nu*(grad(k**0.5))**2 (sloc)
839
rhs = shape_rhs(fshape, detwei_bdy*q_sgin**2.0*visc_sgi*2.0)
840
rhs = rhs/lumpedfmass
843
! high Re shear-stress wall functions for k and epsilon: see e.g. Wilcox (1994), Mathieu p.360
844
else if(wall_fns=="high_Re") then
845
visc_sgi = face_val_at_quad(EV,sele)
846
grad_u = ele_grad_at_quad(u, ele, dshape)
848
! grad(U) at ele_nodes (dim,dim,loc)
849
qq = shape_tensor_rhs(shape, grad_u, detwei)
853
qq(i,j,:) = matmul(invmass,qq(i,j,:))
855
! Pick surface nodes (dim,dim,sloc)
856
qq_s(i,j,:) = qq(i,j,face_local_nodes(field1,sele))
858
! Get values at surface quadrature (dim,dim,sgi)
859
qq_sgi(i,j,:) = matmul(qq_s(i,j,:), fshape%n)
864
!dot with surface normal (dim,sgi)
865
qq_sgin(:,gi) = matmul(qq_sgi(:,:,gi),normal_bdy(:,gi))
867
! Subtract normal component of velocity, leaving tangent components:
868
qq_sgin(:,gi) = qq_sgin(:,gi)-normal_bdy(:,gi)*dot_product(qq_sgin(:,gi),normal_bdy(:,gi))
870
! Get streamwise component by taking sqrt(grad_n.grad_n). Multiply by eddy viscosity.
871
ustar(gi) = norm2(qq_sgin(:,gi)) * visc_sgi(gi)
875
rhs = shape_rhs(fshape, detwei_bdy*ustar/cmu**0.5)
876
rhs = rhs/lumpedfmass
877
else if (index==2) then
878
! calculate wall-normal element size
879
G = matmul(transpose(invJ(:,:,1)), invJ(:,:,1))
880
n(:,1) = normal_bdy(:,1)
881
hb = 1. / sqrt( matmul(matmul(transpose(n), G), n) )
883
! Von Karman's constant
886
rhs = shape_rhs(fshape, detwei_bdy*ustar**1.5/kappa/h)
887
rhs = rhs/lumpedfmass
890
FLAbort("Unknown wall function option for k_epsilon boundary conditions!")
893
end subroutine keps_wall_function
895
!------------------------------------------------------------------------------------------!
896
! Computes the strain rate for the LES model. Double-dot product results in a scalar: !
897
! t:s = [ t11s11 + t12s21 + t13s31 + t21s12 + t22s22 + t23s32 + t31s13 + t32s23 + t33s33 ] !
898
!------------------------------------------------------------------------------------------!
900
function double_dot_product(du_t, nu)
902
real, dimension(:,:,:), intent(in) :: du_t ! derivative of velocity shape function
903
real, dimension(:,:), intent(in) :: nu ! nonlinear velocity
904
real, dimension( size(du_t,2) ) :: double_dot_product
905
real, dimension(size(du_t,3),size(du_t,3)) :: S, T
906
integer :: dim, ngi, gi, i, j
912
S = matmul( nu, du_t(:,gi,:) )
915
double_dot_product(gi) = 0.
919
double_dot_product(gi) = double_dot_product(gi) + T(i,j) * S(j,i)
924
end function double_dot_product
1332
! Check ScalarEddyViscosity is diagnostic
1333
if (have_option(trim(option_path)//'/scalar_field::ScalarEddyViscosity/prescribed')) then
1334
ewrite(0,*) "WARNING: ScalarEddyViscosity field is prescribed"
1339
end subroutine k_epsilon_check_options
926
1341
end module k_epsilon
added
removed
parameterisation/k_epsilon.F90
