Add Frank Herrmann's extra dissipation switches.

Two new parameters, Dissipation::extra_dissipation_in_horizons and
Dissipation::extra_dissipation_at_outerbound, increase the disspation
coefficient inside horizons and near the outer boundary.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Dissipation/trunk@10 850bcc8b-0e4f-0410-8c26-8d28fbf1eda9

6 files changed, 279 insertions, 10 deletions

diff --git a/src/apply_dissipation.F77 b/src/apply_dissipation.F77
index 74cf168..a246774 100644
--- a/src/apply_dissipation.F77
+++ b/src/apply_dissipation.F77
@@ -9,7 +9,7 @@ c     $Header$
       CCTK_REAL var(ni,nj,nk), rhs(ni,nj,nk)
       CCTK_REAL dx(3), dt
       CCTK_INT  order
-      CCTK_REAL epsdis
+      CCTK_REAL epsdis(ni,nj,nk)
       
       integer   i, j, k
       
@@ -18,8 +18,7 @@ c     $Header$
          do k = 2, nk-1
             do j = 2, nj-1
                do i = 2, ni-1
-                  
-                  rhs(i,j,k) = rhs(i,j,k) + epsdis
+                  rhs(i,j,k) = rhs(i,j,k) + epsdis(i,j,k)
      $                 * (+ (var(i-1,j,k) - 2*var(i,j,k) + var(i+1,j,k)) / dx(1)
      $                    + (var(i,j-1,k) - 2*var(i,j,k) + var(i,j+1,k)) / dx(2)
      $                    + (var(i,j,k-1) - 2*var(i,j,k) + var(i,j,k+1)) / dx(3))
@@ -34,7 +33,7 @@ c     $Header$
             do j = 3, nj-2
                do i = 3, ni-2
                   
-                  rhs(i,j,k) = rhs(i,j,k) - epsdis / 16
+                  rhs(i,j,k) = rhs(i,j,k) - epsdis(i,j,k) / 16
      $                 * (+ (var(i-2,j,k) - 4*var(i-1,j,k) + 6*var(i,j,k) - 4*var(i+1,j,k) + var(i+2,j,k)) / dx(1)
      $                    + (var(i,j-2,k) - 4*var(i,j-1,k) + 6*var(i,j,k) - 4*var(i,j+1,k) + var(i,j+2,k)) / dx(2)
      $                    + (var(i,j,k-2) - 4*var(i,j,k-1) + 6*var(i,j,k) - 4*var(i,j,k+1) + var(i,j,k+2)) / dx(3))
diff --git a/src/basegrid.c b/src/basegrid.c
index ad01c26..ad54931 100644
--- a/src/basegrid.c
+++ b/src/basegrid.c
@@ -19,4 +19,8 @@ void dissipation_basegrid (CCTK_ARGUMENTS)
       CCTK_WARN (0, "This thorn requires at least (order+1)/2 ghost zones");
     }
   }
+
+  for (d=0;d<cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2];d++) {
+    epsdisA[d]=epsdis;
+  }
 }
diff --git a/src/dissipation.c b/src/dissipation.c
index 6191cc0..e04c8d5 100644
--- a/src/dissipation.c
+++ b/src/dissipation.c
@@ -2,6 +2,7 @@
 
 #include <assert.h>
 #include <stdlib.h>
+#include <stdio.h>
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
@@ -45,6 +46,9 @@ apply (int const varindex, char const * const optstring, void * const arg)
   int n;
   int d;
   int ierr;
+  int ni,nj,nk;
+  int i;
+  CCTK_REAL minrad,maxrad;
   
   assert (varindex >= 0);
   
@@ -86,25 +90,25 @@ apply (int const varindex, char const * const optstring, void * const arg)
   assert (vargroup >= 0);
   rhsgroup = CCTK_GroupIndexFromVarI (rhsindex);
   assert (rhsgroup >= 0);
-  
+
   ierr = CCTK_GroupData (vargroup, &vardata);
   assert (!ierr);
   ierr = CCTK_GroupData (rhsgroup, &rhsdata);
   assert (!ierr);
-  
+
   assert (vardata.grouptype == CCTK_GF);
   assert (vardata.vartype == CCTK_VARIABLE_REAL);
   assert (vardata.dim == cctk_dim);
   assert (rhsdata.grouptype == CCTK_GF);
   assert (rhsdata.vartype == CCTK_VARIABLE_REAL);
   assert (rhsdata.dim == cctk_dim);
-  
+
   varptr = CCTK_VarDataPtrI (cctkGH, 0, varindex);
   assert (varptr);
   rhsptr = CCTK_VarDataPtrI (cctkGH, 0, rhsindex);
   assert (rhsptr);
-  
+
   CCTK_FNAME(apply_dissipation)
     (varptr, rhsptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
-     dx, &dt, &order, &epsdis);
+     dx, &dt, &order, epsdisA);
 }
diff --git a/src/make.code.defn b/src/make.code.defn
index b1abd5b..a0c8c1f 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = apply_dissipation.F77 basegrid.c dissipation.c
+SRCS = apply_dissipation.F77 basegrid.c dissipation.c paramcheck.c setup_epsdis.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/paramcheck.c b/src/paramcheck.c
new file mode 100644
index 0000000..d3629fe
--- /dev/null
+++ b/src/paramcheck.c
@@ -0,0 +1,29 @@
+/* $Header$ */
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+void
+dissipation_paramcheck (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  int i;
+  int want_horizon;
+  
+  if (extra_dissipation_in_horizons)
+  {
+    want_horizon = 0;
+    for (i=0; i<100; ++i)
+    {
+      want_horizon = want_horizon | (horizon_number[i] >= 0);
+    }
+    
+    if (want_horizon && ! CCTK_IsFunctionAliased ("HorizonRadiusInDirection"))
+    {
+      CCTK_PARAMWARN ("The aliased function \"HorizonRadiusInDirection\" must be defined when the parameter \"extra_dissipation_in_horizons\" is set and one of the sources is AHFinderDirect");
+    }
+  }
+}
diff --git a/src/setup_epsdis.c b/src/setup_epsdis.c
new file mode 100644
index 0000000..9a3b2cf
--- /dev/null
+++ b/src/setup_epsdis.c
@@ -0,0 +1,233 @@
+/* $Header$ */
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#define MAXDIM 3
+
+void
+setup_epsdis (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  int ni,nj,nk;
+  int i,j,k,s;
+  int index,indexP;
+  int ierr;
+  int npts;
+  int *inds;
+  CCTK_REAL *xa,*ya,*za,*rads;
+  CCTK_REAL maxrad;
+  CCTK_REAL xmin,xmax, ymin,ymax, zmin,zmax;
+  CCTK_REAL odx,ody,odz;
+  CCTK_REAL radp;
+  const CCTK_INT  MAXSURFNUM=100; /* XXX hard limit */
+  CCTK_INT doBC[2*MAXDIM],symbnd[2*MAXDIM];  
+  CCTK_INT symtable;
+
+  if (verbose) {
+    CCTK_VInfo(CCTK_THORNSTRING,"Setting up spatially varying dissipation at T=%g",
+               (double)cctk_time);
+  }
+
+  ni=cctk_lsh[0];
+  nj=cctk_lsh[1];
+  nk=cctk_lsh[2];
+
+  for (i=0;i<ni*nj*nk;i++) {
+    epsdisA[i]=epsdis;
+  }
+
+
+
+  if (extra_dissipation_at_outerbound) 
+  {
+    symtable = SymmetryTableHandleForGrid (cctkGH);
+    if (symtable < 0) CCTK_WARN (1, "unable to get symmetry table");
+    ierr=Util_TableGetIntArray(symtable, 6, symbnd, "symmetry_handle");
+    if (ierr != 6) CCTK_WARN (1, "unable to get symmetry handle");
+    for (i = 0; i < 6; i++) {
+      doBC[i] = cctk_bbox[i]!=0 && symbnd[i] < 0;
+    }
+    if(doBC[0]) {
+      for (k=0;k<nk;k++) {
+        for (j=0;j<nj;j++) {
+          for (i=0;i<outer_bound_npoints;i++) {
+            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+            indexP= CCTK_GFINDEX3D(cctkGH,outer_bound_npoints,j,k); 
+            epsdisA[index]=epsdis+epsdis*ob_slope*abs(x[index]-x[indexP]);
+          }
+        }
+      }
+    }
+    if(doBC[1]) {
+      for (k=0;k<nk;k++) {
+        for (j=0;j<nj;j++) {
+          for (i=ni;i>ni-outer_bound_npoints;i--) {
+            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+            indexP= CCTK_GFINDEX3D(cctkGH,ni-outer_bound_npoints,j,k); 
+            epsdisA[index]=epsdis+epsdis*ob_slope*abs(x[index]-x[indexP]);
+          }
+        }
+      }
+    }
+
+    if(doBC[2]) {
+      for (k=0;k<nk;k++) {
+        for (j=0;j<outer_bound_npoints;j++) {
+          for (i=0;i<ni;i++) {
+            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+            indexP= CCTK_GFINDEX3D(cctkGH,i,outer_bound_npoints,k); 
+            epsdisA[index]=epsdis+epsdis*ob_slope*abs(y[index]-y[indexP]);
+          }
+        }
+      }
+    }
+    if(doBC[3]) {
+      for (k=0;k<nk;k++) {
+        for (j=nj;j>nj-outer_bound_npoints;j--) {
+          for (i=0;i<ni;i++) {
+            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+            indexP= CCTK_GFINDEX3D(cctkGH,i,nj-outer_bound_npoints,k); 
+            epsdisA[index]=epsdis+epsdis*ob_slope*abs(y[index]-y[indexP]);
+          }
+        }
+      }
+    }
+
+    if(doBC[4]) {
+      for (k=0;k<outer_bound_npoints;k++) {
+        for (j=0;j<nj;j++) {
+          for (i=0;i<ni;i++) {
+            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+            indexP= CCTK_GFINDEX3D(cctkGH,i,j,outer_bound_npoints); 
+            epsdisA[index]=epsdis+epsdis*ob_slope*abs(z[index]-z[indexP]);
+          }
+        }
+      }
+    }
+    if(doBC[5]) {
+      for (k=nk;k>nk-outer_bound_npoints;k--) {
+        for (j=0;j<nj;j++) {
+          for (i=0;i<ni;i++) {
+            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
+            indexP= CCTK_GFINDEX3D(cctkGH,i,j,nk-outer_bound_npoints); 
+            epsdisA[index]=epsdis+epsdis*ob_slope*abs(z[index]-z[indexP]);
+          }
+        }
+      }
+    }
+  }
+
+
+
+  if (extra_dissipation_in_horizons && cctk_iteration%update_ah_every == 0)
+  {
+    if (verbose) {
+      CCTK_INFO("Linear Interpolation into AH surfaces");
+    }
+
+    for (s=0;s<MAXSURFNUM;s++)
+    {
+      if (surface_number[s]==-1 && horizon_number[s]==-1) {
+        continue;
+      }
+      
+      if (surface_number[s]<0 || horizon_number[s]<=0) {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Invalid specification for horizon %d", s);
+        continue;
+      }
+      
+      if (! sf_valid[surface_number[s]]) {
+        if (verbose) {
+          CCTK_VInfo(CCTK_THORNSTRING,
+                     "Invalid Surface: s=%d, sf_va=%d, surf_no=%d, ah_no=%d",
+                     s,sf_valid[surface_number[s]],surface_number[s],horizon_number[s]);
+        }
+        continue;
+      }
+      
+      if (! CCTK_IsFunctionAliased ("HorizonRadiusInDirection"))
+      {
+        CCTK_WARN (1, "The aliased function \"HorizonRadiusInDirection\" must be defined when the parameter \"extra_dissipation_in_horizons\" is set and one of the sources is AHFinderDirect");
+        continue;
+      }
+
+      maxrad=sf_max_radius[surface_number[s]];
+      odx=sf_origin_x[surface_number[s]];
+      ody=sf_origin_y[surface_number[s]];
+      odz=sf_origin_z[surface_number[s]];
+      xmin=odx-maxrad;
+      xmax=odx+maxrad;
+      ymin=ody-maxrad;
+      ymax=ody+maxrad;
+      zmin=odz-maxrad;
+      zmax=odz+maxrad;
+
+      npts=0;
+      for (i=0;i<ni*nj*nk;i++) {
+        if (  x[i]<=xmax&&x[i]>=xmin
+            &&y[i]<=ymax&&y[i]>=ymin
+            &&z[i]<=zmax&&z[i]>=zmin)
+        {
+          npts++;
+        }
+      }
+
+      xa=(CCTK_REAL *)   malloc(npts*sizeof(CCTK_REAL));
+      ya=(CCTK_REAL *)   malloc(npts*sizeof(CCTK_REAL));
+      za=(CCTK_REAL *)   malloc(npts*sizeof(CCTK_REAL));
+      rads=(CCTK_REAL *) malloc(npts*sizeof(CCTK_REAL));
+      inds=(CCTK_INT *)  malloc(npts*sizeof(CCTK_INT));
+
+      for (i=0;i<npts;i++) {
+        rads[i]=0;
+        inds[i]=0;
+      }
+
+      j=0;
+      for (i=0;i<ni*nj*nk;i++) {
+        if (  x[i]<=xmax&&x[i]>=xmin
+            &&y[i]<=ymax&&y[i]>=ymin
+            &&z[i]<=zmax&&z[i]>=zmin)
+        {
+          xa[j]=x[i];
+          ya[j]=y[i];
+          za[j]=z[i];
+          inds[j]=i;
+          j++;
+        }
+      }
+
+      ierr=HorizonRadiusInDirection(horizon_number[s],
+                                    npts,
+                                    xa, ya, za,  rads);
+      assert(!ierr);
+
+      for (i=0;i<npts;i++) {
+        radp=sqrt(xa[i]*xa[i]+ya[i]*ya[i]+za[i]*za[i]);
+        if (radp<=rads[i]) {
+          epsdisA[inds[i]]=epsdis+ ah_slope*(rads[i]-radp)*epsdis;
+        }
+      }
+
+      free(xa);
+      free(ya);
+      free(za);
+      free(rads);
+      free(inds);
+    }
+  }
+}