From 40174ea0fb74a1addba4f74407cc7af12cc525d0 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Mon, 7 Mar 2022 09:53:31 +0000 Subject: [PATCH 1/5] Begin introducing IO back into x3div --- Makefile | 4 +- decomp2d/io.f90 | 1198 +++++++++++++++++++++++++++++++++ decomp2d/io_read_inflow.f90 | 52 ++ decomp2d/io_read_one.inc | 65 ++ decomp2d/io_read_var.inc | 69 ++ decomp2d/io_write_every.inc | 224 ++++++ decomp2d/io_write_one.inc | 84 +++ decomp2d/io_write_outflow.f90 | 52 ++ decomp2d/io_write_plane.inc | 148 ++++ decomp2d/io_write_var.inc | 69 ++ src/visu.f90 | 541 +++++++++++++++ src/x3d_tools.f90 | 10 +- 12 files changed, 2513 insertions(+), 3 deletions(-) create mode 100644 decomp2d/io.f90 create mode 100644 decomp2d/io_read_inflow.f90 create mode 100644 decomp2d/io_read_one.inc create mode 100644 decomp2d/io_read_var.inc create mode 100644 decomp2d/io_write_every.inc create mode 100644 decomp2d/io_write_one.inc create mode 100644 decomp2d/io_write_outflow.f90 create mode 100644 decomp2d/io_write_plane.inc create mode 100644 decomp2d/io_write_var.inc create mode 100644 src/visu.f90 diff --git a/Makefile b/Makefile index 82e4dca..a156290 100644 --- a/Makefile +++ b/Makefile @@ -59,10 +59,10 @@ DECOMPDIR = ./decomp2d SRCDIR = ./src ### List of files for the main code -SRCDECOMP = $(DECOMPDIR)/decomp_2d.f90 $(DECOMPDIR)/glassman.f90 $(DECOMPDIR)/fft_$(FFT).f90 +SRCDECOMP = $(DECOMPDIR)/decomp_2d.f90 $(DECOMPDIR)/glassman.f90 $(DECOMPDIR)/fft_$(FFT).f90 $(DECOMPDIR)/io.f90 OBJDECOMP = $(SRCDECOMP:%.f90=%.o) OBJ = $(SRC:%.f90=%.o) -SRC = $(SRCDIR)/x3d_precision.f90 $(SRCDIR)/module_param.f90 $(SRCDIR)/time_integrators.f90 $(SRCDIR)/tools.f90 $(SRCDIR)/x3d_transpose.f90 $(SRCDIR)/var.f90 $(SRCDIR)/thomas.f90 $(SRCDIR)/x3d_operator_x_data.f90 $(SRCDIR)/x3d_operator_y_data.f90 $(SRCDIR)/x3d_operator_z_data.f90 $(SRCDIR)/x3d_operator_1d.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/x3d_derive.f90 $(SRCDIR)/x3d_staggered.f90 $(SRCDIR)/x3d_filters.f90 $(SRCDIR)/bc_tgv.f90 $(SRCDIR)/navier.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/bc_tgv2d.f90 $(SRCDIR)/case.f90 $(SRCDIR)/transeq.f90 $(SRCDIR)/x3d_tools.f90 $(SRCDIR)/xcompact3d.f90 +SRC = $(SRCDIR)/x3d_precision.f90 $(SRCDIR)/module_param.f90 $(SRCDIR)/time_integrators.f90 $(SRCDIR)/tools.f90 $(SRCDIR)/x3d_transpose.f90 $(SRCDIR)/var.f90 $(SRCDIR)/thomas.f90 $(SRCDIR)/x3d_operator_x_data.f90 $(SRCDIR)/x3d_operator_y_data.f90 $(SRCDIR)/x3d_operator_z_data.f90 $(SRCDIR)/x3d_operator_1d.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/x3d_derive.f90 $(SRCDIR)/x3d_staggered.f90 $(SRCDIR)/x3d_filters.f90 $(SRCDIR)/bc_tgv.f90 $(SRCDIR)/navier.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/bc_tgv2d.f90 $(SRCDIR)/case.f90 $(SRCDIR)/transeq.f90 $(SRCDIR)/visu.f90 $(SRCDIR)/x3d_tools.f90 $(SRCDIR)/xcompact3d.f90 #######FFT settings########## ifeq ($(FFT),fftw3) diff --git a/decomp2d/io.f90 b/decomp2d/io.f90 new file mode 100644 index 0000000..ae14ad9 --- /dev/null +++ b/decomp2d/io.f90 @@ -0,0 +1,1198 @@ +!======================================================================= +! This is part of the 2DECOMP&FFT library +! +! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) +! decomposition. It also implements a highly scalable distributed +! three-dimensional Fast Fourier Transform (FFT). +! +! Copyright (C) 2009-2013 Ning Li, the Numerical Algorithms Group (NAG) +! Copyright (C) 2021 the University of Edinburgh (UoE) +! +!======================================================================= + +! This module provides parallel IO facilities for applications based on +! 2D decomposition. + +module decomp_2d_io + + use decomp_2d + use MPI +#ifdef T3PIO + use t3pio +#endif + +#ifdef ADIOS2 + use adios2 +#endif + + implicit none + + private ! Make everything private unless declared public + + public :: decomp_2d_write_one, decomp_2d_read_one, & + decomp_2d_write_var, decomp_2d_read_var, & + decomp_2d_write_scalar, decomp_2d_read_scalar, & + decomp_2d_write_plane, decomp_2d_write_every, & + decomp_2d_write_subdomain, & + decomp_2d_write_outflow, decomp_2d_read_inflow +#ifdef ADIOS2 + public :: adios2_register_variable +#endif + + ! Generic interface to handle multiple data types + + interface decomp_2d_write_one + module procedure write_one_real + module procedure write_one_complex + module procedure mpiio_write_real_coarse + module procedure mpiio_write_real_probe +#ifdef ADIOS2 + module procedure adios2_write_real_coarse +#endif + end interface decomp_2d_write_one + + interface decomp_2d_read_one + module procedure read_one_real + module procedure read_one_complex + end interface decomp_2d_read_one + + interface decomp_2d_write_var + module procedure write_var_real + module procedure write_var_complex + end interface decomp_2d_write_var + + interface decomp_2d_read_var + module procedure read_var_real + module procedure read_var_complex + end interface decomp_2d_read_var + + interface decomp_2d_write_scalar + module procedure write_scalar_real + module procedure write_scalar_complex + module procedure write_scalar_integer + module procedure write_scalar_logical + end interface decomp_2d_write_scalar + + interface decomp_2d_read_scalar + module procedure read_scalar_real + module procedure read_scalar_complex + module procedure read_scalar_integer + module procedure read_scalar_logical + end interface decomp_2d_read_scalar + + interface decomp_2d_write_plane + module procedure write_plane_3d_real + module procedure write_plane_3d_complex + ! module procedure write_plane_2d + end interface decomp_2d_write_plane + + interface decomp_2d_write_every + module procedure write_every_real + module procedure write_every_complex + end interface decomp_2d_write_every + + interface decomp_2d_write_subdomain + module procedure write_subdomain + end interface decomp_2d_write_subdomain + + interface decomp_2d_write_outflow + module procedure write_outflow + end interface decomp_2d_write_outflow + + interface decomp_2d_read_inflow + module procedure read_inflow + end interface decomp_2d_read_inflow + +contains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Using MPI-IO library to write a single 3D array to a file +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_one_real(ipencil,var,filename,opt_decomp) + + implicit none + + integer, intent(IN) :: ipencil + real(mytype), dimension(:,:,:), intent(IN) :: var + character(len=*), intent(IN) :: filename + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer(kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, fh, data_type, info, gs + + data_type = real_type + +#include "io_write_one.inc" + + return + end subroutine write_one_real + + + subroutine write_one_complex(ipencil,var,filename,opt_decomp) + + implicit none + + integer, intent(IN) :: ipencil + complex(mytype), dimension(:,:,:), intent(IN) :: var + character(len=*), intent(IN) :: filename + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer(kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, fh, data_type, info, gs + + data_type = complex_type + +#include "io_write_one.inc" + + return + end subroutine write_one_complex + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Using MPI-IO library to read from a file a single 3D array +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine read_one_real(ipencil,var,filename,opt_decomp) + + implicit none + + integer, intent(IN) :: ipencil + real(mytype), dimension(:,:,:), intent(INOUT) :: var + character(len=*), intent(IN) :: filename + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer(kind=MPI_OFFSET_KIND) :: disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, fh, data_type + real(mytype_single), allocatable, dimension(:,:,:) :: varsingle + + data_type = real_type_single + allocate (varsingle(xstV(1):xenV(1),xstV(2):xenV(2),xstV(3):xenV(3))) + + if (present(opt_decomp)) then + decomp = opt_decomp + else + call get_decomp_info(decomp) + end if + + ! determine subarray parameters + sizes(1) = decomp%xsz(1) + sizes(2) = decomp%ysz(2) + sizes(3) = decomp%zsz(3) + + if (ipencil == 1) then + subsizes(1) = decomp%xsz(1) + subsizes(2) = decomp%xsz(2) + subsizes(3) = decomp%xsz(3) + starts(1) = decomp%xst(1)-1 ! 0-based index + starts(2) = decomp%xst(2)-1 + starts(3) = decomp%xst(3)-1 + else if (ipencil == 2) then + subsizes(1) = decomp%ysz(1) + subsizes(2) = decomp%ysz(2) + subsizes(3) = decomp%ysz(3) + starts(1) = decomp%yst(1)-1 + starts(2) = decomp%yst(2)-1 + starts(3) = decomp%yst(3)-1 + else if (ipencil == 3) then + subsizes(1) = decomp%zsz(1) + subsizes(2) = decomp%zsz(2) + subsizes(3) = decomp%zsz(3) + starts(1) = decomp%zst(1)-1 + starts(2) = decomp%zst(2)-1 + starts(3) = decomp%zst(3)-1 + endif + + call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts, & + MPI_ORDER_FORTRAN, data_type, newtype, ierror) + call MPI_TYPE_COMMIT(newtype,ierror) + call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, & + MPI_MODE_RDONLY, MPI_INFO_NULL, & + fh, ierror) + disp = 0_MPI_OFFSET_KIND + call MPI_FILE_SET_VIEW(fh,disp,data_type, & + newtype,'native',MPI_INFO_NULL,ierror) + call MPI_FILE_READ_ALL(fh, varsingle, & + subsizes(1)*subsizes(2)*subsizes(3), & + data_type, MPI_STATUS_IGNORE, ierror) + call MPI_FILE_CLOSE(fh,ierror) + call MPI_TYPE_FREE(newtype,ierror) + var = real(varsingle,mytype) + deallocate(varsingle) + + return + end subroutine read_one_real + + + subroutine read_one_complex(ipencil,var,filename,opt_decomp) + + implicit none + + integer, intent(IN) :: ipencil + complex(mytype), dimension(:,:,:), intent(INOUT) :: var + character(len=*), intent(IN) :: filename + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer(kind=MPI_OFFSET_KIND) :: disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, fh, data_type + + data_type = complex_type + +#include "io_read_one.inc" + + return + end subroutine read_one_complex + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Write a 3D array as part of a big MPI-IO file, starting from + ! displacement 'disp'; 'disp' will be updated after the writing + ! operation to prepare the writing of next chunk of data. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_var_real(fh,disp,ipencil,var,opt_decomp) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: ipencil + real(mytype), dimension(:,:,:), intent(IN) :: var + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, data_type + + data_type = real_type + +#include "io_write_var.inc" + + return + end subroutine write_var_real + + + subroutine write_var_complex(fh,disp,ipencil,var,opt_decomp) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: ipencil + complex(mytype), dimension(:,:,:), intent(IN) :: var + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, data_type + + data_type = complex_type + +#include "io_write_var.inc" + + return + end subroutine write_var_complex + + + subroutine write_outflow(fh,disp,ntimesteps,var,opt_decomp) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: ntimesteps + real(mytype), dimension(:,:,:), intent(IN) :: var + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, data_type + + data_type = real_type + +#include "io_write_outflow.f90" + + return + end subroutine write_outflow + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Read a 3D array as part of a big MPI-IO file, starting from + ! displacement 'disp'; 'disp' will be updated after the reading + ! operation to prepare the reading of next chunk of data. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine read_var_real(fh,disp,ipencil,var,opt_decomp) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: ipencil + real(mytype), dimension(:,:,:), intent(INOUT) :: var + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, data_type + + data_type = real_type + +#include "io_read_var.inc" + + return + end subroutine read_var_real + + + subroutine read_var_complex(fh,disp,ipencil,var,opt_decomp) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: ipencil + complex(mytype), dimension(:,:,:), intent(INOUT) :: var + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, data_type + + data_type = complex_type + +#include "io_read_var.inc" + + return + end subroutine read_var_complex + + + subroutine read_inflow(fh,disp,ntimesteps,var,opt_decomp) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: ntimesteps + real(mytype), dimension(:,:,:), intent(INOUT) :: var + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + TYPE(DECOMP_INFO) :: decomp + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror, newtype, data_type + + data_type = real_type + +#include "io_read_inflow.f90" + + return + end subroutine read_inflow + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Write scalar variables as part of a big MPI-IO file, starting from + ! displacement 'disp'; 'disp' will be updated after the reading + ! operation to prepare the reading of next chunk of data. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_scalar_real(fh,disp,n,var) + + implicit none + + integer, intent(IN) :: fh ! file handle + integer(KIND=MPI_OFFSET_KIND), & + intent(INOUT) :: disp ! displacement + integer, intent(IN) :: n ! number of scalars + real(mytype), dimension(n), & + intent(IN) :: var ! array of scalars + + integer :: m, ierror + + call MPI_FILE_SET_VIEW(fh,disp,real_type, & + real_type,'native',MPI_INFO_NULL,ierror) + if (nrank==0) then + m = n ! only one rank needs to write + else + m = 0 + end if + call MPI_FILE_WRITE_ALL(fh, var, m, real_type, & + MPI_STATUS_IGNORE, ierror) + disp = disp + n*mytype_bytes + + return + end subroutine write_scalar_real + + + subroutine write_scalar_complex(fh,disp,n,var) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: n + complex(mytype), dimension(n), intent(IN) :: var + + integer :: m, ierror + + call MPI_FILE_SET_VIEW(fh,disp,complex_type, & + complex_type,'native',MPI_INFO_NULL,ierror) + if (nrank==0) then + m = n + else + m = 0 + end if + call MPI_FILE_WRITE_ALL(fh, var, m, complex_type, & + MPI_STATUS_IGNORE, ierror) + disp = disp + n*mytype_bytes*2 + + return + end subroutine write_scalar_complex + + + subroutine write_scalar_integer(fh,disp,n,var) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: n + integer, dimension(n), intent(IN) :: var + + integer :: m, ierror + + call MPI_FILE_SET_VIEW(fh,disp,MPI_INTEGER, & + MPI_INTEGER,'native',MPI_INFO_NULL,ierror) + if (nrank==0) then + m = n + else + m = 0 + end if + call MPI_FILE_WRITE_ALL(fh, var, m, MPI_INTEGER, & + MPI_STATUS_IGNORE, ierror) + call MPI_TYPE_SIZE(MPI_INTEGER,m,ierror) + disp = disp + n*m + + return + end subroutine write_scalar_integer + + + subroutine write_scalar_logical(fh,disp,n,var) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: n + logical, dimension(n), intent(IN) :: var + + integer :: m, ierror + + call MPI_FILE_SET_VIEW(fh,disp,MPI_LOGICAL, & + MPI_LOGICAL,'native',MPI_INFO_NULL,ierror) + if (nrank==0) then + m = n + else + m = 0 + end if + call MPI_FILE_WRITE_ALL(fh, var, m, MPI_LOGICAL, & + MPI_STATUS_IGNORE, ierror) + call MPI_TYPE_SIZE(MPI_LOGICAL,m,ierror) + disp = disp + n*m + + return + end subroutine write_scalar_logical + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Read scalar variables as part of a big MPI-IO file, starting from + ! displacement 'disp'; 'disp' will be updated after the reading + ! operation to prepare the reading of next chunk of data. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine read_scalar_real(fh,disp,n,var) + + implicit none + + integer, intent(IN) :: fh ! file handle + integer(KIND=MPI_OFFSET_KIND), & + intent(INOUT) :: disp ! displacement + integer, intent(IN) :: n ! number of scalars + real(mytype), dimension(n), & + intent(INOUT) :: var ! array of scalars + + integer :: ierror + + call MPI_FILE_SET_VIEW(fh,disp,real_type, & + real_type,'native',MPI_INFO_NULL,ierror) + call MPI_FILE_READ_ALL(fh, var, n, real_type, & + MPI_STATUS_IGNORE, ierror) + disp = disp + n*mytype_bytes + + return + end subroutine read_scalar_real + + + subroutine read_scalar_complex(fh,disp,n,var) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: n + complex(mytype), dimension(n), intent(INOUT) :: var + + integer :: ierror + + call MPI_FILE_SET_VIEW(fh,disp,complex_type, & + complex_type,'native',MPI_INFO_NULL,ierror) + call MPI_FILE_READ_ALL(fh, var, n, complex_type, & + MPI_STATUS_IGNORE, ierror) + disp = disp + n*mytype_bytes*2 + + return + end subroutine read_scalar_complex + + + subroutine read_scalar_integer(fh,disp,n,var) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: n + integer, dimension(n), intent(INOUT) :: var + + integer :: m, ierror + + call MPI_FILE_SET_VIEW(fh,disp,MPI_INTEGER, & + MPI_INTEGER,'native',MPI_INFO_NULL,ierror) + call MPI_FILE_READ_ALL(fh, var, n, MPI_INTEGER, & + MPI_STATUS_IGNORE, ierror) + call MPI_TYPE_SIZE(MPI_INTEGER,m,ierror) + disp = disp + n*m + + return + end subroutine read_scalar_integer + + + subroutine read_scalar_logical(fh,disp,n,var) + + implicit none + + integer, intent(IN) :: fh + integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp + integer, intent(IN) :: n + logical, dimension(n), intent(INOUT) :: var + + integer :: m, ierror + + call MPI_FILE_SET_VIEW(fh,disp,MPI_LOGICAL, & + MPI_LOGICAL,'native',MPI_INFO_NULL,ierror) + call MPI_FILE_READ_ALL(fh, var, n, MPI_LOGICAL, & + MPI_STATUS_IGNORE, ierror) + call MPI_TYPE_SIZE(MPI_LOGICAL,m,ierror) + disp = disp + n*m + + return + end subroutine read_scalar_logical + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Write a 2D slice of the 3D data to a file +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_plane_3d_real(ipencil,var,iplane,n,filename, & + opt_decomp) + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + real(mytype), dimension(:,:,:), intent(IN) :: var + integer, intent(IN) :: iplane !(x-plane=1; y-plane=2; z-plane=3) + integer, intent(IN) :: n ! which plane to write (global coordinate) + character(len=*), intent(IN) :: filename + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + real(mytype), allocatable, dimension(:,:,:) :: wk, wk2 + real(mytype), allocatable, dimension(:,:,:) :: wk2d + TYPE(DECOMP_INFO) :: decomp + integer(kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: i,j,k, ierror, newtype, fh, data_type + + data_type = real_type + +#include "io_write_plane.inc" + + return + end subroutine write_plane_3d_real + + + subroutine write_plane_3d_complex(ipencil,var,iplane,n, & + filename,opt_decomp) + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + complex(mytype), dimension(:,:,:), intent(IN) :: var + integer, intent(IN) :: iplane !(x-plane=1; y-plane=2; z-plane=3) + integer, intent(IN) :: n ! which plane to write (global coordinate) + character(len=*), intent(IN) :: filename + TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp + + complex(mytype), allocatable, dimension(:,:,:) :: wk, wk2 + complex(mytype), allocatable, dimension(:,:,:) :: wk2d + TYPE(DECOMP_INFO) :: decomp + integer(kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: i,j,k, ierror, newtype, fh, data_type + + data_type = complex_type + +#include "io_write_plane.inc" + + return + end subroutine write_plane_3d_complex + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Write a 2D array to a file +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !************** TO DO *************** + !* Consider handling distributed 2D data set + ! subroutine write_plane_2d(ipencil,var,filename) + ! integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + ! real(mytype), dimension(:,:), intent(IN) :: var ! 2D array + ! character(len=*), intent(IN) :: filename + ! + ! if (ipencil==1) then + ! ! var should be defined as var(xsize(2) + ! + ! else if (ipencil==2) then + ! + ! else if (ipencil==3) then + ! + ! end if + ! + ! return + ! end subroutine write_plane_2d + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Write 3D array data for every specified mesh point +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_every_real(ipencil,var,iskip,jskip,kskip, & + filename, from1) + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + real(mytype), dimension(:,:,:), intent(IN) :: var + integer, intent(IN) :: iskip,jskip,kskip + character(len=*), intent(IN) :: filename + logical, intent(IN) :: from1 ! .true. - save 1,n+1,2n+1... + ! .false. - save n,2n,3n... + + real(mytype), allocatable, dimension(:,:,:) :: wk, wk2 + integer(kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: i,j,k, ierror, newtype, fh, key,color,newcomm, data_type + integer, dimension(3) :: xsz,ysz,zsz,xst,yst,zst,xen,yen,zen,skip + + data_type = real_type + +#include "io_write_every.inc" + + return + end subroutine write_every_real + + + subroutine write_every_complex(ipencil,var,iskip,jskip,kskip, & + filename, from1) + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + complex(mytype), dimension(:,:,:), intent(IN) :: var + integer, intent(IN) :: iskip,jskip,kskip + character(len=*), intent(IN) :: filename + logical, intent(IN) :: from1 ! .true. - save 1,n+1,2n+1... + ! .false. - save n,2n,3n... + + complex(mytype), allocatable, dimension(:,:,:) :: wk, wk2 + integer(kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: i,j,k, ierror, newtype, fh, key,color,newcomm, data_type + integer, dimension(3) :: xsz,ysz,zsz,xst,yst,zst,xen,yen,zen,skip + + data_type = complex_type + +#include "io_write_every.inc" + + return + end subroutine write_every_complex + + subroutine coarse_extents(ipencil, icoarse, sizes, subsizes, starts) + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + integer, intent(IN) :: icoarse !(nstat=1; nvisu=2) + + integer, dimension(3) :: sizes, subsizes, starts + integer :: ierror + + if ((icoarse.lt.1).or.(icoarse.gt.2)) then + print *, "Error invalid value of icoarse: ", icoarse + call MPI_ABORT(MPI_COMM_WORLD, -1, ierror) + endif + if ((ipencil.lt.1).or.(ipencil.gt.3)) then + print *, "Error invalid value of ipencil: ", ipencil + call MPI_ABORT(MPI_COMM_WORLD, -1, ierror) + endif + + if (icoarse==1) then + sizes(1) = xszS(1) + sizes(2) = yszS(2) + sizes(3) = zszS(3) + + if (ipencil == 1) then + subsizes(1) = xszS(1) + subsizes(2) = xszS(2) + subsizes(3) = xszS(3) + starts(1) = xstS(1)-1 ! 0-based index + starts(2) = xstS(2)-1 + starts(3) = xstS(3)-1 + else if (ipencil == 2) then + subsizes(1) = yszS(1) + subsizes(2) = yszS(2) + subsizes(3) = yszS(3) + starts(1) = ystS(1)-1 + starts(2) = ystS(2)-1 + starts(3) = ystS(3)-1 + else if (ipencil == 3) then + subsizes(1) = zszS(1) + subsizes(2) = zszS(2) + subsizes(3) = zszS(3) + starts(1) = zstS(1)-1 + starts(2) = zstS(2)-1 + starts(3) = zstS(3)-1 + endif + elseif (icoarse==2) then + sizes(1) = xszV(1) + sizes(2) = yszV(2) + sizes(3) = zszV(3) + + if (ipencil == 1) then + subsizes(1) = xszV(1) + subsizes(2) = xszV(2) + subsizes(3) = xszV(3) + starts(1) = xstV(1)-1 ! 0-based index + starts(2) = xstV(2)-1 + starts(3) = xstV(3)-1 + else if (ipencil == 2) then + subsizes(1) = yszV(1) + subsizes(2) = yszV(2) + subsizes(3) = yszV(3) + starts(1) = ystV(1)-1 + starts(2) = ystV(2)-1 + starts(3) = ystV(3)-1 + else if (ipencil == 3) then + subsizes(1) = zszV(1) + subsizes(2) = zszV(2) + subsizes(3) = zszV(3) + starts(1) = zstV(1)-1 + starts(2) = zstV(2)-1 + starts(3) = zstV(3)-1 + endif + endif + + end subroutine coarse_extents + + subroutine mpiio_write_real_coarse(ipencil,var,filename,icoarse) + + ! USE param + ! USE variables + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + integer, intent(IN) :: icoarse !(nstat=1; nvisu=2) + real(mytype), dimension(:,:,:), intent(IN) :: var + real(mytype_single), allocatable, dimension(:,:,:) :: varsingle + character(len=*) :: filename + + integer (kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: i,j,k, ierror, newtype, fh + + call coarse_extents(ipencil, icoarse, sizes, subsizes, starts) + allocate (varsingle(xstV(1):xenV(1),xstV(2):xenV(2),xstV(3):xenV(3))) + varsingle=real(var, mytype_single) + + call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts, & + MPI_ORDER_FORTRAN, real_type_single, newtype, ierror) + call MPI_TYPE_COMMIT(newtype,ierror) + call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, & + MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, & + fh, ierror) + filesize = 0_MPI_OFFSET_KIND + call MPI_FILE_SET_SIZE(fh,filesize,ierror) ! guarantee overwriting + disp = 0_MPI_OFFSET_KIND + call MPI_FILE_SET_VIEW(fh,disp,real_type_single, & + newtype,'native',MPI_INFO_NULL,ierror) + call MPI_FILE_WRITE_ALL(fh, varsingle, & + subsizes(1)*subsizes(2)*subsizes(3), & + real_type_single, MPI_STATUS_IGNORE, ierror) + call MPI_FILE_CLOSE(fh,ierror) + call MPI_TYPE_FREE(newtype,ierror) + + deallocate(varsingle) + + return + end subroutine mpiio_write_real_coarse + +#ifdef ADIOS2 + subroutine adios2_register_variable(io, varname, ipencil, icoarse) + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + integer, intent(IN) :: icoarse !(nstat=1; nvisu=2) + type(adios2_io), intent(in) :: io + + character*(*), intent(in) :: varname + + integer, dimension(3) :: sizes, subsizes, starts + type(adios2_variable) :: var_handle + integer, parameter :: ndims = 3 + logical, parameter :: adios2_constant_dims = .true. + integer :: data_type + integer :: ierror + + call coarse_extents(ipencil, icoarse, sizes, subsizes, starts) + + ! Check if variable already exists, if not create it + call adios2_inquire_variable(var_handle, io, varname, ierror) + if (.not.var_handle % valid) then + !! New variable + if(nrank==0) print *, "Registering variable for IO: ", varname + + ! Need to set the ADIOS2 data type + if (mytype_single.eq.kind(0.0d0)) then + !! Double + data_type = adios2_type_dp + else if (mytype_single.eq.kind(0.0)) then + !! Single + data_type = adios2_type_real + else + print *, "Trying to write unknown data type!" + call MPI_ABORT(MPI_COMM_WORLD, -1, ierror) + endif + + call adios2_define_variable(var_handle, io, varname, data_type, & + ndims, int(sizes, kind=8), int(starts, kind=8), int(subsizes, kind=8), & + adios2_constant_dims, ierror) + endif + + end subroutine adios2_register_variable + subroutine adios2_write_real_coarse(ipencil,var,varname,icoarse,adios,engine,io) + + ! USE param + ! USE variables + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + integer, intent(IN) :: icoarse !(nstat=1; nvisu=2) + real(mytype), dimension(:,:,:), intent(IN) :: var + type(adios2_adios), intent(in) :: adios + type(adios2_engine), intent(in) :: engine + type(adios2_io), intent(in) :: io + character*(*), intent(in) :: varname + + integer (kind=MPI_OFFSET_KIND) :: filesize, disp + integer :: i,j,k, ierror, newtype, fh + type(adios2_variable) :: var_handle + + call adios2_inquire_variable(var_handle, io, varname, ierror) + if (.not.var_handle % valid) then + call adios2_register_variable(io, varname, ipencil, icoarse) + endif + + call adios2_put(engine, var_handle, var, adios2_mode_deferred, ierror) + + return + end subroutine adios2_write_real_coarse +#endif + + subroutine mpiio_write_real_probe(ipencil,var,filename,nlength) + + ! USE param + ! USE variables + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + integer, intent(in) :: nlength + real(mytype), dimension(:,:,:,:), intent(IN) :: var + + character(len=*) :: filename + + integer (kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(4) :: sizes, subsizes, starts + integer :: i,j,k, ierror, newtype, fh + + sizes(1) = xszP(1) + sizes(2) = yszP(2) + sizes(3) = zszP(3) + sizes(4) = nlength + if (ipencil == 1) then + subsizes(1) = xszP(1) + subsizes(2) = xszP(2) + subsizes(3) = xszP(3) + subsizes(4) = nlength + starts(1) = xstP(1)-1 ! 0-based index + starts(2) = xstP(2)-1 + starts(3) = xstP(3)-1 + starts(4) = 0 + else if (ipencil == 2) then + subsizes(1) = yszP(1) + subsizes(2) = yszP(2) + subsizes(3) = yszP(3) + starts(1) = ystP(1)-1 + starts(2) = ystP(2)-1 + starts(3) = ystP(3)-1 + else if (ipencil == 3) then + subsizes(1) = zszP(1) + subsizes(2) = zszP(2) + subsizes(3) = zszP(3) + starts(1) = zstP(1)-1 + starts(2) = zstP(2)-1 + starts(3) = zstP(3)-1 + endif + ! print *,nrank,starts(1),starts(2),starts(3),starts(4) + call MPI_TYPE_CREATE_SUBARRAY(4, sizes, subsizes, starts, & + MPI_ORDER_FORTRAN, real_type, newtype, ierror) + call MPI_TYPE_COMMIT(newtype,ierror) + call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, & + MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, & + fh, ierror) + filesize = 0_MPI_OFFSET_KIND + call MPI_FILE_SET_SIZE(fh,filesize,ierror) ! guarantee overwriting + disp = 0_MPI_OFFSET_KIND + call MPI_FILE_SET_VIEW(fh,disp,real_type, & + newtype,'native',MPI_INFO_NULL,ierror) + call MPI_FILE_WRITE_ALL(fh, var, & + subsizes(1)*subsizes(2)*subsizes(3)*subsizes(4), & + real_type, MPI_STATUS_IGNORE, ierror) + call MPI_FILE_CLOSE(fh,ierror) + call MPI_TYPE_FREE(newtype,ierror) + + + return + end subroutine mpiio_write_real_probe + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Write a 3D data set covering a smaller sub-domain only +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_subdomain(ipencil,var,is,ie,js,je,ks,ke,filename) + + implicit none + + integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3) + real(mytype), dimension(:,:,:), intent(IN) :: var + integer, intent(IN) :: is, ie, js, je, ks, ke + character(len=*), intent(IN) :: filename + + real(mytype), allocatable, dimension(:,:,:) :: wk, wk2 + integer(kind=MPI_OFFSET_KIND) :: filesize, disp + integer, dimension(3) :: sizes, subsizes, starts + integer :: color, key, errorcode, newcomm, ierror + integer :: newtype, fh, data_type, i, j, k + integer :: i1, i2, j1, j2, k1, k2 + + data_type = real_type + + ! validate the input paramters + if (is<1 .OR. ie>nx_global .OR. js<1 .OR. je>ny_global .OR. & + ks<1 .OR. ke>nz_global) then + errorcode = 10 + call decomp_2d_abort(errorcode, & + 'Invalid subdomain specified in I/O') + end if + + ! create a communicator for all those MPI ranks containing the subdomain + color = 1 + key = 1 + if (ipencil==1) then + if (xstart(1)>ie .OR. xend(1)je .OR. xend(2)ke .OR. xend(3)ie .OR. yend(1)je .OR. yend(2)ke .OR. yend(3)ie .OR. zend(1)je .OR. zend(2)ke .OR. zend(3)ie .AND. xstart(1)ie) then + subsizes(1) = ie - xstart(1) + 1 + end if + subsizes(2) = xsize(2) + starts(2) = xstart(2) - js + if (xend(2)>je .AND. xstart(2)je) then + subsizes(2) = je - xstart(2) + 1 + end if + subsizes(3) = xsize(3) + starts(3) = xstart(3) - ks + if (xend(3)>ke .AND. xstart(3)ke) then + subsizes(3) = ke - xstart(3) + 1 + end if + + else if (ipencil==2) then + + ! TODO + + else if (ipencil==3) then + + ! TODO + + end if + + + ! copy data from orginal to a temp array + ! pay attention to blocks only partially cover the sub-domain + if (ipencil==1) then + + if (xend(1)>ie .AND. xstart(1)ie) then + i1 = xstart(1) + i2 = ie + else if (xstart(1)je .AND. xstart(2)je) then + j1 = xstart(2) + j2 = je + else if (xstart(2)ke .AND. xstart(3)ke) then + k1 = xstart(3) + k2 = ke + else if (xstart(3) /dev/null") + end if + end if + + ! HDD usage of visu module + if (nrank==0) then + noutput = (ilast - ifirst +1)/ioutput + + nsnapout = 4 + if (iscalar.ne.0) nsnapout = nsnapout + numscalar + + memout = prec * nsnapout * noutput + if (output2D.eq.0) then + memout = memout * xszV(1) * yszV(2) * zszV(3) + else if (output2D.eq.1) then + memout = memout * yszV(2) * zszV(3) + else if (output2D.eq.2) then + memout = memout * xszV(1) * zszV(3) + else if (output2D.eq.3) then + memout = memout * xszV(1) * yszV(2) + endif + write(*,*)'===========================================================' + write(*,*)'Visu module requires ',real(memout*1e-9,4),'GB' + write(*,*)'===========================================================' + end if + + ! Safety check + if (output2D < 0 .or. output2D > 3 & + .or. (output2d == 2.and.istret /= 0)) then + if (nrank.eq.0) write(*,*) "Visu module: incorrect value for output2D." + call MPI_ABORT(MPI_COMM_WORLD, 0, noutput) + stop + endif + +#ifdef ADIOS2 + !! TODO: make this a runtime-option + adios2_debug_mode = .true. + + call adios2_init(adios, trim(config_file), MPI_COMM_WORLD, adios2_debug_mode, ierror) + if (ierror.ne.0) then + print *, "Error initialising ADIOS2 - is adios2_config.xml present and valid?" + call MPI_ABORT(MPI_COMM_WORLD, -1, ierror) + endif + call adios2_declare_io(io_write_real_coarse, adios, "solution-io", ierror) + if (io_write_real_coarse % engine_type.eq."BP4") then + write(outfile, *) "data.bp4" + else if (io_write_real_coarse % engine_type.eq."HDF5") then + write(outfile, *) "data.hdf5" + else + print *, "Unknown engine!" + call MPI_ABORT(MPI_COMM_WORLD, -1, ierror) + endif + + !! Register variables + call adios2_register_variable(io_write_real_coarse, "ux", 1, 2) + call adios2_register_variable(io_write_real_coarse, "uy", 1, 2) + call adios2_register_variable(io_write_real_coarse, "uz", 1, 2) + call adios2_register_variable(io_write_real_coarse, "pp", 1, 2) + ! if (ilmn) then + ! call adios2_register_variable(io_write_real_coarse, "rho", 1, 2) + ! endif + if (iscalar.ne.0) then + do is = 1, numscalar + call adios2_register_variable(io_write_real_coarse, "phi"//char(48+is), 1, 2) + enddo + endif + + call adios2_open(engine_write_real_coarse, io_write_real_coarse, trim(outfile), adios2_mode_write, ierror) +#endif + + end subroutine visu_init + + ! + ! Finalise the visu module + ! - Currently only needed to clean up ADIOS2 + ! + subroutine visu_finalise() + +#ifdef ADIOS2 + use adios2 +#endif + implicit none + +#ifdef ADIOS2 + integer :: ierr +#endif + +#ifdef ADIOS2 + call adios2_close(engine_write_real_coarse, ierr) + call adios2_finalize(adios, ierr) +#endif + + end subroutine visu_finalise + + ! + ! Write a snapshot + ! + ! subroutine write_snapshot(rho1, ux1, uy1, uz1, pp3, phi1, ep1, itime, num) + subroutine write_snapshot(ux1, uy1, uz1, ep1, itime, num) + + use decomp_2d, only : transpose_z_to_y, transpose_y_to_x + use decomp_2d, only : mytype, xsize, ysize, zsize + use decomp_2d, only : nrank + + use param, only : iscalar, ioutput + + ! use variables, only : sx, cifip6, cisip6, ciwip6, cifx6, cisx6, ciwx6 + ! use variables, only : sy, cifip6y, cisip6y, ciwip6y, cify6, cisy6, ciwy6 + ! use variables, only : sz, cifip6z, cisip6z, ciwip6z, cifz6, cisz6, ciwz6 + use variables, only : numscalar + + ! use var, only : pp1, ta1, di1, nxmsize + ! use var, only : pp2, ppi2, dip2, ph2, nymsize + ! use var, only : ppi3, dip3, ph3, nzmsize + ! use var, only : npress + + ! use tools, only : rescale_pressure + + implicit none + + !! inputs + real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ux1, uy1, uz1 + real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ep1 + ! real(mytype), dimension(xsize(1), xsize(2), xsize(3), nrhotime), intent(in) :: rho1 + ! real(mytype), dimension(ph3%zst(1):ph3%zen(1),ph3%zst(2):ph3%zen(2),nzmsize,npress), intent(in) :: pp3 + ! real(mytype), dimension(xsize(1), xsize(2), xsize(3), numscalar), intent(in) :: phi1 + integer, intent(in) :: itime + character(len=32), intent(out) :: num + + ! Local variables + integer :: is + integer :: ierr + character(len=30) :: scname + + ! Update log file + if (nrank.eq.0) then + call cpu_time(tstart) + print *,'Writing snapshots =>',itime/ioutput + end if +#ifdef ADIOS2 + call adios2_begin_step(engine_write_real_coarse, adios2_step_mode_append, ierr) +#endif + + ! Snapshot number +#ifndef ADIOS2 + if (filenamedigits) then + ! New enumeration system, it works integrated with xcompact3d_toolbox + write(num, ifilenameformat) itime + else + ! Classic enumeration system + write(num, ifilenameformat) itime/ioutput + endif +#else + ! ADIOS2 is zero-indexed + write(num, '(I0)') itime/ioutput - 1 +#endif + + ! Write XDMF header + if (use_xdmf) call write_xdmf_header(".", "snapshot", trim(num)) + + ! Write velocity + call write_field(ux1, ".", "ux", trim(num)) + call write_field(uy1, ".", "uy", trim(num)) + call write_field(uz1, ".", "uz", trim(num)) + + ! ! Interpolate pressure + ! !WORK Z-PENCILS + ! call interzpv(ppi3,pp3(:,:,:,1),dip3,sz,cifip6z,cisip6z,ciwip6z,cifz6,cisz6,ciwz6,& + ! (ph3%zen(1)-ph3%zst(1)+1),(ph3%zen(2)-ph3%zst(2)+1),nzmsize,zsize(3),1) + ! !WORK Y-PENCILS + ! call transpose_z_to_y(ppi3,pp2,ph3) !nxm nym nz + ! call interypv(ppi2,pp2,dip2,sy,cifip6y,cisip6y,ciwip6y,cify6,cisy6,ciwy6,& + ! (ph3%yen(1)-ph3%yst(1)+1),nymsize,ysize(2),ysize(3),1) + ! !WORK X-PENCILS + ! call transpose_y_to_x(ppi2,pp1,ph2) !nxm ny nz + ! call interxpv(ta1,pp1,di1,sx,cifip6,cisip6,ciwip6,cifx6,cisx6,ciwx6,& + ! nxmsize,xsize(1),xsize(2),xsize(3),1) + + + ! ! Rescale pressure + ! call rescale_pressure(ta1) + + ! ! Write pressure + ! call write_field(ta1, ".", "pp", trim(num), .true.) + + ! ! LMN - write density + ! if (ilmn) call write_field(rho1(:,:,:,1), ".", "rho", trim(num)) + + ! ! Write scalars + ! if (iscalar.ne.0) then + ! do is = 1, numscalar + ! write(scname,"('phi',I2.2)") is + ! call write_field(phi1(:,:,:,is), ".", trim(scname), trim(num), .true.) + ! enddo + ! endif + + end subroutine write_snapshot + + subroutine end_snapshot(itime, num) + + use decomp_2d, only : nrank + use param, only : istret, xlx, yly, zlz + use variables, only : nx, ny, nz, beta + + implicit none + + integer, intent(in) :: itime + character(len=32), intent(in) :: num + + character(len=32) :: fmt2, fmt3, fmt4 + integer :: is + integer :: ierr + + ! Write XDMF footer + if (use_xdmf) call write_xdmf_footer() + + ! Add metadata if XDMF is not used + if (.not. use_xdmf) then + if (nrank==0) then + + write(fmt2,'("(A,I16)")') + write(fmt3,'("(A,F16.4)")') + write(fmt4,'("(A,F16.12)")') + + open(newunit=is,file="./data/snap"//trim(num)//".ini",action='write',status='replace') + write(is,'(A)')'[domain]' + write(is,fmt2) 'nx= ',nx + write(is,fmt2) 'ny= ',ny + write(is,fmt2) 'nz= ',nz + write(is,fmt2) 'istret= ',istret + write(is,fmt4) 'beta= ',beta + write(is,fmt3) 'Lx= ',xlx + write(is,fmt3) 'Ly= ',yly + write(is,fmt3) 'Lz= ',zlz + write(is,'(A)')'[time]' + write(is,fmt2) 'itime= ',itime + ! write(is,fmt3) 'dt= ',dt + ! write(is,fmt3) 't = ',t + close(is) + + endif + endif + +#ifdef ADIOS2 + call adios2_end_step(engine_write_real_coarse, ierr) +#endif + + ! Update log file + if (nrank.eq.0) then + call cpu_time(tend) + write(*,'(" Time for writing snapshots (s): ",F12.8)') tend-tstart + endif + + end subroutine end_snapshot + + ! + ! Write the header of the XDMF file + ! Adapted from https://github.com/fschuch/Xcompact3d/blob/master/src/visu.f90 + ! + subroutine write_xdmf_header(pathname, filename, num) + + use variables, only : nvisu, yp + use param, only : dx,dy,dz,istret + use decomp_2d, only : mytype, nrank, xszV, yszV, zszV, ystV + + implicit none + + ! Arguments + character(len=*), intent(in) :: pathname, filename, num + + ! Local variables + integer :: i,k + real(mytype) :: xp(xszV(1)), zp(zszV(3)) + + if (nrank.eq.0) then + OPEN(newunit=ioxdmf,file="./data/"//pathname//"/"//filename//'-'//num//'.xdmf') + + write(ioxdmf,'(A22)')'' + write(ioxdmf,*)'' + write(ioxdmf,*)'' + write(ioxdmf,*)'' + if (istret.ne.0) then + do i=1,xszV(1) + xp(i) = real(i-1,mytype)*dx*nvisu + enddo + do k=1,zszV(3) + zp(k) = real(k-1,mytype)*dz*nvisu + enddo + write(ioxdmf,*)' ' + else if (output2D.eq.1) then + write(ioxdmf,*)' Dimensions="',zszV(3),yszV(2),1,'">' + else if (output2D.eq.2) then + write(ioxdmf,*)' Dimensions="',zszV(3),1,xszV(1),'">' + else if (output2D.eq.3) then + write(ioxdmf,*)' Dimensions="',1,yszV(2),xszV(1),'">' + endif + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + if (output2D.ne.1) then + write(ioxdmf,*)' ' + write(ioxdmf,*)' ',xp(:) + else + write(ioxdmf,*)' ' + write(ioxdmf,*)' ',xp(1) + endif + write(ioxdmf,*)' ' + if (output2D.ne.2) then + write(ioxdmf,*)' ' + write(ioxdmf,*)' ',yp(ystV(1)::nvisu) + else + write(ioxdmf,*)' ' + write(ioxdmf,*)' ',yp(1) + endif + write(ioxdmf,*)' ' + if (output2D.ne.3) then + write(ioxdmf,*)' ' + write(ioxdmf,*)' ',zp(:) + else + write(ioxdmf,*)' ' + write(ioxdmf,*)' ',zp(1) + endif + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + else + write(ioxdmf,*)' ' + else if (output2D.eq.1) then + write(ioxdmf,*)' Dimensions="',zszV(3),yszV(2),1,'">' + else if (output2D.eq.2) then + write(ioxdmf,*)' Dimensions="',zszV(3),1,xszV(1),'">' + else if (output2D.eq.3) then + write(ioxdmf,*)' Dimensions="',1,yszV(2),xszV(1),'">' + endif + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + write(ioxdmf,*)' 0.0 0.0 0.0' + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + if (output2D.eq.0) then + write(ioxdmf,*)' ',nvisu*dz,nvisu*dy,nvisu*dx + else if (output2D.eq.1) then + write(ioxdmf,*)' ',dz,dy,1. + else if (output2D.eq.2) then + write(ioxdmf,*)' ',dz,1.,dx + else if (output2D.eq.3) then + write(ioxdmf,*)' ',1.,dy,dx + endif + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + endif + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + endif + end subroutine write_xdmf_header + + ! + ! Write the footer of the XDMF file + ! Adapted from https://github.com/fschuch/Xcompact3d/blob/master/src/visu.f90 + ! + subroutine write_xdmf_footer() + + use decomp_2d, only : nrank + + implicit none + + if (nrank.eq.0) then + write(ioxdmf,'(/)') + write(ioxdmf,*)' ' + write(ioxdmf,*)'' + write(ioxdmf,'(A7)')'' + close(ioxdmf) + endif + + end subroutine write_xdmf_footer + + ! + ! Write the given field for visualization + ! Adapted from https://github.com/fschuch/Xcompact3d/blob/master/src/visu.f90 + ! + subroutine write_field(f1, pathname, filename, num, skip_ibm) + + use mpi + + use param, only : iibm + use decomp_2d, only : mytype, xsize, xszV, yszV, zszV + use decomp_2d, only : nrank, fine_to_coarseV + use decomp_2d_io, only : decomp_2d_write_one, decomp_2d_write_plane + + implicit none + + real(mytype), intent(in), dimension(xsize(1),xsize(2),xsize(3)) :: f1 + character(len=*), intent(in) :: pathname, filename, num + logical, optional, intent(in) :: skip_ibm + + real(mytype), dimension(xsize(1),xsize(2),xsize(3)) :: uvisu + real(mytype), dimension(xsize(1),xsize(2),xsize(3)) :: local_array + + integer :: ierr + + if (use_xdmf) then + if (nrank.eq.0) then + write(ioxdmf,*)' ' +#ifndef ADIOS2 + write(ioxdmf,*)' ' + else if (output2D.eq.1) then + write(ioxdmf,*)' Dimensions="',zszV(3),yszV(2),1,'">' + else if (output2D.eq.2) then + write(ioxdmf,*)' Dimensions="',zszV(3),1,xszV(1),'">' + else if (output2D.eq.3) then + write(ioxdmf,*)' Dimensions="',1,yszV(2),xszV(1),'">' + endif +#ifndef ADIOS2 + write(ioxdmf,*)' ./'//pathname//"/"//filename//'-'//num//'.bin' +#else + write(ioxdmf,*)' ../data.hdf5:/Step'//num//'/'//filename +#endif + write(ioxdmf,*)' ' + write(ioxdmf,*)' ' + endif + endif + + if (iibm==2 .and. .not.present(skip_ibm)) then + ! local_array(:,:,:) = (one - ep1(:,:,:)) * f1(:,:,:) + local_array(:,:,:) = f1(:,:,:) + else + local_array(:,:,:) = f1(:,:,:) + endif + + if (output2D.eq.0) then +#ifndef ADIOS2 + uvisu = 0.0_mytype + call fine_to_coarseV(1,local_array,uvisu) + call decomp_2d_write_one(1,uvisu,"./data/"//pathname//'/'//filename//'-'//num//'.bin',2) +#else + if (iibm==2 .and. (.not.present(skip_ibm))) then + print *, "Not Implemented: currently ADIOS2 IO doesn't support IBM-blanking" + call MPI_ABORT(MPI_COMM_WORLD, -1, ierr) + endif + call decomp_2d_write_one(1,f1,filename,2,adios,engine_write_real_coarse,io_write_real_coarse) +#endif + else + call decomp_2d_write_plane(1,local_array,output2D,-1,"./data/"//pathname//'/'//filename//'-'//num//'.bin') + endif + + end subroutine write_field + +end module visu diff --git a/src/x3d_tools.f90 b/src/x3d_tools.f90 index edcaadb..81c2a96 100644 --- a/src/x3d_tools.f90 +++ b/src/x3d_tools.f90 @@ -91,6 +91,8 @@ subroutine init_xcompact3d(ndt_max) use var + use visu, only : visu_init + use variables, only : nx, ny, nz, nxm, nym, nzm use variables, only : p_row, p_col use variables, only : test_mode @@ -180,6 +182,8 @@ subroutine init_xcompact3d(ndt_max) call decomp_2d_poisson_init() call decomp_info_init(nxm,nym,nzm,phG) + call visu_init() + endsubroutine init_xcompact3d !######################################################################## !######################################################################## @@ -196,6 +200,8 @@ subroutine finalise_xcompact3d(flag) use x3d_derive, only : x3d_derive_finalize use var, only : var_finalize + use visu, only : visu_finalise + implicit none logical, intent(in) :: flag @@ -214,9 +220,11 @@ subroutine finalise_xcompact3d(flag) call x3d_operator_z_data_finalize() call var_finalize() + call visu_finalise() + call decomp_2d_finalize() if (flag) then CALL MPI_FINALIZE(ierr) endif - + endsubroutine finalise_xcompact3d From 2e2bd84f9ff8d75fe9b21c85d15ac54b2a7a3064 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Mon, 7 Mar 2022 12:53:55 +0000 Subject: [PATCH 2/5] Restore I/O, now outputs but not right freq --- Makefile | 4 +- src/visu.f90 | 98 ++++++++++++++++++++++++++++++++++------------ src/xcompact3d.f90 | 17 +++++++- 3 files changed, 91 insertions(+), 28 deletions(-) diff --git a/Makefile b/Makefile index a156290..aa914e7 100644 --- a/Makefile +++ b/Makefile @@ -20,6 +20,8 @@ LDFLAGS ?= # user can set default linker flags FFLAGS = $(FCFLAGS) LFLAGS = $(LDFLAGS) +NTHREADS ?= 1 # gfortran requires number of threads specified at compile time... + #######CMP settings########### ifeq ($(CMP),intel) FC = mpiifort @@ -38,7 +40,7 @@ else ifeq ($(CMP),gcc) FFLAGS += -ffpe-trap=invalid,zero -fcheck=all -fimplicit-none else FFLAGS += -O3 -march=native - FFLAGS += -fopenmp -ftree-parallelize-loops=12 + FFLAGS += -fopenmp -ftree-parallelize-loops=$(NTHREADS) LFLAGS += -fopenmp endif else ifeq ($(CMP),nagfor) diff --git a/src/visu.f90 b/src/visu.f90 index 16abdc3..efe83a4 100644 --- a/src/visu.f90 +++ b/src/visu.f90 @@ -62,6 +62,12 @@ subroutine visu_init() integer :: is #endif +#ifdef DEBUG + if (nrank == 0) then + print *, "# visu_init start" + endif +#endif + ! Create folder if needed ! XXX: Is this needed for ADIOS2? if (nrank==0) then @@ -71,27 +77,27 @@ subroutine visu_init() end if end if - ! HDD usage of visu module - if (nrank==0) then - noutput = (ilast - ifirst +1)/ioutput - - nsnapout = 4 - if (iscalar.ne.0) nsnapout = nsnapout + numscalar - - memout = prec * nsnapout * noutput - if (output2D.eq.0) then - memout = memout * xszV(1) * yszV(2) * zszV(3) - else if (output2D.eq.1) then - memout = memout * yszV(2) * zszV(3) - else if (output2D.eq.2) then - memout = memout * xszV(1) * zszV(3) - else if (output2D.eq.3) then - memout = memout * xszV(1) * yszV(2) - endif - write(*,*)'===========================================================' - write(*,*)'Visu module requires ',real(memout*1e-9,4),'GB' - write(*,*)'===========================================================' - end if + ! ! HDD usage of visu module + ! if (nrank==0) then + ! noutput = (ilast - ifirst +1)/ioutput + + ! nsnapout = 4 + ! if (iscalar.ne.0) nsnapout = nsnapout + numscalar + + ! memout = prec * nsnapout * noutput + ! if (output2D.eq.0) then + ! memout = memout * xszV(1) * yszV(2) * zszV(3) + ! else if (output2D.eq.1) then + ! memout = memout * yszV(2) * zszV(3) + ! else if (output2D.eq.2) then + ! memout = memout * xszV(1) * zszV(3) + ! else if (output2D.eq.3) then + ! memout = memout * xszV(1) * yszV(2) + ! endif + ! write(*,*)'===========================================================' + ! write(*,*)'Visu module requires ',real(memout*1e-9,4),'GB' + ! write(*,*)'===========================================================' + ! end if ! Safety check if (output2D < 0 .or. output2D > 3 & @@ -137,6 +143,12 @@ subroutine visu_init() call adios2_open(engine_write_real_coarse, io_write_real_coarse, trim(outfile), adios2_mode_write, ierror) #endif +#ifdef DEBUG + if (nrank == 0) then + print *, "# visu_init end" + endif +#endif + end subroutine visu_init ! @@ -165,10 +177,9 @@ end subroutine visu_finalise ! Write a snapshot ! ! subroutine write_snapshot(rho1, ux1, uy1, uz1, pp3, phi1, ep1, itime, num) - subroutine write_snapshot(ux1, uy1, uz1, ep1, itime, num) + subroutine write_snapshot(ux1, uy1, uz1, itime, num) - use decomp_2d, only : transpose_z_to_y, transpose_y_to_x - use decomp_2d, only : mytype, xsize, ysize, zsize + use decomp_2d, only : mytype, xsize use decomp_2d, only : nrank use param, only : iscalar, ioutput @@ -189,7 +200,7 @@ subroutine write_snapshot(ux1, uy1, uz1, ep1, itime, num) !! inputs real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ux1, uy1, uz1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ep1 + ! real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ep1 ! real(mytype), dimension(xsize(1), xsize(2), xsize(3), nrhotime), intent(in) :: rho1 ! real(mytype), dimension(ph3%zst(1):ph3%zen(1),ph3%zst(2):ph3%zen(2),nzmsize,npress), intent(in) :: pp3 ! real(mytype), dimension(xsize(1), xsize(2), xsize(3), numscalar), intent(in) :: phi1 @@ -201,6 +212,12 @@ subroutine write_snapshot(ux1, uy1, uz1, ep1, itime, num) integer :: ierr character(len=30) :: scname +#ifdef DEBUG + if (nrank == 0) then + print *, "# write_snapshot start" + endif +#endif + ! Update log file if (nrank.eq.0) then call cpu_time(tstart) @@ -228,6 +245,7 @@ subroutine write_snapshot(ux1, uy1, uz1, ep1, itime, num) if (use_xdmf) call write_xdmf_header(".", "snapshot", trim(num)) ! Write velocity + print *, "HOLA" call write_field(ux1, ".", "ux", trim(num)) call write_field(uy1, ".", "uy", trim(num)) call write_field(uz1, ".", "uz", trim(num)) @@ -263,6 +281,12 @@ subroutine write_snapshot(ux1, uy1, uz1, ep1, itime, num) ! enddo ! endif +#ifdef DEBUG + if (nrank == 0) then + print *, "# write_snapshot end" + endif +#endif + end subroutine write_snapshot subroutine end_snapshot(itime, num) @@ -341,6 +365,12 @@ subroutine write_xdmf_header(pathname, filename, num) integer :: i,k real(mytype) :: xp(xszV(1)), zp(zszV(3)) +#ifdef DEBUG + if (nrank == 0) then + print *, "## write_xdmf_header start" + end if +#endif + if (nrank.eq.0) then OPEN(newunit=ioxdmf,file="./data/"//pathname//"/"//filename//'-'//num//'.xdmf') @@ -427,6 +457,12 @@ subroutine write_xdmf_header(pathname, filename, num) write(ioxdmf,*)' ' write(ioxdmf,*)' ' endif + +#ifdef DEBUG + if (nrank == 0) then + print *, "## write_xdmf_header end" + end if +#endif end subroutine write_xdmf_header ! @@ -473,6 +509,12 @@ subroutine write_field(f1, pathname, filename, num, skip_ibm) integer :: ierr +#ifdef DEBUG + if (nrank == 0) then + print *, "## write_field -- ", filename, " -- start" + endif +#endif + if (use_xdmf) then if (nrank.eq.0) then write(ioxdmf,*)' ' @@ -536,6 +578,12 @@ subroutine write_field(f1, pathname, filename, num, skip_ibm) call decomp_2d_write_plane(1,local_array,output2D,-1,"./data/"//pathname//'/'//filename//'-'//num//'.bin') endif +#ifdef DEBUG + if (nrank == 0) then + print *, "## write_field -- ", filename, " -- end" + endif +#endif + end subroutine write_field end module visu diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90 index 031fa64..cbdbbcd 100644 --- a/src/xcompact3d.f90 +++ b/src/xcompact3d.f90 @@ -8,12 +8,13 @@ program xcompact3d use var use decomp_2d, only : nrank, xsize, real_type, decomp_2d_warning - use param, only : dt, zero, itr + use param, only : dt, zero, itr, itime, ioutput use transeq, only : calculate_transeq_rhs use navier, only : solve_poisson, cor_vel use case use time_integrators, only : int_time - + use visu, only : write_snapshot + implicit none double precision :: tstart, tend, telapsed, tmin, tmax @@ -21,6 +22,7 @@ program xcompact3d integer :: i, j, k integer :: ndt, ndt_max integer :: code + character(len=32) :: num call boot_xcompact3d() @@ -32,6 +34,12 @@ program xcompact3d tmin = telapsed ndt = 1 + ioutput = 5 + print *, "=== WARNING ===" + print *, " Writing with ioutput=",ioutput,", edit xcompact3d.f90 and recompile to change" + print *, "=== WARNING ===" + + itime = 0 do while(ndt < ndt_max) itr = 1 ! no inner iterations @@ -63,7 +71,12 @@ program xcompact3d print *, "Elapse time min ", tmin, " max ", tmax end if + !! End of timestep + itime = itime + 1 ndt = ndt + 1 + if (mod(itime,ioutput) == 0) then + call write_snapshot(ux1, uy1, uz1, itime, num) + endif call case_postprocess(ux1, uy1, uz1, ndt) end do From 30bbf364e6a21adf82d25ec8d64fcadfa2c1a3b3 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Mon, 7 Mar 2022 13:13:23 +0000 Subject: [PATCH 3/5] Add calls to end_snapshot This "closes" the xdmf file. Also the code was starting at timestep 1 (this should be the value at the END of the first timestep) --- src/xcompact3d.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90 index cbdbbcd..1d3502e 100644 --- a/src/xcompact3d.f90 +++ b/src/xcompact3d.f90 @@ -13,7 +13,7 @@ program xcompact3d use navier, only : solve_poisson, cor_vel use case use time_integrators, only : int_time - use visu, only : write_snapshot + use visu, only : write_snapshot, end_snapshot implicit none @@ -32,7 +32,7 @@ program xcompact3d telapsed = 0 tmin = telapsed - ndt = 1 + ndt = 0 ioutput = 5 print *, "=== WARNING ===" @@ -76,6 +76,7 @@ program xcompact3d ndt = ndt + 1 if (mod(itime,ioutput) == 0) then call write_snapshot(ux1, uy1, uz1, itime, num) + call end_snapshot(itime, num) endif call case_postprocess(ux1, uy1, uz1, ndt) From 3a73e687c275bc6c7b362c81534cae2460ac5d8e Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Mon, 7 Mar 2022 13:38:18 +0000 Subject: [PATCH 4/5] Some changes to visu - I/O now working --- src/visu.f90 | 33 ++++++++++++++++++++++++++++----- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/src/visu.f90 b/src/visu.f90 index efe83a4..70d366d 100644 --- a/src/visu.f90 +++ b/src/visu.f90 @@ -245,7 +245,6 @@ subroutine write_snapshot(ux1, uy1, uz1, itime, num) if (use_xdmf) call write_xdmf_header(".", "snapshot", trim(num)) ! Write velocity - print *, "HOLA" call write_field(ux1, ".", "ux", trim(num)) call write_field(uy1, ".", "uy", trim(num)) call write_field(uz1, ".", "uz", trim(num)) @@ -303,6 +302,12 @@ subroutine end_snapshot(itime, num) character(len=32) :: fmt2, fmt3, fmt4 integer :: is integer :: ierr + +#ifdef DEBUG + if (nrank == 0) then + print *, "# end_snapshot start" + endif +#endif ! Write XDMF footer if (use_xdmf) call write_xdmf_footer() @@ -344,6 +349,12 @@ subroutine end_snapshot(itime, num) write(*,'(" Time for writing snapshots (s): ",F12.8)') tend-tstart endif +#ifdef DEBUG + if (nrank == 0) then + print *, "# end_snapshot end" + endif +#endif + end subroutine end_snapshot ! @@ -374,7 +385,7 @@ subroutine write_xdmf_header(pathname, filename, num) if (nrank.eq.0) then OPEN(newunit=ioxdmf,file="./data/"//pathname//"/"//filename//'-'//num//'.xdmf') - write(ioxdmf,'(A22)')'' + write(ioxdmf,'(A)')'' write(ioxdmf,*)'' write(ioxdmf,*)'' write(ioxdmf,*)'' @@ -475,14 +486,26 @@ subroutine write_xdmf_footer() implicit none +#ifdef DEBUG + if (nrank == 0) then + print *, "## write_xdmf_footer start" + endif +#endif + if (nrank.eq.0) then write(ioxdmf,'(/)') - write(ioxdmf,*)' ' - write(ioxdmf,*)'' - write(ioxdmf,'(A7)')'' + write(ioxdmf,*) ' ' + write(ioxdmf,*) '' + write(ioxdmf,'(A)') '' close(ioxdmf) endif +#ifdef DEBUG + if (nrank == 0) then + print *, "## write_xdmf_footer end" + endif +#endif + end subroutine write_xdmf_footer ! From bd3dd2a5f8be0fa59454f35256246cca5c69323f Mon Sep 17 00:00:00 2001 From: cflag Date: Wed, 9 Mar 2022 09:50:49 +0100 Subject: [PATCH 5/5] Allow case-specific visualization --- Makefile | 2 +- decomp2d/CMakeLists.txt | 1 + decomp2d/Makefile.am | 2 + src/CMakeLists.txt | 1 + src/Makefile.am | 8 ++-- src/bc_tgv.f90 | 94 ++++++++++++++++++++++++++++++++++++++++- src/case.f90 | 36 ++++++++++++---- src/x3d_tools.f90 | 1 + src/xcompact3d.f90 | 8 ++-- 9 files changed, 135 insertions(+), 18 deletions(-) diff --git a/Makefile b/Makefile index aa914e7..b899e11 100644 --- a/Makefile +++ b/Makefile @@ -64,7 +64,7 @@ SRCDIR = ./src SRCDECOMP = $(DECOMPDIR)/decomp_2d.f90 $(DECOMPDIR)/glassman.f90 $(DECOMPDIR)/fft_$(FFT).f90 $(DECOMPDIR)/io.f90 OBJDECOMP = $(SRCDECOMP:%.f90=%.o) OBJ = $(SRC:%.f90=%.o) -SRC = $(SRCDIR)/x3d_precision.f90 $(SRCDIR)/module_param.f90 $(SRCDIR)/time_integrators.f90 $(SRCDIR)/tools.f90 $(SRCDIR)/x3d_transpose.f90 $(SRCDIR)/var.f90 $(SRCDIR)/thomas.f90 $(SRCDIR)/x3d_operator_x_data.f90 $(SRCDIR)/x3d_operator_y_data.f90 $(SRCDIR)/x3d_operator_z_data.f90 $(SRCDIR)/x3d_operator_1d.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/x3d_derive.f90 $(SRCDIR)/x3d_staggered.f90 $(SRCDIR)/x3d_filters.f90 $(SRCDIR)/bc_tgv.f90 $(SRCDIR)/navier.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/bc_tgv2d.f90 $(SRCDIR)/case.f90 $(SRCDIR)/transeq.f90 $(SRCDIR)/visu.f90 $(SRCDIR)/x3d_tools.f90 $(SRCDIR)/xcompact3d.f90 +SRC = $(SRCDIR)/x3d_precision.f90 $(SRCDIR)/module_param.f90 $(SRCDIR)/time_integrators.f90 $(SRCDIR)/tools.f90 $(SRCDIR)/x3d_transpose.f90 $(SRCDIR)/var.f90 $(SRCDIR)/thomas.f90 $(SRCDIR)/x3d_operator_x_data.f90 $(SRCDIR)/x3d_operator_y_data.f90 $(SRCDIR)/x3d_operator_z_data.f90 $(SRCDIR)/x3d_operator_1d.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/x3d_derive.f90 $(SRCDIR)/x3d_staggered.f90 $(SRCDIR)/x3d_filters.f90 $(SRCDIR)/visu.f90 $(SRCDIR)/bc_tgv.f90 $(SRCDIR)/navier.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/bc_tgv2d.f90 $(SRCDIR)/case.f90 $(SRCDIR)/transeq.f90 $(SRCDIR)/x3d_tools.f90 $(SRCDIR)/xcompact3d.f90 #######FFT settings########## ifeq ($(FFT),fftw3) diff --git a/decomp2d/CMakeLists.txt b/decomp2d/CMakeLists.txt index 8ba49d1..44bb5c8 100644 --- a/decomp2d/CMakeLists.txt +++ b/decomp2d/CMakeLists.txt @@ -1,4 +1,5 @@ file(GLOB files_decomp decomp_2d.f90 + io.f90 glassman.f90) include_directories(${CMAKE_SOURCE_DIR}/decomp2d) diff --git a/decomp2d/Makefile.am b/decomp2d/Makefile.am index 2530f42..7983b7b 100644 --- a/decomp2d/Makefile.am +++ b/decomp2d/Makefile.am @@ -9,6 +9,7 @@ endif lib_LTLIBRARIES = libdecomp2d.la libdecomp2d_la_SOURCES = \ decomp_2d.f90 \ +io.f90 \ glassman.f90 \ fft_@FFTENGINE@.f90 libdecomp2d_la_LDFLAGS = -no-undefined @@ -18,6 +19,7 @@ endif # Allow parallel make glassman.lo: decomp_2d.lo +io.lo: decomp_2d.lo if HAVE_MKL fft_@FFTENGINE@.lo: glassman.lo mkl_dfti.lo else diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 38a8ecf..71b5df7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -17,6 +17,7 @@ set(X3DFILES tools.f90 transeq.f90 var.f90 + visu.f90 x3d_derive.f90 x3d_filters.f90 x3d_operator_1d.f90 diff --git a/src/Makefile.am b/src/Makefile.am index b6440c9..5228064 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -8,6 +8,7 @@ time_integrators.f90 \ tools.f90 \ x3d_transpose.f90 \ var.f90 \ +visu.f90 \ thomas.f90 \ x3d_operator_x_data.f90 \ x3d_operator_y_data.f90 \ @@ -31,8 +32,9 @@ time_integrators.lo: module_param.lo tools.lo: module_param.lo x3d_transpose.lo: module_param.lo var.lo: module_param.lo +visu.lo: module_param.lo thomas.lo: module_param.lo -x3d_operator_1d.f90: module_param.lo x3d_operator_x_data.lo x3d_operator_y_data.lo x3d_operator_z_data.lo +x3d_operator_1d.lo: module_param.lo x3d_operator_x_data.lo x3d_operator_y_data.lo x3d_operator_z_data.lo poisson.lo: x3d_operator_x_data.lo x3d_operator_y_data.lo x3d_operator_z_data.lo x3d_derive.lo: x3d_operator_1d.lo x3d_staggered.lo: x3d_operator_1d.lo @@ -40,10 +42,10 @@ x3d_filters.lo: x3d_operator_1d.lo navier.lo: poisson.lo x3d_staggered.lo parameters.lo: module_param.lo bc_tgv2d.lo: x3d_operator_x_data.lo x3d_operator_y_data.lo -bc_tgv.lo: x3d_derive.lo +bc_tgv.lo: x3d_derive.lo visu.lo case.lo: bc_tgv2d.lo bc_tgv.lo transeq.lo: x3d_derive.lo case.lo -x3d_tools.lo: x3d_derive.lo case.lo +x3d_tools.lo: x3d_derive.lo case.lo visu.lo clean-local: -rm -f *.mod diff --git a/src/bc_tgv.f90 b/src/bc_tgv.f90 index ec5b9d4..55f132a 100644 --- a/src/bc_tgv.f90 +++ b/src/bc_tgv.f90 @@ -25,6 +25,8 @@ module bc_tgv tgv_listing, & tgv_init, & tgv_postprocess, & + tgv_visu_init, & + tgv_visu, & tgv_finalize contains @@ -129,7 +131,7 @@ subroutine tgv_init(ux1, uy1, uz1) ! subroutine tgv_postprocess(ux1, uy1, uz1, ndt) - use decomp_2d, only : mytype, real_type, decomp_2d_abort, xsize, ysize, zsize + use decomp_2d, only : decomp_2d_abort, ysize, zsize use param, only : half, two, xnu, dt use variables, only : nx, ny, nz, ppy use var, only : ux2,uy2,uz2,ux3,uy3,uz3 @@ -254,6 +256,96 @@ subroutine tgv_postprocess(ux1, uy1, uz1, ndt) end subroutine tgv_postprocess + ! + ! Register variables for visualization + ! + subroutine tgv_visu_init() + +#ifdef ADIOS2 + use decomp_2d_io, only : adios2_register_variable + use visu, only : io_write_real_coarse +#endif + + implicit none + +#ifdef ADIOS2 + call adios2_register_variable(io_write_real_coarse, "critq", 1, 2) +#endif + + end subroutine tgv_visu_init + + ! + ! Output case-specific variables + ! + subroutine tgv_visu(ux1, uy1, uz1, num) + + use decomp_2d, only : ysize, zsize + use variables, only : ppy + use param, only : half + use var, only : ux2,uy2,uz2,ux3,uy3,uz3 + use var, only : ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1 + use var, only : ta2,tb2,tc2 + use var, only : ta3,tb3,tc3 + use x3d_derive + use x3d_transpose + use x3d_operator_1d + use visu, only : write_field + + implicit none + + ! Arguments + real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: ux1, uy1, uz1 + character(len=32), intent(in) :: num + + ! Local variable + integer :: i, j, k + real(mytype), dimension(xsize(1),xsize(2),xsize(3)) :: array + + ! This is needed to compute derivatives + call x3d_transpose_x_to_y(ux1, ux2) + call x3d_transpose_x_to_y(uy1, uy2) + call x3d_transpose_x_to_y(uz1, uz2) + call x3d_transpose_y_to_z(ux2, ux3) + call x3d_transpose_y_to_z(uy2, uy3) + call x3d_transpose_y_to_z(uz2, uz3) + + ! Compute X derivative + call derx (ta1,ux1,x3d_op_derx, xsize(1),xsize(2),xsize(3)) + call derx (tb1,uy1,x3d_op_derxp,xsize(1),xsize(2),xsize(3)) + call derx (tc1,uz1,x3d_op_derxp,xsize(1),xsize(2),xsize(3)) + + ! Compute Y derivative and transpose back to X + call dery (ta2,ux2,x3d_op_deryp,ppy,ysize(1),ysize(2),ysize(3)) + call dery (tb2,uy2,x3d_op_dery ,ppy,ysize(1),ysize(2),ysize(3)) + call dery (tc2,uz2,x3d_op_deryp,ppy,ysize(1),ysize(2),ysize(3)) + call x3d_transpose_y_to_x(ta2, td1) + call x3d_transpose_y_to_x(tb2, te1) + call x3d_transpose_y_to_x(tc2, tf1) + + ! Compute Z derivative and transpose back to X + call derz (ta3,ux3,x3d_op_derzp,zsize(1),zsize(2),zsize(3)) + call derz (tb3,uy3,x3d_op_derzp,zsize(1),zsize(2),zsize(3)) + call derz (tc3,uz3,x3d_op_derz ,zsize(1),zsize(2),zsize(3)) + call x3d_transpose_z_to_y(ta3, ta2) + call x3d_transpose_z_to_y(tb3, tb2) + call x3d_transpose_z_to_y(tc3, tc2) + call x3d_transpose_y_to_x(ta2, tg1) + call x3d_transpose_y_to_x(tb2, th1) + call x3d_transpose_y_to_x(tc2, ti1) + !du/dx=ta1 du/dy=td1 du/dz=tg1 + !dv/dx=tb1 dv/dy=te1 dv/dz=th1 + !dw/dx=tc1 dw/dy=tf1 dw/dz=ti1 + + ! Compute Q criterion and write it + do concurrent (k=1:xsize(3), j=1:xsize(2), i=1:xsize(1)) + array(i,j,k) = - half*(ta1(i,j,k)**2+te1(i,j,k)**2+ti1(i,j,k)**2) & + - td1(i,j,k)*tb1(i,j,k)& + - tg1(i,j,k)*tc1(i,j,k)& + - th1(i,j,k)*tf1(i,j,k) + enddo + call write_field(array, ".", "critq", trim(num)) + + end subroutine tgv_visu ! ! Finalize case-specific IO ! diff --git a/src/case.f90 b/src/case.f90 index d2d1049..c1f3c11 100644 --- a/src/case.f90 +++ b/src/case.f90 @@ -18,6 +18,7 @@ module case case_init, & case_bc, & case_forcing, & + case_visu_init, & case_visu, & case_postprocess, & case_finalize @@ -94,36 +95,55 @@ subroutine case_forcing(dux1, duy1, duz1) end subroutine case_forcing + ! + ! Register variables for visualization + ! + subroutine case_visu_init() + + implicit none + + if (itype == itype_tgv) call tgv_visu_init() + + end subroutine case_visu_init + ! ! Visualization ! This is called when itime % ioutput = 0 ! - subroutine case_visu() + subroutine case_visu(ux1, uy1, uz1, num) implicit none + ! Arguments + real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: ux1, uy1, uz1 + character(len=32), intent(in) :: num + + if (itype == itype_tgv) call tgv_visu(ux1, uy1, uz1, num) + end subroutine case_visu ! ! Case-specific post-processing ! This is called at the end of each time step ! - subroutine case_postprocess(ux1, uy1, uz1, ndt) + subroutine case_postprocess(ux1, uy1, uz1, it, num) - use param, only : ivisu, ioutput + use param, only : ioutput + use variables, only : nvisu implicit none ! Arguments real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: ux1, uy1, uz1 - integer, intent(in) :: ndt + integer, intent(in) :: it + character(len=32), intent(in) :: num - if (itype == itype_tgv) call tgv_postprocess(ux1, uy1, uz1, ndt) + if (itype == itype_tgv) call tgv_postprocess(ux1, uy1, uz1, it) - if (itype == itype_tgv2d) call tgv2d_postprocess(ux1, uy1, uz1, ndt) + if (itype == itype_tgv2d) call tgv2d_postprocess(ux1, uy1, uz1, it) - if ((ivisu /= 0).and.(ioutput /= 0)) then - if (mod(ndt, ioutput) == 0) call case_visu() + if ((nvisu /= 0).and.(ioutput /= 0)) then + if (mod(it, ioutput) == 0) call case_visu(ux1, uy1, uz1, num) endif end subroutine case_postprocess diff --git a/src/x3d_tools.f90 b/src/x3d_tools.f90 index 81c2a96..cf1cc42 100644 --- a/src/x3d_tools.f90 +++ b/src/x3d_tools.f90 @@ -183,6 +183,7 @@ subroutine init_xcompact3d(ndt_max) call decomp_info_init(nxm,nym,nzm,phG) call visu_init() + call case_visu_init() endsubroutine init_xcompact3d !######################################################################## diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90 index 1d3502e..2858700 100644 --- a/src/xcompact3d.f90 +++ b/src/xcompact3d.f90 @@ -74,11 +74,9 @@ program xcompact3d !! End of timestep itime = itime + 1 ndt = ndt + 1 - if (mod(itime,ioutput) == 0) then - call write_snapshot(ux1, uy1, uz1, itime, num) - call end_snapshot(itime, num) - endif - call case_postprocess(ux1, uy1, uz1, ndt) + if (mod(itime,ioutput) == 0) call write_snapshot(ux1, uy1, uz1, itime, num) + call case_postprocess(ux1, uy1, uz1, itime, num) + if (mod(itime,ioutput) == 0) call end_snapshot(itime, num) end do