Added a new boolean parameter "scale_with_h" that allows the user to

choose to scale the dissipation with the grid spacing h. This results
formally in the loss of one order of convergence. The default value is "no".


git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@52 f69c4107-0314-4c4f-9ad4-17e986b73f4a

7 files changed, 121 insertions, 31 deletions

diff --git a/param.ccl b/param.ccl
index d91cd2a..4e3fda5 100644
--- a/param.ccl
+++ b/param.ccl
@@ -22,6 +22,11 @@ BOOLEAN use_dissipation "Should we add dissipation"
 {
 } "no"
 
+# Note: scaling the dissipation operators with h reduces the order by one.
+BOOLEAN scale_with_h "Should we scale the dissipation with the grid spacing h"
+{
+} "no"
+
 REAL epsdis "Dissipation strength" STEERABLE=always
 {
   *:* :: "Values typical between 0 and 1"
diff --git a/src/Dissipation_2_1.F90 b/src/Dissipation_2_1.F90
index 62ba030..0904b33 100644
--- a/src/Dissipation_2_1.F90
+++ b/src/Dissipation_2_1.F90
@@ -1,18 +1,21 @@
 ! $Header$
 
 #include "cctk.h"
+#include "cctk_Parameters.h"
 
-subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
+subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs)
 
   implicit none
 
+  DECLARE_CCTK_PARAMETERS
+
   integer ::  ni, nj, nk
   CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
   CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
   CCTK_INT, dimension(6), intent(in) :: bb
   CCTK_INT, dimension(3), intent(in) :: gsize
   CCTK_REAL, dimension(3), intent(in) :: delta
-  CCTK_REAL, intent(in) :: epsdis 
+  CCTK_REAL, intent(in) :: epsilon 
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
@@ -24,7 +27,11 @@ subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 4
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 4 * delta(1) )
+  else
+    idel = epsilon / 4
+  end if
 
   if ( bb(1) == 0 ) then
     il = 1 + gsize(1)
@@ -55,7 +62,11 @@ subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 4
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 4 * delta(2) )
+  else
+    idel = epsilon / 4
+  end if
 
   if ( bb(3) == 0 ) then
     jl = 1 + gsize(2)
@@ -86,7 +97,11 @@ subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 4
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 4 * delta(3) )
+  else
+    idel = epsilon / 4
+  end if
 
   if ( bb(5) == 0 ) then
     kl = 1 + gsize(3)
diff --git a/src/Dissipation_4_2.F90 b/src/Dissipation_4_2.F90
index 16a870d..f43073e 100644
--- a/src/Dissipation_4_2.F90
+++ b/src/Dissipation_4_2.F90
@@ -1,18 +1,21 @@
 ! $Header$
 
 #include "cctk.h"
+#include "cctk_Parameters.h"
 
-subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
+subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs)
 
   implicit none
 
+  DECLARE_CCTK_PARAMETERS
+
   integer ::  ni, nj, nk
   CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
   CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
   CCTK_INT, dimension(6), intent(in) :: bb
   CCTK_INT, dimension(3), intent(in) :: gsize
   CCTK_REAL, dimension(3), intent(in) :: delta
-  CCTK_REAL, intent(in) :: epsdis 
+  CCTK_REAL, intent(in) :: epsilon 
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
@@ -24,7 +27,11 @@ subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 16
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 16 * delta(1) )
+  else
+    idel = epsilon / 16
+  end if
 
   if ( bb(1) == 0 ) then
     il = 1 + gsize(1)
@@ -75,7 +82,11 @@ subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 16
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 16 * delta(2) )
+  else
+    idel = epsilon / 16
+  end if
 
   if ( bb(3) == 0 ) then
     jl = 1 + gsize(2)
@@ -126,7 +137,11 @@ subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 16
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 16 * delta(3) )
+  else
+    idel = epsilon / 16
+  end if
 
   if ( bb(5) == 0 ) then
     kl = 1 + gsize(3)
diff --git a/src/Dissipation_4_3_min_err_coeff.F90 b/src/Dissipation_4_3_min_err_coeff.F90
index b4b6ba3..dc4178d 100644
--- a/src/Dissipation_4_3_min_err_coeff.F90
+++ b/src/Dissipation_4_3_min_err_coeff.F90
@@ -4,7 +4,8 @@
 #include "cctk_Parameters.h"
 
 
-subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl, dfl, rhs)
+subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, &
+                                               delta, epsilon, dfl, rhs)
 
   implicit none
 
@@ -17,7 +18,7 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl,
   CCTK_INT, dimension(6), intent(in) :: bb
   CCTK_INT, dimension(3), intent(in) :: gsize
   CCTK_REAL, dimension(3), intent(in) :: delta
-  CCTK_REAL, intent(in) :: epsdisl, dfl
+  CCTK_REAL, intent(in) :: epsilon, dfl
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
@@ -28,6 +29,7 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl,
   CCTK_INT :: il, ir, jl, jr, kl, kr
 
   ni = lsh(1); nj = lsh(2); nk = lsh(3)
+
   allocate ( a(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni) )
   a = zero; d = zero; b = zero; h = zero;
 
@@ -39,7 +41,11 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl,
 
   a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
 
-  idel = epsdisl /  16.0_wp
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 16 * delta(1) )
+  else
+    idel = epsilon / 16
+  end if
 
   if ( bb(1) == 0 ) then
     il = 1 + gsize(1)
@@ -132,7 +138,11 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl,
 
   a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
 
-  idel = epsdisl /  16.0_wp
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 16 * delta(2) )
+  else
+    idel = epsilon / 16
+  end if
 
   if ( bb(3) == 0 ) then
     jl = 1 + gsize(2)
@@ -226,8 +236,11 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl,
 
   a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
 
-    idel = epsdisl /  16.0_wp
-
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 16 * delta(3) )
+  else
+    idel = epsilon / 16
+  end if
 
   if ( bb(5) == 0 ) then
     kl = 1 + gsize(3)
diff --git a/src/Dissipation_6_3.F90 b/src/Dissipation_6_3.F90
index 4f35eec..1f0c74b 100644
--- a/src/Dissipation_6_3.F90
+++ b/src/Dissipation_6_3.F90
@@ -1,18 +1,21 @@
 ! $Header$
 
 #include "cctk.h"
+#include "cctk_Parameters.h"
 
-subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
+subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs)
 
   implicit none
 
+  DECLARE_CCTK_PARAMETERS
+
   integer ::  ni, nj, nk
   CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
   CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
   CCTK_INT, dimension(6), intent(in) :: bb
   CCTK_INT, dimension(3), intent(in) :: gsize
   CCTK_REAL, dimension(3), intent(in) :: delta
-  CCTK_REAL, intent(in) :: epsdis 
+  CCTK_REAL, intent(in) :: epsilon 
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
@@ -24,7 +27,11 @@ subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 64
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 64 * delta(1) )
+  else
+    idel = epsilon / 64
+  end if
 
   if ( bb(1) == 0 ) then
     il = 1 + gsize(1)
@@ -101,7 +108,11 @@ subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 64
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 64 * delta(2) )
+  else
+    idel = epsilon / 64
+  end if
 
   if ( bb(3) == 0 ) then
     jl = 1 + gsize(2)
@@ -178,7 +189,11 @@ subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 64
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 64 * delta(3) )
+  else
+    idel = epsilon / 64
+  end if
 
   if ( bb(5) == 0 ) then
     kl = 1 + gsize(3)
diff --git a/src/Dissipation_6_5_min_err_coeff.F90 b/src/Dissipation_6_5_min_err_coeff.F90
index a88d3ad..312ecef 100644
--- a/src/Dissipation_6_5_min_err_coeff.F90
+++ b/src/Dissipation_6_5_min_err_coeff.F90
@@ -6,7 +6,7 @@
 
 
 subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
-                                               delta, epsdisl, dfl, rhs)
+                                               delta, epsilon, dfl, rhs)
 
   implicit none
 
@@ -20,7 +20,7 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
   CCTK_INT, dimension(6), intent(in) :: bb
   CCTK_INT, dimension(3), intent(in) :: gsize
   CCTK_REAL, dimension(3), intent(in) :: delta
-  CCTK_REAL, intent(in) :: epsdisl, dfl
+  CCTK_REAL, intent(in) :: epsilon, dfl
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
@@ -79,7 +79,11 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
 
   end if
 
-  idel = epsdisl /  64.0_wp
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 64 * delta(1) )
+  else
+    idel = epsilon / 64
+  end if
 
   if ( bb(1) == 0 ) then
 
@@ -275,7 +279,11 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
 
   end if
 
-  idel = epsdisl /  64.0_wp
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 64 * delta(2) )
+  else
+    idel = epsilon / 64
+  end if
  
    if ( bb(3) == 0 ) then
 
@@ -471,7 +479,11 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
 
   end if
 
-  idel = epsdisl /  64.0_wp
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 64 * delta(3) )
+  else
+    idel = epsilon / 64
+  end if
 
   if ( bb(5) == 0 ) then
 
diff --git a/src/Dissipation_8_4.F90 b/src/Dissipation_8_4.F90
index bc225af..ec98407 100644
--- a/src/Dissipation_8_4.F90
+++ b/src/Dissipation_8_4.F90
@@ -1,18 +1,21 @@
 ! $Header$
 
 #include "cctk.h"
+#include "cctk_Parameters.h"
 
-subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
+subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs)
 
   implicit none
 
+  DECLARE_CCTK_PARAMETERS
+
   integer ::  ni, nj, nk
   CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
   CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
   CCTK_INT, dimension(6), intent(in) :: bb
   CCTK_INT, dimension(3), intent(in) :: gsize
   CCTK_REAL, dimension(3), intent(in) :: delta
-  CCTK_REAL, intent(in) :: epsdis 
+  CCTK_REAL, intent(in) :: epsilon 
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
@@ -23,7 +26,11 @@ subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 256
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 256 * delta(1) )
+  else
+    idel = epsilon / 256
+  end if
 
   if ( bb(1) == 0 ) then
     il = 1 + gsize(1)
@@ -134,7 +141,11 @@ subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 256
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 256 * delta(2) )
+  else
+    idel = epsilon / 256
+  end if
 
   if ( bb(3) == 0 ) then
     jl = 1 + gsize(2)
@@ -245,7 +256,11 @@ subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
 
   call set_coeff ( a )
 
-  idel = epsdis / 256
+  if ( scale_with_h > 0 ) then
+    idel = epsilon / ( 256 * delta(3) )
+  else
+    idel = epsilon / 256
+  end if
 
   if ( bb(5) == 0 ) then
     kl = 1 + gsize(3)