added the sor_flat - single proc only at the moment

git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllSOR/trunk@10 fa3da13c-9f13-4301-a575-cf5b8c5e1907

1 files changed, 274 insertions, 0 deletions

diff --git a/src/sor_flat.F b/src/sor_flat.F
new file mode 100644
index 0000000..735916d
--- /dev/null
+++ b/src/sor_flat.F
@@ -0,0 +1,274 @@
+ /*@@
+   @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,
+     $   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))
+      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#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
+c      call CCTK_ReduceLocalScalar(cctkGH, ierr, -1, reduction_handle, 
+c     $        residual, resnorm, CCTK_VARIABLE_REAL) 
+c      if (ierr.ne.0) then
+c         call CCTK_WARN(1,"Reduction of norm failed!");
+c      endif
+
+c#ifdef MPI
+c         call synconefunc(var)
+c#endif
+         residual = resnorm
+
+c         write(*,*)'what is nx0: ',nx0 WHAT IS THIS ??
+c         residual = residual / (nx0 * ny0 * nz0)
+c     Stop?
+         
+         if (residual .lt. tol) then 
+            goto 123
+         endif
+
+c     FIXME: what about boudnary conditions, like robin, etc.
+c     Apply octant Symmetries
+         if (octant) then
+            if (cctk_bbox(1) .eq. 1) then 
+               var(1,:,:) = var(2,:,:)
+            endif
+            if (cctk_bbox(3) .eq. 1) then
+               var(:,1,:) = var(:,2,:)
+            endif
+           if (cctk_bbox(5) .eq. 1) then 
+               var(:,:,1) = var(:,:,2)
+            endif
+       endif
+
+      enddo
+
+      write (*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
+      write (*,*) "!! WARNING: SOR SOLVER DID NOT CONVERGE  !!"
+      write (*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
+
+ 123  continue
+
+      if (verb) write (*,*) "Iteration/Norm",sor_iteration,residual
+
+      return
+      end
+
+
+