module snapshot_m !! Writing and reading snapshot files use MPI use comm_handler_m, only : comm_handler_t use netcdf use precision_grillix_m, only : GP use screen_io_m, only : get_stdout use parallelisation_setup_m, only : is_rank_info_writer use error_handling_grillix_m, only: handle_error_netcdf, handle_error use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use variable_m, only : variable_t use mesh_cart_m, only : mesh_cart_t implicit none type, public :: char_arr_t !! Can be used for character arrays of variable length character(len=:), allocatable, public :: str !! String end type type, public :: snapshot_t !! Datatype for snapshout I/O character(len=:), allocatable, private :: dirpath !! Directory path for snapshot files integer, private :: nvars !! Number of variables stored in snapshotfile type(char_arr_t), allocatable, dimension(:), private :: name_vars !! Names of variables integer, private :: npoints_cano !! Number of canonical mesh points integer, private :: npoints_stag !! Number of staggered mesh points integer, private :: npoints_max !! Maximum number over all mesh points, over phi !! including maximum over canonical and staggered contains procedure, public :: init => init_snapshot procedure, public :: write_snaps procedure, public :: read_snaps procedure, public :: read_snaps_last procedure, public :: display => display_snapshot final :: destructor end type contains subroutine init_snapshot(self, comm_handler, mesh_cano, mesh_stag, dirpath, nvars, name_vars) !! Initialises snapshots, creates dirpath class(snapshot_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! MPI communicator type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) character(len=*), intent(in) :: dirpath !! Directory, where snapshot files are written integer, intent(in) :: nvars !! Number of variables to be written type(char_arr_t), dimension(nvars), intent(in) :: name_vars !! Names of variables integer :: cmd_exit, cmd_err, ierr integer :: npoints_cano_max, npoints_stag_max self%dirpath = dirpath self%nvars = nvars self%npoints_cano = mesh_cano%get_n_points() self%npoints_stag = mesh_stag%get_n_points() call MPI_allreduce(self%npoints_cano, npoints_cano_max, & 1, MPI_INTEGER, MPI_MAX, & comm_handler%get_comm_planes(), ierr) call MPI_allreduce(self%npoints_stag, npoints_stag_max, & 1, MPI_INTEGER, MPI_MAX, & comm_handler%get_comm_planes(), ierr) self%npoints_max = max(npoints_cano_max, npoints_stag_max) allocate(self%name_vars(nvars), source = name_vars) ! Create snapsfile folder if (is_rank_info_writer) then call EXECUTE_COMMAND_LINE('mkdir '//dirpath, & exitstat = cmd_exit, cmdstat = cmd_err) if (cmd_err /= 0) then write(get_stdout(),*)'error(init_snapshot): cmd_err = ', cmd_err endif endif ! Barrier ensuring that folder is created ! before another process possibly writes call MPI_BARRIER(comm_handler%get_comm_cart(), ierr) call self%display() end subroutine subroutine write_snaps(self, comm_handler, isnaps, tau, vars, isubsnaps, allow_overwrite) !! Writes snapshot class(snapshot_t), intent(in) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! MPI communicator integer, intent(in) :: isnaps !! Snapshot number real(GP), intent(in) :: tau !! Time type(variable_t), dimension(self%nvars), intent(in) :: vars !! Variables integer, optional, intent(in) :: isubsnaps !! Optional sub-snapshot number logical, optional, intent(in) :: allow_overwrite !! Whether to allow overwriting snapfiles integer :: ivar character(len=5) :: isnaps_c character(len=3) :: isubsnaps_c character(len=:), allocatable :: snapsfile logical :: allow_overwrite_arg integer :: nf90_stat, nf90_id integer :: id_snaps_dim, id_nplanes_dim, id_species_dim, id_npoints_max_dim integer :: id_tau, id_npoints_cano, id_npoints_stag write(isnaps_c,'(I5.5)')isnaps if (.not. present(isubsnaps)) then snapsfile = self%dirpath // '/snaps_t' // isnaps_c // '.nc' else write(isubsnaps_c,'(I3.3)')isubsnaps snapsfile = self%dirpath // '/fields_d' // isubsnaps_c // '_t' // isnaps_c // '.nc' endif if (present(allow_overwrite)) then allow_overwrite_arg = allow_overwrite else allow_overwrite_arg = .false. endif if (allow_overwrite_arg) then nf90_stat = nf90_create(snapsfile, & NF90_NETCDF4 + NF90_CLOBBER, nf90_id, & comm=comm_handler%get_comm_cart(), info=MPI_INFO_NULL) else nf90_stat = nf90_create(snapsfile, & NF90_NETCDF4 + NF90_NOCLOBBER, nf90_id, & comm=comm_handler%get_comm_cart(), info=MPI_INFO_NULL) endif call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Define dimensions nf90_stat = nf90_def_dim(nf90_id, 'snaps', 1, id_snaps_dim) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_dim(nf90_id, 'planes', & comm_handler%get_nplanes(), id_nplanes_dim) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_dim(nf90_id, 'species', & comm_handler%get_nspecies(), id_species_dim) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_dim(nf90_id, 'npoints_max', & self%npoints_max, id_npoints_max_dim) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Define tau and npoints nf90_stat = nf90_def_var(nf90_id, 'tau', NF90_DOUBLE, & [id_snaps_dim], id_tau) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_var(nf90_id, 'npoints_cano', NF90_DOUBLE, & [id_nplanes_dim], id_npoints_cano) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_var(nf90_id, 'npoints_stag', NF90_DOUBLE, & [id_nplanes_dim], id_npoints_stag) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_enddef(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) if (is_rank_info_writer) then nf90_stat = nf90_put_var(nf90_id, id_tau, [tau]) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) endif nf90_stat = nf90_put_var(nf90_id, id_tau, [tau]) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_put_var(nf90_id, id_npoints_cano, [self%npoints_cano], & start=[comm_handler%get_plane()+1], count=[1]) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_put_var(nf90_id, id_npoints_stag, [self%npoints_stag], & start=[comm_handler%get_plane()+1], count=[1]) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_redef(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Write variables do ivar = 1, self%nvars call vars(ivar)%write_netcdf(comm_handler, nf90_id, & self%name_vars(ivar)%str) enddo nf90_stat = nf90_close(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) end subroutine subroutine read_snaps_last(self, comm_handler, isnaps, tau, vars) !! Reads last available snapshot class(snapshot_t), intent(in) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! MPI communicator integer, intent(out) :: isnaps !! Snapshot number real(GP), intent(out) :: tau !! Time type(variable_t), dimension(self%nvars), intent(inout) :: vars !! Variables logical :: fexist integer :: it character(len=5) :: isnaps_c character(len=:), allocatable :: snapsfile do it = 99998, 1, -1 write(isnaps_c,'(I5.5)')it snapsfile = self%dirpath // '/snaps_t' // isnaps_c // '.nc' inquire(file=snapsfile, exist=fexist) if (fexist) then exit endif enddo isnaps = it ! Check if any snapshotfile was found if (isnaps == 0) then call handle_error('No snapsfile for continuation found', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) endif call self%read_snaps(comm_handler, isnaps, tau, vars) end subroutine subroutine read_snaps(self, comm_handler, isnaps, tau, vars) !! Reads snapshot class(snapshot_t), intent(in) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! MPI communicator integer, intent(in) :: isnaps !! Snapshot number real(GP), intent(out) :: tau !! Time type(variable_t), dimension(self%nvars), intent(inout) :: vars !! Variables integer :: ivar character(len=5) :: isnaps_c character(len=:), allocatable :: snapsfile integer :: nf90_stat, nf90_id, id_tau real(GP), dimension(1) :: wrk write(isnaps_c,'(I5.5)')isnaps snapsfile = self%dirpath // '/snaps_t' // isnaps_c // '.nc' nf90_stat = nf90_open(snapsfile, NF90_NETCDF4, nf90_id ) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! Read tau nf90_stat = nf90_inq_varid(nf90_id, 'tau', id_tau) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_get_var(nf90_id, id_tau, wrk) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) tau = wrk(1) do ivar = 1, self%nvars call vars(ivar)%read_netcdf(comm_handler, nf90_id, & self%name_vars(ivar)%str) enddo nf90_stat = nf90_close(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) end subroutine subroutine display_snapshot(self) !! Displays information class(snapshot_t), intent(inout) :: self !! Instance of the type integer :: ivar if (.not.is_rank_info_writer) return write(get_stdout(),*)'Information on snapshot -------------------------------' write(get_stdout(),101) self%dirpath, self%nvars 101 FORMAT( 3X,'dirpath = ',A32 /, & 3X,'nvars = ',I8) do ivar = 1, self%nvars write(get_stdout(),102)self%name_vars(ivar)%str 102 FORMAT( 3X,'name_var = ',A32) enddo write(get_stdout(),*)'-------------------------------------------------------' end subroutine subroutine destructor(self) !! Destructor type(snapshot_t), intent(inout) :: self !! Instance of the type if (allocated(self%dirpath)) deallocate(self%dirpath) self%nvars = 0 if (allocated(self%name_vars)) deallocate(self%name_vars) end subroutine end module