Skip to content

Custom NVector for GPU #3390

Merged
bendudson merged 36 commits into
nextfrom
custom-nvector
Jun 18, 2026
Merged

Custom NVector for GPU #3390
bendudson merged 36 commits into
nextfrom
custom-nvector

Conversation

@bendudson

@bendudson bendudson commented Jun 16, 2026

Copy link
Copy Markdown
Contributor

Intended to allow SUNDIALS to operate on BOUT++ state data without copying into vectors. The hope is to avoid copying between CPU and GPU.

CVODE, IDA and Arkode can now use this ManyVector backend by setting

solver:nvector=manyvector

The default is still to use the MPI parallel N_Vector packaged with Sundials.

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 42. Check the log or trigger a new build to see more.

Comment thread src/field/field2d.cxx
}
}

void Field2D::swapData(Field2D& other) { std::swap(data, other.data); }

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "std::swap" is directly included [misc-include-cleaner]

void Field2D::swapData(Field2D& other) { std::swap(data, other.data); }
                                              ^

Comment thread src/solver/impls/arkode/arkode.hxx Outdated
N_Vector nvector_from_state(const sundials::Context& ctx) {
std::vector<N_Vector> subvectors;
subvectors.reserve(f2d.size() + f3d.size());
const auto inserter = std::back_inserter(subvectors);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "std::back_inserter" is directly included [misc-include-cleaner]

src/solver/impls/arkode/arkode.hxx:33:

+ #include <iterator>

Comment thread src/solver/impls/arkode/arkode.hxx Outdated
const auto var_str_to_nvector = [&ctx](auto &var_str) {
return BoutNVector::create(ctx, *var_str.var, var_str.evolve_bndry);
};
std::transform(f2d.cbegin(), f2d.cend(), inserter, var_str_to_nvector);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "std::transform" is directly included [misc-include-cleaner]

src/solver/impls/arkode/arkode.hxx:33:

+ #include <algorithm>

Comment thread src/solver/impls/arkode/arkode.hxx Outdated
return BoutNVector::create(ctx, *var_str.var, var_str.evolve_bndry);
};
std::transform(f2d.cbegin(), f2d.cend(), inserter, var_str_to_nvector);
std::transform(f3d.cbegin(), f3d.cend(), inserter, var_str_to_nvector);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "std::transform" is directly included [misc-include-cleaner]

    std::transform(f3d.cbegin(), f3d.cend(), inserter, var_str_to_nvector);
         ^

Comment thread src/solver/impls/arkode/arkode.hxx Outdated
return BoutNVector::create(ctx, subvectors);
}

void swap_state(const N_Vector u) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: 'u' declared with a const-qualified typedef; results in the type being '_generic_N_Vector *const' instead of 'const _generic_N_Vector *' [misc-misplaced-const]

  void swap_state(const N_Vector u) {
                                 ^
Additional context

/usr/include/sundials/sundials_nvector.h:91: typedef declared here

typedef _SUNDIALS_STRUCT_ _generic_N_Vector* N_Vector;
                                             ^

Comment thread src/solver/nvector.hxx
template <typename T>
struct Content {
T& field;
const bool own;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'own' of type 'const bool' is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

    const bool own;
               ^

Comment thread src/solver/nvector.hxx
struct Content {
T& field;
const bool own;
const bool evolve_bndry;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'evolve_bndry' of type 'const bool' is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

    const bool evolve_bndry;
               ^

Comment thread src/solver/nvector.hxx
T& field;
const bool own;
const bool evolve_bndry;
const sunindextype length;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'length' of type 'const sunindextype' (aka 'const int') is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

    const sunindextype length;
                       ^

Comment thread src/solver/nvector.hxx
length(all_reduce(getRegion().size())) {}

auto getRegion() const {
return field.getRegion(evolve_bndry ? RGN_ALL : RGN_NOBNDRY);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "RGN_ALL" is directly included [misc-include-cleaner]

src/solver/nvector.hxx:25:

- #include <bout/field.hxx>
+ #include "bout/bout_types.hxx"
+ #include <bout/field.hxx>

Comment thread src/solver/nvector.hxx
length(all_reduce(getRegion().size())) {}

auto getRegion() const {
return field.getRegion(evolve_bndry ? RGN_ALL : RGN_NOBNDRY);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "RGN_NOBNDRY" is directly included [misc-include-cleaner]

      return field.getRegion(evolve_bndry ? RGN_ALL : RGN_NOBNDRY);
                                                      ^

During configuration catch when `SUNDIALS::nvecmanyvector`
target is missing. If so, fall back to Sundials' nvector.
Defaults to using the Sundials parallel MPI implementation.
If Sundials many-vector is available then the custom N_Vector
is available.
Extends the Sundials ManyVector to Sundials and IDA interfaces.
To use, set
```
solver:nvector=manyvector
```

The default remains the MPI parallel N_Vector.
Backend can select either the Sundials MPI parallel N_Vector,
or a ManyVector backed by BOUT++ fields.
@bendudson bendudson changed the title WIP: Custom NVector for GPU Custom NVector for GPU Jun 17, 2026

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 34. Check the log or trigger a new build to see more.

.withDefault(false)),
nvector_type((*options)["nvector"]
.doc("N_Vector backend to use: sundials or manyvector")
.withDefault(NVectorType::Sundials)),

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "NVectorType" is directly included [misc-include-cleaner]

src/solver/impls/arkode/arkode.cxx:25:

- #include "bout/build_defines.hxx"
+ #include "/github/workspace/src/solver/sundials_nvector_interface.hxx"
+ #include "bout/build_defines.hxx"

use_jacobian((*options)["use_jacobian"].withDefault(false)),
nvector_type((*options)["nvector"]
.doc("N_Vector backend to use: sundials or manyvector")
.withDefault(NVectorType::Sundials)),

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "NVectorType" is directly included [misc-include-cleaner]

src/solver/impls/cvode/cvode.cxx:29:

+ #include "/github/workspace/src/solver/sundials_nvector_interface.hxx"

.withDefault(true)),
nvector_type((*options)["nvector"]
.doc("N_Vector backend to use: sundials or manyvector")
.withDefault(NVectorType::Sundials)),

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "NVectorType" is directly included [misc-include-cleaner]

src/solver/impls/ida/ida.cxx:32:

+ #include "/github/workspace/src/solver/sundials_nvector_interface.hxx"

Comment thread src/solver/nvector.hxx
#include <bout/field.hxx>
#include <bout/field2d.hxx>
#include <bout/field3d.hxx>
#include <functional>

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: included header field3d.hxx is not used directly [misc-include-cleaner]

Suggested change
#include <functional>
#include <functional>

Comment thread src/solver/nvector.hxx

class BoutNVector {
private:
static double all_reduce(const double local, const MPI_Op op = MPI_SUM) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "MPI_Op" is directly included [misc-include-cleaner]

src/solver/nvector.hxx:31:

- #include <nvector/nvector_manyvector.h>
+ #include <mpi.h>
+ #include <nvector/nvector_manyvector.h>

Comment thread src/solver/nvector.hxx

template <typename T, typename = bout::utils::EnableIfField<T>>
static T& get(const N_Vector v, std::size_t subvector) {
return BoutNVector::get<T>(N_VGetSubvector_ManyVector(v, subvector));

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: narrowing conversion from 'std::size_t' (aka 'unsigned long') to signed type 'sunindextype' (aka 'int') is implementation-defined [bugprone-narrowing-conversions]

    return BoutNVector::get<T>(N_VGetSubvector_ManyVector(v, subvector));
                                                             ^

#include "bout/solver.hxx"
#include "bout/sundials_backports.hxx"

#if BOUT_HAS_SUNDIALS_MANYVECTOR

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "BOUT_HAS_SUNDIALS_MANYVECTOR" is directly included [misc-include-cleaner]

src/solver/sundials_nvector_interface.hxx:30:

- #include "bout/solver.hxx"
+ #include "bout/build_defines.hxx"
+ #include "bout/solver.hxx"

#endif
}

auto* vec = callWithSUNContext(N_VNew_Parallel, ctx, BoutComm::get(), local_N, neq);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "BoutComm" is directly included [misc-include-cleaner]

src/solver/sundials_nvector_interface.hxx:30:

- #include "bout/solver.hxx"
+ #include "bout/field2d.hxx"
+ #include "bout/solver.hxx"

}

/// Clone a vector of the same backend/layout, with a descriptive error on failure.
N_Vector clone_vector_like(N_Vector source, const char* description) const {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: method 'clone_vector_like' can be made static [readability-convert-member-functions-to-static]

Suggested change
N_Vector clone_vector_like(N_Vector source, const char* description) const {
static N_Vector clone_vector_like(N_Vector source, const char* description) {


if (use_manyvector()) {
#if BOUT_HAS_SUNDIALS_MANYVECTOR
std::size_t i = 0;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "std::size_t" is directly included [misc-include-cleaner]

src/solver/sundials_nvector_interface.hxx:32:

+ #include <cstddef>

Includes guards to handle when ManyVector is unavailable.

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

clang-tidy made some suggestions

#if BOUT_HAS_SUNDIALS_MANYVECTOR
std::size_t i = 0;
for (std::size_t j = 0; j < f2d_values.size(); ++j, ++i) {
fill_field(BoutNVector::get<Field2D>(vec, i), f2d_values[j],

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "Field2D" is directly included [misc-include-cleaner]

        fill_field(BoutNVector::get<Field2D>(vec, i), f2d_values[j],
                                    ^

solver.f2d[j].evolve_bndry);
}
for (std::size_t j = 0; j < f3d_values.size(); ++j, ++i) {
fill_field(BoutNVector::get<Field3D>(vec, i), f3d_values[j],

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "Field3D" is directly included [misc-include-cleaner]

src/solver/sundials_nvector_interface.hxx:30:

- #include "bout/solver.hxx"
+ #include "bout/field3d.hxx"
+ #include "bout/solver.hxx"


BoutReal* option_data = N_VGetArrayPointer(vec);
int p = 0;
for (const auto& i2d : bout::globals::mesh->getRegion2D("RGN_BNDRY")) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "bout::globals::mesh" is directly included [misc-include-cleaner]

src/solver/sundials_nvector_interface.hxx:30:

- #include "bout/solver.hxx"
+ #include "bout/globals.hxx"
+ #include "bout/solver.hxx"

}

private:
Solver& solver;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'solver' of type 'Solver &' is a reference [cppcoreguidelines-avoid-const-or-ref-data-members]

  Solver& solver;
          ^


private:
Solver& solver;
const sundials::Context& ctx;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'ctx' of type 'const sundials::Context &' is a reference [cppcoreguidelines-avoid-const-or-ref-data-members]

  const sundials::Context& ctx;
                           ^

const sundials::Context& ctx;
NVectorType nvector_type;

void fill_vector_values_op(Ind2D UNUSED(i2d), BoutReal* option_data, int& p,

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "UNUSED" is directly included [misc-include-cleaner]

src/solver/sundials_nvector_interface.hxx:32:

+ #include "bout/unused.hxx"

}

/// Swap the active solver state fields with the subvectors contained in ``u``.
void swap_state(const N_Vector u) const {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: 'u' declared with a const-qualified typedef; results in the type being '_generic_N_Vector *const' instead of 'const _generic_N_Vector *' [misc-misplaced-const]

  void swap_state(const N_Vector u) const {
                                 ^
Additional context

/usr/include/sundials/sundials_nvector.h:91: typedef declared here

typedef _SUNDIALS_STRUCT_ _generic_N_Vector* N_Vector;
                                             ^

}

/// Swap the active solver derivative fields with the subvectors in ``du``.
void swap_deriv(const N_Vector du) const {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: 'du' declared with a const-qualified typedef; results in the type being '_generic_N_Vector *const' instead of 'const _generic_N_Vector *' [misc-misplaced-const]

  void swap_deriv(const N_Vector du) const {
                                 ^
Additional context

/usr/include/sundials/sundials_nvector.h:91: typedef declared here

typedef _SUNDIALS_STRUCT_ _generic_N_Vector* N_Vector;
                                             ^


template <typename FieldType>
static void fill_field(FieldType& field, BoutReal value, bool evolve_bndry) {
BOUT_FOR(i, field.getRegion(evolve_bndry ? RGN_ALL : RGN_NOBNDRY)) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "BOUT_FOR" is directly included [misc-include-cleaner]

src/solver/sundials_nvector_interface.hxx:30:

- #include "bout/solver.hxx"
+ #include "bout/region.hxx"
+ #include "bout/solver.hxx"

Conditionally include `nvector/nvector_manyvector.h`
Missing includes and some `const` annotations.

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

clang-tidy made some suggestions

Comment thread src/solver/nvector.hxx

class BoutNVector {
private:
static double all_reduce(const double local, const MPI_Op op = MPI_SUM) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "MPI_Op" is directly included [misc-include-cleaner]

src/solver/nvector.hxx:32:

- #if BOUT_HAS_SUNDIALS_MANYVECTOR
+ #include <mpi.h>
+ #if BOUT_HAS_SUNDIALS_MANYVECTOR

Comment thread src/solver/nvector.hxx
class BoutNVector {
private:
static double all_reduce(const double local, const MPI_Op op = MPI_SUM) {
double global;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: variable 'global' is not initialized [cppcoreguidelines-init-variables]

src/solver/nvector.hxx:33:

- #include <nvector/nvector_manyvector.h>
+ #include <math.h>
+ #include <nvector/nvector_manyvector.h>
Suggested change
double global;
double global = NAN;

Comment thread src/solver/nvector.hxx
if (x == nullptr) {
return;
}
delete &get_content<T>(x);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: deleting a pointer through a type that is not marked 'gsl::owner<>'; consider using a smart pointer instead [cppcoreguidelines-owning-memory]

      delete &get_content<T>(x);
      ^
Additional context

src/solver/nvector.hxx:125: variable declared here

    v->ops->nvdestroy = [](N_Vector x) {
                           ^

Comment thread src/solver/nvector.hxx
v->ops->nvmaxnorm = [](N_Vector x) { return max(abs(get_field<T>(x)), true); };

v->ops->nvwrmsnorm = [](N_Vector x, N_Vector y) {
return N_VWL2Norm(x, y) / std::sqrt(static_cast<BoutReal>(N_VGetLength(x)));

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "std::sqrt" is directly included [misc-include-cleaner]

src/solver/nvector.hxx:31:

- #include <functional>
+ #include <cmath>
+ #include <functional>

Comment thread src/solver/nvector.hxx
}

template <typename T, typename = bout::utils::EnableIfField<T>>
static void swap(const N_Vector v, T& field, std::size_t subvector) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "std::size_t" is directly included [misc-include-cleaner]

src/solver/nvector.hxx:31:

- #include <functional>
+ #include <cstddef>
+ #include <functional>


if (use_manyvector()) {
#if BOUT_HAS_SUNDIALS_MANYVECTOR
std::size_t i = 0;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "std::size_t" is directly included [misc-include-cleaner]

src/solver/sundials_nvector_interface.hxx:38:

+ #include <cstddef>

@bendudson bendudson merged commit 611ef8d into next Jun 18, 2026
22 of 23 checks passed
@bendudson bendudson deleted the custom-nvector branch June 18, 2026 01:39
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants