Avoid last remaining std::vectors in PolarGrid class

[julianlitz]

@EmilyBourne Let me know if the std::vectors like radii, angles, radial_spacings, angular_spacings need to be changed to our custom Vector data type. I was a bit irritated as you wanted to keep the std::array<int, 9> in the Stencil, but surely the PolarGrid std::vectors need to be replaced right?

[EmilyBourne]

Yes these definitely need replacing. std::vector is not GPU compatible. Whereas std::array is and gets you a stack array instead of a heap array. Everywhere std::vector is used we need a heap array anyway so there is no gain to be had from local stack memory

[julianlitz]

I could implement this on the 18.4 when I come back from vacation.

• Rename lengthSmootherRadial to lengthRadialSmoother
• Swap Vector to std::vector
• In anisotropicRadialRefinement still use std::vector but copy data to a Vector in the end.

@EmilyBourne

[EmilyBourne]

I could implement this on the 18.4 when I come back from vacation.

• Rename lengthSmootherRadial to lengthRadialSmoother
• Swap Vector to std::vector
• In anisotropicRadialRefinement still use std::vector but copy data to a Vector in the end.

@EmilyBourne

Fine for me, or @AbdelhadiKara may be able to do this before then (I am also on holiday until 18.4).

For anisotropicRadialRefinement do you mean constructRadialDivisions ? The size of the outputted Vector seems to be known before you fill it. It is only the allocation that needs doing locally for r_temp. Can't we simply do:

    AllocatableVector<double> r_temp;
    if (anisotropic_factor == 0) {
        // nr = 2**(nr_exp-1) + 1
        int nr                  = (1 << (nr_exp - 1)) + 1;
        double uniform_distance = (R - R0) / (nr - 1);
        assert(uniform_distance > 0.0);
        r_temp = AllocatableVector<double>("r_temp", nr);
        for (int i = 0; i < nr - 1; i++) {
            r_temp[i] = R0 + i * uniform_distance;
        }
        r_temp[nr - 1] = R;
    }

instead of copying?

[julianlitz]

@EmilyBourne Yes I think it will be pretty easy to switch from std::vector to Vector in PolarGrid.
There are just two hurdles:

1. In the Google Tests, there are still many instances of std::vector radii = {...}; that need to be updated. We have two options: either add a function that constructs our Kokkos vector from an initializer list, or use the element-wise assignment pattern we typically follow (e.g., radii(0) = a; radii(1) = b;).

2. The function you mentioned looks fine — the final vector size can be computed without issues. My main concern is PolarGrid::RadialAnisotropicDivision. This is the only function carried over unchanged from v1.0, and it relies heavily on std::vector-specific functionality. My recommendation here would be to leave it using std::vector internally and copy the result into a Kokkos vector at the end.

Note that this function, and more broadly, the entire PolarGrid constructor taking int nr_exp, int nr_theta, is only used in our custom mesh construction, so it's unlikely to be relevant for your purpose.

AI generated solution (only the last few lines are relevant)

#include "../../include/PolarGrid/polargrid.h"

void PolarGrid::RadialAnisotropicDivision(AllocatableVector<double>& r_temp, double R0, double R, const int nr_exp,
                                          double refinement_radius, const int anisotropic_factor) const
{
    // Calculate the percentage of refinement_radius.
    const double percentage = (refinement_radius - R0) / (R - R0);
    assert(percentage >= 0.0 && percentage <= 1.0);

    std::set<double, std::greater<double>>::iterator itr, itr_p1;

    /* uniform division of r in 2^nr_exp - 2^aniso */
    int dummy_lognr = nr_exp;
    int n_elems_equi = (1 << dummy_lognr) - (1 << anisotropic_factor);
    if (anisotropic_factor < 0 || n_elems_equi <= 0) {
        throw std::runtime_error("Please choose anisotropy factor a such that 2^fac_ani < 2^nr_exp.\n");
    }

    if ((anisotropic_factor % 2) == 1)
        n_elems_equi++;
    double uniform_distance      = (R - R0) / n_elems_equi;
    int nr                       = n_elems_equi + 1;
    std::vector<double> r_temp2(nr);
    for (int i = 0; i < nr - 1; i++)
        r_temp2[i] = R0 + i * uniform_distance;
    r_temp2[nr - 1] = R;

    /* refine around refinement_radius */
    int n_elems_refined = (1 << anisotropic_factor);

    int se;
    if (floor(nr * percentage) > nr - (n_elems_refined / 2)) {
        int new_aniso   = static_cast<int>(log2(nr - floor(nr * percentage)) + 1);
        n_elems_refined = (1 << new_aniso);
    }

    se     = static_cast<int>(floor(nr * percentage)) - n_elems_refined / 2;
    int ee = se + n_elems_refined;
    int st = static_cast<int>(ceil((double)n_elems_refined / 4.0 + 1) - 1);
    int et = static_cast<int>(floor(3 * ((double)n_elems_refined / 4.0)));

    std::set<double> r_set;
    std::set<double> r_set_p1;
    int count = 0;
    for (int i = 0; i < n_elems_refined; i++) {
        r_set_p1.insert(r_temp2[se + i]);
        count++;
    }
    double half = uniform_distance / 2.0;
    for (int k = 0; k < anisotropic_factor; k++) {
        std::set<double> r_set_p1_tmp;
        itr_p1     = r_set_p1.begin();
        int r_size = count;
        count      = 0;
        for (int i = 0; i < r_size - 1; i++) {
            r_set.insert((*itr_p1) + half);
            if (k < anisotropic_factor - 1 && i >= st && i < et) {
                r_set_p1_tmp.insert(*(itr_p1));
                r_set_p1_tmp.insert(*(itr_p1) + half);
                count += 2;
            }
            itr_p1++;
        }
        r_set_p1 = r_set_p1_tmp;
        half *= 0.5;
    }

    /* trim r_set so total size is 8*x+1 */
    nr        = nr + r_set.size();
    int shift = std::min(nr % 8 - 1, (int)r_set.size());
    itr       = r_set.begin();
    std::advance(itr, shift);
    r_set.erase(r_set.begin(), itr);
    for (int i = 0; i < n_elems_refined; i++)
        r_set.insert(r_temp2[se + i]);

    /* assemble into a local std::vector */
    int final_nr = n_elems_equi - n_elems_refined + (int)r_set.size() + 1;
    std::vector<double> r_result(final_nr);

    for (int i = 0; i < se; i++)
        r_result[i] = r_temp2[i];
    itr = r_set.begin();
    for (int i = 0; i < (int)r_set.size(); i++) {
        r_result[se + i] = *itr;
        ++itr;
    }
    for (int i = 0; i < n_elems_equi - ee + 1; i++)
        r_result[se + r_set.size() + i] = r_temp2[ee + i];

    /* copy into the Kokkos View — reallocate to the correct size first */
    r_temp = AllocatableVector<double>("r_temp", final_nr);
    for (int i = 0; i < final_nr; i++)
        r_temp(i) = r_result[i];
}

[julianlitz]

@EmilyBourne Third option with the least amount of work: Keep std::vector everywhere except the PolarGrid private members. That would be fine as well

[EmilyBourne]

I actually don't see anything in PolarGrid::RadialAnisotropicDivision stopping you using Vector either. The special functions seem to act on the sets rather than the vectors. But as you say there is not much to be gained from any change here given that the use of sets means that this would run on CPU anyway and given that it is only used for initialisation in tests on your side.

I think adding a constructor would make sense. The only question is which object would own the data.

std::vector<double> my_vec({1,2,3,4});
AllocatableVector<double> owning_vector(my_vec.size());
for (int i(0); i<my_vec.size(); ++i) {
    owning_vector[i] = my_vec[i];
}
AllocatableVector<double> non_owning_vector(my_vec.data(), my_vec.size());

[AbdelhadiKara]

Hello, i begun working on PolarGrid class. front() and back() calls were removed, currently I am trying to remove .begin() and .end() methods since it is not available for kokkos views. For example :

auto it = std::lower_bound(radii_.begin(), radii_.end(), splitting_radius.value()); should be replaced by

double* start=radii_.data(); double* end=radii_.data()+radii_.extent(0); auto it = std::lower_bound(start, end, splitting_radius.value());

If you see drawbacks let me know.

[EmilyBourne]

Hello, i begun working on PolarGrid class. front() and back() calls were removed, currently I am trying to remove .begin() and .end() methods since it is not available for kokkos views. For example :

auto it = std::lower_bound(radii_.begin(), radii_.end(), splitting_radius.value()); should be replaced by

double* start=radii_.data(); double* end=radii_.data()+radii_.extent(0); auto it = std::lower_bound(start, end, splitting_radius.value());

If you see drawbacks let me know.

If you include <Kokkos_StdAlgorithms.hpp> you can use begin and end instead of manipulating pointers:
https://kokkos.org/kokkos-core-wiki/API/algorithms/std-algorithms/Iterators.html

[julianlitz]

@AbdelhadiKara
Instead of

radii_(radii_.extent(0) - 1) + 1.0;

you can use

radii_(nr() - 1) + 1.0; or

radii_[nr() - 1] + 1.0;

[julianlitz]

@AbdelhadiKara Instead of

equals(angles(angles.extent(0) - 1), 2 * M_PI)

we can use

equals(angles(ntheta()), 2 * M_PI)

or

equals(angles[ntheta()], 2 * M_PI)

since angles is of size ntheta + 1.

[AbdelhadiKara]

@AbdelhadiKara Instead of

radii_(radii_.extent(0) - 1) + 1.0;

you can use

radii_(nr() - 1) + 1.0; or

radii_[nr() - 1] + 1.0;

@julianlitz You do not like extent method ?

[julianlitz]

Not a huge fan. @AbdelhadiKara
Something like std::ssize() or nr()-1 is more clear. Extent is a strange kokkos specific function which I haven't seen in any other library.

Also the extent(0) argument 0 seems strange since its just a simple one dimensional vector. @AbdelhadiKara

[EmilyBourne]

Not a huge fan. Something like std::ssize() or nr()-1 is more clear. Extent is a strange kokkos specific function which I haven't seen in any other library.

That is because other libraries don't have multi-D arrays. .extent(dim) is fairly normal in Kokkos. But for a 1D array I agree that just writing .size() would probably be more readable