resc2.f
5.53 KB

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

```
c
c file resc2.f
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c * *
c * copyright (c) 2008 by UCAR *
c * *
c * University Corporation for Atmospheric Research *
c * *
c * all rights reserved *
c * *
c * MUDPACK version 5.0.1 *
c * *
c * A Fortran Package of Multigrid *
c * *
c * Subroutines and Example Programs *
c * *
c * for Solving Elliptic Partial Differential Equations *
c * *
c * by *
c * *
c * John Adams *
c * *
c * of *
c * *
c * the National Center for Atmospheric Research *
c * *
c * Boulder, Colorado (80307) U.S.A. *
c * *
c * which is sponsored by *
c * *
c * the National Science Foundation *
c * *
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c subroutine resc2(nx,ny,work,res)
c
c ... purpose
c
c subroutine resc2 computes the fine grid residual in the nx by ny array
c res after calling cud2 or cuh2 or cud2sa. if
c
c L * p = f
c
c is the n by n (n = nx*ny) block tri-diagonal linear system resulting
c from the pde discretization (done internally in the mudpack solver)
c then resc2 computes the nx by ny residual array
c
c res = f - L * p
c
c one of the vector norms of res,
c
c || res ||
c
c can be computed as a "measure" of how well the approximation satisfies
c the discretization equations.
c
c ... see documentation and test files provided in this distribution
c
c ... notes
c
c let pe be the exact continuous solution to the elliptic pde
c evaluated on the nx by ny discretization grid
c
c let p be the exact solution to the linear discretization
c
c let phi be the approximation to p generated by the mudpack solver
c
c then discretization level error is defined by the condition
c
c || phi - p || < || p - pe ||.
c =
c
c a common measure of multigrid efficieny is that discretization level
c error is reached in one full multigrid cycle (see references [2,9] in
c the mudpack file "readme"). this can happen before the residual is
c reduced to the level of roundoff error. consequently, || res || is
c a conservative measure of accuracy which can be wasteful if multi-
c grid cycles are executed until it reaches the level of roundoff error.
c
c || res || can be used to estimate the convergence rate of multigrid
c iteration. let r(n) be the residual and e(n) be the error after
c executing n cycles. they are related by the residual equation
c
c L * e(n) = r(n).
c
c it follows that the ratio
c
c || r(n+1) || / || r(n) ||
c
c estimates
c
c || e(n+1) || / || e(n) ||
c
c which in turn estimates the convergence rate
c
c c = max || e(k+1) || / || e(k) ||.
c k
c
c notice
c n
c || e(n) || < c || e(0) ||.
c
c
c ... assumptions (see cud2 or cuh2 or cud2sa documentation)
c
c (1) nx,ny have the same values as iparm(10),iparm(11) (used
c to set the fine grid resolution)
c
c (2) work is the same argument used in calling the solver.
c
c (3) work have not changed since the last call to the solver.
c
c (3) assures a copy of the last approximation phi is contained in work.
c if these assumptions are not true then resc2 cannot compute the
c residual in res.
c
subroutine resc2(nx,ny,work,res)
implicit none
integer nx,ny,ic
complex work(*),res(nx,ny)
c
c set approximation and coefficient pointers
c
ic = 1+(nx+2)*(ny+2)
call rec2(nx,ny,work,work(ic),res)
return
end
subroutine rec2(nx,ny,phi,cof,res)
implicit none
integer nx,ny,i,j
complex phi(0:nx+1,0:ny+1),cof(nx,ny,6),res(nx,ny)
c
c compute residual over entire x,y grid
c
do j=1,ny
do i=1,nx
res(i,j) = cof(i,j,6)-(
+ cof(i,j,1)*phi(i-1,j)+
+ cof(i,j,2)*phi(i+1,j)+
+ cof(i,j,3)*phi(i,j-1)+
+ cof(i,j,4)*phi(i,j+1)+
+ cof(i,j,5)*phi(i,j))
end do
end do
return
end
```