[newchem-cpp] Introduce Multi1dView

[mabruzzo]

Overview

This introduces the Multi1DView<T> class template. This construct internally tracks an array of pointers and it provides 2 ways to access the underlying data.

Essentially, this type makes an array of pointers act more like a View<T**> (our type for representing a 2D array).

Motivation

The motivation is to use this as an alternative to grackle_field_data* for accessing field data. Let's sketch out an example. Suppose that we want have a function that needs to access the proton density and neutral Hydrogen density at each index in a provided IndexRange. We'll sketch out out what this looks like right now, and the 2 ways to do it with Multi1DView. For the purpose of this discussion, let's assume that we've taken advantage of the ability to store an internal pointer offset to account for the j,k offsets in an idx_range.

1. Variant with grackle_field_data*

The following logic is representative of the coding idioms currently in use throughout the code base

void fn(IndexRange idx_range, grackle_field_data* field_data, /* <args...> */) {
  View<const gr_float***> HI_density(
      field_data->HII_density, field_data->grid_dimension[0],
      field_data->grid_dimension[1], field_data->grid_dimension[2]);
  View<const gr_float***> HII_density(
      field_data->HII_density, field_data->grid_dimension[0],
      field_data->grid_dimension[1], field_data->grid_dimension[2]);

  for(int i = idx_range.i_start; i < idx_range.i_stop; i++) {
    gr_float rho_proton = HII_density(i, idx_range.j, idx_range.k);
    gr_float rho_neutral_H = HI_density(i, idx_range.j, idx_range.k);
    // < do calculation with rho_proton & rho_neutral_H ...>
  }
}

2. Variant usingMulti1DView::contig1d_ptr method

The following logic should have comparable (maybe marginally better performance than variant 1). The idea would to transition to this style of code throughout the codebase.

void fn(IndexRange idx_range, Multi1DView<gr_float> rho_sp, /*<args...>*/) {
  const gr_float* HI_density = rho_sp.contig1d_ptr(SpLUT::HI);
  const gr_float* HII_density = rho_sp.contig1d_ptr(SpLUT::HII);

  for(int i = idx_range.i_start; i < idx_range.i_stop; i++) {
    gr_float rho_proton = HII_density[i];
    gr_float rho_neutral_H = HI_density[i];
    // < do calculation with rho_proton & rho_neutral_H ...>
  }
}

Note: I'm open to suggestions for better names for this method.

3. Variant using mulit-dimensional indexing

The performance of this approach will generally be comparable or slightly worse than variant 2. I sketch out the reasons for why this is the case in in the docstring for Multi1DView::contig1d_ptr (but I'm happy to provide more details).

void fn(IndexRange idx_range, Multi1DView<gr_float> rho_sp, /*<args...>*/) {
  for(int i = idx_range.i_start; i < idx_range.i_stop; i++) {
    gr_float rho_proton = rho_sp(i, SpLUT::HI);
    gr_float rho_neutral_H = rho_sp(i, SpLUT::HII);
    // < do calculation with rho_proton & rho_neutral_H ...>
  }
}

Note

The fact the order of indices in rho_sp(i, SpLUT::HI), rather than rho_sp(SpLUT::HI, i) is a consequence of the fact that we our regular View object inherit the Fortran multi-dimensional array layouts in the . If we want to change this, now would probably be the time to do so (I may open an issue to discuss this if I have time).

Reviewing Considerations

I suspect that you may want to wait to merge this until you see the followup PR that actually starts using this machinery. We currently don't use this new machinery outside of a few new unit-tests.

The actual code changes in this PR are quite minimal -- a lot of this is docstrings explaining how things work.