/*@@ @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_confmetric_core3d(_CCTK_FARGUMENTS, $ Mlinear_lsh,Mlinear, $ Nsource_lsh,Nsource, $ gxx,gxy,gxz,gyy,gyz,gzz, & psi,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 gxx(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), gxy(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL gxz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), gyy(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL gyz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), gzz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL psi(cctk_lsh(1),cctk_lsh(2),cctk_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 pm4, tmp CCTK_REAL psixp, psiyp, psizp 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 Upper metric and determinant. What a pig. CCTK_REAL uxx(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), uyy(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL uzz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), uxy(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL uxz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), uyz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL det(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) 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 CCTK_REAL finf CCTK_INT npow logical cheb, const, none, verb integer Mlinear_storage,Nsource_storage INTEGER reduction_handle,ierr c stencil size INTEGER sw(3) 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 Set boundary related variables if (CCTK_EQUALS(sor_bound,"robin")) then sw(1)=1 sw(2)=1 sw(2)=1 call Ell_GetRealKey(ierr, finf, "EllLinConfMetric::Bnd::Robin::inf") call Ell_GetIntKey (ierr, npow, "EllLinConfMetric::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 none = .false. const = .false. cheb = .false. if (CCTK_EQUALS(sor_accel,"const")) then const = .true. else if (CCTK_EQUALS(sor_accel,"cheb")) then cheb = .true. else none = .true. endif 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.0d0 c compute determinant det= -(gxz**2*gyy) + 2*gxy*gxz*gyz - gxx*gyz**2 & - gxy**2*gzz + gxx*gyy*gzz c invert metric uxx=(-gyz**2 + gyy*gzz)/det/psi**4 uxy=(gxz*gyz - gxy*gzz)/det/psi**4 uyy=(-gxz**2 + gxx*gzz)/det/psi**4 uxz=(-gxz*gyy + gxy*gyz)/det/psi**4 uyz=(gxy*gxz - gxx*gyz)/det/psi**4 uzz=(-gxy**2 + gxx*gyy)/det/psi**4 det = det * psi**12 if (Mlinear_storage.eq.1) then Mlinear = Mlinear*sqrt(det) endif if (Nsource_storage.eq.1) then Nsource = Nsource*sqrt(det) endif c Re-scaled uppwemetric. PATCHY but effective.... watch out: we get it back after c the solve. uxx=uxx/(2.*dx*dx)*sqrt(det) uyy=uyy/(2.*dy*dy)*sqrt(det) uzz=uzz/(2.*dz*dz)*sqrt(det) uxy=uxy/(4.*dx*dy)*sqrt(det) uxz=uxz/(4.*dx*dz)*sqrt(det) uyz=uyz/(4.*dy*dz)*sqrt(det) 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 c which empirically seems to be pretty good rjacobian = 1. c Set up the omega factor omega = 1. 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 End of first-iteration stuff 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 = -uxx(i+1,j,k) -2.*uxx(i,j,k) -uxx(i-1,j,k) & -uyy(i,j+1,k) -2.*uyy(i,j,k) -uyy(i,j-1,k) & -uzz(i,j,k+1) -2.*uzz(i,j,k) -uzz(i,j,k-1) if (Mlinear_storage.eq.1) then ac = ac + Mlinear(i,j,k) endif ae = uxx(i+1,j,k)+uxx(i,j,k) aw = uxx(i-1,j,k)+uxx(i,j,k) an = uyy(i,j+1,k)+uyy(i,j,k) as = uyy(i,j-1,k)+uyy(i,j,k) at = uzz(i,j,k+1)+uzz(i,j,k) ab = uzz(i,j,k-1)+uzz(i,j,k) ane = uxy(i,j+1,k)+uxy(i+1,j,k) anw = - uxy(i-1,j,k)-uxy(i,j+1,k) ase = - uxy(i+1,j,k)-uxy(i,j-1,k) asw = uxy(i-1,j,k)+uxy(i,j-1,k) ate = uxz(i,j,k+1)+uxz(i+1,j,k) atw = - uxz(i-1,j,k)-uxz(i,j,k+1) abe = - uxz(i+1,j,k)-uxz(i,j,k-1) abw = uxz(i-1,j,k)+uxz(i,j,k-1) atn = uyz(i,j+1,k)+uyz(i,j,k+1) ats = - uyz(i,j,k+1)-uyz(i,j-1,k) abn = - uyz(i,j,k-1)-uyz(i,j+1,k) asb = uyz(i,j,k-1)+uyz(i,j-1,k) 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) & + ane*var(i+1,j+1,k) + anw*var(i-1,j+1,k) & + ase*var(i+1,j-1,k) + asw*var(i-1,j-1,k) & + ate*var(i+1,j,k+1) + atw*var(i-1,j,k+1) & + abe*var(i+1,j,k-1) + abw*var(i-1,j,k-1) & + atn*var(i,j+1,k+1) + ats*var(i,j-1,k+1) & + 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 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 residual = residual / (cctk_gsh(1)*cctk_gsh(2)*cctk_gsh(3)) call CCTK_SyncGroupWithVarI(cctkGH, var_idx) if (residual .lt. tol) then goto 123 endif c Apply Robin boundary 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 write (*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" write (*,*) "!! WARNING: SOR SOLVER DID NOT CONVERGE !!" write (*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" 123 continue if (Mlinear_storage.eq.1) then Mlinear = Mlinear/sqrt(det) endif if (Nsource_storage.eq.1) then Nsource = Nsource/sqrt(det) endif if (verb) write (*,*) "Iteration/Norm",sor_iteration,residual return end