Two part update of the restricted full norm dissipation operators. The first

part is an optimisation. Previously the coefficients were constructed every
time the dissipation was applied. This turned out to have a significant
overhead, especially for large domains. Now the coefficients are constructed
only once for each patch on each processor and stored. The thorn now uses
the aliased functions MultiPatch_GetMap and MultiPatch_GetMaps (provided by
thorn MultiPatch) to figure out how many patches are present and which patch
is currently being worked on. These dissipation operators will not work unless
the MultiPatch thorn is used. Right now this should not be a problem (but can
be changed if it is), unless somebody uses dissipation operators outside of
the multipatch infrastructure.

The second part is an update aimed at making the restricted full norm
dissipation operators better behaved on grids with very different number of
grid points in different directions. For that reason the parameter:

diss_fraction

has been made vector valued (of length 3), so that the width of the transition
region from the boundary to the interior operators can be controlled
separately in the 3 dimensions. Unfortunately this means that any parameter
files, that explicitly sets this parameter, will have to be modified. I don't
think too many people are using the restricted full dissipation operators yet,
so it shouldn't cause too many problems. In additiona new parameter:

h_scaling

has been introduced. This parameter is also vector valued of length 3 and
allows the local grid spacing to be rescaled to more closely resemble the
physical grid spacing in any dimension.

Note, these changes only affect the non-diagonal dissipation operators.


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

9 files changed, 512 insertions, 281 deletions

diff --git a/interface.ccl b/interface.ccl
index 0c85802..e9d1e9e 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -51,6 +51,12 @@ CCTK_INT FUNCTION \
     SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH)
 REQUIRES FUNCTION SymmetryTableHandleForGrid
 
+CCTK_INT FUNCTION MultiPatch_GetMap (CCTK_POINTER_TO_CONST IN cctkGH)
+USES FUNCTION MultiPatch_GetMap
+
+CCTK_INT FUNCTION MultiPatch_GetMaps (CCTK_POINTER_TO_CONST IN cctkGH)
+USES FUNCTION MultiPatch_GetMaps
+
 CCTK_INT FUNCTION SymmetryHandleOfName (CCTK_STRING IN sym_name)
 REQUIRES FUNCTION SymmetryHandleOfName
 
diff --git a/param.ccl b/param.ccl
index 7518ac7..32cb36e 100644
--- a/param.ccl
+++ b/param.ccl
@@ -38,11 +38,16 @@ REAL epsdis "Dissipation strength" STEERABLE=always
   *:* :: "Values typical between 0 and 1"
 } 0.2
 
-REAL diss_fraction "Fractional size of the transition region for the full restricted dissipation operator"
+REAL diss_fraction[3] "Fractional size of the transition region for the full restricted dissipation operator"
 {
   0:0.5 :: ""
 } 0.2
 
+REAL h_scaling[3] "Scaling factor for the local grid spacing in the dissipation operators"
+{
+ 0:* :: "Positive please"
+} 1.0
+
 STRING vars "List of evolved grid functions that should have dissipation added" STEERABLE=always
 {
   .* :: "Must be a valid list of grid functions"
diff --git a/src/Dissipation_4_3_min_err_coeff.F90 b/src/Dissipation_4_3_min_err_coeff.F90
index 1a196ac..7d61389 100644
--- a/src/Dissipation_4_3_min_err_coeff.F90
+++ b/src/Dissipation_4_3_min_err_coeff.F90
@@ -2,14 +2,18 @@
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
+#include "cctk_Functions.h"
 
 
 subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, &
-                                               delta, epsilon, dfl, rhs)
+                                delta, epsilon, dfl, npatches, patch, rhs)
+
+  use dissipation_coeff
 
   implicit none
 
   DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
 
   CCTK_INT, dimension(3), intent(in) ::  lsh, gsh, lbnd
   CCTK_INT :: ni, nj, nk
@@ -18,28 +22,52 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, &
   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, dfl
+  CCTK_REAL, intent(in) :: epsilon
+  CCTK_REAL, dimension(3), intent(in) :: dfl
+  CCTK_INT, intent(in) :: npatches, patch
+
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
   integer :: i, j, k
-  CCTK_REAL, dimension(:,:), allocatable :: a, d, b, h
+  CCTK_REAL, dimension(:,:), allocatable :: atmp, d, b, h
+  CCTK_REAL, dimension(:,:), pointer :: a
   CCTK_REAL :: idel
 
   CCTK_INT :: il, ir, jl, jr, kl, kr
 
   ni = lsh(1); nj = lsh(2); nk = lsh(3)
 
-  allocate ( a(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni) )
-  a = zero; d = zero; b = zero; h = zero;
+  if ( first ) then
+    allocate ( savedx(npatches), savedy(npatches), savedz(npatches) )
+    savedx = .false.; savedy = .false.; savedz = .false.
+    allocate ( xcoeff(npatches), ycoeff(npatches), zcoeff(npatches) )
+    first = .false.
+  end if
+
+  if ( .not. savedx(patch) ) then
+    allocate ( atmp(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni) )
+    atmp = zero; d = zero; b = zero; h = zero
 
-  call set_dmatrix ( d, bb(1:2) )
+    call set_dmatrix ( d, bb(1:2) )
 
-  call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl )
+    call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl(1) )
   
-  call set_hmatrix ( h, bb(1:2) )
+    call set_hmatrix ( h, bb(1:2) )
+
+    atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+
+    allocate ( xcoeff(patch)%coeff(ni,ni) )
+
+    xcoeff(patch)%coeff = atmp
 
-  a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+    savedx(patch) = .true.
+
+    deallocate ( atmp, d, b, h )
+
+  end if
+
+  a => xcoeff(patch)%coeff
 
   if ( scale_with_h > 0 ) then
     idel = epsilon / ( 16 * delta(1) )
@@ -125,20 +153,31 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, &
                      a(i+2,i) * var(i+2,:,:) ) * idel
   end do
 
-  deallocate ( a, d, b, h )
-
   if ( zero_derivs_y == 0 ) then
-    allocate ( a(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj) )
-    a = zero; d = zero; b = zero; h = zero;
+    if ( .not. savedy(patch) ) then
+      allocate ( atmp(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj) )
+      atmp = zero; d = zero; b = zero; h = zero
   
-    call set_dmatrix ( d, bb(3:4) )
+      call set_dmatrix ( d, bb(3:4) )
   
-    call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl )
+      call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl(2) )
     
-    call set_hmatrix ( h, bb(3:4) )
+      call set_hmatrix ( h, bb(3:4) )
   
-    a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+      atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
   
+      allocate ( ycoeff(patch)%coeff(nj,nj) )
+
+      ycoeff(patch)%coeff = atmp
+
+      savedy(patch) = .true.
+
+      deallocate ( atmp, d, b, h )
+
+    end if
+
+    a => ycoeff(patch)%coeff
+
     if ( scale_with_h > 0 ) then
       idel = epsilon / ( 16 * delta(2) )
     else
@@ -224,21 +263,33 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, &
                        a(j+2,j) * var(:,j+2,:) ) * idel
     end do
   
-    deallocate ( a, d, b, h )
   end if
 
   if ( zero_derivs_z == 0 ) then
-    allocate ( a(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk) )
-    a = zero; d = zero; b = zero; h = zero;
+    if ( .not. savedz(patch) ) then
+      allocate ( atmp(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk) )
+      atmp = zero; d = zero; b = zero; h = zero
   
-    call set_dmatrix ( d, bb(5:6) )
+      call set_dmatrix ( d, bb(5:6) )
   
-    call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl )
+      call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl(3) )
   
-    call set_hmatrix ( h, bb(5:6) )
+      call set_hmatrix ( h, bb(5:6) )
   
-    a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+      atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
   
+      allocate ( zcoeff(patch)%coeff(nk,nk) )
+
+      zcoeff(patch)%coeff = atmp
+
+      savedz(patch) = .true.
+
+      deallocate ( atmp, d, b, h )
+
+    end if
+
+    a => zcoeff(patch)%coeff
+
     if ( scale_with_h > 0 ) then
       idel = epsilon / ( 16 * delta(3) )
     else
@@ -324,7 +375,6 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, &
                        a(k+2,k) * var(:,:,k+2) ) * idel
     end do
   
-    deallocate ( a, d, b, h )
   end if
 
 contains
diff --git a/src/Dissipation_4_3_min_err_coeff_alt.F90 b/src/Dissipation_4_3_min_err_coeff_alt.F90
index ef619d0..3310704 100644
--- a/src/Dissipation_4_3_min_err_coeff_alt.F90
+++ b/src/Dissipation_4_3_min_err_coeff_alt.F90
@@ -6,7 +6,9 @@
 
 
 subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, &
-                                               delta, epsilon, dfl, rhs)
+                                delta, epsilon, dfl, npatches, patch, rhs)
+
+  use dissipation_coeff
 
   implicit none
 
@@ -20,65 +22,88 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, &
   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, dfl
+  CCTK_REAL, intent(in) :: epsilon
+  CCTK_REAL, dimension(3), intent(in) :: dfl
+  CCTK_INT, intent(in) :: npatches, patch
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
   integer :: i, j, k, center
-  CCTK_REAL, dimension(:,:), allocatable :: a, d, b, h, tmp
+  CCTK_REAL, dimension(:,:), allocatable :: atmp, d, b, h, tmp
+  CCTK_REAL, dimension(:,:), pointer :: a
   CCTK_REAL :: idel
 
   CCTK_INT :: il, ir, jl, jr, kl, kr
 
   ni = lsh(1); nj = lsh(2); nk = lsh(3)
 
-  allocate ( a(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni), tmp(ni,ni) )
-  a = zero; d = zero; b = zero; h = zero; tmp = zero
+  if ( first ) then
+    allocate ( savedx(npatches), savedy(npatches), savedz(npatches) )
+    savedx = .false.; savedy = .false.; savedz = .false.
+    allocate ( xcoeff(npatches), ycoeff(npatches), zcoeff(npatches) )
+    first = .false.
+  end if
 
-  call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl )
+  if ( .not. savedx(patch) ) then
+    allocate ( atmp(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni), tmp(ni,ni) )
+    atmp = zero; d = zero; b = zero; h = zero; tmp = zero
 
-  call set_hmatrix ( h, bb(1:2) )
+    call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl(1) )
 
-  if ( any ( bb(1:2) == 0 ) ) then
+    call set_hmatrix ( h, bb(1:2) )
 
-    call set_dmatrix ( d, bb(1:2) )
+    if ( any ( bb(1:2) == 0 ) ) then
 
-    a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+      call set_dmatrix ( d, bb(1:2) )
 
-  else
+      atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
 
-    if ( lsh(1) /= gsh(1) ) then
-      call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
-    end if
-
-    center = gsh(1) / 2
-    ir = center + max(3,gsize(1))
+    else
 
-    call set_dmatrix_half ( d(1:ir,1:ir), 1 )
+      if ( lsh(1) /= gsh(1) ) then
+        call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
+      end if
 
-    tmp(1:ir,1:ir) = -transpose ( &
-                        matmul ( h(1:ir,1:ir), &
-                                 matmul ( &
-                                   transpose ( d(1:ir,1:ir) ), &
-                                    matmul ( b(1:ir,1:ir), &
-                                             d(1:ir,1:ir) ) ) ) )
+      center = gsh(1) / 2
+      ir = center + max(3,gsize(1))
 
-    a(1:ir,1:center) = tmp(1:ir,1:center)
+      call set_dmatrix_half ( d(1:ir,1:ir), 1 )
 
-    il = center + 1 - max(3,gsize(1))
-    d = zero
-    call set_dmatrix_half ( d(il:ni,il:ni), 2 )
-    tmp(il:ni,il:ni) = -transpose ( &
-                          matmul ( h(il:ni,il:ni), &
+      tmp(1:ir,1:ir) = -transpose ( &
+                          matmul ( h(1:ir,1:ir), &
                                    matmul ( &
-                                     transpose ( d(il:ni,il:ni) ), &
-                                     matmul ( b(il:ni,il:ni), &
-                                              d(il:ni,il:ni) ) ) ) )
+                                     transpose ( d(1:ir,1:ir) ), &
+                                      matmul ( b(1:ir,1:ir), &
+                                               d(1:ir,1:ir) ) ) ) )
+
+      atmp(1:ir,1:center) = tmp(1:ir,1:center)
+
+      il = center + 1 - max(3,gsize(1))
+      d = zero
+      call set_dmatrix_half ( d(il:ni,il:ni), 2 )
+      tmp(il:ni,il:ni) = -transpose ( &
+                            matmul ( h(il:ni,il:ni), &
+                                     matmul ( &
+                                       transpose ( d(il:ni,il:ni) ), &
+                                       matmul ( b(il:ni,il:ni), &
+                                                d(il:ni,il:ni) ) ) ) )
+
+      atmp(il:ni,center+1:ni) = tmp(il:ni,center+1:ni)
 
-    a(il:ni,center+1:ni) = tmp(il:ni,center+1:ni)
+    end if
+
+    allocate ( xcoeff(patch)%coeff(ni,ni) )
+
+    xcoeff(patch)%coeff = atmp
+
+    savedx(patch) = .true.
+
+    deallocate ( atmp, d, b, h, tmp )
 
   end if
 
+  a => xcoeff(patch)%coeff
+
   idel = epsilon / ( 64 * delta(1) )
 
   if ( bb(1) == 0 ) then
@@ -180,59 +205,69 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, &
 
   end do
 
-  deallocate ( a, d, b, h, tmp )
-
   if ( zero_derivs_y == 0 ) then
-    allocate ( a(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj), tmp(nj,nj) )
-    a = zero; d = zero; b = zero; h = zero;
+    if ( .not. savedy(patch) ) then
+      allocate ( atmp(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj), tmp(nj,nj) )
+      atmp = zero; d = zero; b = zero; h = zero; tmp = zero
   
-    call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl )
+      call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl(2) )
     
-    call set_hmatrix ( h, bb(3:4) )
+      call set_hmatrix ( h, bb(3:4) )
   
-    if ( any ( bb(3:4) == 0 ) ) then
+      if ( any ( bb(3:4) == 0 ) ) then
   
-      call set_dmatrix ( d, bb(3:4) )
+        call set_dmatrix ( d, bb(3:4) )
   
-      a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+        atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
   
-    else
+      else
   
-      if ( lsh(2) /= gsh(2) ) then
-        call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
-      end if
+        if ( lsh(2) /= gsh(2) ) then
+          call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
+        end if
   
-      center = gsh(2) / 2
-      jr = center + max(3,gsize(2))
+        center = gsh(2) / 2
+        jr = center + max(3,gsize(2))
   
-      call set_dmatrix_half ( d(1:jr,1:jr), 1 )
+        call set_dmatrix_half ( d(1:jr,1:jr), 1 )
   
-      tmp(1:jr,1:jr) = -transpose ( &
-                          matmul ( h(1:jr,1:jr), &
-                                   matmul ( &
-                                     transpose ( d(1:jr,1:jr) ), &
-                                      matmul ( b(1:jr,1:jr), &
-                                               d(1:jr,1:jr) ) ) ) )
+        tmp(1:jr,1:jr) = -transpose ( &
+                            matmul ( h(1:jr,1:jr), &
+                                     matmul ( &
+                                       transpose ( d(1:jr,1:jr) ), &
+                                        matmul ( b(1:jr,1:jr), &
+                                                 d(1:jr,1:jr) ) ) ) )
   
-      a(1:jr,1:center) = tmp(1:jr,1:center)
+        atmp(1:jr,1:center) = tmp(1:jr,1:center)
   
-      jl = center + 1 - max(3,gsize(2))
-      d = zero
-      call set_dmatrix_half ( d(jl:nj,jl:nj), 2 )
-      tmp(jl:nj,jl:nj) = -transpose ( &
-                            matmul ( h(jl:nj,jl:nj), &
-                                     matmul ( &
-                                       transpose ( d(jl:nj,jl:nj) ), &
-                                       matmul ( b(jl:nj,jl:nj), &
-                                                d(jl:nj,jl:nj) ) ) ) )
+        jl = center + 1 - max(3,gsize(2))
+        d = zero
+        call set_dmatrix_half ( d(jl:nj,jl:nj), 2 )
+        tmp(jl:nj,jl:nj) = -transpose ( &
+                              matmul ( h(jl:nj,jl:nj), &
+                                       matmul ( &
+                                         transpose ( d(jl:nj,jl:nj) ), &
+                                         matmul ( b(jl:nj,jl:nj), &
+                                                  d(jl:nj,jl:nj) ) ) ) )
   
-      a(jl:nj,center+1:nj) = tmp(jl:nj,center+1:nj)
+        atmp(jl:nj,center+1:nj) = tmp(jl:nj,center+1:nj)
   
+        allocate ( ycoeff(patch)%coeff(nj,nj) )
+
+        ycoeff(patch)%coeff = atmp
+
+        savedy(patch) = .true.
+
+        deallocate ( atmp, d, b, h, tmp )
+
+      end if
+
     end if
   
+    a => ycoeff(patch)%coeff
+
     idel = epsilon / ( 64 * delta(2) )
    
-  
     if ( bb(3) == 0 ) then
   
       jl = 1 + gsize(2)
@@ -332,57 +367,69 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, &
   
     end do
   
-    deallocate ( a, d, b, h, tmp )
   end if
 
   if ( zero_derivs_z == 0 ) then
-    allocate ( a(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk), tmp(nk,nk) )
-    a = zero; d = zero; b = zero; h = zero; tmp = zero
+    if ( .not. savedz(patch) ) then
+      allocate ( atmp(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk), tmp(nk,nk) )
+      atmp = zero; d = zero; b = zero; h = zero; tmp = zero
   
-    call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl )
+      call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl(3) )
   
-    call set_hmatrix ( h, bb(5:6) )
+      call set_hmatrix ( h, bb(5:6) )
   
-    if ( any ( bb(5:6) == 0 ) ) then
+      if ( any ( bb(5:6) == 0 ) ) then
   
-      call set_dmatrix ( d, bb(5:6) )
+        call set_dmatrix ( d, bb(5:6) )
   
-      a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+        atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
   
-    else
+      else
   
-      if ( lsh(3) /= gsh(3) ) then
-        call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
-      end if
+        if ( lsh(3) /= gsh(3) ) then
+          call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
+        end if
   
-      center = gsh(3) / 2
-      kr = center + max(3,gsize(3))
+        center = gsh(3) / 2
+        kr = center + max(3,gsize(3))
   
-      call set_dmatrix_half ( d(1:kr,1:kr), 1 )
+        call set_dmatrix_half ( d(1:kr,1:kr), 1 )
   
-      tmp(1:kr,1:kr) = -transpose ( &
-                          matmul ( h(1:kr,1:kr), &
-                                   matmul ( &
-                                     transpose ( d(1:kr,1:kr) ), &
-                                      matmul ( b(1:kr,1:kr), &
-                                               d(1:kr,1:kr) ) ) ) )
+        tmp(1:kr,1:kr) = -transpose ( &
+                            matmul ( h(1:kr,1:kr), &
+                                     matmul ( &
+                                       transpose ( d(1:kr,1:kr) ), &
+                                        matmul ( b(1:kr,1:kr), &
+                                                 d(1:kr,1:kr) ) ) ) )
   
-      a(1:kr,1:center) = tmp(1:kr,1:center)
+        atmp(1:kr,1:center) = tmp(1:kr,1:center)
   
-      kl = center + 1 - max(3,gsize(3))
-      d = zero
-      call set_dmatrix_half ( d(kl:nk,kl:nk), 2 )
-      tmp(kl:nk,kl:nk) = -transpose ( &
-                            matmul ( h(kl:nk,kl:nk), &
-                                     matmul ( &
-                                       transpose ( d(kl:nk,kl:nk) ), &
-                                       matmul ( b(kl:nk,kl:nk), &
-                                                d(kl:nk,kl:nk) ) ) ) )
+        kl = center + 1 - max(3,gsize(3))
+        d = zero
+        call set_dmatrix_half ( d(kl:nk,kl:nk), 2 )
+        tmp(kl:nk,kl:nk) = -transpose ( &
+                              matmul ( h(kl:nk,kl:nk), &
+                                       matmul ( &
+                                         transpose ( d(kl:nk,kl:nk) ), &
+                                         matmul ( b(kl:nk,kl:nk), &
+                                                  d(kl:nk,kl:nk) ) ) ) )
   
-      a(kl:nk,center+1:nk) = tmp(kl:nk,center+1:nk)
+        atmp(kl:nk,center+1:nk) = tmp(kl:nk,center+1:nk)
   
+        allocate ( zcoeff(patch)%coeff(nk,nk) )
+
+        zcoeff(patch)%coeff = atmp
+
+        savedz(patch) = .true.
+
+        deallocate ( atmp, d, b, h, tmp )
+
+      end if
+
     end if
   
+    a => zcoeff(patch)%coeff
+
     idel = epsilon / ( 64 * delta(3) )
   
     if ( bb(5) == 0 ) then
@@ -484,7 +531,6 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, &
   
     end do
   
-    deallocate ( a, d, b, h, tmp )
   end if
 
 contains
diff --git a/src/Dissipation_6_5_min_err_coeff.F90 b/src/Dissipation_6_5_min_err_coeff.F90
index 63e0fbe..991d1ff 100644
--- a/src/Dissipation_6_5_min_err_coeff.F90
+++ b/src/Dissipation_6_5_min_err_coeff.F90
@@ -6,7 +6,9 @@
 
 
 subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
-                                               delta, epsilon, dfl, rhs)
+                                delta, epsilon, dfl, npatches, patch, rhs)
+
+  use dissipation_coeff
 
   implicit none
 
@@ -20,65 +22,88 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
   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, dfl
+  CCTK_REAL, intent(in) :: epsilon
+  CCTK_REAL, dimension(3), intent(in) :: dfl
+  CCTK_INT, intent(in) :: npatches, patch
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
   integer :: i, j, k, center
-  CCTK_REAL, dimension(:,:), allocatable :: a, d, b, h, tmp
+  CCTK_REAL, dimension(:,:), allocatable :: atmp, d, b, h, tmp
+  CCTK_REAL, dimension(:,:), pointer :: a
   CCTK_REAL :: idel
 
   CCTK_INT :: il, ir, jl, jr, kl, kr
 
   ni = lsh(1); nj = lsh(2); nk = lsh(3)
 
-  allocate ( a(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni), tmp(ni,ni) )
-  a = zero; d = zero; b = zero; h = zero; tmp = zero
+  if ( first ) then
+    allocate ( savedx(npatches), savedy(npatches), savedz(npatches) )
+    savedx = .false.; savedy = .false.; savedz = .false.
+    allocate ( xcoeff(npatches), ycoeff(npatches), zcoeff(npatches) )
+    first = .false.
+  end if
 
-  call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl )
+  if ( .not. savedx(patch) ) then
+    allocate ( atmp(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni), tmp(ni,ni) )
+    atmp = zero; d = zero; b = zero; h = zero; tmp = zero
 
-  call set_hmatrix ( h, bb(1:2) )
+    call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl(1) )
 
-  if ( any ( bb(1:2) == 0 ) ) then
+    call set_hmatrix ( h, bb(1:2) )
 
-    call set_dmatrix ( d, bb(1:2) )
+    if ( any ( bb(1:2) == 0 ) ) then
 
-    a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+      call set_dmatrix ( d, bb(1:2) )
 
-  else
+      atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
 
-    if ( lsh(1) /= gsh(1) ) then
-      call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
-    end if
-
-    center = gsh(1) / 2
-    ir = center + gsize(1)
+    else
 
-    call set_dmatrix_half ( d(1:ir,1:ir), 1 )
+      if ( lsh(1) /= gsh(1) ) then
+        call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
+      end if
 
-    tmp(1:ir,1:ir) = -transpose ( &
-                        matmul ( h(1:ir,1:ir), &
-                                 matmul ( &
-                                   transpose ( d(1:ir,1:ir) ), &
-                                    matmul ( b(1:ir,1:ir), &
-                                             d(1:ir,1:ir) ) ) ) )
+      center = gsh(1) / 2
+      ir = center + gsize(1)
 
-    a(1:ir,1:center) = tmp(1:ir,1:center)
+      call set_dmatrix_half ( d(1:ir,1:ir), 1 )
 
-    il = center + 1 - gsize(1)
-    d = zero
-    call set_dmatrix_half ( d(il:ni,il:ni), 2 )
-    tmp(il:ni,il:ni) = -transpose ( &
-                          matmul ( h(il:ni,il:ni), &
+      tmp(1:ir,1:ir) = -transpose ( &
+                          matmul ( h(1:ir,1:ir), &
                                    matmul ( &
-                                     transpose ( d(il:ni,il:ni) ), &
-                                     matmul ( b(il:ni,il:ni), &
-                                              d(il:ni,il:ni) ) ) ) )
+                                     transpose ( d(1:ir,1:ir) ), &
+                                      matmul ( b(1:ir,1:ir), &
+                                               d(1:ir,1:ir) ) ) ) )
+
+      atmp(1:ir,1:center) = tmp(1:ir,1:center)
+
+      il = center + 1 - gsize(1)
+      d = zero
+      call set_dmatrix_half ( d(il:ni,il:ni), 2 )
+      tmp(il:ni,il:ni) = -transpose ( &
+                            matmul ( h(il:ni,il:ni), &
+                                     matmul ( &
+                                       transpose ( d(il:ni,il:ni) ), &
+                                       matmul ( b(il:ni,il:ni), &
+                                                d(il:ni,il:ni) ) ) ) )
+
+      atmp(il:ni,center+1:ni) = tmp(il:ni,center+1:ni)
+
+    end if
+
+    allocate ( xcoeff(patch)%coeff(ni,ni) )
+
+    xcoeff(patch)%coeff = atmp
+
+    savedx(patch) = .true.
 
-    a(il:ni,center+1:ni) = tmp(il:ni,center+1:ni)
+    deallocate ( atmp, d, b, h, tmp )
 
   end if
 
+  a => xcoeff(patch)%coeff
+
   if ( scale_with_h > 0 ) then
     idel = epsilon / ( 64 * delta(1) )
   else
@@ -230,56 +255,67 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
 
   end do
 
-  deallocate ( a, d, b, h, tmp )
-
   if ( zero_derivs_y == 0 ) then
-    allocate ( a(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj), tmp(nj,nj) )
-    a = zero; d = zero; b = zero; h = zero;
+    if ( .not. savedy(patch) ) then
+      allocate ( atmp(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj), tmp(nj,nj) )
+      atmp = zero; d = zero; b = zero; h = zero; tmp = zero
   
-    call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl )
+      call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl(2) )
     
-    call set_hmatrix ( h, bb(3:4) )
-  
-    if ( any ( bb(3:4) == 0 ) ) then
+      call set_hmatrix ( h, bb(3:4) )
   
-      call set_dmatrix ( d, bb(3:4) )
+      if ( any ( bb(3:4) == 0 ) ) then
   
-      a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+        call set_dmatrix ( d, bb(3:4) )
   
-    else
-  
-      if ( lsh(2) /= gsh(2) ) then
-        call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
-      end if
+        atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
   
-      center = gsh(2) / 2
-      jr = center + gsize(2)
+      else
   
-      call set_dmatrix_half ( d(1:jr,1:jr), 1 )
+        if ( lsh(2) /= gsh(2) ) then
+          call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
+        end if
   
-      tmp(1:jr,1:jr) = -transpose ( &
-                          matmul ( h(1:jr,1:jr), &
-                                   matmul ( &
-                                     transpose ( d(1:jr,1:jr) ), &
-                                      matmul ( b(1:jr,1:jr), &
-                                               d(1:jr,1:jr) ) ) ) )
+        center = gsh(2) / 2
+        jr = center + gsize(2)
   
-      a(1:jr,1:center) = tmp(1:jr,1:center)
+        call set_dmatrix_half ( d(1:jr,1:jr), 1 )
   
-      jl = center + 1 - gsize(2)
-      d = zero
-      call set_dmatrix_half ( d(jl:nj,jl:nj), 2 )
-      tmp(jl:nj,jl:nj) = -transpose ( &
-                            matmul ( h(jl:nj,jl:nj), &
+        tmp(1:jr,1:jr) = -transpose ( &
+                            matmul ( h(1:jr,1:jr), &
                                      matmul ( &
-                                       transpose ( d(jl:nj,jl:nj) ), &
-                                       matmul ( b(jl:nj,jl:nj), &
-                                                d(jl:nj,jl:nj) ) ) ) )
-  
-      a(jl:nj,center+1:nj) = tmp(jl:nj,center+1:nj)
+                                       transpose ( d(1:jr,1:jr) ), &
+                                        matmul ( b(1:jr,1:jr), &
+                                                 d(1:jr,1:jr) ) ) ) )
+  
+        atmp(1:jr,1:center) = tmp(1:jr,1:center)
+  
+        jl = center + 1 - gsize(2)
+        d = zero
+        call set_dmatrix_half ( d(jl:nj,jl:nj), 2 )
+        tmp(jl:nj,jl:nj) = -transpose ( &
+                              matmul ( h(jl:nj,jl:nj), &
+                                       matmul ( &
+                                         transpose ( d(jl:nj,jl:nj) ), &
+                                         matmul ( b(jl:nj,jl:nj), &
+                                                  d(jl:nj,jl:nj) ) ) ) )
+    
+        atmp(jl:nj,center+1:nj) = tmp(jl:nj,center+1:nj)
   
-    end if
+        allocate ( ycoeff(patch)%coeff(nj,nj) )
+
+        ycoeff(patch)%coeff = atmp
+
+        savedy(patch) = .true.
+
+        deallocate ( atmp, d, b, h, tmp )
+
+      end if
+
+    end if 
   
+    a => ycoeff(patch)%coeff
+
     if ( scale_with_h > 0 ) then
       idel = epsilon / ( 64 * delta(2) )
     else
@@ -431,57 +467,69 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
   
     end do
   
-    deallocate ( a, d, b, h, tmp )
   end if
 
   if ( zero_derivs_z == 0 ) then
-    allocate ( a(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk), tmp(nk,nk) )
-    a = zero; d = zero; b = zero; h = zero; tmp = zero
+    if ( .not. savedz(patch) ) then
+      allocate ( atmp(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk), tmp(nk,nk) )
+      atmp = zero; d = zero; b = zero; h = zero; tmp = zero
   
-    call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl )
+      call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl(3) )
   
-    call set_hmatrix ( h, bb(5:6) )
+      call set_hmatrix ( h, bb(5:6) )
   
-    if ( any ( bb(5:6) == 0 ) ) then
+      if ( any ( bb(5:6) == 0 ) ) then
   
-      call set_dmatrix ( d, bb(5:6) )
+        call set_dmatrix ( d, bb(5:6) )
   
-      a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+        atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
   
-    else
+      else
   
-      if ( lsh(3) /= gsh(3) ) then
-        call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
-      end if
+        if ( lsh(3) /= gsh(3) ) then
+          call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' )
+        end if
   
-      center = gsh(3) / 2
-      kr = center + gsize(3)
+        center = gsh(3) / 2
+        kr = center + gsize(3)
   
-      call set_dmatrix_half ( d(1:kr,1:kr), 1 )
+        call set_dmatrix_half ( d(1:kr,1:kr), 1 )
   
-      tmp(1:kr,1:kr) = -transpose ( &
-                          matmul ( h(1:kr,1:kr), &
-                                   matmul ( &
-                                     transpose ( d(1:kr,1:kr) ), &
-                                      matmul ( b(1:kr,1:kr), &
-                                               d(1:kr,1:kr) ) ) ) )
+        tmp(1:kr,1:kr) = -transpose ( &
+                            matmul ( h(1:kr,1:kr), &
+                                     matmul ( &
+                                       transpose ( d(1:kr,1:kr) ), &
+                                        matmul ( b(1:kr,1:kr), &
+                                                 d(1:kr,1:kr) ) ) ) )
   
-      a(1:kr,1:center) = tmp(1:kr,1:center)
+        atmp(1:kr,1:center) = tmp(1:kr,1:center)
   
-      kl = center + 1 - gsize(3)
-      d = zero
-      call set_dmatrix_half ( d(kl:nk,kl:nk), 2 )
-      tmp(kl:nk,kl:nk) = -transpose ( &
-                            matmul ( h(kl:nk,kl:nk), &
-                                     matmul ( &
-                                       transpose ( d(kl:nk,kl:nk) ), &
-                                       matmul ( b(kl:nk,kl:nk), &
-                                                d(kl:nk,kl:nk) ) ) ) )
+        kl = center + 1 - gsize(3)
+        d = zero
+        call set_dmatrix_half ( d(kl:nk,kl:nk), 2 )
+        tmp(kl:nk,kl:nk) = -transpose ( &
+                              matmul ( h(kl:nk,kl:nk), &
+                                       matmul ( &
+                                         transpose ( d(kl:nk,kl:nk) ), &
+                                         matmul ( b(kl:nk,kl:nk), &
+                                                  d(kl:nk,kl:nk) ) ) ) )
   
-      a(kl:nk,center+1:nk) = tmp(kl:nk,center+1:nk)
+        atmp(kl:nk,center+1:nk) = tmp(kl:nk,center+1:nk)
   
-    end if
+        allocate ( zcoeff(patch)%coeff(nk,nk) )
+
+        zcoeff(patch)%coeff = atmp
+
+        savedz(patch) = .true.
+
+        deallocate ( atmp, d, b, h, tmp )
+
+      end if
   
+    end if
+
+    a => zcoeff(patch)%coeff
+
     if ( scale_with_h > 0 ) then
       idel = epsilon / ( 64 * delta(3) )
     else
@@ -633,7 +681,6 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, &
   
     end do
   
-    deallocate ( a, d, b, h, tmp )
   end if
 
 contains
diff --git a/src/Dissipation_6_5_min_err_coeff_alt.F90 b/src/Dissipation_6_5_min_err_coeff_alt.F90
index ac9b19a..e67b2a4 100644
--- a/src/Dissipation_6_5_min_err_coeff_alt.F90
+++ b/src/Dissipation_6_5_min_err_coeff_alt.F90
@@ -6,7 +6,9 @@
 
 
 subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, &
-                                               delta, epsdisl, dfl, rhs)
+                                delta, epsilon, dfl, npatches, patch, rhs)
+
+  use dissipation_coeff
 
   implicit none
 
@@ -20,30 +22,53 @@ subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, &
   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) :: epsdisl, dfl
+  CCTK_REAL, intent(in) :: epsilon
+  CCTK_REAL, dimension(3), intent(in) :: dfl
+  CCTK_INT, intent(in) :: npatches, patch
 
   CCTK_REAL :: zero = 0.0
   integer, parameter :: wp = kind(zero)
   integer :: i, j, k, center
-  CCTK_REAL, dimension(:,:), allocatable :: a, d, b, h, tmp
+  CCTK_REAL, dimension(:,:), allocatable :: atmp, d, b, h
+  CCTK_REAL, dimension(:,:), pointer :: a
   CCTK_REAL :: idel
 
   CCTK_INT :: il, ir, jl, jr, kl, kr
 
   ni = lsh(1); nj = lsh(2); nk = lsh(3)
 
-  allocate ( a(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni), tmp(ni,ni) )
-  a = zero; d = zero; b = zero; h = zero; tmp = zero
+  if ( first ) then
+    allocate ( savedx(npatches), savedy(npatches), savedz(npatches) )
+    savedx = .false.; savedy = .false.; savedz = .false.
+    allocate ( xcoeff(npatches), ycoeff(npatches), zcoeff(npatches) )
+    first = .false.
+  end if
+
+  if ( .not. savedx(patch) ) then
+    allocate ( atmp(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni) )
+    atmp = zero; d = zero; b = zero; h = zero
+
+    call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl(1) )
+
+    call set_hmatrix ( h, bb(1:2) )
 
-  call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl )
+    call set_dmatrix ( d, bb(1:2) )
 
-  call set_hmatrix ( h, bb(1:2) )
+    atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
 
-  call set_dmatrix ( d, bb(1:2) )
+    allocate ( xcoeff(patch)%coeff(ni,ni) )
+
+    xcoeff(patch)%coeff = atmp
+
+    savedx(patch) = .true.
+
+    deallocate ( atmp, d, b, h )
+
+  end if
 
-  a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+  a => xcoeff(patch)%coeff
 
-  idel = epsdisl /  ( 256.0_wp * delta(1) )
+  idel = epsilon /  ( 256.0_wp * delta(1) )
 
   if ( bb(1) == 0 ) then
 
@@ -206,28 +231,39 @@ subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, &
 
   end do
 
-  deallocate ( a, d, b, h, tmp )
-
   if ( zero_derivs_y == 0 ) then
-    allocate ( a(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj), tmp(nj,nj) )
-    a = zero; d = zero; b = zero; h = zero;
+    if ( .not. savedy(patch) ) then
+      allocate ( atmp(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj) )
+      atmp = zero; d = zero; b = zero; h = zero
   
-    call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl )
+      call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl(2) )
     
-    call set_hmatrix ( h, bb(3:4) )
-  
-    call set_dmatrix ( d, bb(3:4) )
+      call set_hmatrix ( h, bb(3:4) )
   
-    a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+      call set_dmatrix ( d, bb(3:4) )
   
-    idel = epsdisl /  ( 256.0_wp * delta(2) )
+      atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+
+      allocate ( ycoeff(patch)%coeff(nj,nj) )
+
+      ycoeff(patch)%coeff = atmp
+ 
+      savedy(patch) = .true.
+
+      deallocate ( atmp, d, b, h )
+
+    end if
+ 
+    a => ycoeff(patch)%coeff
+ 
+    idel = epsilon /  ( 256.0_wp * delta(2) )
    
-     if ( bb(3) == 0 ) then
-  
-       jl = 1 + gsize(2)
+    if ( bb(3) == 0 ) then
   
-     else
+      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,:) + &
@@ -365,11 +401,11 @@ subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, &
                         a(nj,nj) * var(:,nj,:) ) * idel
   
       jr = nj - 7
-  
+   
     end if
-  
+
     do j = jl, jr
-  
+
       rhs(:,j,:) = rhs(:,j,:) + &
                      ( a(j-4,j) * var(:,j-4,:) + &
                        a(j-3,j) * var(:,j-3,:) + &
@@ -383,22 +419,34 @@ subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, &
   
     end do
   
-    deallocate ( a, d, b, h, tmp )
   end if
-
+ 
   if ( zero_derivs_z == 0 ) then
-    allocate ( a(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk), tmp(nk,nk) )
-    a = zero; d = zero; b = zero; h = zero; tmp = zero
-  
-    call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl )
-  
-    call set_hmatrix ( h, bb(5:6) )
-  
-    call set_dmatrix ( d, bb(5:6) )
+    if ( .not. savedz(patch) ) then
+      allocate ( atmp(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk) )
+      atmp = zero; d = zero; b = zero; h = zero
   
-    a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+      call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl(3) )
+
+      call set_hmatrix ( h, bb(5:6) )
+
+      call set_dmatrix ( d, bb(5:6) )
+
+      atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) )
+
+      allocate ( zcoeff(patch)%coeff(nk,nk) )
+
+      zcoeff(patch)%coeff = atmp
+
+      savedz(patch) = .true.
+
+      deallocate ( atmp, d, b, h )
+
+    end if
   
-    idel = epsdisl /  ( 256.0_wp * delta(3) )
+    a => zcoeff(patch)%coeff
+
+    idel = epsilon /  ( 256.0_wp * delta(3) )
   
     if ( bb(5) == 0 ) then
   
@@ -561,7 +609,6 @@ subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, &
   
     end do
   
-    deallocate ( a, d, b, h, tmp )
   end if
 
 contains
@@ -605,7 +652,7 @@ contains
 
     implicit none
 
-    CCTK_REAL, dimension(:,:), intent(out) :: b
+    CCTK_REAL, dimension(:,:), intent(inout) :: b
     CCTK_INT, dimension(2), intent(in) :: bb
     CCTK_INT, intent(in) :: lsh, gsh, lbnd
     CCTK_REAL, intent(in) :: h, dfl
diff --git a/src/dissipation.c b/src/dissipation.c
index 48f7ba4..8f85d32 100644
--- a/src/dissipation.c
+++ b/src/dissipation.c
@@ -6,6 +6,7 @@
 #include "cctk.h"
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
+#include "cctk_Functions.h"
 
 #include "util_ErrorCodes.h"
 #include "util_Table.h"
@@ -19,6 +20,8 @@ void CCTK_FCALL CCTK_FNAME(dissipation_4_3_opt) (const CCTK_REAL *var,
                                              const CCTK_REAL *dx,
                                              const CCTK_REAL *epsdis,
                                              const CCTK_REAL *diss_fraction,
+                                             const CCTK_INT *npatches,
+                                             const CCTK_INT *patch,
                                              CCTK_REAL *rhs);
 void CCTK_FCALL CCTK_FNAME(dissipation_4_3_alt) (const CCTK_REAL *var,
                                              const CCTK_INT *lsh,
@@ -29,6 +32,8 @@ void CCTK_FCALL CCTK_FNAME(dissipation_4_3_alt) (const CCTK_REAL *var,
                                              const CCTK_REAL *dx,
                                              const CCTK_REAL *epsdis,
                                              const CCTK_REAL *diss_fraction,
+                                             const CCTK_INT *npatches,
+                                             const CCTK_INT *patch,
                                              CCTK_REAL *rhs);
 void CCTK_FCALL CCTK_FNAME(dissipation_6_5_opt) (const CCTK_REAL *var,
                                              const CCTK_INT *ni,
@@ -39,6 +44,8 @@ void CCTK_FCALL CCTK_FNAME(dissipation_6_5_opt) (const CCTK_REAL *var,
                                              const CCTK_REAL *dx,
                                              const CCTK_REAL *epsdis,
                                              const CCTK_REAL *diss_fraction,
+                                             const CCTK_INT *npatches,
+                                             const CCTK_INT *patch,
                                              CCTK_REAL *rhs);
 void CCTK_FCALL CCTK_FNAME(dissipation_6_5_alt) (const CCTK_REAL *var,
                                              const CCTK_INT *ni,
@@ -49,6 +56,8 @@ void CCTK_FCALL CCTK_FNAME(dissipation_6_5_alt) (const CCTK_REAL *var,
                                              const CCTK_REAL *dx,
                                              const CCTK_REAL *epsdis,
                                              const CCTK_REAL *diss_fraction,
+                                             const CCTK_INT *npatches,
+                                             const CCTK_INT *patch,
                                              CCTK_REAL *rhs);
 void CCTK_FCALL CCTK_FNAME(dissipation_2_1) (const CCTK_REAL *var,
                                              const CCTK_INT *ni,
@@ -153,11 +162,12 @@ apply (int const varindex, char const * const optstring, void * const arg)
   int ierr;
   CCTK_INT symtable, pen_sym_handle;
   CCTK_INT symbnd[6];
+  CCTK_INT npatches, patch;
   
   assert (varindex >= 0);
   
   for (d=0; d<3; ++d) {
-    dx[d] = CCTK_DELTA_SPACE(d);
+    dx[d] = CCTK_DELTA_SPACE(d)*h_scaling[d];
     gsize[d] = cctk_nghostzones[d];
   }
   
@@ -280,28 +290,30 @@ apply (int const varindex, char const * const optstring, void * const arg)
       assert(0);
     }
   } else {
+    npatches = MultiPatch_GetMaps ( cctkGH );
+    patch = MultiPatch_GetMap ( cctkGH ) + 1;
     switch(order) {
     case 4: {
       if ( CCTK_Equals(dissipation_type,"Mattson-Svard-Nordstrom") ) {
       CCTK_FNAME(dissipation_4_3_opt)
-        (varptr, cctk_lsh, cctk_gsh, cctk_lbnd,
-                            bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr);
+        (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, bbox, gsize,
+                      dx, &epsdis, diss_fraction, &npatches, &patch, rhsptr);
       } else {
       CCTK_FNAME(dissipation_4_3_alt)
-        (varptr, cctk_lsh, cctk_gsh, cctk_lbnd,
-                            bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr);
+        (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, bbox, gsize,
+                      dx, &epsdis, diss_fraction, &npatches, &patch, rhsptr);
       }
       break;
     }
     case 6: {
       if ( CCTK_Equals(dissipation_type,"Mattson-Svard-Nordstrom") ) {
         CCTK_FNAME(dissipation_6_5_opt)
-          (varptr, cctk_lsh, cctk_gsh, cctk_lbnd,
-                            bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr);
+          (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, bbox, gsize,
+                        dx, &epsdis, diss_fraction, &npatches, &patch, rhsptr);
       } else {
         CCTK_FNAME(dissipation_6_5_alt)
-          (varptr, cctk_lsh, cctk_gsh, cctk_lbnd,
-                            bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr);
+          (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, bbox, gsize,
+                        dx, &epsdis, diss_fraction, &npatches, &patch, rhsptr);
       }
       break;
     }
diff --git a/src/dissipation_coeff.F90 b/src/dissipation_coeff.F90
new file mode 100644
index 0000000..dfd3d94
--- /dev/null
+++ b/src/dissipation_coeff.F90
@@ -0,0 +1,17 @@
+! Arrays to store dissipation operator coefficients for the full restricted norm case.
+! $Header$
+
+#include "cctk.h"
+
+module dissipation_coeff
+
+  type coefficients
+    CCTK_REAL, dimension(:,:), pointer :: coeff
+  end type coefficients
+
+  type(coefficients), dimension(:), pointer :: xcoeff, ycoeff, zcoeff
+
+  logical :: first = .true.
+  logical, dimension(:), allocatable :: savedx, savedy, savedz
+
+end module dissipation_coeff
diff --git a/src/make.code.defn b/src/make.code.defn
index 1519773..0371d01 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -31,6 +31,7 @@ SRCS =  call_derivs.c \
         Dissipation_4_3_min_err_coeff_alt.F90 \
         Dissipation_6_5_min_err_coeff.F90 \
         Dissipation_6_5_min_err_coeff_alt.F90 \
+        dissipation_coeff.F90 \
         get_coeffs.c \
         Coefficients_2_1.F90 \
         Coefficients_4_2.F90 \