c----------------------------------------------------------- c 'mgfas2d' and supporting routines constitute a complete c FAS multi-grid solver for c c u_xx + u_yy + G(u) = rhs(x,y) (1) c c on the unit square with homogenous Dirichlet boundary c conditions. Literature (particularly Brandt) should be c consulted for treatment of Neumann and mixed boundary c conditions. c c Currently program solves (1) with c c G(u) = sigma * u*u c c and 'rhs' specified so that c c u = sin(Pi lx x) sin(Pi ly y) c c for some integers, 'lx' and 'ly' defined in c 'mgprob.inc'. c----------------------------------------------------------- program mgfas2d implicit none c----------------------------------------------------------- c Functions called c----------------------------------------------------------- integer iargc, i4arg real*8 r8arg, dmnrm2 c----------------------------------------------------------- c Command-line arguments c----------------------------------------------------------- integer level, nvcycle, preswp, pstswp c----------------------------------------------------------- c 'mg.inc' contains the main 'data structure' used c in this implementation. c----------------------------------------------------------- include 'mg.inc' c----------------------------------------------------------- c 'mgprob.inc' contains problem-specific parameters c including command-line arguments 'sigma' and c 'gft_trace'. c----------------------------------------------------------- include 'mgprob.inc' c----------------------------------------------------------- c Locals. c----------------------------------------------------------- integer l, cycle c----------------------------------------------------------- c Argument parsing. c c level: Finest level solution requested. Full MG c process solves sequence of problems c l = 1 ... level (n_l = 2^l+1 pts per side) c sigma: Controls nonlinearity in PDE (1) c nvcycle: # of vcycles to apply to each problem c preswp: # of relaxation sweeps before a CGC c pstswp: # of relaxation sweeps after a CGC c gft_trace: If present, enables RNPL-style .hdf c output of solution estimate 'u' after c each relaxation sweep. c----------------------------------------------------------- if( iargc() .lt. 2 ) go to 900 level = i4arg(1,-1) if( level .le. 0 ) go to 900 if( level .gt. lmax ) then write(0,*) 'mgfas2d: Specified level is too high' write(0,*) 'mgfas2d: Current maximum level: ', lmax stop end if sigma = r8arg(2,0.0d0) nvcycle = i4arg(3,1) preswp = i4arg(4,1) pstswp = i4arg(5,1) gft_trace = iargc() .gt. 5 if( gft_trace ) write(0,*) 'gft tracing enabled' c----------------------------------------------------------- c Initialize data structure. All grid function arrays c are allocated from global array 'q' (see 'mg.inc') c----------------------------------------------------------- call gsinit(level) c----------------------------------------------------------- c Solve series of problems from l = 1 through c l = level. c----------------------------------------------------------- do l = 1 , level if( l .eq. 1 ) then c----------------------------------------------------------- c Initialize level 1 solution (coarsest) with 0. c----------------------------------------------------------- call dmls ( q(u(l)), 0.0d0, nx(l), ny(l) ) else c----------------------------------------------------------- c Initialize level l soln via linear interpolation c of level l-1 soln. c----------------------------------------------------------- call dmlint ( q(u(l-1)), q(u(l)), & nx(l-1), ny(l-1) ) end if c----------------------------------------------------------- c Initialize 'rhs' (rhs = f). c----------------------------------------------------------- call clcrhs2 ( q(rhs(l)), h(l), nx(l), ny(l) ) c----------------------------------------------------------- c "Solve" level l problem by applying cycle c (preswp,pstswp) V-cycles. c----------------------------------------------------------- do cycle = 1 , nvcycle call vcycle(l,cycle,preswp,pstswp) end do end do c----------------------------------------------------------- c Compute norm of error on finest gird and output to c standard error. c----------------------------------------------------------- l = level call clcu2 ( q(t1(l)), h(l), nx(l), ny(l) ) call dmms ( q(t1(l)), q(u(l)), q(t1(l)), & nx(l), ny(l) ) write(0,1000) dmnrm2( q(t1(l)), nx(l), ny(l) ) 1000 format(' || Error ||: ',1p,e12.4) c----------------------------------------------------------- c Close RNPL-style .hdf files if option enabled. c----------------------------------------------------------- if( gft_trace ) then call gft_close_all() end if c----------------------------------------------------------- c Normal exit. c----------------------------------------------------------- stop c----------------------------------------------------------- c Usage exit. c----------------------------------------------------------- 900 continue write(0,*) 'usage: mgfas2d '// & '[ ]' stop end c----------------------------------------------------------- c Performs one FAS V-cycle of level 'l' problem, using c 'preswp' pre-CGC and 'pstswp' pst-CGC sweeps. c Supply cycle=N for Nth consecutive invocation of this c routine on the same level. c----------------------------------------------------------- subroutine vcycle(l,cycle,preswp,pstswp) implicit none include 'mg.inc' integer l, cycle, preswp, pstswp real*8 dmnrm2, relaxc2 integer m, sweep real*8 rresnrm c----------------------------------------------------------- c Work from fine level to coarse level, performing c pre-CGC sweeps and intitiating coarse grid c corrections. Skipped if on coarsest level, l = 1. c----------------------------------------------------------- do m = l , 2 , -1 c----------------------------------------------------------- c Pre-CGC relaxation sweeps. c----------------------------------------------------------- if( cycle .eq. 1 .or. m .ne. l ) then do sweep = 1 , preswp rresnrm = relaxc2(q(u(m)),q(rhs(m)), & nx(m),ny(m),h(m)) call my_gft_out(m) call trace(m,'Pre-CGC',rresnrm) end do end if c----------------------------------------------------------- c Compute truncation error estimate. c c 2h 2h 2h h 2h h h c tau := L I u - I L u c h h h c c As pointed out in 'mg.inc', this implementation c uses three temporary arrays for clarity (:-)); c production implementation should economize. c----------------------------------------------------------- c----------------------------------------------------------- c 2h 2h 2h h c t2 := L I u c h c----------------------------------------------------------- call dmrshw ( q(u(m)), q(t1(m-1)), & nx(m-1), ny(m-1) ) call lop ( q(t2(m-1)), q(t1(m-1)), & nx(m-1), ny(m-1), h(m-1) ) c----------------------------------------------------------- c 2h 2h h h c t3 := I L u c h c----------------------------------------------------------- call lop ( q(t1(m)), q(u(m)), & nx(m), ny(m), h(m) ) call dmrshw ( q(t1(m)), q(t3(m-1)), & nx(m-1),ny(m-1) ) c----------------------------------------------------------- c 2h 2h 2h c tau := t2 - t3 c h c----------------------------------------------------------- call dmms ( q(t2(m-1)), q(t3(m-1)), q(tau(m-1)), & nx(m-1), ny(m-1) ) c----------------------------------------------------------- c Define right hand side and initialize unknown c for coarse grid problem. c c 2h 2h h 2h c rhs := I rhs + tau c h h c c 2h 2h h c u := I u c h c----------------------------------------------------------- call dmrshw ( q(rhs(m)), q(rhs(m-1)), & nx(m-1),ny(m-1)) call dmma ( q(rhs(m-1)), q(tau(m-1)), q(rhs(m-1)), & nx(m-1), ny(m-1) ) call dmrshw ( q(u(m)), q(u(m-1)), & nx(m-1), ny(m-1) ) c----------------------------------------------------------- c Continue loop to perform coarse-grid c correction. c----------------------------------------------------------- end do c----------------------------------------------------------- c Solve coarse-grid problem using relaxation. c Note use of 'do while()' c----------------------------------------------------------- rresnrm = 1.0d0 + epsi1 sweep = 1 do while( rresnrm .gt. epsi1 .and. & sweep .le. maxsweep1 ) rresnrm = relaxc2( q(u(1)), q(rhs(1)), & nx(1), ny(1), h(1) ) sweep = sweep + 1 call my_gft_out(1) call trace(1,'Coarse grid',rresnrm) end do if( rresnrm .gt. epsi1 ) then write(0,*) '*** WARNING vcycle: Coarse grid '// & 'relaxation failed.' end if c----------------------------------------------------------- c Now work from coarse to fine level, updating c unknowns and performing post-CGC relaxtion c sweeps. Skipped if l = 1. c----------------------------------------------------------- do m = 2 , l c----------------------------------------------------------- c Linearly interpolate correction and update c level m unknown. c c h h h 2h 2h h c u := u + I ( u - I u ) c 2h h c----------------------------------------------------------- call dmrshw ( q(u(m)), q(t1(m-1)), & nx(m-1), ny(m-1) ) call dmms ( q(u(m-1)), q(t1(m-1)), q(t1(m-1)), & nx(m-1), ny(m-1) ) call dmlint ( q(t1(m-1)), q(t2(m)), & nx(m-1),ny(m-1) ) call dmma ( q(u(m)), q(t2(m)), q(u(m)), & nx(m),ny(m) ) c----------------------------------------------------------- c Post-CGC relaxation sweeps. c----------------------------------------------------------- do sweep = 1 , pstswp rresnrm = relaxc2 ( q(u(m)), q(rhs(m)), & nx(m),ny(m),h(m) ) call my_gft_out(m) call trace(m,'Post-CGC',rresnrm) end do end do c----------------------------------------------------------- c End of the V-cycle. c----------------------------------------------------------- write(0,*) return end