state_template_m.f90 Source File

module state_template_m
    !! Holds current state and handles its I/O via netcdf
    use MPI
    use precision_grillix_m, only : GP
    use comm_handler_m, only : comm_handler_t
    use static_data_m, only : mesh_cano, mesh_stag, masks_stag, masks_cano, map
    use params_template_m, only : mms_on
    use mpi_mapping_auxiliaries_m, only : getdata_fwdbwdplane
    implicit none

    integer, public :: tstep
    !! Current time-step of state
    real(GP), public :: tau
    !! Current time of state
    integer, protected :: isnaps = 0
    !! Current snapshot number
    integer, protected :: np_max
    !! Maximum number of points per plane, across all planes and stag/cano

    real(GP), public, allocatable, dimension(:), target :: dens
    !! Density
    real(GP), public, allocatable, dimension(:) :: parmom
    !! Parallel momentum
    real(GP), public, allocatable, dimension(:) :: upar
    !! Parallel ion velocity
    !! Derived field: upar = parmom / dens (on staggered mesh)
    real(GP), public, allocatable, dimension(:) :: pion
    !! Ion pressure
    real(GP), public, allocatable, dimension(:) :: tion
    !! Ion temperature
    !! Derived field: tion = pion / dens

    public :: alloc_state_template
    public :: update_derived_fields
    public :: write_template_state_netcdf
    public :: read_template_state_netcdf

    interface
        module subroutine write_template_state_netcdf(comm_handler, dirname, filename)
            !! Writes state to netcdf file
            !! If filename is not present, it increases isnaps by one and writes to `snaps_<isnaps>.nc`
            !! If filename is present, it writes to filename, without increasing isnaps. 
            !! This is useful for writing e.g. error states
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in) , optional :: dirname
            !! Directory, where snapshots files will be written
            character(len=*), intent(in) , optional :: filename
            !! Filename, where state is written. The filename excludes the directory name,
            !! i.e. the final path is dirname/filename
        end subroutine

        module subroutine read_template_state_netcdf(comm_handler, dirname, filename)
            !! Reads state from netcdf file
            !! If filename is not present, it reads last from `snaps_<isnaps>.nc` with last available isnaps
            !! and sets isnaps to this value
            !! If filename is present it reads from this file and sets isnaps to 0
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communication handler
            character(len=*), intent(in) , optional :: dirname
            !! Directory, where snapshots files will be read
            character(len=*), intent(in) , optional :: filename
            !! Filename, where state is read. The filename excludes the directory name,
            !! i.e. the final path is dirname/filename
        end subroutine
    end interface

contains

    subroutine alloc_state_template(comm_handler)
        !! Allocates memory for state fields
        type(comm_handler_t), intent(in) :: comm_handler
        !! Communication handler

        integer :: ierr, l
        integer :: np_cano, np_stag

        ! Find maximum number of mesh points across all planes
        np_cano = mesh_cano%get_n_points()
        np_stag = mesh_stag%get_n_points()
        np_max = max(np_cano, np_stag)

        call MPI_allreduce(MPI_IN_PLACE, np_max, 1, MPI_INTEGER, MPI_MAX, &
                           comm_handler%get_comm_planes(), ierr)

        allocate(dens(np_max))
        allocate(parmom(np_max))
        allocate(upar(np_max))
        allocate(pion(np_max))
        allocate(tion(np_max))

        ! First touch
        !$omp parallel default(none) private(l) shared(np_max, dens, parmom, upar, pion, tion)
        !$omp do
        do l = 1, np_max
            dens(l) = 0.0_GP
            parmom(l) = 0.0_GP
            upar(l) = 0.0_GP
            pion(l) = 0.0_GP
            tion(l) = 0.0_GP
        enddo
        !$omp end do
        !$omp end parallel

    end subroutine

    subroutine update_derived_fields(comm_handler)
        !! Updates derived fields
        !! upar = parmom /dens_(stag)
        !! tion = pion / dens_(cano)
        type(comm_handler_t), intent(in) :: comm_handler
        !! Communication handler

        integer :: l, nrecv
        real(GP), allocatable, dimension(:) :: halo_buf
        real(GP), dimension(np_max) :: dens_stag

        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, dens, nrecv, halo_buf)
        call map%lift_cano_to_stag(dens, halo_buf, dens_stag)
        deallocate(halo_buf)

        !$omp parallel default(none) private(l) shared(np_max, masks_stag, masks_cano, upar, &
        !$omp                                          parmom, dens_stag, tion, pion, dens)
        !$omp do
        do l = 1, np_max
            if (masks_cano%is_inner(l) .or. masks_cano%is_boundary(l) .or. masks_cano%is_ghost(l)) then
                tion(l) = pion(l) / dens(l)
            endif
            if (masks_stag%is_inner(l) .or. masks_stag%is_boundary(l) .or. masks_stag%is_ghost(l)) then
                upar(l) = parmom(l) / dens_stag(l)
            endif
        enddo
        !$omp end do
        !$omp end parallel

    end subroutine

end module