0

I'm writing a numerical fluid solver in C++ as a hobby project. I will try to explain what I want to accomplish in a simplified manner.

The solver has multiple flow variables (density, velocity, pressure, etc.) stored in each cell in a grid. I would want a convenient way to access the variables and doing computations on them (typically with operator overloading). They are now stored as double* array of size N, where each flow variable belonging to the same cell are stored consecutively as this: density0, u0, v0, w0, pressure0, density1, u1, v1, w1, pressure1 ... density_N-1, u_N-1, v_N-1, w_N-1, pressure_N-1

Keep in mind that I would like to keep everything general; in this specific case there were 5 flow variables, but there might also be a different amount.

What I would ideally like is to have a way to reinterpret my flow variables as a single cell variable without having to copy the memory. In this case the variable in a cell could for instance be a struct like this:

    struct FlowVar{
        double density, u, v, w, p;
    };

I know that there is something called "type-punning" which would allow you to reinterpret memory as a different type. This little example illustrates how the flow variable in cell 10 could be accessed this way:

    double* raw_data = new double[100]; 

    for (int i{0};i<100;i++) raw_data[i] = i;

    FlowVar* flow_var_10 = (FlowVar*)&raw_data[9];

Even though I got the correct variables when running this (9,10,11,12,13) , this is apparently undefined behaviour in C++ https://adriann.github.io/undefined_behavior.html

I have heard about something called std::bit_cast, but my impression is that is can't be used for my kind of purpose. However, please inform me if I'm wrong here.

So at this point I had no defined way to accomplish what I wanted. The next possible solution I checked out was to use the linear algebra library Eigen. I would then use a Eigen::Vector<double, 5> to represent a flow variable. Using Eigen is also convenient in its own right, since it has lots of useful linalg functionality. However, I am not really sure if Eigen is slower or faster than homemade matrix/vector classes for small sizes, so it might be a bad decision Is Eigen slow at multiplying small matrices? .

Eigen has a functionality called Map which allows mapping raw data to vector or matrix types without copying. I'm not sure how this is achieved in a defined and safe way, but I guess it is beyond the level of the average C++ enthusiast.

To map the raw data to a flow variable I could now do something like this:

    using Vec5 = Eigen::Vector<double,5>;
    using FlowVar = Eigen::Map<Vec5>;
    
    double* raw_data = new double[100];

    for (int i{0};i<100;i++) raw_data[i] = i;

    FlowVar flow_var = FlowVar(raw_data + 9);

Now FlowVar shares some of the memory with raw_data, in effect accomplishing the same purpose as the above type punning.

However I fear that this solution might be inefficient as I'm using small vectors and have many grid points and will need to create Maps often. The size of a Eigen::Map (at least on my computer) is 16 bytes, which is more than for instance references and pointers.

I would like some opinions on what design decision would likely be the best here. Where I stand now I have four options:

1: Use the undefined type punning - which seems to work fine for doubles in my case...

2: Use the Eigen::Map solution

3: Simply copy the data to a struct or Eigen::Vector when wanting or needing to view the raw_data as a FlowVar

4: Simply drop the entire FlowVar type and only access the raw_data directly

I would be grateful for some opinions here. Should I pick one of my four options, or are there other possibilities that I'm not aware of?

ander
  • 11
  • 3
  • Why don't you just create an array of `FlowVar` directly? You could simply fill it like `FlowVar data[64]; size_t n = 0; for(auto& d : data) { d.density = n++; d.u = n++; d. [...] }` – Aconcagua Jun 10 '23 at 10:55
  • You could as well provide a constructor for your FlowVar type. – Aconcagua Jun 10 '23 at 10:58
  • @Aconcagua I actually started with this kind of design, I created something like: template struct FlowVars{ double variables [N_VARS]; };to be able to solve different kinds of equations. However, I found this solution inconvenient for various reasons. I found it easier to design the class structure if the FlowVar type of the raw_data don't have to be specified at compile time. – ander Jun 10 '23 at 11:20
  • How about a "view": `struct FlowVarView{ double* data; double& density() { return data[0]; } /* const version and similar for u, v, w, p */ };`? – Jarod42 Jun 10 '23 at 11:29
  • `Eigen::Map` simply holds a pointer to the referenced external data. Basically, `Eigen::Map` is to `Eigen::Matrix` as `std::string_view` is to `std::string`. I don't know where the 16 byte size comes from. For a fixed size vector, it only contains the pointer and this is also the only member that is initialized when you look at the generated assembly code. Overall, `Eigen::Map` can be created very cheaply on-the-fly, basically with zero cost. However, its use should be temporary. Don't store arrays of Maps. – Homer512 Jun 10 '23 at 12:05
  • You could consider `Map>` for an array of variables. However, for your overall memory layout, you currently use an array-of-structs (AoS). A struct of array (SoA) [may have better performance](https://stackoverflow.com/questions/40163722/is-my-understanding-of-aos-vs-soa-advantages-disadvantages-correct) depending on your access pattern. This could be done simply by transposing into `Matrix` (note that Eigen's matrices are [column-major by default](https://eigen.tuxfamily.org/dox/group__TopicStorageOrders.html)) – Homer512 Jun 10 '23 at 12:08
  • @Homer512 , At least when I evaluate for instance `sizeof(Eigen::Matrix>)` I get 16, but it would indeed make more sense that it has the same size as a reference. – ander Jun 10 '23 at 14:23
  • I am aware of the distinction between AoS and SoA. I think it is reasonable to use AoS in my case, since the equations in each cell are typically evaluated together for the type of coupled numerical schemes I am using. At this stage I have just made two custom classes for holding the solution variables, which are not Eigen based. One of these store for instance 5x3 matrices in each cell, representing the gradient of the solution. Not sure if I could create such a construct with eigen directly as it would be a vector of matrices with dimensions `N_CELLS x N_EQUATIONS x 3`. – ander Jun 10 '23 at 14:23
  • I have used the `Eigen::Map` on the fly as you describe for now, hopefully this will turn out to be a decent solution. I added a member function in my solution data class that can return the flow variable in cell i as an eigen map – ander Jun 10 '23 at 14:29
  • 1
    Flexible design and uniquely identifiable variable names contradict each other somehow... If you want to have more variables (dimensions?), how would you want to get additional names for these??? – Aconcagua Jun 10 '23 at 21:18
  • @Aconcagua The way I designed it is that the class holding the solution field is flexible (number of rows are determined at compile time), but when working with a specific equation (for instance Navier Stokes / Euler Equations as described above, the number of rows are known. – ander Jun 13 '23 at 11:23
  • Pretty clear that for a specific problem number of rows are known... But if you design your class for 3 rows and next problem needs four rows, where would you get a new variable or getter name from? Well, you could have a template, provide some *maximum* set of getters and disable the surplus ones by SFINAE. But then again you'd need to extend that template each time you discover yet another row is needed. An index based solution remains more flexible. – Aconcagua Jun 13 '23 at 11:57
  • @Aconcagua I think I perhaps focused on the wrong things when explaining what I wanted to achieve. The important thing for me was not to be able to get for instance have getters for density, velocity etc, the point was to be able to look at the field data as a vector like container with operators defined and range checking in debug mode, etc, instead of just accessing the variables in a cell as a raw pointer, which I guess would be less safe and inconvenient – ander Jun 13 '23 at 12:13
  • I have a class Solver that is abstract, this class contains a SolverData class which again contains all the different fields. I found it very inconvenient to design it this way if I had to specify the number of equations for this generic SolverData class at compile time. A specific Solver type, for instance EulerSolver is derived from Solver, and now the different fields will be constructed with the correct number of rows. When writing the specific methods used by EulerSolver, the number of rows are known (5 rows) and then it could be convenient to access data in each cell as static 5d vectors – ander Jun 13 '23 at 12:15
  • Then why not provide your `FlowVar` with an index operator? Are you relying on a common type for all flow vars or can they be separate types for separate sizes? In the latter case `template struct FlowVar { std::array m_data; };` could be an interesting alternative. Otherwise you'd continue with an array (though you might use a `std::unique_ptr` to avoid all the hassle with memory management). – Aconcagua Jun 13 '23 at 13:29
  • If you additionally provide `begin` and `end` returning some iterators (in most simple case just pointers to the beginning and one past the end of the internal array…) then you could even use range based for loops, that could look like `std::vectordata(whatEverLength, FlowVar(whateverDimensions)); double n = 0.0; for(auto& d : data) { for(auto& v : d) { v = n; n += 1.0; } }` then. – Aconcagua Jun 13 '23 at 13:33

1 Answers1

0

To continue some aspects from the comments in a larger text:

At least when I evaluate for instance sizeof(Eigen::Matrix<Eigen::Vector<double,5>>) I get 16, but it would indeed make more sense that it has the same size as a reference.

Pretty sure you meant to write Eigen::Map<Eigen::Vector<...>>. But yes, that type functionally only needs 8 byte for a single pointer. I have yet to read enough of the code to understand where the second member comes from. If you change the Map to something that needs a second member to store the runtime size, e.g. Map<Vector<double, Eigen::Dynamic>> its size becomes 24. But whatever, if you just use Maps as local variables or the occasional class member, it doesn't really matter.

Not sure if I could create such a construct with eigen directly as it would be a vector of matrices with dimensions N_CELLS x N_EQUATIONS x 3

That is indeed a limitation of Eigen. If you want to look for alternatives that can handle multidimensional data, the keyword is "tensor". There is an unsupported tensor extension for Eigen. Other libraries might include PyTorch's C++ frontend but I cannot vouch for the quality of either library. Personally I just flatten the outer dimensions into extra columns or rows, as appropriate, but that's not really clean.

I added a member function in my solution data class that can return the flow variable in cell i as an eigen map

You mean something like this?

using Vector5d = Eigen::Vector<double, 5>;
struct FlowVar{
    double density, u, v, w, p;

    Eigen::Map<Vector5d> as_vector() noexcept
    { return Eigen::Map<Vector5d>(&u); }
};

I'm not sure this doesn't violate strict aliasing rules. However, I'm not a language lawyer and could be wrong about this. I suggest this alternative that definitely doesn't have that issue:

struct FlowVar{
    Vector5d as_vector;

    double& density() noexcept
    { return as_vector.x(); }

    double& u() noexcept
    { return as_vector.y(); }

    double& v() noexcept
    { return as_vector.z(); }

    double& w() noexcept
    { return as_vector.w(); }

    double& p() noexcept
    { return as_vector[4]; }
};

Similar to the accessors used by std::complex. With an optimizing compiler this should have zero overhead.

Homer512
  • 9,144
  • 2
  • 8
  • 25
  • 1
    Re: "I have yet to read enough of the code to understand where the second member comes from." There are two members for the dimension, which are empty classes if the size is fixed, but the objects still need 1 byte each, so with padding this requires 8 (pointer) + 2 (sizes) + 6 (padding) = 16 bytes – chtz Jun 13 '23 at 09:38
  • @chtz If two sizes need just two bytes sizes would be limited maximally 255 – wouldn't it be more meaningful to have four byte wide integers instead? – Aconcagua Jun 13 '23 at 12:02
  • @Aconcagua If the size is known at compile time, it should actually need no space at all. If it is unknown, then 8 bytes are used (actually the size of `ptrdiff_t` so it depends on what system you are compiling). – chtz Jun 13 '23 at 12:05
  • @chtz That's all clear – solely: 2 bytes only for the sizes and 6 for padding doesn't appear meaningful to me, if we need that space anyway then we'd rather select types that are not limited to such small ranges while at the same time wasting bytes for padding, wouldn't we? – Aconcagua Jun 13 '23 at 13:22
  • @Aconcagua `Eigen::Map` is not intended for long-time storage, but should usually only be a temporary. And with proper inlining the unused two bytes will be optimized away anyways. Otherwise, it is even worse, if one size is dynamic and the other is fixed: 8 byte for the pointer, 8 byte for the dynamic size, 1 byte for the fixed size (never accessed) and 7 bytes padding. – chtz Jun 13 '23 at 14:39
  • @chtz And if both sizes are dynamic? There's a constructor for. Template parameters do not comprise sizes, so we need to be able to store both dimensions then. Why then again fall into the other extremum, 4 + 4 bytes would be much more reasonable (and would result in 12 total size for 32-bit pointers). – Aconcagua Jun 13 '23 at 17:29
  • Pretty unlucky: Inheriting from some type `MapBase` which doesn't seem to be defined? Is [this](https://eigen.tuxfamily.org/dox/Map_8h_source.html) correct documentation at all? – Aconcagua Jun 13 '23 at 17:34
  • @Aconcagua regarding it being 1 byte each: That's not a design choice on Eigen's part. The dimension is stored in an object. Since the size is compile-time, it becomes an empty class. But objects of empty classes still need distinct pointers in C++, so the compiler makes them 1 byte. And then the compiler adds padding because the pointer to the data array needs 8 byte alignment. Therefore `sizeof(Map)` needs to be a multiple of 8 – Homer512 Jun 13 '23 at 18:22
  • @chtz thanks for clarifying. I was suspecting a missing empty base specialization for the same reasoning. It's interesting that `Vector` avoids the extra member but `Map` doesn't. But yeah, for the intended use, it doesn't matter too much – Homer512 Jun 13 '23 at 18:26
  • @Aconcagua Relevant reading: http://www.cantrip.org/emptyopt.html Basically, Eigen doesn't bother with this trick in this particular instance – Homer512 Jun 13 '23 at 18:35
  • @Aconcagua Properly documenting Eigen's template magic is hard with doxygen. Look [here](https://eigen.tuxfamily.org/dox/MapBase_8h_source.html), [here](https://eigen.tuxfamily.org/dox/classEigen_1_1MapBase_3_01Derived_00_01ReadOnlyAccessors_01_4.html), and [here](https://eigen.tuxfamily.org/dox/classEigen_1_1MapBase_3_01Derived_00_01WriteAccessors_01_4.html) for some details. The "magic" with fixed vs dynamic happens in the implementation of `internal::variable_if_dynamic` (which I guess is not in the doxygen docu at all, since it is "internal"). – chtz Jun 14 '23 at 01:39