removing Fortran files

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

2 files changed, 0 insertions, 597 deletions

diff --git a/src/sor_confmetric.F b/src/sor_confmetric.F
deleted file mode 100644
index bee8e8e..0000000
--- a/src/sor_confmetric.F
+++ /dev/null
@@ -1,329 +0,0 @@
- /*@@
-   @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
-
-
-
diff --git a/src/sor_flat.F b/src/sor_flat.F
deleted file mode 100644
index 904f66d..0000000
--- a/src/sor_flat.F
+++ /dev/null
@@ -1,268 +0,0 @@
- /*@@
-   @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
-
-
-