module mms_neutrals_m !! Implementation of Method of Manufactured Solutions for the Neutrals model use MPI use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP, GP_EPS use error_handling_grillix_m, only: handle_error use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST, GRILLIX_ERR_OTHER use variable_m, only : variable_t use screen_io_m, only : get_stdout use parallelisation_setup_m, only : is_rank_info_writer use circular_equilibrium_m, only : circular_t use slab_equilibrium_m, only : slab_t use dommaschk_equilibrium_m, only : dommaschk_t use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use parallel_map_m, only : parallel_map_t use commons_misc_m, only : error_ananum use rate_coeff_neutrals_m, only: rate_coeff_neutrals_t use mms_neutrals_circular_m use mms_neutrals_slab_m use mms_neutrals_dommaschk_m implicit none public :: mms_diagnostics_neutrals public :: mms_sol_neutrals_dens public :: mms_sol_neutrals_parmom public :: mms_sol_neutrals_pressure public :: mms_source_neutrals_dens public :: mms_source_neutrals_parmom public :: mms_source_neutrals_pressure contains subroutine mms_diagnostics_neutrals(comm_handler, equi, & mesh_cano, mesh_stag, map, tau, & neutrals_dens, neutrals_parmom, neutrals_pressure) !! Prints information on numerical errors of numerical solutions type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(in) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) type(parallel_map_t), intent(in) :: map !! Mesh within poloidal plane real(GP), intent(in) :: tau !! Time type(variable_t), intent(in) :: neutrals_dens !! Numerical solution for neutrals density type(variable_t), intent(in) :: neutrals_parmom !! Numerical solution for neutrals parallel momentum type(variable_t), intent(in) :: neutrals_pressure !! Numerical solution for neutrals temp integer :: l real(GP) :: x, y, phi_cano, phi_stag real(GP), dimension(mesh_cano%get_n_points()) :: neutrals_dens_mms real(GP), dimension(mesh_stag%get_n_points()) :: neutrals_parmom_mms real(GP), dimension(mesh_cano%get_n_points()) :: neutrals_pressure_mms real(GP) :: neutrals_dens_nrm2, neutrals_dens_nrmsup, & neutrals_dens_err2, neutrals_dens_errsup real(GP) :: neutrals_parmom_nrm2, neutrals_parmom_nrmsup, & neutrals_parmom_err2, neutrals_parmom_errsup real(GP) :: neutrals_pressure_nrm2, neutrals_pressure_nrmsup, & neutrals_pressure_err2, neutrals_pressure_errsup character(len=22) :: fname logical :: fexist integer :: ioerr character(len=256) :: io_errmsg phi_cano = mesh_cano%get_phi() phi_stag = mesh_stag%get_phi() ! Get MMS solutions ----------------------------------------------------------------------- !$omp parallel default(none) private(l, x, y) & !$omp shared(equi, tau, mesh_cano, phi_cano, mesh_stag, phi_stag, & !$omp neutrals_dens_mms, neutrals_parmom_mms, neutrals_pressure_mms) !$omp do do l = 1, mesh_cano%get_n_points() x = mesh_cano%get_x(l) y = mesh_cano%get_y(l) neutrals_dens_mms(l) = mms_sol_neutrals_dens(equi, x, y, phi_cano, tau) neutrals_pressure_mms(l) = mms_sol_neutrals_pressure(equi, x, y, phi_cano, tau) enddo !$omp end do !$omp do do l = 1, mesh_stag%get_n_points() x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) neutrals_parmom_mms(l) = mms_sol_neutrals_parmom(equi, x, y, phi_stag, tau) enddo !$omp end do !$omp end parallel ! Compute errors between MMS and numeric solution ----------------------------------------- call error_ananum(comm_handler, mesh_cano%get_n_points(), neutrals_dens_mms, neutrals_dens%vals, & neutrals_dens_nrm2, neutrals_dens_nrmsup, neutrals_dens_err2, neutrals_dens_errsup) call error_ananum(comm_handler, mesh_stag%get_n_points(), neutrals_parmom_mms, neutrals_parmom%vals, & neutrals_parmom_nrm2, neutrals_parmom_nrmsup, neutrals_parmom_err2, neutrals_parmom_errsup) call error_ananum(comm_handler, mesh_cano%get_n_points(), neutrals_pressure_mms, neutrals_pressure%vals, & neutrals_pressure_nrm2, neutrals_pressure_nrmsup, neutrals_pressure_err2, neutrals_pressure_errsup) ! Write out results on display ------------------------------------------------------------ if (is_rank_info_writer) then write(get_stdout(),*)'MMS Neutrals diagnostics ------------------------------' write(get_stdout(),202)tau 202 FORMAT(3X,'tau = ',ES14.6E3) write(get_stdout(),203)'Field','norm(L2)','norm(Linf)','err(L2)','errnrm(L2)', 'err(Linf)', 'errnrm(Linf)' 203 FORMAT(A16,6(A19)) write(get_stdout(),204)'neutrals_dens:', neutrals_dens_nrm2, neutrals_dens_nrmsup, & neutrals_dens_err2, neutrals_dens_err2/(neutrals_dens_nrm2+GP_EPS), & neutrals_dens_errsup, neutrals_dens_errsup/(neutrals_dens_nrmsup+GP_EPS) write(get_stdout(),204)'neutrals_parmom:', neutrals_parmom_nrm2, neutrals_parmom_nrmsup, & neutrals_parmom_err2, neutrals_parmom_err2/(neutrals_parmom_nrm2+GP_EPS), & neutrals_parmom_errsup, neutrals_parmom_errsup/(neutrals_parmom_nrmsup+GP_EPS) write(get_stdout(),204)'neutrals_pressure:', neutrals_pressure_nrm2, neutrals_pressure_nrmsup, & neutrals_pressure_err2, neutrals_pressure_err2/(neutrals_pressure_nrm2+GP_EPS), & neutrals_pressure_errsup, neutrals_pressure_errsup/(neutrals_pressure_nrmsup+GP_EPS) 204 FORMAT(3X,A16,6(5X,ES14.6E3)) write(get_stdout(),*)'-------------------------------------------------------' endif ! Write out result to file for integration tests ------------------------------------------ if (is_rank_info_writer) then fname = 'mms_neutrals_job.out' inquire(file=fname, exist = fexist) if (fexist) then open(unit=56, file=fname, status='old', position='append', & action='write', iostat=ioerr, iomsg=io_errmsg) else open(unit=56, 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(56,603)tau write(56,603)neutrals_dens_nrm2 write(56,603)neutrals_dens_nrmsup write(56,603)neutrals_dens_err2 write(56,603)neutrals_dens_errsup write(56,603)neutrals_parmom_nrm2 write(56,603)neutrals_parmom_nrmsup write(56,603)neutrals_parmom_err2 write(56,603)neutrals_parmom_errsup write(56,603)neutrals_pressure_nrm2 write(56,603)neutrals_pressure_nrmsup write(56,603)neutrals_pressure_err2 write(56,603)neutrals_pressure_errsup 603 FORMAT(ES20.12E3) close(56) 301 FORMAT(ES20.12E3) endif end subroutine real(GP) function mms_sol_neutrals_dens(equi, x, y, z, t) !! MMS solution for neutrals density class(equilibrium_t), intent(in) :: 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_neutrals_dens = mms_sol_circular_neutrals_dens(equi, x, y, z, t) type is(slab_t) mms_sol_neutrals_dens = mms_sol_slab_neutrals_dens(equi, x, y, z, t) type is(dommaschk_t) mms_sol_neutrals_dens = mms_sol_dommaschk_neutrals_dens(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_sol_neutrals_parmom(equi, x, y, z, t) !! MMS solution for neutrals parallel momentum class(equilibrium_t), intent(in) :: 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_neutrals_parmom = mms_sol_circular_neutrals_parmom(equi, x, y, z, t) type is(slab_t) mms_sol_neutrals_parmom = mms_sol_slab_neutrals_parmom(equi, x, y, z, t) type is(dommaschk_t) mms_sol_neutrals_parmom = mms_sol_dommaschk_neutrals_parmom(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_sol_neutrals_pressure(equi, x, y, z, t) !! MMS solution for neutrals temperature class(equilibrium_t), intent(in) :: 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_neutrals_pressure = mms_sol_circular_neutrals_pressure(equi, x, y, z, t) type is(slab_t) mms_sol_neutrals_pressure = mms_sol_slab_neutrals_pressure(equi, x, y, z, t) type is(dommaschk_t) mms_sol_neutrals_pressure = mms_sol_dommaschk_neutrals_pressure(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_neutrals_dens(equi, x, y, z, t, chi, epsinv, & neutrals_dens_pen, k_iz_o, k_rec_o) !! MMS source for Neutrals density equation class(equilibrium_t), intent(in) :: 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: epsinv !! Inverse of penalisation epsilon real(GP), intent(in) :: neutrals_dens_pen !! Target value to be penalised towards type(rate_coeff_neutrals_t), intent(in) :: k_iz_o !! Ionization rate coefficient object type(rate_coeff_neutrals_t), intent(in) :: k_rec_o !! Recombination rate coefficient object select type(equi) type is(circular_t) mms_source_neutrals_dens = & mms_source_circular_neutrals_dens(equi, x, y, z, t, & chi, epsinv, neutrals_dens_pen, k_iz_o, k_rec_o) type is(slab_t) mms_source_neutrals_dens = & mms_source_slab_neutrals_dens(equi, x, y, z, t, & k_iz_o, k_rec_o) type is(dommaschk_t) mms_source_neutrals_dens = & mms_source_dommaschk_neutrals_dens(equi, x, y, z, t, & k_iz_o, k_rec_o) class default call handle_error('Equilibrium not valid', GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end function real(GP) function mms_source_neutrals_parmom(equi, x, y, z, t, chi, epsinv, & neutrals_parmom_pen, k_iz_o, k_rec_o) !! MMS source for Neutrals parallel momentum equation class(equilibrium_t), intent(in) :: 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: epsinv !! Inverse of penalisation epsilon real(GP), intent(in) :: neutrals_parmom_pen !! Target value to be penalised towards type(rate_coeff_neutrals_t), intent(in) :: k_iz_o !! Ionization rate coefficient object type(rate_coeff_neutrals_t), intent(in) :: k_rec_o !! Recombination rate coefficient object select type(equi) type is(circular_t) mms_source_neutrals_parmom = & mms_source_circular_neutrals_parmom(equi, x, y, z, t, & chi, epsinv, neutrals_parmom_pen, k_iz_o, k_rec_o) type is(slab_t) mms_source_neutrals_parmom = & mms_source_slab_neutrals_parmom(equi, x, y, z, t, & k_iz_o, k_rec_o) type is(dommaschk_t) mms_source_neutrals_parmom = & mms_source_dommaschk_neutrals_parmom(equi, x, y, z, t, & k_iz_o, k_rec_o) class default call handle_error('Equilibrium not valid', GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end function real(GP) function mms_source_neutrals_pressure(equi, x, y, z, t, chi, epsinv, & neutrals_pressure_pen, k_iz_o, k_rec_o) !! MMS source for Neutrals temperature equation class(equilibrium_t), intent(in) :: 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: epsinv !! Inverse of penalisation epsilon real(GP), intent(in) :: neutrals_pressure_pen !! Target value to be penalised towards type(rate_coeff_neutrals_t), intent(in) :: k_iz_o !! Ionization rate coefficient object type(rate_coeff_neutrals_t), intent(in) :: k_rec_o !! Recombination rate coefficient object select type(equi) type is(circular_t) mms_source_neutrals_pressure = & mms_source_circular_neutrals_pressure(equi, x, y, z, t, & chi, epsinv, neutrals_pressure_pen, k_iz_o, k_rec_o) type is(slab_t) mms_source_neutrals_pressure = & mms_source_slab_neutrals_pressure(equi, x, y, z, t, & k_iz_o, k_rec_o) type is(dommaschk_t) mms_source_neutrals_pressure = & mms_source_dommaschk_neutrals_pressure(equi, x, y, z, t, & k_iz_o, k_rec_o) class default call handle_error('Equilibrium not valid', GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end function end module