Neko  0.9.0
A portable framework for high-order spectral element flow simulations
flow_ic.f90
Go to the documentation of this file.
1 ! Copyright (c) 2021, 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 !
34 module flow_ic
35  use num_types, only : rp
36  use logger, only: neko_log, log_size
37  use gather_scatter, only : gs_t, gs_op_add
38  use neko_config, only : neko_bcknd_device
42  use field, only : field_t
45  use coefs, only : coef_t
46  use math, only : col2, cfill, cfill_mask
48  use user_intf, only : useric
49  use json_module, only : json_file
51  use point_zone, only: point_zone_t
54  use fld_file, only: fld_file_t
55  use file, only: file_t
57  use interpolation, only: interpolator_t
58  use space, only: space_t, gll
59  implicit none
60  private
61 
62  interface set_flow_ic
63  module procedure set_flow_ic_int, set_flow_ic_usr
64  end interface set_flow_ic
65 
66  public :: set_flow_ic
67 
68 contains
69 
71  subroutine set_flow_ic_int(u, v, w, p, coef, gs, type, params)
72  type(field_t), intent(inout) :: u
73  type(field_t), intent(inout) :: v
74  type(field_t), intent(inout) :: w
75  type(field_t), intent(inout) :: p
76  type(coef_t), intent(in) :: coef
77  type(gs_t), intent(inout) :: gs
78  character(len=*) :: type
79  type(json_file), intent(inout) :: params
80  real(kind=rp) :: delta, tol
81  real(kind=rp), allocatable :: uinf(:)
82  real(kind=rp), allocatable :: zone_value(:)
83  character(len=:), allocatable :: read_str
84  character(len=NEKO_FNAME_LEN) :: fname, mesh_fname
85  logical :: interpolate
86 
87 
88  !
89  ! Uniform (Uinf, Vinf, Winf)
90  !
91  if (trim(type) .eq. 'uniform') then
92 
93  call json_get(params, 'case.fluid.initial_condition.value', uinf)
94  call set_flow_ic_uniform(u, v, w, uinf)
95 
96  !
97  ! Blasius boundary layer
98  !
99  else if (trim(type) .eq. 'blasius') then
100 
101  call json_get(params, 'case.fluid.blasius.delta', delta)
102  call json_get(params, 'case.fluid.blasius.approximation', &
103  read_str)
104  call json_get(params, 'case.fluid.blasius.freestream_velocity', uinf)
105 
106  call set_flow_ic_blasius(u, v, w, delta, uinf, read_str)
107 
108  !
109  ! Point zone initial condition
110  !
111  else if (trim(type) .eq. 'point_zone') then
112 
113  call json_get(params, 'case.fluid.initial_condition.base_value', uinf)
114  call json_get(params, 'case.fluid.initial_condition.zone_name', &
115  read_str)
116  call json_get(params, 'case.fluid.initial_condition.zone_value', &
117  zone_value)
118 
119  call set_flow_ic_point_zone(u, v, w, uinf, read_str, zone_value)
120 
121  !
122  ! Field initial condition (from fld file)
123  !
124  else if (trim(type) .eq. 'field') then
125 
126  call json_get(params, 'case.fluid.initial_condition.file_name', &
127  read_str)
128  fname = trim(read_str)
129  call json_get_or_default(params, &
130  'case.fluid.initial_condition.interpolate', interpolate, &
131  .false.)
132  call json_get_or_default(params, &
133  'case.fluid.initial_condition.tolerance', tol, 0.000001_rp)
134  call json_get_or_default(params, &
135  'case.fluid.initial_condition.mesh_file_name', read_str, &
136  "none")
137  mesh_fname = trim(read_str)
138 
139  call set_flow_ic_fld(u, v, w, p, fname, interpolate, tol, mesh_fname)
140 
141  else
142  call neko_error('Invalid initial condition')
143  end if
144 
145  call set_flow_ic_common(u, v, w, p, coef, gs)
146 
147  end subroutine set_flow_ic_int
148 
150  subroutine set_flow_ic_usr(u, v, w, p, coef, gs, usr_ic, params)
151  type(field_t), intent(inout) :: u
152  type(field_t), intent(inout) :: v
153  type(field_t), intent(inout) :: w
154  type(field_t), intent(inout) :: p
155  type(coef_t), intent(in) :: coef
156  type(gs_t), intent(inout) :: gs
157  procedure(useric) :: usr_ic
158  type(json_file), intent(inout) :: params
159 
160 
161  call neko_log%message("Type: user")
162  call usr_ic(u, v, w, p, params)
163 
164  call set_flow_ic_common(u, v, w, p, coef, gs)
165 
166  end subroutine set_flow_ic_usr
167 
168  subroutine set_flow_ic_common(u, v, w, p, coef, gs)
169  type(field_t), intent(inout) :: u
170  type(field_t), intent(inout) :: v
171  type(field_t), intent(inout) :: w
172  type(field_t), intent(inout) :: p
173  type(coef_t), intent(in) :: coef
174  type(gs_t), intent(inout) :: gs
175  integer :: n
176 
177  n = u%dof%size()
178 
179  if (neko_bcknd_device .eq. 1) then
180  call device_memcpy(u%x, u%x_d, n, &
181  host_to_device, sync = .false.)
182  call device_memcpy(v%x, v%x_d, n, &
183  host_to_device, sync = .false.)
184  call device_memcpy(w%x, w%x_d, n, &
185  host_to_device, sync = .false.)
186  end if
187 
188  ! Ensure continuity across elements for initial conditions
189  call gs%op(u%x, u%dof%size(), gs_op_add)
190  call gs%op(v%x, v%dof%size(), gs_op_add)
191  call gs%op(w%x, w%dof%size(), gs_op_add)
192 
193  if (neko_bcknd_device .eq. 1) then
194  call device_col2(u%x_d, coef%mult_d, u%dof%size())
195  call device_col2(v%x_d, coef%mult_d, v%dof%size())
196  call device_col2(w%x_d, coef%mult_d, w%dof%size())
197  else
198  call col2(u%x, coef%mult, u%dof%size())
199  call col2(v%x, coef%mult, v%dof%size())
200  call col2(w%x, coef%mult, w%dof%size())
201  end if
202 
203  end subroutine set_flow_ic_common
204 
206  subroutine set_flow_ic_uniform(u, v, w, uinf)
207  type(field_t), intent(inout) :: u
208  type(field_t), intent(inout) :: v
209  type(field_t), intent(inout) :: w
210  real(kind=rp), intent(in) :: uinf(3)
211  integer :: n, i
212  character(len=LOG_SIZE) :: log_buf
213 
214  call neko_log%message("Type : uniform")
215  write (log_buf, '(A, 3(ES12.6, A))') "Value: [", (uinf(i), ", ", i=1, 2), &
216  uinf(3), "]"
217  call neko_log%message(log_buf)
218 
219  u = uinf(1)
220  v = uinf(2)
221  w = uinf(3)
222  n = u%dof%size()
223  if (neko_bcknd_device .eq. 1) then
224  call cfill(u%x, uinf(1), n)
225  call cfill(v%x, uinf(2), n)
226  call cfill(w%x, uinf(3), n)
227  end if
228 
229  end subroutine set_flow_ic_uniform
230 
233  subroutine set_flow_ic_blasius(u, v, w, delta, uinf, type)
234  type(field_t), intent(inout) :: u
235  type(field_t), intent(inout) :: v
236  type(field_t), intent(inout) :: w
237  real(kind=rp), intent(in) :: delta
238  real(kind=rp), intent(in) :: uinf(3)
239  character(len=*), intent(in) :: type
240  procedure(blasius_profile), pointer :: bla => null()
241  integer :: i
242  character(len=LOG_SIZE) :: log_buf
243 
244  call neko_log%message("Type : blasius")
245  write (log_buf, '(A,ES12.6)') "delta : ", delta
246  call neko_log%message(log_buf)
247  call neko_log%message("Approximation : " // trim(type))
248  write (log_buf, '(A,"[",2(ES12.6,","),ES12.6,"]")') "Value : ", &
249  uinf(1), uinf(2), uinf(3)
250  call neko_log%message(log_buf)
251 
252  select case (trim(type))
253  case ('linear')
254  bla => blasius_linear
255  case ('quadratic')
256  bla => blasius_quadratic
257  case ('cubic')
258  bla => blasius_cubic
259  case ('quartic')
260  bla => blasius_quartic
261  case ('sin')
262  bla => blasius_sin
263  case default
264  call neko_error('Invalid Blasius approximation')
265  end select
266 
267  if ((uinf(1) .gt. 0.0_rp) .and. (uinf(2) .eq. 0.0_rp) &
268  .and. (uinf(3) .eq. 0.0_rp)) then
269  do i = 1, u%dof%size()
270  u%x(i,1,1,1) = bla(u%dof%z(i,1,1,1), delta, uinf(1))
271  v%x(i,1,1,1) = 0.0_rp
272  w%x(i,1,1,1) = 0.0_rp
273  end do
274  else if ((uinf(1) .eq. 0.0_rp) .and. (uinf(2) .gt. 0.0_rp) &
275  .and. (uinf(3) .eq. 0.0_rp)) then
276  do i = 1, u%dof%size()
277  u%x(i,1,1,1) = 0.0_rp
278  v%x(i,1,1,1) = bla(u%dof%x(i,1,1,1), delta, uinf(2))
279  w%x(i,1,1,1) = 0.0_rp
280  end do
281  else if ((uinf(1) .eq. 0.0_rp) .and. (uinf(2) .eq. 0.0_rp) &
282  .and. (uinf(3) .gt. 0.0_rp)) then
283  do i = 1, u%dof%size()
284  u%x(i,1,1,1) = 0.0_rp
285  v%x(i,1,1,1) = 0.0_rp
286  w%x(i,1,1,1) = bla(u%dof%y(i,1,1,1), delta, uinf(3))
287  end do
288  end if
289 
290  end subroutine set_flow_ic_blasius
291 
301  subroutine set_flow_ic_point_zone(u, v, w, base_value, zone_name, zone_value)
302  type(field_t), intent(inout) :: u
303  type(field_t), intent(inout) :: v
304  type(field_t), intent(inout) :: w
305  real(kind=rp), intent(in), dimension(3) :: base_value
306  character(len=*), intent(in) :: zone_name
307  real(kind=rp), intent(in) :: zone_value(:)
308  character(len=LOG_SIZE) :: log_buf
309 
310  ! Internal variables
311  class(point_zone_t), pointer :: zone
312  integer :: size
313 
314  call neko_log%message("Type : point_zone")
315  write (log_buf, '(A,ES12.6)') "Base value : ", base_value
316  call neko_log%message(log_buf)
317  call neko_log%message("Zone name : " // trim(zone_name))
318  write (log_buf, '(A,"[",2(ES12.6,","),ES12.6," ]")') "Value : ", &
319  zone_value(1), zone_value(2), zone_value(3)
320  call neko_log%message(log_buf)
321 
322  call set_flow_ic_uniform(u, v, w, base_value)
323  size = u%dof%size()
324 
325  zone => neko_point_zone_registry%get_point_zone(trim(zone_name))
326 
327  call cfill_mask(u%x, zone_value(1), size, zone%mask, zone%size)
328  call cfill_mask(v%x, zone_value(2), size, zone%mask, zone%size)
329  call cfill_mask(w%x, zone_value(3), size, zone%mask, zone%size)
330 
331  end subroutine set_flow_ic_point_zone
332 
350  subroutine set_flow_ic_fld(u, v, w, p, file_name, &
351  interpolate, tolerance, mesh_file_name)
352  type(field_t), intent(inout) :: u
353  type(field_t), intent(inout) :: v
354  type(field_t), intent(inout) :: w
355  type(field_t), intent(inout) :: p
356  character(len=*), intent(in) :: file_name
357  logical, intent(in) :: interpolate
358  real(kind=rp), intent(in) :: tolerance
359  character(len=*), intent(inout) :: mesh_file_name
360 
361  character(len=LOG_SIZE) :: log_buf
362  integer :: sample_idx, sample_mesh_idx
363  integer :: last_index
364  type(fld_file_data_t) :: fld_data
365  type(file_t) :: f
366  logical :: mesh_mismatch
367 
368  ! ---- For the mesh to mesh interpolation
369  type(global_interpolation_t) :: global_interp
370  ! -----
371 
372  ! ---- For space to space interpolation
373  type(space_t) :: prev_Xh
374  type(interpolator_t) :: space_interp
375  ! ----
376 
377  call neko_log%message("Type : field")
378  call neko_log%message("File name : " // trim(file_name))
379  write (log_buf, '(A,L1)') "Interpolation : ", interpolate
380  call neko_log%message(log_buf)
381  if (interpolate) then
382  end if
383 
384  ! Extract sample index from the file name
385  sample_idx = extract_fld_file_index(file_name, -1)
386 
387  if (sample_idx .eq. -1) &
388  call neko_error("Invalid file name for the initial condition. The&
389  & file format must be e.g. 'mean0.f00001'")
390 
391  ! Change from "field0.f000*" to "field0.fld" for the fld reader
392  call filename_chsuffix(file_name, file_name, 'fld')
393 
394  call fld_data%init
395  f = file_t(trim(file_name))
396 
397  if (interpolate) then
398 
399  ! If no mesh file is specified, use the default file name
400  if (mesh_file_name .eq. "none") then
401  mesh_file_name = trim(file_name)
402  sample_mesh_idx = sample_idx
403  else
404 
405  ! Extract sample index from the mesh file name
406  sample_mesh_idx = extract_fld_file_index(mesh_file_name, -1)
407 
408  if (sample_mesh_idx .eq. -1) &
409  call neko_error("Invalid file name for the initial condition. &
410 &The file format must be e.g. 'mean0.f00001'")
411 
412  write (log_buf, '(A,ES12.6)') "Tolerance : ", tolerance
413  call neko_log%message(log_buf)
414  write (log_buf, '(A,A)') "Mesh file : ", &
415  trim(mesh_file_name)
416  call neko_log%message(log_buf)
417 
418  end if ! if mesh_file_name .eq. none
419 
420  ! Read the mesh coordinates if they are not in our fld file
421  if (sample_mesh_idx .ne. sample_idx) then
422  call f%set_counter(sample_mesh_idx)
423  call f%read(fld_data)
424  end if
425 
426  end if
427 
428  ! Read the field file containing (u,v,w,p)
429  call f%set_counter(sample_idx)
430  call f%read(fld_data)
431 
432  !
433  ! Check that the data in the fld file matches the current case.
434  ! Note that this is a safeguard and there are corner cases where
435  ! two different meshes have the same dimension and same # of elements
436  ! but this should be enough to cover obvious cases.
437  !
438  mesh_mismatch = (fld_data%glb_nelv .ne. u%msh%glb_nelv .or. &
439  fld_data%gdim .ne. u%msh%gdim)
440 
441  if (mesh_mismatch .and. .not. interpolate) then
442  call neko_error("The fld file must match the current mesh! &
443 &Use 'interpolate': 'true' to enable interpolation.")
444  else if (.not. mesh_mismatch .and. interpolate) then
445  call neko_log%warning("You have activated interpolation but you might &
446 &still be using the same mesh.")
447  end if
448 
449  ! Mesh interpolation if specified
450  if (interpolate) then
451 
452  ! Issue a warning if the mesh is in single precision
453  select type (ft => f%file_type)
454  type is (fld_file_t)
455  if (.not. ft%dp_precision) then
456  call neko_warning("The coordinates read from the field file are &
457 &in single precision.")
458  call neko_log%message("It is recommended to use a mesh in double &
459 &precision for better interpolation results.")
460  call neko_log%message("If the interpolation does not work, you&
461 &can try to increase the tolerance.")
462  end if
463  class default
464  end select
465 
466  ! Generates an interpolator object and performs the point search
467  global_interp = fld_data%generate_interpolator(u%dof, u%msh, &
468  tolerance)
469 
470  ! Evaluate velocities and pressure
471  call global_interp%evaluate(u%x, fld_data%u%x)
472  call global_interp%evaluate(v%x, fld_data%v%x)
473  call global_interp%evaluate(w%x, fld_data%w%x)
474  call global_interp%evaluate(p%x, fld_data%p%x)
475 
476  call global_interp%free
477 
478  else ! No interpolation, but potentially just from different spaces
479 
480  ! Build a space_t object from the data in the fld file
481  call prev_xh%init(gll, fld_data%lx, fld_data%ly, fld_data%lz)
482  call space_interp%init(u%Xh, prev_xh)
483 
484  ! Do the space-to-space interpolation
485  call space_interp%map_host(u%x, fld_data%u%x, fld_data%nelv, u%Xh)
486  call space_interp%map_host(v%x, fld_data%v%x, fld_data%nelv, u%Xh)
487  call space_interp%map_host(w%x, fld_data%w%x, fld_data%nelv, u%Xh)
488  call space_interp%map_host(p%x, fld_data%p%x, fld_data%nelv, u%Xh)
489 
490  call space_interp%free
491 
492  end if
493 
494  ! NOTE: we do not copy (u,v,w) to the GPU since `set_flow_ic_common` does
495  ! the copy for us
496  if (neko_bcknd_device .eq. 1) call device_memcpy(p%x, p%x_d, p%dof%size(), &
497  host_to_device, sync = .false.)
498 
499  call fld_data%free
500 
501  end subroutine set_flow_ic_fld
502 
503 end module flow_ic
Copy data between host and device (or device and device)
Definition: device.F90:51
Abstract interface for computing a Blasius flow profile.
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Definition: json_utils.f90:54
Retrieves a parameter by name or throws an error.
Definition: json_utils.f90:45
Abstract interface for user defined initial conditions.
Definition: user_intf.f90:57
Coefficients.
Definition: coef.f90:34
subroutine, public device_col2(a_d, b_d, n)
Vector multiplication .
subroutine, public device_cfill_mask(a_d, c, size, mask_d, mask_size)
Fill a constant to a masked vector. .
subroutine, public device_cfill(a_d, c, n)
Set all elements to a constant c .
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
integer, parameter, public host_to_device
Definition: device.F90:47
Defines a field.
Definition: field.f90:34
Module for file I/O operations.
Definition: file.f90:34
Simple module to handle fld file series. Provides an interface to the different fields sotred in a fl...
NEKTON fld file format.
Definition: fld_file.f90:35
Initial flow condition.
Definition: flow_ic.f90:34
subroutine set_flow_ic_usr(u, v, w, p, coef, gs, usr_ic, params)
Set intial flow condition (user defined)
Definition: flow_ic.f90:151
subroutine set_flow_ic_point_zone(u, v, w, base_value, zone_name, zone_value)
Set the initial condition of the flow based on a point zone.
Definition: flow_ic.f90:302
subroutine set_flow_ic_int(u, v, w, p, coef, gs, type, params)
Set initial flow condition (builtin)
Definition: flow_ic.f90:72
subroutine set_flow_ic_fld(u, v, w, p, file_name, interpolate, tolerance, mesh_file_name)
Set the initial condition of the flow based on a field. @detail The fields are read from an fld file....
Definition: flow_ic.f90:352
subroutine set_flow_ic_uniform(u, v, w, uinf)
Uniform initial condition.
Definition: flow_ic.f90:207
subroutine set_flow_ic_common(u, v, w, p, coef, gs)
Definition: flow_ic.f90:169
subroutine set_flow_ic_blasius(u, v, w, delta, uinf, type)
Set a Blasius profile as initial condition.
Definition: flow_ic.f90:234
Defines a flow profile.
real(kind=rp) function, public blasius_quadratic(y, delta, u)
Quadratic approximate Blasius Profile .
real(kind=rp) function, public blasius_quartic(y, delta, u)
Quartic approximate Blasius Profile .
real(kind=rp) function, public blasius_sin(y, delta, u)
Sinusoidal approximate Blasius Profile .
real(kind=rp) function, public blasius_cubic(y, delta, u)
Cubic approximate Blasius Profile .
real(kind=rp) function, public blasius_linear(y, delta, u)
Linear approximate Blasius profile .
Gather-scatter.
Implements global_interpolation given a dofmap.
Routines to interpolate between different spaces.
Utilities for retrieving parameters from the case files.
Definition: json_utils.f90:34
Logging routines.
Definition: log.f90:34
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
integer, parameter, public log_size
Definition: log.f90:42
Definition: math.f90:60
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
Definition: math.f90:348
subroutine, public col2(a, b, n)
Vector multiplication .
Definition: math.f90:729
subroutine, public cfill_mask(a, c, size, mask, mask_size)
Fill a constant to a masked vector. .
Definition: math.f90:297
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
type(point_zone_registry_t), target, public neko_point_zone_registry
Global point_zone registry.
Defines a function space.
Definition: space.f90:34
integer, parameter, public gll
Definition: space.f90:48
Interfaces for user interaction with NEKO.
Definition: user_intf.f90:34
Utilities.
Definition: utils.f90:35
integer function, public extract_fld_file_index(fld_filename, default_index)
Extracts the index of a field file. For example, "myfield.f00045" will return 45. If the suffix of th...
Definition: utils.f90:94
integer, parameter, public neko_fname_len
Definition: utils.f90:40
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition: utils.f90:245
subroutine, public filename_chsuffix(fname, new_fname, new_suffix)
Change a filename's suffix.
Definition: utils.f90:77
subroutine, public filename_suffix(fname, suffix)
Extract a filename's suffix.
Definition: utils.f90:70
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:54
Interface for NEKTON fld files.
Definition: fld_file.f90:64
Implements global interpolation for arbitrary points in the domain.
Interpolation between two space::space_t.
Base abstract type for point zones.
Definition: point_zone.f90:47
The function space for the SEM solution fields.
Definition: space.f90:62