runge_kutta_s.f90 Source File

submodule (runge_kutta_m) runge_kutta_s
    implicit none
   
contains

    pure integer module function get_ncount(self)
        class(runge_kutta_t), intent(in) :: self
        get_ncount = self%ncount
    end function

    pure integer module function get_ndim(self)
        class(runge_kutta_t), intent(in) :: self
        get_ndim = self%ndim
    end function

    pure integer module function get_order(self)
        class(runge_kutta_t), intent(in) :: self
        get_order = self%order
    end function
       
    module subroutine init_runge_kutta(self, ndim, order)
        class(runge_kutta_t), intent(inout) :: self
        integer, intent(in) :: ndim
        integer, intent(in) :: order
    
        self%ncount = 0
        self%ndim   = ndim
        self%order  = order

    end subroutine

    module subroutine advance(self, dtau, tau, y, rhs)
        class(runge_kutta_t), intent(inout) :: self
        real(GP), intent(in) :: dtau
        real(GP), intent(in) :: tau
        real(GP), dimension(self%ndim), intent(inout) :: y
        interface
            subroutine rhs(ndim, tau, y, dy)
                use precision_grillix_m, only : GP
                integer, intent(in) :: ndim
                real(GP), intent(in) :: tau
                real(GP), dimension(ndim), intent(in) :: y
                real(GP), dimension(ndim), intent(out) :: dy
            end subroutine
        end interface 

        integer :: i
        real(GP) :: tau2, tau3, tau4
        real(GP), dimension(self%ndim) :: wrk, k1, k2, k3, k4

        
        select case(self%order)
        case(1)

            call rhs(self%ndim, tau, y, k1)
            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                y(i) = y(i) + dtau * k1(i)
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 

        case(2)

            call rhs(self%ndim, tau, y, k1)

            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                wrk(i) = y(i) + dtau * k1(i)
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 
            tau2 = tau + dtau         
            call rhs(self%ndim, tau2, wrk, k2)

            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                y(i) = y(i) + dtau / 2.0_GP * ( k1(i) + k2(i) )
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 


        case(3)

            call rhs(self%ndim, tau, y, k1)

            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                wrk(i) = y(i) + dtau/2.0_GP * k1(i)
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 
            tau2 = tau + dtau / 2.0_GP
            call rhs(self%ndim, tau2, wrk, k2)

            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                wrk(i) = y(i) + dtau * ( -k1(i) + 2*k2(i) )
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 
            tau3 = tau + dtau
            call rhs(self%ndim, tau3, wrk, k3)

            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                y(i) = y(i) + dtau / 6.0_GP * ( k1(i) + 4*k2(i) + k3(i) )
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 

        case(4)

            call rhs(self%ndim, tau, y, k1)

            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                wrk(i) = y(i) + dtau/2.0_GP * k1(i)
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 
            tau2 = tau + dtau / 2.0_GP
            call rhs(self%ndim, tau2, wrk, k2)

            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                wrk(i) = y(i) + dtau/2.0_GP * k2(i)
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 
            tau3 = tau + dtau/2.0_GP
            call rhs(self%ndim, tau3, wrk, k3)

            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                wrk(i) = y(i) + dtau * k3(i)
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 
            tau4 = tau + dtau
            call rhs(self%ndim, tau4, wrk, k4)
        
            !$OMP PARALLEL PRIVATE(i)
            !$OMP DO SIMD
            do i = 1, self%ndim
                y(i) = y(i) + dtau / 6.0_GP * ( k1(i) + 2*k2(i) + 2*k3(i) + k4(i) )
            enddo
            !$OMP END DO SIMD
            !$OMP END PARALLEL 
        
        case default
            call handle_error('Runge-Kutta order not valid', &
                              GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                              error_info_t('Order: ',[self%order]))
        end select
   
        self%ncount = self%ncount + 1
    
    end subroutine



    module subroutine display_runge_kutta(self)
        class(runge_kutta_t), intent(in) :: self

        if (.not. is_rank_info_writer) then
            return
        endif
        
        write(get_stdout(),*)'Information on runge_kutta ----------------------------'
        write(get_stdout(),101)self%ncount, self%ndim, self%order                             
101     FORMAT( 3X,'ncount      = ',I10        /, &
                3X,'ndim        = ',I10        /, &
                3X,'order       = ',I3            )  
        write(get_stdout(),*)'-------------------------------------------------------'

    end subroutine 


    module subroutine destructor(self)
        type(runge_kutta_t), intent(inout) :: self

        self%ncount = 0
        self%ndim = 0
        self%order = 0
 
    end subroutine destructor

end submodule