WaveToyMoL: Support 4th order differencing

darcs-hash:20050508171319-891bb-48b570b6edcaa6409e931ef41afbc6049bcbb039.gz

4 files changed, 36 insertions, 13 deletions

diff --git a/CarpetExtra/WaveToyMoL/param.ccl b/CarpetExtra/WaveToyMoL/param.ccl
index fab76bcf2..bb779cca6 100644
--- a/CarpetExtra/WaveToyMoL/param.ccl
+++ b/CarpetExtra/WaveToyMoL/param.ccl
@@ -21,3 +21,9 @@ STRING bound "Type of boundary condition to use"
 BOOLEAN estimate_error "Estimate the truncation error"
 {
 } "no"
+
+INT order "Finite differencing order"
+{
+  2 :: "second order"
+  4 :: "fourth order"
+} 2
diff --git a/CarpetExtra/WaveToyMoL/schedule.ccl b/CarpetExtra/WaveToyMoL/schedule.ccl
index 0122f93a8..c8b81846f 100644
--- a/CarpetExtra/WaveToyMoL/schedule.ccl
+++ b/CarpetExtra/WaveToyMoL/schedule.ccl
@@ -3,10 +3,6 @@
 STORAGE: scalarevolve[3]
 STORAGE: scalarevolvedot
 
-if (estimate_error) {
-  STORAGE: scalarevolveerrorestimate
-}
-
 
 
 SCHEDULE WaveToyMol_Startup AT startup
@@ -47,8 +43,11 @@ SCHEDULE GROUP ApplyBCs IN MoL_PostStep AFTER WaveToyMoL_Boundaries
 
 
 
-if (estimate_error) {
+if (estimate_error)
+{
   
+  STORAGE: scalarevolveerrorestimate
+
   SCHEDULE WaveToyMoL_EstimateError AT poststep
   {
     LANG: Fortran
diff --git a/CarpetExtra/WaveToyMoL/src/calcrhs.F77 b/CarpetExtra/WaveToyMoL/src/calcrhs.F77
index a2fbe5651..371278970 100644
--- a/CarpetExtra/WaveToyMoL/src/calcrhs.F77
+++ b/CarpetExtra/WaveToyMoL/src/calcrhs.F77
@@ -13,10 +13,13 @@
       integer imin(3), imax(3)
       CCTK_REAL dx(3)
       integer i, j, k
+      logical isnearboundary
       integer d
       do d=1,3
-         imin(d) = 1+bndwidth
-         imax(d) = cctk_lsh(d)-bndwidth
+         imin(d) = 1+cctk_nghostzones(d)
+         imax(d) = cctk_lsh(d)-cctk_nghostzones(d)
+         if (cctk_bbox(2*d-1).ne.0) imin(d) = 1+bndwidth
+         if (cctk_bbox(2*d  ).ne.0) imax(d) = cctk_lsh(d)-bndwidth
       end do
       do d=1,3
          dx(d) = CCTK_DELTA_SPACE(3)
@@ -25,10 +28,23 @@
          do j=imin(2),imax(2)
             do i=imin(1),imax(1)
                phidot(i,j,k) = psi(i,j,k)
-               psidot(i,j,k) =
-     $              ( (phi(i-1,j,k) - 2*phi(i,j,k) + phi(i+1,j,k)) / dx(1)**2
-     $              + (phi(i,j-1,k) - 2*phi(i,j,k) + phi(i,j+1,k)) / dx(2)**2
-     $              + (phi(i,j,k-1) - 2*phi(i,j,k) + phi(i,j,k+1)) / dx(3)**2)
+               isnearboundary = .false.
+               isnearboundary = isnearboundary .or. i.eq.imin(1) .or. i.eq.imax(1)
+               isnearboundary = isnearboundary .or. j.eq.imin(2) .or. j.eq.imax(2)
+               isnearboundary = isnearboundary .or. k.eq.imin(3) .or. k.eq.imax(3)
+               if (order.eq.2 .or. isnearboundary) then
+                  psidot(i,j,k) = (
+     $                 + (phi(i-1,j,k) - 2*phi(i,j,k) + phi(i+1,j,k)) / dx(1)**2
+     $                 + (phi(i,j-1,k) - 2*phi(i,j,k) + phi(i,j+1,k)) / dx(2)**2
+     $                 + (phi(i,j,k-1) - 2*phi(i,j,k) + phi(i,j,k+1)) / dx(3)**2)
+               else if (order.eq.4) then
+                  psidot(i,j,k) = (
+     $                 + (- phi(i-2,j,k) + 16*phi(i-1,j,k) - 30*phi(i,j,k) + 16*phi(i+1,j,k) - phi(i+2,j,k)) / (12 * dx(1)**2)
+     $                 + (- phi(i,j-2,k) + 16*phi(i,j-1,k) - 30*phi(i,j,k) + 16*phi(i,j+1,k) - phi(i,j+2,k)) / (12 * dx(2)**2)
+     $                 + (- phi(i,j,k-2) + 16*phi(i,j,k-1) - 30*phi(i,j,k) + 16*phi(i,j,k+1) - phi(i,j,k+2)) / (12 * dx(3)**2))
+               else
+                  call CCTK_WARN (0, "internal error")
+               end if
             end do
          end do
       end do
diff --git a/CarpetExtra/WaveToyMoL/src/estimateerror.F77 b/CarpetExtra/WaveToyMoL/src/estimateerror.F77
index f046183be..ae8d3f0c9 100644
--- a/CarpetExtra/WaveToyMoL/src/estimateerror.F77
+++ b/CarpetExtra/WaveToyMoL/src/estimateerror.F77
@@ -15,8 +15,10 @@
       integer i, j, k
       integer d
       do d=1,3
-         imin(d) = 1+bndwidth
-         imax(d) = cctk_lsh(d)-bndwidth
+         imin(d) = 1+cctk_nghostzones(d)
+         imax(d) = cctk_lsh(d)-cctk_nghostzones(d)
+         if (cctk_bbox(2*d-1).ne.0) imin(d) = 1+bndwidth
+         if (cctk_bbox(2*d  ).ne.0) imax(d) = cctk_lsh(d)-bndwidth
       end do
       do d=1,3
          dx(d) = CCTK_DELTA_SPACE(3)