Implementation of a Kreiss-Oliger type dissipation operator compatible with

the diagonal 2-1 derivative operators.


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

3 files changed, 172 insertions, 2 deletions

diff --git a/src/Dissipation_2_1_alt.F90 b/src/Dissipation_2_1_alt.F90
new file mode 100644
index 0000000..49ee0f0
--- /dev/null
+++ b/src/Dissipation_2_1_alt.F90
@@ -0,0 +1,154 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+subroutine dissipation_2_1_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(4,2) :: a
+  CCTK_REAL :: idel
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 16 * 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,:,:) ) * 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,:,:) ) * idel
+
+    il = 3
+  end if
+  if ( bb(2) == 0 ) then
+    ir = ni - gsize(1)
+  else
+    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,:,:) ) * idel
+    rhs(ni,:,:)   = rhs(ni,:,:) + &
+                 ( a(1,1) * var(ni,:,:) + a(2,1) * var(ni-1,:,:) + &
+                   a(3,1) * var(ni-2,:,:) ) * idel
+
+    ir = ni - 2
+  end if
+  rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+                   ( -6.0_wp * var(il:ir,:,:) + &
+                      4.0_wp * ( var(il-1:ir-1,:,:) + &
+                                 var(il+1:ir+1,:,:) ) - &
+                               ( var(il-2:ir-2,:,:) + &
+                                 var(il+2:ir+2,:,:) ) ) * idel
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 16 * 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,:) ) * 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,:) ) * idel
+
+    jl = 3
+  end if
+  if ( bb(4) == 0 ) then
+    jr = nj - gsize(2)
+  else
+    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,:) ) * idel
+    rhs(:,nj,:)   = rhs(:,nj,:) + &
+                 ( a(1,1) * var(:,nj,:) + a(2,1) * var(:,nj-1,:) + &
+                   a(3,1) * var(:,nj-2,:) ) * idel
+
+    jr = nj - 2
+  end if
+  rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+                   ( -6.0_wp * var(:,jl:jr,:) + &
+                      4.0_wp * ( var(:,jl-1:jr-1,:) + &
+                                 var(:,jl+1:jr+1,:) ) - &
+                               ( var(:,jl-2:jr-2,:) + &
+                                 var(:,jl+2:jr+2,:) ) ) * idel
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 16 * 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) ) * 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) ) * idel
+
+    kl = 3
+  end if
+  if ( bb(6) == 0 ) then
+    kr = nk - gsize(3)
+  else
+    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) ) * idel
+    rhs(:,:,nk)   = rhs(:,:,nk) + &
+                 ( a(1,1) * var(:,:,nk) + a(2,1) * var(:,:,nk-1) + &
+                   a(3,1) * var(:,:,nk-2) ) * idel
+
+    kr = nk - 2
+  end if
+  rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+                   ( -6.0_wp * var(:,:,kl:kr) + &
+                      4.0_wp * ( var(:,:,kl-1:kr-1) + &
+                                 var(:,:,kl+1:kr+1) ) - &
+                               ( var(:,:,kl-2:kr-2) + &
+                                 var(:,:,kl+2:kr+2) ) ) * idel
+
+contains
+
+  subroutine set_coeff ( a )
+
+    implicit none
+
+    CCTK_REAL, dimension(4,2), intent(OUT) :: a
+    CCTK_REAL :: zero = 0.0
+    integer, parameter :: wp = kind(zero)
+
+    a(1,1) = -2.0_wp
+    a(2,1) = 4.0_wp
+    a(3,1) = -2.0_wp
+    a(4,1) = 0.0_wp
+    a(1,2) = 2.0_wp
+    a(2,2) = -5.0_wp
+    a(3,2) = 4.0_wp
+    a(4,2) = -1.0_wp
+
+  end subroutine set_coeff
+
+end subroutine dissipation_2_1_alt
diff --git a/src/dissipation.c b/src/dissipation.c
index a1557cc..29dff9f 100644
--- a/src/dissipation.c
+++ b/src/dissipation.c
@@ -46,6 +46,15 @@ void CCTK_FCALL CCTK_FNAME(dissipation_2_1) (const CCTK_REAL *var,
                                              const CCTK_REAL *dx,
                                              const CCTK_REAL *epsdis,
                                              CCTK_REAL *rhs);
+void CCTK_FCALL CCTK_FNAME(dissipation_2_1_alt) (const CCTK_REAL *var,
+                                             const CCTK_INT *ni,
+                                             const CCTK_INT *nj,
+                                             const CCTK_INT *nk,
+                                             const CCTK_INT *bbox,
+                                             const CCTK_INT *gsize,
+                                             const CCTK_REAL *dx,
+                                             const CCTK_REAL *epsdis,
+                                             CCTK_REAL *rhs);
 void CCTK_FCALL CCTK_FNAME(dissipation_4_2) (const CCTK_REAL *var,
                                              const CCTK_INT *ni,
                                              const CCTK_INT *nj,
@@ -171,9 +180,15 @@ apply (int const varindex, char const * const optstring, void * const arg)
   if ( CCTK_Equals(norm_type,"Diagonal") ) {
     switch(order) {
     case 2: {
-      CCTK_FNAME(dissipation_2_1)
-        (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+      if ( CCTK_Equals(dissipation_type,"Mattson-Svard-Nordstrom") ) {
+        CCTK_FNAME(dissipation_2_1)
+          (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
                             bbox, gsize, dx, &epsdis, rhsptr);
+      } else {
+        CCTK_FNAME(dissipation_2_1_alt)
+          (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+                            bbox, gsize, dx, &epsdis, rhsptr);
+      }
       break;
     }
     case 4: {
diff --git a/src/make.code.defn b/src/make.code.defn
index 6ff753b..87f521e 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -19,6 +19,7 @@ SRCS =  call_derivs.c \
         CheckGridSizes.F90 \
         dissipation.c \
         Dissipation_2_1.F90 \
+        Dissipation_2_1_alt.F90 \
         Dissipation_4_2.F90 \
         Dissipation_6_3.F90 \
         Dissipation_8_4.F90 \