-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathscrapcode-ins-2015.f
109 lines (85 loc) · 3.22 KB
/
scrapcode-ins-2015.f
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
!! Code scraps: incompressible Navier-Stokes solver: 2015
!! Dr Jack R. Edwards Jr.
! ------------------------------------------------------------
! "j" direction flux balance (inviscid only)
!
! Note: beta2 stores beta**2
! Note: const scales dissipation terms
! Note: rho is density
! ------------------------------------------------------------
do i=2,ii
do j=1,jj
area = sqrt(del(i,j,2,1)**2 + del(i,j,2,2)**2)
xnx = del(i,j,2,1)/area
xny = del(i,j,2,2)/area
uave = 0.5*(u(i,j)+u(i,j+1))
vave = 0.5*(v(i,j)+v(i,j+1))
pave = 0.5*(p(i,j)+p(i,j+1))
udotn = uave*xnx + vave*xny
btave2 = 0.5*(beta2(i,j)+beta2(i,j+1))
eiga = 0.5*(abs(udotn) + sqrt(udotn**2 + 4.0*btave)
delp = p(i,j)-p(i,j+1) ! Note: first order model
c delp = pleft-pright ! higher-order model
delu = u(i,j)-u(i,j+1)
delv = v(i,j)-v(i,j+1)
dissip = 0.5*area*const*eiga*delp/btave2 ! multiply by 0.5 to make const < 1.0
xmassflux = area*rho*udotn
h(j,1) = xmassflux + dissip
h(j,2) = xmassflux*uave + xnx*area*pave
h(j,3) = xmassflux*vave + xny*area*pave
enddo
c add viscous flux calculation here
do j=1,jj
h(j,2) = h(j,2) -
h(j,3) = h(j,3) -
enddo
do j=2,jj
res(i,j,1) = res(i,j,1) + h(j,1) - h(j-1,1)
res(i,j,2) = res(i,j,2) + h(j,2) - h(j-1,2)
res(i,j,3) = res(i,j,3) + h(j,3) - h(j-1,3)
enddo
enddo
c----- minimum modulus averaging
function xmd(a,b) = 0.5*(sign(1.0,a)+sign(1.0,b))*min(abs(a),abs(b))
c----- left and right states using minimum modulus averaging ('j' direction)
do j=2,jj
du(j) = xmd(p(i,j+1)-p(i,j),p(i,j)-p(i,j-1))
enddo
du(1) = du(2)
du(jj+1) = du(jj)
pleft = p(i,j) + 0.5*du(j)
pright = p(i,j+1) - 0.5*du(j+1)
c ----- local time step (average areas to cell centers)
do j=2,jj
do i=2,ii
a11ave = 0.5*(del(i,j,1,1) + del(i-1,j,1,1))
a12ave = 0.5*(del(i,j,1,2) + del(i-1,j,1,2))
b11ave = 0.5*(del(i,j,2,1) + del(i,j-1,2,1))
b12ave = 0.5*(del(i,j,2,2) + del(i,j-1,2,2))
areaaave = sqrt(a11ave**2 + a12ave**2)
areabave = sqrt(b11ave**2 + b12ave**2)
vdotna = (u(i,j)*a11ave + v(i,j)*a12ave)/areaaave
vdotnb = (u(i,j)*b11ave + v(i,j)*b12ave)/areabave
eiga = 0.5*areaaave*(abs(vdotna) + sqrt(vdotna**2 + 4.0*beta2(i,j)))
eigb = 0.5*areabave*(abs(vdotnb) + sqrt(vdotnb**2 + 4.0*beta2(i,j)))
voldt = eiga + eigb
dt(i,j) = cfl*vol(i,j)/voldt !cfl is cfl number
enddo
enddo
c ---- solution update
do j=2,jj
do i=2,ii
dtv = dt(i,j)/vol(i,j)
deltap = -beta2(i,j)*dtv*res(i,j,1)
deltau = -dtv*res(i,j,2)/rho
deltav = -dtv*res(i,j,3)/rho
p(i,j) = p(i,j) + deltap
u(i,j) = u(i,j) + deltau
v(i,j) = v(i,j) + deltav
enddo
enddo
! Don't have here:
! - Zeroing out of mass-flux terms in inviscid flux
! - Boundary conditions, ie, ghost cell values
! - Calculation of beta
! - Viscous fluxes