-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathscrapcodeb-2015.f
112 lines (85 loc) · 3.62 KB
/
scrapcodeb-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
110
111
! A few code scraps (PLEASE CHECK FOR ERRORS!!!!)
! ------------------------------------------------------------
! "j" direction flux balance (parallelogram CV)
! ------------------------------------------------------------
do i=2,ii
do j=1,jj
f1 = u(i,j+1)
f3 = u(i,j)
f2 = 0.25*(u(i,j)+u(i,j+1)+u(i+1,j+1)+u(i+1,j))
f4 = 0.25*(u(i,j)+u(i,j+1)+u(i-1,j+1)+u(i-1,j))
volave = 0.5*(vol(i,j) + vol(i,j+1))
deln(1) = yc(i,j+1) - yc(i,j)
deln(2) = -(xc(i,j+1)-xc(i,j))
gradf(1:2) = (f1-f3)*area(i,j,2,1:2) + (f2-f4)*deln(1:2) ! Green's theorem
gradf(1:2) = gradf(1:2)/volave !divide by average volume
h(j) = gradf(1)*area(i,j,2,1) + gradf(2)*area(i,j,2,2) !grad u * n * area
enddo
do j=2,jj
res(i,j) = res(i,j) + h(j)-h(j-1) ! note augmentation of residual
enddo
enddo
! ---------------------------------------------------------------------------
! "j" direction flux balance (normal / tangential gradient decomposition
! ---------------------------------------------------------------------------
do i=2,ii
do j=1,jj
f1 = u(i,j+1)
f3 = u(i,j)
delt(1) = xc(i,j+1) - xc(i,j)
delt(2) = yc(i,j+1) - yc(i,j)
dels = sqrt(delt(1)**2 + delt(2)**2)
tau(1) = delt(1)/dels !unit vector pointing from j to j+1
tau(2) = delt(2)/dels
areaface = sqrt(area(i,j,2,1)**2+area(i,j,2,2)**2)
xn(1) = area(i,j,2,1)/areaface !normal vector to mesh face
xn(2) = area(i,j,2,2)/areaface
taudotn = tau(1)*xn(1) + tau(2)*xn(2) !dot product of tau and n
gradave(1:2) = 0.5*(grad(i,j,1:2) + grad(i,j+1,1:2)) !average gradient vector
graddottau = (gradave(1)*tau(1) + gradave(2)*tau(2))*taudotn
gradf(1:2) = (f1-f3)*taudotn*xn(1:2)/dels + gradave(1:2) - graddottau*xn(1:2) !gradient
h(j) = gradf(1)*area(i,j,2,1) + gradf(2)*area(i,j,2,2) !grad u * n * area
enddo
do j=2,jj
res(i,j) = res(i,j) + h(j)-h(j-1) ! note augmentation of residual
enddo
enddo
! -----------------------------------------------------------------------------
! Matrix calculations (based on 'thin layer form' of parallel CV model
! Note that you only have to do this once (before the iterations)
! ----------------------------------------------------------------------------
do j=2,jj
do i=1,ii
volave = 0.5*(vol(i,j) + vol(i+1,j))
ahalf(i) = (area(i,j,1,1)**2 + area(i,j,1,2)**2)/volave
enddo
do i=2,ii
a(i,j,1) = ahalf(i-1)
a(i,j,5) = ahalf(i)
a(i,j,3) = -(ahalf(i)+ahalf(i-1))
enddo
enddo
do i=2,ii
do j=1,jj
volave = 0.5*(vol(i,j) + vol(i,j+1))
bhalf(j) = (area(i,j,2,1)**2 + area(i,j,2,2)**2)/volave
enddo
do j=2,jj
a(i,j,2) = bhalf(j-1)
a(i,j,4) = bhalf(j)
a(i,j,3) = a(i,j,3) -(bhalf(j)+bhalf(j-1))
enddo
enddo
! ------------------------------------------------------------
! Cell centered gradients using Green's theorem
! ------------------------------------------------------------
do j=2,jj
do i=2,ii
f1 = 0.5*(u(i,j)+u(i+1,j))
f2 = 0.5*(u(i,j)+u(i,j+1))
f3 = 0.5*(u(i,j)+u(i-1,j))
f4 = 0.5*(u(i,j)+u(i,j-1))
grad(i,j,:) = f1*area(i,j,1,:) - f3*area(i-1,j,1,:) + f2*area(i,j,2,:) - f4*area(i,j-1,2,:)
grad(i,j,:) = grad(i,j,:)/vol(i,j)
enddo
enddo