Neko  0.9.0
A portable framework for high-order spectral element flow simulations
re2_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 re2_file
36  use generic_file, only: generic_file_t
37  use num_types, only: rp
38  use utils
39  use mesh, only: mesh_t
40  use point, only: point_t
41  use comm
42  use mpi_f08
43  use neko_mpi_types
44  use datadist
45  use re2
46  use map
47  use map_file
48  use htable
49  use logger, only: neko_log, log_size
50  implicit none
51  private
52 
53  ! Defines the conventions for conversion of re2 labels to labeled zones.
54  integer, public :: neko_w_bc_label = -1
55  integer, public :: neko_v_bc_label = -1
56  integer, public :: neko_o_bc_label = -1
57  integer, public :: neko_sym_bc_label = -1
58  integer, public :: neko_on_bc_label = -1
59  integer, public :: neko_shl_bc_label = -1
60 
62  type, public, extends(generic_file_t) :: re2_file_t
63  contains
64  procedure :: read => re2_file_read
65  procedure :: write => re2_file_write
66  end type re2_file_t
67 
68 contains
69 
71  subroutine re2_file_read(this, data)
72  class(re2_file_t) :: this
73  class(*), target, intent(inout) :: data
74  type(mesh_t), pointer :: msh
75  character(len=5) :: hdr_ver
76  character(len=54) :: hdr_str
77  character(len=80) :: hdr_full
78  integer :: nel, ndim, nelv, ierr, nBCre2
79  type(mpi_status) :: status
80  type(mpi_file) :: fh
81  integer (kind=MPI_OFFSET_KIND) :: mpi_offset
82  real(kind=sp) :: test
83  real(kind=dp) :: t2
84  integer :: ncurv, nbcs
85  type(linear_dist_t) :: dist
86  type(map_t) :: nm
87  type(map_file_t) :: map_file
88  character(len=1024) :: map_fname
89  logical :: read_map
90  integer :: re2_data_xy_size
91  integer :: re2_data_xyz_size
92  integer :: re2_data_cv_size
93  integer :: re2_data_bc_size
94  logical :: v2_format
95  character(len=LOG_SIZE) :: log_buf
96 
97  call this%check_exists()
98 
99  select type(data)
100  type is (mesh_t)
101  msh => data
102  class default
103  call neko_error('Invalid output data')
104  end select
105 
106  v2_format = .false.
107  open(unit=9,file=trim(this%fname), status='old', iostat=ierr)
108  call neko_log%message('Reading binary NEKTON file ' // this%fname)
109 
110  read(9,'(a80)') hdr_full
111  read(hdr_full, '(a5)') hdr_ver
112 
113  if (hdr_ver .eq. '#v004') then
114  read(hdr_full, '(a5,i16,i3,i16,i4,a36)') hdr_ver, nel, ndim, nelv, nbcre2, hdr_str
115  v2_format = .true.
116  else if (hdr_ver .eq. '#v002' .or. hdr_ver .eq. '#v003') then
117  read(hdr_full, '(a5,i9,i3,i9,a54)') hdr_ver, nel, ndim, nelv, hdr_str
118  v2_format = .true.
119  else if (hdr_ver .eq. '#v001') then
120  read(hdr_full, '(a5,i9,i3,i9,a54)') hdr_ver, nel, ndim, nelv, hdr_str
121  end if
122 
123  if (v2_format) then
124  call mpi_type_size(mpi_re2v2_data_xy, re2_data_xy_size, ierr)
125  call mpi_type_size(mpi_re2v2_data_xyz, re2_data_xyz_size, ierr)
126  call mpi_type_size(mpi_re2v2_data_cv, re2_data_cv_size, ierr)
127  call mpi_type_size(mpi_re2v2_data_bc, re2_data_bc_size, ierr)
128  else
129  call mpi_type_size(mpi_re2v1_data_xy, re2_data_xy_size, ierr)
130  call mpi_type_size(mpi_re2v1_data_xyz, re2_data_xyz_size, ierr)
131  call mpi_type_size(mpi_re2v1_data_cv, re2_data_cv_size, ierr)
132  call mpi_type_size(mpi_re2v1_data_bc, re2_data_bc_size, ierr)
133  end if
134 
135  write(log_buf,1) ndim, nelv
136 1 format('gdim = ', i1, ', nelements =', i7)
137  call neko_log%message(log_buf)
138  close(9)
139 
140  call filename_chsuffix(this%fname, map_fname,'map')
141 
142  inquire(file=map_fname, exist=read_map)
143  if (read_map) then
144  call map_init(nm, nelv, 2**ndim)
145  call map_file%init(map_fname)
146  call map_file%read(nm)
147  else
148  call neko_log%warning('No NEKTON map file found')
149  end if
150 
151  call mpi_file_open(neko_comm, trim(this%fname), &
152  mpi_mode_rdonly, mpi_info_null, fh, ierr)
153 
154  if (ierr .ne. 0) then
155  call neko_log%error("Can't open binary NEKTON file ")
156  end if
157  dist = linear_dist_t(nelv, pe_rank, pe_size, neko_comm)
158 
159 
160  call msh%init(ndim, dist)
161 
162  ! Set offset (header)
163  mpi_offset = re2_hdr_size * mpi_character_size
164 
165  call mpi_file_read_at_all(fh, mpi_offset, test, 1, mpi_real, status, ierr)
166  mpi_offset = mpi_offset + mpi_real_size
167 
168  if (abs(re2_endian_test - test) .gt. 1e-4) then
169  call neko_error('Invalid endian of re2 file, byte swap not implemented yet')
170  end if
171 
172  call re2_file_read_points(msh, ndim, nel, dist, fh, &
173  mpi_offset, re2_data_xy_size, re2_data_xyz_size, v2_format)
174 
175 
176  ! Set offset to start of curved side data
178  if (ndim .eq. 2) then
179  mpi_offset = mpi_offset + int(dist%num_global(),i8) * int(re2_data_xy_size,i8)
180  else
181  mpi_offset = mpi_offset + int(dist%num_global(),i8) * int(re2_data_xyz_size,i8)
182  end if
183 
186  if (v2_format) then
187  call mpi_file_read_at_all(fh, mpi_offset, t2, 1, mpi_double_precision, status, ierr)
188  ncurv = int(t2)
189  mpi_offset = mpi_offset + mpi_double_precision_size
190  call re2_file_read_curve(msh, ncurv, dist, fh, mpi_offset, v2_format)
191  mpi_offset = mpi_offset + int(ncurv,i8) * int(re2_data_cv_size,i8)
192  call mpi_file_read_at_all(fh, mpi_offset, t2, 1, mpi_double_precision, status, ierr)
193  nbcs = int(t2)
194  mpi_offset = mpi_offset + mpi_double_precision_size
195 
196  call re2_file_read_bcs(msh, nbcs, dist, fh, mpi_offset, v2_format)
197  else
198  call mpi_file_read_at_all(fh, mpi_offset, ncurv, 1, mpi_integer, status, ierr)
199  mpi_offset = mpi_offset + mpi_integer_size
200  call re2_file_read_curve(msh, ncurv, dist, fh, mpi_offset, v2_format)
201  mpi_offset = mpi_offset + int(ncurv,i8) * int(re2_data_cv_size,i8)
202  call mpi_file_read_at_all(fh, mpi_offset, nbcs, 1, mpi_integer, status, ierr)
203  mpi_offset = mpi_offset + mpi_integer_size
204 
205  call re2_file_read_bcs(msh, nbcs, dist, fh, mpi_offset, v2_format)
206 
207  end if
208 
209  call mpi_file_close(fh, ierr)
210  call msh%finalize()
211 
212 
213  call neko_log%message('Done')
214 
215 
216  end subroutine re2_file_read
217 
218  subroutine re2_file_write(this, data, t)
219  class(re2_file_t), intent(inout) :: this
220  class(*), target, intent(in) :: data
221  real(kind=rp), intent(in), optional :: t
222  type(re2v1_xy_t), allocatable :: re2_data_xy(:)
223  type(re2v1_xyz_t), allocatable :: re2_data_xyz(:)
224  type(mesh_t), pointer :: msh
225  character(len=5), parameter :: RE2_HDR_VER = '#v001'
226  character(len=54), parameter :: RE2_HDR_STR = 'RE2 exported by NEKO'
227  integer :: i, j, ierr, nelgv
228  type(mpi_status) :: status
229  type(mpi_file) :: fh
230  integer (kind=MPI_OFFSET_KIND) :: mpi_offset
231  integer :: element_offset
232  integer :: re2_data_xy_size
233  integer :: re2_data_xyz_size
234 
235  select type(data)
236  type is (mesh_t)
237  msh => data
238  class default
239  call neko_error('Invalid output data')
240  end select
241 
242  call mpi_type_size(mpi_re2v1_data_xy, re2_data_xy_size, ierr)
243  call mpi_type_size(mpi_re2v1_data_xyz, re2_data_xyz_size, ierr)
244  call mpi_reduce(msh%nelv, nelgv, 1, &
245  mpi_integer, mpi_sum, 0, neko_comm, ierr)
246  element_offset = 0
247  call mpi_exscan(msh%nelv, element_offset, 1, &
248  mpi_integer, mpi_sum, neko_comm, ierr)
249 
250  call neko_log%message('Writing data as a binary NEKTON file ' // this%fname)
251 
252  if (pe_rank .eq. 0) then
253  open(unit=9,file=trim(this%fname), status='new', iostat=ierr)
254  write(9, '(a5,i9,i3,i9,a54)') re2_hdr_ver, nelgv, msh%gdim,&
255  nelgv, re2_hdr_str
256  close(9)
257  end if
258 
259  call mpi_barrier(neko_comm, ierr)
260  call mpi_file_open(neko_comm, trim(this%fname), &
261  mpi_mode_wronly + mpi_mode_create, mpi_info_null, fh, ierr)
262  mpi_offset = re2_hdr_size * mpi_character_size
263 
264  call mpi_file_write_at(fh, mpi_offset, re2_endian_test, 1, &
265  mpi_real, status, ierr)
266  mpi_offset = mpi_offset + mpi_real_size
267 
268  if (msh%gdim .eq. 2) then
269  allocate(re2_data_xy(msh%nelv))
270  do i = 1, msh%nelv
271  re2_data_xy(i)%rgroup = 1.0 ! Not used
272  do j = 1, 4
273  re2_data_xy(i)%x(j) = real(msh%elements(i)%e%pts(j)%p%x(1))
274  re2_data_xy(i)%y(j) = real(msh%elements(i)%e%pts(j)%p%x(2))
275  end do
276  end do
277  mpi_offset = mpi_offset + element_offset * re2_data_xy_size
278  call mpi_file_write_at(fh, mpi_offset, &
279  re2_data_xy, msh%nelv, mpi_re2v1_data_xy, status, ierr)
280 
281  deallocate(re2_data_xy)
282 
283  else if (msh%gdim .eq. 3) then
284  allocate(re2_data_xyz(msh%nelv))
285  do i = 1, msh%nelv
286  re2_data_xyz(i)%rgroup = 1.0 ! Not used
287  do j = 1, 8
288  re2_data_xyz(i)%x(j) = real(msh%elements(i)%e%pts(j)%p%x(1))
289  re2_data_xyz(i)%y(j) = real(msh%elements(i)%e%pts(j)%p%x(2))
290  re2_data_xyz(i)%z(j) = real(msh%elements(i)%e%pts(j)%p%x(3))
291  end do
292  end do
293  mpi_offset = mpi_offset + element_offset * re2_data_xyz_size
294  call mpi_file_write_at(fh, mpi_offset, &
295  re2_data_xyz, msh%nelv, mpi_re2v1_data_xyz, status, ierr)
296 
297  deallocate(re2_data_xyz)
298  else
299  call neko_error("Invalid dimension of mesh")
300  end if
301 
302  call mpi_file_close(fh, ierr)
303  call neko_log%message('Done')
304 
306 
307  end subroutine re2_file_write
308 
309  subroutine re2_file_read_points(msh, ndim, nel, dist, fh, &
310  mpi_offset, re2_data_xy_size, re2_data_xyz_size, v2_format)
311  type(mesh_t), intent(inout) :: msh
312  integer (kind=MPI_OFFSET_KIND) :: mpi_offset
313  integer, intent(inout) :: ndim
314  integer, intent(inout) :: nel
315  type(mpi_file), intent(inout) :: fh
316  integer, intent(in) :: re2_data_xy_size
317  integer, intent(in) :: re2_data_xyz_size
318  logical, intent(in) :: v2_format
319  type(linear_dist_t) :: dist
320  integer :: element_offset
321  type(re2v1_xy_t), allocatable :: re2v1_data_xy(:)
322  type(re2v1_xyz_t), allocatable :: re2v1_data_xyz(:)
323  type(re2v2_xy_t), allocatable :: re2v2_data_xy(:)
324  type(re2v2_xyz_t), allocatable :: re2v2_data_xyz(:)
325  type(mpi_status) :: status
326  type(htable_pt_t) :: htp
327  type(point_t) :: p(8)
328  integer :: pt_idx, nelv
329  integer :: i, j, ierr
330 
331  nelv = dist%num_local()
332  element_offset = dist%start_idx()
333 
334  call htp%init(2*nel, ndim)
335  pt_idx = 0
336  if (ndim .eq. 2) then
337  mpi_offset = mpi_offset + element_offset * re2_data_xy_size
338  if (.not. v2_format) then
339  allocate(re2v1_data_xy(nelv))
340  call mpi_file_read_at_all(fh, mpi_offset, &
341  re2v1_data_xy, nelv, mpi_re2v1_data_xy, status, ierr)
342  do i = 1, nelv
343  do j = 1, 4
344  p(j) = point_t(real(re2v1_data_xy(i)%x(j),dp), &
345  real(re2v1_data_xy(i)%y(j),dp), 0.0d0)
346  call re2_file_add_point(htp, p(j), pt_idx)
347  end do
348  if (nelv > 10) then
349  if(mod(i,nelv/10) .eq. 0) write(*,*) i, 'elements read'
350  end if
351  ! swap vertices to keep symmetric vertex numbering in neko
352  call msh%add_element(i, p(1), p(2), p(4), p(3))
353  end do
354  deallocate(re2v1_data_xy)
355  else
356  allocate(re2v2_data_xy(nelv))
357  call mpi_file_read_at_all(fh, mpi_offset, &
358  re2v2_data_xy, nelv, mpi_re2v2_data_xy, status, ierr)
359  do i = 1, nelv
360  do j = 1, 4
361  p(j) = point_t(re2v2_data_xy(i)%x(j), &
362  re2v2_data_xy(i)%y(j), 0.0d0)
363  call re2_file_add_point(htp, p(j), pt_idx)
364  end do
365  if (nelv > 10) then
366  if(mod(i,nelv/10) .eq. 0) write(*,*) i, 'elements read'
367  end if
368  ! swap vertices to keep symmetric vertex numbering in neko
369  call msh%add_element(i, p(1), p(2), p(4), p(3))
370  end do
371  deallocate(re2v2_data_xy)
372  end if
373  else if (ndim .eq. 3) then
374  mpi_offset = mpi_offset + element_offset * re2_data_xyz_size
375  if (.not. v2_format) then
376  allocate(re2v1_data_xyz(nelv))
377  call mpi_file_read_at_all(fh, mpi_offset, &
378  re2v1_data_xyz, nelv, mpi_re2v1_data_xyz, status, ierr)
379  do i = 1, nelv
380  do j = 1, 8
381  p(j) = point_t(real(re2v1_data_xyz(i)%x(j),dp), &
382  real(re2v1_data_xyz(i)%y(j),dp),&
383  real(re2v1_data_xyz(i)%z(j),dp))
384  call re2_file_add_point(htp, p(j), pt_idx)
385  end do
386  if (nelv > 100) then
387  if(mod(i,nelv/100) .eq. 0) write(*,*) i, 'elements read'
388  end if
389  ! swap vertices to keep symmetric vertex numbering in neko
390  call msh%add_element(i, &
391  p(1), p(2), p(4), p(3), p(5), p(6), p(8), p(7))
392  end do
393  deallocate(re2v1_data_xyz)
394  else
395  allocate(re2v2_data_xyz(nelv))
396  call mpi_file_read_at_all(fh, mpi_offset, &
397  re2v2_data_xyz, nelv, mpi_re2v2_data_xyz, status, ierr)
398  do i = 1, nelv
399  do j = 1, 8
400  p(j) = point_t(re2v2_data_xyz(i)%x(j), &
401  re2v2_data_xyz(i)%y(j),&
402  re2v2_data_xyz(i)%z(j))
403  call re2_file_add_point(htp, p(j), pt_idx)
404  end do
405  if (nelv > 100) then
406  if(mod(i,nelv/100) .eq. 0) write(*,*) i, 'elements read'
407  end if
408  ! swap vertices to keep symmetric vertex numbering in neko
409  call msh%add_element(i, &
410  p(1), p(2), p(4), p(3), p(5), p(6), p(8), p(7))
411  end do
412  deallocate(re2v2_data_xyz)
413  end if
414  end if
415 
416  call htp%free()
417  end subroutine re2_file_read_points
418 
419  subroutine re2_file_read_curve(msh, ncurve, dist, fh, mpi_offset, v2_format)
420  type(mesh_t), intent(inout) :: msh
421  integer (kind=MPI_OFFSET_KIND) :: mpi_offset
422  integer, intent(inout) :: ncurve
423  type(linear_dist_t) :: dist
424  type(mpi_file), intent(inout) :: fh
425  logical, intent(in) :: v2_format
426  type(mpi_status) :: status
427  integer :: i, j, l, ierr, el_idx, id
428  type(re2v1_curve_t), allocatable :: re2v1_data_curve(:)
429  type(re2v2_curve_t), allocatable :: re2v2_data_curve(:)
430  real(kind=dp), allocatable :: curve_data(:,:,:)
431  integer, allocatable :: curve_type(:,:)
432  logical, allocatable :: curve_element(:)
433  character(len=1) :: chtemp
434  logical :: curve_skip = .false.
435 
436  allocate(curve_data(5,12,msh%nelv))
437  allocate(curve_element(msh%nelv))
438  allocate(curve_type(12,msh%nelv))
439  do i = 1, msh%nelv
440  curve_element(i) = .false.
441  do j = 1, 12
442  curve_type(j,i) = 0
443  do l = 1, 5
444  curve_data(l,j,i) = 0d0
445  end do
446  end do
447  end do
448 
449  call neko_log%message('Reading curved data')
450  if (.not. v2_format) then
451  allocate(re2v1_data_curve(ncurve))
452  call mpi_file_read_at_all(fh, mpi_offset, re2v1_data_curve, ncurve, &
453  mpi_re2v1_data_cv, status, ierr)
454  else
455  allocate(re2v2_data_curve(ncurve))
456  call mpi_file_read_at_all(fh, mpi_offset, re2v2_data_curve, ncurve, &
457  mpi_re2v2_data_cv, status, ierr)
458  end if
459  !This can probably be made nicer...
460  do i = 1, ncurve
461  if(v2_format) then
462  el_idx = re2v2_data_curve(i)%elem - dist%start_idx()
463  id = re2v2_data_curve(i)%zone
464  chtemp = re2v2_data_curve(i)%type
465  do j = 1, 5
466  curve_data(j,id, el_idx) = re2v2_data_curve(i)%point(j)
467  enddo
468  else
469  el_idx = re2v1_data_curve(i)%elem - dist%start_idx()
470  id = re2v1_data_curve(i)%zone
471  chtemp = re2v1_data_curve(i)%type
472  do j = 1, 5
473  curve_data(j,id, el_idx) = real(re2v1_data_curve(i)%point(j),dp)
474  enddo
475  end if
476 
477  curve_element(el_idx) = .true.
478  !This might need to be extended
479  select case(trim(chtemp))
480  case ('s')
481  curve_type(id,el_idx) = 1
482  call neko_log%warning('curve type s not supported, treating mesh as non-curved')
483  curve_skip = .true.
484  exit
485  case ('e')
486  curve_type(id,el_idx) = 2
487  call neko_log%warning('curve type e not supported, treating mesh as non-curved')
488  curve_skip = .true.
489  exit
490  case ('C')
491  curve_type(id,el_idx) = 3
492  case ('m')
493  curve_type(id,el_idx) = 4
494  case default
495  write(*,*) chtemp, 'curve type not supported yet, treating mesh as non-curved',id, el_idx
496 
497  curve_skip = .true.
498  end select
499  end do
500 
501  if( v2_format) then
502  deallocate(re2v2_data_curve)
503  else
504  deallocate(re2v1_data_curve)
505  end if
506  if (.not. curve_skip) then
507  do el_idx = 1, msh%nelv
508  if (curve_element(el_idx)) then
509  call msh%mark_curve_element(el_idx, &
510  curve_data(1,1,el_idx), curve_type(1,el_idx))
511  end if
512  end do
513  end if
514 
515  deallocate(curve_data)
516  deallocate(curve_element)
517  deallocate(curve_type)
518  end subroutine re2_file_read_curve
519 
520 
521  subroutine re2_file_read_bcs(msh, nbcs, dist, fh, mpi_offset, v2_format)
522  type(mesh_t), intent(inout) :: msh
523  integer (kind=MPI_OFFSET_KIND) :: mpi_offset
524  integer, intent(inout) :: nbcs
525  type(linear_dist_t) :: dist
526  type(mpi_file), intent(inout) :: fh
527  logical, intent(in) :: v2_format
528  type(mpi_status) :: status
529  integer :: pids(4)
530  integer :: sym_facet, label
531  integer :: p_el_idx, p_facet
532  integer :: i, j, ierr, el_idx
533  integer, parameter, dimension(6) :: facet_map = (/3, 2, 4, 1, 5, 6/)
534  logical :: periodic
535  type(re2v1_bc_t), allocatable :: re2v1_data_bc(:)
536  type(re2v2_bc_t), allocatable :: re2v2_data_bc(:)
537  character(len=LOG_SIZE) :: log_buf
538 
539  ! ---- Offsets for conversion from internal zones to labeled zones
540  integer :: user_labeled_zones(NEKO_MSH_MAX_ZLBLS)
541  integer :: labeled_zone_offsets(NEKO_MSH_MAX_ZLBLS)
542  integer :: total_labeled_zone_offset, current_internal_zone
543  total_labeled_zone_offset = 0
544  labeled_zone_offsets = 0
545  user_labeled_zones = 0
546  current_internal_zone = 1
547  ! ----
548 
549  call neko_log%message("Reading boundary conditions")
550  if (.not. v2_format) then
551  allocate(re2v1_data_bc(nbcs))
552  call mpi_file_read_at_all(fh, mpi_offset, re2v1_data_bc, nbcs, &
553  mpi_re2v1_data_bc, status, ierr)
554  else
555  allocate(re2v2_data_bc(nbcs))
556  call mpi_file_read_at_all(fh, mpi_offset, re2v2_data_bc, nbcs, &
557  mpi_re2v2_data_bc, status, ierr)
558  end if
559 
560  periodic = .false.
561 
563  if (v2_format) then ! V2 format
564 
565  !
566  ! Search for user labeled zones
567  !
568  do i = 1, nbcs
569 
570  el_idx = int(re2v2_data_bc(i)%elem) - dist%start_idx()
571  sym_facet = facet_map(int(re2v2_data_bc(i)%face))
572 
573  select case(trim(re2v2_data_bc(i)%type))
574  case ('MSH', 'msh', 'EXO','exo')
575 
576  label = int(re2v2_data_bc(i)%bc_data(5))
577 
578  if (label .lt. 1 .or. label .gt. neko_msh_max_zlbls) then
579  call neko_error('Invalid label id (valid range [1,...,20])')
580  end if
581 
582  if (user_labeled_zones(label) .eq. 0) then
583  write (log_buf, "(A,I2,A)") "Labeled zone ", label, &
584  " found."
585  call neko_log%message(log_buf)
586  user_labeled_zones(label) = 1
587  end if
588 
589  call msh%mark_labeled_facet(sym_facet, el_idx, label)
590 
591  ! Get the largest label possible in case we need to convert
592  ! re2 labels to labeled zones (see below).
593  total_labeled_zone_offset = max(total_labeled_zone_offset, label)
594 
595 
596  end select
597  end do
598 
599  do i = 1, nbcs
600  el_idx = int(re2v2_data_bc(i)%elem) - dist%start_idx()
601  sym_facet = facet_map(int(re2v2_data_bc(i)%face))
602  select case(trim(re2v2_data_bc(i)%type))
603  case ('MSH', 'msh', 'EXO','exo')
604  !Do nothing, already handled
605  case ('W')
606  if (neko_w_bc_label .eq. -1) then
607  neko_w_bc_label = current_internal_zone
608  current_internal_zone = current_internal_zone + 1
609  end if
610 
611  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
612  re2v2_data_bc(i)%type, &
613  neko_w_bc_label, total_labeled_zone_offset, &
614  labeled_zone_offsets(neko_w_bc_label) .eq. 0)
615 
616  labeled_zone_offsets(neko_w_bc_label) = 1
617 
618  case ('v', 'V')
619  if (neko_v_bc_label .eq. -1) then
620  neko_v_bc_label = current_internal_zone
621  current_internal_zone = current_internal_zone + 1
622  end if
623 
624  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
625  re2v2_data_bc(i)%type, &
626  neko_v_bc_label, total_labeled_zone_offset, &
627  labeled_zone_offsets(neko_v_bc_label) .eq. 0)
628 
629  labeled_zone_offsets(neko_v_bc_label) = 1
630  case ('O', 'o')
631  if (neko_o_bc_label .eq. -1) then
632  neko_o_bc_label = current_internal_zone
633  current_internal_zone = current_internal_zone + 1
634  end if
635 
636  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
637  re2v2_data_bc(i)%type, &
638  neko_o_bc_label, total_labeled_zone_offset, &
639  labeled_zone_offsets(neko_o_bc_label) .eq. 0)
640 
641  labeled_zone_offsets(neko_o_bc_label) = 1
642  case ('SYM','sym')
643  if (neko_sym_bc_label .eq. -1) then
644  neko_sym_bc_label = current_internal_zone
645  current_internal_zone = current_internal_zone + 1
646  end if
647 
648  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
649  re2v2_data_bc(i)%type, &
650  neko_sym_bc_label, total_labeled_zone_offset, &
651  labeled_zone_offsets(neko_sym_bc_label) .eq. 0)
652 
653  labeled_zone_offsets(neko_sym_bc_label) = 1
654  case ('ON', 'on')
655  if (neko_on_bc_label .eq. -1) then
656  neko_on_bc_label = current_internal_zone
657  current_internal_zone = current_internal_zone + 1
658  end if
659 
660  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
661  re2v2_data_bc(i)%type, &
662  neko_on_bc_label, total_labeled_zone_offset, &
663  labeled_zone_offsets(neko_on_bc_label) .eq. 0)
664 
665  labeled_zone_offsets(neko_on_bc_label) = 1
666  case ('s', 'sl', 'sh', 'shl', 'S', 'SL', 'SH', 'SHL')
667  if (neko_shl_bc_label .eq. -1) then
668  neko_shl_bc_label = current_internal_zone
669  current_internal_zone = current_internal_zone + 1
670  end if
671 
672  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
673  re2v2_data_bc(i)%type, &
674  neko_shl_bc_label, total_labeled_zone_offset, &
675  labeled_zone_offsets(neko_shl_bc_label) .eq. 0)
676 
677  labeled_zone_offsets(neko_shl_bc_label) = 1
678  case ('P')
679  periodic = .true.
680  p_el_idx = int(re2v2_data_bc(i)%bc_data(1))
681  p_facet = facet_map(int(re2v2_data_bc(i)%bc_data(2)))
682  call msh%get_facet_ids(sym_facet, el_idx, pids)
683  call msh%mark_periodic_facet(sym_facet, el_idx, &
684  p_facet, p_el_idx, pids)
685  case default
686  write (*,*) trim(re2v2_data_bc(i)%type), 'bc type not supported yet'
687  write (*,*) re2v2_data_bc(i)%bc_data
688  end select
689  end do
690 
691  !
692  ! Fix periodic condition for shared nodes
693  !
694  if (periodic) then
695  do j = 1, 3
696  do i = 1, nbcs
697  el_idx = re2v2_data_bc(i)%elem - dist%start_idx()
698  sym_facet = facet_map(int(re2v2_data_bc(i)%face))
699  select case(trim(re2v2_data_bc(i)%type))
700  case ('P')
701  p_el_idx = int(re2v2_data_bc(i)%bc_data(1))
702  p_facet = facet_map(int(re2v2_data_bc(i)%bc_data(2)))
703  call msh%create_periodic_ids(sym_facet, el_idx, &
704  p_facet, p_el_idx)
705  end select
706  end do
707  end do
708  end if
709  deallocate(re2v2_data_bc)
710 
711  else ! V! format
712 
713  !
714  ! Search for user labeled zones
715  !
716  do i = 1, nbcs
717 
718  el_idx = re2v1_data_bc(i)%elem - dist%start_idx()
719  sym_facet = facet_map(re2v1_data_bc(i)%face)
720 
721  select case(trim(re2v1_data_bc(i)%type))
722  case ('MSH', 'msh', 'EXO','exo')
723 
724  label = int(re2v1_data_bc(i)%bc_data(5))
725 
726  if (label .lt. 1 .or. label .gt. neko_msh_max_zlbls) then
727  call neko_error('Invalid label id (valid range [1,...,20])')
728  end if
729 
730  if (user_labeled_zones(label) .eq. 0) then
731  write (log_buf, "(A,I2,A,I3)") "Labeled zone ", label, &
732  " found."
733  call neko_log%message(log_buf)
734  user_labeled_zones(label) = 1
735  end if
736 
737  ! Get the largest label possible in case we need to convert
738  ! re2 labels to labeled zones (see below).
739  total_labeled_zone_offset = max(total_labeled_zone_offset, label)
740 
741  call msh%mark_labeled_facet(sym_facet, el_idx, label)
742 
743  end select
744  end do
745 
746  do i = 1, nbcs
747 
748  el_idx = re2v1_data_bc(i)%elem - dist%start_idx()
749  sym_facet = facet_map(re2v1_data_bc(i)%face)
750  select case(trim(re2v1_data_bc(i)%type))
751  case ('MSH', 'msh', 'EXO','exo')
752  !Do nothing, already handled
753  case ('W')
754  if (neko_w_bc_label .eq. -1) then
755  neko_w_bc_label = current_internal_zone
756  current_internal_zone = current_internal_zone + 1
757  end if
758 
759  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
760  re2v1_data_bc(i)%type, &
761  neko_w_bc_label, total_labeled_zone_offset, &
762  labeled_zone_offsets(neko_w_bc_label) .eq. 0)
763 
764  labeled_zone_offsets(neko_w_bc_label) = 1
765 
766  case ('v', 'V')
767  if (neko_v_bc_label .eq. -1) then
768  neko_v_bc_label = current_internal_zone
769  current_internal_zone = current_internal_zone + 1
770  end if
771 
772  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
773  re2v1_data_bc(i)%type, &
774  neko_v_bc_label, total_labeled_zone_offset, &
775  labeled_zone_offsets(neko_v_bc_label) .eq. 0)
776 
777  labeled_zone_offsets(neko_v_bc_label) = 1
778  case ('O', 'o')
779  if (neko_o_bc_label .eq. -1) then
780  neko_o_bc_label = current_internal_zone
781  current_internal_zone = current_internal_zone + 1
782  end if
783 
784  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
785  re2v1_data_bc(i)%type, &
786  neko_o_bc_label, total_labeled_zone_offset, &
787  labeled_zone_offsets(neko_o_bc_label) .eq. 0)
788 
789  labeled_zone_offsets(neko_o_bc_label) = 1
790  case ('SYM','sym')
791  if (neko_sym_bc_label .eq. -1) then
792  neko_sym_bc_label = current_internal_zone
793  current_internal_zone = current_internal_zone + 1
794  end if
795 
796  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
797  re2v1_data_bc(i)%type, &
798  neko_sym_bc_label, total_labeled_zone_offset, &
799  labeled_zone_offsets(neko_sym_bc_label) .eq. 0)
800 
801  labeled_zone_offsets(neko_sym_bc_label) = 1
802  case ('ON', 'on')
803  if (neko_on_bc_label .eq. -1) then
804  neko_on_bc_label = current_internal_zone
805  current_internal_zone = current_internal_zone + 1
806  end if
807 
808  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
809  re2v1_data_bc(i)%type, &
810  neko_on_bc_label, total_labeled_zone_offset, &
811  labeled_zone_offsets(neko_on_bc_label) .eq. 0)
812 
813  labeled_zone_offsets(neko_on_bc_label) = 1
814  case ('s', 'sl', 'sh', 'shl', 'S', 'SL', 'SH', 'SHL')
815  if (neko_shl_bc_label .eq. -1) then
816  neko_shl_bc_label = current_internal_zone
817  current_internal_zone = current_internal_zone + 1
818  end if
819 
820  call re2_file_mark_labeled_bc(msh, el_idx, sym_facet, &
821  re2v1_data_bc(i)%type, &
822  neko_shl_bc_label, total_labeled_zone_offset, &
823  labeled_zone_offsets(neko_shl_bc_label) .eq. 0)
824 
825  labeled_zone_offsets(neko_shl_bc_label) = 1
826  case ('P')
827  periodic = .true.
828  p_el_idx = int(re2v1_data_bc(i)%bc_data(1))
829  p_facet = facet_map(int(re2v1_data_bc(i)%bc_data(2)))
830  call msh%get_facet_ids(sym_facet, el_idx, pids)
831  call msh%mark_periodic_facet(sym_facet, el_idx, &
832  p_facet, p_el_idx, pids)
833  case default
834  write (*,*) re2v1_data_bc(i)%type, 'bc type not supported yet'
835  write (*,*) re2v1_data_bc(i)%bc_data
836  end select
837  end do
838 
839  !
840  ! Fix periodic condition for shared nodes
841  !
842  if (periodic) then
843  do j = 1, 3
844  do i = 1, nbcs
845  el_idx = re2v1_data_bc(i)%elem - dist%start_idx()
846  sym_facet = facet_map(re2v1_data_bc(i)%face)
847  select case(trim(re2v1_data_bc(i)%type))
848  case ('P')
849  p_el_idx = int(re2v1_data_bc(i)%bc_data(1))
850  p_facet = facet_map(int(re2v1_data_bc(i)%bc_data(2)))
851  call msh%create_periodic_ids(sym_facet, el_idx, &
852  p_facet, p_el_idx)
853  end select
854  end do
855  end do
856  end if
857 
858  deallocate(re2v1_data_bc)
859  end if
860 
861  end subroutine re2_file_read_bcs
862 
873  subroutine re2_file_mark_labeled_bc(msh, el_idx, facet, type, label, offset, print_info)
874  type(mesh_t), intent(inout) :: msh
875  integer, intent(in) :: el_idx
876  integer, intent(in) :: facet
877  character(len=*), intent(in) :: type
878  integer, intent(in) :: label
879  integer, intent(in) :: offset
880  logical, intent(in) :: print_info
881 
882  integer :: mark_label
883  character(len=LOG_SIZE) :: log_buf
884 
885  mark_label = offset + label
886 
887  if (mark_label .lt. 1 .or. mark_label .gt. neko_msh_max_zlbls) then
888  call neko_error("You have reached the maximum amount of allowed labeled&
889 & zones (max allowed: 20). This happened when converting re2 internal labels&
890 & like e.g. 'w', 'V' or 'o' to labeled zones. Please reduce the number of&
891 & labeled zones that you have defined or make sure that they are labeled&
892 & from [1,...,20].")
893  end if
894 
895  if (print_info) then
896  write (log_buf, "(A3,A,I2)") trim(type), " => Labeled index ", mark_label
897  call neko_log%message(log_buf)
898  end if
899  call msh%mark_labeled_facet(facet, el_idx, mark_label)
900 
901  end subroutine re2_file_mark_labeled_bc
902 
903  subroutine re2_file_add_point(htp, p, idx)
904  type(htable_pt_t), intent(inout) :: htp
905  type(point_t), intent(inout) :: p
906  integer, intent(inout) :: idx
907  integer :: tmp
908 
909  if (htp%get(p, tmp) .gt. 0) then
910  idx = idx + 1
911  call htp%set(p, idx)
912  call p%set_id(idx)
913  else
914  call p%set_id(tmp)
915  end if
916 
917  end subroutine re2_file_add_point
918 
919 end module re2_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
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
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
integer, parameter, public log_size
Definition: log.f90:42
NEKTON map file.
Definition: map_file.f90:35
NEKTON map.
Definition: map.f90:3
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
type(mpi_datatype), public mpi_re2v2_data_xyz
MPI derived type for 3D NEKTON re2 data.
Definition: mpi_types.f90:52
type(mpi_datatype), public mpi_re2v2_data_xy
MPI derived type for 2D NEKTON re2 data.
Definition: mpi_types.f90:53
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
type(mpi_datatype), public mpi_re2v1_data_cv
MPI derived type for NEKTON re2 cv data.
Definition: mpi_types.f90:49
type(mpi_datatype), public mpi_re2v2_data_cv
MPI derived type for NEKTON re2 cv data.
Definition: mpi_types.f90:54
type(mpi_datatype), public mpi_re2v1_data_xy
MPI derived type for 2D NEKTON re2 data.
Definition: mpi_types.f90:48
type(mpi_datatype), public mpi_re2v1_data_bc
MPI derived type for NEKTON re2 bc data.
Definition: mpi_types.f90:50
type(mpi_datatype), public mpi_re2v2_data_bc
MPI derived type for NEKTON re2 bc data.
Definition: mpi_types.f90:55
type(mpi_datatype), public mpi_re2v1_data_xyz
MPI derived type for 3D NEKTON re2 data.
Definition: mpi_types.f90:47
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
subroutine re2_file_mark_labeled_bc(msh, el_idx, facet, type, label, offset, print_info)
Mark a labeled zone based on an offset.
Definition: re2_file.f90:874
subroutine re2_file_read(this, data)
Load a binary NEKTON mesh from a re2 file.
Definition: re2_file.f90:72
subroutine re2_file_write(this, data, t)
Definition: re2_file.f90:219
integer, public neko_on_bc_label
Definition: re2_file.f90:58
integer, public neko_w_bc_label
Definition: re2_file.f90:54
integer, public neko_o_bc_label
Definition: re2_file.f90:56
subroutine re2_file_read_bcs(msh, nbcs, dist, fh, mpi_offset, v2_format)
Definition: re2_file.f90:522
subroutine re2_file_add_point(htp, p, idx)
Definition: re2_file.f90:904
integer, public neko_v_bc_label
Definition: re2_file.f90:55
integer, public neko_shl_bc_label
Definition: re2_file.f90:59
subroutine re2_file_read_curve(msh, ncurve, dist, fh, mpi_offset, v2_format)
Definition: re2_file.f90:420
integer, public neko_sym_bc_label
Definition: re2_file.f90:57
subroutine re2_file_read_points(msh, ndim, nel, dist, fh, mpi_offset, re2_data_xy_size, re2_data_xyz_size, v2_format)
Definition: re2_file.f90:311
NEKTON re2 format.
Definition: re2.f90:2
real(kind=sp), parameter re2_endian_test
NEKTION re2 endian test.
Definition: re2.f90:10
integer, parameter re2_hdr_size
NEKTON re2 header size.
Definition: re2.f90:7
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
NEKTON re2 bc data (version 1)
Definition: re2.f90:39
NEKTON re2 curve data (version 1)
Definition: re2.f90:31
NEKTON re2 element data (2d) (version 1)
Definition: re2.f90:25
NEKTON re2 element data (3d) (version 1)
Definition: re2.f90:18
NEKTON re2 bc data (version 2)
Definition: re2.f90:73
NEKTON re2 curve data (version 2)
Definition: re2.f90:65
NEKTON re2 element data (2d) (version 2)
Definition: re2.f90:59
NEKTON re2 element data (3d) (version 2)
Definition: re2.f90:52
Interface for NEKTON re2 files.
Definition: re2_file.f90:62
#define max(a, b)
Definition: tensor.cu:40