Add Peter's new higher order dissipation operators.

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

2 files changed, 32 insertions, 0 deletions

diff --git a/param.ccl b/param.ccl
index 00ebb03..686f17f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -11,6 +11,8 @@ INT order "Dissipation order" STEERABLE=always
 {
   1 :: "first order accurate dissipation (using a second derivative)"
   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 eight derivative)"
 } 3
 
 REAL epsdis "Dissipation strength" STEERABLE=always
diff --git a/src/apply_dissipation.F77 b/src/apply_dissipation.F77
index a246774..d67a02e 100644
--- a/src/apply_dissipation.F77
+++ b/src/apply_dissipation.F77
@@ -42,6 +42,36 @@ c     $Header$
             end do
          end do
          
+      else if (order .eq. 5) then
+         
+        do k = 4, nk-3
+          do j = 4, nj-3
+            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))
+              
+            end do
+          end do
+        end do
+        
+      else if (order .eq. 7) then
+         
+        do k = 5, nk-4
+          do j = 5, nj-4
+            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))
+              
+            end do
+          end do
+        end do
+        
       else
          
          call CCTK_WARN (0, "internal error")