Add ninth order dissipation operator

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

2 files changed, 34 insertions, 13 deletions

diff --git a/param.ccl b/param.ccl
index 81798b5..fe6f50f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -13,6 +13,7 @@ INT order "Dissipation order" STEERABLE=always
   3 :: "third order accurate dissipation (using a fourth derivative)"
   5 :: "fifth order accurate dissipation (using a sixth derivative)"
   7 :: "seventh order accurate dissipation (using an eighth derivative)"
+  9 :: "ninth order accurate dissipation (using a tenth derivative)"
 } 3
 
 REAL epsdis "Dissipation strength" STEERABLE=always
diff --git a/src/apply_dissipation.F77 b/src/apply_dissipation.F77
index d4a9519..a54b80c 100644
--- a/src/apply_dissipation.F77
+++ b/src/apply_dissipation.F77
@@ -7,12 +7,16 @@ c     $Header$
       
       integer   ni, nj, nk
       CCTK_REAL var(ni,nj,nk), rhs(ni,nj,nk)
-      CCTK_REAL dx(3)
+      CCTK_REAL dx(3), idx(3)
       CCTK_INT  order
       CCTK_REAL epsdis(ni,nj,nk)
       
       integer   i, j, k
       
+      idx(1) = 1 / dx(1)
+      idx(2) = 1 / dx(2)
+      idx(3) = 1 / dx(3)
+      
       if (order .eq. 1) then
          
 c$omp parallel do private(i,j,k)
@@ -20,9 +24,9 @@ c$omp parallel do private(i,j,k)
             do j = 2, nj-1
                do i = 2, ni-1
                   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))
+     $                 * (+ (var(i-1,j,k) - 2*var(i,j,k) + var(i+1,j,k)) * idx(1)
+     $                    + (var(i,j-1,k) - 2*var(i,j,k) + var(i,j+1,k)) * idx(2)
+     $                    + (var(i,j,k-1) - 2*var(i,j,k) + var(i,j,k+1)) * idx(3))
                   
                end do
             end do
@@ -36,9 +40,9 @@ c$omp parallel do private(i,j,k)
                do i = 3, ni-2
                   
                   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))
+     $                 * (+ (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)) * idx(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)) * idx(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)) * idx(3))
                   
                end do
             end do
@@ -52,9 +56,9 @@ c$omp parallel do private(i,j,k)
             do i = 4, ni-3
                
               rhs(i,j,k) = rhs(i,j,k) + epsdis(i,j,k) / 64
-     $           * (+ (var(i-3,j,k) - 6*var(i-2,j,k) + 15*var(i-1,j,k) - 20*var(i,j,k) + 15*var(i+1,j,k) - 6*var(i+2,j,k) + var(i+3,j,k)) / dx(1)
-     $              + (var(i,j-3,k) - 6*var(i,j-2,k) + 15*var(i,j-1,k) - 20*var(i,j,k) + 15*var(i,j+1,k) - 6*var(i,j+2,k) + var(i,j+3,k) ) / dx(2)
-     $              + (var(i,j,k-3) - 6*var(i,j,k-2) + 15*var(i,j,k-1) - 20*var(i,j,k) + 15*var(i,j,k+1) - 6*var(i,j,k+2) + var(i,j,k+3) ) / dx(3))
+     $           * (+ (var(i-3,j,k) - 6*var(i-2,j,k) + 15*var(i-1,j,k) - 20*var(i,j,k) + 15*var(i+1,j,k) - 6*var(i+2,j,k) + var(i+3,j,k)) * idx(1)
+     $              + (var(i,j-3,k) - 6*var(i,j-2,k) + 15*var(i,j-1,k) - 20*var(i,j,k) + 15*var(i,j+1,k) - 6*var(i,j+2,k) + var(i,j+3,k) ) * idx(2)
+     $              + (var(i,j,k-3) - 6*var(i,j,k-2) + 15*var(i,j,k-1) - 20*var(i,j,k) + 15*var(i,j,k+1) - 6*var(i,j,k+2) + var(i,j,k+3) ) * idx(3))
               
             end do
           end do
@@ -68,9 +72,25 @@ c$omp parallel do private(i,j,k)
             do i = 5, ni-4
                
               rhs(i,j,k) = rhs(i,j,k) - epsdis(i,j,k) / 256
-     $           * (+ (var(i-4,j,k) - 8*var(i-3,j,k) + 28*var(i-2,j,k) - 56*var(i-1,j,k) + 70*var(i,j,k) - 56*var(i+1,j,k) + 28*var(i+2,j,k) - 8*var(i+3,j,k) + var(i+4,j,k)) / dx(1)
-     $              + (var(i,j-4,k) - 8*var(i,j-3,k) + 28*var(i,j-2,k) - 56*var(i,j-1,k) + 70*var(i,j,k) - 56*var(i,j+1,k) + 28*var(i,j+2,k) - 8*var(i,j+3,k) + var(i,j+4,k)) / dx(2)
-     $              + (var(i,j,k-4) - 8*var(i,j,k-3) + 28*var(i,j,k-2) - 56*var(i,j,k-1) + 70*var(i,j,k) - 56*var(i,j,k+1) + 28*var(i,j,k+2) - 8*var(i,j,k+3) + var(i,j,k+4)) / dx(3))
+     $           * (+ (var(i-4,j,k) - 8*var(i-3,j,k) + 28*var(i-2,j,k) - 56*var(i-1,j,k) + 70*var(i,j,k) - 56*var(i+1,j,k) + 28*var(i+2,j,k) - 8*var(i+3,j,k) + var(i+4,j,k)) * idx(1)
+     $              + (var(i,j-4,k) - 8*var(i,j-3,k) + 28*var(i,j-2,k) - 56*var(i,j-1,k) + 70*var(i,j,k) - 56*var(i,j+1,k) + 28*var(i,j+2,k) - 8*var(i,j+3,k) + var(i,j+4,k)) * idx(2)
+     $              + (var(i,j,k-4) - 8*var(i,j,k-3) + 28*var(i,j,k-2) - 56*var(i,j,k-1) + 70*var(i,j,k) - 56*var(i,j,k+1) + 28*var(i,j,k+2) - 8*var(i,j,k+3) + var(i,j,k+4)) * idx(3))
+              
+            end do
+          end do
+        end do
+        
+      else if (order .eq. 9) then
+         
+c$omp parallel do private(i,j,k)
+        do k = 6, nk-5
+          do j = 6, nj-5
+            do i = 6, ni-5
+               
+              rhs(i,j,k) = rhs(i,j,k) + epsdis(i,j,k) / 1024
+     $           * (+ (var(i-5,j,k) - 10*var(i-4,j,k) + 45*var(i-3,j,k) - 120*var(i-2,j,k) + 210*var(i-1,j,k) - 252*var(i,j,k) + 210*var(i+1,j,k) - 120*var(i+2,j,k) + 45*var(i+3,j,k) - 10*var(i+4,j,k) + var(i+5,j,k)) * idx(1)
+     $              + (var(i,j-5,k) - 10*var(i,j-4,k) + 45*var(i,j-3,k) - 120*var(i,j-2,k) + 210*var(i,j-1,k) - 252*var(i,j,k) + 210*var(i,j+1,k) - 120*var(i,j+2,k) + 45*var(i,j+3,k) - 10*var(i,j+4,k) + var(i,j+5,k)) * idx(2)
+     $              + (var(i,j,k-5) - 10*var(i,j,k-4) + 45*var(i,j,k-3) - 120*var(i,j,k-2) + 210*var(i,j,k-1) - 252*var(i,j,k) + 210*var(i,j,k+1) - 120*var(i,j,k+2) + 45*var(i,j,k+3) - 10*var(i,j,k+4) + var(i,j,k+5)) * idx(3))
               
             end do
           end do