Neko  0.9.0
A portable framework for high-order spectral element flow simulations
fld_file.f90
Go to the documentation of this file.
1 ! Copyright (c) 2020-2023, 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 !
35 module fld_file
36  use num_types, only: rp, dp, sp, i8
37  use generic_file, only: generic_file_t
38  use field, only: field_t
39  use field_list, only: field_list_t
40  use dofmap, only: dofmap_t
41  use space, only: space_t
42  use structs, only: array_ptr_t
43  use vector, only: vector_t
45  use mean_flow, only: mean_flow_t
47  use vector, only : vector_t
48  use space, only : space_t
49  use mesh, only : mesh_t
51  use utils, only: neko_error
52  use comm
53  use datadist, only: linear_dist_t
54  use math, only: vlmin, vlmax
57  implicit none
58  private
59 
60  real(kind=dp), private, allocatable :: tmp_dp(:)
61  real(kind=sp), private, allocatable :: tmp_sp(:)
62 
64  type, public, extends(generic_file_t) :: fld_file_t
65  logical :: dp_precision = .false.
66  contains
67  procedure :: read => fld_file_read
68  procedure :: write => fld_file_write
69  procedure :: set_precision => fld_file_set_precision
70  end type fld_file_t
71 
72 
73 contains
74 
77  subroutine fld_file_write(this, data, t)
78  class(fld_file_t), intent(inout) :: this
79  class(*), target, intent(in) :: data
80  real(kind=rp), intent(in), optional :: t
81  type(array_ptr_t) :: x, y, z, u, v, w, p, tem
82  real(kind=rp), allocatable, target :: tempo(:)
83  type(mesh_t), pointer :: msh
84  type(dofmap_t), pointer :: dof
85  type(space_t), pointer :: Xh
86  real(kind=dp) :: time
87  character(len= 132) :: hdr
88  character :: rdcode(10)
89  character(len=6) :: id_str
90  character(len= 1024) :: fname
91  character(len= 1024) :: start_field
92  integer :: i, ierr, n, suffix_pos, tslash_pos
93  integer :: lx, ly, lz, lxyz, gdim, glb_nelv, nelv, offset_el
94  integer, allocatable :: idx(:)
95  type(mpi_status) :: status
96  type(mpi_file) :: fh
97  integer (kind=MPI_OFFSET_KIND) :: mpi_offset, byte_offset, temp_offset
98  real(kind=sp), parameter :: test_pattern = 6.54321
99  type(array_ptr_t), allocatable :: scalar_fields(:)
100  logical :: write_mesh, write_velocity, write_pressure, write_temperature
101  integer :: fld_data_size, n_scalar_fields
102 
103  if (present(t)) then
104  time = real(t, dp)
105  else
106  time = 0d0
107  end if
108 
109  nullify(msh)
110  nullify(dof)
111  nullify(xh)
112  n_scalar_fields = 0
113  write_pressure = .false.
114  write_velocity = .false.
115  write_temperature = .false.
116 
117  select type (data)
118  type is (fld_file_data_t)
119  nelv = data%nelv
120  lx = data%lx
121  ly = data%ly
122  lz = data%lz
123  gdim = data%gdim
124  glb_nelv = data%glb_nelv
125  offset_el = data%offset_el
126 
127  if (data%x%n .gt. 0) x%ptr => data%x%x
128  if (data%y%n .gt. 0) y%ptr => data%y%x
129  if (data%z%n .gt. 0) z%ptr => data%z%x
130  if (gdim .eq. 2) z%ptr => data%y%x
131  if (data%u%n .gt. 0) then
132  u%ptr => data%u%x
133  write_velocity = .true.
134  end if
135  if (data%v%n .gt. 0) v%ptr => data%v%x
136  if (data%w%n .gt. 0) w%ptr => data%w%x
137  if (data%p%n .gt. 0) then
138  p%ptr => data%p%x
139  write_pressure = .true.
140  end if
141  if (data%t%n .gt. 0) then
142  write_temperature = .true.
143  tem%ptr => data%t%x
144  end if
145  ! If gdim = 2 and Z-velocity component exists,
146  ! it is stored in last scalar field
147  if (gdim .eq. 2 .and. data%w%n .gt. 0) then
148  n_scalar_fields = data%n_scalars + 1
149  allocate(scalar_fields(n_scalar_fields))
150  do i = 1, n_scalar_fields -1
151  scalar_fields(i)%ptr => data%s(i)%x
152  end do
153  scalar_fields(n_scalar_fields)%ptr => data%w%x
154  else
155  n_scalar_fields = data%n_scalars
156  allocate(scalar_fields(n_scalar_fields+1))
157  do i = 1, n_scalar_fields
158  scalar_fields(i)%ptr => data%s(i)%x
159  end do
160  scalar_fields(n_scalar_fields+1)%ptr => data%w%x
161  end if
162  ! This is very stupid...
163  ! Some compilers cannot handle that these pointers dont point to anything
164  ! (although they are not used) this fixes a segfault due to this.
165  if (nelv .eq. 0) then
166  allocate(tempo(1))
167  x%ptr => tempo
168  y%ptr => tempo
169  z%ptr => tempo
170  u%ptr => tempo
171  v%ptr => tempo
172  w%ptr => tempo
173  p%ptr => tempo
174  tem%ptr => tempo
175  end if
176 
177  allocate(idx(nelv))
178  do i = 1, nelv
179  idx(i) = data%idx(i)
180  end do
181  type is (field_t)
182  p%ptr => data%x(:,1,1,1)
183  dof => data%dof
184  write_pressure = .true.
185  write_velocity = .false.
186  type is (field_list_t)
187  select case (data%size())
188  case (1)
189  p%ptr => data%items(1)%ptr%x(:,1,1,1)
190  write_pressure = .true.
191  write_velocity = .false.
192  case (2)
193  p%ptr => data%items(1)%ptr%x(:,1,1,1)
194  tem%ptr => data%items(2)%ptr%x(:,1,1,1)
195  write_pressure = .true.
196  write_temperature = .true.
197  case (3)
198  u%ptr => data%items(1)%ptr%x(:,1,1,1)
199  v%ptr => data%items(2)%ptr%x(:,1,1,1)
200  w%ptr => data%items(3)%ptr%x(:,1,1,1)
201  write_velocity = .true.
202  case (4)
203  p%ptr => data%items(1)%ptr%x(:,1,1,1)
204  u%ptr => data%items(2)%ptr%x(:,1,1,1)
205  v%ptr => data%items(3)%ptr%x(:,1,1,1)
206  w%ptr => data%items(4)%ptr%x(:,1,1,1)
207  write_pressure = .true.
208  write_velocity = .true.
209  case (5:99)
210  p%ptr => data%items(1)%ptr%x(:,1,1,1)
211  u%ptr => data%items(2)%ptr%x(:,1,1,1)
212  v%ptr => data%items(3)%ptr%x(:,1,1,1)
213  w%ptr => data%items(4)%ptr%x(:,1,1,1)
214  tem%ptr => data%items(5)%ptr%x(:,1,1,1)
215  n_scalar_fields = data%size() - 5
216  allocate(scalar_fields(n_scalar_fields))
217  do i = 1, n_scalar_fields
218  scalar_fields(i)%ptr => data%items(i+5)%ptr%x(:,1,1,1)
219  end do
220  write_pressure = .true.
221  write_velocity = .true.
222  write_temperature = .true.
223  case default
224  call neko_error('This many fields not supported yet, fld_file')
225  end select
226  dof => data%dof(1)
227 
228  type is (mean_flow_t)
229  u%ptr => data%u%mf%x(:,1,1,1)
230  v%ptr => data%v%mf%x(:,1,1,1)
231  w%ptr => data%w%mf%x(:,1,1,1)
232  p%ptr => data%p%mf%x(:,1,1,1)
233  dof => data%u%mf%dof
234  write_pressure = .true.
235  write_velocity = .true.
236  type is (mean_sqr_flow_t)
237  u%ptr => data%uu%mf%x(:,1,1,1)
238  v%ptr => data%vv%mf%x(:,1,1,1)
239  w%ptr => data%ww%mf%x(:,1,1,1)
240  p%ptr => data%pp%mf%x(:,1,1,1)
241  dof => data%pp%mf%dof
242  write_pressure = .true.
243  write_velocity = .true.
244  class default
245  call neko_error('Invalid data')
246  end select
247  ! Fix things for pointers that do not exist in all data types...
248  if (associated(dof)) then
249  x%ptr => dof%x(:,1,1,1)
250  y%ptr => dof%y(:,1,1,1)
251  z%ptr => dof%z(:,1,1,1)
252  msh => dof%msh
253  xh => dof%Xh
254  end if
255 
256  if (associated(msh)) then
257  nelv = msh%nelv
258  glb_nelv = msh%glb_nelv
259  offset_el = msh%offset_el
260  gdim = msh%gdim
261  ! Store global idx of each element
262  allocate(idx(msh%nelv))
263  do i = 1, msh%nelv
264  idx(i) = msh%elements(i)%e%id()
265  end do
266  end if
267 
268  if (associated(xh)) then
269  lx = xh%lx
270  ly = xh%ly
271  lz = xh%lz
272  end if
273 
274  lxyz = lx*ly*lz
275  n = nelv*lxyz
276 
277  if (this%dp_precision) then
278  fld_data_size = mpi_double_precision_size
279  else
280  fld_data_size = mpi_real_size
281  end if
282  if (this%dp_precision) then
283  allocate(tmp_dp(gdim*n))
284  else
285  allocate(tmp_sp(gdim*n))
286  end if
287 
288 
289  !
290  ! Create fld header for NEKTON's multifile output
291  !
292 
293  write_mesh = (this%counter .eq. this%start_counter)
294  call mpi_allreduce(mpi_in_place, write_mesh, 1, &
295  mpi_logical, mpi_lor, neko_comm)
296  call mpi_allreduce(mpi_in_place, write_velocity, 1, &
297  mpi_logical, mpi_lor, neko_comm)
298  call mpi_allreduce(mpi_in_place, write_pressure, 1, &
299  mpi_logical, mpi_lor, neko_comm)
300  call mpi_allreduce(mpi_in_place, write_temperature, 1, &
301  mpi_logical, mpi_lor, neko_comm)
302  call mpi_allreduce(mpi_in_place, n_scalar_fields, 1, &
303  mpi_integer, mpi_max, neko_comm)
304 
305  ! Build rdcode note that for field_t, we only support scalar
306  ! fields at the moment
307  rdcode = ' '
308  i = 1
309  if (write_mesh) then
310  rdcode(i) = 'X'
311  i = i + 1
312  end if
313  if (write_velocity) then
314  rdcode(i) = 'U'
315  i = i + 1
316  end if
317  if (write_pressure) then
318  rdcode(i) = 'P'
319  i = i + 1
320  end if
321  if (write_temperature) then
322  rdcode(i) = 'T'
323  i = i + 1
324  end if
325  if (n_scalar_fields .gt. 0 ) then
326  rdcode(i) = 'S'
327  i = i + 1
328  write(rdcode(i), '(i1)') (n_scalar_fields)/10
329  i = i + 1
330  write(rdcode(i), '(i1)') (n_scalar_fields) - 10*((n_scalar_fields)/10)
331  i = i + 1
332  end if
333 
335  write(hdr, 1) fld_data_size, lx, ly, lz, glb_nelv, glb_nelv,&
336  time, this%counter, 1, 1, (rdcode(i),i = 1, 10)
337 1 format('#std', 1x, i1, 1x, i2, 1x, i2, 1x, i2, 1x, i10, 1x, i10, &
338  1x, e20.13, 1x, i9, 1x, i6, 1x, i6, 1x, 10a)
339 
340  ! Change to NEKTON's fld file format
341  suffix_pos = filename_suffix_pos(this%fname)
342  write(id_str, '(a,i5.5)') 'f', this%counter
343  fname = trim(this%fname(1:suffix_pos-1))//'0.'//id_str
344 
345  call mpi_file_open(neko_comm, trim(fname), &
346  mpi_mode_wronly + mpi_mode_create, mpi_info_null, fh, &
347  ierr)
348 
349  call mpi_file_write_all(fh, hdr, 132, mpi_character, status, ierr)
350  mpi_offset = 132 * mpi_character_size
351 
352  call mpi_file_write_all(fh, test_pattern, 1, mpi_real, status, ierr)
353  mpi_offset = mpi_offset + mpi_real_size
354 
355  byte_offset = mpi_offset + &
356  int(offset_el, i8) * int(mpi_integer_size, i8)
357  call mpi_file_write_at_all(fh, byte_offset, idx, nelv, &
358  mpi_integer, status, ierr)
359  mpi_offset = mpi_offset + int(glb_nelv, i8) * int(mpi_integer_size, i8)
360  deallocate(idx)
361  if (write_mesh) then
362 
363  byte_offset = mpi_offset + int(offset_el, i8) * &
364  (int(gdim*lxyz, i8) * &
365  int(fld_data_size, i8))
366  call fld_file_write_vector_field(this, fh, byte_offset, &
367  x%ptr, y%ptr, z%ptr, &
368  n, gdim, lxyz, nelv)
369  mpi_offset = mpi_offset + int(glb_nelv, i8) * &
370  (int(gdim *lxyz, i8) * &
371  int(fld_data_size, i8))
372  end if
373  if (write_velocity) then
374  byte_offset = mpi_offset + int(offset_el, i8) * &
375  (int(gdim * (lxyz), i8) * int(fld_data_size, i8))
376  call fld_file_write_vector_field(this, fh, byte_offset, &
377  u%ptr, v%ptr, w%ptr, n, gdim, lxyz, nelv)
378 
379  mpi_offset = mpi_offset + int(glb_nelv, i8) * &
380  (int(gdim * (lxyz), i8) * &
381  int(fld_data_size, i8))
382 
383  end if
384 
385  if (write_pressure) then
386  byte_offset = mpi_offset + int(offset_el, i8) * &
387  (int((lxyz), i8) * int(fld_data_size, i8))
388  call fld_file_write_field(this, fh, byte_offset, p%ptr, n)
389  mpi_offset = mpi_offset + int(glb_nelv, i8) * &
390  (int((lxyz), i8) * int(fld_data_size, i8))
391  end if
392 
393  if (write_temperature) then
394  byte_offset = mpi_offset + int(offset_el, i8) * &
395  (int((lxyz), i8) * &
396  int(fld_data_size, i8))
397  call fld_file_write_field(this, fh, byte_offset, tem%ptr, n)
398  mpi_offset = mpi_offset + int(glb_nelv, i8) * &
399  (int((lxyz), i8) * &
400  int(fld_data_size, i8))
401  end if
402 
403  temp_offset = mpi_offset
404 
405  do i = 1, n_scalar_fields
406  ! Without this redundant if statement Cray optimizes this loop to
407  ! Oblivion
408  if (i .eq. 2) then
409  mpi_offset = int(temp_offset, i8) + int(1_i8*glb_nelv, i8) * &
410  (int(lxyz, i8) * int(fld_data_size, i8))
411  end if
412  byte_offset = int(mpi_offset, i8) + int(offset_el, i8) * &
413  (int((lxyz), i8) * &
414  int(fld_data_size, i8))
415  call fld_file_write_field(this, fh, byte_offset, scalar_fields(i)%ptr, n)
416  mpi_offset = int(mpi_offset, i8) + int(glb_nelv, i8) * &
417  (int(lxyz, i8) * &
418  int(fld_data_size, i8))
419  end do
420 
421  if (gdim .eq. 3) then
422 
424  if (write_mesh) then
425  !The offset is:
426  ! mpioff + element_off * 2(min max value)
427  ! * 4(single precision) * gdim(dimensions)
428  byte_offset = int(mpi_offset, i8) + &
429  int(offset_el, i8) * &
430  int(2, i8) * &
431  int(mpi_real_size, i8) * &
432  int(gdim, i8)
433  call fld_file_write_metadata_vector(this, fh, byte_offset, &
434  x%ptr, y%ptr, z%ptr, gdim, lxyz, nelv)
435  mpi_offset = int(mpi_offset, i8) + &
436  int(glb_nelv, i8) * &
437  int(2, i8) * &
438  int(mpi_real_size, i8) * &
439  int(gdim, i8)
440  end if
441 
442  if (write_velocity) then
443  byte_offset = int(mpi_offset, i8) + &
444  int(offset_el, i8) * &
445  int(2, i8) * &
446  int(mpi_real_size, i8) * &
447  int(gdim, i8)
448  call fld_file_write_metadata_vector(this, fh, byte_offset, &
449  u%ptr, v%ptr, w%ptr, gdim, lxyz, nelv)
450  mpi_offset = int(mpi_offset, i8) + &
451  int(glb_nelv, i8) * &
452  int(2, i8) * &
453  int(mpi_real_size, i8) * &
454  int(gdim, i8)
455 
456  end if
457 
458  if (write_pressure) then
459  byte_offset = int(mpi_offset, i8) + &
460  int(offset_el, i8) * &
461  int(2, i8) * &
462  int(mpi_real_size, i8)
463  call fld_file_write_metadata_scalar(this, fh, byte_offset, &
464  p%ptr, lxyz, nelv)
465  mpi_offset = int(mpi_offset, i8) + &
466  int(glb_nelv, i8) * &
467  int(2, i8) * &
468  int(mpi_real_size, i8)
469 
470  end if
471 
472  if (write_temperature) then
473  byte_offset = int(mpi_offset, i8) + &
474  int(offset_el, i8) * &
475  int(2, i8) * &
476  int(mpi_real_size, i8)
477  call fld_file_write_metadata_scalar(this, fh, byte_offset, &
478  tem%ptr, lxyz, nelv)
479  mpi_offset = int(mpi_offset, i8) + &
480  int(glb_nelv, i8) * &
481  int(2, i8) * &
482  int(mpi_real_size, i8)
483 
484  end if
485 
486 
487 
488  temp_offset = mpi_offset
489 
490  do i = 1, n_scalar_fields
491  ! Without this redundant if statement, Cray optimizes this loop to
492  ! Oblivion
493  if (i .eq. 2) then
494  mpi_offset = int(temp_offset, i8) + &
495  int(1_i8*glb_nelv, i8) * &
496  int(2, i8) * &
497  int(mpi_real_size, i8)
498  end if
499 
500  byte_offset = int(mpi_offset, i8) + &
501  int(offset_el, i8) * &
502  int(2, i8) * &
503  int(mpi_real_size, i8)
504  call fld_file_write_metadata_scalar(this, fh, byte_offset, &
505  scalar_fields(i)%ptr, lxyz, nelv)
506  mpi_offset = int(mpi_offset, i8) + &
507  int(glb_nelv, i8) * &
508  int(2, i8) * &
509  int(mpi_real_size, i8)
510  end do
511  end if
512 
513 
514  call mpi_file_sync(fh, ierr)
515  call mpi_file_close(fh, ierr)
516  ! Write metadata file
517  if (pe_rank .eq. 0) then
518  tslash_pos = filename_tslash_pos(this%fname)
519  write(start_field, "(I5,A8)") this%start_counter, '.nek5000'
520  open(unit = 9, &
521  file = trim(this%fname(1:suffix_pos-1)) // &
522  trim(adjustl(start_field)), status = 'replace')
523  write(9, fmt = '(A,A,A)') 'filetemplate: ', &
524  this%fname(tslash_pos+1:suffix_pos-1), '%01d.f%05d'
525  write(9, fmt = '(A,i5)') 'firsttimestep: ', this%start_counter
526  write(9, fmt = '(A,i5)') 'numtimesteps: ', &
527  (this%counter + 1) - this%start_counter
528  write(9, fmt = '(A)') 'type: binary'
529  close(9)
530  end if
531 
532  this%counter = this%counter + 1
533 
534  if (allocated(tmp_dp)) deallocate(tmp_dp)
535  if (allocated(tmp_sp)) deallocate(tmp_sp)
536  end subroutine fld_file_write
537 
538  subroutine fld_file_write_metadata_vector(this, fh, byte_offset, x, y, z, &
539  gdim, lxyz, nelv)
540  class(fld_file_t), intent(inout) :: this
541  type(mpi_file), intent(inout) :: fh
542  integer, intent(in) :: gdim, lxyz, nelv
543  real(kind=rp), intent(in) :: x(lxyz, nelv), y(lxyz, nelv), z(lxyz, nelv)
544  integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
545  integer :: el, j, ierr, nout
546  type(mpi_status) :: status
547  real(kind=sp) :: buffer(2*gdim*nelv)
548 
549  j = 1
550  do el = 1, nelv
551  buffer(j+0) = real(vlmin(x(1, el), lxyz), sp)
552  buffer(j+1) = real(vlmax(x(1, el), lxyz), sp)
553  buffer(j+2) = real(vlmin(y(1, el), lxyz), sp)
554  buffer(j+3) = real(vlmax(y(1, el), lxyz), sp)
555  j = j + 4
556  if (gdim .eq. 3) then
557  buffer(j+0) = real(vlmin(z(1, el), lxyz), sp)
558  buffer(j+1) = real(vlmax(z(1, el), lxyz), sp)
559  j = j + 2
560  end if
561  end do
562 
563  ! write out data
564  nout = 2*gdim*nelv
565 
566  call mpi_file_write_at_all(fh, byte_offset, buffer, nout, &
567  mpi_real, status, ierr)
568 
569  end subroutine fld_file_write_metadata_vector
570 
571  subroutine fld_file_write_metadata_scalar(this, fh, byte_offset, x, lxyz, &
572  nelv)
573  class(fld_file_t), intent(inout) :: this
574  type(mpi_file), intent(inout) :: fh
575  integer, intent(in) :: lxyz, nelv
576  real(kind=rp), intent(in) :: x(lxyz, nelv)
577  integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
578  integer :: el, j, ierr, nout
579  type(mpi_status) :: status
580  real(kind=sp) :: buffer(2*nelv)
581 
582  j = 1
583  do el = 1, nelv
584  buffer(j+0) = real(vlmin(x(1, el), lxyz), sp)
585  buffer(j+1) = real(vlmax(x(1, el), lxyz), sp)
586  j = j + 2
587  end do
588 
589  ! write out data
590  nout = 2*nelv
591 
592  call mpi_file_write_at_all(fh, byte_offset, buffer, nout, &
593  mpi_real, status, ierr)
594 
595  end subroutine fld_file_write_metadata_scalar
596 
597  subroutine fld_file_write_field(this, fh, byte_offset, p, n)
598  class(fld_file_t), intent(inout) :: this
599  type(mpi_file), intent(inout) :: fh
600  integer, intent(inout) :: n
601  real(kind=rp), intent(inout) :: p(n)
602  integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
603  integer :: i, ierr
604  type(mpi_status) :: status
605 
606  if ( this%dp_precision) then
607  do i = 1, n
608  tmp_dp(i) = real(p(i), dp)
609  end do
610 
611  call mpi_file_write_at_all(fh, byte_offset, tmp_dp, n, &
612  mpi_double_precision, status, ierr)
613  else
614  do i = 1, n
615  tmp_sp(i) = real(p(i), sp)
616  end do
617  call mpi_file_write_at_all(fh, byte_offset, tmp_sp, n, &
618  mpi_real, status, ierr)
619  end if
620 
621  end subroutine fld_file_write_field
622 
623  subroutine fld_file_write_vector_field(this, fh, byte_offset, x, y, z, n, &
624  gdim, lxyz, nelv)
625  class(fld_file_t), intent(inout) :: this
626  type(mpi_file), intent(inout) :: fh
627  integer, intent(in) :: n, gdim, lxyz, nelv
628  real(kind=rp), intent(in) :: x(lxyz, nelv), y(lxyz, nelv), z(lxyz, nelv)
629  integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
630  integer :: i, el, j, ierr
631  type(mpi_status) :: status
632 
633  if (this%dp_precision) then
634  i = 1
635  do el = 1, nelv
636  do j = 1, lxyz
637  tmp_dp(i) = real(x(j, el), dp)
638  i = i +1
639  end do
640  do j = 1, lxyz
641  tmp_dp(i) = real(y(j, el), dp)
642  i = i +1
643  end do
644  if (gdim .eq. 3) then
645  do j = 1, lxyz
646  tmp_dp(i) = real(z(j, el), dp)
647  i = i +1
648  end do
649  end if
650  end do
651  call mpi_file_write_at_all(fh, byte_offset, tmp_dp, gdim*n, &
652  mpi_double_precision, status, ierr)
653  else
654  i = 1
655  do el = 1, nelv
656  do j = 1, lxyz
657  tmp_sp(i) = real(x(j, el), sp)
658  i = i +1
659  end do
660  do j = 1, lxyz
661  tmp_sp(i) = real(y(j, el), sp)
662  i = i +1
663  end do
664  if (gdim .eq. 3) then
665  do j = 1, lxyz
666  tmp_sp(i) = real(z(j, el), sp)
667  i = i +1
668  end do
669  end if
670  end do
671  call mpi_file_write_at_all(fh, byte_offset, tmp_sp, gdim*n, &
672  mpi_real, status, ierr)
673  end if
674 
675 
676  end subroutine fld_file_write_vector_field
677 
679  subroutine fld_file_read(this, data)
680  class(fld_file_t) :: this
681  class(*), target, intent(inout) :: data
682  character(len= 132) :: hdr
683  integer :: ierr, suffix_pos, i, j
684  type(mpi_file) :: fh
685  type(mpi_status) :: status
686  character(len= 1024) :: fname, meta_fname, string, path
687  logical :: meta_file, read_mesh, read_velocity, read_pressure
688  logical :: read_temp
689  character(len=6) :: id_str
690  integer (kind=MPI_OFFSET_KIND) :: mpi_offset, byte_offset
691  integer :: lx, ly, lz, glb_nelv, counter, lxyz
692  integer :: FLD_DATA_SIZE, n_scalars, n
693  real(kind=rp) :: time
694  real(kind=sp) :: temp
695  type(linear_dist_t) :: dist
696  real(kind=sp), parameter :: test_pattern = 6.54321
697  character :: rdcode(10), temp_str(4)
698 
699  select type (data)
700  type is (fld_file_data_t)
701  call filename_chsuffix(this%fname, meta_fname, 'nek5000')
702 
703  inquire(file = trim(meta_fname), exist = meta_file)
704  if (meta_file .and. data%meta_nsamples .eq. 0) then
705  if (pe_rank .eq. 0) then
706  open(unit = 9, file = trim(meta_fname))
707  read(9, fmt = '(A)') string
708  read(string(14:), fmt = '(A)') string
709  string = trim(string)
710  data%fld_series_fname = string(:scan(trim(string), '%')-1)
711  data%fld_series_fname = adjustl(data%fld_series_fname)
712  data%fld_series_fname = trim(data%fld_series_fname)//'0'
713  read(9, fmt = '(A)') string
714  read(string(scan(string, ':')+1:), *) data%meta_start_counter
715  read(9, fmt = '(A)') string
716  read(string(scan(string, ':')+1:), *) data%meta_nsamples
717 
718  close(9)
719  write(*,*) 'Reading meta file for fld series'
720  write(*,*) 'Name: ', trim(data%fld_series_fname)
721  write(*,*) 'Start counter: ', data%meta_start_counter, &
722  'Nsamples: ', data%meta_nsamples
723  end if
724  call mpi_bcast(data%fld_series_fname, 1024, mpi_character, 0, &
725  neko_comm, ierr)
726  call mpi_bcast(data%meta_start_counter, 1, mpi_integer, 0, &
727  neko_comm, ierr)
728  call mpi_bcast(data%meta_nsamples, 1, mpi_integer, 0, &
729  neko_comm, ierr)
730  if (this%counter .eq. 0) this%counter = data%meta_start_counter
731  end if
732 
733  if (meta_file) then
734  write(id_str, '(a,i5.5)') 'f', this%counter
735  path = trim(meta_fname(1:scan(meta_fname, '/', .true. )))
736  fname = trim(path)//trim(data%fld_series_fname)//'.'//id_str
737  if (this%counter .ge. data%meta_nsamples+data%meta_start_counter) then
738  call neko_error('Trying to read more fld files than exist')
739  end if
740  else
741  suffix_pos = filename_suffix_pos(this%fname)
742  write(id_str, '(a,i5.5)') 'f', this%counter
743  fname = trim(this%fname(1:suffix_pos-1))//'.'//id_str
744  end if
745  call mpi_file_open(neko_comm, trim(fname), &
746  mpi_mode_rdonly, mpi_info_null, fh, ierr)
747 
748  if (ierr .ne. 0) call neko_error("Could not read "//trim(fname))
749 
750  call mpi_file_read_all(fh, hdr, 132, mpi_character, status, ierr)
751  ! This read can prorbably be done wihtout the temp variables,
752  ! temp_str, i, j
753 
754  read(hdr, 1) temp_str, fld_data_size, lx, ly, lz, glb_nelv, glb_nelv, &
755  time, counter, i, j, (rdcode(i), i = 1, 10)
756 1 format(4a, 1x, i1, 1x, i2, 1x, i2, 1x, i2, 1x, i10, 1x, i10, &
757  1x, e20.13, 1x, i9, 1x, i6, 1x, i6, 1x, 10a)
758  if (data%nelv .eq. 0) then
759  dist = linear_dist_t(glb_nelv, pe_rank, pe_size, neko_comm)
760  data%nelv = dist%num_local()
761  data%offset_el = dist%start_idx()
762  end if
763  data%lx = lx
764  data%ly = ly
765  data%lz = lz
766  data%glb_nelv = glb_nelv
767  data%t_counter = counter
768  data%time = time
769  lxyz = lx * ly * lz
770  n = lxyz * data%nelv
771 
772  if (lz .eq. 1) then
773  data%gdim = 2
774  else
775  data%gdim = 3
776  end if
777 
778 
779  if (fld_data_size .eq. mpi_double_precision_size) then
780  this%dp_precision = .true.
781  else
782  this%dp_precision = .false.
783  end if
784  if (this%dp_precision) then
785  allocate(tmp_dp(data%gdim*n))
786  else
787  allocate(tmp_sp(data%gdim*n))
788  end if
789 
790 
791  i = 1
792  read_mesh = .false.
793  read_velocity = .false.
794  read_pressure = .false.
795  read_temp = .false.
796  if (rdcode(i) .eq. 'X') then
797  read_mesh = .true.
798  if (data%x%n .ne. n) call data%x%init(n)
799  if (data%y%n .ne. n) call data%y%init(n)
800  if (data%z%n .ne. n) call data%z%init(n)
801  i = i + 1
802  end if
803  if (rdcode(i) .eq. 'U') then
804  read_velocity = .true.
805  if (data%u%n .ne. n) call data%u%init(n)
806  if (data%v%n .ne. n) call data%v%init(n)
807  if (data%w%n .ne. n) call data%w%init(n)
808  i = i + 1
809  end if
810  if (rdcode(i) .eq. 'P') then
811  read_pressure = .true.
812  if (data%p%n .ne. n) call data%p%init(n)
813  i = i + 1
814  end if
815  if (rdcode(i) .eq. 'T') then
816  read_temp = .true.
817  if (data%t%n .ne. n) call data%t%init(n)
818  i = i + 1
819  end if
820  n_scalars = 0
821  if (rdcode(i) .eq. 'S') then
822  i = i + 1
823  read(rdcode(i),*) n_scalars
824  n_scalars = n_scalars*10
825  i = i + 1
826  read(rdcode(i),*) j
827  n_scalars = n_scalars+j
828  i = i + 1
829  if (allocated(data%s)) then
830  if (data%n_scalars .ne. n_scalars) then
831  do j = 1, data%n_scalars
832  call data%s(j)%free()
833  end do
834  deallocate(data%s)
835  data%n_scalars = n_scalars
836  allocate(data%s(n_scalars))
837  do j = 1, data%n_scalars
838  call data%s(j)%init(n)
839  end do
840  end if
841  else
842  data%n_scalars = n_scalars
843  allocate(data%s(data%n_scalars))
844  do j = 1, data%n_scalars
845  call data%s(j)%init(n)
846  end do
847  end if
848  i = i + 1
849  end if
850 
851  mpi_offset = 132 * mpi_character_size
852  call mpi_file_read_at_all(fh, mpi_offset, temp, 1, &
853  mpi_real, status, ierr)
854  if (temp .ne. test_pattern) then
855  call neko_error('Incorrect format for fld file, &
856  &test pattern does not match.')
857  end if
858  mpi_offset = mpi_offset + mpi_real_size
859 
860 
861  if (allocated(data%idx)) then
862  if (size(data%idx) .ne. data%nelv) then
863  deallocate(data%idx)
864  allocate(data%idx(data%nelv))
865  end if
866  else
867  allocate(data%idx(data%nelv))
868  end if
869 
870  byte_offset = mpi_offset + &
871  int(data%offset_el, i8) * int(mpi_integer_size, i8)
872 
873  call mpi_file_read_at_all(fh, byte_offset, data%idx, data%nelv, &
874  mpi_integer, status, ierr)
875 
876  mpi_offset = mpi_offset + &
877  int(data%glb_nelv, i8) * int(mpi_integer_size, i8)
878 
879  if (read_mesh) then
880  byte_offset = mpi_offset + int(data%offset_el, i8) * &
881  (int(data%gdim*lxyz, i8) * &
882  int(fld_data_size, i8))
883  call fld_file_read_vector_field(this, fh, byte_offset, &
884  data%x, data%y, data%z, data)
885  mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
886  (int(data%gdim *lxyz, i8) * &
887  int(fld_data_size, i8))
888  end if
889 
890  if (read_velocity) then
891  byte_offset = mpi_offset + int(data%offset_el, i8) * &
892  (int(data%gdim*lxyz, i8) * &
893  int(fld_data_size, i8))
894  call fld_file_read_vector_field(this, fh, byte_offset, &
895  data%u, data%v, data%w, data)
896  mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
897  (int(data%gdim *lxyz, i8) * &
898  int(fld_data_size, i8))
899  end if
900 
901  if (read_pressure) then
902  byte_offset = mpi_offset + int(data%offset_el, i8) * &
903  (int(lxyz, i8) * &
904  int(fld_data_size, i8))
905  call fld_file_read_field(this, fh, byte_offset, data%p, data)
906  mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
907  (int(lxyz, i8) * &
908  int(fld_data_size, i8))
909  end if
910 
911  if (read_temp) then
912  byte_offset = mpi_offset + int(data%offset_el, i8) * &
913  (int(lxyz, i8) * &
914  int(fld_data_size, i8))
915  call fld_file_read_field(this, fh, byte_offset, data%t, data)
916  mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
917  (int(lxyz, i8) * &
918  int(fld_data_size, i8))
919  end if
920 
921  do i = 1, n_scalars
922  byte_offset = mpi_offset + int(data%offset_el, i8) * &
923  (int(lxyz, i8) * &
924  int(fld_data_size, i8))
925  call fld_file_read_field(this, fh, byte_offset, data%s(i), data)
926  mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
927  (int(lxyz, i8) * &
928  int(fld_data_size, i8))
929  end do
930 
931  this%counter = this%counter + 1
932 
933  if (allocated(tmp_dp)) deallocate(tmp_dp)
934  if (allocated(tmp_sp)) deallocate(tmp_sp)
935  class default
936  call neko_error('Currently we only read into fld_file_data_t, &
937  &please use that data structure instead.')
938  end select
939 
940  end subroutine fld_file_read
941 
942  subroutine fld_file_read_field(this, fh, byte_offset, x, fld_data)
943  class(fld_file_t), intent(inout) :: this
944  type(vector_t), intent(inout) :: x
945  type(fld_file_data_t) :: fld_data
946  integer(kind=MPI_OFFSET_KIND) :: byte_offset
947  type(mpi_file) :: fh
948  type(mpi_status) :: status
949  integer :: n, ierr, lxyz, i
950 
951  n = x%n
952  lxyz = fld_data%lx*fld_data%ly*fld_data%lz
953 
954  if (this%dp_precision) then
955  call mpi_file_read_at_all(fh, byte_offset, tmp_dp, n, &
956  mpi_double_precision, status, ierr)
957  else
958  call mpi_file_read_at_all(fh, byte_offset, tmp_sp, n, &
959  mpi_real, status, ierr)
960  end if
961 
962  if (this%dp_precision) then
963  do i = 1, n
964  x%x(i) = tmp_dp(i)
965  end do
966  else
967  do i = 1, n
968  x%x(i) = tmp_sp(i)
969  end do
970  end if
971 
972 
973  end subroutine fld_file_read_field
974 
975 
976  subroutine fld_file_read_vector_field(this, fh, byte_offset, &
977  x, y, z, fld_data)
978  class(fld_file_t), intent(inout) :: this
979  type(vector_t), intent(inout) :: x, y, z
980  type(fld_file_data_t) :: fld_data
981  integer(kind=MPI_OFFSET_KIND) :: byte_offset
982  type(mpi_file) :: fh
983  type(mpi_status) :: status
984  integer :: n, ierr, lxyz, i, j, e, nd
985 
986  n = x%n
987  nd = n*fld_data%gdim
988  lxyz = fld_data%lx*fld_data%ly*fld_data%lz
989 
990  if (this%dp_precision) then
991  call mpi_file_read_at_all(fh, byte_offset, tmp_dp, nd, &
992  mpi_double_precision, status, ierr)
993  else
994  call mpi_file_read_at_all(fh, byte_offset, tmp_sp, nd, &
995  mpi_real, status, ierr)
996  end if
997 
998 
999  if (this%dp_precision) then
1000  i = 1
1001  do e = 1, fld_data%nelv
1002  do j = 1, lxyz
1003  x%x((e-1)*lxyz+j) = tmp_dp(i)
1004  i = i +1
1005  end do
1006  do j = 1, lxyz
1007  y%x((e-1)*lxyz+j) = tmp_dp(i)
1008  i = i +1
1009  end do
1010  if (fld_data%gdim .eq. 3) then
1011  do j = 1, lxyz
1012  z%x((e-1)*lxyz+j) = tmp_dp(i)
1013  i = i +1
1014  end do
1015  end if
1016  end do
1017  else
1018  i = 1
1019  do e = 1, fld_data%nelv
1020  do j = 1, lxyz
1021  x%x((e-1)*lxyz+j) = tmp_sp(i)
1022  i = i +1
1023  end do
1024  do j = 1, lxyz
1025  y%x((e-1)*lxyz+j) = tmp_sp(i)
1026  i = i +1
1027  end do
1028  if (fld_data%gdim .eq. 3) then
1029  do j = 1, lxyz
1030  z%x((e-1)*lxyz+j) = tmp_sp(i)
1031  i = i +1
1032  end do
1033  end if
1034  end do
1035  end if
1036 
1037  end subroutine fld_file_read_vector_field
1038 
1039  subroutine fld_file_set_precision(this, precision)
1040  class(fld_file_t) :: this
1041  integer, intent(in) :: precision
1042 
1043  if (precision .eq. dp) then
1044  this%dp_precision = .true.
1045  else if (precision .eq. sp) then
1046  this%dp_precision = .false.
1047  else
1048  call neko_error('Invalid precision')
1049  end if
1050 
1051  end subroutine fld_file_set_precision
1052 
1053 
1054 end module fld_file
double real
Definition: device_config.h:12
Definition: comm.F90:1
integer pe_rank
MPI rank.
Definition: comm.F90:28
type(mpi_comm) neko_comm
MPI communicator.
Definition: comm.F90:16
integer pe_size
MPI size of communicator.
Definition: comm.F90:31
Defines practical data distributions.
Definition: datadist.f90:34
Defines a mapping of the degrees of freedom.
Definition: dofmap.f90:35
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
real(kind=dp), dimension(:), allocatable, private tmp_dp
Definition: fld_file.f90:60
subroutine fld_file_set_precision(this, precision)
Definition: fld_file.f90:1040
subroutine fld_file_read_vector_field(this, fh, byte_offset, x, y, z, fld_data)
Definition: fld_file.f90:978
subroutine fld_file_read(this, data)
Load a field from a NEKTON fld file.
Definition: fld_file.f90:680
subroutine fld_file_write_vector_field(this, fh, byte_offset, x, y, z, n, gdim, lxyz, nelv)
Definition: fld_file.f90:625
subroutine fld_file_write_metadata_vector(this, fh, byte_offset, x, y, z, gdim, lxyz, nelv)
Definition: fld_file.f90:540
real(kind=sp), dimension(:), allocatable, private tmp_sp
Definition: fld_file.f90:61
subroutine fld_file_read_field(this, fh, byte_offset, x, fld_data)
Definition: fld_file.f90:943
subroutine fld_file_write_field(this, fh, byte_offset, p, n)
Definition: fld_file.f90:598
subroutine fld_file_write_metadata_scalar(this, fh, byte_offset, x, lxyz, nelv)
Definition: fld_file.f90:573
subroutine fld_file_write(this, data, t)
Write fields to a NEKTON fld file.
Definition: fld_file.f90:78
Definition: math.f90:60
real(kind=rp) function, public vlmin(vec, n)
minimun value of a vector of length n
Definition: math.f90:463
real(kind=rp) function, public vlmax(vec, n)
maximum value of a vector of length n
Definition: math.f90:452
Defines a mean flow field.
Definition: mean_flow.f90:34
Defines a mean squared flow field.
Defines a mesh.
Definition: mesh.f90:34
MPI derived types.
Definition: mpi_types.f90:34
integer, public mpi_double_precision_size
Size of MPI type double precision.
Definition: mpi_types.f90:61
integer, public mpi_character_size
Size of MPI type character.
Definition: mpi_types.f90:62
integer, public mpi_real_size
Size of MPI type real.
Definition: mpi_types.f90:60
integer, public mpi_integer_size
Size of MPI type integer.
Definition: mpi_types.f90:63
integer, parameter, public i8
Definition: num_types.f90:7
integer, parameter, public dp
Definition: num_types.f90:9
integer, parameter, public sp
Definition: num_types.f90:8
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Defines a function space.
Definition: space.f90:34
Defines structs that are used... Dont know if we should keep it though.
Definition: structs.f90:2
Utilities.
Definition: utils.f90:35
subroutine, public filename_chsuffix(fname, new_fname, new_suffix)
Change a filename's suffix.
Definition: utils.f90:77
pure integer function, public filename_tslash_pos(fname)
Find position (in the string) of a filename's trailing slash.
Definition: utils.f90:63
pure integer function, public filename_suffix_pos(fname)
Find position (in the string) of a filename's suffix.
Definition: utils.f90:56
Defines a vector.
Definition: vector.f90:34
Load-balanced linear distribution .
Definition: datadist.f90:50
field_list_t, To be able to group fields together
Definition: field_list.f90:13
Interface for NEKTON fld files.
Definition: fld_file.f90:64
A generic file handler.
The function space for the SEM solution fields.
Definition: space.f90:62
Pointer to array.
Definition: structs.f90:14