module mms_diffusion_m !! Method of Manufactured Solutions for diffusion model use MPI use comm_handler_m, only : comm_handler_t use parallelisation_setup_m, only : is_rank_info_writer use precision_grillix_m, only : GP, GP_EPS, MPI_GP 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 parallel_map_m, only : parallel_map_t use commons_misc_m, only : error_ananum use mms_diffusion_circular_m, only : mms_sol_diffusion_circular, mms_source_diffusion_circular use mms_diffusion_slab_m, only : mms_sol_diffusion_slab, mms_source_diffusion_slab use mms_diffusion_dommaschk_m, only : mms_sol_diffusion_dommaschk, mms_source_diffusion_dommaschk 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 dommaschk_equilibrium_m, only : dommaschk_t use mesh_cart_m, only : mesh_cart_t implicit none public :: mms_sol_diffusion public :: mms_source_diffusion public :: mms_diagnostics_diffusion contains subroutine mms_diagnostics_diffusion(comm_handler, equi, mesh, map, tau, num_sol) !! Prints information on numerical errors of numerical solution type(comm_handler_t), intent(in) :: comm_handler class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Mesh within poloidal plane real(GP), intent(in) :: tau !! Time type(variable_t), intent(in) :: num_sol !! Numerical solution integer :: l real(GP) :: x, y, z, nrm2, nrmsup, err2, errsup, err2nrm real(GP), dimension(mesh%get_n_points()) :: mms_sol logical :: fexist integer :: ioerr character(len=256) :: io_errmsg z = map%get_phi_cano() !$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(l) = mms_sol_diffusion(equi, x, y, z, tau) enddo !$OMP END DO !$OMP END PARALLEL call error_ananum(comm_handler, mesh%get_n_points(), mms_sol, num_sol%vals, nrm2, nrmsup, err2, errsup) err2nrm = err2/(nrm2+GP_EPS) if (is_rank_info_writer) then write(get_stdout(),*)'MMS diagnostics ---------------------------------------' write(get_stdout(),204)tau, nrm2, nrmsup, err2, err2nrm, errsup 204 FORMAT(3X,'tau = ',ES14.6E3 /, & 3X,'nrm2 = ',ES14.6E3 /, & 3X,'nrmsup = ',ES14.6E3 /, & 3X,'err2 = ',ES14.6E3 /, & 3X,'err2nrm = ',ES14.6E3 /, & 3X,'errsup = ',ES14.6E3 ) write(get_stdout(),*)'-------------------------------------------------------' endif ! Write out result to file for integration tests ------------------------------------------ if (is_rank_info_writer) then inquire(file='mms_diffusion_job.out', exist = fexist) if (fexist) then open(unit = 56, file = 'mms_diffusion_job.out', & status ='old', & position = 'append', & action = 'write', & iostat = ioerr, & iomsg=io_errmsg) else open(unit = 56, & file = 'mms_diffusion_job.out', & 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 write(56,603)nrm2 write(56,603)err2 write(56,603)err2nrm write(56,603)errsup 603 FORMAT(ES20.12E3) close(56) 301 FORMAT(ES20.12E3) endif end subroutine real(GP) function mms_sol_diffusion(equi, x, y, z, t) !! MMS solution class(equilibrium_t) :: equi !! Equilibrium 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) mms_sol_diffusion = mms_sol_diffusion_circular(equi, x, y, z, t) type is(slab_t) mms_sol_diffusion = mms_sol_diffusion_slab(equi, x, y, z, t) type is(dommaschk_t) mms_sol_diffusion = mms_sol_diffusion_dommaschk(equi, x, y, z, t) class default call handle_error('Equilibrium not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end function real(GP) function mms_source_diffusion(equi, x, y, z, t, charfun, penval) !! Source corresponding to MMS solution class(equilibrium_t), intent(inout) :: equi !! Equilibrium 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 real(GP), intent(in) :: charfun !! Value of characteristic function (penalisation) real(GP), intent(in) :: penval !! Value that solution shall be penalised to select type(equi) type is(circular_t) mms_source_diffusion = mms_source_diffusion_circular(equi, x, y, z, t, charfun, penval) type is(slab_t) mms_source_diffusion = mms_source_diffusion_slab(equi, x, y, z, t) type is(dommaschk_t) mms_source_diffusion = mms_source_diffusion_dommaschk(equi, x, y, z, t) class default call handle_error('Equilibrium not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end function end module