Writing and Reading variables of type 'complex' should now work. Please test.

darcs-hash:20040412205904-19929-21858f81ba0fbfc5618ec47c5083fc3b146ddaa8.gz

4 files changed, 74 insertions, 16 deletions

diff --git a/Carpet/CarpetIOHDF5/src/iohdf5.cc b/Carpet/CarpetIOHDF5/src/iohdf5.cc
index 83891a809..35f0ae08a 100644
--- a/Carpet/CarpetIOHDF5/src/iohdf5.cc
+++ b/Carpet/CarpetIOHDF5/src/iohdf5.cc
@@ -17,7 +17,7 @@
 #include "cctk_Parameters.h"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.cc,v 1.26 2004/04/07 14:10:23 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.cc,v 1.27 2004/04/12 22:59:04 cott Exp $";
   CCTK_FILEVERSION(Carpet_CarpetIOHDF5_iohdf5_cc);
 }
 
@@ -126,11 +126,41 @@ namespace CarpetIOHDF5 {
 
     myGH->open_output_files = NULL;
 
+    // Now set hdf5 complex datatypes
+    // Stolen from Thomas Radke
+    HDF5_ERROR (myGH->HDF5_COMPLEX =
+              H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX)));
+    HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX, "real",
+                         offsetof (CCTK_COMPLEX, Re), HDF5_REAL));
+    HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX, "imag",
+                         offsetof (CCTK_COMPLEX, Im), HDF5_REAL));
+#ifdef CCTK_REAL4
+    HDF5_ERROR (myGH->HDF5_COMPLEX8 =
+		H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX8)));
+    HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX8, "real",
+			   offsetof (CCTK_COMPLEX8, Re), HDF5_REAL4));
+    HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX8, "imag",
+			   offsetof (CCTK_COMPLEX8, Im), HDF5_REAL4));
+#endif
+#ifdef CCTK_REAL8
+    HDF5_ERROR (myGH->HDF5_COMPLEX16 =
+		H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX16)));
+    HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX16, "real",
+			   offsetof (CCTK_COMPLEX16, Re), HDF5_REAL8));
+    HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX16, "imag",
+			   offsetof (CCTK_COMPLEX16, Im), HDF5_REAL8));
+#endif
+#ifdef CCTK_REAL16
+    HDF5_ERROR (myGH->HDF5_COMPLEX32 =
+		H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX32)));
+    HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX32, "real",
+			   offsetof (CCTK_COMPLEX32, Re), HDF5_REAL16));
+    HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX32, "imag",
+			   offsetof (CCTK_COMPLEX32, Im), HDF5_REAL16));
+#endif
     return (myGH);
   }
   
-  
-  
   int OutputGH (const cGH* const cctkGH) {
     for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) {
       if (TimeToOutput(cctkGH, vindex)) {
@@ -355,7 +385,7 @@ namespace CarpetIOHDF5 {
 
 	      // Select datatype
 
-	      const hid_t datatype = h5DataType(CCTK_VarTypeI(n));
+	      const hid_t datatype = h5DataType(cctkGH,CCTK_VarTypeI(n));
           
 	      ostringstream datasetnamebuf;
 	      datasetnamebuf << varname
@@ -630,7 +660,7 @@ namespace CarpetIOHDF5 {
     int recovery_rl = -1;
     int recovery_comp = -1;
 
-    CCTK_REAL *h5data;
+    void * h5data;
 
     // Check for storage
     if (! CCTK_QueryGroupStorageI(cctkGH, group)) {
@@ -685,12 +715,16 @@ namespace CarpetIOHDF5 {
 	   }
 
 	 }
-	 const hid_t datatype = H5T_NATIVE_DOUBLE;
+         const int cctkDataType = CCTK_VarTypeI(n);
+	 const hid_t datatype = h5DataType(cctkGH,cctkDataType);
 	 
          //cout << "datalength: " << datalength << " rank: " << rank << "\n";
          //cout << shape[0] << " " << shape[1] << " " << shape[2] << "\n";
 	 
-	 h5data = (CCTK_REAL*) malloc(sizeof(double)*datalength);
+	 // to do: read in an allocate with correct datatype
+
+	 h5data = (void*) malloc( CCTK_VarTypeSize(cctkDataType) *datalength);
+
 	 herr = H5Dread(dataset,datatype,H5S_ALL, H5S_ALL, H5P_DEFAULT,(void*)h5data);
 	 assert(!herr);
 	 
diff --git a/Carpet/CarpetIOHDF5/src/iohdf5.hh b/Carpet/CarpetIOHDF5/src/iohdf5.hh
index 8b2fa6eb0..21e5b512f 100644
--- a/Carpet/CarpetIOHDF5/src/iohdf5.hh
+++ b/Carpet/CarpetIOHDF5/src/iohdf5.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.hh,v 1.9 2004/04/03 12:40:21 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.hh,v 1.10 2004/04/12 22:59:04 cott Exp $
 
 #ifndef CARPETIOHDF5_HH
 #define CARPETIOHDF5_HH
@@ -52,7 +52,6 @@
 
 /*** Define the different datatypes used for HDF5 I/O
      NOTE: the complex datatype SHOULD be [is] defined dynamically at runtime in Startup.c
-
      100% of the definitions below were taken from Thomas Radke's IOHDF5Util thorn for PUGH
  ***/
 /* char type is easy */
@@ -118,6 +117,23 @@
 #define HDF5_INT    HDF5_INT1
 #endif
 
+/* Nice error handling. Stolen from Thomas Radke */
+/* check return code of HDF5 call and print a warning in case of an error */
+#define HDF5_ERROR(fn_call)                                                   \
+        {                                                                     \
+          int _error_code = fn_call;                                          \
+                                                                              \
+                                                                              \
+          if (_error_code < 0)                                                \
+          {                                                                   \
+            CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,              \
+                        "HDF5 call '%s' returned error code %d",              \
+                        #fn_call, _error_code);                               \
+          }                                                                   \
+        }
+
+
+
 
 namespace CarpetIOHDF5 {
   
@@ -180,7 +196,7 @@ namespace CarpetIOHDF5 {
   int GetnDatasets (const hid_t reader);
   void GetDatasetName (const hid_t reader, const int _index, char* name);
 
-  hid_t h5DataType(int cctk_type);
+  hid_t h5DataType(const cGH* const cctkGH, int cctk_type);
 
 } // namespace CarpetIOHDF5
 
diff --git a/Carpet/CarpetIOHDF5/src/iohdf5GH.h b/Carpet/CarpetIOHDF5/src/iohdf5GH.h
index 21463ac44..0516201c2 100644
--- a/Carpet/CarpetIOHDF5/src/iohdf5GH.h
+++ b/Carpet/CarpetIOHDF5/src/iohdf5GH.h
@@ -46,6 +46,10 @@ typedef struct
   /* iteration number of the last checkpoint */
   int last_checkpoint_iteration;
 
+  /* hdf5 datatype for stupid complex variables; to be set at run time */
+  hid_t HDF5_COMPLEX, HDF5_COMPLEX8, HDF5_COMPLEX16, HDF5_COMPLEX32;
+  
+
 } CarpetIOHDF5GH;
 
 
diff --git a/Carpet/CarpetIOHDF5/src/iohdf5utils.cc b/Carpet/CarpetIOHDF5/src/iohdf5utils.cc
index 49fecffe1..085870743 100644
--- a/Carpet/CarpetIOHDF5/src/iohdf5utils.cc
+++ b/Carpet/CarpetIOHDF5/src/iohdf5utils.cc
@@ -17,7 +17,7 @@
 #include "cctk_Parameters.h"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5utils.cc,v 1.5 2004/04/03 12:40:21 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5utils.cc,v 1.6 2004/04/12 22:59:04 cott Exp $";
   CCTK_FILEVERSION(Carpet_CarpetIOHDF5_iohdf5utils_cc);
 }
 
@@ -33,6 +33,7 @@ extern "C" {
 #include "carpet.hh"
 
 #include "iohdf5.hh"
+#include "iohdf5GH.h"
 
 
 
@@ -460,9 +461,12 @@ namespace CarpetIOHDF5 {
 		     &getn)<0){}
   }
 
-  hid_t h5DataType (int cctk_type) {
+  hid_t h5DataType (const cGH* const cctkGH, int cctk_type) {
 
     hid_t retval;  
+
+    CarpetIOHDF5GH *myGH;
+    myGH = (CarpetIOHDF5GH *) CCTK_GHExtension (cctkGH, "CarpetIOHDF5");
   
     // this is adapted from Thomas Radke's IOHDF5Util. Thanks, Thomas!
 
@@ -471,7 +475,7 @@ namespace CarpetIOHDF5 {
       case CCTK_VARIABLE_CHAR:      retval = HDF5_CHAR; break;
       case CCTK_VARIABLE_INT:       retval = HDF5_INT; break;
       case CCTK_VARIABLE_REAL:      retval = HDF5_REAL; break;
-	//      case CCTK_VARIABLE_COMPLEX:   retval = myGH->HDF5_COMPLEX; break;
+      case CCTK_VARIABLE_COMPLEX:   retval = myGH->HDF5_COMPLEX; break;
 #ifdef CCTK_INT1
       case CCTK_VARIABLE_INT1:      retval = HDF5_INT1; break;
 #endif
@@ -486,15 +490,15 @@ namespace CarpetIOHDF5 {
 #endif
 #ifdef CCTK_REAL4
       case CCTK_VARIABLE_REAL4:     retval = HDF5_REAL4; break;
-	//      case CCTK_VARIABLE_COMPLEX8:  retval = myGH->HDF5_COMPLEX8; break;
+      case CCTK_VARIABLE_COMPLEX8:  retval = myGH->HDF5_COMPLEX8; break;
 #endif
 #ifdef CCTK_REAL8
       case CCTK_VARIABLE_REAL8:     retval = HDF5_REAL8; break;
-	//case CCTK_VARIABLE_COMPLEX16: retval = myGH->HDF5_COMPLEX16; break;
+      case CCTK_VARIABLE_COMPLEX16: retval = myGH->HDF5_COMPLEX16; break;
 #endif
 #ifdef CCTK_REAL16
       case CCTK_VARIABLE_REAL16:    retval = HDF5_REAL16; break;
-	//case CCTK_VARIABLE_COMPLEX32: retval = myGH->HDF5_COMPLEX32; break;
+      case CCTK_VARIABLE_COMPLEX32: retval = myGH->HDF5_COMPLEX32; break;
 #endif
   
       default: CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,