The second derivative in the 4-2 case is not unique, so I added the minimal

bandwidth operator in addition to the ABTE operator already there.


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

6 files changed, 257 insertions, 15 deletions

diff --git a/src/All_Coeffs_mod.F90 b/src/All_Coeffs_mod.F90
index ef965d8..76dfcc9 100644
--- a/src/All_Coeffs_mod.F90
+++ b/src/All_Coeffs_mod.F90
@@ -194,6 +194,45 @@ module All_Coeffs_mod
       CCTK_REAL, dimension(3), intent(OUT) :: a
       CCTK_REAL, dimension(6,4), intent(OUT) :: q
 
+      a(1) = -2.500000000000000000000000000000000000000_wp
+      a(2) = 1.333333333333333333333333333333333333333_wp
+      a(3) = -0.08333333333333333333333333333333333333333_wp
+
+      q(1,1) = 2.000000000000000000000000000000000000000_wp
+      q(2,1) = -5.000000000000000000000000000000000000000_wp
+      q(3,1) = 4.000000000000000000000000000000000000000_wp
+      q(4,1) = -1.000000000000000000000000000000000000000_wp
+      q(5,1) = 0_wp
+      q(6,1) = 0_wp
+
+      q(1,2) = 1.000000000000000000000000000000000000000_wp
+      q(2,2) = -2.000000000000000000000000000000000000000_wp
+      q(3,2) = 1.000000000000000000000000000000000000000_wp
+      q(4,2) = 0_wp
+      q(5,2) = 0_wp
+      q(6,2) = 0_wp
+
+      q(1,3) = -0.09302325581395348837209302325581395348837_wp
+      q(2,3) = 1.372093023255813953488372093023255813953_wp
+      q(3,3) = -2.558139534883720930232558139534883720930_wp
+      q(4,3) = 1.372093023255813953488372093023255813953_wp
+      q(5,3) = -0.09302325581395348837209302325581395348837_wp
+      q(6,3) = 0_wp
+
+      q(1,4) = -0.02040816326530612244897959183673469387755_wp
+      q(2,4) = 0_wp
+      q(3,4) = 1.204081632653061224489795918367346938776_wp
+      q(4,4) = -2.408163265306122448979591836734693877551_wp
+      q(5,4) = 1.306122448979591836734693877551020408163_wp
+      q(6,4) = -0.08163265306122448979591836734693877551020_wp
+
+    end subroutine coeffs_2_4_2
+
+    subroutine coeffs_2_4_2_opt ( a, q )
+
+      CCTK_REAL, dimension(3), intent(OUT) :: a
+      CCTK_REAL, dimension(6,4), intent(OUT) :: q
+
       a(1) = -2.5_wp
       a(2) = 1.333333333333333333333333333333333333333_wp
       a(3) = -0.08333333333333333333333333333333333333333_wp
@@ -226,7 +265,7 @@ module All_Coeffs_mod
       q(5,4) = -0.9258553965216912759079296546576971585129_wp
       q(6,4) = 0.4150007042443253806547133319378374518503_wp
 
-    end subroutine coeffs_2_4_2
+    end subroutine coeffs_2_4_2_opt
 
     subroutine coeffs_1_6_3 ( a, q )
 
diff --git a/src/Derivatives2_4_2.F90 b/src/Derivatives2_4_2.F90
index 1827966..443775f 100644
--- a/src/Derivatives2_4_2.F90
+++ b/src/Derivatives2_4_2.F90
@@ -179,10 +179,10 @@ subroutine deriv2_gf_4_2 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar2 )
                             q(5,3) * var(:,:,nk-4) + &
                             q(6,3) * var(:,:,nk-5) ) * idel2
         dvar2(:,:,nk-1) = ( q(1,2) * var(:,:,nk) + &
-                            q(2,2) * var(:,:,nk-1) * &
-                            q(3,2) * var(:,:,nk-2) * &
-                            q(4,2) * var(:,:,nk-3) * &
-                            q(5,2) * var(:,:,nk-4) * &
+                            q(2,2) * var(:,:,nk-1) + &
+                            q(3,2) * var(:,:,nk-2) + &
+                            q(4,2) * var(:,:,nk-3) + &
+                            q(5,2) * var(:,:,nk-4) + &
                             q(6,2) * var(:,:,nk-5) ) * idel2
         dvar2(:,:,nk)   = ( q(1,1) * var(:,:,nk) + &
                             q(2,1) * var(:,:,nk-1) + &
diff --git a/src/Derivatives2_4_2_min_err_coeff.F90 b/src/Derivatives2_4_2_min_err_coeff.F90
new file mode 100644
index 0000000..ed3d81f
--- /dev/null
+++ b/src/Derivatives2_4_2_min_err_coeff.F90
@@ -0,0 +1,177 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine deriv2_gf_4_2_opt ( var, ni, nj, nk, dir, bb, gsize, delta, dvar2 )
+
+  use All_Coeffs_mod
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(IN) :: ni, nj, nk
+  CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var
+  CCTK_INT, intent(IN) :: dir
+  CCTK_INT, intent(IN) :: bb(2)
+  CCTK_INT, intent(IN) :: gsize
+  CCTK_REAL, intent(IN) :: delta
+  CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar2
+
+  CCTK_REAL, dimension(3), save :: a
+  CCTK_REAL, dimension(6,4), save :: q 
+  CCTK_REAL :: idel2
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    call coeffs_2_4_2_opt ( a, q )
+    first = .false.
+  end if
+
+  idel2 = 1.0_wp / delta**2
+
+  direction: select case (dir)
+  case (0) direction
+    if ( bb(1) == 0 ) then
+      il = 1 + gsize
+    else
+      dvar2(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:)  + &
+                       q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) ) * idel2
+      dvar2(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + &
+                       q(3,2) * var(3,:,:) ) * idel2
+      dvar2(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + &
+                       q(3,3) * var(3,:,:) + q(4,3) * var(4,:,:) + &
+                       q(5,3) * var(5,:,:) ) * idel2
+      dvar2(4,:,:) = ( q(1,4) * var(1,:,:) + q(3,4) * var(3,:,:) + &
+                       q(4,4) * var(4,:,:) + q(5,4) * var(5,:,:) + &
+                       q(6,4) * var(6,:,:) ) * idel2
+      il = 5
+    end if
+    if ( bb(2) == 0 ) then
+      ir = ni - gsize
+    else
+      dvar2(ni-3,:,:) = ( q(1,4) * var(ni,:,:) + &
+                          q(3,4) * var(ni-2,:,:) + &
+                          q(4,4) * var(ni-3,:,:) + &
+                          q(5,4) * var(ni-4,:,:) + &
+                          q(6,4) * var(ni-5,:,:) ) * idel2
+      dvar2(ni-2,:,:) = ( q(1,3) * var(ni,:,:) + &
+                          q(2,3) * var(ni-1,:,:) + &
+                          q(3,3) * var(ni-2,:,:) + &
+                          q(4,3) * var(ni-3,:,:) + &
+                          q(5,3) * var(ni-4,:,:) ) * idel2
+      dvar2(ni-1,:,:) = ( q(1,2) * var(ni,:,:) + &
+                          q(2,2) * var(ni-1,:,:) + &
+                          q(3,2) * var(ni-2,:,:) ) * idel2
+      dvar2(ni,:,:)   = ( q(1,1) * var(ni,:,:) + &
+                          q(2,1) * var(ni-1,:,:) + &
+                          q(3,1) * var(ni-2,:,:) + &
+                          q(4,1) * var(ni-3,:,:) ) * idel2
+      ir = ni - 4
+    end if
+    dvar2(il:ir,:,:) = ( a(1) * var(il:ir,:,:) + &
+                         a(2) * ( var(il+1:ir+1,:,:) + &
+                                  var(il-1:ir-1,:,:) ) + &
+                         a(3) * ( var(il+2:ir+2,:,:) + &
+                                  var(il-2:ir-2,:,:) ) ) * idel2
+  case (1) direction
+    if ( zero_derivs_y /= 0 ) then
+      dvar2 = zero
+    else
+      if ( bb(1) == 0 ) then
+        jl = 1 + gsize
+      else
+        dvar2(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + &
+                         q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) ) * idel2
+        dvar2(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + &
+                         q(3,2) * var(:,3,:) ) * idel2
+        dvar2(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + &
+                         q(3,3) * var(:,3,:) + q(4,3) * var(:,4,:) + &
+                         q(5,3) * var(:,5,:) ) * idel2
+        dvar2(:,4,:) = ( q(1,4) * var(:,1,:) + q(3,4) * var(:,3,:) + &
+                         q(4,4) * var(:,4,:) + q(5,4) * var(:,5,:) + &
+                         q(6,4) * var(:,6,:) ) * idel2
+        jl = 5
+      end if
+      if ( bb(2) == 0 ) then
+        jr = nj - gsize
+      else
+        dvar2(:,nj-3,:) = ( q(1,4) * var(:,nj,:) + &
+                            q(3,4) * var(:,nj-2,:) + &
+                            q(4,4) * var(:,nj-3,:) + &
+                            q(5,4) * var(:,nj-4,:) + &
+                            q(6,4) * var(:,nj-5,:) ) * idel2
+        dvar2(:,nj-2,:) = ( q(1,3) * var(:,nj,:) + &
+                            q(2,3) * var(:,nj-1,:) + &
+                            q(3,3) * var(:,nj-2,:) + &
+                            q(4,3) * var(:,nj-3,:) + &
+                            q(5,3) * var(:,nj-4,:) ) * idel2
+        dvar2(:,nj-1,:) = ( q(1,2) * var(:,nj,:) + &
+                            q(2,2) * var(:,nj-1,:) + &
+                            q(3,2) * var(:,nj-2,:) ) * idel2
+        dvar2(:,nj,:)   = ( q(1,1) * var(:,nj,:) + &
+                            q(2,1) * var(:,nj-1,:) + &
+                            q(3,1) * var(:,nj-2,:) + &
+                            q(4,1) * var(:,nj-3,:) ) * idel2
+        jr = nj - 4
+      end if
+      dvar2(:,jl:jr,:) = ( a(1) * var(:,jl:jr,:) + &
+                           a(2) * ( var(:,jl+1:jr+1,:) + &
+                                    var(:,jl-1:jr-1,:) ) + &
+                           a(3) * ( var(:,jl+2:jr+2,:) + &
+                                    var(:,jl-2:jr-2,:) ) ) * idel2
+
+    end if
+  case (2) direction
+    if ( zero_derivs_z /= 0 ) then
+      dvar2 = zero
+    else
+      if ( bb(1) == 0 ) then
+        kl = 1 + gsize
+      else
+        dvar2(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + &
+                         q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) ) * idel2
+        dvar2(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + &
+                         q(3,2) * var(:,:,3) ) * idel2
+        dvar2(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + &
+                         q(3,3) * var(:,:,3) + q(4,3) * var(:,:,4) + &
+                         q(5,3) * var(:,:,5) ) * idel2
+        dvar2(:,:,4) = ( q(1,4) * var(:,:,1) + q(3,4) * var(:,:,3) + &
+                         q(4,4) * var(:,:,4) + q(5,4) * var(:,:,5) + &
+                         q(6,4) * var(:,:,6) ) * idel2
+        kl = 5
+      end if
+      if ( bb(2) == 0 ) then
+        kr = nk - gsize
+      else 
+        dvar2(:,:,nk-3) = ( q(1,4) * var(:,:,nk) + &
+                            q(3,4) * var(:,:,nk-2) + &
+                            q(4,4) * var(:,:,nk-3) + &
+                            q(5,4) * var(:,:,nk-4) + &
+                            q(6,4) * var(:,:,nk-5) ) * idel2
+        dvar2(:,:,nk-2) = ( q(1,3) * var(:,:,nk) + &
+                            q(2,3) * var(:,:,nk-1) + &
+                            q(3,3) * var(:,:,nk-2) + &
+                            q(4,3) * var(:,:,nk-3) + &
+                            q(5,3) * var(:,:,nk-4) ) * idel2
+        dvar2(:,:,nk-1) = ( q(1,2) * var(:,:,nk) + &
+                            q(2,2) * var(:,:,nk-1) + &
+                            q(3,2) * var(:,:,nk-2) ) * idel2
+        dvar2(:,:,nk)   = ( q(1,1) * var(:,:,nk) + &
+                            q(2,1) * var(:,:,nk-1) + &
+                            q(3,1) * var(:,:,nk-2) + &
+                            q(4,1) * var(:,:,nk-3) ) * idel2
+        kr = nk - 4
+      end if
+      dvar2(:,:,kl:kr) = ( a(1) * var(:,:,kl:kr) + &
+                           a(2) * ( var(:,:,kl+1:kr+1) + &
+                                    var(:,:,kl-1:kr-1) ) + &
+                           a(3) * ( var(:,:,kl+2:kr+2) + &
+                                    var(:,:,kl-2:kr-2) ) ) * idel2
+    end if
+  end select direction
+end subroutine deriv2_gf_4_2_opt
diff --git a/src/call_derivs2.c b/src/call_derivs2.c
index 23a5899..cf193c6 100644
--- a/src/call_derivs2.c
+++ b/src/call_derivs2.c
@@ -41,6 +41,15 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
                                        const CCTK_INT *gsize,
                                        const CCTK_REAL *delta,
                                        CCTK_REAL *dvar2);
+  void CCTK_FCALL CCTK_FNAME(deriv2_gf_4_2_opt)(const CCTK_REAL *var,
+                                       const CCTK_INT *ni,
+                                       const CCTK_INT *nj,
+                                       const CCTK_INT *nk,
+                                       const CCTK_INT *dir,
+                                       const CCTK_INT *bb,
+                                       const CCTK_INT *gsize,
+                                       const CCTK_REAL *delta,
+                                       CCTK_REAL *dvar2);
 
   
   ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2];
@@ -85,17 +94,32 @@ void DiffGv2 ( const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT dir,
   }
 
   if ( CCTK_Equals(norm_type,"Diagonal") ) {
-    switch(loc_order) {
-    case 2: {
-      CCTK_FNAME(deriv2_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2);
-      break;
-    }
-    case 4: {
-      CCTK_FNAME(deriv2_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2);
-      break;
-    }
-    default:
+    if ( CCTK_Equals(operator_type,"Minimal Bandwidth") ) {
+      switch(loc_order) {
+      case 2: {
+        CCTK_FNAME(deriv2_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2);
+        break;
+      }
+      case 4: {
+        CCTK_FNAME(deriv2_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2);
+        break;
+      }
+      default:
       CCTK_WARN (0, "Unknown 2nd derivative stencil specified");
+      }
+    } else {
+      switch(loc_order) {
+      case 2: {
+        CCTK_FNAME(deriv2_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2);
+        break;
+      }
+      case 4: {
+        CCTK_FNAME(deriv2_gf_4_2_opt)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar2);
+        break;
+      }
+      default:
+        CCTK_WARN (0, "Unknown stencil specified");
+      }
     }
   } else {
     switch(loc_order) {
diff --git a/src/make.code.defn b/src/make.code.defn
index b5bb131..4357264 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -11,6 +11,7 @@ SRCS =  call_derivs.c \
         Derivatives2_2_1.F90 \
         Derivatives_4_2.F90 \
         Derivatives2_4_2.F90 \
+        Derivatives2_4_2_min_err_coeff.F90 \
         Derivatives_6_3.F90 \
         Derivatives_6_3_min_err_coeff.F90 \
         Derivatives_8_4.F90 \
diff --git a/src/make.code.deps b/src/make.code.deps
index 23655a1..6091656 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -5,6 +5,7 @@ Coefficients2_2_1.o:	All_Coeffs_mod.F90.o
 Derivatives_4_2.o:	All_Coeffs_mod.F90.o
 Coefficients_4_2.o:	All_Coeffs_mod.F90.o
 Derivatives2_4_2.o:	All_Coeffs_mod.F90.o
+Derivatives2_4_2_min_err_coeff.o:	All_Coeffs_mod.F90.o
 Coefficients2_4_2.o:	All_Coeffs_mod.F90.o
 Derivatives_4_3.o:	All_Coeffs_mod.F90.o
 Coefficients_4_3.o:	All_Coeffs_mod.F90.o