First attempt at a summation by parts finite differencing thorn.

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

1 files changed, 239 insertions, 0 deletions

diff --git a/src/Derivatives_6_3.F90 b/src/Derivatives_6_3.F90
new file mode 100644
index 0000000..3da882d
--- /dev/null
+++ b/src/Derivatives_6_3.F90
@@ -0,0 +1,239 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine deriv_gf_6_3 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar )
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL, parameter :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  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) :: dvar
+
+  CCTK_REAL, dimension(3), save :: a
+  CCTK_REAL, dimension(9,6), save :: q 
+  CCTK_REAL :: idel
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    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) = -21600.0_wp/13649.0_wp; q(2,1) = 81763.0_wp/40947.0_wp
+    q(3,1) = 131.0_wp/27298.0_wp; q(4,1) = -9143.0_wp/13649.0_wp
+    q(5,1) = 20539.0_wp/81894.0_wp; q(6,1) = zero;
+    q(7,1) = zero; q(8,1) = zero; q(9,1) = zero
+    q(1,2) = -81763.0_wp/180195.0_wp; q(2,2) = zero
+    q(3,2) = 7357.0_wp/36039.0_wp; q(4,2) = 30637.0_wp/72078.0_wp
+    q(5,2) = -2328.0_wp/12013.0_wp; q(6,2) = 6611.0_wp/360390.0_wp
+    q(7,2) = zero; q(8,2) = zero; q(9,2) = zero
+    q(1,3) = -131.0_wp/54220.0_wp; q(2,3) = -7357.0_wp/16266.0_wp
+    q(3,3) = zero; q(4,3) = 645.0_wp/2711.0_wp
+    q(5,3) = 11237.0_wp/32532.0_wp; q(6,3) = -3487.0_wp/27110.0_wp
+    q(7,3) = zero; q(8,3) = zero; q(9,3) = zero
+    q(1,4) = 9143.0_wp/53590.0_wp; q(2,4) = -30637.0_wp/64308.0_wp
+    q(3,4) =  -645.0_wp/5359.0_wp; q(4,4) = zero
+    q(5,4) = 13733.0_wp/32154.0_wp; q(6,4) = -67.0_wp/4660.0_wp
+    q(7,4) = 72.0_wp/5359.0_wp; q(8,4) = zero; q(9,4) = zero;
+    q(1,5) = -20539.0_wp/236310.0_wp; q(2,5) = 2328.0_wp/7877.0_wp
+    q(3,5) = -11237.0_wp/47262.0_wp; q(4,5) = -13733.0_wp/23631.0_wp
+    q(5,5) = zero; q(6,5) = 89387.0_wp/118155.0_wp
+    q(7,5) = -1296.0_wp/7877.0_wp; q(8,5) = 144.0_wp/7877.0_wp; q(9,5) = zero
+    q(1,6) = zero; q(2,6) = -6611.0_wp/262806.0_wp
+    q(3,6) = 3487.0_wp/43801.0_wp; q(4,6) = 1541.0_wp/87602.0_wp
+    q(5,6) = -89387.0_wp/131403.0_wp; q(6,6) = zero
+    q(7,6) = 32400.0_wp/43801.0_wp; q(8,6) = -6480.0_wp/43801.0_wp
+    q(9,6) = 720.0_wp/43801.0_wp
+    first = .false. 
+  end if
+
+  idel = 1.0_wp / delta
+
+  direction: select case (dir)
+  case (0) direction
+    if ( bb(1) == 0 ) then
+      il = 1 + gsize
+    else
+      dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + &
+                      q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) + &
+                      q(5,1) * var(5,:,:) ) * idel
+      dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(3,2) * var(3,:,:) + &
+                      q(4,2) * var(4,:,:) + q(5,2) * var(5,:,:) + &
+                      q(6,2) * var(6,:,:) ) * idel
+      dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + &
+                      q(4,3) * var(4,:,:) + q(5,3) * var(5,:,:) + &
+                      q(6,3) * var(6,:,:) ) * idel
+      dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + &
+                      q(3,4) * var(3,:,:) + q(5,4) * var(5,:,:) + &
+                      q(6,4) * var(6,:,:) + q(7,4) * var(7,:,:) ) * idel
+      dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + &
+                      q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + &
+                      q(6,5) * var(6,:,:) + q(7,5) * var(7,:,:) + &
+                      q(8,5) * var(8,:,:) ) * idel
+      dvar(6,:,:) = ( q(2,6) * var(2,:,:) + q(3,6) * var(3,:,:) + &
+                      q(4,6) * var(4,:,:) + q(5,6) * var(5,:,:) + &
+                      q(7,6) * var(7,:,:) + q(8,6) * var(8,:,:) + &
+                      q(9,6) * var(9,:,:) ) * idel
+      il = 7
+    end if
+    if ( bb(2) == 0 ) then
+      ir = ni - gsize
+    else
+      dvar(ni,:,:) =   - ( q(1,1) * var(ni,:,:) + q(2,1) * var(ni-1,:,:) + &
+                           q(3,1) * var(ni-2,:,:) + q(4,1) * var(ni-3,:,:) + &
+                           q(5,1) * var(ni-4,:,:) ) * idel
+      dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(3,2) * var(ni-2,:,:) + &
+                           q(4,2) * var(ni-3,:,:) + q(5,2) * var(ni-4,:,:) + &
+                           q(6,2) * var(ni-5,:,:) ) * idel
+      dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + &
+                           q(4,3) * var(ni-3,:,:) + q(5,3) * var(ni-4,:,:) + &
+                           q(6,3) * var(ni-5,:,:) ) * idel
+      dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + &
+                           q(3,4) * var(ni-2,:,:) + q(5,4) * var(ni-4,:,:) + &
+                           q(6,4) * var(ni-5,:,:) + &
+                           q(7,4) * var(ni-6,:,:) ) * idel
+      dvar(ni-4,:,:) = - ( q(1,5) * var(ni,:,:) + q(2,5) * var(ni-1,:,:) + &
+                           q(3,5) * var(ni-2,:,:) + q(4,5) * var(ni-3,:,:) + &
+                           q(6,5) * var(ni-5,:,:) + q(7,5) * var(ni-6,:,:) + &
+                           q(8,5) * var(ni-7,:,:) ) * idel
+      dvar(ni-5,:,:) = - ( q(2,6) * var(ni-1,:,:) + q(3,6) * var(ni-2,:,:) + &
+                           q(4,6) * var(ni-3,:,:) + q(5,6) * var(ni-4,:,:) + &
+                           q(7,6) * var(ni-6,:,:) + q(8,6) * var(ni-7,:,:) + &
+                           q(9,6) * var(ni-8,:,:) ) * idel
+      ir = ni - 6
+    end if
+    dvar(il:ir,:,:) = ( a(1) * ( var(il+1:ir+1,:,:) - &
+                                 var(il-1:ir-1,:,:) ) + &
+                        a(2) * ( var(il+2:ir+2,:,:) - &
+                                 var(il-2:ir-2,:,:) ) + &
+                        a(3) * ( var(il+3:ir+3,:,:) - &
+                                 var(il-3:ir-3,:,:) ) ) * idel
+  case (1) direction
+    if ( bb(1) == 0 ) then
+      jl = 1 + gsize
+    else
+      dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + &
+                      q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) + &
+                      q(5,1) * var(:,5,:) ) * idel
+      dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(3,2) * var(:,3,:) + &
+                      q(4,2) * var(:,4,:) + q(5,2) * var(:,5,:) + &
+                      q(6,2) * var(:,6,:) ) * idel
+      dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + &
+                      q(4,3) * var(:,4,:) + q(5,3) * var(:,5,:) + &
+                      q(6,3) * var(:,6,:) ) * idel
+      dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + &
+                      q(3,4) * var(:,3,:) + q(5,4) * var(:,5,:) + &
+                      q(6,4) * var(:,6,:) + q(7,4) * var(:,7,:) ) * idel
+      dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + &
+                      q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + &
+                      q(6,5) * var(:,6,:) + q(7,5) * var(:,7,:) + &
+                      q(8,5) * var(:,8,:) ) * idel
+      dvar(:,6,:) = ( q(2,6) * var(:,2,:) + q(3,6) * var(:,3,:) + &
+                      q(4,6) * var(:,4,:) + q(5,6) * var(:,5,:) + &
+                      q(7,6) * var(:,7,:) + q(8,6) * var(:,8,:) + &
+                      q(9,6) * var(:,9,:) ) * idel
+      jl = 7
+    end if
+    if ( bb(2) == 0 ) then
+      jr = nj - gsize
+    else
+      dvar(:,nj,:) =   - ( q(1,1) * var(:,nj,:) + q(2,1) * var(:,nj-1,:) + &
+                           q(3,1) * var(:,nj-2,:) + q(4,1) * var(:,nj-3,:) + &
+                           q(5,1) * var(:,nj-4,:) ) * idel
+      dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(3,2) * var(:,nj-2,:) + &
+                           q(4,2) * var(:,nj-3,:) + q(5,2) * var(:,nj-4,:) + &
+                           q(6,2) * var(:,nj-5,:) ) * idel
+      dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + &
+                           q(4,3) * var(:,nj-3,:) + q(5,3) * var(:,nj-4,:) + &
+                           q(6,3) * var(:,nj-5,:) ) * idel
+      dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + &
+                           q(3,4) * var(:,nj-2,:) + q(5,4) * var(:,nj-4,:) + &
+                           q(6,4) * var(:,nj-5,:) + &
+                           q(7,4) * var(:,nj-6,:) ) * idel
+      dvar(:,nj-4,:) = - ( q(1,5) * var(:,nj,:) + q(2,5) * var(:,nj-1,:) + &
+                           q(3,5) * var(:,nj-2,:) + q(4,5) * var(:,nj-3,:) + &
+                           q(6,5) * var(:,nj-5,:) + q(7,5) * var(:,nj-6,:) + &
+                           q(8,5) * var(:,nj-7,:) ) * idel
+      dvar(:,nj-5,:) = - ( q(2,6) * var(:,nj-1,:) + q(3,6) * var(:,nj-2,:) + &
+                           q(4,6) * var(:,nj-3,:) + q(5,6) * var(:,nj-4,:) + &
+                           q(7,6) * var(:,nj-6,:) + q(8,6) * var(:,nj-7,:) + &
+                           q(9,6) * var(:,nj-8,:) ) * idel
+      jr = nj - 6
+    end if
+    dvar(:,jl:jr,:) = ( a(1) * ( var(:,jl+1:jr+1,:) - &
+                                 var(:,jl-1:jr-1,:) ) + &
+                        a(2) * ( var(:,jl+2:jr+2,:) - &
+                                 var(:,jl-2:jr-2,:) ) + &
+                        a(3) * ( var(:,jl+3:jr+3,:) - &
+                                 var(:,jl-3:jr-3,:) ) ) * idel
+  case (2) direction
+    if ( bb(1) == 0 ) then
+      kl = 1 + gsize
+    else
+      dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + &
+                      q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) + &
+                      q(5,1) * var(:,:,5) ) * idel
+      dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(3,2) * var(:,:,3) + &
+                      q(4,2) * var(:,:,4) + q(5,2) * var(:,:,5) + &
+                      q(6,2) * var(:,:,6) ) * idel
+      dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + &
+                      q(4,3) * var(:,:,4) + q(5,3) * var(:,:,5) + &
+                      q(6,3) * var(:,:,6) ) * idel
+      dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + &
+                      q(3,4) * var(:,:,3) + q(5,4) * var(:,:,5) + &
+                      q(6,4) * var(:,:,6) + q(7,4) * var(:,:,7) ) * idel
+      dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + &
+                      q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + &
+                      q(6,5) * var(:,:,6) + q(7,5) * var(:,:,7) + &
+                      q(8,5) * var(:,:,8) ) * idel
+      dvar(:,:,6) = ( q(2,6) * var(:,:,2) + q(3,6) * var(:,:,3) + &
+                      q(4,6) * var(:,:,4) + q(5,6) * var(:,:,5) + &
+                      q(7,6) * var(:,:,7) + q(8,6) * var(:,:,8) + &
+                      q(9,6) * var(:,:,9) ) * idel
+      kl = 7
+    end if
+    if ( bb(2) == 0 ) then
+      kr = nk - gsize
+    else 
+      dvar(:,:,nk) =   - ( q(1,1) * var(:,:,nk) + q(2,1) * var(:,:,nk-1) + &
+                           q(3,1) * var(:,:,nk-2) + q(4,1) * var(:,:,nk-3) + &
+                           q(5,1) * var(:,:,nk-4) ) * idel
+      dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(3,2) * var(:,:,nk-2) + &
+                           q(4,2) * var(:,:,nk-3) + q(5,2) * var(:,:,nk-4) + &
+                           q(6,2) * var(:,:,nk-5) ) * idel
+      dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + &
+                           q(4,3) * var(:,:,nk-3) + q(5,3) * var(:,:,nk-4) + &
+                           q(6,3) * var(:,:,nk-5) ) * idel
+      dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + &
+                           q(3,4) * var(:,:,nk-2) + q(5,4) * var(:,:,nk-4) + &
+                           q(6,4) * var(:,:,nk-5) + &
+                           q(7,4) * var(:,:,nk-6) ) * idel
+      dvar(:,:,nk-4) = - ( q(1,5) * var(:,:,nk) + q(2,5) * var(:,:,nk-1) + &
+                           q(3,5) * var(:,:,nk-2) + q(4,5) * var(:,:,nk-3) + &
+                           q(6,5) * var(:,:,nk-5) + q(7,5) * var(:,:,nk-6) + &
+                           q(8,5) * var(:,:,nk-7) ) * idel
+      dvar(:,:,nk-5) = - ( q(2,6) * var(:,:,nk-1) + q(3,6) * var(:,:,nk-2) + &
+                           q(4,6) * var(:,:,nk-3) + q(5,6) * var(:,:,nk-4) + &
+                           q(7,6) * var(:,:,nk-6) + q(8,6) * var(:,:,nk-7) + &
+                           q(9,6) * var(:,:,nk-8) ) * idel
+      kr = nk - 6
+    end if
+    dvar(:,:,kl:kr) = ( a(1) * ( var(:,:,kl+1:kr+1) - &
+                                 var(:,:,kl-1:kr-1) ) + &
+                        a(2) * ( var(:,:,kl+2:kr+2) - &
+                                 var(:,:,kl-2:kr-2) ) + &
+                        a(3) * ( var(:,:,kl+3:kr+3) - &
+                                 var(:,:,kl-3:kr-3) ) ) * idel
+  end select direction
+end subroutine deriv_gf_6_3