-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrk2.f90
150 lines (127 loc) · 5.84 KB
/
rk2.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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
module rk2_module
implicit none
contains
subroutine RK2_SSP(dens, velx, vely, pres, ener, dt)
#include "header.h"
use riemann_module, only: hllc, hlle
use recon_module, only: recon_getcellfaces
use eos_module, only: eos_getp
use sim_data, only: grav, usegrav, ilo, ihi, jlo, jhi, &
gamma, dx, dy, yTpts, xTpts, smallf, &
flux_solver
use applyBC_module, only: applyBC_all
use misc_module, only: to_upper
implicit none
real(8), dimension(xTpts, yTpts), intent(inout) :: dens, velx, vely, pres, ener
real(8), intent(in) :: dt
real(8), dimension(xTpts, yTpts) :: momx, momy, dens_n, &
momx_n, momy_n, ener_n
real(8), dimension(xTpts, yTpts) :: xr_plus, xu_plus, xv_plus, xp_plus
real(8), dimension(xTpts, yTpts) :: xr_minus, xu_minus, xv_minus, xp_minus
real(8), dimension(xTpts, yTpts) :: yr_plus, yu_plus, yv_plus, yp_plus
real(8), dimension(xTpts, yTpts) :: yr_minus, yu_minus, yv_minus, yp_minus
real(8), dimension(xTpts, yTpts) :: xrF, xruF, xrvF, xrwF, xeF
real(8), dimension(xTpts, yTpts) :: yrF, yruF, yrvF, yrwF, yeF
real(8), dimension(5) :: Uleft, Uright, Vleft, Vright, xF, yF
real(8), dimension(4) :: sgrav, sterm
integer :: i, j, k, l, m, n
#ifdef DEBUG_PAR
print *, "usegrav, grav = ", usegrav, grav
print *, "ilo, ihi, jlo, jhi = ", ilo, ihi, jlo, jhi
print *, "xTpts, yTpts, dx, dy = ", xTpts, yTpts, dx, dy
print *, "gamma = ", gamma
#endif
if (usegrav) then
sterm(:) = (/0.0e0, 0.0e0, grav, grav/)
else
sterm(:) = 0.0e0
end if
do j = 1, yTpts
do i = 1, xTpts
momx(i,j) = velx(i,j) * dens(i,j)
momy(i,j) = vely(i,j) * dens(i,j)
dens_n(i,j) = dens(i,j)
momx_n(i,j) = momx(i,j)
momy_n(i,j) = momy(i,j)
ener_n(i,j) = ener(i,j)
end do
end do
! rk2 has 2 steps
do k = 1, 2
call recon_getcellfaces(dens, velx, vely, pres, dt, &
!method, slimiter, &
xr_plus, xu_plus, xv_plus, xp_plus, &
xr_minus, xu_minus, xv_minus, xp_minus, &
yr_plus, yu_plus, yv_plus, yp_plus, &
yr_minus, yu_minus, yv_minus, yp_minus)
do j=jlo, jhi + 1
do i = ilo, ihi + 1
! initialize them to 0
xF(:) = 0.0
yF(:) = 0.0
Uleft = (/xr_plus(i-1,j), xu_plus(i-1,j), &
xv_plus(i-1,j), 0.0e0, xp_plus(i-1,j)/)
Uright = (/xr_minus(i,j), xu_minus(i,j), &
xv_minus(i,j), 0.0e0, xp_minus(i,j)/)
Vleft = (/yr_plus(i,j-1), yu_plus(i,j-1), &
yv_plus(i,j-1), 0.0e0, yp_plus(i,j-1)/)
Vright = (/yr_minus(i,j), yu_minus(i,j), &
yv_minus(i,j), 0.0e0, yp_minus(i,j)/)
if (to_upper(trim(flux_solver)) == "HLLC") then
!print *, "using HLLC"
call hllc(Uleft, Uright, "x", xF)
call hllc(Vleft, Vright, "y", yF)
else if (to_upper(trim(flux_solver)) == "HLLE") then
!print *, "using HLLE"
call hlle(Uleft, Uright, "x", xF)
call hlle(Vleft, Vright, "y", yF)
end if
xrF(i,j) = xF(1); xruF(i,j) = xF(2); xrvF(i,j) = xF(3); xeF(i,j) = xF(5)
yrF(i,j) = yF(1); yruF(i,j) = yF(2); yrvF(i,j) = yF(3); yeF(i,j) = yF(5)
end do
end do
do j=jlo, jhi
do i = ilo, ihi
sgrav = sterm * (/1.0e0, dens(i,j), dens(i,j), momy(i,j)/)
dens(i,j) = dens(i,j) &
+ (dt / dx) * (xrF(i,j) - xrF(i+1, j)) &
+ (dt / dy) * (yrF(i,j) - yrF(i, j+1)) &
+ dt * sgrav(1)
momx(i,j) = momx(i,j) &
+ (dt / dx) * (xruF(i,j) - xruF(i+1, j)) &
+ (dt / dy) * (yruF(i,j) - yruF(i, j+1)) &
+ dt * sgrav(2)
momy(i,j) = momy(i,j) &
+ (dt / dx) * (xrvF(i,j) - xrvF(i+1, j)) &
+ (dt / dy) * (yrvF(i,j) - yrvF(i, j+1)) &
+ dt * sgrav(3)
ener(i,j) = ener(i,j) &
+ (dt / dx) * (xeF(i,j) - xeF(i+1, j)) &
+ (dt / dy) * (yeF(i,j) - yeF(i, j+1)) &
+ dt * sgrav(4)
velx(i,j) = momx(i,j)/dens(i,j)
vely(i,j) = momy(i,j)/dens(i,j)
pres(i,j) = eos_getp((/dens(i,j), velx(i,j), vely(i,j), 0.0, ener(i,j)/))
! floor values for minimum is smallf
!dens(i,j) = max(dens(i,j), smallf)
!momx(i,j) = max(abs(momx(i,j)), smallf)
!momy(i,j) = max(abs(momy(i,j)), smallf)
!ener(i,j) = max(ener(i,j), smallf)
end do
end do
call applyBC_all(dens, velx, vely, pres, ener)
end do
do j=jlo, jhi
do i = ilo, ihi
dens(i,j) = 0.5e0 * (dens_n(i,j) + dens(i,j))
momx(i,j) = 0.5e0 * (momx_n(i,j) + momx(i,j))
momy(i,j) = 0.5e0 * (momy_n(i,j) + momy(i,j))
ener(i,j) = 0.5e0 * (ener_n(i,j) + ener(i,j))
! get primitive quantities
velx(i,j) = momx(i,j) / dens(i,j)
vely(i,j) = momy(i,j) / dens(i,j)
pres(i,j) = eos_getp((/dens(i,j), velx(i,j), vely(i,j), 0.0, ener(i,j)/))
end do
end do
end subroutine RK2_SSP
end module rk2_module