add support for different data types in arrays, also on GPU #324

rhaas80 · 2024-12-06T14:34:44Z

This adds proper support for data types other than double to CarpetX, also for GPUs.

It introduces a helper class AnyVectorthat can hold any data type (similar to gdata in CarpetLib), as well as a number of switch statements to handle openPMD's typed output routines.

I am told that very new versions of openPMD-api lets one call the internal un-typed output routine, but at the time this code was written, that was not yet possible, and thus may not be possible on all cluster installed copies of openPMD.

Python code to dump data to compare:

#!/usr/bin/env python3

import openpmd_api as io
import sys
import numpy

series = io.Series(sys.argv[1], io.Access.read_only)

print("Read a Series with openPMD standard version %s" %
      series.openPMD)

print("The Series contains {0} iterations:".format(len(series.iterations)))
for i in series.iterations:
    print("\t {0}".format(i))
print("")

it = 0
i = series.iterations[it]
print("Iteration {0} contains {1} meshes:".format(it, len(i.meshes)))
for m in i.meshes:
    print("\t {0}".format(m))
print("")

gf = i.meshes["testoutput_gf_patch00_lev00"]["testoutput_gf"]
shape = gf.shape
print("Field gf has shape {0} and datatype {1}".format(
      shape, gf.dtype))
arrdata = gf.load_chunk()
series.flush()
numpy.savetxt("testoutput_gf.asc", arrdata.flatten())

for var in ["sc", "a1", "a2", "a3"]:
    for suffix in [""]: #, "_int", "_complex"]:
        data = i.meshes["testoutput_"+var+suffix]["testoutput_"+var+suffix]
        print("Field {0} has shape {1} and datatype {2}".format("testoutput_"+var+suffix, data.shape, data.dtype))
        arrdata = data.load_chunk()
        series.flush()
        numpy.savetxt("testoutput_"+var+suffix+".asc", arrdata.flatten())

rhaas80 · 2024-12-06T14:36:30Z

Ticket is here: https://bitbucket.org/einsteintoolkit/tickets/issues/2836/add-support-for-different-data-types-to

rhaas80 · 2024-12-06T16:34:23Z

Fixing the various CUDA issues right now.

CarpetX/src/io_openpmd.cxx

eschnett · 2024-12-06T23:46:44Z

CarpetX/src/driver.hxx

+      AnyTypeVector(const AnyTypeVector &) = delete;
+      AnyTypeVector &operator=(const AnyTypeVector &) = delete;
+      AnyTypeVector &operator=(const AnyTypeVector &&other) {
+        this->_type = other._type;


This is incorrect. The moved-from object will still have its destructor called, which would deallocate the memory. The moved-from object needs to be put into an "inactive" state.

Fixed in 617d0d9

eschnett · 2024-12-06T23:46:53Z

CarpetX/src/driver.hxx

+        return *this;
+      }
+      AnyTypeVector(AnyTypeVector &&other)
+          : _type(other._type), _data(other._data), _count(other._count) {}


Fixed in 617d0d9

eschnett · 2024-12-06T23:48:57Z

CarpetX/src/driver.hxx

+      AnyTypeVector() : _type(-1), _data(nullptr), _count(0){};
+      AnyTypeVector(int type_, size_t count_)
+          : _type(type_),
+            _data(amrex::The_Arena()->alloc(CCTK_VarTypeSize(type_) * count_)),


This looks dangerous. Passing an invalid type makes CCTK_VarTypeSize return a negative number, and C++ will convert this into a very large unsigned number for alloc. Better to check and call CCTK_ERROR.

Moved into constructor code and / or added an assert for the correct values. Those who compile with NDEBUG are on their own. In 4fc6ceb

eschnett · 2024-12-06T23:49:55Z

CarpetX/src/driver.hxx

+
+      int type() const { return _type; };
+
+      void *data_at(size_t i) const {


There should probably be two versions of this function, one const and returning a const void*, and the other non-const and returning a void *.

Done in 3605564. Note that differs from AMReX's Array4D and `GridPtrDesc1` class (and thus requires a const-cast somehwere in schedule.cxx)

For grid-functions:
```
auto &restrict groupdata = *leveldata.groupdata.at(gi);
const GridPtrDesc1 grid1(leveldata, groupdata, mfp);
for (int tl = 0; tl < int(groupdata.mfab.size()); ++tl) {
const amrex::Array4<CCTK_REAL> vars =
groupdata.mfab.at(tl)->array(mfp.index());
for (int vi = 0; vi < groupdata.numvars; ++vi)
cctkGH->data[groupdata.firstvarindex + vi][tl] = grid1.ptr(vars, vi);
}
```

For grid-arrays (and scalars):
```
auto &restrict arraygroupdata = *globaldata.arraygroupdata.at(gi);
for (int tl = 0; tl < int(arraygroupdata.data.size()); ++tl) {
const auto &restrict vars = arraygroupdata.data.at(tl);
for (int vi = 0; vi < arraygroupdata.numvars; ++vi)
cctkGH->data[arraygroupdata.firstvarindex + vi][tl] =
const_cast<void*>(vars.data_at(vi * arraygroupdata.array_size));
```

eschnett · 2024-12-06T23:51:19Z

CarpetX/src/io_openpmd.cxx

+                  static_cast<CCTK_COMPLEX *>(cactus_var_ptr), start, count);
+              break;
+            default:
+              assert(0 && "Unexpected variable type");


I'd call CCTK_ERROR here because it can't be disabled by compiler options.

this is not user triggerable (via parfile input) and the list is fixed by Cactus (not user or CarpetX code), so I'd rather keep the assert, if possible.

eschnett · 2024-12-06T23:54:55Z

CarpetX/src/io_openpmd.cxx

+                std::vector<char> poison(typesize);
+                // for int: deadbeef for little endian machines
+                for (size_t i = 0; i < typesize; i += sizeof poison)
+                  std::memcpy(&poison[i], &ipoison, std::min(poison.size(), typesize - i));


I would define three different poison values here (for int, real, and complex). Using just some of the real bits for int might give values that aren't very distinctive. I found that large positive values make for good integer poison. Negative values don't work that well because they are more likely to be clipped off.

done in 49de23b

eschnett · 2024-12-06T23:55:41Z

CarpetX/src/io_openpmd.cxx

+                  start, count);
+              break;
+            default:
+              assert(0 && "Unexpected variable type");


I would call CCTK_ERROR here.

#324 (comment)

eschnett · 2024-12-06T23:57:28Z

CarpetX/src/io_openpmd.cxx

+          case CCTK_VARIABLE_COMPLEX:
+            return openPMD::determineDatatype<CCTK_COMPLEX>();
+          default:
+            assert(0 && "Unexpected varType");


CCTK_ERROR

#324 (comment)

eschnett · 2024-12-06T23:58:16Z

CarpetX/src/io_openpmd.cxx

+            switch (cgroup.vartype) {
+            case CCTK_VARIABLE_REAL:
+              record_components.at(vi).storeChunkRaw(
+                  static_cast<CCTK_REAL const *const>(var_ptr), start, count);


This should be static_cast<CCTK_REAL const *>, the additional const doesn't have any effect since the pointer is an rvalue.

Done in 016fec7

eschnett · 2024-12-06T23:59:56Z

CarpetX/src/io_tsv.cxx

-    file << sep << arraygroupdata.data.at(tl).at(vi);
+    switch (cgroup.vartype) {
+    case CCTK_VARIABLE_REAL:
+      file << sep << *(CCTK_REAL *)arraygroupdata.data.at(tl).data_at(vi);


Should be (const CCTK_REAL *)

Done in 4704af8

eschnett · 2024-12-07T00:00:36Z

CarpetX/src/driver.cxx

+  for (size_t i = 0; i < anytypevector.size(); ++i) {
+    switch (anytypevector.type()) {
+    case CCTK_VARIABLE_COMPLEX:
+      yaml << *(CCTK_COMPLEX *)anytypevector.data_at(i);


This should be (const CCTK_COMPLEX *).

Done in 8f8fbad

eschnett · 2024-12-07T00:01:21Z

CarpetX/src/io_tsv.cxx

+      file << sep << col++ << ":" << varname << ".real";
+      file << sep << col++ << ":" << varname << ".imag";
+    } else
+      assert(0 && "Unexpected variable type");


Better call CCTK_ERROR

Since this is not user (parfile or other runtime input) triggerable and the list of allowed options is determined by the flesh and not the thorn being written (CarpetX) I'd rather keep the assert (since it's a bug) rather than the CCTK_Error.

eschnett · 2024-12-07T00:01:56Z

CarpetX/src/io_tsv.cxx

+      switch (cgroup.vartype) {
+      case CCTK_VARIABLE_REAL:
+        file << sep
+             << *(CCTK_REAL *)arraygroupdata.data.at(tl).data_at(


(const CCTK_REAL *)

Done in 4704af8

eschnett · 2024-12-07T00:02:32Z

CarpetX/src/io_tsv.cxx

+        break;
+      case CCTK_VARIABLE_INT:
+        file << sep
+             << *(CCTK_INT *)arraygroupdata.data.at(tl).data_at(


The expression arraygroupdata.data.at(tl).data_at(np * vi + DI[out_dir] * i) could be calculated before the switch statement to reduce code duplication.

It could, something like this:
```
for (int vi = 0; vi < arraygroupdata.numvars; ++vi)
file << sep;
const void *ptr = arraygroupdata.data.at(tl).data_at(np * vi + DI[out_dir] * i);
switch (cgroup.vartype) {
case CCTK_VARIABLE_REAL:
file << *(const CCTK_REAL *)ptr;
break;
case CCTK_VARIABLE_INT:
file << *(const CCTK_INT *)ptr;
break;
case CCTK_VARIABLE_COMPLEX: {
CCTK_COMPLEX value = *(const CCTK_COMPLEX *)ptr;
file << value.real() << sep << value.imag();
} break;
default:
assert(0 && "Unexpected variable type");
break;
}
file << "\n";
}
```
though it comes at the price of information about what is output being less local and requiring backtracking in the code when reading it.

In addition to AnyVector, you could introduce an AnyScalar, and add a virtual function that outputs to a stream. (Not sure whether that makes the code simpler or more complex in the end.)

Hmm, maybe. Not having all those switch scattered around would be nice.

The tricky one is the CCTK_COMPLEX output since it currently uses sep which would not be known to AnyScalar::operator<< .

Still an AnyScalar and a AnyScalar AnyVector::operator[](size_t i) const might be useful. Would be a return by value (at least unless I want to really complicate my life and have the AnyScalar carry a back-reference to its parent AnyVector so that assignment to it would change the AnyVector element) for now.

Something like 98ab37c ? So far only used for output. Could possibly also be used for setting poison values (though there performance issues may show up).

eschnett · 2024-12-07T00:04:02Z

CarpetX/src/valid.cxx

-  CCTK_REAL poison;
-  std::memcpy(&poison, &ipoison, sizeof poison);
+  const size_t typesize = size_t(CCTK_VarTypeSize(group.vartype));
+  char poison[typesize];


Didn't you switch to std::vector elsewhere for a similar construct?

Yes, and here it does not complain about the varuable sized array, but it does in io_openpmd.cxx. I have no clue why.

I will redo this all anyway, to have one poison code for both places where it is used.

eschnett · 2024-12-07T00:04:38Z

CarpetX/src/valid.cxx

-  CCTK_REAL poison;
-  std::memcpy(&poison, &ipoison, sizeof poison);
+  const size_t typesize = size_t(CCTK_VarTypeSize(group.vartype));
+  char poison[typesize];


std::vector?

removed in 1834168

eschnett · 2024-12-07T00:05:53Z

CarpetX/src/valid.cxx

+        using std::isnan;
+        if (CCTK_BUILTIN_EXPECT(nan_handling == nan_handling_t::allow_nans
+                                    ? std::memcmp(ptr, poison, sizeof *ptr) == 0
+                                    : isnan(*ptr),


You could call isnan(real) || isnan(imag) for complex numbers.

I guess, yes.

done in 49de23b

eschnett

I've left a few comments. Some are important and need to be addressed. Most are suggestions.

used clang-format-16 -style file -i

rhaas80 requested review from shankar-1729 and eschnett December 6, 2024 14:36

rhaas80 force-pushed the rhaas/anyvector branch 3 times, most recently from 8e0e23b to a665253 Compare December 6, 2024 15:10

eschnett reviewed Dec 6, 2024

View reviewed changes

CarpetX/src/io_openpmd.cxx Outdated Show resolved Hide resolved

rhaas80 force-pushed the rhaas/anyvector branch from 0a2d057 to a722e7d Compare December 6, 2024 17:19

rhaas80 requested a review from eschnett December 6, 2024 21:59

eschnett reviewed Dec 6, 2024

View reviewed changes

eschnett reviewed Dec 7, 2024

View reviewed changes

eschnett requested changes Dec 7, 2024

View reviewed changes

rhaas80 added 2 commits January 1, 2025 14:42

CarpetX: add support for any type of grid arrays and scalars

152f739

TestArrayGroup: add tests for integer and complex grid arrays

73c6d3d

rhaas80 added 17 commits January 1, 2025 14:42

TestOutput: add interger and complex grid array and grid scalar

c4de77c

CarpetX: avoid varsize array

dfa47b2

CarpetX: add missing debug output for COMPLEX

24adbe1

CarpetX: correct move operators for AnyVector

7315838

CarpetX: keep track of type size internally

10bed84

CarpetX: correct count in AnyTypeVector::alloc()

f73cdf3

CarpetX: re-use alloc() function in allocating constructor

43a69e5

CarpetX: split AnyTypeVector::data into const and non-const

06f0c23

CarpetX: check array access in AnyTypeVector

96c4c50

CarpetX: consolidate poison values into a single class

0b401c5

CarpetX: use common poison code in openPMD output

5f41762

CarpetX: remove redundant const in cast

286e9e4

CarpetX: add const in cast for output

6b7d888

CarpetX: add const in cast for debug output

00264ef

CarpetX: clang-format code

dc86fe5

used clang-format-16 -style file -i

CarpetX: introduce AnyTypeScalarRef

9dc3680

CarpetX: output invalid array indices in debug mode

b5f228f

rhaas80 force-pushed the rhaas/anyvector branch from 3ab11f1 to b5f228f Compare January 1, 2025 20:42


		int type() const { return _type; };

		void *data_at(size_t i) const {

add support for different data types in arrays, also on GPU #324

Are you sure you want to change the base?

add support for different data types in arrays, also on GPU #324

Conversation

rhaas80 commented Dec 6, 2024

rhaas80 commented Dec 6, 2024

rhaas80 commented Dec 6, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rhaas80 Dec 7, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

eschnett left a comment

Choose a reason for hiding this comment

rhaas80 Dec 7, 2024 •

edited

Loading