/*@@ @file sor.F @date Thu Mar 13 07:46:55 1997 @author Joan Masso + Paul Walker @desc The SOR solver @enddesc @version $Header$ @@*/ #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" #include "EllBase.h" /*@@ @routine sor @date Thu Apr 24 13:29:52 1997 @author Joan Masso @desc This is a standalone sor solver which does all the MPI itself. It is a pretty good example of doing pugh-based communications in FORTRAN. @enddesc @@*/ subroutine sor_flat_core3d( & ierr,_CCTK_FARGUMENTS, $ Mlinear_lsh,Mlinear, $ Nsource_lsh,Nsource, $ var,var_idx, $ abstol,reltol) implicit none _DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer ierr INTEGER Mlinear_lsh(3) CCTK_REAL Mlinear(Mlinear_lsh(1),Mlinear_lsh(2),Mlinear_lsh(3)) INTEGER Nsource_lsh(3) CCTK_REAL Nsource(Nsource_lsh(1),Nsource_lsh(2),Nsource_lsh(3)) CCTK_REAL var(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) INTEGER var_idx CCTK_REAL abstol(3),reltol(3) CCTK_REAL tol INTEGER toltype CCTK_REAL dx,dy,dz c Temporaries INTEGER sor_iteration CCTK_REAL tmp c Numbers... CCTK_REAL two, four c Total residual CCTK_REAL resnorm, residual c Stencil CCTK_REAL a(-1:1,-1:1,-1:1) c Loopers INTEGER i,j,k c Acceleration factor CCTK_REAL omega, rjacobian c conformal test logical conformal logical octant c Loop bounds. Starts, Ends, and deltas (steps) INTEGER is, js, ks, ie, je, ke, di, dj, dk, kstep c Coeeficients for the solver: 19 point stencil... CCTK_REAL ac,ac_orig,aw,ae,an,as,at,ab CCTK_REAL ane,anw,ase,asw,ate,atw CCTK_REAL abe,abw,atn,ats,abn,asb CCTK_REAL finf CCTK_INT npow logical cheb, const, none, verb integer Mlinear_storage,Nsource_storage INTEGER sum_handle c stencil size INTEGER sw(3) ierr = ELL_SUCCESS tol = AbsTol(1) c Get the reduction handel for the sum operation call CCTK_ReductionArrayHandle(sum_handle,"sum"); if (sum_handle.lt.0) then call CCTK_WARN(1,"Cannot get reduction handle.") endif c Set boundary related variables if (CCTK_EQUALS(sor_bound,"robin")) then sw(1)=1 sw(2)=1 sw(3)=1 call Ell_GetRealKey(ierr, finf, "EllLinFlat::Bnd::Robin::inf") call Ell_GetIntKey (ierr, npow, "EllLinFlat::Bnd::Robin::falloff") end if c We have no storage for M/N if they are of size one in each direction if ((Mlinear_lsh(1).eq.1).and.(Mlinear_lsh(2).eq.1).and.(Mlinear_lsh(3).eq.1)) then Mlinear_storage=0 else Mlinear_storage=1 endif if ((Nsource_lsh(1).eq.1).and.(Nsource_lsh(2).eq.1).and.(Nsource_lsh(3).eq.1)) then Nsource_storage=0 else Nsource_storage=1 endif c Set up shorthand for the grid spacings dx=cctk_delta_space(1) dy=cctk_delta_space(2) dz=cctk_delta_space(3) verb = CCTK_Equals(elliptic_verbose,"yes").eq.1 octant = CCTK_Equals(domain,"octant").eq.1 c cheb = contains("sor_accel","cheb").ne.0 c const = contains("sor_accel","const").ne.0 c none = contains("sor_accel","none").ne.0 verb = .true. cheb = .false. none = .false. const = .false. if (verb .and. cheb) $ print *,"Chebyshev Acceleration with radius of 1" if (verb .and. const) $ print *,"SOR with omega = 1.8" if (verb .and. none) $ print *,"Un-accelearted relaxation (omega = 1)" two = 2.0D0 four = 4.0D0 resnorm = 0 ae = 1.0d0/dx**2. aw = 1.0d0/dx**2. an = 1.0d0/dy**2. as = 1.0d0/dy**2. at = 1.0d0/dz**2. ab = 1.0d0/dz**2. ac_orig = -2.0d0/dx**2. - 2.0d0/dy**2. - 2.0d0/dz**2. do sor_iteration=1,maxit c We do not know the spectral radius of the Jacobi iteration, c so we will take it to be one which empirically seems to be pretty good rjacobian = 1.0D0 c Set up the omega factor omega = 1.0D0 if (cheb) then do i=2,sor_iteration omega = 1.0D0/(1.0D0 - .25D0*rjacobian**2*omega) enddo endif if (const) then omega = 1.8 endif c Total norm of the residual zero resnorm = 0. c Start loop with Red Black ks = mod(sor_iteration,2)+2 if (cctk_lsh(3) .eq. 3) then ks = 2 endif kstep = 2 do k=ks,cctk_lsh(3)-1,kstep do j=2,cctk_lsh(2)-1 do i=2,cctk_lsh(1)-1 ac = ac_orig if (Mlinear_storage.eq.1) then ac = ac + Mlinear(i,j,k) endif residual = ac*var(i,j,k) & + ae*var(i+1,j,k) + aw*var(i-1,j,k) & + an*var(i,j+1,k) + as*var(i,j-1,k) & + at*var(i,j,k+1) + ab*var(i,j,k-1) if (Nsource_storage.eq.1) then residual = residual + Nsource(i,j,k) endif c Accumulate to the total Norm of the residual resnorm = resnorm + abs(residual) c Update var(i,j,k) = var(i,j,k) - omega*residual/ac end do end do end do c Reduce the norm call CCTK_ReduceLocalScalar(ierr, cctkGH, -1, sum_handle, $ resnorm, residual, CCTK_VARIABLE_REAL) if (ierr.ne.0) then call CCTK_WARN(1,"Reduction of norm failed!"); endif residual = residual / (cctk_gsh(1)*cctk_gsh(2)*cctk_gsh(3)) c write (*,*) "Iteration/Norm",sor_iteration,residual c Synchronize the variables call CCTK_SyncGroupWithVarI(cctkGH, var_idx) if (residual .lt. tol) then goto 123 endif c Apply boundary conditions c call Ell_GetStringKey(nchar, mybound,"EllLinFlat::Bnd") if (CCTK_EQUALS(sor_bound,"robin")) then call RobinBCVarI(ierr, cctkGH, finf, npow, sw, var_idx); if (ierr.ne.0) then call CCTK_WARN(1,"Could not apply Robin BC !") endif endif c Apply octant Symmetries call CartSymVI(ierr, cctkGH, var_idx) enddo call CCTK_INFO("SOR solver did not converge") ierr = ELL_NOCONVERGENCE 123 continue if (verb) write (*,*) "Iteration/Norm",sor_iteration,residual return end