Line data Source code
1 : module edyn_mudmod
2 : use shr_kind_mod, only: r8 => shr_kind_r8
3 : use cam_logfile, only: iulog
4 : use edyn_mud, only: dismd2cr, mud2cr1, adjmd2cr, kcymd2cr, relmd2cr, resmd2cr
5 : use edyn_mudcom, only: swk2, trsfc2, prolon2, cor2, res2
6 :
7 : implicit none
8 :
9 : private
10 :
11 : public :: mudmod
12 :
13 : contains
14 : !-----------------------------------------------------------------------
15 19 : subroutine mudmod(pe,phi_out,jntl,isolve,nlev,ier)
16 : use edyn_solver_coefs,only: cee
17 : use edyn_params, only: pi
18 :
19 : integer,intent(in) :: jntl, isolve, nlev
20 : integer,intent(out) :: ier ! output: not converged ier < 0
21 : !
22 : ! set grid size params
23 : !
24 : integer iixp,jjyq,iiex,jjey,nnx,nny,llwork
25 : parameter (iixp = 5 , jjyq = 3)
26 : !
27 : ! estimate work space for point relaxation (see mud2cr.d)
28 : !
29 : real(r8) :: phi_out(0:iixp*2**(nlev-1)+1+1,0:jjyq*2**(nlev-1)+1+1)
30 19 : real(r8),allocatable :: phi(:,:),rhs(:,:),work(:)
31 : !
32 : ! put integer and floating point argument names in contiguous
33 : ! storage for labelling in vectors iprm,fprm
34 : !
35 : integer iprm(17),mgopt(4)
36 : real(r8) :: fprm(6)
37 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nx,ny,&
38 : iguess,maxcy,method,nwork,lwrkqd,itero
39 : common/itmud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nx,ny, &
40 : iguess,maxcy,method,nwork,lwrkqd,itero
41 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
42 : common/ftmud2cr/xa,xb,yc,yd,tolmax,relmax
43 : equivalence(intl,iprm)
44 : equivalence(xa,fprm)
45 : integer i,j,ierror
46 : real(r8) :: PE(iixp*2**(nlev-1)+1,*)
47 : integer, parameter :: maxcya=50
48 : integer jj,jjj,ij
49 :
50 19 : iiex = nlev
51 19 : jjey = nlev
52 19 : nnx=iixp*2**(iiex-1)+1
53 19 : nny=jjyq*2**(jjey-1)+1
54 19 : llwork=(7*(nnx+2)*(nny+2)+76*nnx*nny)/3
55 :
56 152 : allocate(phi(nnx,nny),rhs(nnx,nny),work(llwork))
57 : !
58 : ! SET INPUT INTEGER PARAMETERS
59 : !
60 19 : INTL = JNTL
61 : !
62 : ! set boundary condition flags
63 : !
64 19 : nxa = 0
65 19 : nxb = 0
66 19 : nyc = 2
67 19 : nyd = 1
68 : !
69 : ! set grid sizes from parameter statements
70 : !
71 19 : ixp = iixp
72 19 : jyq = jjyq
73 19 : iex = iiex
74 19 : jey = jjey
75 19 : nx = nnx
76 19 : ny = nny
77 : !
78 : ! set multigrid arguments (w(2,1) cycling with fully weighted
79 : ! residual restriction and cubic prolongation)
80 : !
81 19 : mgopt(1) = 2
82 19 : mgopt(2) = 3
83 19 : mgopt(3) = 2
84 19 : mgopt(4) = 3
85 : !
86 : ! set for one cycle
87 : !
88 19 : maxcy = maxcya
89 : !
90 : ! set no initial guess forcing full multigrid cycling
91 : !
92 19 : iguess = 0
93 : !
94 : ! set work space length approximation from parameter statement
95 : !
96 19 : nwork = llwork
97 : !
98 : ! set line z relaxation
99 : !
100 19 : method = 3
101 : !
102 : ! set end points of solution rectangle in (x,y) space
103 : !
104 19 : xa = -pi
105 19 : xb = pi
106 19 : yc = 0.0_r8
107 19 : yd = 0.5_r8*pi
108 : !
109 : ! set error control flag
110 : !
111 19 : tolmax = 0.001_r8
112 : !
113 : ! set right hand side in rhs
114 : ! initialize phi to zero
115 : !
116 1558 : do i=1,nx
117 76969 : do j=1,ny
118 150822 : rhs(i,j) = cee(i+(j-1)*nx+9*nx*ny)
119 76950 : phi(i,j) = 0.0_r8
120 : end do
121 : end do
122 : !
123 : ! set specified boundaries in phi
124 : !
125 1558 : do i=1,nx
126 3097 : phi(i,ny) = rhs(i,ny)/cee(i+(ny-1)*nx+8*nx*ny)
127 : end do
128 :
129 19 : call mud2cm(iprm,fprm,work,rhs,phi,mgopt,ierror,isolve)
130 :
131 : !
132 : ! attempt solution
133 : !
134 19 : intl = 1
135 :
136 19 : call mud2cm(iprm,fprm,work,rhs,phi,mgopt,ierror,isolve)
137 19 : ier = ierror ! ier < 0 not converged
138 19 : if(ier < 0 ) goto 108
139 : !
140 : ! COPY PHI TO PE
141 : !
142 950 : DO J = 1,NY
143 931 : JJ = NY+J-1
144 931 : JJJ = NY+1-J
145 76361 : DO I = 1,NX
146 75411 : PE(I,JJ) = PHI(I,J)
147 76342 : PE(I,JJJ) = PHI(I,J)
148 : END DO
149 : END DO
150 :
151 : ! am 8/10 for calculating residual: convert work array (solution) into array
152 : ! sized as coefficient stencil (c0, cofum) including values at index 0, nmlon0+1
153 : ! and nmlat0+1
154 :
155 988 : do j=0,ny+1
156 969 : jj = j*(nx+2)
157 81415 : do i=0,nx+1
158 80427 : ij = jj+i+1
159 81396 : phi_out(i,j) = work(ij)
160 : end do
161 : end do
162 :
163 : 108 continue
164 :
165 19 : deallocate(phi,rhs,work)
166 :
167 19 : end subroutine mudmod
168 : !-------------------------------------------------------------------
169 38 : subroutine mud2cm(iparm,fparm,work,rhs,phi,mgopt,ierror,isolve)
170 :
171 : integer,intent(in) :: isolve
172 : integer iparm,mgopt,ierror
173 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
174 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
175 : kcycle,iprer,ipost,intpol,kps
176 : real(r8) :: fparm,xa,xb,yc,yd,tolmax,relmax
177 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
178 : integer int,iw,k,kb,nx,ny,ic,itx,ity
179 : dimension iparm(17),fparm(6),mgopt(4)
180 : real(r8) :: work(*),phi(*),rhs(*)
181 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
182 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
183 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
184 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
185 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), &
186 : nxk(50),nyk(50),isx,jsy
187 :
188 : data int / 0 /
189 : save int
190 :
191 38 : ierror = 1
192 38 : intl = iparm(1) ! set and check intl on all calls
193 38 : if (intl*(intl-1).ne.0) return
194 38 : if (int.eq.0) then
195 2 : int = 1
196 2 : if (intl.ne.0) return ! very first call is not intl=0
197 : end if
198 38 : ierror = 0
199 : !
200 : ! set arguments internally
201 : ! these will not be rechecked if intl=1!
202 : !
203 38 : nxa = iparm(2)
204 38 : nxb = iparm(3)
205 38 : nyc = iparm(4)
206 38 : nyd = iparm(5)
207 38 : ixp = iparm(6)
208 38 : jyq = iparm(7)
209 38 : iex = iparm(8)
210 38 : jey = iparm(9)
211 38 : ngrid = max0(iex,jey)
212 38 : nfx = iparm(10)
213 38 : nfy = iparm(11)
214 38 : iguess = iparm(12)
215 38 : maxcy = iparm(13)
216 38 : method = iparm(14)
217 38 : nwork = iparm(15)
218 38 : kcycle = mgopt(1)
219 38 : if (kcycle .eq. 0) then
220 : ! set defaults
221 0 : kcycle = 2
222 0 : iprer = 2
223 0 : ipost = 1
224 0 : intpol = 3
225 : else
226 38 : iprer = mgopt(2)
227 38 : ipost = mgopt(3)
228 38 : intpol = mgopt(4)
229 : end if
230 38 : xa = fparm(1)
231 38 : xb = fparm(2)
232 38 : yc = fparm(3)
233 38 : yd = fparm(4)
234 38 : tolmax = fparm(5)
235 38 : if (intl .eq. 0) then ! intialization call
236 : !
237 : ! check input arguments
238 : !
239 19 : ierror = 2 ! check boundary condition flags
240 19 : if (max0(nxa,nxb,nyc,nyd).gt.2) return
241 19 : if (min0(nxa,nxb,nyc,nyd).lt.0) return
242 19 : if (nxa.eq.0.and.nxb.ne.0) return
243 19 : if (nxa.ne.0.and.nxb.eq.0) return
244 19 : if (nyc.eq.0.and.nyd.ne.0) return
245 19 : if (nyc.ne.0.and.nyd.eq.0) return
246 19 : ierror = 3 ! check grid sizes
247 19 : if (ixp.lt.2) return
248 19 : if (jyq.lt.2) return
249 19 : ierror = 4
250 19 : ngrid = max0(iex,jey)
251 19 : if (iex.lt.1) return
252 19 : if (jey.lt.1) return
253 19 : if (ngrid.gt.50) return
254 19 : ierror = 5
255 19 : if (nfx.ne.ixp*2**(iex-1)+1) return
256 19 : if (nfy.ne.jyq*2**(jey-1)+1) return
257 19 : ierror = 6
258 19 : if (iguess*(iguess-1).ne.0) return
259 19 : ierror = 7
260 19 : if (maxcy.lt.1) return
261 19 : ierror = 8
262 19 : if (method.lt.0 .or. method.gt.3) return
263 19 : ierror = 9
264 : ! compute and test minimum work space
265 19 : isx = 0
266 19 : if (method.eq.1 .or. method.eq.3) then
267 0 : if (nxa.ne.0) isx = 3
268 19 : if (nxa.eq.0) isx = 5
269 : end if
270 19 : jsy = 0
271 19 : if (method.eq.2 .or. method.eq.3) then
272 19 : if (nyc.ne.0) jsy = 3
273 19 : if (nyc.eq.0) jsy = 5
274 : end if
275 19 : kps = 1
276 114 : do k=1,ngrid
277 : ! set subgrid sizes
278 190 : nxk(k) = ixp*2**(max0(k+iex-ngrid,1)-1)+1
279 95 : nyk(k) = jyq*2**(max0(k+jey-ngrid,1)-1)+1
280 95 : nx = nxk(k)
281 95 : ny = nyk(k)
282 19 : kps = kps+(nx+2)*(ny+2)+nx*ny*(10+isx+jsy)
283 : end do
284 19 : iparm(16) = kps+(nfx+2)*(nfy+2) ! exact minimum work space
285 19 : lwork = iparm(16)
286 19 : if (lwork .gt. nwork) return
287 19 : ierror = 10 ! check solution region
288 19 : if (xb.le.xa .or. yd.le.yc) return
289 19 : ierror = 11
290 19 : if (tolmax .lt. 0.0_r8) return
291 19 : ierror = 12 ! multigrid parameters
292 19 : if (kcycle.lt.0) return
293 19 : if (min0(iprer,ipost).lt.1) return
294 19 : if ((intpol-1)*(intpol-3).ne.0) return
295 19 : if (max0(kcycle,iprer,ipost).gt.2) then
296 19 : ierror = -5 ! inefficient multigrid cycling
297 : end if
298 19 : if (ierror .gt. 0) ierror = 0 ! no fatal errors
299 : !
300 : ! set work space pointers and discretize pde at each grid level
301 : !
302 19 : iw = 1
303 114 : do kb=1,ngrid
304 95 : k = ngrid-kb+1
305 95 : nx = nxk(k)
306 95 : ny = nyk(k)
307 95 : kpbgn(k) = iw
308 95 : kcbgn(k) = kpbgn(k)+(nx+2)*(ny+2)
309 95 : ktxbgn(k) = kcbgn(k)+10*nx*ny
310 95 : ktybgn(k) = ktxbgn(k)+isx*nx*ny
311 95 : iw = ktybgn(k)+jsy*nx*ny
312 95 : ic = kcbgn(k)
313 95 : itx = ktxbgn(k)
314 95 : ity = ktybgn(k)
315 95 : klevel = k
316 0 : call dismd2cr(nx,ny,work(ic),work(itx),work(ity), &
317 114 : work,ierror,isolve)
318 : end do
319 : return
320 : end if ! end of intl=0 initialization call block
321 19 : nx = nfx
322 19 : ny = nfy
323 19 : call mud2c1m(nx,ny,rhs,phi,work)
324 19 : iparm(17) = itero
325 19 : if (tolmax.gt.0.0_r8) then ! check for convergence
326 19 : fparm(6) = relmax
327 19 : if (relmax.gt.tolmax) then
328 :
329 : ! ierror = -1 ! flag convergenc failure
330 0 : write(iulog,*) "no convergence with mudmod"
331 : !
332 0 : iguess = 1
333 0 : iparm(12)= iguess
334 0 : call mud2cr1(nx,ny,rhs,phi,work) ! solve with modified stencils
335 :
336 0 : fparm(6) = relmax
337 0 : if (relmax.gt.tolmax) then
338 0 : write(iulog,*) "no convergence with mud"
339 0 : ierror = -1 ! flag convergenc failure
340 : end if
341 :
342 : end if
343 : end if
344 :
345 : return
346 80408 : end subroutine mud2cm
347 : !------------------------------------------------------------------------
348 38 : subroutine mud2c1m(nx,ny,rhsf,phif,wk)
349 :
350 : integer nx,ny
351 : real(r8) :: phif(nx,ny),rhsf(nx,ny),wk(*)
352 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
353 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
354 : kcycle,iprer,ipost,intpol,kps
355 : real(r8) :: xa,xb,yc,yd,tolmax,relmax,phmax
356 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
357 : integer k,kb,ip,ic,ir,ipc,irc,icc
358 : integer ncx,ncy,jj,ij,i,j,iter
359 :
360 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
361 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
362 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
363 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
364 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), &
365 : nxk(50),nyk(50),isx,jsy
366 :
367 19 : nx = nxk(ngrid)
368 19 : ny = nyk(ngrid)
369 19 : ip = kpbgn(ngrid)
370 19 : ic = kcbgn(ngrid)
371 19 : ir = ic+9*nx*ny
372 : !
373 : ! set phif,rhsf in wk and adjust right hand side
374 : !
375 19 : call swk2(nx,ny,phif,rhsf,wk(ip),wk(ir))
376 19 : if (iguess.eq.0) then
377 : !
378 : ! no initial guess at finest grid level!
379 : !
380 95 : do kb=2,ngrid
381 76 : k = ngrid-kb+1
382 76 : nx = nxk(k+1)
383 76 : ny = nyk(k+1)
384 76 : ip = kpbgn(k+1)
385 76 : ir = kcbgn(k+1)+9*nx*ny
386 76 : ncx = nxk(k)
387 76 : ncy = nyk(k)
388 76 : ipc = kpbgn(k)
389 76 : icc = kcbgn(k)
390 76 : irc = icc+9*ncx*ncy
391 : !
392 : ! transfer down to all grid levels
393 : !
394 95 : call trsfc2(nx,ny,wk(ip),wk(ir),ncx,ncy,wk(ipc),wk(irc))
395 : end do
396 : !
397 : ! adjust right hand side at all grid levels in case
398 : ! rhs or specified b.c. in phi or gbdy changed
399 : !
400 114 : do k=1,ngrid
401 95 : nx = nxk(k)
402 95 : ny = nyk(k)
403 95 : ip = kpbgn(k)
404 95 : ic = kcbgn(k)
405 114 : call adjmd2cr(nx,ny,wk(ip),wk(ic))
406 : end do
407 : !
408 : ! execute one full multigrid cycle
409 : !
410 95 : do k=1,ngrid-1
411 76 : kcur = k
412 76 : call kcymd2cr(wk)
413 76 : nx = nxk(k+1)
414 76 : ny = nyk(k+1)
415 76 : ip = kpbgn(k+1)
416 76 : ipc = kpbgn(k)
417 76 : ncx = nxk(k)
418 76 : ncy = nyk(k)
419 : !
420 : ! lift or prolong approximation from k to k+1
421 : !
422 95 : call prolon2(ncx,ncy,wk(ipc),nx,ny,wk(ip),nxa,nxb,nyc,nyd,intpol)
423 : end do
424 : else
425 : !
426 : ! adjust rhs at finest grid level only
427 : !
428 0 : nx = nxk(ngrid)
429 0 : ny = nyk(ngrid)
430 0 : ip = kpbgn(ngrid)
431 0 : ic = kcbgn(ngrid)
432 0 : call adjmd2cr(nx,ny,wk(ip),wk(ic))
433 : end if
434 : !
435 : ! execute maxcy more multigrid k cycles from finest level
436 : !
437 19 : kcur = ngrid
438 76 : do iter=1,maxcy
439 76 : itero = iter
440 76 : call kcym2cm(wk)
441 0 : if (tolmax.gt.0.0_r8) then
442 : !
443 : ! error control
444 : !
445 76 : relmax = 0.0_r8
446 76 : phmax = 0.0_r8
447 :
448 3800 : do j=1,nfy
449 3724 : jj = j*(nfx+2)
450 305444 : do i=1,nfx
451 301644 : ij = jj+i+1
452 301644 : phmax = max(phmax,abs(wk(ij)))
453 301644 : relmax = max(relmax,abs(wk(ij)-phif(i,j)))
454 :
455 305368 : phif(i,j) = wk(ij)
456 : end do
457 : end do
458 : !
459 : ! set maximum relative difference and check for convergence
460 : !
461 76 : if (phmax.gt.0.0_r8) relmax = relmax/phmax
462 76 : if (relmax.le.tolmax) return
463 : end if
464 : end do
465 : !
466 : ! set final interate after maxcy cycles in phif
467 : !
468 0 : do j=1,nfy
469 0 : jj = j*(nfx+2)
470 0 : do i=1,nfx
471 0 : ij = jj+i+1
472 0 : phif(i,j) = wk(ij)
473 : end do
474 : end do
475 : return
476 646 : end subroutine mud2c1m
477 :
478 : !------------------------------------------------------------------------
479 152 : subroutine kcym2cm(wk)
480 4066 : use edyn_solver_coefs,only: cofum
481 : !
482 : ! execute multigrid k cycle from kcur grid level
483 : ! kcycle=1 for v cycles, kcycle=2 for w cycles
484 : !
485 : real(r8) :: wk(*)
486 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
487 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
488 : kcycle,iprer,ipost,intpol,kps
489 : integer nx,ny,ip,ic,ipc,irc,itx,ity,ncx,ncy,l,nrel
490 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
491 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
492 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
493 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
494 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
495 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
496 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), &
497 : nxk(50),nyk(50),isx,jsy
498 : integer kount(50)
499 : ! real(r8) :: :: cofum
500 : ! common/mudmd/cofum(1)
501 :
502 76 : klevel = kcur
503 76 : nx = nxk(klevel)
504 76 : ny = nyk(klevel)
505 76 : ip = kpbgn(klevel)
506 76 : ic = kcbgn(klevel)
507 76 : itx = ktxbgn(klevel)
508 76 : ity = ktybgn(klevel)
509 : !
510 : ! prerelax at current finest grid level
511 : !
512 304 : do l=1,iprer
513 532 : call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
514 : end do
515 76 : if (kcur .eq. 1) go to 5
516 : !
517 : ! restrict residual to kcur-1 level
518 : !
519 76 : ipc = kpbgn(klevel-1)
520 76 : ncx = nxk(klevel-1)
521 76 : ncy = nyk(klevel-1)
522 76 : irc = kcbgn(klevel-1)+9*ncx*ncy
523 : ! call resmd2cr(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic),wk(kps))
524 :
525 0 : call resm2cm(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic), &
526 152 : wk(kps),cofum)
527 : !
528 : ! set counter for grid levels to zero
529 : !
530 456 : do l = 1,kcur
531 456 : kount(l) = 0
532 : end do
533 : !
534 : ! set new grid level and continue k-cycling
535 : !
536 76 : klevel = kcur-1
537 76 : nrel = iprer
538 : !
539 : ! kcycle control point
540 : !
541 : 10 continue
542 : !
543 : ! post relax when kcur revisited
544 : !
545 2280 : if (klevel .eq. kcur) go to 5
546 : !
547 : ! count hit at current level
548 : !
549 2204 : kount(klevel) = kount(klevel)+1
550 : !
551 : ! relax at current level
552 : !
553 2204 : nx = nxk(klevel)
554 2204 : ny = nyk(klevel)
555 2204 : ip = kpbgn(klevel)
556 2204 : ic = kcbgn(klevel)
557 2204 : itx = ktxbgn(klevel)
558 2204 : ity = ktybgn(klevel)
559 7752 : do l=1,nrel
560 13300 : call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
561 : end do
562 2204 : if (kount(klevel) .eq. kcycle+1) then
563 : !
564 : ! kcycle complete at klevel
565 : !
566 684 : ipc = ip
567 684 : ip = kpbgn(klevel+1)
568 684 : ncx = nxk(klevel)
569 684 : ncy = nyk(klevel)
570 684 : nx = nxk(klevel+1)
571 684 : ny = nyk(klevel+1)
572 : !
573 : ! inject correction to finer grid
574 : !
575 0 : call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd, &
576 1368 : intpol,wk(kps))
577 : !
578 : ! reset counter to zero
579 : !
580 684 : kount(klevel) = 0
581 : !
582 : ! ascend to next higher level and set to postrelax there
583 : !
584 684 : klevel = klevel+1
585 684 : nrel = ipost
586 684 : go to 10
587 : else
588 1520 : if (klevel .gt. 1) then
589 : !
590 : ! kcycle not complete so descend unless at coarsest grid
591 : !
592 1064 : ipc = kpbgn(klevel-1)
593 1064 : ncx = nxk(klevel-1)
594 1064 : ncy = nyk(klevel-1)
595 1064 : irc = kcbgn(klevel-1)+9*ncx*ncy
596 0 : call resmd2cr(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic), &
597 2128 : wk(kps))
598 : !
599 : ! prerelax at next coarser level
600 : !
601 1064 : klevel = klevel-1
602 1064 : nrel = iprer
603 1064 : go to 10
604 : else
605 : !
606 : ! postrelax at coarsest level
607 : !
608 1368 : do l=1,ipost
609 2280 : call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
610 : end do
611 456 : ipc = ip
612 456 : ip = kpbgn(2)
613 456 : ncx = nxk(1)
614 456 : ncy = nyk(1)
615 456 : nx = nxk(2)
616 456 : ny = nyk(2)
617 : !
618 : ! inject correction to level 2
619 : !
620 0 : call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd, &
621 912 : intpol,wk(kps))
622 : !
623 : ! set to postrelax at level 2
624 : !
625 456 : nrel = ipost
626 456 : klevel = 2
627 532 : go to 10
628 : end if
629 : end if
630 : 5 continue
631 : !
632 : ! post relax at current finest grid level
633 : !
634 76 : nx = nxk(kcur)
635 76 : ny = nyk(kcur)
636 76 : ip = kpbgn(kcur)
637 76 : ic = kcbgn(kcur)
638 76 : itx = ktxbgn(kcur)
639 76 : ity = ktybgn(kcur)
640 228 : do l=1,ipost
641 380 : call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
642 : end do
643 76 : return
644 : end subroutine kcym2cm
645 : !----------------------------------------------------------------------
646 76 : subroutine resm2cm(nx,ny,phi,ncx,ncy,phic,rhsc,cof,resf,cofum)
647 : !
648 : ! restrict residual from fine to coarse mesh using fully weighted
649 : ! residual restriction
650 : !
651 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
652 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
653 : kcycle,iprer,ipost,intpol,kps
654 : integer nx,ny,ncx,ncy,i,j,ic,jc
655 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
656 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
657 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
658 : real(r8) :: rhsc(ncx,ncy),resf(nx,ny)
659 : real(r8) :: phi(0:nx+1,0:ny+1),phic(0:ncx+1,0:ncy+1)
660 : real(r8) :: cof(nx,ny,10),cofum(nx,ny,9)
661 : real(r8) :: l2norm
662 : !
663 : ! set phic zero
664 : !
665 2128 : do jc=0,ncy+1
666 90364 : do ic=0,ncx+1
667 90288 : phic(ic,jc) = 0.0_r8
668 : end do
669 : end do
670 :
671 76 : call bnd2cm(nx,ny,cofum)
672 : !
673 : ! compute residual on fine mesh in resf
674 : !
675 76 : l2norm = 0._r8
676 : !$OMP PARALLEL DO SHARED(resf,cof,phi,nx,ny) PRIVATE(i,j)
677 3800 : do j=1,ny
678 305444 : do i=1,nx
679 603288 : resf(i,j) = cof(i,j,10)-( &
680 603288 : cofum(i,j,1)*phi(i+1,j)+ &
681 301644 : cofum(i,j,2)*phi(i+1,j+1)+ &
682 301644 : cofum(i,j,3)*phi(i,j+1)+ &
683 301644 : cofum(i,j,4)*phi(i-1,j+1)+ &
684 : cofum(i,j,5)*phi(i-1,j)+ &
685 301644 : cofum(i,j,6)*phi(i-1,j-1)+ &
686 : cofum(i,j,7)*phi(i,j-1)+ &
687 : cofum(i,j,8)*phi(i+1,j-1)+ &
688 1809864 : cofum(i,j,9)*phi(i,j))
689 :
690 305368 : l2norm = l2norm + resf(i,j)*resf(i,j)
691 : end do
692 : end do
693 : !
694 : ! restrict resf to coarse mesh in rhsc
695 : !
696 76 : call res2(nx,ny,resf,ncx,ncy,rhsc,nxa,nxb,nyc,nyd)
697 76 : return
698 17784 : end subroutine resm2cm
699 :
700 : !-----------------------------------------------------------------------
701 76 : subroutine bnd2cm(nx,ny,cf)
702 : !
703 : ! set stencil & boundary condition for finest stencil
704 : !
705 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
706 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
707 : kcycle,iprer,ipost,intpol,kps
708 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
709 : integer nx,ny,i,j,l
710 : real(r8) :: cf(nx,ny,*)
711 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
712 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
713 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
714 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
715 :
716 : !
717 : ! set coefficient for specified boundaries
718 : !
719 76 : if (nxa.eq.1) then
720 0 : i = 1
721 0 : do j=1,ny
722 0 : do l=1,9
723 0 : cf(i,j,l) = 0.0_r8
724 : end do
725 0 : cf(i,j,9) = 1.0_r8
726 : end do
727 : end if
728 76 : if (nxb.eq.1) then
729 0 : i = nx
730 0 : do j=1,ny
731 0 : do l=1,9
732 0 : cf(i,j,l) = 0.0_r8
733 : end do
734 0 : cf(i,j,9) = 1.0_r8
735 : end do
736 : end if
737 76 : if (nyc.eq.1) then
738 0 : j = 1
739 0 : do i=1,nx
740 0 : do l=1,9
741 0 : cf(i,j,l) = 0.0_r8
742 : end do
743 0 : cf(i,j,9) = 1.0_r8
744 : end do
745 : end if
746 76 : if (nyd.eq.1) then
747 76 : j = ny
748 6232 : do i=1,nx
749 61560 : do l=1,9
750 61560 : cf(i,j,l) = 0.0_r8
751 : end do
752 6232 : cf(i,j,9) = 1.0_r8
753 : end do
754 : end if
755 : !
756 76 : return
757 : end subroutine bnd2cm
758 : !-----------------------------------------------------------------------
759 : end module edyn_mudmod
|
