-new parameter ah_max_epsdis

-bugfix for ah not at (0,0,0)


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

2 files changed, 12 insertions, 2 deletions

diff --git a/param.ccl b/param.ccl
index 30578a0..29b8bc0 100644
--- a/param.ccl
+++ b/param.ccl
@@ -39,6 +39,11 @@ REAL ah_slope "Slope inside AH" STEERABLE=always
   *:* :: "0 for no dissipation.  Unstable when total epsdis<0 or total epsdis>1/3"
 } 0.2
 
+REAL ah_max_epsdis "maximal epsdis" STEERABLE=always
+{
+  *:* :: "<0 for 'off', >=0 for maximal epsdis in horizon"
+} -1.0
+
 INT surface_number[100] "horizon number for extra dissipation in horizons -- SphericalSurface number" STEERABLE=always
 {
   -1  :: "do not use a spherical surface"
diff --git a/src/setup_epsdis.c b/src/setup_epsdis.c
index 5acb649..6239e6c 100644
--- a/src/setup_epsdis.c
+++ b/src/setup_epsdis.c
@@ -154,7 +154,8 @@ setup_epsdis (CCTK_ARGUMENTS)
         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]);
+                     s,sf_valid[surface_number[s]],
+                     surface_number[s],horizon_number[s]);
         }
         continue;
       }
@@ -217,9 +218,13 @@ setup_epsdis (CCTK_ARGUMENTS)
       assert(!ierr);
 
       for (i=0;i<npts;i++) {
-        radp=sqrt(xa[i]*xa[i]+ya[i]*ya[i]+za[i]*za[i]);
+        radp=sqrt((xa[i]-odx)*(xa[i]-odx)+(ya[i]-ody)*(ya[i]-ody)+
+                  (za[i]-odz)*(za[i]-odz));
         if (radp<=rads[i]) {
           epsdisA[inds[i]]=epsdis+ ah_slope*(rads[i]-radp);
+          if (epsdisA[inds[i]] > ah_max_epsdis) {
+            epsdisA[inds[i]] = ah_max_epsdis;
+          }
         }
       }