Fix a severe bug: forgot to divide by the grid spacing.

Remove spurious factor of 1/16.

Add second order dissipation, which is first order correct.

Checked stability interval.  It should now be as described.


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

5 files changed, 65 insertions, 18 deletions

diff --git a/interface.ccl b/interface.ccl
index 80a3e75..6e26e97 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -4,4 +4,4 @@
 IMPLEMENTS: Dissipation
 
 CCTK_INT FUNCTION MoLQueryEvolvedRHS (CCTK_INT IN EvolvedIndex)
-USES FUNCTION MoLQueryEvolvedRHS
+REQUIRES FUNCTION MoLQueryEvolvedRHS
diff --git a/param.ccl b/param.ccl
index b50b7da..5640ea1 100644
--- a/param.ccl
+++ b/param.ccl
@@ -5,9 +5,15 @@ BOOLEAN verbose "produce log output" STEERABLE=always
 {
 } "no"
 
+INT order "Dissipation order" STEERABLE=always
+{
+  2 :: "second order dissipation (first order accurate)"
+  4 :: "fourth order dissipation (third order accurate)"
+} 4
+
 REAL epsdis "Dissipation strength" STEERABLE=always
 {
-  *:* :: "0 for no dissipation.  Unstable for epsdis<0 and epsdis>1/"
+  *:* :: "0 for no dissipation.  Unstable for epsdis<0 and epsdis>1/3"
 } 0.2
 
 STRING vars "List of evolved grid functions that should have dissipation added" STEERABLE=always
diff --git a/src/apply_dissipation.F77 b/src/apply_dissipation.F77
index 9363824..ab78713 100644
--- a/src/apply_dissipation.F77
+++ b/src/apply_dissipation.F77
@@ -2,26 +2,51 @@ c     $Header$
 
 #include "cctk.h"
 
-      subroutine apply_dissipation (var, rhs, ni, nj, nk, epsdis)
+      subroutine apply_dissipation (var, rhs, ni, nj, nk, dx, order, epsdis)
       implicit none
       
       integer   ni, nj, nk
       CCTK_REAL var(ni,nj,nk), rhs(ni,nj,nk)
+      CCTK_REAL dx(3)
+      CCTK_INT  order
       CCTK_REAL epsdis
       
       integer   i, j, k
       
-      do k = 3, nk-2
-         do j = 3, nj-2
-            do i = 3, ni-2
-               
-               rhs(i,j,k) = rhs(i,j,k) - epsdis / 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)
-     $                 + 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)
-     $                 + 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))
-               
+      if (order .eq. 2) then
+         
+         do k = 2, nk-1
+            do j = 2, nj-1
+               do i = 2, ni-1
+                  
+                  rhs(i,j,k) = rhs(i,j,k) - epsdis
+     $                 * (+ (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))
+                  
+               end do
             end do
          end do
-      end do
+         
+      else if (order .eq. 4) then
+         
+         do k = 3, nk-2
+            do j = 3, nj-2
+               do i = 3, ni-2
+                  
+                  rhs(i,j,k) = rhs(i,j,k) - epsdis
+     $                 * (+ (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))
+                  
+               end do
+            end do
+         end do
+         
+      else
+         
+         call CCTK_WARN (0, "internal error")
+         
+      end if
       
       end
diff --git a/src/dissipation.c b/src/dissipation.c
index 362bfbb..fd59e31 100644
--- a/src/dissipation.c
+++ b/src/dissipation.c
@@ -13,6 +13,8 @@ CCTK_FNAME(apply_dissipation) (CCTK_REAL const * const var,
                                int       const * const ni,
                                int       const * const nj,
                                int       const * const nk,
+                               CCTK_REAL const * const dx,
+                               CCTK_INT  const * const order,
                                CCTK_REAL const * const epsdis);
 
 static void
@@ -38,11 +40,23 @@ apply (int const varindex, char const * const optstring, void * const arg)
   cGroup vardata, rhsdata;
   CCTK_REAL const * varptr;
   CCTK_REAL       * rhsptr;
+  CCTK_REAL dx[3];
   int n;
+  int d;
   int ierr;
   
   assert (varindex >= 0);
   
+  for (d=0; d<cctk_dim; ++d) {
+    if (cctk_nghostzones[d] < order/2) {
+      CCTK_WARN (0, "This thorn requires at least order/2 ghost zones");
+    }
+  }
+  
+  for (d=0; d<3; ++d) {
+    dx[d] = CCTK_DELTA_SPACE(d);
+  }
+  
   rhsindex = MoLQueryEvolvedRHS (varindex);
   if (rhsindex < 0) {
     char * const fullvarname = CCTK_FullName (varindex);
@@ -89,5 +103,6 @@ apply (int const varindex, char const * const optstring, void * const arg)
   assert (rhsptr);
   
   CCTK_FNAME(apply_dissipation)
-    (varptr, rhsptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2], &epsdis);
+    (varptr, rhsptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+     dx, &order, &epsdis);
 }
diff --git a/src/paramcheck.c b/src/paramcheck.c
index 3cce974..47d541a 100644
--- a/src/paramcheck.c
+++ b/src/paramcheck.c
@@ -8,14 +8,15 @@ void dissipation_paramcheck (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
+  int d;
   
   if (cctk_dim != 3) {
     CCTK_PARAMWARN ("This thorn supports only dim=3");
-    CCTK_WARN (0, "This thorn supports only dim=3");
   }
   
-  if (cctk_nghostzones[0] < 2 || cctk_nghostzones[1] < 2
-      || cctk_nghostzones[2] < 2) {
-    CCTK_PARAMWARN ("This thorn requires at least 2 ghost zones");
+  for (d=0; d<cctk_dim; ++d) {
+    if (cctk_nghostzones[d] < order/2) {
+      CCTK_PARAMWARN ("This thorn requires at least order/2 ghost zones");
+    }
   }
 }