-
Notifications
You must be signed in to change notification settings - Fork 0
/
soln_RK4.F90
56 lines (43 loc) · 1.36 KB
/
soln_RK4.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
subroutine soln_RK4(dt, j, Vj, Flux)
!performs the jth step of the RK4 algorithm
!Total Flux = (k1 + 2k2 + 2k3 + k4)/6
#include "definition.h"
use grid_data
use primconsflux
implicit none
real, intent(IN) :: dt
integer, intent(IN) :: j
real, dimension(NUMB_VAR,gr_imax), intent(INOUT) :: Vj
real, dimension(NSYS_VAR,gr_imax), intent(INOUT) :: Flux
real :: A, F, dtx
integer :: i
real, dimension(NUMB_VAR,gr_imax) :: Uj
if (j == 1 .OR. j == 2) then
A = 0.5
else
A = 1.
end if
if (j == 1 .OR. j == 4) then
F = 1./6.
else
F = 1./3.
end if
dtx = A*dt/gr_dx
call soln_reconstruct(dt, Vj)
call soln_getFlux()
if (j .NE. 4) then
!update cons variables to jth step only if not 4th step
do i = gr_ibeg, gr_iend
Uj(DENS_VAR:ENER_VAR,i) = gr_U(DENS_VAR:ENER_VAR,i) - &
dtx*(gr_flux(DENS_VAR:ENER_VAR,i+1) - gr_flux(DENS_VAR:ENER_VAR,i))
end do
!! get updated primitive vars from the updated conservative vars
do i = gr_ibeg, gr_iend
! Eos is automatically callled inside cons2prim
call cons2prim(Uj(DENS_VAR:ENER_VAR,i),Vj(DENS_VAR:GAME_VAR,i))
end do
end if
!finally add to total flux
Flux(DENS_VAR:ENER_VAR,gr_i0:gr_imax) = Flux(DENS_VAR:ENER_VAR,gr_i0:gr_imax) + &
F*gr_flux(DENS_VAR:ENER_VAR,gr_i0:gr_imax)
end subroutine soln_RK4