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 / &
317 (this%coef%Xh%lx - 1.0_rp) / &
318 (this%coef%Xh%ly - 1.0_rp) / &
319 (this%coef%Xh%lz - 1.0_rp) ) ** (1.0_rp / 3.0_rp)
321 else if (this%delta_type .eq.
"pointwise")
then
322 do e = 1, this%coef%msh%nelv
323 do k = 1, this%coef%Xh%lz
325 kp = min(this%coef%Xh%lz, k+1)
327 do j = 1, this%coef%Xh%ly
329 jp = min(this%coef%Xh%ly, j+1)
331 do i = 1, this%coef%Xh%lx
333 ip = min(this%coef%Xh%lx, i+1)
335 di = (this%coef%dof%x(ip, j, k, e) - &
336 this%coef%dof%x(im, j, k, e))**2 &
337 + (this%coef%dof%y(ip, j, k, e) - &
338 this%coef%dof%y(im, j, k, e))**2 &
339 + (this%coef%dof%z(ip, j, k, e) - &
340 this%coef%dof%z(im, j, k, e))**2
342 dj = (this%coef%dof%x(i, jp, k, e) - &
343 this%coef%dof%x(i, jm, k, e))**2 &
344 + (this%coef%dof%y(i, jp, k, e) - &
345 this%coef%dof%y(i, jm, k, e))**2 &
346 + (this%coef%dof%z(i, jp, k, e) - &
347 this%coef%dof%z(i, jm, k, e))**2
349 dk = (this%coef%dof%x(i, j, kp, e) - &
350 this%coef%dof%x(i, j, km, e))**2 &
351 + (this%coef%dof%y(i, j, kp, e) - &
352 this%coef%dof%y(i, j, km, e))**2 &
353 + (this%coef%dof%z(i, j, kp, e) - &
354 this%coef%dof%z(i, j, km, e))**2
356 di = sqrt(di) / (ip - im)
357 dj = sqrt(dj) / (jp - jm)
358 dk = sqrt(dk) / (kp - km)
359 this%delta%x(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
366 call neko_type_error(
"delta_type for LES model", &
367 this%delta_type, delta_known_types)
371 if (neko_bcknd_device .eq. 1)
then
372 call device_memcpy(this%delta%x, this%delta%x_d, this%delta%dof%size(),&
373 host_to_device, sync = .false.)
374 call this%coef%gs_h%op(this%delta%x, this%delta%dof%size(), gs_op_add)
375 call device_col2(this%delta%x_d, this%coef%mult_d, this%delta%dof%size())
377 call this%coef%gs_h%op(this%delta%x, this%delta%dof%size(), gs_op_add)
378 call col2(this%delta%x, this%coef%mult, this%delta%dof%size())