snapshots_helpers_m.f90 Source File


Contents


Source Code

module snapshots_helpers_m
    !! Module containing useful routines for handling snapshots
    use netcdf
    use precision_grillix_m, only : GP
    use comm_handler_m, only : comm_handler_t
    use error_handling_grillix_m, only: handle_error, error_info_t, handle_error_netcdf
    use mesh_cart_m, only: mesh_cart_t
    use variable_m, only : variable_t
    implicit none
    public :: snaps_average
contains 

    subroutine snaps_average(comm_handler, mesh, dirpath, nsnaps_start, nsnaps_end, var_name, res)
        !! Averages an variable in the snapshots from nsnaps_start until nsnaps_end
        type(comm_handler_t), intent(in) :: comm_handler
        !! MPI communication handler
        type(mesh_cart_t), intent(in) :: mesh
        !! Mesh within poloidal plane
        !! Canonical or staggered. Should agree with var_name and res
        character(len=*), intent(in)  :: dirpath
        !! Directory, where snapshot files are read
        integer, intent(in) :: nsnaps_start
        !! Snapshot number where the average period starts
        integer, intent(in) :: nsnaps_end
        !! Snapshot number where the average period ends
        character(len=*), intent(in) :: var_name
        !! Name of the variable to be averaged, must be contained in dirpath
        real(GP), dimension(mesh%get_n_points()), intent(inout) :: res
        !! Stores results of average

        type(variable_t) :: var_temp
        real(GP) :: step
        integer :: isnaps, l
        character(len=5) :: isnaps_c
        character(len=:), allocatable :: snapsfile
        integer :: nf90_stat, nf90_id
        
        res = 0.0_GP

        do isnaps = nsnaps_start, nsnaps_end 
            step = real(isnaps - nsnaps_start + 1, GP)
            ! Read variable from snapshots
            write(isnaps_c,'(I5.5)')isnaps
            snapsfile = dirpath // '/snaps_t' // isnaps_c // '.nc' 
            nf90_stat = nf90_open(snapsfile, NF90_NETCDF4, nf90_id )
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            call var_temp%read_netcdf(comm_handler, nf90_id, var_name)   
            ! Time average
            !$omp parallel default(none) private(l) &
            !$omp shared(res, var_temp, step, mesh)
            !$omp do
            do l = 1, mesh%get_n_points()
                res(l) = res(l) * (step - 1.0_GP) / step + var_temp%vals(l) / step
            enddo
            !$omp end do
            !$omp end parallel
        enddo

    end subroutine

end module