-
Notifications
You must be signed in to change notification settings - Fork 1
/
grid_reconstruct.f90
143 lines (118 loc) · 5.41 KB
/
grid_reconstruct.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
subroutine grid_reconstruct(nx,ny,rho_floor,&
rho,rho_left_x,rho_right_x,rho_left_y,rho_right_y,&
tau,tau_left_x,tau_right_x,tau_left_y,tau_right_y,&
etot,etot_left_x,etot_right_x,etot_left_y,etot_right_y,&
mom_A,mom_A_left_x,mom_A_right_x,mom_A_left_y,mom_A_right_y,&
mom_B,mom_B_left_x,mom_B_right_x,mom_B_left_y,mom_B_right_y,&
mom_x,mom_x_left_x,mom_x_right_x,mom_x_left_y,mom_x_right_y,&
mom_y,mom_y_left_x,mom_y_right_x,mom_y_left_y,mom_y_right_y,&
x,x_left_x,x_right_x,x_left_y,x_right_y,&
y,y_left_x,y_right_x,y_left_y,y_right_y,&
phi,phi_left_x,phi_right_x,phi_left_y,phi_right_y,recons)
implicit none
integer :: recons
include 'variables.h'
include 'grid.h'
double precision :: mom_A_2(nx,ny),mom_B_2(nx,ny),mom_x_2(nx,ny),mom_y_2(nx,ny)
call getrtheta(nx,ny,x,y,r,theta)
!divide conserved quantities by rho
mom_A_2=mom_A/rho
mom_B_2=mom_B/rho
mom_x_2=mom_x/rho
mom_y_2=mom_y/rho
if (recons .eq. 1) then
! print*, 'using ppm!'
call reconstruct_ppm(nx,ny,rho,rho_left_x,rho_right_x,rho_left_y,rho_right_y)
call reconstruct_ppm(nx,ny,tau,tau_left_x,tau_right_x,tau_left_y,tau_right_y)
call reconstruct_ppm(nx,ny,etot,etot_left_x,etot_right_x,etot_left_y,etot_right_y)
call reconstruct_ppm(nx,ny,mom_A_2,mom_A_left_x,mom_A_right_x,mom_A_left_y,mom_A_right_y)
call reconstruct_ppm(nx,ny,mom_B_2,mom_B_left_x,mom_B_right_x,mom_B_left_y,mom_B_right_y)
call reconstruct_ppm(nx,ny,mom_x_2,mom_x_left_x,mom_x_right_x,mom_x_left_y,mom_x_right_y)
call reconstruct_ppm(nx,ny,mom_y_2,mom_y_left_x,mom_y_right_x,mom_y_left_y,mom_y_right_y)
else if (recons .eq. 0) then
! print*, 'using minmod!'
call reconstruct(nx,ny,rho,rho_left_x,rho_right_x,rho_left_y,rho_right_y)
call reconstruct(nx,ny,tau,tau_left_x,tau_right_x,tau_left_y,tau_right_y)
call reconstruct(nx,ny,etot,etot_left_x,etot_right_x,etot_left_y,etot_right_y)
call reconstruct(nx,ny,mom_A_2,mom_A_left_x,mom_A_right_x,mom_A_left_y,mom_A_right_y)
call reconstruct(nx,ny,mom_B_2,mom_B_left_x,mom_B_right_x,mom_B_left_y,mom_B_right_y)
call reconstruct(nx,ny,mom_x_2,mom_x_left_x,mom_x_right_x,mom_x_left_y,mom_x_right_y)
call reconstruct(nx,ny,mom_y_2,mom_y_left_x,mom_y_right_x,mom_y_left_y,mom_y_right_y)
else
print*, 'error invalid reconstruction!'
call exit(1)
end if
do j=2,ny
do i=2,nx
x_left_x(i,j) = ( x(i,j) + x(i-1,j) )/2d0
x_right_x(i,j) = ( x(i,j) + x(i-1,j) )/2d0
x_left_y(i,j) = ( x(i,j) + x(i,j-1) )/2d0
x_right_y(i,j) = ( x(i,j) + x(i,j-1) )/2d0
y_left_x(i,j) = ( y(i,j) + y(i-1,j) )/2d0
y_right_x(i,j) = ( y(i,j) + y(i-1,j) )/2d0
y_left_y(i,j) = ( y(i,j) + y(i,j-1) )/2d0
y_right_y(i,j) = ( y(i,j) + y(i,j-1) )/2d0
end do
end do
do i=1,nx
j=1
x_left_x(i,j) = ( x(i,j) + x(i,j) )/2d0
x_right_x(i,j) = ( x(i,j) + x(i,j) )/2d0
x_left_y(i,j) = ( x(i,j) + x(i,j) )/2d0
x_right_y(i,j) = ( x(i,j) + x(i,j) )/2d0
y_left_x(i,j) = ( y(i,j) + y(i,j) )/2d0
y_right_x(i,j) = ( y(i,j) + y(i,j) )/2d0
y_left_y(i,j) = ( y(i,j) + y(i,j) )/2d0
y_right_y(i,j) = ( y(i,j) + y(i,j) )/2d0
end do
do j=1,ny
i=1
x_left_x(i,j) = ( x(i,j) + x(i,j) )/2d0
x_right_x(i,j) = ( x(i,j) + x(i,j) )/2d0
x_left_y(i,j) = ( x(i,j) + x(i,j) )/2d0
x_right_y(i,j) = ( x(i,j) + x(i,j) )/2d0
y_left_x(i,j) = ( y(i,j) + y(i,j) )/2d0
y_right_x(i,j) = ( y(i,j) + y(i,j) )/2d0
y_left_y(i,j) = ( y(i,j) + y(i,j) )/2d0
y_right_y(i,j) = ( y(i,j) + y(i,j) )/2d0
end do
!multiplying reconstructed values by rho to return to conserved variables
mom_A_left_x = mom_A_left_x*rho_left_x
mom_A_right_x = mom_A_right_x*rho_right_x
mom_A_left_y = mom_A_left_y*rho_left_y
mom_A_right_y = mom_A_right_y*rho_right_y
call getrtheta(nx,ny,x_left_x,y_left_x,r_left_x,theta)
call getrtheta(nx,ny,x_right_x,y_right_x,r_right_x,theta)
call getrtheta(nx,ny,x_left_y,y_left_y,r_left_y,theta)
call getrtheta(nx,ny,x_right_y,y_right_y,r_right_y,theta)
mom_B_left_x = mom_B_left_x * rho_left_x
mom_B_right_x = mom_B_right_x * rho_right_x
mom_B_left_y = mom_B_left_y * rho_left_y
mom_B_right_y = mom_B_right_y * rho_right_y
!!$ mom_B_left_x = mom_B_left_x * rho_left_x * r_left_x
!!$ mom_B_right_x = mom_B_right_x * rho_right_x * r_right_x
!!$ mom_B_left_y = mom_B_left_y * rho_left_y * r_left_y
!!$ mom_B_right_y = mom_B_right_y * rho_right_y * r_right_y
mom_x_left_x = mom_x_left_x*rho_left_x
mom_x_right_x = mom_x_right_x*rho_right_x
mom_x_left_y = mom_x_left_y*rho_left_y
mom_x_right_y = mom_x_right_y*rho_right_y
mom_y_left_x = mom_y_left_x*rho_left_x
mom_y_right_x = mom_y_right_x*rho_right_x
mom_y_left_y = mom_y_left_y*rho_left_y
mom_y_right_y = mom_y_right_y*rho_right_y
! call reconstruct(nx,ny,x,x_left_x,x_right_x,x_left_y,x_right_y)
! call reconstruct(nx,ny,y,y_left_x,y_right_x,y_left_y,y_right_y)
! call reconstruct(nx,ny,phi,phi_left_x,phi_right_x,phi_left_y,phi_right_y)
!density floor
!!!! !$OMP PARALLEL DO PRIVATE(i,j)
do j=1,ny
do i=1,nx
rho_left_x(i,j) = max(rho_left_x(i,j),rho_floor)
rho_left_y(i,j) = max(rho_left_y(i,j),rho_floor)
rho_right_x(i,j) = max(rho_right_x(i,j),rho_floor)
rho_right_y(i,j) = max(rho_right_y(i,j),rho_floor)
end do
end do
!!!!! !$OMP END PARALLEL DO
end subroutine grid_reconstruct