267 subroutine les_model_compute_delta(this)
268 class(les_model_t),
intent(inout) :: this
269 integer :: e, i, j, k
270 integer :: im, ip, jm, jp, km, kp
271 real(kind=rp) :: di, dj, dk, ndim_inv, volume_element
272 integer :: lx_half, ly_half, lz_half
273 character(len=:),
allocatable :: type_string
275 lx_half = this%coef%Xh%lx / 2
276 ly_half = this%coef%Xh%ly / 2
277 lz_half = this%coef%Xh%lz / 2
279 if (this%delta_type .eq.
"elementwise_max")
then
282 do e = 1, this%coef%msh%nelv
283 di = (this%coef%dof%x(lx_half, 1, 1, e) &
284 - this%coef%dof%x(lx_half + 1, 1, 1, e))**2 &
285 + (this%coef%dof%y(lx_half, 1, 1, e) &
286 - this%coef%dof%y(lx_half + 1, 1, 1, e))**2 &
287 + (this%coef%dof%z(lx_half, 1, 1, e) &
288 - this%coef%dof%z(lx_half + 1, 1, 1, e))**2
290 dj = (this%coef%dof%x(1, ly_half, 1, e) &
291 - this%coef%dof%x(1, ly_half + 1, 1, e))**2 &
292 + (this%coef%dof%y(1, ly_half, 1, e) &
293 - this%coef%dof%y(1, ly_half + 1, 1, e))**2 &
294 + (this%coef%dof%z(1, ly_half, 1, e) &
295 - this%coef%dof%z(1, ly_half + 1, 1, e))**2
297 dk = (this%coef%dof%x(1, 1, lz_half, e) &
298 - this%coef%dof%x(1, 1, lz_half + 1, e))**2 &
299 + (this%coef%dof%y(1, 1, lz_half, e) &
300 - this%coef%dof%y(1, 1, lz_half + 1, e))**2 &
301 + (this%coef%dof%z(1, 1, lz_half, e) &
302 - this%coef%dof%z(1, 1, lz_half + 1, e))**2
306 this%delta%x(:,:,:,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
308 else if (this%delta_type .eq.
"elementwise_average")
then
311 do e = 1, this%coef%msh%nelv
312 volume_element = 0.0_rp
313 do k = 1, this%coef%Xh%lx * this%coef%Xh%ly * this%coef%Xh%lz
314 volume_element = volume_element + this%coef%B(k, 1, 1, e)
316 this%delta%x(:,:,:,e) = (volume_element / this%coef%Xh%lx &
317 / this%coef%Xh%ly / this%coef%Xh%lz)**(1.0_rp / 3.0_rp)
319 else if (this%delta_type .eq.
"pointwise")
then
320 do e = 1, this%coef%msh%nelv
321 do k = 1, this%coef%Xh%lz
323 kp = min(this%coef%Xh%lz, k+1)
325 do j = 1, this%coef%Xh%ly
327 jp = min(this%coef%Xh%ly, j+1)
329 do i = 1, this%coef%Xh%lx
331 ip = min(this%coef%Xh%lx, i+1)
333 di = (this%coef%dof%x(ip, j, k, e) - &
334 this%coef%dof%x(im, j, k, e))**2 &
335 + (this%coef%dof%y(ip, j, k, e) - &
336 this%coef%dof%y(im, j, k, e))**2 &
337 + (this%coef%dof%z(ip, j, k, e) - &
338 this%coef%dof%z(im, j, k, e))**2
340 dj = (this%coef%dof%x(i, jp, k, e) - &
341 this%coef%dof%x(i, jm, k, e))**2 &
342 + (this%coef%dof%y(i, jp, k, e) - &
343 this%coef%dof%y(i, jm, k, e))**2 &
344 + (this%coef%dof%z(i, jp, k, e) - &
345 this%coef%dof%z(i, jm, k, e))**2
347 dk = (this%coef%dof%x(i, j, kp, e) - &
348 this%coef%dof%x(i, j, km, e))**2 &
349 + (this%coef%dof%y(i, j, kp, e) - &
350 this%coef%dof%y(i, j, km, e))**2 &
351 + (this%coef%dof%z(i, j, kp, e) - &
352 this%coef%dof%z(i, j, km, e))**2
354 di = sqrt(di) / (ip - im)
355 dj = sqrt(dj) / (jp - jm)
356 dk = sqrt(dk) / (kp - km)
357 this%delta%x(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
364 call neko_type_error(
"delta_type for LES model", &
365 this%delta_type, delta_known_types)
369 if (neko_bcknd_device .eq. 1)
then
370 call device_memcpy(this%delta%x, this%delta%x_d, this%delta%dof%size(),&
371 host_to_device, sync = .false.)
372 call this%coef%gs_h%op(this%delta%x, this%delta%dof%size(), gs_op_add)
373 call device_col2(this%delta%x_d, this%coef%mult_d, this%delta%dof%size())
375 call this%coef%gs_h%op(this%delta%x, this%delta%dof%size(), gs_op_add)
376 call col2(this%delta%x, this%coef%mult, this%delta%dof%size())