The CactusElliptic arrangement

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

8 files changed, 463 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..e90ec41
--- /dev/null
+++ b/README
@@ -0,0 +1,20 @@
+Cactus Code Thorn EllSOR
+Authors    : ...
+Managed by : ... <...@...........>
+Version    : ...  
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+1. Purpose of the thorn
+
+This thorn does ...
+
+2. Dependencies of the thorn
+
+This thorn additionally requires implementations and thorns ...
+
+3. Thorn distribution
+
+This thorn is available to ...
+
+4. Additional information
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..ce0621b
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,4 @@
+# Interface definition for thorn EllSOR
+# $Header$
+
+implements: ellsor
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..abaf803
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,8 @@
+# Parameter definitions for thorn EllSOR
+# $Header$
+
+shares:grid
+
+USES KEYWORD domain ""
+{
+}
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..4ac4c6b
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,7 @@
+# Schedule definitions for thorn EllSOR
+# $Header$
+
+schedule EllSOR_Register at CCTK_INITIAL
+{
+  LANG:C
+} "Register the SOR solvers"
diff --git a/src/Startup.c b/src/Startup.c
new file mode 100644
index 0000000..a7d3c22
--- /dev/null
+++ b/src/Startup.c
@@ -0,0 +1,16 @@
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+
+int EllSOR_Register(cGH *GH) {
+  void sor_confmetric(cGH *GH, int *MetricPsiI, int *FieldI, int *MI,
+                      int *NI, int *TolArray);
+
+  printf("SOR: registering sor...");
+  Ell_RegisterSolver(sor_confmetric,"sor","Ell_LinConfMetric");
+  printf("...done\n");
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..9e920e0
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn EllSOR
+# $Header$
+
+# Source files in this directory
+SRCS = Startup.c sor_wrapper.c sor_confmetric.F
+
+# Subdirectories containing source files
+SUBDIRS = 
+
diff --git a/src/sor_confmetric.F b/src/sor_confmetric.F
new file mode 100644
index 0000000..5937f87
--- /dev/null
+++ b/src/sor_confmetric.F
@@ -0,0 +1,311 @@
+ /*@@
+   @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,
+     $   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))
+      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
+
+      INTEGER maxit
+      
+      logical cheb, const, none, verb
+      integer Mlinear_storage,Nsource_storage
+c#ifdef MPI
+c      include 'mpif.h'
+c      INTEGER ierror
+c#endif
+            
+      tol = AbsTol(1)
+
+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
+
+      print *,"NStor ",Nsource_storage
+      print *,"Mlin",Mlinear_storage
+
+c     Set up shorthand for the grid spacings
+      dx=cctk_delta_space(1)
+      dy=cctk_delta_space(2)
+      dz=cctk_delta_space(3)
+
+c      write(*,*)'supper: still need to do something with params'
+c      verb = contains("elliptic_verbose","yes").ne.0
+      octant = CCTK_Equals(domain,"octant").eq.1
+c      maxit = geti("sor_maxit")
+      maxit = 100
+
+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
+
+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#ifdef MPI
+c         call mpi_allreduce(resnorm,residual,1,MPI_DOUBLE_PRECISION,
+c     $        MPI_SUM,MPI_COMM_WORLD,ierror)
+c         call synconefunc(var)
+c#else
+         residual = resnorm
+c#endif
+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 (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_wrapper.c b/src/sor_wrapper.c
new file mode 100644
index 0000000..989229a
--- /dev/null
+++ b/src/sor_wrapper.c
@@ -0,0 +1,88 @@
+ /*@@
+   @file      jacobi_wrapper.c
+   @date      Tue Aug 24 12:50:07 1999
+   @author    Gerd Lanfermann
+   @desc 
+     Provides the C wrapper to the different elliptic FORTRAN routines in this thorn
+     * These C routines are registered in the Startup routine.
+     * The LINELL_*_ARGS macro used here is provided by LinearElliptic.h
+       This cannot be modified, since this argument structure is the same 
+       across all solves which will be registered as *_hhgr3d. See LinearElliptic.c
+   @enddesc 
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+#include "cctk_FortranString.h"
+
+#include "CactusElliptic/LinearElliptic/src/LinearElliptic.h"
+
+ 
+void sor_confmetric(cGH *GH, int *MetricPsiI, int *FieldIndex, 
+		    int *MIndex, int *NIndex, int *AbsTol,int *RelTol) {
+
+  CCTK_REAL *gxx=NULL, *gxy=NULL, *gxz=NULL;
+  CCTK_REAL *gyy=NULL, *gyz=NULL, *gzz=NULL;
+  CCTK_REAL *psi=NULL;
+  CCTK_REAL *Mlinear=NULL, *Nsources=NULL;
+  CCTK_REAL *var=NULL;
+  
+  CCTK_REAL tolerance;
+  int i;
+  int toltype;
+
+  int Mlinear_lsh[3], Nsource_lsh[3];
+
+  int retcode;
+  
+  /* derive the metric data pointer from the index array. Note the ordering. */
+  gxx = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[0]);
+  gxy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[1]);
+  gxz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[2]);
+  gyy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[3]);
+  gyz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[4]);
+  gzz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[5]);
+  psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, MetricPsiI[6]);
+  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*FieldIndex);
+
+  if ((!gxx)||(!gxy)||(!gxz)||(!gyy)||(!gyz)||(!gzz)||(!psi)||(!var)) 
+    CCTK_WARN(0,"SOR_WRAPPER: One of the metric data fields, or the GF to solve could not be found!");
+  
+  /* derive the data pointer for the fields. the M/N fields are not 
+     allocated (better: are of size 1), if the passed index is negative, 
+     or we get back an empty GF of size 1 */
+  if (*MIndex>0) Mlinear   = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*MIndex);
+  if (*NIndex>0) Nsources  = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,*NIndex);
+
+
+  /* we pass the size of M/N through to frotran, so F can tell the difference 
+     between an allocated GF (Mlinear_lsh=cctk_lsh) or unallocated GF (Mlinear_lsh=1) 
+     maximal dimension is three. */
+  /*$FIXME: We need to get the group index first! 
+    printf("XXXX %d \n",CCTK_QueryGroupStorage(GH,*MIndex,NULL));$*/
+  if (GH->cctk_dim>3) 
+    CCTK_WARN(0,"This elliptic solver implementation does not do dimension>3!");
+  for (i=0;i<GH->cctk_dim;i++) {
+    if((*MIndex<0)) /*$|| (CCTK_QueryGroupStorage(GH,*MIndex,NULL)==0))$*/ 
+          Mlinear_lsh[i]=1;
+    else  Mlinear_lsh[i]=GH->cctk_lsh[i];
+    if((*NIndex<0))  /*$|| (CCTK_QueryGroupStorage(GH,*NIndex,NULL)==0))$*/ 
+          Nsource_lsh[i]=1;
+    else  Nsource_lsh[i]=GH->cctk_lsh[i];
+    printf("%d %d \n",Nsource_lsh[i],Mlinear_lsh[i]);
+  }
+ 
+
+  /* call the fortran routine */
+  FORTRAN_NAME(sor_confmetric_core3d)(_PASS_CCTK_C2F(GH),
+	 Mlinear_lsh, Mlinear,
+         Nsource_lsh, Nsources,
+	 gxx,gxy,gxz,gyy,gyz,gzz,psi, 
+	 var, AbsTol, RelTol);
+
+}
+