-
Notifications
You must be signed in to change notification settings - Fork 0
/
hsolve.f90
129 lines (115 loc) · 4.34 KB
/
hsolve.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
subroutine hsolve(h,oldh,hdt,dtime)
! --------------
USE header, only : NI,NJ,rc_kind
! modified for periodicew boundaries
! ----------------------------------
integer iter,i,j,kount,l
REAL(kind=rc_kind) :: ch(9,NI,NJ),rhs(NI,NJ), &
& oldh(0:NI+1,0:NJ+1),h(0:NI+1,0:NJ+1),res, &
& hdt(0:NI+1,0:NJ+1),dti,dtime,tol,maxres,hstar,rlx
! & dti,dtime,tol,maxres,hstar,rlx
tol= 1.d-12
dti= 1.d0/dtime
kount=0
! compute fine grid coefficients and rhs on fine grid
call chfine(dtime,ch,rhs)
rlx= 1.72
do iter=1,100
do j=1,NJ
i=1
hstar= ( ch(2,i,j)*h(i+1,j) &
& +ch(3,i,j)*h(NI,j) &
& +ch(4,i,j)*h(i,j+1) &
& +ch(5,i,j)*h(i,j-1) &
& +ch(6,i,j)*h(i+1,j+1) &
& +ch(7,i,j)*h(NI,j+1) &
& +ch(8,i,j)*h(i+1,j-1) &
& +ch(9,i,j)*h(NI,j-1) &
& - rhs(i,j) )/(-ch(1,i,j))
h(i,j)= (1.d0 -rlx)*h(i,j) +rlx*hstar
do i=2,NI-1
hstar= ( ch(2,i,j)*h(i+1,j) &
& +ch(3,i,j)*h(i-1,j) &
& +ch(4,i,j)*h(i,j+1) &
& +ch(5,i,j)*h(i,j-1) &
& +ch(6,i,j)*h(i+1,j+1) &
& +ch(7,i,j)*h(i-1,j+1) &
& +ch(8,i,j)*h(i+1,j-1) &
& +ch(9,i,j)*h(i-1,j-1) &
& - rhs(i,j) )/(-ch(1,i,j))
h(i,j)= (1.d0 -rlx)*h(i,j) +rlx*hstar
enddo ! i
i=NI
hstar= ( ch(2,i,j)*h(1,j) &
& +ch(3,i,j)*h(i-1,j) &
& +ch(4,i,j)*h(i,j+1) &
& +ch(5,i,j)*h(i,j-1) &
& +ch(6,i,j)*h(1,j+1) &
& +ch(7,i,j)*h(i-1,j+1) &
& +ch(8,i,j)*h(1,j-1) &
& +ch(9,i,j)*h(i-1,j-1) &
& - rhs(i,j) )/(-ch(1,i,j))
h(i,j)= (1.d0 -rlx)*h(i,j) +rlx*hstar
enddo ! j
do l=1,3
call hfill(dtime,h)
enddo
maxres= 0.d0
do j=1,NJ
i=1
res= rhs(i,j)- &
& (ch(1,i,j)*h(i,j) &
& +ch(2,i,j)*h(i+1,j) &
& +ch(3,i,j)*h(NI,j) &
& +ch(4,i,j)*h(i,j+1) &
& +ch(5,i,j)*h(i,j-1) &
& +ch(6,i,j)*h(i+1,j+1) &
& +ch(7,i,j)*h(NI,j+1) &
& +ch(8,i,j)*h(i+1,j-1) &
& +ch(9,i,j)*h(NI,j-1))
if (abs(res).gt.maxres) then
maxres= abs(res)
end if
do i=2,NI-1 ! res = b - A x'
res= rhs(i,j)- &
& (ch(1,i,j)*h(i,j) &
& +ch(2,i,j)*h(i+1,j) &
& +ch(3,i,j)*h(i-1,j) &
& +ch(4,i,j)*h(i,j+1) &
& +ch(5,i,j)*h(i,j-1) &
& +ch(6,i,j)*h(i+1,j+1) &
& +ch(7,i,j)*h(i-1,j+1) &
& +ch(8,i,j)*h(i+1,j-1) &
& +ch(9,i,j)*h(i-1,j-1))
if (abs(res).gt.maxres) then
maxres= abs(res)
end if
enddo !i
i=NI
res= rhs(i,j)- &
& (ch(1,i,j)*h(i,j) &
& +ch(2,i,j)*h(1,j) &
& +ch(3,i,j)*h(i-1,j) &
& +ch(4,i,j)*h(i,j+1) &
& +ch(5,i,j)*h(i,j-1) &
& +ch(6,i,j)*h(1,j+1) &
& +ch(7,i,j)*h(i-1,j+1) &
& +ch(8,i,j)*h(1,j-1) &
& +ch(9,i,j)*h(i-1,j-1))
if (abs(res).gt.maxres) then
maxres= abs(res)
end if
enddo !j
! write(6,*) 'iter,maxres',iter,maxres
if (maxres.gt.3000.d0) then
print*, 'hsolve.f90: STOP. res too large, i,j,maxres=',i,j,maxres
stop
end if
if (maxres.lt.tol) goto 13
enddo !iter
13 continue
do l=1,1
call mprove(h,ch,rhs,dtime)
end do
return
END