robin for sor

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

4 files changed, 65 insertions, 42 deletions

diff --git a/param.ccl b/param.ccl
index e71e986..df52af5 100644
--- a/param.ccl
+++ b/param.ccl
@@ -19,3 +19,10 @@ CCTK_INT maxit "Maximal number of iteration in the SOR loop"
 {
 0:* :: ""
 } 100
+
+KEYWORD sor_bound "Boundary condition to apply in SOR loop"
+{
+  "flat"     :: "flat / copying boundary conditions"
+  "constant" :: "constant boundary conditions"
+  "robin"    :: "Robin type boundary conditions"
+} "flat"
diff --git a/src/sor_confmetric.F b/src/sor_confmetric.F
index fafc8a5..8d10c79 100644
--- a/src/sor_confmetric.F
+++ b/src/sor_confmetric.F
@@ -26,7 +26,7 @@
      $   Mlinear_lsh,Mlinear,
      $   Nsource_lsh,Nsource,
      $   gxx,gxy,gxz,gyy,gyz,gzz,
-     &   psi,var,
+     &   psi,var, var_idx,
      $   abstol,reltol)
 
       implicit none
@@ -46,6 +46,8 @@
       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
@@ -94,9 +96,19 @@ c     Coeeficients for the solver: 19 point stencil...
 
       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     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
@@ -253,33 +265,36 @@ c     Update
             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?
          
+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
-         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)
+
+c     Apply Robin boundary 
+c         write (*,*) "Robin ?"
+         if (CCTK_EQUALS(sor_bound,"robin")) then
+c            write (*,*) "Robin"
+            sw(1)=1
+            sw(2)=1
+            sw(2)=1
+            call RobinBCVarI(ierr, cctkGH, 1.0, 1, sw, var_idx);
+            if (ierr.ne.0) then
+               call CCTK_WARN(1,"Could not Robin BC !")
             endif
-       endif
+         endif
+
+c     Apply octant Symmetries         
+         call CartSymBCVarI(ierr, cctkGH, var_idx)
 
       enddo
 
diff --git a/src/sor_flat.F b/src/sor_flat.F
index a6d8b23..98781e1 100644
--- a/src/sor_flat.F
+++ b/src/sor_flat.F
@@ -228,25 +228,27 @@ c     Reduce the norm
             call CCTK_WARN(1,"Reduction of norm failed!");
          endif
          residual = resnorm
-         
+      
+c     Synchronize the variables
          call CCTK_SyncGroupWithVarI(cctkGH, var_idx)
          
          if (residual .lt. tol) then 
             goto 123
          endif
 
-c     FIXME: what about boundary conditions, like robin, etc.
-c     Apply octant Symmetries
+c     Apply Robin boundary 
+         if (CCTK_EQUALS(sor_bound,"robin")) then
+            sw(1)=1
+            sw(2)=1
+            sw(2)=1
+            call RobinBCVarI(ierr, cctkGH, 1.0, 1, sw, var_idx);
+            if (ierr.ne.0) then
+               call CCTK_WARN(1,"Could not apply Robin BC !")
+            endif
+         endif
 
-c         sw(1)=1
-c         sw(2)=1
-c         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 CartSymBCVarI(ierr,cctkGH, var_idx)
+c     Apply octant Symmetries         
+         call CartSymBCVarI(ierr, cctkGH, var_idx)
 
       enddo
 
diff --git a/src/sor_wrapper.c b/src/sor_wrapper.c
index 93057cf..956a6cd 100644
--- a/src/sor_wrapper.c
+++ b/src/sor_wrapper.c
@@ -34,7 +34,7 @@ void sor_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex,
   CCTK_REAL *gyy=NULL, *gyz=NULL, *gzz=NULL;
   CCTK_REAL *psi=NULL;
   CCTK_REAL *Mlinear=NULL, *Nsources=NULL;
-  CCTK_REAL *var=NULL;
+  CCTK_REAL *Field  =NULL;
   
   CCTK_REAL tolerance;
   int i;
@@ -54,9 +54,10 @@ void sor_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex,
   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)) 
+  Field = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,FieldIndex);
+
+  if ((!gxx)||(!gxy)||(!gxz)||(!gyy)||(!gyz)||(!gzz)||(!psi)||(!Field)) 
     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 
@@ -86,7 +87,7 @@ void sor_confmetric(cGH *GH, int *MetricPsiI, int FieldIndex,
 	 Mlinear_lsh, Mlinear,
          Nsource_lsh, Nsources,
 	 gxx,gxy,gxz,gyy,gyz,gzz,psi, 
-	 var, AbsTol, RelTol);
+	 Field, &FieldIndex, AbsTol, RelTol);
 
 }
  
@@ -103,8 +104,6 @@ void sor_flat(cGH *GH, int FieldIndex, int MIndex,
   int Mlinear_lsh[3], Nsource_lsh[3];
 
   int retcode;
-  
-  
 
   Field = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,FieldIndex); 
   if (MIndex>0) Mlinear   = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex);