Neko  0.9.0
A portable framework for high-order spectral element flow simulations
rea_file.f90
Go to the documentation of this file.
1 ! Copyright (c) 2019-2022, 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 rea_file
36  use generic_file
37  use num_types
38  use utils
39  use mesh
40  use point
41  use map
42  use rea
43  use re2_file, only: re2_file_t
44  use map_file
45  use comm
46  use datadist
47  use htable
48  use logger
49  implicit none
50  private
51 
52  ! Defines the conventions for conversion of rea labels to labeled zones.
53  integer, public :: neko_w_bc_label = -1
54  integer, public :: neko_v_bc_label = -1
55  integer, public :: neko_o_bc_label = -1
56  integer, public :: neko_sym_bc_label = -1
57  integer, public :: neko_on_bc_label = -1
58  integer, public :: neko_shl_bc_label = -1
59 
61  type, public, extends(generic_file_t) :: rea_file_t
62  contains
63  procedure :: read => rea_file_read
64  procedure :: write => rea_file_write
65  end type rea_file_t
66 
67 contains
68 
70  subroutine rea_file_read(this, data)
71  class(rea_file_t) :: this
72  class(*), target, intent(inout) :: data
73  type(mesh_t), pointer :: msh
74  real(kind=dp), pointer :: params(:)
75  character(len=3), pointer :: cbc(:,:)
76  integer, allocatable :: curve_type(:,:)
77  logical, allocatable :: curve_element(:)
78  character(len=1) :: chtemp
79  integer :: ndim, nparam, nskip, nlogic, ncurve
80  integer :: nelgs, nelgv, i, j, ierr, l
81  integer :: el_idx, pt_idx
82  logical :: read_param, read_bcs, read_map
83  real(kind=dp) :: xc(8), yc(8), zc(8), curve(5)
84  real(kind=dp), allocatable :: bc_data(:,:,:)
85  real(kind=dp), allocatable :: curve_data(:,:,:)
86  type(point_t) :: p(8)
87  type(re2_file_t) :: re2_file
88  type(map_file_t) :: map_file
89  character(len=1024) :: re2_fname, map_fname
90  integer :: start_el, end_el, nel, edge
91  type(linear_dist_t) :: dist
92  type(map_t) :: nm
93  type(htable_pt_t) :: htp
94  integer :: sym_facet, pids(4), p_el_idx, p_facet
95  integer :: off
96  integer, parameter, dimension(6) :: facet_map = (/3, 2, 4, 1, 5, 6/)
97  logical :: curve_skip = .false.
98  character(len=LOG_SIZE) :: log_buf
99  ! ---- Offsets for conversion from internal zones to labeled zones
100  integer :: user_labeled_zones(NEKO_MSH_MAX_ZLBLS)
101  integer :: labeled_zone_offsets(NEKO_MSH_MAX_ZLBLS)
102  integer :: total_labeled_zone_offset, current_internal_zone
103  total_labeled_zone_offset = 0
104  labeled_zone_offsets = 0
105  user_labeled_zones = 0
106  current_internal_zone = 1
107  ! ----
108 
109  select type(data)
110  type is (rea_t)
111  call rea_free(data)
112  msh => data%msh
113  params => data%params
114  cbc => data%cbc
115  read_param = .true.
116  read_bcs = .true.
117  type is (mesh_t)
118  msh => data
119  read_param = .false.
120  read_bcs = .false.
121  class default
122  call neko_error('Invalid output data')
123  end select
124 
125  if (read_param .and. read_bcs .and. pe_size .gt. 1) then
126  call neko_error('Reading NEKTON session data only implemented in serial')
127  end if
128 
129 
130  open(unit=9,file=trim(this%fname), status='old', iostat=ierr)
131  call neko_log%message('Reading NEKTON file ' // this%fname)
132 
133  read(9, *)
134  read(9, *)
135  read(9, *) ndim
136  read(9, *) nparam
137 
138  if (.not. read_param) then
139  ! Skip parameters
140  do i = 1, nparam
141  read(9, *)
142  end do
143  else
144  allocate(params(nparam))
145  do i = 1, nparam
146  read(9, *) params(i)
147  end do
148  end if
149 
150  ! Skip passive scalars
151  read(9, *) nskip
152  do i = 1, nskip
153  read(9, *)
154  end do
155 
156  ! Skip logic switches
157  read(9, *) nlogic
158  do i = 1, nlogic
159  read(9, *)
160  end do
161 
162  ! Read mesh info
163  read(9, *)
164  read(9, *)
165  read(9, *) nelgs,ndim, nelgv
166  if (nelgs .lt. 0) then
167  re2_fname = trim(this%fname(1:scan(trim(this%fname), &
168  '.', back=.true.)))//'re2'
169  call re2_file%init(re2_fname)
170  call re2_file%read(msh)
171  else
172  write(log_buf,1) ndim, nelgv
173 1 format('gdim = ', i1, ', nelements =', i7)
174  call neko_log%message(log_buf)
175 
176  call filename_chsuffix(this%fname, map_fname, 'map')
177  inquire(file=map_fname, exist=read_map)
178  if (read_map) then
179  call map_init(nm, nelgv, 2**ndim)
180  call map_file%init(map_fname)
181  call map_file%read(nm)
182  else
183  call neko_log%warning('No NEKTON map file found')
184  end if
185 
186  ! Use a load-balanced linear distribution
187  dist = linear_dist_t(nelgv, pe_rank, pe_size, neko_comm)
188  nel = dist%num_local()
189  start_el = dist%start_idx() + 1
190  end_el = dist%end_idx() + 1
191 
192  call msh%init(ndim, dist)
193 
194  call htp%init(2**ndim * nel, ndim)
195 
196  el_idx = 1
197  pt_idx = 0
198  do i = 1, nelgv
199  read(9, *)
200  if (ndim .eq. 2) then
201  read(9, *) (xc(j),j=1,4)
202  read(9, *) (yc(j),j=1,4)
203  if (i .ge. start_el .and. i .le. end_el) then
204  do j = 1, 4
205  p(j) = point_t(real(xc(j),dp), real(yc(j),dp),real(0d0,dp))
206  call rea_file_add_point(htp, p(j), pt_idx)
207  end do
208  ! swap vertices to keep symmetric vertex numbering in neko
209  call msh%add_element(el_idx, p(1), p(2), p(4), p(3))
210  end if
211  else if (ndim .eq. 3) then
212  read(9, *) (xc(j),j=1,4)
213  read(9, *) (yc(j),j=1,4)
214  read(9, *) (zc(j),j=1,4)
215  read(9, *) (xc(j),j=5,8)
216  read(9, *) (yc(j),j=5,8)
217  read(9, *) (zc(j),j=5,8)
218  if (i .ge. start_el .and. i .le. end_el) then
219  do j = 1, 8
220  p(j) = point_t(real(xc(j),dp), real(yc(j),dp), real(zc(j),dp))
221  call rea_file_add_point(htp, p(j), pt_idx)
222  end do
223  ! swap vertices to keep symmetric vertex numbering in neko
224  call msh%add_element(el_idx, &
225  p(1), p(2), p(4), p(3), p(5), p(6), p(8), p(7))
226  end if
227  end if
228  if (i .ge. start_el .and. i .le. end_el) then
229  el_idx = el_idx + 1
230  end if
231  end do
232 
233  call htp%free()
234 
235  read(9, *)
236  read(9, *) ncurve
237  allocate(curve_data(5,8,nelgv))
238  allocate(curve_element(nelgv))
239  allocate(curve_type(8,nelgv))
240  do i = 1, nelgv
241  curve_element(i) = .false.
242  do j = 1, 8
243  curve_type(j,i) = 0
244  do l = 1, 5
245  curve_data(l,j,i) = 0d0
246  end do
247  end do
248  end do
249  do i = 1, ncurve
250  read(9, *) edge, el_idx, (curve(j),j=1,5), chtemp
251  do j = 1, 5
252  curve_data(j,edge,el_idx) = curve(j)
253  end do
254  curve_element(el_idx) = .true.
255  select case(trim(chtemp))
256  case ('s')
257  curve_type(edge,el_idx) = 1
258  curve_skip = .true.
259  case ('e')
260  curve_type(edge,el_idx) = 2
261  curve_skip = .true.
262  case ('C')
263  curve_type(edge,el_idx) = 3
264  case ('m')
265  curve_type(edge,el_idx) = 4
266  end select
267  end do
268  if (curve_skip) then
269  call neko_log%warning('Curve type: s, e are not supported, treating mesh as non-curved.')
270  else
271  do el_idx = 1, nelgv
272  if (curve_element(el_idx)) then
273  call msh%mark_curve_element(el_idx, &
274  curve_data(1,1,el_idx), curve_type(1,el_idx))
275  end if
276  end do
277  end if
278  deallocate(curve_data)
279  deallocate(curve_element)
280  deallocate(curve_type)
281 
282  ! Read fluid boundary conditions
283  read(9,*)
284  read(9,*)
285  if (.not. read_bcs) then ! Mark zones in the mesh
286  call neko_log%message("Reading boundary conditions", neko_log_debug)
287  allocate(cbc(6,nelgv))
288  allocate(bc_data(6,2*ndim,nelgv))
289  off = 0
290  !Fix for different horrible .rea periodic bc formats.
291  if (nelgv .lt. 1000) off = 1
292  do i = 1, nelgv
293  if (i .ge. start_el .and. i .le. end_el) then
294  el_idx = i - start_el + 1
295  do j = 1, 2*ndim
296  read(9, *) cbc(j, i), (bc_data(l,j,i),l=1,6)
297  sym_facet = facet_map(j)
298 
299  select case(trim(cbc(j,i)))
300  case ('W')
301  if (neko_w_bc_label .eq. -1) then
302  neko_w_bc_label = current_internal_zone
303  current_internal_zone = current_internal_zone + 1
304  end if
305 
306  call rea_file_mark_labeled_bc(msh, el_idx, sym_facet, &
307  cbc(j,i), &
308  neko_w_bc_label, total_labeled_zone_offset, &
309  labeled_zone_offsets(neko_w_bc_label) .eq. 0)
310 
311  labeled_zone_offsets(neko_w_bc_label) = 1
312 
313  case ('v', 'V')
314  if (neko_v_bc_label .eq. -1) then
315  neko_v_bc_label = current_internal_zone
316  current_internal_zone = current_internal_zone + 1
317  end if
318 
319  call rea_file_mark_labeled_bc(msh, el_idx, sym_facet, &
320  cbc(j,i), &
321  neko_v_bc_label, total_labeled_zone_offset, &
322  labeled_zone_offsets(neko_v_bc_label) .eq. 0)
323 
324  labeled_zone_offsets(neko_v_bc_label) = 1
325  case ('O', 'o')
326  if (neko_o_bc_label .eq. -1) then
327  neko_o_bc_label = current_internal_zone
328  current_internal_zone = current_internal_zone + 1
329  end if
330 
331  call rea_file_mark_labeled_bc(msh, el_idx, sym_facet, &
332  cbc(j,i), &
333  neko_o_bc_label, total_labeled_zone_offset, &
334  labeled_zone_offsets(neko_o_bc_label) .eq. 0)
335 
336  labeled_zone_offsets(neko_o_bc_label) = 1
337  case ('SYM')
338  if (neko_sym_bc_label .eq. -1) then
339  neko_sym_bc_label = current_internal_zone
340  current_internal_zone = current_internal_zone + 1
341  end if
342 
343  call rea_file_mark_labeled_bc(msh, el_idx, sym_facet, &
344  cbc(j,i), &
345  neko_sym_bc_label, total_labeled_zone_offset, &
346  labeled_zone_offsets(neko_sym_bc_label) .eq. 0)
347 
348  labeled_zone_offsets(neko_sym_bc_label) = 1
349  case ('ON', 'on')
350  if (neko_on_bc_label .eq. -1) then
351  neko_on_bc_label = current_internal_zone
352  current_internal_zone = current_internal_zone + 1
353  end if
354 
355  call rea_file_mark_labeled_bc(msh, el_idx, sym_facet, &
356  cbc(j,i), &
357  neko_on_bc_label, total_labeled_zone_offset, &
358  labeled_zone_offsets(neko_on_bc_label) .eq. 0)
359 
360  labeled_zone_offsets(neko_on_bc_label) = 1
361  case ('s', 'sl', 'sh', 'shl', 'S', 'SL', 'SH', 'SHL')
362  if (neko_shl_bc_label .eq. -1) then
363  neko_shl_bc_label = current_internal_zone
364  current_internal_zone = current_internal_zone + 1
365  end if
366 
367  call rea_file_mark_labeled_bc(msh, el_idx, sym_facet, &
368  cbc(j,i), &
369  neko_shl_bc_label, total_labeled_zone_offset, &
370  labeled_zone_offsets(neko_shl_bc_label) .eq. 0)
371 
372  labeled_zone_offsets(neko_shl_bc_label) = 1
373  case ('P')
374  p_el_idx = int(bc_data(2+off,j,i))
375  p_facet = facet_map(int(bc_data(3+off,j,i)))
376  call msh%get_facet_ids(sym_facet, el_idx, pids)
377  call msh%mark_periodic_facet(sym_facet, el_idx, &
378  p_facet, p_el_idx, pids)
379  end select
380  end do
381  end if
382  end do
383 
384  do i = 1, nelgv
385  if (i .ge. start_el .and. i .le. end_el) then
386  el_idx = i - start_el + 1
387  do j = 1, 2*ndim
388  sym_facet = facet_map(j)
389  select case(trim(cbc(j,i)))
390  case ('P')
391  p_el_idx = int(bc_data(2+off,j,i))
392  p_facet = facet_map(int(bc_data(3+off,j,i)))
393  call msh%create_periodic_ids(sym_facet, el_idx, &
394  p_facet, p_el_idx)
395  end select
396  end do
397  end if
398  end do
399  do i = 1, nelgv
400  if (i .ge. start_el .and. i .le. end_el) then
401  el_idx = i - start_el + 1
402  do j = 1, 2*ndim
403  sym_facet = facet_map(j)
404  select case(trim(cbc(j,i)))
405  case ('P')
406  p_el_idx = int(bc_data(2+off,j,i))
407  p_facet = facet_map(int(bc_data(3+off,j,i)))
408  call msh%create_periodic_ids(sym_facet, el_idx, &
409  p_facet, p_el_idx)
410  end select
411  end do
412  end if
413  end do
414  do i = 1, nelgv
415  if (i .ge. start_el .and. i .le. end_el) then
416  el_idx = i - start_el + 1
417  do j = 1, 2*ndim
418  sym_facet = facet_map(j)
419  select case(trim(cbc(j,i)))
420  case ('P')
421  p_el_idx = int(bc_data(2+off,j,i))
422  p_facet = facet_map(int(bc_data(3+off,j,i)))
423  call msh%create_periodic_ids(sym_facet, el_idx, &
424  p_facet, p_el_idx)
425  end select
426  end do
427  end if
428  end do
429  deallocate(cbc)
430  deallocate(bc_data)
431  else ! Store bcs in a NEKTON session structure
432  allocate(cbc(6,nelgv))
433  do i = 1, nelgv
434  do j = 1, 2*ndim
435  read(9,'(a1, a3)') chtemp, cbc(j, i)
436  end do
437  end do
438  end if
439 
440  call msh%finalize()
441 
442  call neko_log%message('Done')
443  close(9)
444  endif
445 
446  end subroutine rea_file_read
447 
448  subroutine rea_file_write(this, data, t)
449  class(rea_file_t), intent(inout) :: this
450  class(*), target, intent(in) :: data
451  real(kind=rp), intent(in), optional :: t
452  end subroutine rea_file_write
453 
454  subroutine rea_file_add_point(htp, p, idx)
455  type(htable_pt_t), intent(inout) :: htp
456  type(point_t), intent(inout) :: p
457  integer, intent(inout) :: idx
458  integer :: tmp
459 
460  if (htp%get(p, tmp) .gt. 0) then
461  idx = idx + 1
462  call htp%set(p, idx)
463  call p%set_id(idx)
464  else
465  call p%set_id(tmp)
466  end if
467 
468  end subroutine rea_file_add_point
469 
480  subroutine rea_file_mark_labeled_bc(msh, el_idx, facet, type, label, offset, print_info)
481  type(mesh_t), intent(inout) :: msh
482  integer, intent(in) :: el_idx
483  integer, intent(in) :: facet
484  character(len=*), intent(in) :: type
485  integer, intent(in) :: label
486  integer, intent(in) :: offset
487  logical, intent(in) :: print_info
488 
489  integer :: mark_label
490  character(len=LOG_SIZE) :: log_buf
491 
492  mark_label = offset + label
493 
494  if (mark_label .lt. 1 .or. mark_label .gt. neko_msh_max_zlbls) then
495  call neko_error("You have reached the maximum amount of allowed labeled&
496 & zones (max allowed: 20). This happened when converting re2 internal labels&
497 & like e.g. 'w', 'V' or 'o' to labeled zones. Please reduce the number of&
498 & labeled zones that you have defined or make sure that they are labeled&
499 & from [1,...,20].")
500  end if
501 
502  if (print_info) then
503  write (log_buf, "(A3,A,I2)") trim(type), " => Labeled index ", mark_label
504  call neko_log%message(log_buf)
505  end if
506  call msh%mark_labeled_facet(facet, el_idx, mark_label)
507 
508  end subroutine rea_file_mark_labeled_bc
509 
510 end module rea_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 a domain as a subset of facets in a mesh.
Definition: curve.f90:2
Defines practical data distributions.
Definition: datadist.f90:34
Module for file I/O operations.
Definition: file.f90:34
Implements a hash table ADT.
Definition: htable.f90:36
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
NEKTON map file.
Definition: map_file.f90:35
NEKTON map.
Definition: map.f90:3
Defines a mesh.
Definition: mesh.f90:34
integer, parameter, public neko_msh_max_zlbls
Max num. zone labels.
Definition: mesh.f90:56
integer, parameter, public dp
Definition: num_types.f90:9
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Implements a point.
Definition: point.f90:35
NEKTON mesh data in re2 format.
Definition: re2_file.f90:35
NEKTON session data reader.
Definition: rea_file.f90:35
integer, public neko_shl_bc_label
Definition: rea_file.f90:58
integer, public neko_o_bc_label
Definition: rea_file.f90:55
subroutine rea_file_read(this, data)
Load NEKTON session data from an ascii file.
Definition: rea_file.f90:71
integer, public neko_w_bc_label
Definition: rea_file.f90:53
integer, public neko_v_bc_label
Definition: rea_file.f90:54
integer, public neko_on_bc_label
Definition: rea_file.f90:57
subroutine rea_file_add_point(htp, p, idx)
Definition: rea_file.f90:455
integer, public neko_sym_bc_label
Definition: rea_file.f90:56
subroutine rea_file_write(this, data, t)
Definition: rea_file.f90:449
subroutine rea_file_mark_labeled_bc(msh, el_idx, facet, type, label, offset, print_info)
Mark a labeled zone based on an offset.
Definition: rea_file.f90:481
NEKTON session data.
Definition: rea.f90:4
subroutine, public rea_free(r)
Free a NEKTON session data.
Definition: rea.f90:26
Utilities.
Definition: utils.f90:35
subroutine, public filename_chsuffix(fname, new_fname, new_suffix)
Change a filename's suffix.
Definition: utils.f90:77
Load-balanced linear distribution .
Definition: datadist.f90:50
A generic file handler.
Point based hash table.
Definition: htable.f90:112
NEKTON vertex mapping.
Definition: map.f90:8
Interface for NEKTON map files.
Definition: map_file.f90:44
A point in with coordinates .
Definition: point.f90:43
Interface for NEKTON re2 files.
Definition: re2_file.f90:62
NEKTON session data struct.
Definition: rea.f90:13
Interface for NEKTON ascii files.
Definition: rea_file.f90:61