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
81 integer (kind=MPI_OFFSET_KIND) :: mpi_offset
84 integer :: ncurv, nbcs
88 character(len=1024) :: map_fname
90 integer :: re2_data_xy_size
91 integer :: re2_data_xyz_size
92 integer :: re2_data_cv_size
93 integer :: re2_data_bc_size
95 character(len=LOG_SIZE) :: log_buf
97 call this%check_exists()
107 open(unit=9,
file=trim(this%fname), status=
'old', iostat=ierr)
108 call neko_log%message(
'Reading binary NEKTON file ' // this%fname)
110 read(9,
'(a80)') hdr_full
111 read(hdr_full,
'(a5)') hdr_ver
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
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
119 else if (hdr_ver .eq.
'#v001')
then
120 read(hdr_full,
'(a5,i9,i3,i9,a54)') hdr_ver, nel, ndim, nelv, hdr_str
135 write(log_buf,1) ndim, nelv
136 1
format(
'gdim = ', i1,
', nelements =', i7)
142 inquire(
file=map_fname, exist=read_map)
148 call neko_log%warning(
'No NEKTON map file found')
151 call mpi_file_open(
neko_comm, trim(this%fname), &
152 mpi_mode_rdonly, mpi_info_null, fh, ierr)
154 if (ierr .ne. 0)
then
155 call neko_log%error(
"Can't open binary NEKTON file ")
160 call msh%init(ndim, dist)
165 call mpi_file_read_at_all(fh, mpi_offset, test, 1, mpi_real, status, ierr)
169 call neko_error(
'Invalid endian of re2 file, byte swap not implemented yet')
173 mpi_offset, re2_data_xy_size, re2_data_xyz_size, v2_format)
178 if (ndim .eq. 2)
then
179 mpi_offset = mpi_offset + int(dist%num_global(),i8) * int(re2_data_xy_size,i8)
181 mpi_offset = mpi_offset + int(dist%num_global(),i8) * int(re2_data_xyz_size,i8)
187 call mpi_file_read_at_all(fh, mpi_offset, t2, 1, mpi_double_precision, status, ierr)
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)
198 call mpi_file_read_at_all(fh, mpi_offset, ncurv, 1, mpi_integer, status, ierr)
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)
209 call mpi_file_close(fh, ierr)
220 class(*),
target,
intent(in) :: data
221 real(kind=
rp),
intent(in),
optional :: t
222 type(
re2v1_xy_t),
allocatable :: re2_data_xy(:)
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
230 integer (kind=MPI_OFFSET_KIND) :: mpi_offset
231 integer :: element_offset
232 integer :: re2_data_xy_size
233 integer :: re2_data_xyz_size
244 call mpi_reduce(msh%nelv, nelgv, 1, &
245 mpi_integer, mpi_sum, 0,
neko_comm, ierr)
247 call mpi_exscan(msh%nelv, element_offset, 1, &
250 call neko_log%message(
'Writing data as a binary NEKTON file ' // this%fname)
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,&
260 call mpi_file_open(
neko_comm, trim(this%fname), &
261 mpi_mode_wronly + mpi_mode_create, mpi_info_null, fh, ierr)
265 mpi_real, status, ierr)
268 if (msh%gdim .eq. 2)
then
269 allocate(re2_data_xy(msh%nelv))
271 re2_data_xy(i)%rgroup = 1.0
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))
277 mpi_offset = mpi_offset + element_offset * re2_data_xy_size
278 call mpi_file_write_at(fh, mpi_offset, &
281 deallocate(re2_data_xy)
283 else if (msh%gdim .eq. 3)
then
284 allocate(re2_data_xyz(msh%nelv))
286 re2_data_xyz(i)%rgroup = 1.0
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))
293 mpi_offset = mpi_offset + element_offset * re2_data_xyz_size
294 call mpi_file_write_at(fh, mpi_offset, &
297 deallocate(re2_data_xyz)
302 call mpi_file_close(fh, ierr)
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
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
328 integer :: pt_idx, nelv
329 integer :: i, j, ierr
331 nelv = dist%num_local()
332 element_offset = dist%start_idx()
334 call htp%init(2*nel, ndim)
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, &
345 real(re2v1_data_xy(i)%y(j),dp), 0.0d0)
349 if(mod(i,nelv/10) .eq. 0)
write(*,*) i,
'elements read'
352 call msh%add_element(i, p(1), p(2), p(4), p(3))
354 deallocate(re2v1_data_xy)
356 allocate(re2v2_data_xy(nelv))
357 call mpi_file_read_at_all(fh, mpi_offset, &
361 p(j) =
point_t(re2v2_data_xy(i)%x(j), &
362 re2v2_data_xy(i)%y(j), 0.0d0)
366 if(mod(i,nelv/10) .eq. 0)
write(*,*) i,
'elements read'
369 call msh%add_element(i, p(1), p(2), p(4), p(3))
371 deallocate(re2v2_data_xy)
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, &
382 real(re2v1_data_xyz(i)%y(j),dp),&
383 real(re2v1_data_xyz(i)%z(j),dp))
387 if(mod(i,nelv/100) .eq. 0)
write(*,*) i,
'elements read'
390 call msh%add_element(i, &
391 p(1), p(2), p(4), p(3), p(5), p(6), p(8), p(7))
393 deallocate(re2v1_data_xyz)
395 allocate(re2v2_data_xyz(nelv))
396 call mpi_file_read_at_all(fh, mpi_offset, &
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))
406 if(mod(i,nelv/100) .eq. 0)
write(*,*) i,
'elements read'
409 call msh%add_element(i, &
410 p(1), p(2), p(4), p(3), p(5), p(6), p(8), p(7))
412 deallocate(re2v2_data_xyz)
420 type(
mesh_t),
intent(inout) :: msh
421 integer (kind=MPI_OFFSET_KIND) :: mpi_offset
422 integer,
intent(inout) :: ncurve
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
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.
436 allocate(curve_data(5,12,msh%nelv))
437 allocate(curve_element(msh%nelv))
438 allocate(curve_type(12,msh%nelv))
440 curve_element(i) = .false.
444 curve_data(l,j,i) = 0d0
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, &
455 allocate(re2v2_data_curve(ncurve))
456 call mpi_file_read_at_all(fh, mpi_offset, re2v2_data_curve, ncurve, &
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
466 curve_data(j,id, el_idx) = re2v2_data_curve(i)%point(j)
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
473 curve_data(j,id, el_idx) =
real(re2v1_data_curve(i)%point(j),dp)
477 curve_element(el_idx) = .true.
479 select case(trim(chtemp))
481 curve_type(id,el_idx) = 1
482 call neko_log%warning(
'curve type s not supported, treating mesh as non-curved')
486 curve_type(id,el_idx) = 2
487 call neko_log%warning(
'curve type e not supported, treating mesh as non-curved')
491 curve_type(id,el_idx) = 3
493 curve_type(id,el_idx) = 4
495 write(*,*) chtemp,
'curve type not supported yet, treating mesh as non-curved',id, el_idx
502 deallocate(re2v2_data_curve)
504 deallocate(re2v1_data_curve)
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))
515 deallocate(curve_data)
516 deallocate(curve_element)
517 deallocate(curve_type)
522 type(
mesh_t),
intent(inout) :: msh
523 integer (kind=MPI_OFFSET_KIND) :: mpi_offset
524 integer,
intent(inout) :: nbcs
526 type(mpi_file),
intent(inout) :: fh
527 logical,
intent(in) :: v2_format
528 type(mpi_status) :: status
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/)
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
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
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, &
555 allocate(re2v2_data_bc(nbcs))
556 call mpi_file_read_at_all(fh, mpi_offset, re2v2_data_bc, nbcs, &
570 el_idx = int(re2v2_data_bc(i)%elem) - dist%start_idx()
571 sym_facet = facet_map(int(re2v2_data_bc(i)%face))
573 select case(trim(re2v2_data_bc(i)%type))
576 label = int(re2v2_data_bc(i)%bc_data(5))
578 if (label .lt. 1 .or. label .gt. neko_msh_max_zlbls)
then
579 call neko_error(
'Invalid label id (valid range [1,...,20])')
582 if (user_labeled_zones(label) .eq. 0)
then
583 write (log_buf,
"(A,I2,A)")
"Labeled zone ", label, &
586 user_labeled_zones(label) = 1
589 call msh%mark_labeled_facet(sym_facet, el_idx, label)
593 total_labeled_zone_offset =
max(total_labeled_zone_offset, label)
600 el_idx = int(re2v2_data_bc(i)%elem) - dist%start_idx()
601 sym_facet = facet_map(int(re2v2_data_bc(i)%face))
603 select case(trim(re2v2_data_bc(i)%type))
607 current_internal_zone = current_internal_zone + 1
611 re2v2_data_bc(i)%type, &
620 current_internal_zone = current_internal_zone + 1
624 re2v2_data_bc(i)%type, &
632 current_internal_zone = current_internal_zone + 1
636 re2v2_data_bc(i)%type, &
644 current_internal_zone = current_internal_zone + 1
648 re2v2_data_bc(i)%type, &
656 current_internal_zone = current_internal_zone + 1
660 re2v2_data_bc(i)%type, &
665 case (
's',
'sl',
'sh',
'shl',
'S',
'SL',
'SH',
'SHL')
668 current_internal_zone = current_internal_zone + 1
672 re2v2_data_bc(i)%type, &
679 p_el_idx = int(re2v2_data_bc(i)%bc_data(1))
680 p_facet = facet_map(int(re2v2_data_bc(i)%bc_data(2)))
681 call msh%get_facet_ids(sym_facet, el_idx, pids)
682 call msh%mark_periodic_facet(sym_facet, el_idx, &
683 p_facet, p_el_idx, pids)
685 write (*,*) re2v1_data_bc(i)%type,
'bc type not supported yet'
686 write (*,*) re2v1_data_bc(i)%bc_data
696 el_idx = re2v2_data_bc(i)%elem - dist%start_idx()
697 sym_facet = facet_map(int(re2v2_data_bc(i)%face))
698 select case(trim(re2v2_data_bc(i)%type))
700 p_el_idx = int(re2v2_data_bc(i)%bc_data(1))
701 p_facet = facet_map(int(re2v2_data_bc(i)%bc_data(2)))
702 call msh%create_periodic_ids(sym_facet, el_idx, &
708 deallocate(re2v2_data_bc)
717 el_idx = re2v1_data_bc(i)%elem - dist%start_idx()
718 sym_facet = facet_map(re2v1_data_bc(i)%face)
720 select case(trim(re2v1_data_bc(i)%type))
723 label = int(re2v1_data_bc(i)%bc_data(5))
725 if (label .lt. 1 .or. label .gt. neko_msh_max_zlbls)
then
726 call neko_error(
'Invalid label id (valid range [1,...,20])')
729 if (user_labeled_zones(label) .eq. 0)
then
730 write (log_buf,
"(A,I2,A,I3)")
"Labeled zone ", label, &
733 user_labeled_zones(label) = 1
738 total_labeled_zone_offset =
max(total_labeled_zone_offset, label)
740 call msh%mark_labeled_facet(sym_facet, el_idx, label)
747 el_idx = re2v1_data_bc(i)%elem - dist%start_idx()
748 sym_facet = facet_map(re2v1_data_bc(i)%face)
749 select case(trim(re2v1_data_bc(i)%type))
753 current_internal_zone = current_internal_zone + 1
757 re2v1_data_bc(i)%type, &
766 current_internal_zone = current_internal_zone + 1
770 re2v1_data_bc(i)%type, &
778 current_internal_zone = current_internal_zone + 1
782 re2v1_data_bc(i)%type, &
790 current_internal_zone = current_internal_zone + 1
794 re2v1_data_bc(i)%type, &
802 current_internal_zone = current_internal_zone + 1
806 re2v1_data_bc(i)%type, &
811 case (
's',
'sl',
'sh',
'shl',
'S',
'SL',
'SH',
'SHL')
814 current_internal_zone = current_internal_zone + 1
818 re2v1_data_bc(i)%type, &
825 p_el_idx = int(re2v1_data_bc(i)%bc_data(1))
826 p_facet = facet_map(int(re2v1_data_bc(i)%bc_data(2)))
827 call msh%get_facet_ids(sym_facet, el_idx, pids)
828 call msh%mark_periodic_facet(sym_facet, el_idx, &
829 p_facet, p_el_idx, pids)
831 write (*,*) re2v1_data_bc(i)%type,
'bc type not supported yet'
832 write (*,*) re2v1_data_bc(i)%bc_data
842 el_idx = re2v1_data_bc(i)%elem - dist%start_idx()
843 sym_facet = facet_map(re2v1_data_bc(i)%face)
844 select case(trim(re2v1_data_bc(i)%type))
846 p_el_idx = int(re2v1_data_bc(i)%bc_data(1))
847 p_facet = facet_map(int(re2v1_data_bc(i)%bc_data(2)))
848 call msh%create_periodic_ids(sym_facet, el_idx, &
855 deallocate(re2v1_data_bc)
871 type(
mesh_t),
intent(inout) :: msh
872 integer,
intent(in) :: el_idx
873 integer,
intent(in) :: facet
874 character(len=*),
intent(in) :: type
875 integer,
intent(in) :: label
876 integer,
intent(in) :: offset
877 logical,
intent(in) :: print_info
879 integer :: mark_label
880 character(len=LOG_SIZE) :: log_buf
882 mark_label = offset + label
884 if (mark_label .lt. 1 .or. mark_label .gt. neko_msh_max_zlbls)
then
885 call neko_error(
"You have reached the maximum amount of allowed labeled&
886 & zones (max allowed: 20). This happened when converting re2 internal labels&
887 & like e.g. 'w', 'V' or 'o' to labeled zones. Please reduce the number of&
888 & labeled zones that you have defined or make sure that they are labeled&
893 write (log_buf,
"(A3,A,I2)") trim(type),
" => Labeled index ", mark_label
896 call msh%mark_labeled_facet(facet, el_idx, mark_label)
902 type(
point_t),
intent(inout) :: p
903 integer,
intent(inout) :: idx
906 if (htp%get(p, tmp) .gt. 0)
then
type(mpi_comm) neko_comm
MPI communicator.
integer pe_size
MPI size of communicator.
Defines practical data distributions.
Module for file I/O operations.
Implements a hash table ADT.
type(log_t), public neko_log
Global log stream.
integer, parameter, public log_size
integer, public mpi_double_precision_size
Size of MPI type double precision.
type(mpi_datatype), public mpi_re2v2_data_xyz
MPI derived type for 3D NEKTON re2 data.
type(mpi_datatype), public mpi_re2v2_data_xy
MPI derived type for 2D NEKTON re2 data.
integer, public mpi_character_size
Size of MPI type character.
integer, public mpi_real_size
Size of MPI type real.
integer, public mpi_integer_size
Size of MPI type integer.
type(mpi_datatype), public mpi_re2v1_data_cv
MPI derived type for NEKTON re2 cv data.
type(mpi_datatype), public mpi_re2v2_data_cv
MPI derived type for NEKTON re2 cv data.
type(mpi_datatype), public mpi_re2v1_data_xy
MPI derived type for 2D NEKTON re2 data.
type(mpi_datatype), public mpi_re2v1_data_bc
MPI derived type for NEKTON re2 bc data.
type(mpi_datatype), public mpi_re2v2_data_bc
MPI derived type for NEKTON re2 bc data.
type(mpi_datatype), public mpi_re2v1_data_xyz
MPI derived type for 3D NEKTON re2 data.
integer, parameter, public rp
Global precision used in computations.
NEKTON mesh data in re2 format.
subroutine re2_file_mark_labeled_bc(msh, el_idx, facet, type, label, offset, print_info)
Mark a labeled zone based on an offset.
subroutine re2_file_read(this, data)
Load a binary NEKTON mesh from a re2 file.
subroutine re2_file_write(this, data, t)
integer, public neko_on_bc_label
integer, public neko_w_bc_label
integer, public neko_o_bc_label
subroutine re2_file_read_bcs(msh, nbcs, dist, fh, mpi_offset, v2_format)
subroutine re2_file_add_point(htp, p, idx)
integer, public neko_v_bc_label
integer, public neko_shl_bc_label
subroutine re2_file_read_curve(msh, ncurve, dist, fh, mpi_offset, v2_format)
integer, public neko_sym_bc_label
subroutine re2_file_read_points(msh, ndim, nel, dist, fh, mpi_offset, re2_data_xy_size, re2_data_xyz_size, v2_format)
real(kind=sp), parameter re2_endian_test
NEKTION re2 endian test.
integer, parameter re2_hdr_size
NEKTON re2 header size.
subroutine, public filename_chsuffix(fname, new_fname, new_suffix)
Change a filename's suffix.
Load-balanced linear distribution .
Interface for NEKTON map files.
A point in with coordinates .
NEKTON re2 bc data (version 1)
NEKTON re2 curve data (version 1)
NEKTON re2 element data (2d) (version 1)
NEKTON re2 element data (3d) (version 1)
NEKTON re2 bc data (version 2)
NEKTON re2 curve data (version 2)
NEKTON re2 element data (2d) (version 2)
NEKTON re2 element data (3d) (version 2)
Interface for NEKTON re2 files.