!----------------------------------------------------------------------- ! ! (c) 2000-2001 Dale Choi ! ! HISTORY: mg_main.F (J. Gao's Unigrid MG solver) ! ! Main Driver of MultiGrid Elliptic Solver ! this code solves a poisson-type elliptic equation using LCS (Linear ! Correction Scheme), with octant symmentry and robin boundary condition ! ! (1) Tested in parADM Version 1 ! ! Input GF (gup,gam,Acoeff,\rho) and Output GF (\psi) defined at node-centers ! ! (2) Further modification to solve: ! ! gup11 * u,xx + gup22 * u,yy + gup33 * u,zz + acoeff * u^2 = rhs ! ! (3) Modification (being made) for AMR (non-uniform) grids ! ! Mar 14 2001 DONE: [1] reworked mg_gsrb_level.F / mg_gsrb_relax.F ! [2] OBC (mg_relaxb,mg_octant,mg_relaxb_robin) ! [3] check mg_lop.F (flux_matching) ! ! Thu Sep 6 2001 Stand Alone version (LCS) created ! Thu Sep 27 2001 2nd order convergence against analytic data verified! ! ! Sat Sep 22 16:08:38 EDT 2001 ! FAS (Full Approximation Storage scheme) -- implemented & being tested now ! !----------------------------------------------------------------------- program mg_main implicit none #include "physicaldata.fh" ! Use this for PARAMESH V.1 ! include 'mpp/shmem.fh' ! Use this for PARAMESH V.2 #include "amr_shmem.fh" include 'tree.fh' include 'shmem_reduce.fh' include 'mg_data.fh' include 'mg_node.fh' include 'mg_oct.fh' integer mype, i, ir, iv, level, ix,iy,iz,lb, nprocs, maxprocs integer nrelax_up integer shmem_my_pe, shmem_n_pes save mype, nprocs !----------------------------------------------------------------------- ! initial paramesh !----------------------------------------------------------------------- ! call comm_start(maxprocs,nprocs,mype) mype = shmem_my_pe() nprocs= shmem_n_pes() call amr_initialize call shmem_barrier_all() !----------------------------------------------------------------------- ! read-in from parameter file and set up initial refinement mesh !----------------------------------------------------------------------- call mg_init_param(mype) call shmem_barrier_all() call mg_init_grid(mype) call shmem_barrier_all() !----------------------------------------------------------------------- ! set up initial data ! set RHS, RES here! (calls mg_lop_nf) !----------------------------------------------------------------------- call mg_init_soln(lrefine_max) call shmem_barrier_all() !----------------------------------------------------------------------- ! Special treatment for LMAX = 1 (It *is* a simple relaxation!) !----------------------------------------------------------------------- if( lrefine_max .eq. 1 ) then ! redistribute blocks among processors call amr_morton_order(nprocs,lnblocks) write(0,*)'mype,lnblocks = ',mype,lnblocks if( mg_type .eq. 0 .or. mg_type .eq. 2 ) then do i = 1, ns do ir = 0, 1 call shmem_barrier_all() call mg_guardcell(iu,1) call shmem_barrier_all() !sep26 2001 #if OCTANT_TEST == 0 call mg_octant(iu,1) #else if OCTANT_TEST == 1 call mg_octant2(iu,1) #endif call shmem_barrier_all() !sep26 2001 call mg_relax_swid(ir,1) !swid end do do ir = 0, 1 call shmem_barrier_all() call mg_guardcell(iu,1) call shmem_barrier_all() !sep26 2001 #if OCTANT_TEST == 0 call mg_octant(iu,1) #else if OCTANT_TEST == 1 call mg_octant2(iu,1) #endif !call mg_relaxb(ir,1) ! original (DONOT USE!) call mg_fillc(itmp,1.0e0,1) call shmem_barrier_all() !sep26 2001 call mg_relaxb_robin(iu,itmp,ir,1) !to control 'rhs' for robin boundary eqn. end do call shmem_barrier_all() !call mg_comp_res(itmp,1,i,0) ! use itmp as temp. storage call mg_comp_res(itmp,1,-100,i) end do else if( mg_type .eq. 1 ) then do i=1,ns do ir=0,1 call shmem_barrier_all() call mg_guardcell1(iu,1) ! apr 3 call mg_octant(iu,1) call mg_lapl_0(iu,irhs,1,ir) ! simple relaxation end do do ir=0,1 call shmem_barrier_all() call mg_guardcell1(iu,1) ! apr 3 call mg_octant(iu,1) call mg_relaxb_robin(iu,irhs,ir,1) end do call shmem_barrier_all() ! Mar 20, 2001 call mg_comp_res(itmp,1,i,0) ! use itmp as temp. storage enddo end if call mg_dump_1lev(iu,1,'u_') call shmem_barrier_all() stop endif !----------------------------------------------------------------------- ! ! V-Cycle nvcycle := # of V-Cycle ! !----------------------------------------------------------------------- do iv = 1, nvcycle if( mype .eq. 0 ) write(0,*)'V-cycle #= ',iv if( mype .eq. 0 ) write(59,*)'V-cycle #= ',iv !---------------------------------------------------------------- ! Down leg of V-cycle !---------------------------------------------------------------- do level = lrefine_max, 2, -1 if( mype .eq. 0 ) write(0,*)' Down Leg: level = ',level if( mype .eq. 0 ) write(59,*)' Down Leg: level = ',level !-------------------------------------------------- ! Relaxation on 'level' nrelax := # of relaxation !-------------------------------------------------- if( mg_type .eq. 0 ) then do i = 1, nrelax do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,level) call shmem_barrier_all() !sep26 2001 #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif call shmem_barrier_all() !sep26 2001 call mg_relax_swid(ir,level) end do do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,level) #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif !call mg_relaxb(ir,level) !original (DONOT USE!) call mg_fill0(itmp,level) !apr 21 2001 if(level.eq.lrefine_max) & call mg_fillc(itmp,1.0e0,level) !apr 21 2001 call shmem_barrier_all() !sep26 2001 call mg_relaxb_robin(iu,itmp,ir,level) !apr 21 2001 end do call shmem_barrier_all() !sep26 2001 call mg_comp_res(itmp,level,0,iv) call shmem_barrier_all() !sep26 2001 end do else if( mg_type .eq. 1 ) then call shmem_barrier_all() ! Feb 27 2001 call mg_copy(iusave,iu,level) call mg_guardcell1(iusave,level) ! Apr 10 2001 call mg_octant(iusave,level) ! Mar 14 2001 call set_interface_from_c(iusave,level) !Apr 2 2001 call mg_fill0(ierr,level-1) call mg_guardcell1(ires,level) ! Apr 10 2001 call mg_octant(ires,level) ! Mar 5 2001 call set_interface_from_c(ires,level) !Apr 2 2001 call shmem_barrier_all() ! Mar 14 2001 call mg_gsrb_level(ierr,level) call shmem_barrier_all() ! Mar 14 2001 call mg_guardcell1(ierr,level) ! Apr 10 2001 call mg_octant(ierr,level) ! Mar 5 2001 call set_interface_from_c(ierr,level) !Apr 2 2001 call shmem_barrier_all() ! Mar 14 2001 call mg_add(iu,iu,ierr,level) !ok call shmem_barrier_all() ! Mar 14 2001 call mg_guardcell1(iu,level) ! Apr 10 2001 call mg_octant(iu,level) ! Mar 5 2001 call set_interface_from_c(iu,level) !Apr 2 2001 else if( mg_type .eq. 2 ) then do i = 1, nrelax do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,level) #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif call mg_relax_swid(ir,level) end do do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,level) #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif call mg_fillc(itmp,1.0e0,level) call mg_relaxb_robin(iu,itmp,ir,level) end do call mg_comp_res(itmp,level,0,iv) end do end if call shmem_barrier_all() !-------------------------------------------------- ! Set up the coarse grid problem on 'level-1' !-------------------------------------------------- if( mg_type .eq. 0 ) then call mg_guardcell(iu,level) call shmem_barrier_all() #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif call shmem_barrier_all() !sep26 2001 call mg_residue_swid(level) ! Modified for swid(011901) & AMR(062101) !check mg_residue_swid again (might want to use ivar to compue residue?!?! ! call set_interface_from_c(ires,level) !Jul 9, 2001 call shmem_barrier_all() call mg_guardcell(ires,level) #if OCTANT_TEST == 0 call mg_octant(ires,level) #else if OCTANT_TEST == 1 call mg_octant2(ires,level) #endif call shmem_barrier_all() call mg_restrict(ires,level-1) call shmem_barrier_all() ! sep26 2001 call mg_guardcell(ires,level-1) !? need guardcell1 call for level-1 data? (Jun 21 2001) ! !--------------------------------------------------- ! ! trying to restrict coefficients ! !--------------------------------------------------- ! call mg_restrict(iacoeff,level-1) ! added on 020101 ! call mg_restrict(igup11c,level-1) ! added on 020101 ! call mg_restrict(igup22c,level-1) ! added on 020101 ! call mg_restrict(igup33c,level-1) ! added on 020101 ! ok, but restriction is done only for non-guardcells (Jun 19, 2001) ! also over-write some data that causes problems! ! ! Apr 19 2001 guardcell call really helps! ! call mg_guardcell(iacoeff,level-1) !no '1's here (jun 20, 2001) ! call mg_guardcell(igup11c,level-1) ! call mg_guardcell(igup22c,level-1) ! call mg_guardcell(igup33c,level-1) #if OCTANT_TEST == 0 ! call mg_octant(iacoeff,level-1) ! call mg_octant(igup11c,level-1) ! call mg_octant(igup22c,level-1) ! call mg_octant(igup33c,level-1) #else if OCTANT_TEST == 1 ! call mg_octant2(iacoeff,level-1) ! call mg_octant2(igup11c,level-1) ! call mg_octant2(igup22c,level-1) ! call mg_octant2(igup33c,level-1) #endif call mg_copy(irhs,ires,level-1) call mg_fill0(iu,level-1) else if( mg_type .eq. 1 ) then call mg_residual(level-1) ! to test (mg_lop) only for AMR call shmem_barrier_all() call mg_restrict(iacoeff,level-1) ! 022201 call mg_restrict(igup11c,level-1) ! 022201 call mg_restrict(igup22c,level-1) ! 022201 call mg_restrict(igup33c,level-1) ! 022201 call mg_restrict(irhs,level-1) ! 022201 call shmem_barrier_all() call mg_guardcell1(iacoeff,level-1) ! Mar 19 2001 call mg_guardcell1(igup11c,level-1) call mg_guardcell1(igup22c,level-1) call mg_guardcell1(igup33c,level-1) call mg_guardcell1(irhs,level-1) call mg_octant(iacoeff,level-1) ! Mar 5 2001 call mg_octant(igup11c,level-1) ! Mar 5 2001 call mg_octant(igup22c,level-1) ! Mar 5 2001 call mg_octant(igup33c,level-1) ! Mar 5 2001 call mg_octant(irhs,level-1) ! Mar 5 2001 else if( mg_type .eq. 2 ) then !----------------------------------------------------------- ! Compute truncation error estimate. ! ! 2h 2h 2h h 2h h h ! tau := L I u - I L u ! h h h ! ! this implementation uses three temporary arrays for clarity (:-)); ! production implementation should economize. !----------------------------------------------------------- ! restriction to compute I^{2h}_{h} u^{h} part ! u(l-1) contains I^{2h}_{h} u^{h} call shmem_barrier_all() ! call mg_restrict(iu,level-1) ! dale's this gives increase of residual vs. # of relax call inject(mype,iu,level-1) ! jimin's inject gives decrease of resi vs. # of relax call shmem_barrier_all() call mg_guardcell(iu,level-1) call mg_octant(iu,level-1) ! compute L^{2h} (I^{2h}_{h} u^{h}) from u(l-1) = I^{2h}_{h} u^{h}, ! and store it in tmp(l-1) call shmem_barrier_all() call mg_fill0(itmp,level-1) call mg_lop2(iu,level-1) call shmem_barrier_all() call mg_guardcell(itmp,level-1) call mg_octant(itmp,level-1) ! compute L^{h} u^{h} part and store in tmp(l) call shmem_barrier_all() call mg_fill0(itmp,level) call mg_lop2(iu,level) call shmem_barrier_all() call mg_guardcell(itmp,level) call mg_octant(itmp,level) ! call mg_fill_bnd(itmp,level,0.0e0) ! restriction of L^{h} u^{h} and store in tau(l-1)= I^{2h}_{h}(L^{h} u^{h}) call shmem_barrier_all() call mg_fill0(itau,level-1) call mg_restrict2(itau,itmp,level-1) call shmem_barrier_all() call mg_guardcell(itau,level-1) call mg_octant(itau,level-1) ! compute truncation error and store in q(tau(l-1)) call shmem_barrier_all() call mg_sub(itau,itmp,itau,level-1) call shmem_barrier_all() call mg_guardcell(itau,level-1) call mg_octant(itau,level-1) ! at this point tau contains truncation error and tmp L^{h} u^{h} ! call mg_fill_bnd(itau,level-1,0.0e0) call mg_anorm(mype,itau,level-1) call mg_dump_1lev(itau,level-1,'lte_') !----------------------------------------------------------- ! Define right hand side and initialize unknown ! for coarse grid problem. ! ! 2h 2h h 2h ! rhs := I rhs + tau ! h h ! ! 2h 2h h ! u := I u ! h !----------------------------------------------------------- ! restrict rhs call shmem_barrier_all() call mg_restrict(irhs,level-1) call shmem_barrier_all() call mg_guardcell(irhs,level-1) call mg_octant(irhs,level-1) ! add truncation error to the rhs call shmem_barrier_all() call mg_add(irhs,irhs,itau,level-1) call shmem_barrier_all() call mg_guardcell(irhs,level-1) call mg_octant(irhs,level-1) ! call mg_fill_bnd(irhs,level-1,0.0e0) call mg_dump_1lev(irhs,level-1,'rhs_r') ! initial guess for the coarse grid problem ! restriction from the fine grid solution call shmem_barrier_all() call mg_restrict(iu,level-1) call shmem_barrier_all() call mg_guardcell(iu,level-1) call mg_octant(iu,level-1) call mg_copy(iuold,iu,level-1) ! jimin's endif call shmem_barrier_all() end do !---------------------------------------------------------------- ! solve on the coarest grid !---------------------------------------------------------------- if( mg_type .eq. 0 ) then do i=1,ns do ir=0,1 call shmem_barrier_all() call mg_guardcell1(iu,1) !apr 19 2001 #if OCTANT_TEST == 0 call mg_octant(iu,1) #else if OCTANT_TEST == 1 call mg_octant2(iu,1) #endif call mg_relax_swid(ir,1) !swid end do do ir=0,1 call shmem_barrier_all() call mg_guardcell1(iu,1) !apr 19 2001 #if OCTANT_TEST == 0 call mg_octant(iu,1) #else if OCTANT_TEST == 1 call mg_octant2(iu,1) #endif !call mg_relaxb(ir,1) !original (DONOT USE!) call mg_fill0(itmp,1) !apr 21 2001 call mg_relaxb_robin(iu,itmp,ir,1) !apr 21 2001 end do call shmem_barrier_all() call mg_comp_res(itmp,1,i,iv) ! use itmp as temp. storage end do else if( mg_type .eq. 1 ) then call mg_guardcell1(ires,1) ! apr 10 call mg_octant(ires,1) ! NEW Mar 14 call set_interface_from_c(ires,level) !Apr 2 2001 call shmem_barrier_all() ! NEW Mar 14 !call mg_dump_1lev(ierr,1,'error_') do i=1,ns do ir=0,1 call shmem_barrier_all() call mg_guardcell1(ierr,1) ! Apr 20 2001 call mg_octant(ierr,1) ! call set_interface_from_c(ierr,1) !Apr 10 2001 call mg_lapl_0(ierr,ires,1,ir) ! simple relaxation end do do ir=0,1 call shmem_barrier_all() call mg_guardcell1(ierr,1) ! Apr 10 2001 call mg_octant(ierr,1) ! call set_interface_from_c(ierr,1) !Apr 10 2001 call mg_relaxb_robin(ierr,ires,ir,1) ! NEW Feb 27 end do call shmem_barrier_all() call mg_add(iu,iu,ierr,1) ! what iu?? call shmem_barrier_all() call mg_comp_res(itmp,1,i,iv) !can use itmp ! call mg_anorm(mype,ierr,1) ! apr 11 2001 call shmem_barrier_all() call mg_sub(iu,iu,ierr,1) enddo call shmem_barrier_all() call mg_guardcell1(ierr,1) ! NEW Apr 10 call mg_octant(ierr,1) ! NEw Mar 14 call shmem_barrier_all() else if( mg_type .eq. 2 ) then do i= 1, ns do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,1) #if OCTANT_TEST == 0 call mg_octant(iu,1) #else if OCTANT_TEST == 1 call mg_octant2(iu,1) #endif call mg_relax_swid(ir,1) end do do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,1) #if OCTANT_TEST == 0 call mg_octant(iu,1) #else if OCTANT_TEST == 1 call mg_octant2(iu,1) #endif call mg_fillc(itmp,1.0e0,1) call mg_relaxb_robin(iu,itmp,ir,1) end do call shmem_barrier_all() call mg_guardcell(iu,1) ! NEW Oct 3 2001 call mg_octant(iu,1) ! NEw Oct 3 2001 call mg_comp_res(itmp,1,i,iv) ! use itmp as temp. storage end do end if call shmem_barrier_all() !---------------------------------------------------------------- ! Up leg of V-cycle !---------------------------------------------------------------- do level = 2, lrefine_max if( mype .eq. 0 ) write(0,*)' Up leg: level =', level if( mype .eq. 0 ) write(59,*)' Up leg: level =',level if( mg_type .eq. 0 ) then call shmem_barrier_all() !sep26 2001 call mg_copy(itmp,iu,level-1) call shmem_barrier_all() call mg_prolong(mype,itmp,level) !check this again! call shmem_barrier_all() !sep26 2001 call mg_add(iu,iu,itmp,level) ! call set_interface_from_c(iu,level) !Jul 9, 2001 call shmem_barrier_all() nrelax_up = nrelax if( i .eq. lrefine_max ) nrelax_up = nrelax_1 do i = 1, nrelax_up do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,level) #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif call mg_relax_swid(ir,level) !swid end do do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,level) #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif !call mg_relaxb(ir,level) !original (DONOT USE!) call mg_fill0(itmp,level) !apr 21 2001 if(level.eq.lrefine_max) & call mg_fillc(itmp,1.0e0,level) !apr 21 2001 call mg_relaxb_robin(iu,itmp,ir,level) !apr 21 2001 end do call mg_comp_res(itmp,level,0,iv) end do else if( mg_type .eq. 1 ) then call shmem_barrier_all() call mg_guardcell1(ierr,level-1) ! NEW Apr 10 call mg_octant(ierr,level-1) ! NEW Mar 14 call shmem_barrier_all() call mg_prolong1(mype,itmp,ierr,level) !itmp=interp(err{l-1}) call shmem_barrier_all() call mg_guardcell1(itmp,level) ! NEW Apr 10 call mg_octant(itmp,level) ! NEW Mar 5 call shmem_barrier_all() call mg_add(ierr,ierr,itmp,level) !ierr=ierr+itmp call shmem_barrier_all() call mg_guardcell1(ierr,level) ! NEW Apr 10 call mg_octant(ierr,level) ! NEW Mar 5 ! should we do smoother here for ierr?? apr 11 2001 call mg_guardcell1(ires,level) call mg_octant(ires,level) call mg_gsrb_level(ierr,level) ! apr 11 2001 call shmem_barrier_all() call mg_lop_nf(ierr,level) !itmp=lop_nf(ierr) call shmem_barrier_all() call mg_guardcell1(itmp,level) ! NEW Apr 10 call mg_octant(itmp,level) ! NEW Mar 14 call shmem_barrier_all() call mg_sub(ires,ires,itmp,level) !ires=ires-lop_nf call shmem_barrier_all() call mg_guardcell1(ires,level) ! NEW Apr 10 call mg_octant(ires,level) ! NEW Mar 5 call mg_fill0(iderr,level) ! call mg_fill0(iderr,level-1) ! NEW Apr 10 call shmem_barrier_all() call mg_gsrb_level(iderr,level) call shmem_barrier_all() call mg_guardcell1(iderr,level) ! apr 11 call mg_octant(iderr,level) ! Mar 14 call shmem_barrier_all() call mg_add(ierr,ierr,iderr,level) call shmem_barrier_all() call mg_guardcell1(ierr,level) ! apr 11 call mg_octant(ierr,level) ! Mar 14 call shmem_barrier_all() call mg_guardcell1(iusave,level) ! apr 11 call mg_octant(iusave,level) ! Mar 14 call shmem_barrier_all() call mg_add(iu,iusave,ierr,level) call shmem_barrier_all() call mg_guardcell1(iu,level) ! apr 11 call mg_octant(iu,level) ! NEW Mar 5 call mg_comp_res(itmp,level,0,iv) else if( mg_type .eq. 2 ) then 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----------------------------------------------------------- ! restrict previous fine grid solution and store in q(tau(l-1)) = I^{2h}_{h} u^{h} call shmem_barrier_all() call mg_fill0(itau,level-1) call mg_restrict2(itau,iu,level-1) call shmem_barrier_all() call mg_guardcell(itau,level-1) call mg_octant(itau,level-1) ! compute u^{2h} - I^{2h}_{h} u^{h} and store in q(tau(l-1)) call shmem_barrier_all() call mg_sub(itau,iu,itau,level-1) !dale's ! call mg_sub(itau,iu,iuold,level-1) !jimin's call shmem_barrier_all() call mg_guardcell(itau,level-1) call mg_octant(itau,level-1) ! interpolate to get correction to the previous solution; store in q(tmp(l)) call shmem_barrier_all() call mg_fill0(itmp,level) call mg_prolong1(mype,itmp,itau,level) call shmem_barrier_all() call mg_guardcell(itmp,level) call mg_octant(itmp,level) ! add the correction to the previous solution call shmem_barrier_all() call mg_add(iu,iu,itmp,level) call mg_guardcell(iu,level) call mg_octant(iu,level) ! relax to remove high freq. components do i = 1, nrelax do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,level) #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif call mg_relax_swid(ir,level) end do do ir=0,1 call shmem_barrier_all() call mg_guardcell(iu,level) #if OCTANT_TEST == 0 call mg_octant(iu,level) #else if OCTANT_TEST == 1 call mg_octant2(iu,level) #endif call mg_guardcell(iu,level) !oct 3 2001 call mg_fillc(itmp,1.0e0,level) call mg_relaxb_robin(iu,itmp,ir,level) call mg_guardcell(iu,level) !oct 3 2001 end do call mg_comp_res(itmp,level,0,iv) end do end if ! mg_type end do !up leg of V-cycle call shmem_barrier_all() !sep26 2001 call mg_comp_res(itmp,lrefine_max,-100,iv) end do !of V-Cycle call mg_dump_1lev(iu,lrefine_max,'u_') call mg_test(itmp,lrefine_max) call mg_dump_1lev(itmp,lrefine_max,'uan_') call mg_dump_1lev(irhs,lrefine_max,'rhs_') call mg_comp_admmass(iu,lrefine_max,0) call shmem_barrier_all() stop end