Introduce a boolean parameter to switch of SBP near outer boundary. Instead

switch to progressively lower order centered stencil as the boundery is
approached, similar to what is done with the second order and upwind stencils.


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

8 files changed, 450 insertions, 26 deletions

diff --git a/param.ccl b/param.ccl
index d4314cf..d382d36 100644
--- a/param.ccl
+++ b/param.ccl
@@ -15,6 +15,10 @@ KEYWORD operator_type "Type of operator" STEERABLE=always
   "Optimized" :: "Optimized for performance"
 } "Optimized"
 
+BOOLEAN sbp_1st_deriv "Should the 1st derivative operator be SBP" STEERABLE=always
+{
+} "yes"
+  
 BOOLEAN sbp_2nd_deriv "Should the 2nd derivative operator be SBP" STEERABLE=always
 {
 } "yes"
diff --git a/src/All_Coeffs_mod.F90 b/src/All_Coeffs_mod.F90
index 1886497..55a76ab 100644
--- a/src/All_Coeffs_mod.F90
+++ b/src/All_Coeffs_mod.F90
@@ -28,6 +28,17 @@ module All_Coeffs_mod
 
     end subroutine coeffs_1_2_1
 
+    subroutine coeffs_1_2 ( a, q )
+
+      CCTK_REAL, dimension(1), intent(OUT) :: a
+      CCTK_REAL, dimension(1,1), intent(OUT) :: q
+
+      a(1) = 1.0_wp/2.0_wp
+
+      q(1,1) = zero
+
+    end subroutine coeffs_1_2
+
     subroutine coeffs_up_2_1 ( a1, q1, a2, q2 )
 
       CCTK_REAL, dimension(-1:1), intent(OUT) :: a1, a2
@@ -161,6 +172,24 @@ module All_Coeffs_mod
 
     end subroutine coeffs_1_4_2
 
+    subroutine coeffs_1_4 ( a, q )
+
+      CCTK_REAL, dimension(2), intent(OUT) :: a
+      CCTK_REAL, dimension(3,2), intent(OUT) :: q
+
+      a(1) = 2.0_wp/3.0_wp
+      a(2) = -1.0_wp/12.0_wp
+
+      q(1,1) = zero
+      q(2,1) = zero
+      q(3,1) = zero
+
+      q(1,2) = -1.0_wp/2.0_wp
+      q(2,2) = zero
+      q(3,2) = 1.0_wp/2.0_wp
+
+    end subroutine coeffs_1_4
+
     subroutine coeffs_up_4_2 ( a1, q1, a2, q2 )
 
       CCTK_REAL, dimension(-2:2), intent(OUT) :: a1, a2
@@ -526,6 +555,35 @@ module All_Coeffs_mod
 
     end subroutine coeffs_1_6_3_opt
 
+    subroutine coeffs_1_6 ( a, q )
+
+      CCTK_REAL, dimension(3), intent(OUT) :: a
+      CCTK_REAL, dimension(5,3), intent(OUT) :: q
+
+      a(1) = 3.0_wp/4.0_wp
+      a(2) = -3.0_wp/20.0_wp
+      a(3) = 1.0_wp/60.0_wp
+
+      q(1,1) = zero
+      q(2,1) = zero
+      q(3,1) = zero
+      q(4,1) = zero
+      q(5,1) = zero
+
+      q(1,2) = -1.0_wp/2.0_wp
+      q(2,2) = zero
+      q(3,2) = 1.0_wp/2.0_wp
+      q(4,2) = zero
+      q(5,2) = zero
+
+      q(1,3) = 1.0_wp/12.0_wp
+      q(2,3) = -2.0_wp/3.0_wp
+      q(3,3) = zero
+      q(4,3) = 2.0_wp/3.0_wp
+      q(5,3) = -1.0_wp/12.0_wp
+
+    end subroutine coeffs_1_6
+
     subroutine coeffs_up_6_3_opt ( a1, q1, a2, q2 )
 
       CCTK_REAL, dimension(-3:3), intent(OUT) :: a1, a2
@@ -1129,6 +1187,50 @@ module All_Coeffs_mod
 
     end subroutine coeffs_1_8_4_opt
 
+    subroutine coeffs_1_8 ( a, q )
+
+      CCTK_REAL, dimension(4), intent(OUT) :: a
+      CCTK_REAL, dimension(7,4), intent(OUT) :: q
+
+      a(1) = 4.0_wp/5.0_wp
+      a(2) = -1.0_wp/5.0_wp
+      a(3) = 4.0_wp/105.0_wp
+      a(4) = -1.0_wp/280.0_wp
+
+      q(1,1) = zero
+      q(2,1) = zero
+      q(3,1) = zero
+      q(4,1) = zero
+      q(5,1) = zero
+      q(6,1) = zero
+      q(7,1) = zero
+
+      q(1,2) = -1.0_wp/2.0_wp
+      q(2,2) = zero
+      q(3,2) = 1.0_wp/2.0_wp
+      q(4,2) = zero
+      q(5,2) = zero
+      q(6,2) = zero
+      q(7,1) = zero
+
+      q(1,3) = 1.0_wp/12.0_wp
+      q(2,3) = -2.0_wp/3.0_wp
+      q(3,3) = zero
+      q(4,3) = 2.0_wp/3.0_wp
+      q(5,3) = -1.0_wp/12.0_wp
+      q(6,3) = zero
+      q(7,3) = zero
+
+      q(1,4) = -1.0_wp/60.0_wp
+      q(2,4) = 3.0_wp/20.0_wp
+      q(3,4) = -3.0_wp/4.0_wp
+      q(4,4) = zero
+      q(5,4) = 3.0_wp/4.0_wp
+      q(6,4) = -3.0_wp/20.0_wp
+      q(7,4) = 1.0_wp/60.0_wp
+
+    end subroutine coeffs_1_8
+
     subroutine coeffs_up_8_4_opt ( a1, q1, a2, q2 )
 
       CCTK_REAL, dimension(-4:4), intent(OUT) :: a1, a2
diff --git a/src/Coefficients_2.F90 b/src/Coefficients_2.F90
new file mode 100644
index 0000000..122961b
--- /dev/null
+++ b/src/Coefficients_2.F90
@@ -0,0 +1,56 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine set_coeff_2 ( nsize, loc_order, bb, gsize, imin, imax, dd )
+
+  use All_Coeffs_mod
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(IN) :: nsize, loc_order
+  CCTK_INT, dimension(2), intent(IN) :: bb
+  CCTK_INT, intent(IN) :: gsize
+  CCTK_INT, dimension(nsize), intent(OUT) :: imin, imax
+  CCTK_REAL, dimension(nsize,nsize), intent(OUT) :: dd
+
+  CCTK_REAL, dimension(1), save :: a
+  CCTK_REAL, dimension(1,1), save :: q 
+
+  CCTK_INT :: i, il, ir
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    call coeffs_1_2 ( a, q )
+    first = .false.
+  end if
+
+  dd = zero
+  imin = 0
+  imax = -1
+
+  if ( bb(1) == 0 ) then
+    il = 1 + gsize
+  else
+    dd(1,1) = q(1,1);
+    imin(1) = 1; imax(1) = 1
+    il = 2
+  end if
+  if ( bb(2) == 0 ) then
+    ir = nsize - gsize
+  else
+    dd(nsize,nsize) = q(1,1)
+    imin(nsize) = nsize; imax(nsize) = nsize
+    ir = nsize - 1
+  end if
+  do i = il, ir
+    dd(i-1,i) = -a(1); 
+    dd(i+1,i) = a(1)
+    imin(i) = i-1; imax(i) = i+1
+  end do
+
+end subroutine set_coeff_2
diff --git a/src/Coefficients_4.F90 b/src/Coefficients_4.F90
new file mode 100644
index 0000000..2211559
--- /dev/null
+++ b/src/Coefficients_4.F90
@@ -0,0 +1,63 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine set_coeff_4 ( nsize, loc_order, bb, gsize, imin, imax, dd )
+
+  use All_Coeffs_mod
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(IN) :: nsize, loc_order
+  CCTK_INT, dimension(2), intent(IN) :: bb
+  CCTK_INT, intent(IN) :: gsize
+  CCTK_INT, dimension(nsize), intent(OUT) :: imin, imax
+  CCTK_REAL, dimension(nsize,nsize), intent(OUT) :: dd
+
+  CCTK_REAL, dimension(2), save :: a
+  CCTK_REAL, dimension(3,2), save :: q 
+
+  CCTK_INT :: i, il, ir
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    call coeffs_1_4 ( a, q )
+    first = .false.
+  end if
+
+  dd = zero
+  imin = 0
+  imax = -1
+
+  if ( bb(1) == 0 ) then
+    il = 1 + gsize
+  else
+    dd(1:3,1) = q(1:3,1)
+    imin(1) = 1; imax(1) = 1
+
+    dd(1:3,2) = q(1:3,2)
+    imin(2) = 1; imax(2) = 3
+
+    il = 3
+  end if
+  if ( bb(2) == 0 ) then
+    ir = nsize - gsize
+  else
+    dd(nsize-2:nsize,nsize-1) = -q(3:1:-1,2)
+    imin(nsize-1) = nsize-2; imax(nsize-1) = nsize
+
+    dd(nsize-2:nsize,nsize) = -q(3:1:-1,1)
+    imin(nsize) = nsize; imax(nsize) = nsize
+
+    ir = nsize - 2
+  end if
+  do i = il, ir
+    dd(i-2:i-1,i) = -a(2:1:-1);
+    dd(i+1:i+2,i) = a(1:2)
+    imin(i) = i-2; imax(i) = i+2
+  end do
+end subroutine set_coeff_4
diff --git a/src/Coefficients_6.F90 b/src/Coefficients_6.F90
new file mode 100644
index 0000000..a27089f
--- /dev/null
+++ b/src/Coefficients_6.F90
@@ -0,0 +1,69 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine set_coeff_6 ( nsize, loc_order, bb, gsize, imin, imax, dd )
+
+  use All_Coeffs_mod
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(IN) :: nsize, loc_order
+  CCTK_INT, dimension(2), intent(IN) :: bb
+  CCTK_INT, intent(IN) :: gsize
+  CCTK_INT, dimension(nsize), intent(OUT) :: imin, imax
+  CCTK_REAL, dimension(nsize,nsize), intent(OUT) :: dd
+
+  CCTK_REAL, dimension(3), save :: a
+  CCTK_REAL, dimension(5,3), save :: q 
+
+  CCTK_INT :: i, il, ir
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    call coeffs_1_6 ( a, q )
+    first = .false.
+  end if
+
+  dd = zero
+  imin = 0
+  imax = -1
+
+  if ( bb(1) == 0 ) then
+    il = 1 + gsize
+  else
+    dd(1:5,1) = q(1:5,1)
+    imin(1) = 1; imax(1) = 1
+
+    dd(1:5,2) = q(1:5,2)
+    imin(2) = 1; imax(2) = 3
+
+    dd(1:5,3) = q(1:5,3)
+    imin(3) = 1; imax(3) = 5
+
+    il = 4
+  end if
+  if ( bb(2) == 0 ) then
+    ir = nsize - gsize
+  else
+    dd(nsize-4:nsize,nsize-2) = -q(5:1:-1,3)
+    imin(nsize-2) = nsize-4; imax(nsize-2) = nsize
+
+    dd(nsize-4:nsize,nsize-1) = -q(5:1:-1,2)
+    imin(nsize-1) = nsize-2; imax(nsize-1) = nsize
+
+    dd(nsize-4:nsize,nsize) = -q(5:1:-1,1)
+    imin(nsize) = nsize; imax(nsize) = nsize
+
+    ir = nsize - 3
+  end if
+  do i = il, ir
+    dd(i-3:i-1,i) = -a(3:1:-1);
+    dd(i+1:i+3,i) = a(1:3)
+    imin(i) = i-3; imax(i) = i+3
+  end do
+end subroutine set_coeff_6
diff --git a/src/Coefficients_8.F90 b/src/Coefficients_8.F90
new file mode 100644
index 0000000..a9fdeb6
--- /dev/null
+++ b/src/Coefficients_8.F90
@@ -0,0 +1,75 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine set_coeff_8 ( nsize, loc_order, bb, gsize, imin, imax, dd )
+
+  use All_Coeffs_mod
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(IN) :: nsize, loc_order
+  CCTK_INT, dimension(2), intent(IN) :: bb
+  CCTK_INT, intent(IN) :: gsize
+  CCTK_INT, dimension(nsize), intent(OUT) :: imin, imax
+  CCTK_REAL, dimension(nsize,nsize), intent(OUT) :: dd
+
+  CCTK_REAL, dimension(4), save :: a
+  CCTK_REAL, dimension(7,4), save :: q 
+
+  CCTK_INT :: i, il, ir
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    call coeffs_1_8 ( a, q )
+    first = .false.
+  end if
+
+  dd = zero
+  imin = 0
+  imax = -1
+
+  if ( bb(1) == 0 ) then
+    il = 1 + gsize
+  else
+    dd(1:7,1) = q(1:7,1)
+    imin(1) = 1; imax(1) = 1
+
+    dd(1:7,2) = q(1:7,2)
+    imin(2) = 1; imax(2) = 3
+
+    dd(1:7,3) = q(1:7,3)
+    imin(3) = 1; imax(3) = 5
+
+    dd(1:7,4) = q(1:7,4)
+    imin(4) = 1; imax(4) = 7
+
+    il = 5
+  end if
+  if ( bb(2) == 0 ) then
+    ir = nsize - gsize
+  else
+    dd(nsize-6:nsize,nsize-3) = -q(7:1:-1,4)
+    imin(nsize-3) = nsize-6; imax(nsize-3) = nsize
+
+    dd(nsize-6:nsize,nsize-2) = -q(7:1:-1,3)
+    imin(nsize-2) = nsize-4; imax(nsize-2) = nsize
+
+    dd(nsize-6:nsize,nsize-1) = -q(7:1:-1,2)
+    imin(nsize-1) = nsize-2; imax(nsize-1) = nsize
+
+    dd(nsize-6:nsize,nsize) = -q(7:1:-1,1)
+    imin(nsize) = nsize; imax(nsize) = nsize
+
+    ir = nsize - 4
+  end if
+  do i = il, ir
+    dd(i-4:i-1,i) = -a(4:1:-1);
+    dd(i+1:i+4,i) = a(1:4)
+    imin(i) = i-4; imax(i) = i+4
+  end do
+end subroutine set_coeff_8
diff --git a/src/get_coeffs.c b/src/get_coeffs.c
index 836ff52..31b3601 100644
--- a/src/get_coeffs.c
+++ b/src/get_coeffs.c
@@ -32,6 +32,13 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
                                             CCTK_INT *imin,
                                             CCTK_INT *imax,
                                             CCTK_REAL *q);
+  void CCTK_FCALL CCTK_FNAME(set_coeff_2)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
   void CCTK_FCALL CCTK_FNAME(set_coeff_4_2)(const CCTK_INT *nsize,
                                             const CCTK_INT *loc_order,
                                             const CCTK_INT *bb,
@@ -39,6 +46,13 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
                                             CCTK_INT *imin,
                                             CCTK_INT *imax,
                                             CCTK_REAL *q);
+  void CCTK_FCALL CCTK_FNAME(set_coeff_4)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
   void CCTK_FCALL CCTK_FNAME(set_coeff_4_3)(const CCTK_INT *nsize,
                                             const CCTK_INT *loc_order,
                                             const CCTK_INT *bb,
@@ -67,6 +81,13 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
                                             CCTK_INT *imin,
                                             CCTK_INT *imax,
                                             CCTK_REAL *q);
+  void CCTK_FCALL CCTK_FNAME(set_coeff_6)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
   void CCTK_FCALL CCTK_FNAME(set_coeff_6_5)(const CCTK_INT *nsize,
                                             const CCTK_INT *loc_order,
                                             const CCTK_INT *bb,
@@ -95,6 +116,13 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
                                             CCTK_INT *imin,
                                             CCTK_INT *imax,
                                             CCTK_REAL *q);
+  void CCTK_FCALL CCTK_FNAME(set_coeff_8)(const CCTK_INT *nsize,
+                                          const CCTK_INT *loc_order,
+                                          const CCTK_INT *bb,
+                                          const CCTK_INT *gsize,
+                                          CCTK_INT *imin,
+                                          CCTK_INT *imax,
+                                          CCTK_REAL *q);
 
   ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2];
 
@@ -135,47 +163,70 @@ void DiffCoeff ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
   }
 
   if ( CCTK_Equals(norm_type,"Diagonal") ) {
-    if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) {
-      switch(loc_order) {
-      case 2: {
-        CCTK_FNAME(set_coeff_2_1)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
-        break;
-      }
-      case 4: {
-        CCTK_FNAME(set_coeff_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
-        break;
-      }
-      case 6: {
-        CCTK_FNAME(set_coeff_6_3)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
-        break;
-      }
-      case 8: {
-        CCTK_FNAME(set_coeff_8_4)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
-        break;
-      }
-      default:
-        CCTK_WARN (0, "Unknown stencil specified");
+    if ( sbp_1st_deriv ) {
+      if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) {
+        switch(loc_order) {
+        case 2: {
+          CCTK_FNAME(set_coeff_2_1)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+          break;
+        }
+        case 4: {
+          CCTK_FNAME(set_coeff_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+          break;
+        }
+        case 6: {
+          CCTK_FNAME(set_coeff_6_3)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+          break;
+        }
+        case 8: {
+          CCTK_FNAME(set_coeff_8_4)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+          break;
+        }
+        default:
+          CCTK_WARN (0, "Unknown stencil specified");
+        }
+      } else {
+        switch(loc_order) {
+        case 2: {
+          CCTK_FNAME(set_coeff_2_1)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+          break;
+        }
+        case 4: {
+          CCTK_FNAME(set_coeff_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+          break;
+        }
+        case 6: {
+          CCTK_FNAME(set_coeff_6_3_opt)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+          break;
+        }
+        case 8: {
+          CCTK_FNAME(set_coeff_8_4_opt)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+          break;
+        }
+        default:
+          CCTK_WARN (0, "Unknown stencil specified");
+        }
       }
     } else {
-      switch(loc_order) {
+     switch(loc_order) {
       case 2: {
-        CCTK_FNAME(set_coeff_2_1)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+        CCTK_FNAME(set_coeff_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
         break;
       }
       case 4: {
-        CCTK_FNAME(set_coeff_4_2)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+        CCTK_FNAME(set_coeff_4)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
         break;
       }
       case 6: {
-        CCTK_FNAME(set_coeff_6_3_opt)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+        CCTK_FNAME(set_coeff_6)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
         break;
       }
       case 8: {
-        CCTK_FNAME(set_coeff_8_4_opt)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
+        CCTK_FNAME(set_coeff_8)(&nsize,&loc_order,bb,&gsize,imin,imax,q);
         break;
       }
       default:
-        CCTK_WARN (0, "Unknown stencil specified");
+        CCTK_WARN (0, "Unknown 1st derivative stencil specified");
       }
     }
   } else {
diff --git a/src/make.code.defn b/src/make.code.defn
index 8eaaf16..7fa5829 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -46,19 +46,23 @@ SRCS =  call_derivs.c \
         get_up_coeffs.c \
         get_coeffs2.c \
         Coefficients_2_1.F90 \
+        Coefficients_2.F90 \
         Coefficients2_2_1.F90 \
         Coefficients2_2.F90 \
         Coefficients_4_2.F90 \
+        Coefficients_4.F90 \
         Coefficients2_4_2.F90 \
         Coefficients2_4_2_min_err_coeff.F90 \
         Coefficients2_4.F90 \
         Coefficients_6_3.F90 \
         Coefficients2_6_3.F90 \
         Coefficients_6_3_min_err_coeff.F90 \
+        Coefficients_6.F90 \
         Coefficients2_6.F90 \
         Coefficients_8_4.F90 \
         Coefficients2_8_4.F90 \
         Coefficients_8_4_min_err_coeff.F90 \
+        Coefficients_8.F90 \
         Coefficients2_8.F90 \
         Coefficients_4_3.F90 \
         Coefficients_4_3_min_err_coeff.F90 \