mms_zhdanov_circular_s.f90 Source File

submodule(mms_zhdanov_m) mms_zhdanov_circular_s 
    !! MMS Solutions and Sources for Zhdanov model in circular geometry
    implicit none
    
    real(GP), parameter :: ZaDe = 1.0_GP, &
                           ZaTr = 1.0_GP, &
                           ZaHe = 2.0_GP
       
    real(GP), parameter :: AmpDensDe = 1.0_GP, &
                           AmpDensTr = 0.3_GP, &
                           AmpDensHe = 0.1_GP
                           
    real(GP), parameter :: KrDensDe = 1.0_GP, &
                           KrDensTr = 2.0_GP, &
                           KrDensHe = 1.0_GP
                           
    real(GP), parameter :: KpDensDe = 2.0_GP, &
                           KpDensTr = 1.0_GP, &
                           KpDensHe = 3.0_GP
                        
    real(GP), parameter :: PPhDensDe = 0.5_GP, &
                           PPhDensTr = 0.9_GP, &
                           PPhDensHe = 1.6_GP

    real(GP), parameter :: KzDensDe = 2.0_GP, &
                           KzDensTr = 1.0_GP, &
                           KzDensHe = 1.0_GP
                           
    real(GP), parameter :: ZPhDensDe = 0.8_GP, &
                           ZPhDensTr = 1.5_GP, &
                           ZPhDensHe = 1.9_GP
                           
    real(GP), parameter :: OmegaDensDe = 50.0_GP, &
                           OmegaDensTr = 96.0_GP, &
                           OmegaDensHe = 73.0_GP
                           
    real(GP), parameter :: TPhDensDe = 2.0_GP, &
                           TPhDensTr = 1.2_GP, &
                           TPhDensHe = 0.9_GP
                                                  
    !! Offset of MMS solutions (for positively defined quantites, i.e. density and temperatures)
    real(GP), parameter :: OffsetDensDe = AmpDensDe + 1.5E-1_GP, & 
                           OffsetDensTr = AmpDensTr + 1.3E-1_GP, &
                           OffsetDensHe = AmpDensHe + 1.2E-1_GP
contains

    real(GP) module function MmsSolDensCircular(equi, species, x, y, z, t)
        class(equilibrium_t) :: equi
        integer, intent(in) :: species
        real(GP), intent(in) :: x
        real(GP), intent(in) :: y
        real(GP), intent(in) :: z
        real(GP), intent(in) :: t
        
        real(GP) ::  rmin, rmax, q, r, p, twopi
                          
        twopi = TWO_PI     
        ! Convert Cartesian to polar coordinates and safety factor
        call cart_to_circ(equi, x, y, rho = r, theta = p, rhomin = rmin,&
            rhomax = rmax, qval = q)  
        
        select case(species)
        case(0)
            ! Electrons
            MmsSolDensCircular = MmsSolDensCircularEl()
        case(1)
            ! Deuterium
            MmsSolDensCircular = MmsSolDensCircularDe()
        case(2)
            ! Tritium
            MmsSolDensCircular = MmsSolDensCircularTr()
        case(3)
            ! Helium
            MmsSolDensCircular = MmsSolDensCircularHe()
        case default
            call handle_error('Species not valid', &
                              GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                              error_info_t('Species: ',[species]))
        end select

    contains
        
        real(GP) function MmsSolDensCircularEl() result(res)
            include 'MmsSol/MmsSolDensEl_header.txt'
            include 'MmsSol/MmsSolDensEl.txt'
        end
        
        real(GP) function MmsSolDensCircularDe() result(res)
            include 'MmsSol/MmsSolDensDe_header.txt'
            include 'MmsSol/MmsSolDensDe.txt'
        end
        
        real(GP) function MmsSolDensCircularTr() result(res)
            include 'MmsSol/MmsSolDensTr_header.txt'
            include 'MmsSol/MmsSolDensTr.txt'
        end
        
        real(GP) function MmsSolDensCircularHe() result(res)
            include 'MmsSol/MmsSolDensHe_header.txt'
            include 'MmsSol/MmsSolDensHe.txt'
        end

    end function  
    
    real(GP) module function MmsSourceDensCircular(equi, species,&
        x, y, z, t, params_species)
        class(equilibrium_t) :: equi
        integer, intent(in) :: species
        real(GP), intent(in) :: x
        real(GP), intent(in) :: y
        real(GP), intent(in) :: z
        real(GP), intent(in) :: t
        type(params_zhdanov_species_t) :: params_species
        
        real(GP) :: rmin, rmax, q, r, p, twopi, SourceParam
                          
        twopi = TWO_PI     
        ! Convert Cartesian to polar coordinates and safety factor
        call cart_to_circ(equi, x, y, rho = r, theta = p, rhomin = rmin,&
            rhomax = rmax, qval = q)  
        
        SourceParam = params_species%source_param
        
        select case(species)
        case(0)
            ! Electrons
            MmsSourceDensCircular = 0.0_GP
        case(1)
            ! Deuterium
            MmsSourceDensCircular = MmsSourceDensCircularDe()
        case(2)
            ! Tritium
            MmsSourceDensCircular = MmsSourceDensCircularTr()
        case(3)
            ! Helium
            MmsSourceDensCircular = MmsSourceDensCircularHe()
        case default
            call handle_error('Species not valid', &
                              GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                              error_info_t('Species: ',[species]))
        end select      
        
    contains
                
        real(GP) function MmsSourceDensCircularDe() result(res)
            include 'MmsSource/MmsSourceDensDe_header.txt'
            include 'MmsSource/MmsSourceDensDe.txt'
        end
        
        real(GP) function MmsSourceDensCircularTr() result(res)
            include 'MmsSource/MmsSourceDensTr_header.txt'
            include 'MmsSource/MmsSourceDensTr.txt'
        end
        
        real(GP) function MmsSourceDensCircularHe() result(res)
            include 'MmsSource/MmsSourceDensHe_header.txt'
            include 'MmsSource/MmsSourceDensHe.txt'
        end

    end function  
       
    subroutine cart_to_circ(equi, x, y, rho, theta, rhomin, rhomax, qval) 
        !! Converts Cartesian (x,y) coordinates to (rho, theta) coordinates and returns equilibrium specific quantities
        class(equilibrium_t), intent(inout) :: equi
        !! Equilibrium
        real(GP), intent(in) :: x
        !! x-coordinate 
        real(GP), intent(in) :: y
        !! y-coordinate
        real(GP), intent(out) :: rho
        !! radial coordinate 
        real(GP), intent(out) :: theta
        !! poloidal angle
        real(GP), intent(out) :: rhomin
        !! Inner limiting flux surface
        real(GP), intent(out) :: rhomax
        !! Outer limiting flux surface
        real(GP), intent(out) :: qval
        !! safety factor
        
        real(GP) :: phi

        select type(equi)
            type is(circular_t)
                phi = 0.0_GP ! Axisymmetry
                rho     = equi%rho(x, y, phi)
                theta   = equi%theta(x, y)
                rhomin  = equi%rhomin
                rhomax  = equi%rhomax
                qval    = equi%qval(rho)
            class default
                call handle_error('Equilibrium is not circular', &
                                  GRILLIX_ERR_OTHER, __LINE__, __FILE__)
        end select

    end subroutine
    

end submodule