module mms_template_m !! Method of Manufactured Solutions for template model use MPI use precision_grillix_m, only : GP, GP_EPS, MPI_GP use comm_handler_m, only : comm_handler_t 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 variable_m, only: variable_t use commons_misc_m, only : error_ananum use static_data_m, only : equi, mesh_cano, masks_cano, mesh_stag, masks_stag, map use state_template_m, only : np_max, dens, parmom, pion use params_template_m, only : path_mmsjobout use screen_io_m, only : get_stdout use equilibrium_m, only : equilibrium_t use circular_toroidal_equilibrium_m, only : circular_toroidal_t use slab_equilibrium_m, only : slab_t use mesh_cart_m, only : mesh_cart_t implicit none private public :: init_mms_template public :: display_mms_params_template public :: mms_diagnostics_template public :: mms_sol_dens_template public :: mms_source_dens_template public :: mms_sol_parmom_template public :: mms_source_parmom_template public :: mms_sol_pion_template public :: mms_source_pion_template ! Interface for parameter initialization and display interface module subroutine init_mms_slab_template() !! Initializes MMS parameters for slab geometry end subroutine module subroutine init_mms_circular_toroidal_template() !! Initializes MMS parameters for circular toroidal geometry end subroutine module subroutine display_mms_params_template_slab() !! Displays MMS parameters for slab geometry end subroutine module subroutine display_mms_params_template_circular_toroidal() !! Displays MMS parameters for circular toroidal geometry end subroutine end interface ! Abstract interface for MMS pointers abstract interface real(GP) function mms_func_interface(x, y, phi, t) !! MMS functions interface import GP real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: phi !! Axial coordinate phi (toroidal coordinate) real(GP), intent(in) :: t !! Time end function end interface ! Declare a procedure pointer procedure(mms_func_interface), pointer :: mms_sol_dens_template => null() !! MMS solution function pointer for density procedure(mms_func_interface), pointer :: mms_source_dens_template => null() !! MMS source function pointer for density procedure(mms_func_interface), pointer :: mms_sol_parmom_template => null() !! MMS solution function pointer for parallel momentum procedure(mms_func_interface), pointer :: mms_source_parmom_template => null() !! MMS source function pointer for parallel momentum procedure(mms_func_interface), pointer :: mms_sol_pion_template => null() !! MMS solution function pointer for ion pressure procedure(mms_func_interface), pointer :: mms_source_pion_template => null() !! MMS source function pointer for ion pressure contains subroutine init_mms_template() !! Initializes MMS parameters for template model select type(equi) type is(circular_toroidal_t) call init_mms_circular_toroidal_template() type is(slab_t) call init_mms_slab_template() class default call handle_error('Equilibrium not valid', GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end subroutine subroutine display_mms_params_template() !! Displays MMS parameters for template model select type(equi) type is(circular_toroidal_t) call display_mms_params_template_circular_toroidal() type is(slab_t) call display_mms_params_template_slab() class default call handle_error('Equilibrium not valid', GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end subroutine subroutine mms_diagnostics_template(comm_handler, tau) !! Prints information on numerical errors of numerical solution type(comm_handler_t), intent(in) :: comm_handler !! Communication handler real(GP), intent(in) :: tau !! Time integer :: l real(GP) :: x, y, phi_cano, phi_stag, & nrm2_dens, nrmsup_dens, err2_dens, errsup_dens, err2nrm_dens, errsupnrm_dens, & nrm2_parmom, nrmsup_parmom, err2_parmom, errsup_parmom, err2nrm_parmom, errsupnrm_parmom, & nrm2_pion, nrmsup_pion, err2_pion, errsup_pion, err2nrm_pion, errsupnrm_pion real(GP), dimension(np_max) :: mms_sol_dens, mms_sol_parmom, mms_sol_pion logical :: fexist integer :: ioerr character(len=256) :: io_errmsg mms_sol_dens = 0.0_GP mms_sol_parmom = 0.0_GP mms_sol_pion = 0.0_GP phi_cano = mesh_cano%get_phi() phi_stag = mesh_stag%get_phi() !$omp parallel default(none) & !$omp private(l, x, y) & !$omp shared(np_max, mesh_cano, masks_cano, phi_cano, mesh_stag, masks_stag, phi_stag, & !$omp tau, mms_sol_dens, mms_sol_dens_template, mms_sol_parmom, mms_sol_parmom_template, mms_sol_pion, mms_sol_pion_template) !$omp do do l = 1, np_max if (masks_cano%is_inner(l) .or. masks_cano%is_ghost(l) .or. masks_cano%is_boundary(l)) then x = mesh_cano%get_x(l) y = mesh_cano%get_y(l) mms_sol_dens(l) = mms_sol_dens_template(x, y, phi_cano, tau) mms_sol_pion(l) = mms_sol_pion_template(x, y, phi_cano, tau) endif if (masks_stag%is_inner(l) .or. masks_stag%is_ghost(l) .or. masks_stag%is_boundary(l)) then x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) mms_sol_parmom(l) = mms_sol_parmom_template(x, y, phi_stag, tau) endif enddo !$omp end do !$omp end parallel call error_ananum(comm_handler, np_max, mms_sol_dens, dens, & nrm2_dens, nrmsup_dens, err2_dens, errsup_dens) err2nrm_dens = err2_dens / (nrm2_dens + GP_EPS) errsupnrm_dens = errsup_dens / (nrmsup_dens + GP_EPS) call error_ananum(comm_handler, np_max, mms_sol_parmom, parmom, & nrm2_parmom, nrmsup_parmom, err2_parmom, errsup_parmom) err2nrm_parmom = err2_parmom / (nrm2_parmom + GP_EPS) errsupnrm_parmom = errsup_parmom / (nrmsup_parmom + GP_EPS) call error_ananum(comm_handler, np_max, mms_sol_pion, pion, & nrm2_pion, nrmsup_pion, err2_pion, errsup_pion) err2nrm_pion = err2_pion / (nrm2_pion + GP_EPS) errsupnrm_pion = errsup_pion / (nrmsup_pion + GP_EPS) if (is_rank_info_writer) then write(get_stdout(),*)'MMS diagnostics ' // repeat('-',64) write(get_stdout(),203)tau 203 FORMAT(3X,'at tau = ',ES14.6E3) write(get_stdout(),204)'nrm2', 'nrmsup', 'err2nrm', 'errsupnrm' 204 FORMAT(17X,A4,11X,A6,9X,A7,8X,A9) write(get_stdout(),205)'dens', nrm2_dens, nrmsup_dens, err2nrm_dens, errsupnrm_dens write(get_stdout(),205)'parmom', nrm2_parmom, nrmsup_parmom, err2nrm_parmom, errsupnrm_parmom write(get_stdout(),205)'pion', nrm2_pion, nrmsup_pion, err2nrm_pion, errsupnrm_pion 205 FORMAT(3X, A10, 3X,4(ES14.6E3,X)) write(get_stdout(),*) repeat('-',80) // new_line('a') endif ! Write out result to file for integration tests ------------------------------------------ if (is_rank_info_writer) then inquire(file=path_mmsjobout, exist=fexist) open(unit = 56, file = path_mmsjobout, & status ='unknown', & position = 'append', & action = 'write', & iostat = ioerr, & iomsg=io_errmsg) if (ioerr /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif write(56, 603) tau, & nrm2_dens, nrmsup_dens, err2nrm_dens, errsupnrm_dens, & nrm2_parmom, nrmsup_parmom, err2nrm_parmom, errsupnrm_parmom, & nrm2_pion, nrmsup_pion, err2nrm_pion, errsupnrm_pion 603 FORMAT(13(ES20.12E3, :, ', ')) close(56) endif end subroutine end module