Add double precision option for HDF5 output about athena

[felker]

This is still true for most data, correct? What would need to be done to extend output to double precision? I could use this capability for some linear wave tests that I am currently running.

[c-white]

For a quick hack, just changing H5T_IEEE_F32BE to H5T_IEEE_F64BE should do the trick, especially in the line beginning datasets_celldata[n]. (The Python reader actually creates double precision arrays in any case, which it probably shouldn't do.) I'll work on smoothly doing this with a configure option. For small but precise tests, another workaround is to use the tab format with format specifier %24.16e.

[felker]

• The dtype option of athdf() is defunct; should it be removed from the wiki, or will it be added soon?
• After changing H5T_IEEE_F32BE to H5T_IEEE_F64BE in the H5Dcreate() calls, the resulting dataset type still appears as a 4-byte float in viewers such as HDFCompass. Should H5T_NATIVE_FLOAT also be replaced with H5T_NATIVE_DOUBLE in the H5Dwrite() calls, @c-white?

[felker]

Also, the suggested quick hack for outputs/athena_hdf5.cpp doesn't seem to produce the correct values to double precision when imported into Python with the athdf() function.

E.g., the Orszag-Tang test's initial uniform density profile only agrees to 8 digits with the values I am observing internally in the code and from tab format output generated with a 16 digit format specifier.

[c-white]

• dtype is working in the development branch. Are you using the public version?

• Changing the format in H5Dcreate() calls changes how the file is written, no matter what else is happening. I'm not familiar with HDFCompass, but in Python you can check

  import h5py
  with h5py.File('output.athdf', 'r') as f:
    print(f['prim'].dtype)

  This should give >f4 or >f8 for single or double precision.

• You are right that the H5Dwrite() calls need to be updated, but not because they specify how the data is to look in the written file, confusingly. That second argument specifies the format of the data in memory (pointed to by the last argument). In particular, one must also change the declaration of data_buffers (and the corresponding coordinate buffers if one desires) from float to double . The unaltered code copies the doubles into float arrays (the copy is done anyway to make the active zone contiguous in memory), and it is these that are written to file. This explains why the values were only agreeing to single precision.

[felker]

The original issue here ended up being due to yt's anti-aliasing causing graphical artifacts. But I believe that the discrepancy between HDF5 files and internal code representation was still there. Also, I recommend against using HDFCompass.

I am going to close this, but I believe we still should update the H5Dwrite() calls, so I will mark it wontfix.

[felker]

Ok, I actually need this functionality right now, so I reopened the issue.

Do we want a user switch for switching between:

• H5T_NATIVE_FLOAT --> H5T_NATIVE_DOUBLE
• H5T_IEEE_F32BE --> H5T_IEEE_F64BE
• data_buffers allocation to doubles

and how should it be implemented? Runtime or configure-time option?

I thought about using the existing typedef float Real; functionality, but I assume users may want to toggle the output precision independently of the solver's precision. Although it wouldnt make sense to write out double precision values if the solver is single precision.... should those settings be mutually exclusive?

I think I will go ahead with duplicating that configure-time functionality with a new H5Real type for now.