module diagnostics_template_m !! Performs diagnostics on state use MPI use parallelisation_setup_m, only : is_rank_info_writer use screen_io_m, only : get_stdout use precision_grillix_m, only : GP, GP_LARGE, MPI_GP use error_handling_grillix_m, only: handle_error use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use comm_handler_m, only : comm_handler_t use static_data_m, only : masks_cano, masks_stag, map use state_template_m, only : np_max, tau, dens, parmom implicit none public :: diagnostics_template private :: state_norms contains subroutine diagnostics_template(comm_handler, diag_file) !! Diagnostics for template model !! Yet just the L1, L2 global minimum and global maximum values of dynamical fields are computed type(comm_handler_t), intent(in) :: comm_handler !! Communication handler character(len=*), intent(in), optional :: diag_file !! File where diagnostics can be written to real(GP) :: dens_avg, dens_l2, dens_min, dens_max, parmom_avg, parmom_l2, parmom_min, parmom_max logical :: fexist character(:), allocatable :: fstat integer :: io_unit, io_error character(len=256) :: io_errmsg integer :: len_diag_file call state_norms(comm_handler, dens, .false., dens_avg, dens_l2, dens_min, dens_max) call state_norms(comm_handler, parmom, .true., parmom_avg, parmom_l2, parmom_min, parmom_max) ! Write diagnostics to screen if (is_rank_info_writer) then write(get_stdout(),701)tau, dens_avg, dens_l2, dens_min, dens_max, parmom_avg, parmom_l2, parmom_min, parmom_max 701 FORMAT(3X,'tau = ',ES20.12E3, / & 3X,'dens_avg = ',ES20.12E3, / & 3X,'dens_l2 = ',ES20.12E3, / & 3X,'dens_min = ',ES20.12E3, / & 3X,'dens_max = ',ES20.12E3, / & 3X,'parmom_avg = ',ES20.12E3, / & 3X,'parmom_l2 = ',ES20.12E3, / & 3X,'parmom_min = ',ES20.12E3, / & 3X,'parmom_max = ',ES20.12E3 ) endif ! Write diagnostics to file len_diag_file = 0 if (present(diag_file)) then len_diag_file = len(trim(diag_file)) endif if (len_diag_file > 0) then if (is_rank_info_writer) then write(get_stdout(),*)'Appending diagnostics to file: ', diag_file inquire(file=diag_file, exist=fexist) if (fexist) then fstat = 'old' else fstat = 'new' endif open(newunit=io_unit, file=trim(diag_file), status=fstat, action='write', position="append", iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_OTHER, __LINE__, __FILE__) endif write(io_unit,702) tau, dens_avg, dens_l2, dens_min, dens_max, parmom_avg, parmom_l2, parmom_min, parmom_max 702 FORMAT(9(ES20.12E3,X)) close(io_unit) endif endif end subroutine subroutine state_norms(comm_handler, var, is_stag, nrm_avg, nrm_l2, nrm_min, nrm_max) type(comm_handler_t), intent(in) :: comm_handler !! Communication handler real(GP), dimension(np_max), intent(in) :: var !! Dynamical field logical, intent(in) :: is_stag !! If false, var is located on canonical mesh, !! if true, var is located on staggered mesh real(GP), intent(out) :: nrm_avg !! Average value of field real(GP), intent(out) :: nrm_l2 !! L2 norm of field real(GP), intent(out) :: nrm_min !! Global minimum of field real(GP), intent(out) :: nrm_max !! Global maximum of field integer :: l, ierr real(GP) :: vol nrm_avg = 0.0_GP nrm_l2 = 0.0_GP vol = 0.0_GP nrm_min = GP_LARGE nrm_max = -GP_LARGE if (is_stag) then !$omp parallel default(none) private(l) shared(np_max, masks_stag, map, var) & !$omp reduction(+:nrm_avg, nrm_l2, vol) reduction(min:nrm_min) reduction(max:nrm_max) !$omp do do l = 1, np_max if (masks_stag%is_inner(l) .or. masks_stag%is_boundary(l) .or. masks_stag%is_ghost(l)) then nrm_avg = nrm_avg + var(l) * map%fluxbox_vol_stag(l) nrm_l2 = nrm_l2 + var(l)**2 * map%fluxbox_vol_stag(l) nrm_min = min(nrm_min, var(l)) nrm_max = max(nrm_max, var(l)) vol = vol + map%fluxbox_vol_stag(l) endif enddo !$omp end do !$omp end parallel else !$omp parallel default(none) private(l) shared(np_max, masks_cano, map, var) & !$omp reduction(+:nrm_avg, nrm_l2, vol) reduction(min:nrm_min) reduction(max:nrm_max) !$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 nrm_avg = nrm_avg + var(l) * map%fluxbox_vol_cano(l) nrm_l2 = nrm_l2 + var(l)**2 * map%fluxbox_vol_cano(l) nrm_min = min(nrm_min, var(l)) nrm_max = max(nrm_max, var(l)) vol = vol + map%fluxbox_vol_cano(l) endif enddo !$omp end do !$omp end parallel endif call MPI_allreduce(MPI_IN_PLACE, nrm_avg, 1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr) call MPI_allreduce(MPI_IN_PLACE, nrm_l2, 1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr) call MPI_allreduce(MPI_IN_PLACE, vol, 1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr) call MPI_allreduce(MPI_IN_PLACE, nrm_min, 1, MPI_GP, MPI_MIN, comm_handler%get_comm_planes(), ierr) call MPI_allreduce(MPI_IN_PLACE, nrm_max, 1, MPI_GP, MPI_MAX, comm_handler%get_comm_planes(), ierr) nrm_avg = nrm_avg / vol nrm_l2 = sqrt(nrm_l2 / vol) end subroutine end module