MAKE_BUTCHERARRAY.f90

Go to the documentation of this file.

!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.
 
 
!  where p is the order of the scheme and s is the number of stages
!  Schemes Implemented: RK(2,2), RK(3,3), RK(4,4)  
 
 
 
     subroutine make_butcherarray(order, stages, A_low, B_low, c)
          
     implicit none
     
     integer*4:: order, stages
     real*8, dimension(stages,stages) :: a
     real*8, dimension(stages) :: b, c, a_low, b_low
 
     a = 0; b = 0; c = 0; a_low = 0; b_low = 0;
     
     select case(order)   
 
        case(2)               
           ! RK(2,2)
           select case(stages)
              case(2)
                 c(2) = 1
                 b(1) = 0.5; b(2)=0.5
                 a(2,1) = 1
                 
                 call make_low_storage_coefficients(stages,a,b,c,a_low,b_low)
                 
              case default   
                  write(*,'(A)') 'RK scheme not implemented'
           end select               
        
        case(3)
           ! RK(3,3)
           select case(stages)
              case(3)
                 c(2) = 0.5; c(3) = 1
                 b(1) = 1.d0/6; b(2) = 2.d0/3; b(3) = 1.d0/6
                 a(2,1) = 0.5; 
                 a(3,1) = -1; a(3,2) = 2
                 
                 call make_low_storage_coefficients(stages,a,b,c,a_low,b_low)
                                  
 
              case default   
                  write(*,'(A)') 'RK scheme not implemented'
           end select                       
        
        case(4)
 
           select case(stages)
              ! RK(4,4)
              case(4)
                 c(2) = 0.5; c(3) = 0.5; c(4) = 1
                 b(1) = 1.d0/6; b(2) = 1.d0/3; b(3) = 1.d0/3; b(4) = 1.d0/6
                 a(2,1) = 0.5; 
                 a(3,2) = 0.5;
                 a(4,3) = 1;
                 
                 call make_low_storage_coefficients(stages,a,b,c,a_low,b_low)
               
                 
              case default   
                  write(*,'(A)') 'RK scheme not implemented'
                  
           end select                
       
       
      
       
       
        case default
           write(*,'(A)') 'RK scheme not implemented'
           
      end select    
           
     
     end subroutine make_butcherarray     
     
!**************************************************************************************************
 
 
 
 
     
     subroutine make_low_storage_coefficients(stages,A,b,c,A_low,B_low)
     
     implicit none
     
     integer*4:: stages
     real*8, dimension(stages,stages) :: a
     real*8, dimension(stages) :: b, c, a_low, b_low
     integer*4 :: i
     
 
     ! Compute Low storage coefficients according to the
     ! paper Toulorge and Desmet: "Optimal Runge-Kutta schemes for discontinuous
     ! Galerkin space discretizations applied to wave propagation problems"
     
     a_low(1) = 0;
     
     do i = 1, stages -1
        b_low(i) = a(i+1,i)
     enddo
     b_low(stages) = b(stages)
    
     do i = 2, stages 
        if (b(i) .ne. 0.d0) then
           a_low(i) = (b(i-1) - b_low(i-1))/b(i);
        else
           a_low(i) = (a(i+1,i-1) - c(i))/b_low(i);
        endif
     enddo   
        
     
     end subroutine make_low_storage_coefficients