Add support for reading and writing units

The full dUnits object is not recovered, just the scaling at this point.
The vector itself is redimensionalized. There is short-circuit in case the
stored vector has exactly the same scaling factor as the dFS reading it, so
the current implementation may not have bitwise reproducibility.

include/dohpfs.h

@@ -12,6 +12,7 @@
 
 #include "dohpmesh.h"
 #include "dohpjacobi.h"
+#include "dohpunits.h"
 
 dEXTERN_C_BEGIN
 
@@ -65,6 +66,9 @@ extern dErr dFSSetDegree(dFS,dJacobi,dMeshTag);
 extern dErr dFSSetBlockSize(dFS,dInt);
 extern dErr dFSGetBlockSize(dFS,dInt*);
 extern dErr dFSSetFieldName(dFS,dInt,const char*);
+extern dErr dFSGetFieldName(dFS,dInt,const char**);
+extern dErr dFSSetFieldUnit(dFS fs,dInt fn,dUnit unit);
+extern dErr dFSGetFieldUnit(dFS fs,dInt fn,dUnit *unit);
 extern dErr dFSRegisterBoundary(dFS,dInt,dFSBStatus,dFSConstraintFunction,void*);
 extern dErr dFSRegisterBoundarySet(dFS,dMeshESH,dFSBStatus,dFSConstraintFunction,void*);
 extern dErr dFSSetFromOptions(dFS);

include/dohpfsimpl.h

@@ -126,6 +126,7 @@ struct _p_dFS {
   dFSRotation  rot;             /**< Rotation for local vector */
   char         orderingtype[256];
   char       **fieldname;
+  dUnit       *fieldunit;
   void        *data;
 };
 

include/dohpviewerdhm.h

@@ -91,6 +91,7 @@ typedef struct {
 typedef struct {
   hdset_reg_ref_t fs;
   dReal           time;
+  dht_Units       units;
   dInt            internal_state;
 } dht_Vec;
 

src/fs/impls/cont/contviewdhm.c

@@ -87,8 +87,8 @@ dErr dFSView_Cont_DHM(dFS fs,dViewer viewer)
     err = dMallocA(fs5.fields.len,&field5);dCHK(err);
     for (i=0; i<bs; i++) {
       field5[i].name = fs->fieldname[i];
-      field5[i].units.dimensions = (char*)"m s-1"; /* we only use it as \c const */
-      field5[i].units.scale = exp(1.0);
+      field5[i].units.dimensions = (char*)dUnitName(fs->fieldunit[i]); /* we only use it as \c const */
+      field5[i].units.scale = dUnitDimensionalize(fs->fieldunit[i],1.0);
     }
     fs5.fields.p = field5;
     err = dFSGetDHMLink(fs,viewer,&fsdset,&fsspace);dCHK(err); /* Get location to write this FS */
@@ -133,6 +133,7 @@ static dErr VecView_Dohp_FSCont_DHM(Vec X,PetscViewer viewer)
   const char  *xname;
   dFS          fs;
   Vec          Xclosure;
+  dUnit        fieldunit;
   dScalar     *x;
   hid_t        fsdset,fsspace,curstep,dset,fspace,mspace,plist;
   hsize_t      gdim[2],offset[2],count[2];
@@ -142,8 +143,9 @@ static dErr VecView_Dohp_FSCont_DHM(Vec X,PetscViewer viewer)
   dFunctionBegin;
   err = dViewerDHMSetUp(viewer);dCHK(err);
   err = PetscObjectGetName((PetscObject)X,&xname);dCHK(err);
-  err = PetscObjectQuery((PetscObject)X,"dFS",(PetscObject*)&fs);dCHK(err);
+  err = VecDohpGetFS(X,&fs);dCHK(err);
   if (!fs) dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Vector not generated from a FS");
+  err = dFSGetFieldUnit(fs,0,&fieldunit);dCHK(err);
   err = dFSGetDHMLink(fs,viewer,&fsdset,&fsspace);dCHK(err); /* we are responsible for closing */
   err = dViewerDHMGetStep(viewer,&curstep);dCHK(err);        /* leave curstep open */
   err = dFSView_Cont_DHM(fs,viewer);dCHK(err);
@@ -180,6 +182,8 @@ static dErr VecView_Dohp_FSCont_DHM(Vec X,PetscViewer viewer)
       /* Initialize data vecatt[0] */
       herr = H5Rcreate(&vecatt[0].fs,dhm->file,fsname,H5R_DATASET_REGION,fsspace);dH5CHK(herr,H5Rcreate);
       vecatt[0].time = dhm->time;
+      vecatt[0].units.dimensions = NULL;
+      vecatt[0].units.scale = dUnitDimensionalize(fieldunit,1.0);
       err = PetscObjectStateQuery((PetscObject)X,&vecatt[0].internal_state);dCHK(err);
       /* Write attribute */
       aspace = H5Screate_simple(1,dims,NULL);dH5CHK(aspace,H5Screate_simple);
@@ -206,8 +210,7 @@ dErr VecView_Dohp_FSCont(Vec x,PetscViewer viewer)
   dErr err;
 
   dFunctionBegin;
-  err = PetscObjectQuery((PetscObject)x,"dFS",(PetscObject*)&fs);dCHK(err);
-  if (!fs) dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Vector not generated from a FS");
+  err = VecDohpGetFS(x,&fs);dCHK(err);
   err = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDHM,&isdhm);dCHK(err);
   if (isdhm) {
     err = VecView_Dohp_FSCont_DHM(x,viewer);dCHK(err);
@@ -249,7 +252,7 @@ dErr dFSLoadIntoFS_Cont_DHM(PetscViewer viewer,const char name[],dFS fs)
     herr = H5Sclose(memspace);
 
     if (debug) {
-      dPrintf(PETSC_COMM_SELF,"degree = %s\nglobal_offset = %s\npartition = %s\nordered_subdomain = %s\n",fs5.degree,fs5.global_offset,fs5.partition,fs5.ordered_subdomain);dCHK(err);
+      err = dPrintf(PETSC_COMM_SELF,"degree = %s\nglobal_offset = %s\npartition = %s\nordered_subdomain = %s\n",fs5.degree,fs5.global_offset,fs5.partition,fs5.ordered_subdomain);dCHK(err);
     }
     meshobj = H5Rdereference(dhm->meshroot,H5R_OBJECT,&fs5.mesh);dH5CHK(meshobj,H5Rdereference);
     {
@@ -346,11 +349,15 @@ dErr VecDohpLoadIntoVector(PetscViewer viewer,const char fieldname[],Vec X)
 {
   dErr    err;
   dBool   match;
-  hid_t   curstep,dset,memspace,filespace;
+  dht_Vec vecmeta;
+  hid_t   curstep,dset,memspace,filespace,vattr,vectype;
   hsize_t gdim[2],offset[2],count[2];
   herr_t  herr;
   Vec     Xclosure;
   dScalar *x;
+  dReal   scale;
+  dFS     fs;
+  dUnit   unit;
 
   dFunctionBegin;
   err = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDHM,&match);dCHK(err);
@@ -360,7 +367,11 @@ dErr VecDohpLoadIntoVector(PetscViewer viewer,const char fieldname[],Vec X)
 
   err = dViewerDHMSetUp(viewer);dCHK(err);
   err = dViewerDHMGetStep(viewer,&curstep);dCHK(err);
+  err = dViewerDHMGetVecType(viewer,&vectype);dCHK(err);
   err = dH5Dopen(curstep,fieldname,H5P_DEFAULT,&dset);dCHK(err);
+  err = dH5Aopen(dset,"meta",H5P_DEFAULT,&vattr);dCHK(err);
+  herr = H5Aread(vattr,vectype,&vecmeta);dH5CHK(herr,H5Aread);
+  herr = H5Aclose(vattr);dH5CHK(herr,H5Aclose);
 
   err = VecDohpGetClosure(X,&Xclosure);dCHK(err);
 
@@ -377,6 +388,12 @@ dErr VecDohpLoadIntoVector(PetscViewer viewer,const char fieldname[],Vec X)
   herr = H5Sclose(memspace);dH5CHK(herr,H5Sclose);
   herr = H5Sclose(filespace);dH5CHK(herr,H5Sclose);
   herr = H5Dclose(dset);dH5CHK(herr,H5Dclose);
+
+  err = VecDohpGetFS(X,&fs);dCHK(err);
+  err = dFSGetFieldUnit(fs,0,&unit);dCHK(err);
+  scale = dUnitNonDimensionalize(unit,vecmeta.units.scale);
+  err = VecScale(Xclosure,scale);dCHK(err);
+
   err = VecDohpRestoreClosure(X,&Xclosure);dCHK(err);
   dFunctionReturn(0);
 }

src/fs/interface/fs.c

@@ -103,6 +103,8 @@ dErr dFSSetBlockSize(dFS fs,dInt bs)
   for (dInt i=0; i<obs; i++) {err = dFree(fs->fieldname[i]);dCHK(err);}
   err = dFree(fs->fieldname);dCHK(err);
   err = dCallocA(bs,&fs->fieldname);dCHK(err);
+  err = dFree(fs->fieldunit);dCHK(err);
+  err = dCallocA(bs,&fs->fieldunit);dCHK(err);
   fs->dm.bs = bs;
   dFunctionReturn(0);
 }
@@ -138,6 +140,44 @@ dErr dFSSetFieldName(dFS fs,dInt fn,const char *fname)
   err = PetscStrallocpy(fname,&fs->fieldname[fn]);dCHK(err);
   dFunctionReturn(0);
 }
+dErr dFSGetFieldName(dFS fs,dInt fn,const char **fname)
+{
+  dErr err;
+  dInt bs;
+
+  dFunctionBegin;
+  dValidHeader(fs,DM_CLASSID,1);
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
+  if (fn < 0 || bs <= fn) dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field number %d out of range",fn);
+  *fname = fs->fieldname[fn];
+  dFunctionReturn(0);
+}
+
+// Set the physical units to use for the field
+dErr dFSSetFieldUnit(dFS fs,dInt fn,dUnit unit)
+{
+  dErr err;
+  dInt bs;
+
+  dFunctionBegin;
+  dValidHeader(fs,DM_CLASSID,1);
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
+  if (fn < 0 || bs <= fn) dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field number %d out of range",fn);
+  fs->fieldunit[fn] = unit;
+  dFunctionReturn(0);
+}
+dErr dFSGetFieldUnit(dFS fs,dInt fn,dUnit *unit)
+{
+  dErr err;
+  dInt bs;
+
+  dFunctionBegin;
+  dValidHeader(fs,DM_CLASSID,1);
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
+  if (fn < 0 || bs <= fn) dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field number %d out of range",fn);
+  *unit = fs->fieldunit[fn];
+  dFunctionReturn(0);
+}
 
 static dErr dFSBoundarySetClosure_Private(dFS fs,dMeshESH bset)
 {
@@ -333,6 +373,7 @@ dErr DMDestroy_dFS(DM dm)
   err = dFSGetBlockSize(fs,&bs);dCHK(err);
   for (dInt i=0; i<bs; i++) {err = dFree(fs->fieldname[i]);dCHK(err);}
   err = dFree(fs->fieldname);dCHK(err);
+  err = dFree(fs->fieldunit);dCHK(err);
   err = VecDestroy(&fs->gvec);dCHK(err);
   err = VecDestroy(&fs->dcache);dCHK(err);
   err = VecScatterDestroy(&fs->dscat);dCHK(err);

src/fs/interface/fsreg.c

@@ -61,6 +61,7 @@ dErr dFSCreate(MPI_Comm comm,dFS *infs)
 /* Defaults */
   fs->dm.bs = 1;
   err = dCallocA(1,&fs->fieldname);dCHK(err);
+  err = dCallocA(1,&fs->fieldunit);dCHK(err);
 
   *infs = fs;
   dFunctionReturn(0);

src/viewer/dhm.c

@@ -190,11 +190,13 @@ dErr dViewerDHMGetVecType(PetscViewer viewer,hid_t *intype)
 
   dFunctionBegin;
   if (dhm->h5t_vec < 0) {
-    hid_t vectype;
+    hid_t vectype,unitstype;
     herr_t herr;
+    err = dViewerDHMGetUnitsType(viewer,&unitstype);dCHK(err);
     vectype = H5Tcreate(H5T_COMPOUND,sizeof(dht_Vec));dH5CHK(vectype,H5Tcreate);
     herr = H5Tinsert(vectype,"fs",offsetof(dht_Vec,fs),H5T_STD_REF_DSETREG);dH5CHK(herr,H5Tinsert);
     herr = H5Tinsert(vectype,"time",offsetof(dht_Vec,time),dH5T_REAL);dH5CHK(herr,H5Tinsert);
+    herr = H5Tinsert(vectype,"units",offsetof(dht_Vec,units),unitstype);dH5CHK(herr,H5Tinsert);
     herr = H5Tinsert(vectype,"internal_state",offsetof(dht_Vec,internal_state),dH5T_INT);dH5CHK(herr,H5Tinsert);
     err = dViewerDHM_H5Tcommit(viewer,"Vec_type",vectype);dCHK(err);
     dhm->h5t_vec = vectype;