/*@@ @file sor.F @date Thu Mar 13 07:46:55 1997 @author Joan Masso + Paul Walker @desc The SOR solver @enddesc @@*/ #include "cctk.h" #include "cctk_arguments.h" #include "cctk_parameters.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(_CCTK_FARGUMENTS, $ Mlinear_lsh,Mlinear, $ Nsource_lsh,Nsource, $ var,var_idx, $ abstol,reltol) implicit none _DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS INTEGER CCTK_Equals 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,aw,ae,an,as,at,ab CCTK_REAL ane,anw,ase,asw,ate,atw CCTK_REAL abe,abw,atn,ats,abn,asb logical cheb, const, none, verb integer Mlinear_storage,Nsource_storage INTEGER reduction_handle,ierr c stencil size INTEGER sw(3) c#ifdef MPI c include 'mpif.h' c INTEGER ierror c#endif tol = AbsTol(1) c Get the reduction handel for the sum operation call CCTK_ReductionArrayHandle(reduction_handle,"sum"); if (reduction_handle.lt.0) then call CCTK_WARN(1,"Cannot get reduction handle.") endif 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. aw = 1. an = 1. as = 1. at = 1. ab = 1. ac =-6. ane = 0. anw = 0. ase = 0. asw = 0. ate = 0. atw = 0. abe = 0. abw = 0. atn = 0. ats = 0. abn = 0. asb = 0. do sor_iteration=1,sor_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./(1. - .25*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 = -10.38 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) c & + ane*var(i+1,j+1,k) + anw*var(i-1,j+1,k) c & + ase*var(i+1,j-1,k) + asw*var(i-1,j-1,k) c & + ate*var(i+1,j,k+1) + atw*var(i-1,j,k+1) c & + abe*var(i+1,j,k-1) + abw*var(i-1,j,k-1) c & + atn*var(i,j+1,k+1) + ats*var(i,j-1,k+1) c & + abn*var(i,j+1,k-1) + asb*var(i,j-1,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, reduction_handle, $ resnorm, residual, CCTK_VARIABLE_REAL) if (ierr.ne.0) then call CCTK_WARN(1,"Reduction of norm failed!"); endif residual = resnorm call CCTK_SyncGroupWithVarI(cctkGH, var_idx) if (residual .lt. tol) then goto 123 endif c FIXME: what about boudnary conditions, like robin, etc. c Apply octant Symmetries sw(1)=1 sw(2)=1 sw(2)=1 c call FlatBCVarI(ierr,cctkGH,sw,var_idx); c if (ierr.ne.0) then c call CCTK_WARN(1,"Could not apply BC !") c endif call ApplySymmetryVarI(cctkGH, var_idx) c if (octant) then c if (cctk_bbox(1) .eq. 1) then c var(1,:,:) = var(2,:,:) c endif c if (cctk_bbox(3) .eq. 1) then c var(:,1,:) = var(:,2,:) c endif c if (cctk_bbox(5) .eq. 1) then c var(:,:,1) = var(:,:,2) c endif c endif enddo write (*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" write (*,*) "!! WARNING: SOR SOLVER DID NOT CONVERGE !!" write (*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" 123 continue if (verb) write (*,*) "Iteration/Norm",sor_iteration,residual return end