Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
wall_model.f90
Go to the documentation of this file.
1! Copyright (c) 2024, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
33!
36 use num_types, only : rp
37 use field, only : field_t
38 use json_module, only : json_file
40 use dofmap, only : dofmap_t
41 use coefs, only : coef_t
44 use vector, only : vector_t
46 use math, only : glmin, glmax
47 use comm, only : pe_rank
49 use file, only : file_t
51 use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr, c_associated
54 use json_utils, only : json_get
55 implicit none
56 private
57
59 type, abstract, public :: wall_model_t
61 type(coef_t), pointer :: coef => null()
63 type(dofmap_t), pointer :: dof => null()
65 type(field_t), pointer :: mu => null()
67 type(field_t), pointer :: rho => null()
71 character(len=:), allocatable :: scheme_name
73 integer, pointer :: msk(:) => null()
74 type(c_ptr) :: msk_d = c_null_ptr
76 integer, pointer :: facet(:) => null()
78 type(vector_t) :: tau_x
80 type(vector_t) :: tau_y
82 type(vector_t) :: tau_z
84 type(vector_t) :: n_x
86 type(vector_t) :: n_y
88 type(vector_t) :: n_z
90 integer, allocatable :: ind_r(:)
91 type(c_ptr) :: ind_r_d = c_null_ptr
93 integer, allocatable :: ind_s(:)
94 type(c_ptr) :: ind_s_d = c_null_ptr
96 integer, allocatable :: ind_t(:)
97 type(c_ptr) :: ind_t_d = c_null_ptr
99 integer, allocatable :: ind_e(:)
100 type(c_ptr) :: ind_e_d = c_null_ptr
102 type(vector_t) :: h
104 integer :: h_index = 0
106 integer :: n_nodes = 0
108 type(field_t), pointer :: tau_field => null()
109 contains
111 procedure, pass(this) :: init_base => wall_model_init_base
113 procedure, pass(this) :: partial_init_base => wall_model_partial_init_base
115 procedure, pass(this) :: finalize_base => wall_model_finalize_base
117 procedure, pass(this) :: free_base => wall_model_free_base
119 procedure, pass(this) :: compute_mag_field => wall_model_compute_mag_field
121 procedure(wall_model_init), pass(this), deferred :: init
129 procedure(wall_model_partial_init), pass(this), deferred :: partial_init
131 procedure(wall_model_finalize), pass(this), deferred :: finalize
133 procedure(wall_model_free), pass(this), deferred :: free
135 procedure(wall_model_compute), pass(this), deferred :: compute
137 procedure, pass(this) :: find_points => wall_model_find_points
138 end type wall_model_t
139
140 abstract interface
141
144 subroutine wall_model_compute(this, t, tstep)
145 import wall_model_t, rp
146 class(wall_model_t), intent(inout) :: this
147 real(kind=rp), intent(in) :: t
148 integer, intent(in) :: tstep
149 end subroutine wall_model_compute
150 end interface
151
152 abstract interface
153
161 subroutine wall_model_init(this, scheme_name, coef, msk, facet, &
162 h_index, json)
163 import wall_model_t, json_file, dofmap_t, coef_t, rp
164 class(wall_model_t), intent(inout) :: this
165 character(len=*), intent(in) :: scheme_name
166 type(coef_t), intent(in) :: coef
167 integer, intent(in) :: msk(:)
168 integer, intent(in) :: facet(:)
169 integer, intent(in) :: h_index
170 type(json_file), intent(inout) :: json
171 end subroutine wall_model_init
172 end interface
173
174 abstract interface
175
179 subroutine wall_model_partial_init(this, coef, json)
180 import wall_model_t, json_file, dofmap_t, coef_t, rp
181 class(wall_model_t), intent(inout) :: this
182 type(coef_t), intent(in) :: coef
183 type(json_file), intent(inout) :: json
184 end subroutine wall_model_partial_init
185 end interface
186
187 abstract interface
188
191 subroutine wall_model_finalize(this, msk, facet)
192 import wall_model_t
193 class(wall_model_t), intent(inout) :: this
194 integer, intent(in) :: msk(:)
195 integer, intent(in) :: facet(:)
196 end subroutine wall_model_finalize
197 end interface
198
199 abstract interface
200
201 subroutine wall_model_free(this)
202 import wall_model_t
203 class(wall_model_t), intent(inout) :: this
204 end subroutine wall_model_free
205 end interface
206
207 interface
208
217 module subroutine wall_model_factory(object, scheme_name, coef, msk, &
218 facet, json)
219 class(wall_model_t), allocatable, intent(inout) :: object
220 character(len=*), intent(in) :: scheme_name
221 type(coef_t), intent(in) :: coef
222 integer, intent(in) :: msk(:)
223 integer, intent(in) :: facet(:)
224 type(json_file), intent(inout) :: json
225 end subroutine wall_model_factory
226 end interface
227
228 interface
229
232 module subroutine wall_model_allocator(object, type_name)
233 class(wall_model_t), allocatable, intent(inout) :: object
234 character(len=*), intent(in) :: type_name
235 end subroutine wall_model_allocator
236 end interface
237
238 !
239 ! Machinery for injecting user-defined types
240 !
241
245 abstract interface
246 subroutine wall_model_allocate(obj)
247 import wall_model_t
248 class(wall_model_t), allocatable, intent(inout) :: obj
249 end subroutine wall_model_allocate
250 end interface
251
252 interface
253
254 module subroutine register_wall_model(type_name, allocator)
255 character(len=*), intent(in) :: type_name
256 procedure(wall_model_allocate), pointer, intent(in) :: allocator
257 end subroutine register_wall_model
258 end interface
259
260 ! A name-allocator pair for user-defined types. A helper type to define a
261 ! registry of custom allocators.
262 type allocator_entry
263 character(len=20) :: type_name
264 procedure(wall_model_allocate), pointer, nopass :: allocator
265 end type allocator_entry
266
268 type(allocator_entry), allocatable :: wall_model_registry(:)
269
271 integer :: wall_model_registry_size = 0
272
273 public :: wall_model_factory, wall_model_allocator, register_wall_model, &
274 wall_model_allocate
275
276contains
283 subroutine wall_model_init_base(this, scheme_name, coef, msk, facet, index)
284 class(wall_model_t), intent(inout) :: this
285 type(coef_t), target, intent(in) :: coef
286 integer, target, intent(in) :: msk(0:)
287 integer, target, intent(in) :: facet(0:)
288 character(len=*) :: scheme_name
289 integer, intent(in) :: index
290
291 call this%free_base
292
293 this%coef => coef
294 this%dof => coef%dof
295 this%h_index = index
296 this%scheme_name = trim(scheme_name)
297 this%mu => neko_field_registry%get_field_by_name(this%scheme_name // "_mu")
298 this%rho => neko_field_registry%get_field_by_name(this%scheme_name // &
299 "_rho")
300
301 call neko_field_registry%add_field(this%dof, "tau", &
302 ignore_existing = .true.)
303 this%tau_field => neko_field_registry%get_field("tau")
304
305 call this%finalize_base(msk, facet)
306 end subroutine wall_model_init_base
307
312 subroutine wall_model_partial_init_base(this, coef, json)
313 class(wall_model_t), intent(inout) :: this
314 type(coef_t), target, intent(in) :: coef
315 type(json_file), intent(inout) :: json
316
317 call this%free_base()
318
319 this%coef => coef
320 this%dof => coef%dof
321 call json_get(json, "h_index", this%h_index)
322 call json_get(json, "scheme_name", this%scheme_name)
323
324 this%mu => neko_field_registry%get_field_by_name(this%scheme_name // "_mu")
325 this%rho => neko_field_registry%get_field_by_name(this%scheme_name // &
326 "_rho")
327
328 call neko_field_registry%add_field(this%dof, "tau", &
329 ignore_existing = .true.)
330 this%tau_field => neko_field_registry%get_field("tau")
331 end subroutine wall_model_partial_init_base
332
333 subroutine wall_model_finalize_base(this, msk, facet)
334 class(wall_model_t), intent(inout) :: this
335 integer, target, intent(in) :: msk(0:)
336 integer, target, intent(in) :: facet(:)
337
338 this%msk(0:msk(0)) => msk
339 if (neko_bcknd_device .eq. 1) this%msk_d = device_get_ptr(msk)
340 this%facet(0:msk(0)) => facet
341
342 call this%tau_x%init(this%msk(0))
343 call this%tau_y%init(this%msk(0))
344 call this%tau_z%init(this%msk(0))
345
346 allocate(this%ind_r(this%msk(0)))
347 allocate(this%ind_s(this%msk(0)))
348 allocate(this%ind_t(this%msk(0)))
349 allocate(this%ind_e(this%msk(0)))
350
351 call this%h%init(this%msk(0))
352 call this%n_x%init(this%msk(0))
353 call this%n_y%init(this%msk(0))
354 call this%n_z%init(this%msk(0))
355
356 call this%find_points()
357
358 ! Initialize pointers for device
359 if (neko_bcknd_device .eq. 1) then
360 call device_map(this%ind_r, this%ind_r_d, this%n_nodes)
361 call device_map(this%ind_s, this%ind_s_d, this%n_nodes)
362 call device_map(this%ind_t, this%ind_t_d, this%n_nodes)
363 call device_map(this%ind_e, this%ind_e_d, this%n_nodes)
364 call device_memcpy(this%ind_r, this%ind_r_d, this%n_nodes, &
365 host_to_device, sync = .false.)
366 call device_memcpy(this%ind_s, this%ind_s_d, this%n_nodes, &
367 host_to_device, sync = .false.)
368 call device_memcpy(this%ind_t, this%ind_t_d, this%n_nodes, &
369 host_to_device, sync = .false.)
370 call device_memcpy(this%ind_e, this%ind_e_d, this%n_nodes, &
371 host_to_device, sync = .false.)
372 end if
373
374 end subroutine wall_model_finalize_base
375
377 subroutine wall_model_free_base(this)
378 class(wall_model_t), intent(inout) :: this
379
380 nullify(this%coef)
381 nullify(this%msk)
382 nullify(this%facet)
383 nullify(this%tau_field)
384 nullify(this%mu)
385 nullify(this%rho)
386
387 call this%tau_x%free()
388 call this%tau_y%free()
389 call this%tau_z%free()
390
391 if (allocated(this%ind_r)) then
392 deallocate(this%ind_r)
393 end if
394 if (allocated(this%ind_s)) then
395 deallocate(this%ind_s)
396 end if
397 if (allocated(this%ind_t)) then
398 deallocate(this%ind_t)
399 end if
400
401 if (c_associated(this%msk_d)) then
402 call device_free(this%msk_d)
403 end if
404 if (c_associated(this%ind_r_d)) then
405 call device_free(this%ind_r_d)
406 end if
407 if (c_associated(this%ind_s_d)) then
408 call device_free(this%ind_s_d)
409 end if
410 if (c_associated(this%ind_t_d)) then
411 call device_free(this%ind_t_d)
412 end if
413 if (c_associated(this%ind_e_d)) then
414 call device_free(this%ind_e_d)
415 end if
416
417 call this%h%free()
418 call this%n_x%free()
419 call this%n_y%free()
420 call this%n_z%free()
421 end subroutine wall_model_free_base
422
424 subroutine wall_model_find_points(this)
425 class(wall_model_t), intent(inout) :: this
426 integer :: n_nodes, fid, idx(4), i, linear
427 real(kind=rp) :: normal(3), p(3), x, y, z, xw, yw, zw, magp
428 real(kind=rp) :: hmin, hmax
429 type(field_t), pointer :: h_field
430 type(file_t) :: h_file
431 character(len=LOG_SIZE), allocatable :: log_msg
432
433 n_nodes = this%msk(0)
434 this%n_nodes = n_nodes
435
436 call neko_field_registry%add_field(this%coef%dof, "sampling_height", &
437 ignore_existing=.true.)
438
439 h_field => neko_field_registry%get_field_by_name("sampling_height")
440
441 do i = 1, n_nodes
442 linear = this%msk(i)
443 fid = this%facet(i)
444 idx = nonlinear_index(linear, this%coef%Xh%lx, this%coef%Xh%ly,&
445 this%coef%Xh%lz)
446 normal = this%coef%get_normal(idx(1), idx(2), idx(3), idx(4), fid)
447
448 this%n_x%x(i) = normal(1)
449 this%n_y%x(i) = normal(2)
450 this%n_z%x(i) = normal(3)
451
452 ! inward normal
453 normal = -normal
454
455 select case (fid)
456 case (1)
457 this%ind_r(i) = idx(1) + this%h_index
458 this%ind_s(i) = idx(2)
459 this%ind_t(i) = idx(3)
460 case (2)
461 this%ind_r(i) = idx(1) - this%h_index
462 this%ind_s(i) = idx(2)
463 this%ind_t(i) = idx(3)
464 case (3)
465 this%ind_r(i) = idx(1)
466 this%ind_s(i) = idx(2) + this%h_index
467 this%ind_t(i) = idx(3)
468 case (4)
469 this%ind_r(i) = idx(1)
470 this%ind_s(i) = idx(2) - this%h_index
471 this%ind_t(i) = idx(3)
472 case (5)
473 this%ind_r(i) = idx(1)
474 this%ind_s(i) = idx(2)
475 this%ind_t(i) = idx(3) + this%h_index
476 case (6)
477 this%ind_r(i) = idx(1)
478 this%ind_s(i) = idx(2)
479 this%ind_t(i) = idx(3) - this%h_index
480 case default
481 call neko_error("The face index is not correct ")
482 end select
483 this%ind_e(i) = idx(4)
484
485 ! Location of the wall node
486 xw = this%dof%x(idx(1), idx(2), idx(3), idx(4))
487 yw = this%dof%y(idx(1), idx(2), idx(3), idx(4))
488 zw = this%dof%z(idx(1), idx(2), idx(3), idx(4))
489
490 ! Location of the sampling point
491 x = this%dof%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
492 this%ind_e(i))
493 y = this%dof%y(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
494 this%ind_e(i))
495 z = this%dof%z(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
496 this%ind_e(i))
497
498 ! Vector from the sampling point to the wall
499 p(1) = x - xw
500 p(2) = y - yw
501 p(3) = z - zw
502
503 ! Total distance to the sampling point
504 magp = sqrt(p(1)**2 + p(2)**2 + p(3)**2)
505
506 ! Project on the normal direction to get h
507 this%h%x(i) = p(1)*normal(1) + p(2)*normal(2) + p(3)*normal(3)
508
509 h_field%x(linear,1,1,1) = this%h%x(i)
510
511 ! Look at how much the total distance distance from the normal and warn
512 ! if significant
513 if ((this%h%x(i) - magp) / magp > 0.1) then
514 write(log_msg,*) "Significant misalignment between wall normal and"
515 call neko_log%message(log_msg, neko_log_debug)
516 write(log_msg,*) "sampling point direction at wall node", xw, yw, zw
517 call neko_log%message(log_msg, neko_log_debug)
518 end if
519 end do
520
521! hmin = glmin(this%h%x, n_nodes)
522! hmax = glmax(this%h%x, n_nodes)
523! if (pe_rank .eq. 0) then
524! write(*, "(A, F10.4, F10.4)") " h min / max:", hmin, hmax
525! end if
526
527 if (neko_bcknd_device .eq. 1) then
528 call device_memcpy(this%h%x, this%h%x_d, n_nodes, host_to_device,&
529 sync = .false.)
530 call device_memcpy(this%n_x%x, this%n_x%x_d, n_nodes, host_to_device, &
531 sync = .false.)
532 call device_memcpy(this%n_y%x, this%n_y%x_d, n_nodes, host_to_device, &
533 sync = .false.)
534 call device_memcpy(this%n_z%x, this%n_z%x_d, n_nodes, host_to_device, &
535 sync = .true.)
536 end if
537
538 ! Each wall_model bc will do a write unfortunately... But very helpful
539 ! for setup debugging.
540 call h_file%init("sampling_height.fld")
541 call h_file%write(h_field)
542 end subroutine wall_model_find_points
543
545 class(wall_model_t), intent(inout) :: this
546 integer :: i, m
547 real(kind=rp) :: magtau
548
549 m = this%msk(0)
550 if (m > 0) then
551 if (neko_bcknd_device .eq. 1) then
552 call wall_model_compute_mag_field_device(this%tau_x%x_d, &
553 this%tau_y%x_d, &
554 this%tau_z%x_d, &
555 this%tau_field%x_d, &
556 this%msk_d, m)
557 else
558 do i = 1, m
559 magtau = sqrt(this%tau_x%x(i)**2 + &
560 this%tau_y%x(i)**2 + &
561 this%tau_z%x(i)**2)
562 this%tau_field%x(this%msk(i),1,1,1) = magtau
563 end do
564 end if
565 end if
566
567 end subroutine wall_model_compute_mag_field
568
569end module wall_model
__inline__ __device__ void nonlinear_index(const int idx, const int lx, int *index)
Definition bc_utils.h:44
__global__ void wall_model_compute_mag_field(const T *__restrict__ tau_x_d, const T *__restrict__ tau_y_d, const T *__restrict__ tau_z_d, T *__restrict__ tau_field_d, const int *__restrict__ msk_d, const int m)
Return the device pointer for an associated Fortran array.
Definition device.F90:96
Map a Fortran array to a device (allocate and associate)
Definition device.F90:72
Copy data between host and device (or device and device)
Definition device.F90:66
Retrieves a parameter by name or throws an error.
Compute wall shear stress.
Finilzation of partial construction, similar to bc_t
Common constructor.
Partial constructor from JSON, meant to work as the first stage of initialization before the finalize...
Coefficients.
Definition coef.f90:34
Definition comm.F90:1
integer, public pe_rank
MPI rank.
Definition comm.F90:55
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:47
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:214
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Defines a field.
Definition field.f90:34
Module for file I/O operations.
Definition file.f90:34
Utilities for retrieving parameters from the case files.
Logging routines.
Definition log.f90:34
integer, parameter, public neko_log_debug
Debug log level.
Definition log.f90:78
type(log_t), public neko_log
Global log stream.
Definition log.f90:70
integer, parameter, public log_size
Definition log.f90:46
Definition math.f90:60
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
Definition math.f90:522
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Definition math.f90:552
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Utilities.
Definition utils.f90:35
Defines a vector.
Definition vector.f90:34
Implements the device kernel for the wall_model_t type.
subroutine, public wall_model_compute_mag_field_device(tau_x_d, tau_y_d, tau_z_d, tau_field_d, msk_d, m)
Compute the wall shear stress's magnitude on device.
Implements wall_model_t.
subroutine wall_model_finalize_base(this, msk, facet)
subroutine wall_model_init_base(this, scheme_name, coef, msk, facet, index)
Constructor for the wall_model_t (base) class.
subroutine wall_model_partial_init_base(this, coef, json)
Partial initialization based on JSON, prior to knowing the mask and facets.
subroutine wall_model_find_points(this)
Find sampling points based on the requested index.
subroutine wall_model_free_base(this)
Destructor for the wall_model_t (base) class.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
A wrapper around a polymorphic generic_file_t that handles its init. This is essentially a factory fo...
Definition file.f90:55
Base abstract type for wall-stress models for wall-modelled LES.