Fix a bug in setting extra dissipation at the outer boundary.

Add a parameter to specify the maximum epsdis for the extra dissipation at the outer boundary, if wanted.


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

2 files changed, 33 insertions, 14 deletions

diff --git a/param.ccl b/param.ccl
index 9b7dcd0..2ecde5b 100644
--- a/param.ccl
+++ b/param.ccl
@@ -87,6 +87,10 @@ REAL ob_slope "slope at outer boundary" STEERABLE=always
   0:* :: "increase dissipation"
 } 5
 
+REAL outer_boundary_max_epsdis "maximal epsdis" STEERABLE=always
+{
+  *:* :: "<0 for 'off', >=0 for maximal epsdis at the outer boundary"
+} -1.0
 
 
 
diff --git a/src/setup_epsdis.c b/src/setup_epsdis.c
index 4be7a2c..b3a3b9c 100644
--- a/src/setup_epsdis.c
+++ b/src/setup_epsdis.c
@@ -49,8 +49,6 @@ setup_epsdis (CCTK_ARGUMENTS)
     epsdisA[i]=epsdis;
   }
 
-
-
   if (extra_dissipation_at_outerbound) 
   {
     symtable = SymmetryTableHandleForGrid (cctkGH);
@@ -60,6 +58,7 @@ setup_epsdis (CCTK_ARGUMENTS)
     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++) {
@@ -67,6 +66,9 @@ setup_epsdis (CCTK_ARGUMENTS)
             index = CCTK_GFINDEX3D(cctkGH,i,j,k);
             indexP= CCTK_GFINDEX3D(cctkGH,outer_bound_npoints,j,k); 
             epsdisA[index]=epsdis+ob_slope*fabs(x[index]-x[indexP]);
+	    if (epsdisA[index] > outer_boundary_max_epsdis) {
+	      epsdisA[index] = outer_boundary_max_epsdis;
+	    }
           }
         }
       }
@@ -74,15 +76,17 @@ setup_epsdis (CCTK_ARGUMENTS)
     if(doBC[1]) {
       for (k=0;k<nk;k++) {
         for (j=0;j<nj;j++) {
-          for (i=ni;i>ni-outer_bound_npoints;i--) {
+          for (i=ni-1;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+ob_slope*fabs(x[index]-x[indexP]);
+            indexP= CCTK_GFINDEX3D(cctkGH,ni-outer_bound_npoints-1,j,k); 
+	    epsdisA[index]=epsdis+ob_slope*fabs(x[index]-x[indexP]);
+	    if (epsdisA[index] > outer_boundary_max_epsdis) {
+	      epsdisA[index] = outer_boundary_max_epsdis;
+	    }
           }
         }
       }
     }
-
     if(doBC[2]) {
       for (k=0;k<nk;k++) {
         for (j=0;j<outer_bound_npoints;j++) {
@@ -90,19 +94,25 @@ setup_epsdis (CCTK_ARGUMENTS)
             index = CCTK_GFINDEX3D(cctkGH,i,j,k);
             indexP= CCTK_GFINDEX3D(cctkGH,i,outer_bound_npoints,k); 
             epsdisA[index]=epsdis+ob_slope*fabs(y[index]-y[indexP]);
+	    if (epsdisA[index] > outer_boundary_max_epsdis) {
+	      epsdisA[index] = outer_boundary_max_epsdis;
+	    }
           }
         }
       }
     }
     if(doBC[3]) {
       for (k=0;k<nk;k++) {
-        for (j=nj;j>nj-outer_bound_npoints;j--) {
+        for (j=nj-1;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); 
+            indexP= CCTK_GFINDEX3D(cctkGH,i,nj-outer_bound_npoints-1,k); 
             epsdisA[index]=epsdis+ob_slope*fabs(y[index]-y[indexP]);
-          }
-        }
+	    if (epsdisA[index] > outer_boundary_max_epsdis) {
+	      epsdisA[index] = outer_boundary_max_epsdis;
+	    }
+	  }
+	}
       }
     }
 
@@ -113,25 +123,29 @@ setup_epsdis (CCTK_ARGUMENTS)
             index = CCTK_GFINDEX3D(cctkGH,i,j,k);
             indexP= CCTK_GFINDEX3D(cctkGH,i,j,outer_bound_npoints); 
             epsdisA[index]=epsdis+ob_slope*fabs(z[index]-z[indexP]);
+	    if (epsdisA[index] > outer_boundary_max_epsdis) {
+	      epsdisA[index] = outer_boundary_max_epsdis;
+	    }
           }
         }
       }
     }
     if(doBC[5]) {
-      for (k=nk;k>nk-outer_bound_npoints;k--) {
+      for (k=nk-1;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); 
+            indexP= CCTK_GFINDEX3D(cctkGH,i,j,nk-outer_bound_npoints-1); 
             epsdisA[index]=epsdis+ob_slope*fabs(z[index]-z[indexP]);
+	    if (epsdisA[index] > outer_boundary_max_epsdis) {
+	      epsdisA[index] = outer_boundary_max_epsdis;
+	    }
           }
         }
       }
     }
   }
 
-
-
   if (extra_dissipation_in_horizons && cctk_iteration%update_ah_every == 0)
   {
     if (verbose) {
@@ -271,4 +285,5 @@ setup_epsdis (CCTK_ARGUMENTS)
       free(inds);
     }
   }
+
 }