variable_set_s.f90 Source File

submodule (variable_m) variable_set_s
    !! Routines related with setting variable values
    implicit none
contains
      
    module subroutine set_blob_toroidal(self, mesh, x0, y0, phi0, wx, wy, wphi, amp)
        class(variable_t), intent(inout) :: self  
        type(mesh_cart_t), intent(in) :: mesh
        real(GP), intent(in) :: x0
        real(GP), intent(in) :: y0
        real(GP), intent(in) :: phi0
        real(GP), intent(in) :: wx
        real(GP), intent(in) :: wy
        real(GP), intent(in) :: wphi
        real(GP), intent(in) :: amp
        
        real(GP) :: phi, gauss_tor, x, y
        integer :: l
        
        phi = mesh%get_phi()
    
        ! Gaussian argument in toroidal direction
        if (wphi <= GP_EPS) then
            ! Delta peak 
            if (abs(phi - phi0) <= GP_EPS) then
                gauss_tor = amp
            else
                gauss_tor = 0.0_GP
            endif           
        else
            gauss_tor = amp * gaussian(phi0, wphi, phi)
        endif

        !$OMP PARALLEL DO PRIVATE(l, x, y)
        do l = 1,mesh%get_n_points()
            x = mesh%get_x(l)
            y = mesh%get_y(l)
            self%vals(l) = gauss_tor * gaussian(x0, wx, x) * gaussian(y0, wy, y)
        enddo
        !$OMP END PARALLEL DO
        
    end subroutine
    
    module subroutine set_blob_aligned(self, equi, mesh, x0, y0, phi0, wx, wy, wpar, amp)
        class(variable_t), intent(inout) :: self  
        class(equilibrium_t), intent(inout) :: equi
        type(mesh_cart_t), intent(in) :: mesh
        real(GP), intent(in) :: x0
        real(GP), intent(in) :: y0
        real(GP), intent(in) :: phi0
        real(GP), intent(in) :: wx
        real(GP), intent(in) :: wy
        real(GP), intent(in) :: wpar
        real(GP), intent(in) :: amp
        
        real(GP) :: phi, xs, ys, dphi_trace, xf, yf, dpar
        integer :: l, nturns, iturns
        
        phi = mesh%get_phi()
        
        ! Number of toroidal turns to be traced, captures roughly 8*parallel width of blob
        nturns = max(1, ceiling(8 * wpar/TWO_PI))

        !$omp parallel private(l, xs, ys, dphi_trace, xf, yf, dpar, iturns)
        !$omp do
        do l = 1, mesh%get_n_points()
            xs = mesh%get_x(l)
            ys = mesh%get_y(l)

            self%vals(l) = 0.0_GP
            do iturns = -nturns, nturns
    
                dphi_trace = phi0 - phi + iturns * TWO_PI
                                         
                call trace(x_init=xs, y_init=ys, phi_init=phi, &
                           dphi=dphi_trace, equi=equi, &
                           x_end=xf, y_end=yf, arclen=dpar)

                self%vals(l) = self%vals(l) +  amp * &
                              gaussian(x0, wx, xf) * &
                              gaussian(y0, wy, yf) * &
                              gaussian(0.0_GP, wpar, dpar)
            enddo

        enddo
        !$omp end do
        !$omp end parallel
        
    end subroutine

end submodule