#include "unused_dummy.inc" !> FIXME : Add documentation module gs2_save use constants, only: run_name_size implicit none private public :: gs2_restore, gs2_save_for_restart, finish_gs2_save public :: read_many, save_many, gs2_save_response, gs2_restore_response public :: init_gs2_save, init_dt, read_t0_from_restart_file, init_ant_amp public :: set_restart_file, include_explicit_source_in_restart, proc_to_save_fields public :: init_vnm, restart_writable, EigNetcdfID public :: init_eigenfunc_file, finish_eigenfunc_file, add_eigenpair_to_file public :: gs2_save_slice interface gs2_restore module procedure gs2_restore_many end interface !> A custom type to look after the netcdf ids for the eigenvalue file type EigNetcdfID integer :: ncid !< File handle integer :: nconv_count !< Current size of `conv` dimension end type EigNetcdfID !> Read and write single or multiple restart files. !> Effectively ignored if not built with parallel netcdf logical :: read_many = .false., save_many = .false. !> Do we want to include the explicit source terms and !> related timesteps when saving/restoring from the restart file. logical :: include_explicit_source_in_restart = .true. !> Which processor should save the potentials. If <0 (default) then !> all processors save their potentials. If >= nproc then will be !> set to mod(proc_to_save_fields, nproc) integer :: proc_to_save_fields !> Used to hold the base of the restart file name. character (run_name_size) :: restart_file #ifdef NETCDF !> The NETCDF_PARALLEL directives include code for parallel !> netcdf using HDF5 to write the output to a single restart file !> The save_many/read_many flag allows the old style multiple file output !> Here we set a run time flag to enable us to handle some of the different !> paths through the code at run time rather than with ifdef mazes #ifdef NETCDF_PARALLEL logical, parameter :: has_netcdf_parallel = .true. #else logical, parameter :: has_netcdf_parallel = .false. #endif #endif contains #ifdef NETCDF integer function get_dim_size(ncid, dim_name) result(dim_size) use neasyf, only: neasyf_dim use netcdf, only: nf90_inquire_dimension implicit none integer, intent(in) :: ncid character(len = *), intent(in) :: dim_name integer :: id, istatus call neasyf_dim(ncid, dim_name, dimid = id) istatus = nf90_inquire_dimension(ncid, id, len=dim_size) end function get_dim_size subroutine verify_compatible_dim_sizes(ncid, dim_name, expected_size) use mp, only: mp_abort, iproc integer, intent(in) :: ncid, expected_size character(len = *), intent(in) :: dim_name integer :: actual_size actual_size = get_dim_size(ncid, dim_name) if (actual_size /= expected_size) then write(*,'("Restart error (",A,"): got/expected = ",I0,"/",I0," : ",I0)') & dim_name, actual_size, expected_size, iproc call mp_abort('Incompatible '//dim_name//' in restart', .true.) end if end subroutine verify_compatible_dim_sizes !> Returns the file corresponding to restart file in current setup function get_file_proc(is_one_file_per_processor, fileopt) result(file_proc) use mp, only: iproc use optionals, only: get_option_with_default implicit none logical, intent(in) :: is_one_file_per_processor character(len=*), intent(in), optional :: fileopt character(run_name_size) :: file_proc file_proc = trim(restart_file) // trim(adjustl(get_option_with_default(fileopt, ''))) if (is_one_file_per_processor) write(file_proc, '(A,".",I0)') trim(file_proc), iproc end function get_file_proc logical function variable_exists(ncid, var_name) use netcdf, only: nf90_inq_varid, NF90_NOERR implicit none integer, intent(in) :: ncid character(len = *), intent(in) :: var_name integer :: id variable_exists = nf90_inq_varid(ncid, var_name, id) == NF90_NOERR end function variable_exists #endif !> Create and fill the restart files subroutine gs2_save_for_restart & (g, t0, vnm, has_phi, has_apar, has_bpar, & code_dt, code_dt_prev1, code_dt_prev2, code_dt_max, & fileopt, & save_glo_info_and_grids, save_velocities, header) #ifdef NETCDF use fields_arrays, only: phinew, aparnew, bparnew use dist_fn_arrays, only: gexp_1, gexp_2, gexp_3, kx_shift, theta0_shift use kt_grids, only: naky, ntheta0, jtwist_out, akx, aky use antenna_data, only: a_ant, b_ant, ant_on use gs2_metadata, only: create_metadata use netcdf, only: nf90_collective, nf90_independent use neasyf, only: neasyf_open, neasyf_dim, neasyf_write, neasyf_close use mp, only: iproc, barrier, proc0, mp_comm, mp_info use theta_grid, only: theta use gs2_layouts, only: proc_id,ik_idx,it_idx,is_idx,ie_idx,il_idx, layout use layouts_type, only: g_layout_type use le_grids, only: energy, al, negrid, nlambda use species, only: nspec use dist_fn_arrays, only: vpa, vperp2 use unit_tests, only: debug_message use optionals, only: get_option_with_default #else use file_utils, only: error_unit use mp, only: proc0 #endif use theta_grid, only: ntgrid use gs2_layouts, only: g_lo use standard_header, only : standard_header_type implicit none !> The distribution function \(g\) complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in) :: g !> Current simulation time real, intent (in) :: t0 !> FIXME: Collisionality multiplier? real, dimension (2), intent (in) :: vnm !> The current time step, the two previous time steps and the !> maximum allowed time step. real, intent(in) :: code_dt, code_dt_prev1, code_dt_prev2, code_dt_max !> Do we include electrostatic potential, \(\phi\) logical, intent (in) :: has_phi !> Do we include vector potential, \(A_\parallel\) logical, intent (in) :: has_apar !> Do we include parallel magnetic field, \(B_\parallel\) logical, intent (in) :: has_bpar !> Optional extra filename infix character (len=*), intent (in), optional :: fileopt !> Control if layout and diagnostic information is saved logical, intent (in), optional :: save_glo_info_and_grids !> Control if \(v_\parallel, v_\perp\) information is saved logical, intent (in), optional :: save_velocities !> Header for files with build and run information type(standard_header_type), optional, intent(in) :: header #ifdef NETCDF !> Names of dimensions for distribution-sized variables in netCDF files character(len=5), dimension(*), parameter :: g_dim_names = ["theta", "sign ", "glo "] character(len=5), dimension(*), parameter :: field_dim_names = ["theta", "akx ", "aky "] !> Arrays holding information about the layout for plotting integer, dimension (:,:,:,:,:), allocatable :: procid_kykxsel, iglo_kykxsel character (run_name_size) :: file_proc real, dimension(:), allocatable :: atmp real, dimension(naky) :: stmp integer, dimension(5) :: glo_counts integer, dimension(3) :: start_pos, counts, field_counts integer, dimension(1) :: scalar_count, exb_count logical :: is_one_file_per_processor, this_proc_saves_fields, has_gexp_1 logical :: save_glo_info_and_grids_local, save_velocities_info logical :: should_define_single_entry_variables integer, parameter :: verb = 3 integer :: ncid, iglo, n_elements ! Whether to use collective (parallel) or independent access to netCDF variables integer :: par_access type(standard_header_type) :: local_header if (present(header)) then local_header = header else local_header = standard_header_type() end if call debug_message(verb, 'gs2_save::gs2_save_for_restart initialising') ! Handle some optional inputs save_glo_info_and_grids_local = get_option_with_default(save_glo_info_and_grids, .false.) save_velocities_info = get_option_with_default(save_velocities, .false.) ! Decide if we want one file per processor or not is_one_file_per_processor = save_many .or. (.not. has_netcdf_parallel) ! Open the file and define dimensions and variables file_proc = get_file_proc(is_one_file_per_processor, fileopt) call debug_message(verb, 'gs2_save::gs2_save_for_restart opening file') call debug_message(verb, 'file proc is') call debug_message(verb, trim(adjustl(file_proc))) if (is_one_file_per_processor) then ncid = neasyf_open(file_proc, "w") else call barrier call debug_message(verb, & 'gs2_save::gs2_save_for_restart called barrier before opening') ncid = neasyf_open(file_proc, "w", comm=mp_comm, info=mp_info) end if call create_metadata(ncid, "GS2 restart file", local_header) call debug_message(verb, 'gs2_save::gs2_save_for_restart defining dimensions') ! Work out the problem size if (is_one_file_per_processor) then n_elements = max(g_lo%ulim_proc-g_lo%llim_proc+1, 0) else n_elements = g_lo%ulim_world+1 end if ! Note when processors have no local data we find n_elements = 0 so here we ! would be trying to create a zero size dimension. Netcdf interprets this as ! an unlimited dimension so neasyf requires unlimited = .true. in this case, ! hence the unlimited argument below. This will lead to empty gr/gi data being ! written on such processors. call neasyf_dim(ncid, "glo", dim_size=n_elements, unlimited = n_elements == 0) call neasyf_dim(ncid, "theta", theta) call neasyf_dim(ncid, "sign", dim_size=2) call neasyf_dim(ncid, "aky", aky) call neasyf_dim(ncid, "akx", akx) ! Decide if this processor should take part in defining items ! which should only appear once in outputs. should_define_single_entry_variables = proc0 .or. (.not. is_one_file_per_processor) !For saving distribution function if (save_glo_info_and_grids_local .and. should_define_single_entry_variables) then call neasyf_dim(ncid, "lambda", values=al, long_name="Energy/magnetic moment", units="1/(2 B_a)") call neasyf_dim(ncid, "species", dim_size=nspec) ! Energy coordinate is a function of species, but the dimension must be 1D call neasyf_dim(ncid, "egrid", dim_size=negrid, long_name="See 'energy' for energy grid values") call neasyf_write(ncid, "energy", energy, dim_names=["egrid ", "species"], & long_name="Energy grid values") end if call debug_message(verb, 'gs2_save::gs2_save_for_restart defining variables') if (ant_on) call neasyf_dim(ncid, "nk_stir", dim_size=size(a_ant)) call debug_message(verb, 'gs2_save::gs2_save_for_restart writing scalars') if (is_one_file_per_processor .or. proc0) then scalar_count = [1] else scalar_count = [0] end if call neasyf_write(ncid, "delt0", code_dt, count=scalar_count) call neasyf_write(ncid, "delt1", code_dt_prev1, count=scalar_count) call neasyf_write(ncid, "delt2", code_dt_prev2, count=scalar_count) call neasyf_write(ncid, "delt_max", code_dt_max, count=scalar_count) call neasyf_write(ncid, "t0", t0, count=scalar_count) call neasyf_write(ncid, "layout", layout, count=scalar_count) call neasyf_write(ncid, "ntheta0", ntheta0, count=scalar_count) call neasyf_write(ncid, "naky", naky, count=scalar_count) call neasyf_write(ncid, "nlambda", nlambda, count=scalar_count) call neasyf_write(ncid, "negrid", negrid, count=scalar_count) call neasyf_write(ncid, "nspec", nspec, count=scalar_count) call neasyf_write(ncid, "vnm1", vnm(1), count=scalar_count) call neasyf_write(ncid, "vnm2", vnm(2), count=scalar_count) if (ant_on) then atmp = real(a_ant) call neasyf_write(ncid, "a_ant_r", atmp, dim_names=["nk_stir"]) atmp = aimag(a_ant) call neasyf_write(ncid, "a_ant_i", atmp, dim_names=["nk_stir"]) atmp = real(b_ant) call neasyf_write(ncid, "b_ant_r", atmp, dim_names=["nk_stir"]) atmp = aimag(b_ant) call neasyf_write(ncid, "b_ant_i", atmp, dim_names=["nk_stir"]) end if call debug_message(verb, 'gs2_save::gs2_save_for_restart writing dist fn') if (is_one_file_per_processor) then start_pos = [1, 1, 1] par_access = nf90_independent else start_pos = [1, 1, g_lo%llim_proc+1] par_access = nf90_collective end if counts = shape(g) call neasyf_write(ncid, "gr", real(g), dim_names=g_dim_names, & start=start_pos, count=counts, par_access=par_access) call neasyf_write(ncid, "gi", aimag(g), dim_names=g_dim_names, & start=start_pos, count=counts, par_access=par_access) if (include_explicit_source_in_restart) then !Explicit source term gexp_1 #ifndef SHMEM has_gexp_1 = allocated(gexp_1) #else has_gexp_1 = associated(gexp_1) #endif if (has_gexp_1) then call neasyf_write(ncid, "gexp_1_r", real(gexp_1), dim_names=g_dim_names, & start=start_pos, count=counts, par_access=par_access) call neasyf_write(ncid, "gexp_1_i", aimag(gexp_1), dim_names=g_dim_names, & start=start_pos, count=counts, par_access=par_access) end if !Explicit source term gexp_2 if(allocated(gexp_2)) then call neasyf_write(ncid, "gexp_2_r", real(gexp_2), dim_names=g_dim_names, & start=start_pos, count=counts, par_access=par_access) call neasyf_write(ncid, "gexp_2_i", aimag(gexp_2), dim_names=g_dim_names, & start=start_pos, count=counts, par_access=par_access) end if !Explicit source term gexp_3 if(allocated(gexp_3)) then call neasyf_write(ncid, "gexp_3_r", real(gexp_3), dim_names=g_dim_names, & start=start_pos, count=counts, par_access=par_access) call neasyf_write(ncid, "gexp_3_i", aimag(gexp_3), dim_names=g_dim_names, & start=start_pos, count=counts, par_access=par_access) end if end if ! For saving distribution function if (save_glo_info_and_grids_local) then if (proc0 .or. .not. is_one_file_per_processor) then if(.not. allocated(iglo_kykxsel)) allocate(iglo_kykxsel(naky,ntheta0,nspec,negrid,nlambda)) if(.not. allocated(procid_kykxsel)) allocate(procid_kykxsel(naky,ntheta0,nspec,negrid,nlambda)) end if if (proc0) then do iglo=g_lo%llim_world,g_lo%ulim_world iglo_kykxsel(ik_idx(g_lo,iglo),it_idx(g_lo,iglo),is_idx(g_lo,iglo), & ie_idx(g_lo,iglo),il_idx(g_lo,iglo)) = iglo end do do iglo=g_lo%llim_world,g_lo%ulim_world procid_kykxsel(ik_idx(g_lo,iglo),it_idx(g_lo,iglo),is_idx(g_lo,iglo), & ie_idx(g_lo,iglo),il_idx(g_lo,iglo)) = proc_id(g_lo,iglo) end do glo_counts = shape(iglo_kykxsel) else glo_counts = [0, 0, 0, 0, 0] end if ! Either only write to proc0's serial file, or all processes need to call ! these functions in parallel if (proc0 .or. .not. is_one_file_per_processor) then call neasyf_write(ncid, "iglo_llim_proc", g_lo%llim_proc, count=scalar_count) ! Fill processor and layout information for plotting purposes call neasyf_write(ncid, "iglo_kykxsel", iglo_kykxsel, & dim_names=[character(len=7)::"aky", "akx", "species", "egrid", "lambda"], & count=glo_counts) call neasyf_write(ncid, "procid_kykxsel", procid_kykxsel, & dim_names=[character(len=7)::"aky", "akx", "species", "egrid", "lambda"], & count=glo_counts) call neasyf_write(ncid, "jtwist", jtwist_out, count=scalar_count) end if end if if (save_velocities_info) then call neasyf_write(ncid, "vpa", vpa, dim_names=["theta", "sign ", "glo "], & start=start_pos, count=counts) call neasyf_write(ncid, "vperp2", vperp2, dim_names=["theta", "glo "], & start=start_pos(1:3:2), count=counts(1:3:2)) end if ! Decide if we are saving the fields this_proc_saves_fields = (iproc == proc_to_save_fields) .or. (proc_to_save_fields < 0) if (this_proc_saves_fields .or. .not. is_one_file_per_processor) then ! In parallel, all processes need to call `neasyf_write`, but ! only one should actually write any data if (this_proc_saves_fields) then field_counts = shape(phinew) else field_counts = [0, 0, 0] end if if (has_phi) then call neasyf_write(ncid, "phi_r", real(phinew), dim_names=field_dim_names, count=field_counts) call neasyf_write(ncid, "phi_i", aimag(phinew), dim_names=field_dim_names, count=field_counts) end if if (has_apar) then call neasyf_write(ncid, "apar_r", real(aparnew), dim_names=field_dim_names, count=field_counts) call neasyf_write(ncid, "apar_i", aimag(aparnew), dim_names=field_dim_names, count=field_counts) end if if (has_bpar) then call neasyf_write(ncid, "bpar_r", real(bparnew), dim_names=field_dim_names, count=field_counts) call neasyf_write(ncid, "bpar_i", aimag(bparnew), dim_names=field_dim_names, count=field_counts) end if end if ! Because we want to be able to restart using exb shear from a ! case which does not have exb shear we always add kx_shift and ! theta0_shift even if no exb shear present in simulation) if (is_one_file_per_processor .or. proc0) then exb_count = [naky] else exb_count = [0] end if if (allocated(kx_shift)) then stmp = kx_shift else stmp = 0. end if call neasyf_write(ncid, "kx_shift", stmp, dim_names=["aky"], count=exb_count) if (allocated(theta0_shift)) then stmp = theta0_shift else stmp = 0. end if call neasyf_write(ncid, "theta0_shift", stmp, dim_names=["aky"], count=exb_count) ! Close the file call neasyf_close(ncid) ! Deallocate local arrays if (allocated(iglo_kykxsel)) deallocate(iglo_kykxsel) if (allocated(procid_kykxsel)) deallocate(procid_kykxsel) #else if (proc0) write (error_unit(),*) 'WARNING: gs2_save_for_restart is called without netcdf library' UNUSED_DUMMY(g); UNUSED_DUMMY(t0); UNUSED_DUMMY(vnm); UNUSED_DUMMY(has_phi) UNUSED_DUMMY(has_apar); UNUSED_DUMMY(has_bpar); UNUSED_DUMMY(code_dt) UNUSED_DUMMY(code_dt_prev1); UNUSED_DUMMY(code_dt_prev2); UNUSED_DUMMY(code_dt_max) UNUSED_DUMMY(fileopt); UNUSED_DUMMY(save_glo_info_and_grids) UNUSED_DUMMY(save_velocities); UNUSED_DUMMY(header) #endif end subroutine gs2_save_for_restart !> Save g slices subroutine gs2_save_slice(g_in, it, ik, il, ie, is, time, nout, fileopt, header) # ifdef NETCDF use neasyf, only: neasyf_open, neasyf_dim, neasyf_write, neasyf_close use gs2_metadata, only: create_metadata use unit_tests, only: debug_message use mp, only: iproc, proc0 use optionals, only: get_option_with_default use convert, only: c2r use kt_grids, only: aky, akx, naky, ntheta0 use le_grids, only: energy, al, negrid, nlambda use species, only: nspec use theta_grid, only: theta use gs2_layouts, only: proc_id, idx # endif use theta_grid, only: ntgrid use gs2_layouts, only: g_lo use standard_header, only : standard_header_type implicit none !> The distribution function \(g\) complex, dimension (-ntgrid:, :, g_lo%llim_proc:), intent (in) :: g_in !> The indices of the slice to save integer, intent(in) :: it, ik, il, ie, is !> The current time real, intent(in) :: time !> The current time output index integer, intent(in), optional :: nout !> Optional extra filename infix character (len=*), intent (in), optional :: fileopt !> Header for files with build and run information type(standard_header_type), optional, intent(in) :: header # ifdef NETCDF character (run_name_size) :: file_proc integer :: the_iglo, the_proc integer, parameter :: verb = 3 type(standard_header_type) :: local_header character (len=*), parameter :: file_opt_base = '.slice' logical :: appending, previous_file_exists real, dimension(:, :, :), allocatable :: ri_gslice integer, save :: nout_local = 1 logical :: error_exit ! File handle integer :: ncid call debug_message(verb, 'gs2_save::gs2_save_slice initialising') ! Sanity check the inputs error_exit = .false. if (it < 1 .or. it > ntheta0) error_exit = .true. if (ik < 1 .or. ik > naky) error_exit = .true. if (il < 1 .or. il > nlambda) error_exit = .true. if (ie < 1 .or. ie > negrid) error_exit = .true. if (is < 1 .or. is > nspec) error_exit = .true. if (proc0 .and. error_exit) write(6,'("Error: Invalid index based to gs2_save_slice.")') if (error_exit) return the_iglo = idx(g_lo, ik, it, il, ie, is) the_proc = proc_id(g_lo, the_iglo) ! The proc which owns this slice is responsible for writing and ! does the rest of the work so let us return if that's not this ! proc if (iproc /= the_proc) return if (present(header)) then local_header = header else local_header = standard_header_type() end if ! Open the file and define dimensions and variables file_proc = get_file_proc(.false., file_opt_base // fileopt) call debug_message(verb, 'gs2_save::gs2_save_slice opening file') call debug_message(verb, 'file proc is') call debug_message(verb, trim(adjustl(file_proc))) nout_local = get_option_with_default(nout, nout_local) ! See if the output file already exists, so we know how to open ! the file. This would be simpler if netCDF had a "open if it ! exists, create if it doesn't" mode, but it doesn't so we have ! to do this instead inquire(file=trim(file_proc), exist=previous_file_exists) appending = (previous_file_exists .and. nout_local > 1) if (appending) then ncid = neasyf_open(trim(file_proc), 'rw') else ncid = neasyf_open(trim(file_proc), 'w') call create_metadata(ncid, "GS2 slice file", local_header) call debug_message(verb, 'gs2_save::gs2_save_slice defining dimensions') call neasyf_dim(ncid, "ri", dim_size = 2, long_name = "real/imaginary components") call neasyf_dim(ncid, "theta", values = theta, long_name="Poloidal angle", units="rad") call neasyf_dim(ncid, "sign", dim_size = 2, long_name = "sign of v||") call neasyf_dim(ncid, "t", unlimited = .true., long_name = "Time", units = "L/vt") call debug_message(verb, 'gs2_save::gs2_save_slice defining constant variables') call neasyf_write(ncid, "it", it) call neasyf_write(ncid, "ik", ik) call neasyf_write(ncid, "il", il) call neasyf_write(ncid, "ie", ie) call neasyf_write(ncid, "is", is) call neasyf_write(ncid, "iglo", the_iglo) call neasyf_write(ncid, "iproc", the_proc) call neasyf_write(ncid, "kx", akx(it)) call neasyf_write(ncid, "ky", aky(ik)) call neasyf_write(ncid, "lambda", al(il)) call neasyf_write(ncid, "energy", energy(ie, is)) end if call neasyf_write(ncid, "t", time, start = [nout_local]) ! Convert from complex to real allocate(ri_gslice(2, size(theta), 2)) call c2r(g_in(:, :, the_iglo), ri_gslice) call neasyf_write(ncid, "g_slice", ri_gslice, start = [1, 1, 1, nout_local]) ! Close the file call neasyf_close(ncid) nout_local = nout_local + 1 # else UNUSED_DUMMY(g_in); UNUSED_DUMMY(it); UNUSED_DUMMY(ik) ; UNUSED_DUMMY(il) UNUSED_DUMMY(ie); UNUSED_DUMMY(is); UNUSED_DUMMY(time); UNUSED_DUMMY(nout) UNUSED_DUMMY(fileopt); UNUSED_DUMMY(header) # endif end subroutine gs2_save_slice !> FIXME : Add documentation subroutine gs2_restore_many (g, scale, has_phi, has_apar, has_bpar, fileopt) !MR, 2007: restore kx_shift array if already allocated # ifdef NETCDF use array_utils, only: zero_array use mp, only: iproc, proc0, broadcast, mp_comm, mp_info use fields_arrays, only: phinew, aparnew, bparnew use fields_arrays, only: phi, apar, bpar use dist_fn_arrays, only: gexp_1, gexp_2, gexp_3, kx_shift, theta0_shift use kt_grids, only: naky, ntheta0 use gs2_layouts, only: layout use neasyf, only: neasyf_open, neasyf_read, neasyf_close use netcdf_utils, only: check_netcdf_file_precision, kind_nf # endif use mp, only: mp_abort use theta_grid, only: ntgrid use gs2_layouts, only: g_lo implicit none complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (out) :: g real, intent (in) :: scale logical, intent (in) :: has_phi, has_apar, has_bpar character (len=*), intent(in), optional :: fileopt # ifdef NETCDF logical :: this_proc_reads_fields, need_to_broadcast_fields integer, dimension(3) :: counts, start_pos real, dimension(naky) :: stmp character (run_name_size) :: file_proc character (len = 5) :: restart_layout real :: fac logical :: has_explicit_source_terms, is_one_file_per_processor, has_gexp_1 logical, save :: explicit_warning_given = .false. integer(kind_nf) :: ncid real, allocatable, dimension(:,:,:) :: tmpr, tmpi, ftmpr, ftmpi is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel) file_proc = get_file_proc(is_one_file_per_processor, fileopt) if (is_one_file_per_processor) then ncid = neasyf_open(file_proc, "r") else ncid = neasyf_open(file_proc, "r", comm=mp_comm, info=mp_info) end if call check_netcdf_file_precision (ncid) ! Check dimension sizes are compatible with the current simulation inputs call verify_compatible_dim_sizes(ncid, "theta", 2 * ntgrid + 1) call verify_compatible_dim_sizes(ncid, "sign", 2) call verify_compatible_dim_sizes(ncid, "aky", naky) call verify_compatible_dim_sizes(ncid, "akx", ntheta0) if (is_one_file_per_processor) then call verify_compatible_dim_sizes(ncid, "glo", g_lo%ulim_proc - g_lo%llim_proc + 1) else call verify_compatible_dim_sizes(ncid, "glo", g_lo%ulim_world + 1) end if if (include_explicit_source_in_restart) then has_explicit_source_terms = .true. #ifndef SHMEM has_gexp_1 = allocated(gexp_1) #else has_gexp_1 = associated(gexp_1) #endif if (has_gexp_1) & has_explicit_source_terms = has_explicit_source_terms .and. & (variable_exists(ncid, "gexp_1_r") .and. variable_exists(ncid, "gexp_1_i")) if (allocated(gexp_2)) & has_explicit_source_terms = has_explicit_source_terms .and. & (variable_exists(ncid, "gexp_2_r") .and. variable_exists(ncid, "gexp_2_i")) if (allocated(gexp_3)) & has_explicit_source_terms = has_explicit_source_terms .and. & (variable_exists(ncid, "gexp_3_r") .and. variable_exists(ncid, "gexp_3_i")) if ((.not. explicit_warning_given).and. proc0 .and. .not. has_explicit_source_terms) then write(*, '("Warning: At least one explicit source term absent in restart file.")') write(*, '(" Will not load these terms. This is probably OK unless you expect the restart file to contain these terms.")') write(*, '(" This warning will not be repeated.")') ! Update this flag to ensure we don't give this warning again explicit_warning_given = .true. end if else ! Always skip trying to get the explicit source terms if we're not ! including them has_explicit_source_terms = .false. end if call neasyf_read(ncid, "layout", restart_layout) ! Abort if the layouts don't match if(restart_layout /= layout) then call mp_abort("Incompatible layouts in restart file ("//restart_layout//") and simulation ("//layout//")",.true.) end if allocate (tmpr(2*ntgrid+1,2,g_lo%llim_proc:g_lo%ulim_alloc)) allocate (tmpi(2*ntgrid+1,2,g_lo%llim_proc:g_lo%ulim_alloc)) call zero_array(tmpr); call zero_array(tmpi) if (is_one_file_per_processor) then start_pos = [1, 1, 1] else start_pos = [1, 1, g_lo%llim_proc+1] end if counts = shape(tmpr) call neasyf_read(ncid, "gr", tmpr, start=start_pos, count=counts) call neasyf_read(ncid, "gi", tmpi, start=start_pos, count=counts) g = cmplx(tmpr, tmpi) if (include_explicit_source_in_restart .and. has_explicit_source_terms) then ! Explicit source term #ifndef SHMEM has_gexp_1 = allocated(gexp_1) #else has_gexp_1 = associated(gexp_1) #endif if (has_gexp_1) then call neasyf_read(ncid, "gexp_1_r", tmpr, start=start_pos, count=counts) call neasyf_read(ncid, "gexp_1_i", tmpi, start=start_pos, count=counts) gexp_1 = cmplx(tmpr, tmpi) endif ! Explicit source term if(allocated(gexp_2)) then call neasyf_read(ncid, "gexp_2_r", tmpr, start=start_pos, count=counts) call neasyf_read(ncid, "gexp_2_i", tmpi, start=start_pos, count=counts) gexp_2 = cmplx(tmpr, tmpi) endif ! Explicit source term if(allocated(gexp_3)) then call neasyf_read(ncid, "gexp_3_r", tmpr, start=start_pos, count=counts) call neasyf_read(ncid, "gexp_3_i", tmpi, start=start_pos, count=counts) gexp_3 = cmplx(tmpr, tmpi) endif end if deallocate(tmpr, tmpi) if (allocated(kx_shift)) then ! MR begin call neasyf_read(ncid, "kx_shift", stmp) kx_shift = stmp endif ! MR end if (allocated(theta0_shift)) then call neasyf_read(ncid, "theta0_shift", stmp) theta0_shift = stmp endif ! Decide if we are reading the fields need_to_broadcast_fields = proc_to_save_fields >= 0 this_proc_reads_fields = (iproc == proc_to_save_fields) .or. (.not. need_to_broadcast_fields) if (this_proc_reads_fields) then allocate (ftmpr(2*ntgrid+1,ntheta0,naky)) allocate (ftmpi(2*ntgrid+1,ntheta0,naky)) if (has_phi) then call neasyf_read(ncid, "phi_r", ftmpr) call neasyf_read(ncid, "phi_i", ftmpi) phinew = cmplx(ftmpr, ftmpi) end if if (has_apar) then call neasyf_read(ncid, "apar_r", ftmpr) call neasyf_read(ncid, "apar_i", ftmpi) aparnew = cmplx(ftmpr, ftmpi) end if if (has_bpar) then call neasyf_read(ncid, "bpar_r", ftmpr) call neasyf_read(ncid, "bpar_i", ftmpi) bparnew = cmplx(ftmpr, ftmpi) end if deallocate(ftmpr, ftmpi) end if call neasyf_close(ncid) if (need_to_broadcast_fields) then if (has_phi) call broadcast(phinew, proc_to_save_fields) if (has_apar) call broadcast(aparnew, proc_to_save_fields) if (has_bpar) call broadcast(bparnew, proc_to_save_fields) end if if (has_phi) call zero_array(phi) if (has_apar) call zero_array(apar) if (has_bpar) call zero_array(bpar) if (scale > 0.) then g = g*scale phinew = phinew*scale aparnew = aparnew*scale bparnew = bparnew*scale else fac = - scale/(maxval(abs(phinew))) g = g*fac phinew = phinew*fac aparnew = aparnew*fac bparnew = bparnew*fac end if # else call mp_abort("Cannot restore from restart file without netcdf.", .true.) UNUSED_DUMMY(scale); UNUSED_DUMMY(has_phi); UNUSED_DUMMY(has_apar) UNUSED_DUMMY(has_bpar); UNUSED_DUMMY(fileopt) ; g = 0.0 # endif end subroutine gs2_restore_many !> This routine writes a passed square complex array to a file !> with passed name subroutine gs2_save_response(resp, fname, code_dt, condition, header) use file_utils, only: error_unit use standard_header, only: standard_header_type #ifdef NETCDF use optionals, only: get_option_with_default use convert, only: c2r use gs2_metadata, only: create_metadata use neasyf, only: neasyf_open, neasyf_dim, neasyf_write, neasyf_close #else use file_utils, only: get_unused_unit #endif implicit none complex,dimension(:,:), intent(in) :: resp character(len=*), intent(in) :: fname real, intent(in) :: code_dt real, intent(in), optional :: condition !> Header for files with build and run information type(standard_header_type), optional, intent(in) :: header integer :: sz #ifdef NETCDF integer :: ncid real, dimension(:,:,:), allocatable :: ri_resp ! Actual value for optional `header` input type(standard_header_type) :: local_header #else integer :: unit #endif !Currently only support serial writing, but could be by any proc !so we have to make sure only one proc calls this routine !Verify we have a square array sz=size(resp(:,1)) if(sz/=size(resp(1,:))) then write(error_unit(),'("Error: gs2_save_response expects a square array input.")') return endif #ifdef NETCDF if (present(header)) then local_header = header else local_header = standard_header_type() end if ncid = neasyf_open(fname, "w") call create_metadata(ncid, "GS2 response matrix", local_header) call neasyf_dim(ncid, "ri", dim_size=2) call neasyf_dim(ncid, "ax1", dim_size=sz) call neasyf_dim(ncid, "ax2", dim_size=sz) call neasyf_write(ncid, "dt", code_dt) call neasyf_write(ncid, "condition", get_option_with_default(condition, -1.0)) !/Convert complex to ri and write allocate(ri_resp(2,sz,sz)) call c2r(resp,ri_resp) call neasyf_write(ncid, "response", ri_resp, dim_names=["ri ", "ax1", "ax2"]) deallocate(ri_resp) !/Now close the file call neasyf_close(ncid) #else !Fall back on binary output if no NETCDF !Get a free unit call get_unused_unit(unit) !Open file and write open(unit=unit,file=fname,form="unformatted") write(unit) resp close(unit) UNUSED_DUMMY(code_dt); UNUSED_DUMMY(condition); UNUSED_DUMMY(header) #endif end subroutine gs2_save_response !> This routine reads a square complex array from a file !! with passed name subroutine gs2_restore_response(resp, fname, code_dt, condition) use file_utils, only: error_unit #ifdef NETCDF use neasyf, only: neasyf_open, neasyf_close, neasyf_read use convert, only: r2c #else use file_utils, only: get_unused_unit #endif implicit none complex, dimension(:,:), intent(out) :: resp character(len=*), intent(in) :: fname real, intent(out) :: code_dt real, intent(out), optional :: condition integer :: sz #ifdef NETCDF integer :: ncid real, dimension(:,:,:), allocatable :: ri_resp #else integer :: unit #endif !Currently only support serial reading, but could be by any proc !so we have to make sure only one proc calls this routine !Verify we have a square array sz=size(resp(:,1)) if(sz/=size(resp(1,:))) then write(error_unit(),'("Error: gs2_restore_response expects a square array output.")') return endif #ifdef NETCDF ncid = neasyf_open(fname, "r") !/Read and convert ri to complex allocate(ri_resp(2,sz,sz)) call neasyf_read(ncid, "response", ri_resp) call r2c(resp,ri_resp) deallocate(ri_resp) !/Get the time step call neasyf_read(ncid, "dt", code_dt) !/Get the condition number if (present(condition)) call neasyf_read(ncid, "condition", condition) !/Now close the file call neasyf_close(ncid) #else !Fall back on binary output if no NETCDF !Get a free unit call get_unused_unit(unit) !Open file and write open(unit=unit,file=fname,form="unformatted") read(unit) resp close(unit) !Set outputs to something invalid code_dt = -1 ; condition = 0.0 #endif end subroutine gs2_restore_response !> Initialises a file for saving output of eigensolver to netcdf subroutine init_eigenfunc_file(fname, IDs) #ifdef NETCDF use theta_grid, only: theta use neasyf, only: neasyf_open, neasyf_dim #endif implicit none character(len=*), intent(in) :: fname type(EigNetcdfID), intent(inout) :: IDs !Set nconv counter to 0 IDs%nconv_count=0 #ifdef NETCDF IDs%ncid = neasyf_open(fname, "w") call neasyf_dim(IDs%ncid, "ri", dim_size=2) call neasyf_dim(IDs%ncid, "theta", theta) call neasyf_dim(IDs%ncid, "conv", unlimited=.true.) #else UNUSED_DUMMY(fname) #endif end subroutine init_eigenfunc_file !> Add an eigenpairs data to file subroutine add_eigenpair_to_file(eval, has_phi, has_apar, has_bpar, IDs, my_conv) #ifdef NETCDF use fields_arrays, only: phinew, aparnew, bparnew use convert, only: c2r use theta_grid, only: ntgrid use neasyf, only: neasyf_write use optionals, only: get_option_with_default #endif complex, intent(in) :: eval !Note just use fields to get eigenvectors real, intent(in), optional :: my_conv type(EigNetcdfID), intent(inout) :: IDs logical, intent(in) :: has_phi, has_apar, has_bpar #ifdef NETCDF real, dimension(2) :: ri_omega real, dimension(:,:), allocatable :: ri_field integer, dimension(3) :: start_field, count_field real :: local_conv character(*), dimension(*), parameter :: eigen_field_dim_names = [character(len=5)::"ri", "theta", "conv"] !First increment counter IDs%nconv_count=IDs%nconv_count+1 start_field = [1, 1, IDs%nconv_count] count_field = [2, 2*ntgrid+1, 1] local_conv = get_option_with_default(my_conv, real(IDs%nconv_count)) call neasyf_write(IDs%ncid, "conv", local_conv, dim_names=["conv"], & start=[IDs%nconv_count]) !/Omega ri_omega(1)=real(eval) ri_omega(2)=aimag(eval) call neasyf_write(IDs%ncid, "omega", ri_omega, dim_names=["ri ", "conv"], & start=[1, IDs%nconv_count], count=[2,1]) !/Fields allocate(ri_field(2,2*ntgrid+1)) if(has_phi)then call c2r(phinew(:,1,1),ri_field) call neasyf_write(IDs%ncid, "phi", ri_field, dim_names=eigen_field_dim_names, & start=start_field, count=count_field) endif if(has_apar)then call c2r(aparnew(:,1,1),ri_field) call neasyf_write(IDs%ncid, "apar", ri_field, dim_names=eigen_field_dim_names, & start=start_field, count=count_field) endif if(has_bpar)then call c2r(bparnew(:,1,1),ri_field) call neasyf_write(IDs%ncid, "bpar", ri_field, dim_names=eigen_field_dim_names, & start=start_field, count=count_field) endif deallocate(ri_field) #else UNUSED_DUMMY(eval); UNUSED_DUMMY(has_phi); UNUSED_DUMMY(has_apar); UNUSED_DUMMY(has_bpar) UNUSED_DUMMY(IDs); UNUSED_DUMMY(my_conv) #endif end subroutine add_eigenpair_to_file !> Close the eigenfunction file subroutine finish_eigenfunc_file(IDs) #ifdef NETCDF use neasyf, only: neasyf_close #endif implicit none type(EigNetcdfID), intent(inout) :: IDs #ifdef NETCDF call neasyf_close(IDs%ncid) #else UNUSED_DUMMY(IDs) #endif end subroutine finish_eigenfunc_file !> This function checks to see if we can create a file with name !! <restart_file>//<SomeSuffix> if not then our restarts are not !! going to be possible and we return false. Can also be used to check !! that we can read from the restart file (which assumes it exists). function restart_writable(read_only, my_file, error_message) use mp, only: proc0, broadcast use file_utils, only: get_unused_unit use optionals, only: get_option_with_default implicit none !> If present and true, only check that files can be read logical, intent(in), optional :: read_only !> An optional specific filename to check character(len=*),intent(in),optional::my_file !> Error message returned from `open` if there was a problem character(len=:), allocatable, optional, intent(out) :: error_message character(len=*), parameter :: SuffixTmp = '.ThisIsATestFile' character(9) :: open_mode character(6) :: close_mode character(run_name_size) :: local_file, filename logical :: restart_writable integer :: unit, ierr character(len=200) :: io_err_msg ierr=-200 local_file = trim(get_option_with_default(my_file, restart_file)) ! On proc0 try to open tmp file for writing if (proc0)then call get_unused_unit(unit) ! Default open and close modes open_mode = 'readwrite' close_mode = 'delete' ! If we want to test write capability then do it with an unusual ! file name to prevent clobber filename = trim(local_file) // SuffixTmp ! Set filemode to READ if read_only=T if (get_option_with_default(read_only, .false.)) then ! If we're only checking a file can be read, then don't ! delete the file after we're done open_mode = 'read' close_mode = 'keep' ! If checking readonly then we need to make sure we try to ! read from an existing file filename = local_file end if open(unit=unit, file=trim(filename), & iostat=ierr, action=open_mode, & iomsg=io_err_msg) ! If open was successful then we can close the file and delete it if (ierr == 0) close(unit=unit, status=close_mode) endif ! Now make sure everyone knows the answer call broadcast(ierr) restart_writable = (ierr == 0) if (.not. present(error_message)) return if (restart_writable) then error_message = "" else error_message = trim(io_err_msg) end if end function restart_writable !> FIXME : Add documentation subroutine init_gs2_save end subroutine init_gs2_save !> Sets the base of the restart file to be used when !> writing/reading files. subroutine set_restart_file (file) use mp, only: proc0 use file_utils, only: error_unit implicit none character(len=*), intent (in) :: file if (proc0 .and. len_trim(file) > run_name_size) then write( error_unit(), '("Argument to set_restart_file exceeds restart_file size so truncating")') end if restart_file = file(1:min(len_trim(file), run_name_size)) end subroutine set_restart_file !> FIXME : Add documentation subroutine finish_gs2_save end subroutine finish_gs2_save !> FIXME : Add documentation subroutine init_dt (delt0, delt1, delt2, delt_max, not_set_value) # ifdef NETCDF use mp, only: proc0, broadcast use optionals, only: get_option_with_default use neasyf, only: neasyf_open, neasyf_read, neasyf_close use netcdf_utils, only: kind_nf # else use mp, only: mp_abort # endif implicit none real, intent (in out) :: delt0, delt1, delt2, delt_max real, intent(in), optional :: not_set_value # ifdef NETCDF character (run_name_size) :: file_proc real :: not_set_value_to_use logical :: is_one_file_per_processor ! NetCDF handles integer(kind_nf) :: ncid real, dimension(0:3) :: time_values if (proc0) then is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel) file_proc = get_file_proc(is_one_file_per_processor) ncid = neasyf_open(file_proc, "r") not_set_value_to_use = get_option_with_default(not_set_value, -1.0) time_values = not_set_value_to_use ! Note unlike the explicit source terms all three time steps should always ! be available in the restart file so we don't silence the error messages here. ! The only situation we are likely to come across where the delt1 and delt2 ! values aren't available is where we are trying to read an old restart file. ! This will then lead to the error/warning being displayed but the code should ! carry on as intended and the missing steps will be set to a special value to ! indicate they have not been set yet. if (variable_exists(ncid, "delt0")) call neasyf_read(ncid, "delt0", time_values(0)) if (variable_exists(ncid, "delt1")) call neasyf_read(ncid, "delt1", time_values(1)) if (variable_exists(ncid, "delt2")) call neasyf_read(ncid, "delt2", time_values(2)) if (variable_exists(ncid, "delt_max")) call neasyf_read(ncid, "delt_max", time_values(3)) call neasyf_close(ncid) endif call broadcast(time_values) delt0 = time_values(0) ; delt1 = time_values(1) ; delt2 = time_values(2) delt_max = time_values(3) # else call mp_abort("Cannot load initial dt from restart without netcdf.", .true.) UNUSED_DUMMY(delt0); UNUSED_DUMMY(delt1); UNUSED_DUMMY(delt2); UNUSED_DUMMY(delt_max) UNUSED_DUMMY(not_set_value) # endif end subroutine init_dt !> FIXME : Add documentation subroutine init_vnm (vnm) # ifdef NETCDF use mp, only: proc0, broadcast use neasyf, only: neasyf_open, neasyf_read, neasyf_close use netcdf_utils, only: kind_nf # else use mp, only: mp_abort # endif implicit none real, dimension(2), intent (in out) :: vnm # ifdef NETCDF character (run_name_size) :: file_proc logical :: is_one_file_per_processor ! NetCDF handles integer(kind_nf) :: ncid if (proc0) then is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel) file_proc = get_file_proc(is_one_file_per_processor) ncid = neasyf_open (file_proc, "r") call neasyf_read(ncid, "vnm1", vnm(1)) call neasyf_read(ncid, "vnm2", vnm(2)) call neasyf_close(ncid) endif call broadcast (vnm) # else call mp_abort("Cannot load vnm from restart without netcdf.", .true.) UNUSED_DUMMY(vnm) # endif end subroutine init_vnm !> FIXME : Add documentation !! !! @note This routine gets a_ant and b_ant for proc 0 only! subroutine init_ant_amp (a_ant, b_ant, nk_stir) # ifdef NETCDF use constants, only: zi use neasyf, only: neasyf_open, neasyf_read, neasyf_close use mp, only: proc0, broadcast, mp_abort use netcdf_utils, only: kind_nf # else use mp, only: mp_abort # endif implicit none complex, dimension(:), intent (in out) :: a_ant, b_ant integer, intent (in) :: nk_stir # ifdef NETCDF character (run_name_size) :: file_proc logical :: is_one_file_per_processor ! NetCDF handles integer(kind_nf) :: ncid real, dimension(nk_stir) :: atmp if (proc0) then is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel) file_proc = get_file_proc(is_one_file_per_processor) ncid = neasyf_open (file_proc, "r") call verify_compatible_dim_sizes(ncid, "nk_stir", nk_stir) atmp = 0. call neasyf_read(ncid, "a_ant_r", atmp) a_ant = atmp call neasyf_read(ncid, "a_ant_i", atmp) a_ant = a_ant + zi * atmp call neasyf_read(ncid, "b_ant_r", atmp) b_ant = atmp call neasyf_read(ncid, "b_ant_i", atmp) b_ant = b_ant + zi * atmp call neasyf_close(ncid) end if call broadcast(a_ant) call broadcast(b_ant) # else call mp_abort("Cannot load antenna amplitudes from restart without netcdf.", .true.) UNUSED_DUMMY(a_ant); UNUSED_DUMMY(b_ant); UNUSED_DUMMY(nk_stir) # endif end subroutine init_ant_amp !> FIXME : Add documentation subroutine read_t0_from_restart_file (tstart) # ifdef NETCDF use mp, only: proc0, broadcast use neasyf, only: neasyf_open, neasyf_read, neasyf_close use netcdf_utils, only: kind_nf # else use mp, only: mp_abort # endif implicit none real, intent (in out) :: tstart # ifdef NETCDF character (run_name_size) :: file_proc logical :: is_one_file_per_processor ! NetCDF handles integer(kind_nf) :: ncid if (proc0) then is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel) file_proc = get_file_proc(is_one_file_per_processor) ncid = neasyf_open(file_proc, "r") call neasyf_read(ncid, "t0", tstart) call neasyf_close(ncid) end if call broadcast (tstart) # else call mp_abort("Cannot read t0 from restart file in build with no netcdf.", .true.) UNUSED_DUMMY(tstart) # endif end subroutine read_t0_from_restart_file end module gs2_save