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))
574 case (
'MSH',
'msh',
'EXO',
'exo')
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))
602 select case(trim(re2v2_data_bc(i)%type))
603 case (
'MSH',
'msh',
'EXO',
'exo')
608 current_internal_zone = current_internal_zone + 1
612 re2v2_data_bc(i)%type, &
621 current_internal_zone = current_internal_zone + 1
625 re2v2_data_bc(i)%type, &
633 current_internal_zone = current_internal_zone + 1
637 re2v2_data_bc(i)%type, &
645 current_internal_zone = current_internal_zone + 1
649 re2v2_data_bc(i)%type, &
657 current_internal_zone = current_internal_zone + 1
661 re2v2_data_bc(i)%type, &
666 case (
's',
'sl',
'sh',
'shl',
'S',
'SL',
'SH',
'SHL')
669 current_internal_zone = current_internal_zone + 1
673 re2v2_data_bc(i)%type, &
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)
686 write (*,*) trim(re2v2_data_bc(i)%type),
'bc type not supported yet'
687 write (*,*) re2v2_data_bc(i)%bc_data
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))
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, &
709 deallocate(re2v2_data_bc)
718 el_idx = re2v1_data_bc(i)%elem - dist%start_idx()
719 sym_facet = facet_map(re2v1_data_bc(i)%face)
721 select case(trim(re2v1_data_bc(i)%type))
722 case (
'MSH',
'msh',
'EXO',
'exo')
724 label = int(re2v1_data_bc(i)%bc_data(5))
726 if (label .lt. 1 .or. label .gt. neko_msh_max_zlbls)
then
727 call neko_error(
'Invalid label id (valid range [1,...,20])')
730 if (user_labeled_zones(label) .eq. 0)
then
731 write (log_buf,
"(A,I2,A,I3)")
"Labeled zone ", label, &
734 user_labeled_zones(label) = 1
739 total_labeled_zone_offset =
max(total_labeled_zone_offset, label)
741 call msh%mark_labeled_facet(sym_facet, el_idx, label)
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')
756 current_internal_zone = current_internal_zone + 1
760 re2v1_data_bc(i)%type, &
769 current_internal_zone = current_internal_zone + 1
773 re2v1_data_bc(i)%type, &
781 current_internal_zone = current_internal_zone + 1
785 re2v1_data_bc(i)%type, &
793 current_internal_zone = current_internal_zone + 1
797 re2v1_data_bc(i)%type, &
805 current_internal_zone = current_internal_zone + 1
809 re2v1_data_bc(i)%type, &
814 case (
's',
'sl',
'sh',
'shl',
'S',
'SL',
'SH',
'SHL')
817 current_internal_zone = current_internal_zone + 1
821 re2v1_data_bc(i)%type, &
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)
834 write (*,*) re2v1_data_bc(i)%type,
'bc type not supported yet'
835 write (*,*) re2v1_data_bc(i)%bc_data
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))
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, &
858 deallocate(re2v1_data_bc)
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
882 integer :: mark_label
883 character(len=LOG_SIZE) :: log_buf
885 mark_label = offset + label
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&
896 write (log_buf,
"(A3,A,I2)") trim(type),
" => Labeled index ", mark_label
899 call msh%mark_labeled_facet(facet, el_idx, mark_label)
905 type(
point_t),
intent(inout) :: p
906 integer,
intent(inout) :: idx
909 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.