Neko  0.9.99
A portable framework for high-order spectral element flow simulations
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 !
35 module wall_model
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
42  use neko_config, only : neko_bcknd_device
44  use vector, only : vector_t
45  use utils, only : neko_error, nonlinear_index
46  use math, only : glmin, glmax
47  use comm, only : pe_rank
48  use logger, only : neko_log, neko_log_debug
49 
50  implicit none
51  private
52 
54  type, abstract, public :: wall_model_t
56  type(coef_t), pointer :: coef => null()
58  type(dofmap_t), pointer :: dof => null()
60  integer, pointer :: msk(:) => null()
62  integer, pointer :: facet(:) => null()
64  real(kind=rp), allocatable :: tau_x(:)
66  real(kind=rp), allocatable :: tau_y(:)
68  real(kind=rp), allocatable :: tau_z(:)
70  type(vector_t) :: n_x
72  type(vector_t) :: n_y
74  type(vector_t) :: n_z
76  integer, allocatable :: ind_r(:)
78  integer, allocatable :: ind_s(:)
80  integer, allocatable :: ind_t(:)
82  integer, allocatable :: ind_e(:)
84  type(vector_t) :: h
86  integer :: h_index = 0
88  integer :: n_nodes = 0
90  real(kind=rp) :: nu = 0_rp
92  type(field_t), pointer :: tau_field => null()
93  contains
95  procedure, pass(this) :: init_base => wall_model_init_base
97  procedure, pass(this) :: free_base => wall_model_free_base
99  procedure(wall_model_init), pass(this), deferred :: init
101  procedure(wall_model_free), pass(this), deferred :: free
103  procedure(wall_model_compute), pass(this), deferred :: compute
105  procedure, pass(this) :: find_points => wall_model_find_points
106  end type wall_model_t
107 
108  abstract interface
109 
112  subroutine wall_model_compute(this, t, tstep)
113  import wall_model_t, rp
114  class(wall_model_t), intent(inout) :: this
115  real(kind=rp), intent(in) :: t
116  integer, intent(in) :: tstep
117  end subroutine wall_model_compute
118  end interface
119 
120  abstract interface
121 
128  subroutine wall_model_init(this, coef, msk, facet, nu, h_index, json)
129  import wall_model_t, json_file, dofmap_t, coef_t, rp
130  class(wall_model_t), intent(inout) :: this
131  type(coef_t), intent(in) :: coef
132  integer, intent(in) :: msk(:)
133  integer, intent(in) :: facet(:)
134  real(kind=rp), intent(in) :: nu
135  integer, intent(in) :: h_index
136  type(json_file), intent(inout) :: json
137  end subroutine wall_model_init
138  end interface
139 
140  abstract interface
141 
142  subroutine wall_model_free(this)
143  import wall_model_t
144  class(wall_model_t), intent(inout) :: this
145  end subroutine wall_model_free
146  end interface
147 
148  interface
149 
157  module subroutine wall_model_factory(object, coef, msk, facet, nu, &
158  json)
159  class(wall_model_t), allocatable, intent(inout) :: object
160  type(coef_t), intent(in) :: coef
161  integer, intent(in) :: msk(:)
162  integer, intent(in) :: facet(:)
163  real(kind=rp), intent(in) :: nu
164  type(json_file), intent(inout) :: json
165  end subroutine wall_model_factory
166  end interface
167 
168  public :: wall_model_factory
169 
170 contains
177  subroutine wall_model_init_base(this, coef, msk, facet, nu, index)
178  class(wall_model_t), intent(inout) :: this
179  type(coef_t), target, intent(in) :: coef
180  integer, target, intent(in) :: msk(0:)
181  integer, target, intent(in) :: facet(0:)
182  real(kind=rp), intent(in) :: nu
183  integer, intent(in) :: index
184 
185  call this%free_base
186 
187  this%coef => coef
188  this%dof => coef%dof
189  this%msk(0:msk(0)) => msk
190  this%facet(0:msk(0)) => facet
191  this%nu = nu
192  this%h_index = index
193 
194  call neko_field_registry%add_field(this%dof, "tau", &
195  ignore_existing = .true.)
196 
197  this%tau_field => neko_field_registry%get_field("tau")
198 
199  allocate(this%tau_x(this%msk(0)))
200  allocate(this%tau_y(this%msk(0)))
201  allocate(this%tau_z(this%msk(0)))
202 
203  allocate(this%ind_r(this%msk(0)))
204  allocate(this%ind_s(this%msk(0)))
205  allocate(this%ind_t(this%msk(0)))
206  allocate(this%ind_e(this%msk(0)))
207 
208  call this%h%init(this%msk(0))
209  call this%n_x%init(this%msk(0))
210  call this%n_y%init(this%msk(0))
211  call this%n_z%init(this%msk(0))
212 
213  call this%find_points
214 
215  end subroutine wall_model_init_base
216 
218  subroutine wall_model_free_base(this)
219  class(wall_model_t), intent(inout) :: this
220 
221  nullify(this%coef)
222  nullify(this%msk)
223  nullify(this%facet)
224  nullify(this%tau_field)
225 
226  if (allocated(this%tau_x)) then
227  deallocate(this%tau_x)
228  end if
229  if (allocated(this%tau_y)) then
230  deallocate(this%tau_y)
231  end if
232  if (allocated(this%tau_z)) then
233  deallocate(this%tau_z)
234  end if
235  if (allocated(this%ind_r)) then
236  deallocate(this%ind_r)
237  end if
238 
239  call this%h%free()
240  call this%n_x%free()
241  call this%n_y%free()
242  call this%n_z%free()
243  end subroutine wall_model_free_base
244 
246  subroutine wall_model_find_points(this)
247  class(wall_model_t), intent(inout) :: this
248  integer :: n_nodes, fid, idx(4), i, linear
249  real(kind=rp) :: normal(3), p(3), x, y, z, xw, yw, zw, magp
250  real(kind=rp) :: hmin, hmax
251 
252  n_nodes = this%msk(0)
253  this%n_nodes = n_nodes
254 
255  do i = 1, n_nodes
256  linear = this%msk(i)
257  fid = this%facet(i)
258  idx = nonlinear_index(linear, this%coef%Xh%lx, this%coef%Xh%lx,&
259  this%coef%Xh%lx)
260  normal = this%coef%get_normal(idx(1), idx(2), idx(3), idx(4), fid)
261 
262  this%n_x%x(i) = normal(1)
263  this%n_y%x(i) = normal(2)
264  this%n_z%x(i) = normal(3)
265 
266  ! inward normal
267  normal = -normal
268 
269  select case (fid)
270  case (1)
271  this%ind_r(i) = idx(1) + this%h_index
272  this%ind_s(i) = idx(2)
273  this%ind_t(i) = idx(3)
274  case (2)
275  this%ind_r(i) = idx(1) - this%h_index
276  this%ind_s(i) = idx(2)
277  this%ind_t(i) = idx(3)
278  case (3)
279  this%ind_r(i) = idx(1)
280  this%ind_s(i) = idx(2) + this%h_index
281  this%ind_t(i) = idx(3)
282  case (4)
283  this%ind_r(i) = idx(1)
284  this%ind_s(i) = idx(2) - this%h_index
285  this%ind_t(i) = idx(3)
286  case (5)
287  this%ind_r(i) = idx(1)
288  this%ind_s(i) = idx(2)
289  this%ind_t(i) = idx(3) + this%h_index
290  case (6)
291  this%ind_r(i) = idx(1)
292  this%ind_s(i) = idx(2)
293  this%ind_t(i) = idx(3) - this%h_index
294  case default
295  call neko_error("The face index is not correct ")
296  end select
297  this%ind_e(i) = idx(4)
298 
299  ! Location of the wall node
300  xw = this%dof%x(idx(1), idx(2), idx(3), idx(4))
301  yw = this%dof%y(idx(1), idx(2), idx(3), idx(4))
302  zw = this%dof%z(idx(1), idx(2), idx(3), idx(4))
303 
304  ! Location of the sampling point
305  x = this%dof%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
306  this%ind_e(i))
307  y = this%dof%y(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
308  this%ind_e(i))
309  z = this%dof%z(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
310  this%ind_e(i))
311 
312 
313  ! Vector from the sampling point to the wall
314  p(1) = x - xw
315  p(2) = y - yw
316  p(3) = z - zw
317 
318  ! Total distance to the sampling point
319  magp = sqrt(p(1)**2 + p(2)**2 + p(3)**2)
320 
321  ! Project on the normal direction to get h
322  this%h%x(i) = p(1)*normal(1) + p(2)*normal(2) + p(3)*normal(3)
323 
324  ! Look at how much the total distance distance from the normal and warn
325  ! if significant
326  if ((this%h%x(i) - magp) / magp > 0.1 &
327  .and. (neko_log%level_ .eq. neko_log_debug)) then
328  write(*,*) "Significant missalignment between wall normal and &
329  & sampling point direction at wall node", xw, yw, zw
330  end if
331  end do
332 
333  hmin = glmin(this%h%x, n_nodes)
334 ! hmax = glmax(this%h%x, n_nodes)
335 !
336 ! if (pe_rank .eq. 0) then
337 ! write(*, "(A, F10.4, F10.4)") " h min / max:", hmin, hmax
338 ! end if
339 
340 
341 
342  if (neko_bcknd_device .eq. 1) then
343  call device_memcpy(this%h%x, this%h%x_d, n_nodes, host_to_device,&
344  sync = .false.)
345  call device_memcpy(this%n_x%x, this%n_x%x_d, n_nodes, host_to_device, &
346  sync = .false.)
347  call device_memcpy(this%n_y%x, this%n_y%x_d, n_nodes, host_to_device, &
348  sync = .false.)
349  call device_memcpy(this%n_z%x, this%n_z%x_d, n_nodes, host_to_device, &
350  sync = .true.)
351  end if
352  end subroutine wall_model_find_points
353 
354 end module wall_model
__device__ void nonlinear_index(const int idx, const int lx, int *index)
Copy data between host and device (or device and device)
Definition: device.F90:51
Compute wall shear stress.
Definition: wall_model.f90:112
Common constructor.
Definition: wall_model.f90:128
Coefficients.
Definition: coef.f90:34
Definition: comm.F90:1
integer pe_rank
MPI rank.
Definition: comm.F90:28
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
integer, parameter, public host_to_device
Definition: device.F90:47
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
Logging routines.
Definition: log.f90:34
integer, parameter, public neko_log_debug
Debug log level.
Definition: log.f90:73
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
Definition: math.f90:60
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
Definition: math.f90:377
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Definition: math.f90:407
Build configurations.
Definition: neko_config.f90:34
integer, parameter neko_bcknd_device
Definition: neko_config.f90:44
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 wall_model_t.
Definition: wall_model.f90:35
subroutine wall_model_init_base(this, coef, msk, facet, nu, index)
Wall model factory. Both constructs and initializes the object.
Definition: wall_model.f90:178
subroutine wall_model_find_points(this)
Find sampling points based on the requested index.
Definition: wall_model.f90:247
subroutine wall_model_free_base(this)
Destructor for the wall_model_t (base) class.
Definition: wall_model.f90:219
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
Base abstract type for wall-stress models for wall-modelled LES.
Definition: wall_model.f90:54