Add 5th order spatial interpolation.

darcs-hash:20030624120058-07bb3-01d8634b9ffbc710b6102af4d3128f32b73f21d2.gz

5 files changed, 67 insertions, 22 deletions

diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 0f76c7830..155e220db 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.25 2003/02/11 10:54:16 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.26 2003/06/24 14:00:58 schnetter Exp $
 
 #include <assert.h>
 
@@ -314,6 +314,14 @@ extern "C" {
      const int srcbbox[3][3],
      const int dstbbox[3][3],
      const int regbbox[3][3]);
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_o5)
+    (const CCTK_REAL8* src,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
   
   void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl)
     (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
@@ -333,6 +341,15 @@ extern "C" {
      const int srcbbox[3][3],
      const int dstbbox[3][3],
      const int regbbox[3][3]);
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o5)
+    (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
+     const CCTK_REAL8* src2, const CCTK_REAL8& t2,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst, const CCTK_REAL8& t,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
   
   void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl)
     (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
@@ -354,6 +371,16 @@ extern "C" {
      const int srcbbox[3][3],
      const int dstbbox[3][3],
      const int regbbox[3][3]);
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_o5)
+    (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
+     const CCTK_REAL8* src2, const CCTK_REAL8& t2,
+     const CCTK_REAL8* src3, const CCTK_REAL8& t3,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst, const CCTK_REAL8& t,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
   
 }
 
@@ -444,6 +471,15 @@ void data<CCTK_REAL8,3>
 	   dstshp[0], dstshp[1], dstshp[2],
 	   srcbbox, dstbbox, regbbox);
 	break;
+      case 4:
+      case 5:
+	CCTK_FNAME(prolongate_3d_real8_o5)
+	  ((const CCTK_REAL8*)srcs[0]->storage(),
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(),
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
       default:
 	assert (0);
       }
@@ -472,6 +508,16 @@ void data<CCTK_REAL8,3>
 	   dstshp[0], dstshp[1], dstshp[2],
 	   srcbbox, dstbbox, regbbox);
 	break;
+      case 4:
+      case 5:
+	CCTK_FNAME(prolongate_3d_real8_2tl_o5)
+	  ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+	   (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(), time,
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
       default:
 	assert (0);
       }
@@ -502,6 +548,17 @@ void data<CCTK_REAL8,3>
 	   dstshp[0], dstshp[1], dstshp[2],
 	   srcbbox, dstbbox, regbbox);
 	break;
+      case 4:
+      case 5:
+	CCTK_FNAME(prolongate_3d_real8_3tl_o5)
+	  ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+	   (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+	   (const CCTK_REAL8*)srcs[2]->storage(), times[2],
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(), time,
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
       default:
 	assert (0);
       }
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index c6e1b111e..6459559fa 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -1,5 +1,5 @@
 # Main make.code.defn file for thorn CarpetLib   -*-Makefile-*-
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.7 2003/01/03 15:49:36 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.8 2003/06/24 14:00:58 schnetter Exp $
 
 # Source files in this directory
 SRCS =	bbox.cc					\
@@ -17,10 +17,13 @@ SRCS =	bbox.cc					\
 	copy_3d_real8.F77			\
 	prolongate_3d_real8.F77			\
 	prolongate_3d_real8_o3.F77		\
+	prolongate_3d_real8_o5.F77		\
 	prolongate_3d_real8_2tl.F77		\
 	prolongate_3d_real8_2tl_o3.F77		\
+	prolongate_3d_real8_2tl_o5.F77		\
 	prolongate_3d_real8_3tl.F77		\
 	prolongate_3d_real8_3tl_o3.F77		\
+	prolongate_3d_real8_3tl_o5.F77		\
 	restrict_3d_real8.F77
 
 # Subdirectories containing source files
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o5.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o5.F77
index e28354048..0eb705b88 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o5.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o5.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o5.F77,v 1.3 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o5.F77,v 1.1 2003/06/24 14:00:58 schnetter Exp $
       
 #include "cctk.h"
 
@@ -88,11 +88,6 @@ c     bbox(:,3) is stride
 c     This could be handled, but is likely to point to an error elsewhere
             call CCTK_WARN (0, "Internal error: region extent is empty")
          end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
          regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1
          dstkfac = srcbbox(d,3) / dstbbox(d,3)
          srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)
@@ -152,7 +147,7 @@ c     Linear (first order) interpolation
          call CCTK_WARN (0, "Internal error: arrays have same time")
       end if
       if (t.lt.min(t1,t2)-eps .or. t.gt.max(t1,t2)+eps) then
-         call CCTK_WARN (0, "Internal error: extrapolation in time")
+         call CCTK_WARN (0, "Internal error: extrapolation")
       end if
       
       s1fac = (t - t2) / (t1 - t2)
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o5.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o5.F77
index 28f8c155f..e664ac47e 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o5.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o5.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o5.F77,v 1.3 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o5.F77,v 1.1 2003/06/24 14:00:58 schnetter Exp $
       
 #include "cctk.h"
 
@@ -90,11 +90,6 @@ c     bbox(:,3) is stride
 c     This could be handled, but is likely to point to an error elsewhere
             call CCTK_WARN (0, "Internal error: region extent is empty")
          end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
          regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1
          dstkfac = srcbbox(d,3) / dstbbox(d,3)
          srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)
@@ -154,7 +149,7 @@ c     Quadratic (second order) interpolation
          call CCTK_WARN (0, "Internal error: arrays have same time")
       end if
       if (t.lt.min(t1,t2,t3)-eps .or. t.gt.max(t1,t2,t3)+eps) then
-         call CCTK_WARN (0, "Internal error: extrapolation in time")
+         call CCTK_WARN (0, "Internal error: extrapolation")
       end if
       
       s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3))
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o5.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o5.F77
index 8c86178ad..ce9e4f5c5 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_o5.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o5.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o5.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o5.F77,v 1.1 2003/06/24 14:00:58 schnetter Exp $
       
 #include "cctk.h"
 
@@ -82,11 +82,6 @@ c     bbox(:,3) is stride
 c     This could be handled, but is likely to point to an error elsewhere
             call CCTK_WARN (0, "Internal error: region extent is empty")
          end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
          regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1
          dstkfac = srcbbox(d,3) / dstbbox(d,3)
          srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)