mms_template_circular_toroidal_s.f90 Source File

submodule (mms_template_m) mms_template_circular_toroidal_s
    !! MMS Solutions and Sources for template model in circular toroidal geometry
    use precision_grillix_m, only : GP
    use error_handling_grillix_m, only: handle_error
    use status_codes_grillix_m, only : GRILLIX_ERR_OTHER
    use params_template_m, only : dperp_coeff_dens, path_mms
    use static_data_m, only : equi

    use constants_m, only : TWO_PI
    use equilibrium_m, only : equilibrium_t
    use circular_toroidal_equilibrium_m, only : circular_toroidal_t
    implicit none

    ! Default values for circular toroidal MMS parameters
    real(GP) :: bckg_n = 1.1_GP
    !! Constant background of MMS solution
    real(GP) :: amp_n = 1.0_GP
    !! Amplitude of MMS solution
    real(GP) :: kr_n = 1.0_GP
    !! Radial mode number of MMS Solution
    real(GP) :: kp_n = 2.0_GP
    !! Poloidal mode number of MMS Solution
    real(GP) :: php_n = 0.2_GP
    !! Poloidal phase shift of MMS Solution
    real(GP) :: kz_n = 1.0_GP
    !! Axial mode number of MMS Solution
    real(GP) :: phz_n = 1.1_GP
    !! Axial phase shift of MMS Solution
    real(GP) :: omega_n = 0.5_GP
    !! Angular time frequency of MMS Solution
    real(GP) :: pht_n = 0.4_GP
    !! Temporal phase shift of MMS Solution

    namelist / mms_circtor_template / bckg_n, amp_n, kr_n, kp_n, php_n, kz_n, phz_n, omega_n, pht_n

contains
    
    module subroutine init_mms_circular_toroidal_template()
        !! Initializes MMS parameters for circular toroidal geometry
        integer :: io_unit, io_error
        character(len=256) :: io_errmsg

        ! Read circular toroidal MMS parameters
        open(newunit=io_unit, file=path_mms, status='old', action='read', iostat=io_error, iomsg=io_errmsg)
        if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
        read(io_unit, nml=mms_circtor_template, iostat=io_error, iomsg=io_errmsg)
        if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
        close(io_unit)

        ! Set procedure pointers
        mms_sol_dens_template => mms_sol_dens_template_circular_toroidal
        mms_source_dens_template => mms_source_dens_template_circular_toroidal

    end subroutine

    module subroutine display_mms_params_template_circular_toroidal()

        if (is_rank_info_writer) then
            write(get_stdout(), nml=mms_circtor_template)
        endif

    end subroutine
    
    function mms_sol_dens_template_circular_toroidal(x, y, z, t) result(res)
        !! MMS solution
        real(GP), intent(in) :: x
        !! x-coordinate 
        real(GP), intent(in) :: y
        !! y-coordinate
        real(GP), intent(in) :: z
        !! Axial coordinate phi (toroidal coordinate)
        real(GP), intent(in) :: t
        !! Time

        real(GP) :: q, rmin, rmax

        include 'generated/mms_sol_Dens_circular_toroidal_head.inc'

        select type(equi)
            type is(circular_toroidal_t)
                q     = equi%qval(equi%rho(x, y, z))
                rmin  = equi%rhomin
                rmax  = equi%rhomax
            class default
                call handle_error('Equilibrium is not circular', GRILLIX_ERR_OTHER, __LINE__, __FILE__)
        end select

        include 'generated/mms_sol_Dens_circular_toroidal.inc'
                    
    end function
    
    function mms_source_dens_template_circular_toroidal(x, y, z, t) result(res)
        !! Source corresponding to MMS solution
        real(GP), intent(in) :: x
        !! x-coordinate 
        real(GP), intent(in) :: y
        !! y-coordinate
        real(GP), intent(in) :: z
        !! Axial coordinate phi (toroidal coordinate)
        real(GP), intent(in) :: t
        !! Time
        
        real(GP) :: q, rmin, rmax
        
        include 'generated/mms_source_Dens_circular_toroidal_head.inc'

        select type(equi)
            type is(circular_toroidal_t)
                q     = equi%qval(equi%rho(x, y, z))
                rmin  = equi%rhomin
                rmax  = equi%rhomax
            class default
                call handle_error('Equilibrium is not circular', GRILLIX_ERR_OTHER, __LINE__, __FILE__)
        end select

        include 'generated/mms_source_Dens_circular_toroidal.inc'

    end function
  
end submodule