Blame view

src/resm3sp.f 6.92 KB
af19620a   Kévin Duquette   Add all files
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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
c
c     file resm3sp.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 resm3sp(nx,ny,nz,nxa,nxb,nyc,nyd,nze,nzf,work,res)
c
c
c ... purpose
c
c
c     subroutine resm3sp computes the fine grid residual in the nx by ny by nz
c     array res after calling mud3sp.  if
c
c          l * p = f
c
c     is the n by n (n = nx*ny*nz) block tri-diagonal linear system resulting
c     from the pde discretization (done internally in mud3sp) and phi is the
c     approximation to p obtained by calling mud3sp, then resm3sp computes
c     the nx by ny by nz residual array
c
c          res = f - l * phi.
c
c     one of the vector norms of res,
c
c          || res ||
c
c     can be computed as a "measure" of how well phi satisfies the
c     discretization equations.  for example, the following statements
c     will compute the location and size of the maximum residual in res
c     on cray computers:
c
c          ijk = isamax(nx*ny*nz,res,1)
c
c          kmax = (ijk-1)/(nx*ny) + 1
c
c          jmax = (ijk-(kmax-1)*nx*ny-1)/nx + 1
c
c          imax = ij - nx*((kmax-1)*ny-jmax+1)
c
c          resmax = abs(res(imax,jmax,kmax))
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 by nz 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 mud3sp documentation)
c
c     (1) nx,ny,nz have the same values as iparm(14),iparm(15),iparm(16)
c         (used to set the fine grid resolution when calling mud3sp)
c
c     (2) nxa,nxb,nyc,nyd,nze,nzf have the same values as iparm(2),
c         iparm(3),iparm(4),iparm(5),iparm(6),iparm(7) (used to flag
c         boundary conditions when calling mud3sp).
c
c     (3) work,phi are the same parameters used in calling mud3sp.
c
c     (4) work,phi have not changed since the last call to mud3sp.
c
c     (3) assures a copy of the last approximation phi is contained in work.
c     if (1)-(4) are not true then resm3sp cannot compute the residual.
c
      subroutine resm3sp(nx,ny,nz,nxa,nxb,nyc,nyd,nze,nzf,wk,res)
      implicit none
      integer nx,ny,nz,nxa,nxb,nyc,nyd,nze,nzf,ir,ix,iy,iz
      real wk(*),res(*)
c
c     set pointers
c
      ir = 1+(nx+2)*(ny+2)*(nz+2)
      ix = ir+nx*ny*nz
      iy = ix+3*nx
      iz = iy+3*ny
      call rem3sp(nx,ny,nz,nxa,nxb,nyc,nyd,nze,nzf,wk,wk(ir),wk(ix),
     +            wk(iy),wk(iz),res)
      return
      end

      subroutine rem3sp(nx,ny,nz,nxa,nxb,nyc,nyd,nze,nzf,phi,
     +                  rhs,cofx,cofy,cofz,resf)
      implicit none
      integer nx,ny,nz,nxa,nxb,nyc,nyd,nze,nzf
      integer i,j,k,ist,ifn,jst,jfn,kst,kfn
      real cofx(nx,3),cofy(ny,3),cofz(nz,3)
      real rhs(nx,ny,nz),phi(0:nx+1,0:ny+1,0:nz+1),resf(nx,ny,nz)
c
c     intialize residual to zero and set loop limits
c
      do k=1,nz
	do j=1,ny
	  do i=1,nx
	    resf(i,j,k) = 0.0
	  end do
	end do
      end do
      ist = 1
      if (nxa.eq.1) ist = 2
      ifn = nx
      if (nxb.eq.1) ifn = nx-1
      jst = 1
      if (nyc.eq.1) jst = 2
      jfn = ny
      if (nyd.eq.1) jfn = ny-1
      kst = 1
      if (nze.eq.1) kst = 2
      kfn = nz
      if (nzf.eq.1) kfn = nz-1
c
c     compute residual
c
      do k=kst,kfn
	do j=jst,jfn
	  do i=ist,ifn
	    resf(i,j,k) =  rhs(i,j,k)-(
     +      cofx(i,1)*phi(i-1,j,k)+cofx(i,2)*phi(i+1,j,k) +
     +      cofy(j,1)*phi(i,j-1,k)+cofy(j,2)*phi(i,j+1,k) +
     +      cofz(k,1)*phi(i,j,k-1)+cofz(k,2)*phi(i,j,k+1) +
     +      (cofx(i,3)+cofy(j,3)+cofz(k,3))*phi(i,j,k)  )
	  end do
	end do
      end do
      return
      end