module checksum_braginskii_m !! Some basic diagnostics, execudeted during runtime use MPI use comm_handler_m, only : comm_handler_t use screen_io_m, only : get_stdout use parallelisation_setup_m, only : is_rank_info_writer use error_handling_grillix_m, only: handle_error use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST use precision_grillix_m, only : GP, GP_LARGEST, MPI_GP use variable_m, only : variable_t use mesh_cart_m, only : mesh_cart_t implicit none public :: checksum_braginskii private :: compute_checksum_var_local contains subroutine checksum_braginskii(comm_handler, mesh, tau, ne, pot, vort, upar, jpar, apar, te, ti) !! Computes and writes checksum of a plasma state !! The checksum does not have any physical meaning, but is used e.g in CI/CD to ensure code integrity type(comm_handler_t), intent(in) :: comm_handler type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane real(GP), intent(in) :: tau !! Time type(variable_t), intent(in) :: ne !! Electron density type(variable_t), intent(in) :: pot !! Electrostatic potential type(variable_t), intent(in) :: vort !! Generalised vorticity type(variable_t), intent(in) :: upar !! Parallel ion velocity type(variable_t), intent(in) :: jpar !! Parallel current type(variable_t), intent(in) :: apar !! Electromagnetic potential type(variable_t), intent(in) :: te !! Electron temperature type(variable_t), intent(in) :: ti !! Ion temperature logical :: fexist integer :: plane, ioerr character(len=256) :: io_errmsg character(len=13) :: fname real(GP) :: mean_ne, l1nrm_ne, maxv_ne, minv_ne real(GP) :: mean_pot, l1nrm_pot, maxv_pot, minv_pot real(GP) :: mean_vort, l1nrm_vort, maxv_vort, minv_vort real(GP) :: mean_upar, l1nrm_upar, maxv_upar, minv_upar real(GP) :: mean_jpar, l1nrm_jpar, maxv_jpar, minv_jpar real(GP) :: mean_apar, l1nrm_apar, maxv_apar, minv_apar real(GP) :: mean_te, l1nrm_te, maxv_te, minv_te real(GP) :: mean_ti, l1nrm_ti, maxv_ti, minv_ti real(GP) :: plane_fac ! Checksum within plane ----------------------------------------------- call compute_checksum_var_local(mesh, ne, mean_ne, l1nrm_ne, maxv_ne, minv_ne) call compute_checksum_var_local(mesh, pot, mean_pot, l1nrm_pot, maxv_pot, minv_pot) call compute_checksum_var_local(mesh, vort, mean_vort, l1nrm_vort, maxv_vort, minv_vort) call compute_checksum_var_local(mesh, upar, mean_upar, l1nrm_upar, maxv_upar, minv_upar) call compute_checksum_var_local(mesh, jpar, mean_jpar, l1nrm_jpar, maxv_jpar, minv_jpar) call compute_checksum_var_local(mesh, apar, mean_apar, l1nrm_apar, maxv_apar, minv_apar) call compute_checksum_var_local(mesh, te, mean_te, l1nrm_te, maxv_te, minv_te) call compute_checksum_var_local(mesh, ti, mean_ti, l1nrm_ti, maxv_ti, minv_ti) ! Reduce over planes with plane dependent factor ---------------------- plane = comm_handler%get_plane() plane_fac = 1.0_GP*plane**2 call planes_reduction(comm_handler, plane_fac, mean_ne, l1nrm_ne, maxv_ne, maxv_ne) call planes_reduction(comm_handler, plane_fac, mean_pot, l1nrm_pot, maxv_pot, minv_pot) call planes_reduction(comm_handler, plane_fac, mean_vort, l1nrm_vort, maxv_vort, minv_vort) call planes_reduction(comm_handler, plane_fac, mean_upar, l1nrm_upar, maxv_upar, minv_upar) call planes_reduction(comm_handler, plane_fac, mean_jpar, l1nrm_jpar, maxv_jpar, minv_jpar) call planes_reduction(comm_handler, plane_fac, mean_apar, l1nrm_apar, maxv_apar, minv_apar) call planes_reduction(comm_handler, plane_fac, mean_te, l1nrm_te, maxv_te, minv_te) call planes_reduction(comm_handler, plane_fac, mean_ti, l1nrm_ti, maxv_ti, minv_ti) ! Write out result to file for integration tests ------------------------------------------ if (is_rank_info_writer) then fname = 'checksums.out' inquire(file=fname, exist = fexist) if (fexist) then open(unit=58, file=fname, status='old', position='append', & action='write', iostat=ioerr, iomsg=io_errmsg) else open(unit=58, file=fname, status='new', & action='write', iostat=ioerr, iomsg=io_errmsg) endif if (ioerr /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif write(58,605)tau write(58,605)mean_ne write(58,605)l1nrm_ne write(58,605)maxv_ne write(58,605)minv_ne write(58,605)mean_pot write(58,605)l1nrm_pot write(58,605)maxv_pot write(58,605)minv_pot write(58,605)mean_vort write(58,605)l1nrm_vort write(58,605)maxv_vort write(58,605)minv_vort write(58,605)mean_upar write(58,605)l1nrm_upar write(58,605)maxv_upar write(58,605)minv_upar write(58,605)mean_jpar write(58,605)l1nrm_jpar write(58,605)maxv_jpar write(58,605)minv_jpar write(58,605)mean_apar write(58,605)l1nrm_apar write(58,605)maxv_apar write(58,605)minv_apar write(58,605)mean_te write(58,605)l1nrm_te write(58,605)maxv_te write(58,605)minv_te write(58,605)mean_ti write(58,605)l1nrm_ti write(58,605)maxv_ti write(58,605)minv_ti 605 FORMAT(ES20.12E3) close(58) endif end subroutine subroutine compute_checksum_var_local(mesh, var, mean, l1nrm, maxv, minv) !! Computes some checksums per plane for a given variable type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane type(variable_t), intent(in) :: var !! Variable real(GP), intent(out) :: mean !! Mean real(GP), intent(out) :: l1nrm !! L1-norm real(GP), intent(out) :: maxv !! Maximum real(GP), intent(out) :: minv !! Minimum integer :: l mean = 0.0_GP l1nrm = 0.0_GP maxv = -GP_LARGEST minv = GP_LARGEST !$OMP PARALLEL REDUCTION(+: mean, l1nrm) Reduction(max: maxv) Reduction(min: minv) !$OMP DO do l = 1, mesh%get_n_points() mean = mean + var%vals(l) l1nrm = l1nrm + abs(var%vals(l)) maxv = max(maxv, var%vals(l)) minv = min(minv, var%vals(l)) enddo !$OMP END DO !$OMP END PARALLEL mean = mean/(1.0_GP*mesh%get_n_points()) l1nrm = l1nrm/(1.0_GP*mesh%get_n_points()) end subroutine subroutine planes_reduction(comm_handler, plane_fac, mean, l1nrm, maxv, minv) !! Performs reduction over planes with plane dependent factor type(comm_handler_t), intent(in) :: comm_handler !! Communicators real(GP), intent(in) :: plane_fac !! plane dependent factor real(GP), intent(inout) :: mean !! Mean (on output summed over planes) real(GP), intent(inout) :: l1nrm !! L1-norm (on output summed over planes) real(GP), intent(inout) :: maxv !! Maximum (on output summed over planes, not maximum over planes!!!) real(GP), intent(inout) :: minv !! Minimum (on output summed over planes, not maximum over planes!!!) integer :: ierr mean = mean*plane_fac call MPI_Allreduce(MPI_IN_PLACE, mean, 1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr) l1nrm = l1nrm*plane_fac call MPI_Allreduce(MPI_IN_PLACE, l1nrm, 1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr) maxv = maxv*plane_fac call MPI_Allreduce(MPI_IN_PLACE, maxv, 1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr) minv = minv*plane_fac call MPI_Allreduce(MPI_IN_PLACE, minv, 1, MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr) end subroutine end module