NeoN
A framework for CFD software
Loading...
Searching...
No Matches
sundials.hpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2023 - 2025 NeoN authors
2//
3// SPDX-License-Identifier: MIT
4
5#pragma once
6
7#if NN_WITH_SUNDIALS
8
9#include <concepts>
10#include <functional>
11#include <memory>
12
13#include <sundials/sundials_nvector.h>
14#include <sundials/sundials_core.hpp>
15#include <nvector/nvector_serial.h>
16#include <nvector/nvector_kokkos.hpp>
17#include <arkode/arkode_arkstep.h>
18#include <arkode/arkode_erkstep.h>
19
20#include "NeoN/core/error.hpp"
22#include "NeoN/fields/field.hpp"
23
24namespace NeoN::sundials
25{
26
32inline auto SUN_CONTEXT_DELETER = [](SUNContext* ctx)
33{
34 if (ctx != nullptr)
35 {
36 SUNContext_Free(ctx);
37 }
38};
39
45inline auto SUN_ARK_DELETER = [](char* ark)
46{
47 if (ark != nullptr)
48 {
49 void* arkodMem = reinterpret_cast<void*>(ark);
50 ARKodeFree(&arkodMem);
51 }
52};
53
60inline ARKODE_ERKTableID stringToERKTable(const std::string& key)
61{
62 if (key == "Forward-Euler") return ARKODE_FORWARD_EULER_1_1;
63 if (key == "Heun")
64 {
65 NF_ERROR_EXIT("Currently unsupported until field time step-stage indexing resolved.");
66 return ARKODE_HEUN_EULER_2_1_2;
67 }
68 if (key == "Midpoint")
69 {
70 NF_ERROR_EXIT("Currently unsupported until field time step-stage indexing resolved.");
71 return ARKODE_EXPLICIT_MIDPOINT_EULER_2_1_2;
72 }
74 "Unsupported Runge-Kutta time integration method selectied: " + key + ".\n"
75 + "Supported methods are: Forward-Euler, Heun, Midpoint."
76 );
77 return ARKODE_ERK_NONE; // avoids compiler warnings.
78}
79
88template<typename SKVectorType, typename ValueType>
89void fieldToSunNVectorImpl(const NeoN::Vector<ValueType>& field, N_Vector& vector)
90{
91 auto view = ::sundials::kokkos::GetVec<SKVectorType>(vector)->View();
92 auto fieldView = field.view();
94 field.exec(), field.range(), KOKKOS_LAMBDA(const localIdx i) { view(i) = fieldView[i]; }
95 );
96};
97
105template<typename ValueType>
106void fieldToSunNVector(const NeoN::Vector<ValueType>& field, N_Vector& vector)
107{
108 // CHECK FOR N_Vector on correct space in DEBUG
109 if (std::holds_alternative<NeoN::GPUExecutor>(field.exec()))
110 {
111 fieldToSunNVectorImpl<::sundials::kokkos::Vector<Kokkos::DefaultExecutionSpace>>(
112 field, vector
113 );
114 return;
115 }
116 if (std::holds_alternative<NeoN::CPUExecutor>(field.exec()))
117 {
118 fieldToSunNVectorImpl<::sundials::kokkos::Vector<Kokkos::DefaultHostExecutionSpace>>(
119 field, vector
120 );
121 return;
122 }
123 if (std::holds_alternative<NeoN::SerialExecutor>(field.exec()))
124 {
125 fieldToSunNVectorImpl<::sundials::kokkos::Vector<Kokkos::Serial>>(field, vector);
126 return;
127 }
128 NF_ERROR_EXIT("Unsupported NeoN executor for field.");
129};
130
139template<typename SKVectorType, typename ValueType>
140void sunNVectorToVectorImpl(const N_Vector& vector, NeoN::Vector<ValueType>& field)
141{
142 auto view = ::sundials::kokkos::GetVec<SKVectorType>(vector)->View();
143 ValueType* fieldData = field.data();
145 field.exec(), field.range(), KOKKOS_LAMBDA(const localIdx i) { fieldData[i] = view(i); }
146 );
147};
148
155template<typename ValueType>
156void sunNVectorToVector(const N_Vector& vector, NeoN::Vector<ValueType>& field)
157{
158 if (std::holds_alternative<NeoN::GPUExecutor>(field.exec()))
159 {
160 sunNVectorToVectorImpl<::sundials::kokkos::Vector<Kokkos::DefaultExecutionSpace>>(
161 vector, field
162 );
163 return;
164 }
165 if (std::holds_alternative<NeoN::CPUExecutor>(field.exec()))
166 {
167 sunNVectorToVectorImpl<::sundials::kokkos::Vector<Kokkos::DefaultHostExecutionSpace>>(
168 vector, field
169 );
170 return;
171 }
172 if (std::holds_alternative<NeoN::SerialExecutor>(field.exec()))
173 {
174 sunNVectorToVectorImpl<::sundials::kokkos::Vector<Kokkos::Serial>>(vector, field);
175 return;
176 }
177 NF_ERROR_EXIT("Unsupported NeoN executor for field.");
178};
179
194template<typename SolutionVectorType>
195int explicitRKSolve([[maybe_unused]] sunrealtype t, N_Vector y, N_Vector ydot, void* userData)
196{
197 // Pointer wrangling
198 using ValueType = typename SolutionVectorType::VectorValueType;
200 reinterpret_cast<NeoN::dsl::Expression<ValueType>*>(userData);
201 sunrealtype* yDotArray = N_VGetArrayPointer(ydot);
202 sunrealtype* yArray = N_VGetArrayPointer(y);
203
204 NF_ASSERT(
205 yDotArray != nullptr && yArray != nullptr && pdeExpre != nullptr,
206 "Failed to dereference pointers in sundails."
207 );
208
209 auto size = static_cast<localIdx>(N_VGetLength(y));
210 // Copy initial value from y to source.
211 NeoN::Vector<NeoN::scalar> source = pdeExpre->explicitOperation(size) * -1.0; // compute spatial
212 if (std::holds_alternative<NeoN::GPUExecutor>(pdeExpre->exec()))
213 {
214 Kokkos::fence();
215 }
216 NeoN::sundials::fieldToSunNVector(source, ydot); // assign rhs to ydot.
217 return 0;
218}
219
220namespace detail
221{
222
230template<typename Vector>
231void initNVector(size_t size, std::shared_ptr<SUNContext> context, Vector& vec)
232{
233 vec.initNVector(size, context);
234}
235
242template<typename Vector>
243const N_Vector& sunNVector(const Vector& vec)
244{
245 return vec.sunNVector();
246}
247
254template<typename Vector>
255N_Vector& sunNVector(Vector& vec)
256{
257 return vec.sunNVector();
258}
259}
260
266template<typename ValueType>
267class SKVectorSerial
268{
269public:
270
271 SKVectorSerial() {};
272 ~SKVectorSerial() = default;
273 SKVectorSerial(const SKVectorSerial& other)
274 : kvector_(other.kvector_), svector_(other.kvector_) {};
275 SKVectorSerial(SKVectorSerial&& other) noexcept
276 : kvector_(std::move(other.kvector_)), svector_(std::move(other.svector_)) {};
277 SKVectorSerial& operator=(const SKVectorSerial& other) = delete;
278 SKVectorSerial& operator=(SKVectorSerial&& other) = delete;
279
280
281 using KVector = ::sundials::kokkos::Vector<Kokkos::Serial>;
282 void initNVector(size_t size, std::shared_ptr<SUNContext> context)
283 {
284 kvector_ = KVector(size, *context);
285 svector_ = kvector_;
286 };
287 const N_Vector& sunNVector() const { return svector_; };
288 N_Vector& sunNVector() { return svector_; };
289
290private:
291
292 KVector kvector_ {};
293 N_Vector svector_ {nullptr};
294};
295
301template<typename ValueType>
302class SKVectorHostDefault
303{
304public:
305
306 using KVector = ::sundials::kokkos::Vector<Kokkos::DefaultHostExecutionSpace>;
307
308 SKVectorHostDefault() = default;
309 ~SKVectorHostDefault() = default;
310 SKVectorHostDefault(const SKVectorHostDefault& other)
311 : kvector_(other.kvector_), svector_(other.kvector_) {};
312 SKVectorHostDefault(SKVectorHostDefault&& other) noexcept
313 : kvector_(std::move(other.kvector_)), svector_(std::move(other.svector_)) {};
314 SKVectorHostDefault& operator=(const SKVectorHostDefault& other) = delete;
315 SKVectorHostDefault& operator=(SKVectorHostDefault&& other) = delete;
316
317 void initNVector(size_t size, std::shared_ptr<SUNContext> context)
318 {
319 kvector_ = KVector(size, *context);
320 svector_ = kvector_;
321 };
322 const N_Vector& sunNVector() const { return svector_; };
323 N_Vector& sunNVector() { return svector_; };
324
325private:
326
327 KVector kvector_ {};
328 N_Vector svector_ {nullptr};
329};
330
336template<typename ValueType>
337class SKVectorDefault
338{
339public:
340
341 using KVector = ::sundials::kokkos::Vector<Kokkos::DefaultExecutionSpace>;
342
343 SKVectorDefault() = default;
344 ~SKVectorDefault() = default;
345 SKVectorDefault(const SKVectorDefault& other)
346 : kvector_(other.kvector_), svector_(other.kvector_) {};
347 SKVectorDefault(SKVectorDefault&& other) noexcept
348 : kvector_(std::move(other.kvector_)), svector_(std::move(other.svector_)) {};
349 SKVectorDefault& operator=(const SKVectorDefault& other) = delete;
350 SKVectorDefault& operator=(SKVectorDefault&& other) = delete;
351
352 void initNVector(size_t size, std::shared_ptr<SUNContext> context)
353 {
354 kvector_ = KVector(size, *context);
355 svector_ = kvector_;
356 };
357
358 const N_Vector& sunNVector() const { return svector_; };
359
360 N_Vector& sunNVector() { return svector_; };
361
362private:
363
364 KVector kvector_ {};
365 N_Vector svector_ {nullptr};
366};
367
374template<typename ValueType>
375class SKVector
376{
377public:
378
379 using SKVectorSerialV = SKVectorSerial<ValueType>;
380 using SKVectorHostDefaultV = SKVectorHostDefault<ValueType>;
381 using SKDefaultVectorV = SKVectorDefault<ValueType>;
382 using SKVectorVariant = std::variant<SKVectorSerialV, SKVectorHostDefaultV, SKDefaultVectorV>;
383
387 SKVector() { vector_.template emplace<SKVectorHostDefaultV>(); };
388
392 ~SKVector() = default;
393
398 SKVector(const SKVector&) = default;
399
403 SKVector& operator=(const SKVector&) = delete;
404
409 SKVector(SKVector&&) noexcept = default;
410
414 SKVector& operator=(SKVector&&) noexcept = delete;
415
420 void setExecutor(const NeoN::Executor& exec)
421 {
422 if (std::holds_alternative<NeoN::GPUExecutor>(exec))
423 {
424 vector_.template emplace<SKDefaultVectorV>();
425 return;
426 }
427 if (std::holds_alternative<NeoN::CPUExecutor>(exec))
428 {
429 vector_.template emplace<SKVectorHostDefaultV>();
430 return;
431 }
432 if (std::holds_alternative<NeoN::SerialExecutor>(exec))
433 {
434 vector_.template emplace<SKVectorSerialV>();
435 return;
436 }
437
439 "Unsupported NeoN executor: "
440 << std::visit([](const auto& e) { return e.name(); }, exec) << "."
441 );
442 }
443
449 void initNVector(size_t size, std::shared_ptr<SUNContext> context)
450 {
451 std::visit(
452 [size, &context](auto& vec) { detail::initNVector(size, context, vec); }, vector_
453 );
454 }
455
460 const N_Vector& sunNVector() const
461 {
462 return std::visit(
463 [](const auto& vec) -> const N_Vector& { return detail::sunNVector(vec); }, vector_
464 );
465 }
466
471 N_Vector& sunNVector()
472 {
473 return std::visit([](auto& vec) -> N_Vector& { return detail::sunNVector(vec); }, vector_);
474 }
475
480 const SKVectorVariant& variant() const { return vector_; }
481
486 SKVectorVariant& variant() { return vector_; }
487
488private:
489
490 SKVectorVariant vector_;
491};
492}
493
494#endif
A class to contain the data and executors for a field and define some basic operations.
Definition vector.hpp:30
ValueType * data()
Direct access to the underlying field data.
Definition vector.hpp:217
std::pair< localIdx, localIdx > range() const
Gets the range of the field.
Definition vector.hpp:305
const Executor & exec() const
Gets the executor associated with the field.
Definition vector.hpp:229
View< ValueType > view() &&=delete
Vector< ValueType > explicitOperation(localIdx nCells) const
const Executor & exec() const
#define NF_ERROR_EXIT(message)
Macro for printing an error message and aborting the program.
Definition error.hpp:110
#define NF_ASSERT(condition, message)
Macro for asserting a condition and printing an error message if the condition is false.
Definition error.hpp:144
SpatialOperator< scalar > source(fvcc::VolumeField< scalar > &coeff, fvcc::VolumeField< scalar > &phi)
Definition array.hpp:20
void parallelFor(const Executor &exec, std::pair< localIdx, localIdx > range, Kernel kernel, std::string name="parallelFor")
int32_t localIdx
Definition label.hpp:32