Neko  0.9.0
A portable framework for high-order spectral element flow simulations
scalar_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 scalar_ic
35  use gather_scatter, only : gs_t, gs_op_add
36  use neko_config, only : neko_bcknd_device
37  use num_types, only : rp
38  use device_math, only : device_col2
40  use field, only : field_t
43  use coefs, only : coef_t
44  use math, only : col2, cfill, cfill_mask
45  use user_intf, only : useric_scalar
46  use json_module, only : json_file
48  use point_zone, only: point_zone_t
51  use logger, only: neko_log, log_size
53  use fld_file, only: fld_file_t
54  use checkpoint, only: chkp_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_scalar_ic
63  module procedure set_scalar_ic_int, set_scalar_ic_usr
64  end interface set_scalar_ic
65 
66  public :: set_scalar_ic
67 
68 contains
69 
81  subroutine set_scalar_ic_int(s, coef, gs, type, params)
82  type(field_t), intent(inout) :: s
83  type(coef_t), intent(in) :: coef
84  type(gs_t), intent(inout) :: gs
85  character(len=*) :: type
86  type(json_file), intent(inout) :: params
87 
88  ! Variables for retrieving JSON parameters
89  real(kind=rp) :: ic_value
90  character(len=:), allocatable :: read_str
91  character(len=NEKO_FNAME_LEN) :: fname, mesh_fname
92  real(kind=rp) :: zone_value, tol
93  logical :: interpolate
94 
95  if (trim(type) .eq. 'uniform') then
96 
97  call json_get(params, 'case.scalar.initial_condition.value', ic_value)
98  call set_scalar_ic_uniform(s, ic_value)
99 
100  else if (trim(type) .eq. 'point_zone') then
101 
102  call json_get(params, 'case.scalar.initial_condition.base_value', &
103  ic_value)
104  call json_get(params, 'case.scalar.initial_condition.zone_name', &
105  read_str)
106  call json_get(params, 'case.scalar.initial_condition.zone_value', &
107  zone_value)
108 
109  call set_scalar_ic_point_zone(s, ic_value, read_str, zone_value)
110 
111  else if (trim(type) .eq. 'field') then
112 
113  call json_get(params, 'case.scalar.initial_condition.file_name', &
114  read_str)
115  fname = trim(read_str)
116  call json_get_or_default(params, &
117  'case.scalar.initial_condition.interpolate', interpolate, &
118  .false.)
119  call json_get_or_default(params, &
120  'case.scalar.initial_condition.tolerance', tol, 0.000001_rp)
121  call json_get_or_default(params, &
122  'case.scalar.initial_condition.mesh_file_name', read_str, &
123  "none")
124  mesh_fname = trim(read_str)
125 
126  call set_scalar_ic_fld(s, fname, interpolate, tol, mesh_fname)
127 
128  else
129  call neko_error('Invalid initial condition')
130  end if
131 
132  call set_scalar_ic_common(s, coef, gs)
133 
134  end subroutine set_scalar_ic_int
135 
143  subroutine set_scalar_ic_usr(s, coef, gs, usr_ic, params)
144  type(field_t), intent(inout) :: s
145  type(coef_t), intent(in) :: coef
146  type(gs_t), intent(inout) :: gs
147  procedure(useric_scalar) :: usr_ic
148  type(json_file), intent(inout) :: params
149 
150  call neko_log%message("Type: user")
151  call usr_ic(s, params)
152 
153  call set_scalar_ic_common(s, coef, gs)
154 
155  end subroutine set_scalar_ic_usr
156 
163  subroutine set_scalar_ic_common(s, coef, gs)
164  type(field_t), intent(inout) :: s
165  type(coef_t), intent(in) :: coef
166  type(gs_t), intent(inout) :: gs
167  integer :: n
168 
169  n = s%dof%size()
170  if (neko_bcknd_device .eq. 1) then
171  call device_memcpy(s%x, s%x_d, n, &
172  host_to_device, sync = .false.)
173  end if
174 
175  ! Ensure continuity across elements for initial conditions
176  call gs%op(s%x, n, gs_op_add)
177 
178  if (neko_bcknd_device .eq. 1) then
179  call device_col2(s%x_d, coef%mult_d, n)
180  else
181  call col2(s%x, coef%mult, n)
182  end if
183 
184  end subroutine set_scalar_ic_common
185 
190  subroutine set_scalar_ic_uniform(s, ic_value)
191  type(field_t), intent(inout) :: s
192  real(kind=rp), intent(in) :: ic_value
193  integer :: n
194  character(len=LOG_SIZE) :: log_buf
195 
196  call neko_log%message("Type : uniform")
197  write (log_buf, '(A,ES12.6)') "Value: ", ic_value
198  call neko_log%message(log_buf)
199 
200  s = ic_value
201  n = s%dof%size()
202  if (neko_bcknd_device .eq. 1) then
203  call cfill(s%x, ic_value, n)
204  end if
205 
206  end subroutine set_scalar_ic_uniform
207 
215  subroutine set_scalar_ic_point_zone(s, base_value, zone_name, zone_value)
216  type(field_t), intent(inout) :: s
217  real(kind=rp), intent(in) :: base_value
218  character(len=*), intent(in) :: zone_name
219  real(kind=rp), intent(in) :: zone_value
220 
221  ! Internal variables
222  character(len=LOG_SIZE) :: log_buf
223  class(point_zone_t), pointer :: zone
224  integer :: size
225 
226  call neko_log%message("Type : point_zone")
227  write (log_buf, '(A,ES12.6)') "Base value: ", base_value
228  call neko_log%message(log_buf)
229  call neko_log%message("Zone name : " // trim(zone_name))
230  write (log_buf, '(A,ES12.6)') "Zone value: ", zone_value
231  call neko_log%message(log_buf)
232 
233  size = s%dof%size()
234  zone => neko_point_zone_registry%get_point_zone(trim(zone_name))
235 
236  call set_scalar_ic_uniform(s, base_value)
237  call cfill_mask(s%x, zone_value, size, zone%mask, zone%size)
238 
239  end subroutine set_scalar_ic_point_zone
240 
255  subroutine set_scalar_ic_fld(s, file_name, &
256  interpolate, tolerance, mesh_file_name)
257  type(field_t), intent(inout) :: s
258  character(len=*), intent(in) :: file_name
259  logical, intent(in) :: interpolate
260  real(kind=rp), intent(in) :: tolerance
261  character(len=*), intent(inout) :: mesh_file_name
262 
263  character(len=LOG_SIZE) :: log_buf
264  integer :: sample_idx, sample_mesh_idx
265  integer :: last_index
266  type(fld_file_data_t) :: fld_data
267  type(file_t) :: f
268  logical :: mesh_mismatch
269 
270  ! ---- For the mesh to mesh interpolation
271  type(global_interpolation_t) :: global_interp
272  ! -----
273 
274  ! ---- For space to space interpolation
275  type(space_t) :: prev_Xh
276  type(interpolator_t) :: space_interp
277  ! ----
278 
279  call neko_log%message("Type : field")
280  call neko_log%message("File name : " // trim(file_name))
281  write (log_buf, '(A,L1)') "Interpolation : ", interpolate
282  call neko_log%message(log_buf)
283  if (interpolate) then
284  end if
285 
286  ! Extract sample index from the file name
287  sample_idx = extract_fld_file_index(file_name, -1)
288 
289  if (sample_idx .eq. -1) &
290  call neko_error("Invalid file name for the initial condition. The&
291  & file format must be e.g. 'mean0.f00001'")
292 
293  ! Change from "field0.f000*" to "field0.fld" for the fld reader
294  call filename_chsuffix(file_name, file_name, 'fld')
295 
296  call fld_data%init
297  f = file_t(trim(file_name))
298 
299  if (interpolate) then
300 
301  ! If no mesh file is specified, use the default file name
302  if (mesh_file_name .eq. "none") then
303  mesh_file_name = trim(file_name)
304  sample_mesh_idx = sample_idx
305  else
306 
307  ! Extract sample index from the mesh file name
308  sample_mesh_idx = extract_fld_file_index(mesh_file_name, -1)
309 
310  if (sample_mesh_idx .eq. -1) &
311  call neko_error("Invalid file name for the initial condition. &
312 &The file format must be e.g. 'mean0.f00001'")
313 
314  write (log_buf, '(A,ES12.6)') "Tolerance : ", tolerance
315  call neko_log%message(log_buf)
316  write (log_buf, '(A,A)') "Mesh file : ", &
317  trim(mesh_file_name)
318  call neko_log%message(log_buf)
319 
320  end if ! if mesh_file_name .eq. none
321 
322  ! Read the mesh coordinates if they are not in our fld file
323  if (sample_mesh_idx .ne. sample_idx) then
324  call f%set_counter(sample_mesh_idx)
325  call f%read(fld_data)
326  end if
327 
328  end if
329 
330  ! Read the field file containing (u,v,w,p)
331  call f%set_counter(sample_idx)
332  call f%read(fld_data)
333 
334  !
335  ! Check that the data in the fld file matches the current case.
336  ! Note that this is a safeguard and there are corner cases where
337  ! two different meshes have the same dimension and same # of elements
338  ! but this should be enough to cover obvious cases.
339  !
340  mesh_mismatch = (fld_data%glb_nelv .ne. s%msh%glb_nelv .or. &
341  fld_data%gdim .ne. s%msh%gdim)
342 
343  if (mesh_mismatch .and. .not. interpolate) then
344  call neko_error("The fld file must match the current mesh! &
345 &Use 'interpolate': 'true' to enable interpolation.")
346  else if (.not. mesh_mismatch .and. interpolate) then
347  call neko_log%warning("You have activated interpolation but you might &
348 &still be using the same mesh.")
349  end if
350 
351 
352  ! Mesh interpolation if specified
353  if (interpolate) then
354  ! Issue a warning if the mesh is in single precision
355  select type (ft => f%file_type)
356  type is (fld_file_t)
357  if (.not. ft%dp_precision) then
358  call neko_warning("The coordinates read from the field file are &
359 &in single precision.")
360  call neko_log%message("It is recommended to use a mesh in double &
361 &precision for better interpolation results.")
362  call neko_log%message("If the interpolation does not work, you&
363 &can try to increase the tolerance.")
364  end if
365  class default
366  end select
367 
368  ! Generates an interpolator object and performs the point search
369  global_interp = fld_data%generate_interpolator(s%dof, s%msh, tolerance)
370 
371  ! Evaluate scalar
372  call global_interp%evaluate(s%x, fld_data%t%x)
373  call global_interp%free
374 
375  else ! No interpolation, just potentially from different spaces
376 
377  ! Build a space_t object from the data in the fld file
378  call prev_xh%init(gll, fld_data%lx, fld_data%ly, fld_data%lz)
379  call space_interp%init(s%Xh, prev_xh)
380 
381  ! Do the space-to-space interpolation
382  call space_interp%map_host(s%x, fld_data%t%x, fld_data%nelv, s%Xh)
383 
384  call space_interp%free
385 
386  end if
387 
388  call fld_data%free
389 
390  end subroutine set_scalar_ic_fld
391 
392 end module scalar_ic
Copy data between host and device (or device and device)
Definition: device.F90:51
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 scalar initial conditions.
Definition: user_intf.f90:70
Defines a checkpoint.
Definition: checkpoint.f90:34
Coefficients.
Definition: coef.f90:34
subroutine, public device_col2(a_d, b_d, n)
Vector multiplication .
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
integer, parameter, public host_to_device
Definition: device.F90:47
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
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
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.
Scalar initial condition.
Definition: scalar_ic.f90:34
subroutine set_scalar_ic_uniform(s, ic_value)
Uniform initial condition.
Definition: scalar_ic.f90:191
subroutine set_scalar_ic_point_zone(s, base_value, zone_name, zone_value)
Point zone initial condition.
Definition: scalar_ic.f90:216
subroutine set_scalar_ic_common(s, coef, gs)
Set scalar initial condition (common)
Definition: scalar_ic.f90:164
subroutine set_scalar_ic_usr(s, coef, gs, usr_ic, params)
Set scalar intial condition (user defined)
Definition: scalar_ic.f90:144
subroutine set_scalar_ic_fld(s, file_name, interpolate, tolerance, mesh_file_name)
Set the initial condition of the scalar based on a field. @detail The field is read from an fld file....
Definition: scalar_ic.f90:257
subroutine set_scalar_ic_int(s, coef, gs, type, params)
Set scalar initial condition (builtin)
Definition: scalar_ic.f90:82
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