Add new parameters derivatives_zero_y and derivatives_zero_z in order to

ecplictly setting any derivatives in the y- and/or z-direction to zero
when simulating 1- or 2-D problems with a 3D code. This avoids issues with
errors at roundoff level in the y- and z-directions. Should be used with care
and together with the parameter check_grid_sizes to avoid complaints about
too small grids in the y- and z-directions.

Has only been tested in a few cases. If problems are caused by this commit,
please report them immediately.


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

1 files changed, 70 insertions, 66 deletions

diff --git a/src/Dissipation_2_1.F90 b/src/Dissipation_2_1.F90
index 0904b33..0b625f4 100644
--- a/src/Dissipation_2_1.F90
+++ b/src/Dissipation_2_1.F90
@@ -60,75 +60,79 @@ subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs)
                              ( var(il-1:ir-1,:,:) + &
                                var(il+1:ir+1,:,:) ) ) * idel
 
-  call set_coeff ( a )
-
-  if ( scale_with_h > 0 ) then
-    idel = epsilon / ( 4 * delta(2) )
-  else
-    idel = epsilon / 4
-  end if
-
-  if ( bb(3) == 0 ) then
-    jl = 1 + gsize(2)
-  else
-    rhs(:,1,:) = rhs(:,1,:) + &
-                 ( a(1,1) * var(:,1,:) + a(2,1) * var(:,2,:) ) * idel
-    rhs(:,2,:) = rhs(:,2,:) + &
-                 ( a(1,2) * var(:,1,:) + a(2,2) * var(:,2,:) + &
-                   a(3,2) * var(:,3,:) ) * 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,:) ) * idel
-    rhs(:,nj,:)   = rhs(:,nj,:) + &
-                 ( a(1,1) * var(:,nj,:) + a(2,1) * var(:,nj-1,:) ) * idel
-
-    jr = nj - 2
+  if ( zero_derivs_y == 0 ) then
+    call set_coeff ( a )
+  
+    if ( scale_with_h > 0 ) then
+      idel = epsilon / ( 4 * delta(2) )
+    else
+      idel = epsilon / 4
+    end if
+  
+    if ( bb(3) == 0 ) then
+      jl = 1 + gsize(2)
+    else
+      rhs(:,1,:) = rhs(:,1,:) + &
+                   ( a(1,1) * var(:,1,:) + a(2,1) * var(:,2,:) ) * idel
+      rhs(:,2,:) = rhs(:,2,:) + &
+                   ( a(1,2) * var(:,1,:) + a(2,2) * var(:,2,:) + &
+                     a(3,2) * var(:,3,:) ) * 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,:) ) * idel
+      rhs(:,nj,:)   = rhs(:,nj,:) + &
+                   ( a(1,1) * var(:,nj,:) + a(2,1) * var(:,nj-1,:) ) * idel
+  
+      jr = nj - 2
+    end if
+    rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+                     ( -2.0_wp * var(:,jl:jr,:) + &
+                               ( var(:,jl-1:jr-1,:) + &
+                                 var(:,jl+1:jr+1,:) ) ) * idel
   end if
-  rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
-                   ( -2.0_wp * var(:,jl:jr,:) + &
-                             ( var(:,jl-1:jr-1,:) + &
-                               var(:,jl+1:jr+1,:) ) ) * idel
-
-  call set_coeff ( a )
-
-  if ( scale_with_h > 0 ) then
-    idel = epsilon / ( 4 * delta(3) )
-  else
-    idel = epsilon / 4
-  end if
-
-  if ( bb(5) == 0 ) then
-    kl = 1 + gsize(3)
-  else
-    rhs(:,:,1) = rhs(:,:,1) + &
-                 ( a(1,1) * var(:,:,1) + a(2,1) * var(:,:,2) ) * idel
-    rhs(:,:,2) = rhs(:,:,2) + &
-                 ( a(1,2) * var(:,:,1) + a(2,2) * var(:,:,2) + &
-                   a(3,2) * var(:,:,3) ) * 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) ) * idel
-    rhs(:,:,nk)   = rhs(:,:,nk) + &
-                 ( a(1,1) * var(:,:,nk) + a(2,1) * var(:,:,nk-1) ) * idel
 
-    kr = nk - 2
+  if ( zero_derivs_z == 0 ) then
+    call set_coeff ( a )
+  
+    if ( scale_with_h > 0 ) then
+      idel = epsilon / ( 4 * delta(3) )
+    else
+      idel = epsilon / 4
+    end if
+  
+    if ( bb(5) == 0 ) then
+      kl = 1 + gsize(3)
+    else
+      rhs(:,:,1) = rhs(:,:,1) + &
+                   ( a(1,1) * var(:,:,1) + a(2,1) * var(:,:,2) ) * idel
+      rhs(:,:,2) = rhs(:,:,2) + &
+                   ( a(1,2) * var(:,:,1) + a(2,2) * var(:,:,2) + &
+                     a(3,2) * var(:,:,3) ) * 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) ) * idel
+      rhs(:,:,nk)   = rhs(:,:,nk) + &
+                   ( a(1,1) * var(:,:,nk) + a(2,1) * var(:,:,nk-1) ) * idel
+  
+      kr = nk - 2
+    end if
+    rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+                     ( -2.0_wp * var(:,:,kl:kr) + &
+                               ( var(:,:,kl-1:kr-1) + &
+                                 var(:,:,kl+1:kr+1) ) ) * idel
   end if
-  rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
-                   ( -2.0_wp * var(:,:,kl:kr) + &
-                             ( var(:,:,kl-1:kr-1) + &
-                               var(:,:,kl+1:kr+1) ) ) * idel
 
 contains