module mms_zhdanov_m !! Method of Manufactured Solutions for Zhdanov model !! Note: Usage of camel case conventiuon for consistency in naming with mathematica use MPI use precision_grillix_m, only : GP, GP_EPS, MPI_GP use constants_m, only : TWO_PI use parallelisation_setup_m, only : is_rank_info_writer use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST, GRILLIX_ERR_OTHER use comm_handler_m, only : comm_handler_t use screen_io_m, only : get_stdout use equilibrium_m, only : equilibrium_t use circular_equilibrium_m, only : circular_t use slab_equilibrium_m, only : slab_t use commons_misc_m, only : error_ananum use mesh_cart_m, only : mesh_cart_t use variable_m, only : variable_t use params_zhdanov_species_m, only : params_zhdanov_species_t implicit none public :: mms_diagnostics_zhdanov public :: MmsSolDens public :: MmsSourceDens private :: MmsSolDensCircular private :: MmsSolDensSlab private :: MmsSourceDensCircular private :: MmsSourceDensSlab interface real(GP) module function MmsSolDensCircular(equi, species, x, y, z, t) !! MMS solution of density for specific species in circular geometry class(equilibrium_t) :: equi !! Equilibrium integer, intent(in) :: species !! Species number real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: z !! Axial coordinate z (represented by toroidal coordinate) real(GP), intent(in) :: t !! Time end function real(GP) module function MmsSolDensSlab(equi, species, x, y, z, t) !! MMS solution of density for specific species in circular geometry class(equilibrium_t) :: equi !! Equilibrium integer, intent(in) :: species !! Species number real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: z !! Axial coordinate z (represented by toroidal coordinate) real(GP), intent(in) :: t !! Time end function real(GP) module function MmsSourceDensCircular(equi, species,& x, y, z, t, params_species) !! MMS solution of density for specific species class(equilibrium_t) :: equi !! Equilibrium integer, intent(in) :: species !! Species number real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: z !! Axial coordinate z (represented by toroidal coordinate) real(GP), intent(in) :: t !! Time type(params_zhdanov_species_t) :: params_species !! Species parameters end function real(GP) module function MmsSourceDensSlab(equi, species,& x, y, z, t, params_species) !! MMS solution of density for specific species class(equilibrium_t) :: equi !! Equilibrium integer, intent(in) :: species !! Species number real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: z !! Axial coordinate z (represented by toroidal coordinate) real(GP), intent(in) :: t !! Time type(params_zhdanov_species_t) :: params_species !! Species parameters end function end interface contains subroutine mms_diagnostics_zhdanov(comm_handler, equi, mesh,& tau, dens_species) !! Prints information on numerical errors of numerical solution type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane real(GP), intent(in) :: tau !! Time type(variable_t), intent(in) :: dens_species !! Density of species integer :: l, s, ierr, ioerr real(GP) :: x, y, phi real(GP), dimension(mesh%get_n_points()) :: mms_sol_dens real(GP) :: dens_species_nrm2, dens_species_nrmsup,& dens_species_err2, dens_species_errsup, dens_species_err2nrm integer, dimension(comm_handler%get_nspecies()) :: gath_species real(GP), dimension(comm_handler%get_nspecies()) ::& gath_nrm2, gath_err2nrm, gath_errsup character(len=19) :: jobfile character(len=256) :: io_errmsg logical :: fexist ! Get MMS solution phi = mesh%get_phi() !$OMP PARALLEL PRIVATE(l, x, y) !$OMP DO do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) mms_sol_dens(l) = MmsSolDens(equi, comm_handler%get_species(),& x, y, phi, tau) enddo !$OMP END DO !$OMP END PARALLEL call error_ananum(comm_handler, mesh%get_n_points(),& mms_sol_dens, dens_species%vals,& dens_species_nrm2, dens_species_nrmsup,& dens_species_err2, dens_species_errsup) dens_species_err2nrm = dens_species_err2 / (dens_species_nrm2 + GP_EPS) if (is_rank_info_writer) then write(get_stdout(),*)'MMS diagnostics ------------------------------------------' endif call MPI_barrier(comm_handler%get_comm_cart(), ierr) if (comm_handler%get_plane() == 0) then write(get_stdout(),204)comm_handler%get_species(), tau, & dens_species_nrm2, dens_species_err2nrm, dens_species_errsup 204 FORMAT('species = ',I3,X,'tau = ',ES14.6E3,X,'nrm2 = ',ES14.6E3,X,& 'err2nrm = ',ES14.6E3,X,'errsup',ES14.6E3) endif call MPI_barrier(comm_handler%get_comm_cart(), ierr) if (is_rank_info_writer) then write(get_stdout(),*)'----------------------------------------------------------' endif ! Gather results on master rank call MPI_gather(comm_handler%get_species(), 1, MPI_INTEGER, & gath_species, 1, MPI_INTEGER, 0, & comm_handler%get_comm_species(), ierr) call MPI_gather(dens_species_nrm2, 1, MPI_GP, & gath_nrm2, 1, MPI_GP, 0, & comm_handler%get_comm_species(), ierr) call MPI_gather(dens_species_err2nrm, 1, MPI_GP, & gath_err2nrm, 1, MPI_GP, 0, & comm_handler%get_comm_species(), ierr) call MPI_gather(dens_species_errsup, 1, MPI_GP, & gath_errsup, 1, MPI_GP, 0, & comm_handler%get_comm_species(), ierr) ! Write results to mms_zhdanov_job.out file if (is_rank_info_writer) then jobfile = 'mms_zhdanov_job.out' inquire(file=jobfile, exist=fexist) if (fexist) then open(unit=56, file=jobfile, & status='old', & position='append', & action='write', & iostat=ioerr, & iomsg=io_errmsg) else open(unit=56, file=jobfile, & 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(56,603)tau do s = 1, comm_handler%get_nspecies() write(56,*)gath_species(s) write(56,603)gath_nrm2(s) write(56,603)gath_err2nrm(s) write(56,603)gath_errsup(s) 603 FORMAT(ES20.12E3) enddo close(56) endif end subroutine real(GP) function MmsSolDens(equi, species, x, y, z, t) !! MMS solution of density for specific species class(equilibrium_t) :: equi !! Equilibrium integer, intent(in) :: species !! Species number real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: z !! Axial coordinate z (represented by toroidal coordinate) real(GP), intent(in) :: t !! Time select type(equi) type is(circular_t) MmsSolDens = MmsSolDensCircular(equi, species, x, y, z, t) type is(slab_t) MmsSolDens = MmsSolDensSlab(equi, species, x, y, z, t) class default call handle_error('Equilibrium not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end function real(GP) function MmsSourceDens(equi, species, x, y, z, t, params_species) !! MMS solution of density for specific species class(equilibrium_t) :: equi !! Equilibrium integer, intent(in) :: species !! Species number real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: z !! Axial coordinate z (represented by toroidal coordinate) real(GP), intent(in) :: t !! Time type(params_zhdanov_species_t) :: params_species !! Species parameters select type(equi) type is(circular_t) MmsSourceDens = MmsSourceDensCircular(equi, species,& x, y, z, t, params_species) type is(slab_t) MmsSourceDens = MmsSourceDensSlab(equi, species,& x, y, z, t, params_species) class default call handle_error('Equilibrium not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end function end module