Kreiss-Oliger type dissipation operators compatible with the diagonal 4-2 and

6-3 derivative operators.


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

1 files changed, 240 insertions, 0 deletions

diff --git a/src/Dissipation_4_2_alt.F90 b/src/Dissipation_4_2_alt.F90
new file mode 100644
index 0000000..c0a68ba
--- /dev/null
+++ b/src/Dissipation_4_2_alt.F90
@@ -0,0 +1,240 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+subroutine dissipation_4_2_alt (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) :: epsilon 
+
+  CCTK_REAL :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  integer :: i, j, k
+  CCTK_REAL, dimension(7,4) :: a
+  CCTK_REAL :: idel
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 64 * delta(1) )
+
+  if ( bb(1) == 0 ) then
+    il = 1 + gsize(1)
+  else
+    rhs(1,:,:) = rhs(1,:,:) + &
+                 ( a(1,1) * var(1,:,:) + a(2,1) * var(2,:,:) + &
+                   a(3,1) * var(3,:,:) + a(4,1) * var(4,:,:) ) * idel
+    rhs(2,:,:) = rhs(2,:,:) + &
+                 ( a(1,2) * var(1,:,:) + a(2,2) * var(2,:,:) + &
+                   a(3,2) * var(3,:,:) + a(4,2) * var(4,:,:) + &
+                   a(5,2) * var(5,:,:) ) * idel
+    rhs(3,:,:) = rhs(3,:,:) + &
+                 ( a(1,3) * var(1,:,:) + a(2,3) * var(2,:,:) + &
+                   a(3,3) * var(3,:,:) + a(4,3) * var(4,:,:) + &
+                   a(5,3) * var(5,:,:) + a(6,3) * var(6,:,:) ) * idel
+    rhs(4,:,:) = rhs(4,:,:) + &
+                 ( a(1,4) * var(1,:,:) + a(2,4) * var(2,:,:) + &
+                   a(3,4) * var(3,:,:) + a(4,4) * var(4,:,:) + &
+                   a(5,4) * var(5,:,:) + a(6,4) * var(6,:,:) + &
+                   a(7,4) * var(7,:,:) ) * idel
+
+    il = 5
+  end if
+  if ( bb(2) == 0 ) then
+    ir = ni - gsize(1)
+  else
+    rhs(ni-3,:,:) = rhs(ni-3,:,:) + &
+                 ( a(1,4) * var(ni,:,:) + a(2,4) * var(ni-1,:,:) + &
+                   a(3,4) * var(ni-2,:,:) + a(4,4) * var(ni-3,:,:) + &
+                   a(5,4) * var(ni-4,:,:) + a(6,4) * var(ni-5,:,:) + &
+                   a(7,4) * var(ni-6,:,:) ) * idel
+    rhs(ni-2,:,:) = rhs(ni-2,:,:) + &
+                 ( a(1,3) * var(ni,:,:) + a(2,3) * var(ni-1,:,:) + &
+                   a(3,3) * var(ni-2,:,:) + a(4,3) * var(ni-3,:,:) + &
+                   a(5,3) * var(ni-4,:,:) + a(6,3) * var(ni-5,:,:) ) * idel
+    rhs(ni-1,:,:) = rhs(ni-1,:,:) + &
+                 ( a(1,2) * var(ni,:,:) + a(2,2) * var(ni-1,:,:) + &
+                   a(3,2) * var(ni-2,:,:) + a(4,2) * var(ni-3,:,:) + &
+                   a(5,2) * var(ni-4,:,:) ) * idel
+    rhs(ni,:,:)   = rhs(ni,:,:) + &
+                 ( a(1,1) * var(ni,:,:) + a(2,1) * var(ni-1,:,:) + &
+                   a(3,1) * var(ni-2,:,:) + a(4,1) * var(ni-3,:,:) ) * idel
+
+    ir = ni - 4
+  end if
+  rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+                   ( -20.0_wp * var(il:ir,:,:) + &
+                      15.0_wp * ( var(il-1:ir-1,:,:) + &
+                                  var(il+1:ir+1,:,:) ) - &
+                       6.0_wp * ( var(il-2:ir-2,:,:) + &
+                                  var(il+2:ir+2,:,:) ) + &
+                                ( var(il-3:ir-3,:,:) + &
+                                  var(il+3:ir+3,:,:) ) ) * idel
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 64 * delta(2) )
+
+  if ( bb(3) == 0 ) then
+    jl = 1 + gsize(2)
+  else
+    rhs(:,1,:) = rhs(:,1,:) + &
+                 ( a(1,1) * var(:,1,:) + a(2,1) * var(:,2,:) + &
+                   a(3,1) * var(:,3,:) + a(4,1) * var(:,4,:) ) * idel
+    rhs(:,2,:) = rhs(:,2,:) + &
+                 ( a(1,2) * var(:,1,:) + a(2,2) * var(:,2,:) + &
+                   a(3,2) * var(:,3,:) + a(4,2) * var(:,4,:) + &
+                   a(5,2) * var(:,5,:) ) * idel
+    rhs(:,3,:) = rhs(:,3,:) + &
+                 ( a(1,3) * var(:,1,:) + a(2,3) * var(:,2,:) + &
+                   a(3,3) * var(:,3,:) + a(4,3) * var(:,4,:) + &
+                   a(5,3) * var(:,5,:) + a(6,3) * var(:,6,:) ) * idel
+    rhs(:,4,:) = rhs(:,4,:) + &
+                 ( a(1,4) * var(:,1,:) + a(2,4) * var(:,2,:) + &
+                   a(3,4) * var(:,3,:) + a(4,4) * var(:,4,:) + &
+                   a(5,4) * var(:,5,:) + a(6,4) * var(:,6,:) + &
+                   a(7,4) * var(:,7,:) ) * idel
+
+    jl = 5
+  end if
+  if ( bb(4) == 0 ) then
+    jr = nj - gsize(2)
+  else
+    rhs(:,nj-3,:) = rhs(:,nj-3,:) + &
+                 ( a(1,4) * var(:,nj,:) + a(2,4) * var(:,nj-1,:) + &
+                   a(3,4) * var(:,nj-2,:) + a(4,4) * var(:,nj-3,:) + &
+                   a(5,4) * var(:,nj-4,:) + a(6,4) * var(:,nj-5,:) + &
+                   a(7,4) * var(:,nj-6,:) ) * idel
+    rhs(:,nj-2,:) = rhs(:,nj-2,:) + &
+                 ( a(1,3) * var(:,nj,:) + a(2,3) * var(:,nj-1,:) + &
+                   a(3,3) * var(:,nj-2,:) + a(4,3) * var(:,nj-3,:) + &
+                   a(5,3) * var(:,nj-4,:) + a(6,3) * var(:,nj-5,:) ) * idel
+    rhs(:,nj-1,:) = rhs(:,nj-1,:) + &
+                 ( a(1,2) * var(:,nj,:) + a(2,2) * var(:,nj-1,:) + &
+                   a(3,2) * var(:,nj-2,:) + a(4,2) * var(:,nj-3,:) + &
+                   a(5,2) * var(:,nj-4,:) ) * idel
+    rhs(:,nj,:)   = rhs(:,nj,:) + &
+                 ( a(1,1) * var(:,nj,:) + a(2,1) * var(:,nj-1,:) + &
+                   a(3,1) * var(:,nj-2,:) + a(4,1) * var(:,nj-3,:) ) * idel
+
+    jr = nj - 4
+  end if
+  rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+                   ( -20.0_wp * var(:,jl:jr,:) + &
+                      15.0_wp * ( var(:,jl-1:jr-1,:) + &
+                                  var(:,jl+1:jr+1,:) ) - &
+                       6.0_wp * ( var(:,jl-2:jr-2,:) + &
+                                  var(:,jl+2:jr+2,:) ) + &
+                                ( var(:,jl-3:jr-3,:) + &
+                                  var(:,jl+3:jr+3,:) ) ) * idel
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 64 * delta(3) )
+
+  if ( bb(5) == 0 ) then
+    kl = 1 + gsize(3)
+  else
+    rhs(:,:,1) = rhs(:,:,1) + &
+                 ( a(1,1) * var(:,:,1) + a(2,1) * var(:,:,2) + &
+                   a(3,1) * var(:,:,3) + a(4,1) * var(:,:,4) ) * idel
+    rhs(:,:,2) = rhs(:,:,2) + &
+                 ( a(1,2) * var(:,:,1) + a(2,2) * var(:,:,2) + &
+                   a(3,2) * var(:,:,3) + a(4,2) * var(:,:,4) + &
+                   a(5,2) * var(:,:,5) ) * idel
+    rhs(:,:,3) = rhs(:,:,3) + &
+                 ( a(1,3) * var(:,:,1) + a(2,3) * var(:,:,2) + &
+                   a(3,3) * var(:,:,3) + a(4,3) * var(:,:,4) + &
+                   a(5,3) * var(:,:,5) + a(6,3) * var(:,:,6) ) * idel
+    rhs(:,:,4) = rhs(:,:,4) + &
+                 ( a(1,4) * var(:,:,1) + a(2,4) * var(:,:,2) + &
+                   a(3,4) * var(:,:,3) + a(4,4) * var(:,:,4) + &
+                   a(5,4) * var(:,:,5) + a(6,4) * var(:,:,6) + &
+                   a(7,4) * var(:,:,7) ) * idel
+
+    kl = 5
+  end if
+  if ( bb(6) == 0 ) then
+    kr = nk - gsize(3)
+  else
+    rhs(:,:,nk-3) = rhs(:,:,nk-3) + &
+                 ( a(1,4) * var(:,:,nk) + a(2,4) * var(:,:,nk-1) + &
+                   a(3,4) * var(:,:,nk-2) + a(4,4) * var(:,:,nk-3) + &
+                   a(5,4) * var(:,:,nk-4) + a(6,4) * var(:,:,nk-5) + &
+                   a(7,4) * var(:,:,nk-6) ) * idel
+    rhs(:,:,nk-2) = rhs(:,:,nk-2) + &
+                 ( a(1,3) * var(:,:,nk) + a(2,3) * var(:,:,nk-1) + &
+                   a(3,3) * var(:,:,nk-2) + a(4,3) * var(:,:,nk-3) + &
+                   a(5,3) * var(:,:,nk-4) + a(6,3) * var(:,:,nk-5) ) * idel
+    rhs(:,:,nk-1) = rhs(:,:,nk-1) + &
+                 ( a(1,2) * var(:,:,nk) + a(2,2) * var(:,:,nk-1) + &
+                   a(3,2) * var(:,:,nk-2) + a(4,2) * var(:,:,nk-3) + &
+                   a(5,2) * var(:,:,nk-4) ) * idel
+    rhs(:,:,nk)   = rhs(:,:,nk) + &
+                 ( a(1,1) * var(:,:,nk) + a(2,1) * var(:,:,nk-1) + &
+                   a(3,1) * var(:,:,nk-2) + a(4,1) * var(:,:,nk-3) ) * idel
+
+    kr = nk - 4
+  end if
+  rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+                   ( -20.0_wp * var(:,:,kl:kr) + &
+                      15.0_wp * ( var(:,:,kl-1:kr-1) + &
+                                  var(:,:,kl+1:kr+1) ) - &
+                       6.0_wp * ( var(:,:,kl-2:kr-2) + &
+                                  var(:,:,kl+2:kr+2) ) + &
+                                ( var(:,:,kl-3:kr-3) + &
+                                  var(:,:,kl+3:kr+3) ) ) * idel
+
+contains
+
+  subroutine set_coeff ( a )
+
+    implicit none
+
+    CCTK_REAL, dimension(7,4), intent(OUT) :: a
+    CCTK_REAL :: zero = 0.0
+    integer, parameter :: wp = kind(zero)
+
+    a(1,1) = -2.8235294117647058823529411764705882352941176470588_wp
+    a(2,1) = 8.4705882352941176470588235294117647058823529411765_wp
+    a(3,1) = -8.4705882352941176470588235294117647058823529411765_wp
+    a(4,1) = 2.8235294117647058823529411764705882352941176470588_wp
+    a(5,1) = 0.0_wp
+    a(6,1) = 0.0_wp
+    a(7,1) = 0.0_wp
+    a(1,2) = 2.4406779661016949152542372881355932203389830508475_wp
+    a(2,2) = -8.1355932203389830508474576271186440677966101694915_wp
+    a(3,2) = 9.7627118644067796610169491525423728813559322033898_wp
+    a(4,2) = -4.8813559322033898305084745762711864406779661016949_wp
+    a(5,2) = 0.81355932203389830508474576271186440677966101694915_wp
+    a(6,2) = 0.0_wp
+    a(7,2) = 0.0_wp
+    a(1,3) = -3.3488372093023255813953488372093023255813953488372_wp
+    a(2,3) = 13.395348837209302325581395348837209302325581395349_wp
+    a(3,3) = -21.209302325581395348837209302325581395348837209302_wp
+    a(4,3) = 16.744186046511627906976744186046511627906976744186_wp
+    a(5,3) = -6.6976744186046511627906976744186046511627906976744_wp
+    a(6,3) = 1.1162790697674418604651162790697674418604651162791_wp
+    a(7,3) = 0.0_wp
+    a(1,4) = 0.97959183673469387755102040816326530612244897959184_wp
+    a(2,4) = -5.8775510204081632653061224489795918367346938775510_wp
+    a(3,4) = 14.693877551020408163265306122448979591836734693878_wp
+    a(4,4) = -19.591836734693877551020408163265306122448979591837_wp
+    a(5,4) = 14.693877551020408163265306122448979591836734693878_wp
+    a(6,4) = -5.8775510204081632653061224489795918367346938775510_wp
+    a(7,4) = 0.97959183673469387755102040816326530612244897959184_wp
+
+  end subroutine set_coeff
+
+end subroutine dissipation_4_2_alt