deal.II version GIT relicensing-2645-g9d42f38b55 2025-02-15 04:20:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
sunlinsol_wrapper.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2021 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#include <deal.II/base/config.h>
17
19
20#ifdef DEAL_II_WITH_SUNDIALS
21
23
27# include <deal.II/lac/vector.h>
28# ifdef DEAL_II_WITH_TRILINOS
31# endif
32# ifdef DEAL_II_WITH_PETSC
35# endif
36
39
40# include <sundials/sundials_iterative.h>
41# include <sundials/sundials_linearsolver.h>
42
44
45
46
47namespace SUNDIALS
48{
50 int,
51 << "One of the SUNDIALS linear solver internal"
52 << " functions returned a negative error code: " << arg1
53 << ". Please consult SUNDIALS manual.");
54
55# define AssertSundialsSolver(code) \
56 Assert(code >= 0, ExcSundialsSolverError(code))
57
58 namespace
59 {
72 template <typename F, typename... Args>
73 int
74 call_and_possibly_capture_exception(const F &f,
75 std::exception_ptr &eptr,
76 Args &&...args)
77 {
78 // See whether there is already something in the exception pointer
79 // variable. This can only happen if we had previously had
80 // a recoverable exception, and the underlying library actually
81 // did recover successfully. In that case, we can abandon the
82 // exception previously thrown. If eptr contains anything other,
83 // then we really don't know how that could have happened, and
84 // should probably bail out:
85 if (eptr)
86 {
87 try
88 {
89 std::rethrow_exception(eptr);
90 }
91 catch (const RecoverableUserCallbackError &)
92 {
93 // ok, ignore, but reset the pointer
94 eptr = nullptr;
95 }
96 catch (...)
97 {
98 // uh oh:
100 }
101 }
102
103 // Call the function and if that succeeds, return zero:
104 try
105 {
106 f(std::forward<Args>(args)...);
107 eptr = nullptr;
108 return 0;
109 }
110 // If the call failed with a recoverable error, then
111 // ignore the exception for now (but store a pointer to it)
112 // and return a positive return value (+1). If the underlying
113 // implementation manages to recover
114 catch (const RecoverableUserCallbackError &)
115 {
116 eptr = std::current_exception();
117 return 1;
118 }
119 // For any other exception, capture the exception and
120 // return -1:
121 catch (const std::exception &)
122 {
123 eptr = std::current_exception();
124 return -1;
125 }
126 }
127 } // namespace
128
129
130 namespace internal
131 {
135 template <typename VectorType>
137 {
139 : a_times_fn(nullptr)
140 , preconditioner_setup(nullptr)
141 , preconditioner_solve(nullptr)
143 , linsol_ctx(nullptr)
144# endif
145 , P_data(nullptr)
146 , A_data(nullptr)
148 {}
149
150# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
154
156# else
157 ATimesFn a_times_fn;
158 PSetupFn preconditioner_setup;
159 PSolveFn preconditioner_solve;
160# endif
161
163
164 void *P_data;
165 void *A_data;
166
171 std::exception_ptr &pending_exception;
172 };
173 } // namespace internal
174
175 namespace
176 {
177 using namespace SUNDIALS::internal;
178
183 template <typename VectorType>
185 access_content(SUNLinearSolver ls)
186 {
187 Assert(ls->content != nullptr, ExcInternalError());
188 return static_cast<LinearSolverContent<VectorType> *>(ls->content);
189 }
190
191
192
193 SUNLinearSolver_Type
194 arkode_linsol_get_type(SUNLinearSolver)
195 {
196 return SUNLINEARSOLVER_ITERATIVE;
197 }
198
199
200
201 template <typename VectorType>
202 int
203 arkode_linsol_solve(SUNLinearSolver LS,
204 SUNMatrix /*ignored*/,
205 N_Vector x,
206 N_Vector b,
208 {
209 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
210
211 auto *src_b = unwrap_nvector_const<VectorType>(b);
212 auto *dst_x = unwrap_nvector<VectorType>(x);
213
215 content->a_times_fn
217 ,
218 content->linsol_ctx
219# endif
220 );
221
223 content->P_data,
224 content->preconditioner_solve,
226 content->linsol_ctx,
227# endif
228 tol);
229
231 content->pending_exception,
232 op,
233 preconditioner,
234 *dst_x,
235 *src_b,
236 tol);
237 }
238
239
240
241 template <typename VectorType>
242 int
243 arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix /*ignored*/)
244 {
245 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
246
248 return content->preconditioner_setup(content->P_data);
249 return 0;
250 }
251
252
253
254 template <typename VectorType>
255 int
256 arkode_linsol_initialize(SUNLinearSolver)
257 {
258 // this method is currently only provided because SUNDIALS 4.0.0 requires
259 // it - no user-set action is implemented so far
260 return 0;
261 }
262
263
264 template <typename VectorType>
265 int
266 arkode_linsol_set_a_times(SUNLinearSolver LS,
267 void *A_data,
269 SUNATimesFn ATimes
270# else
271 ATimesFn ATimes
272# endif
273 )
274 {
275 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
276
277 content->A_data = A_data;
278 content->a_times_fn = ATimes;
279 return 0;
280 }
281
282
283
284 template <typename VectorType>
285 int
286 arkode_linsol_set_preconditioner(SUNLinearSolver LS,
287 void *P_data,
289 SUNPSetupFn p_setup,
290 SUNPSolveFn p_solve
291# else
292 PSetupFn p_setup,
293 PSolveFn p_solve
294# endif
295 )
296 {
297 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
298
299 content->P_data = P_data;
300 content->preconditioner_setup = p_setup;
301 content->preconditioner_solve = p_solve;
302 return 0;
303 }
304 } // namespace
305
306
307
308 template <typename VectorType>
311 std::exception_ptr &pending_exception
313 ,
314 SUNContext linsol_ctx
315# endif
316 )
317 : content(
318 std::make_unique<LinearSolverContent<VectorType>>(pending_exception))
319 {
320# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
321 sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx);
322# else
323 sun_linear_solver = SUNLinSolNewEmpty();
324# endif
325
326 sun_linear_solver->ops->gettype = arkode_linsol_get_type;
327 sun_linear_solver->ops->solve = arkode_linsol_solve<VectorType>;
328 sun_linear_solver->ops->setup = arkode_linsol_setup<VectorType>;
329 sun_linear_solver->ops->initialize = arkode_linsol_initialize<VectorType>;
330 sun_linear_solver->ops->setatimes = arkode_linsol_set_a_times<VectorType>;
331 sun_linear_solver->ops->setpreconditioner =
332 arkode_linsol_set_preconditioner<VectorType>;
333
334 content->lsolve = lsolve;
335# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
336 content->linsol_ctx = linsol_ctx;
337# endif
338 sun_linear_solver->content = content.get();
339 }
340
341
342
343 namespace internal
344 {
345 template <typename VectorType>
347 {
348 SUNLinSolFreeEmpty(sun_linear_solver);
349 }
350
351
352
353 template <typename VectorType>
355 {
356 return sun_linear_solver;
357 }
358 } // namespace internal
359
360
361
362# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
363 template <typename VectorType>
365 SUNATimesFn a_times_fn,
366 SUNContext linsol_ctx)
367 : A_data(A_data)
368 , a_times_fn(a_times_fn)
369 , linsol_ctx(linsol_ctx)
370
371 {
372 Assert(a_times_fn != nullptr, ExcInternalError());
373 Assert(linsol_ctx != nullptr, ExcInternalError());
374 }
375# else
376 template <typename VectorType>
378 ATimesFn a_times_fn)
379 : A_data(A_data)
380 , a_times_fn(a_times_fn)
381
382 {
383 Assert(a_times_fn != nullptr, ExcInternalError());
384 }
385# endif
386
387
388
389 template <typename VectorType>
390 void
392 const VectorType &src) const
393 {
394 auto sun_dst = internal::make_nvector_view(dst
396 ,
397 linsol_ctx
398# endif
399 );
400 auto sun_src = internal::make_nvector_view(src
402 ,
403 linsol_ctx
404# endif
405 );
406 int status = a_times_fn(A_data, sun_src, sun_dst);
407 (void)status;
408 AssertSundialsSolver(status);
409 }
410
411
412
413# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
414 template <typename VectorType>
416 void *P_data,
417 SUNPSolveFn p_solve_fn,
418 SUNContext linsol_ctx,
419 double tol)
420 : P_data(P_data)
421 , p_solve_fn(p_solve_fn)
422 , linsol_ctx(linsol_ctx)
423 , tol(tol)
424 {}
425# else
426 template <typename VectorType>
428 void *P_data,
429 PSolveFn p_solve_fn,
430 double tol)
431 : P_data(P_data)
432 , p_solve_fn(p_solve_fn)
433 , tol(tol)
434 {}
435# endif
436
437
438
439 template <typename VectorType>
440 void
442 const VectorType &src) const
443 {
444 // apply identity preconditioner if nothing else specified
445 if (!p_solve_fn)
446 {
447 dst = src;
448 return;
449 }
450
451 auto sun_dst = internal::make_nvector_view(dst
453 ,
454 linsol_ctx
455# endif
456 );
457 auto sun_src = internal::make_nvector_view(src
459 ,
460 linsol_ctx
461# endif
462 );
463 // for custom preconditioners no distinction between left and right
464 // preconditioning is made
465 int status =
466 p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
467 (void)status;
468 AssertSundialsSolver(status);
469 }
470
471 template struct SundialsOperator<Vector<double>>;
474 template struct SundialsOperator<
476
479 template struct SundialsPreconditioner<
481 template struct SundialsPreconditioner<
483
486 template class internal::LinearSolverWrapper<
488 template class internal::LinearSolverWrapper<
490
491# ifdef DEAL_II_WITH_MPI
492# ifdef DEAL_II_WITH_TRILINOS
495
498
500 template class internal::LinearSolverWrapper<
502# endif // DEAL_II_WITH_TRILINOS
503
504# ifdef DEAL_II_WITH_PETSC
505# ifndef PETSC_USE_COMPLEX
506
509
512
515# endif // PETSC_USE_COMPLEX
516# endif // DEAL_II_WITH_PETSC
517
518# endif // DEAL_II_WITH_MPI
519} // namespace SUNDIALS
521
522#endif
std::unique_ptr< LinearSolverContent< VectorType > > content
LinearSolverWrapper(const LinearSolveFunction< VectorType > &lsolve, std::exception_ptr &pending_exception, SUNContext linsol_ctx)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition config.h:422
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:513
static ::ExceptionBase & RecoverableUserCallbackError()
static ::ExceptionBase & ExcSundialsSolverError(int arg1)
#define AssertThrow(cond, exc)
int call_and_possibly_capture_exception(const F &f, std::exception_ptr &eptr, Args &&...args)
Definition utilities.h:48
std::function< void(SundialsOperator< VectorType > &op, SundialsPreconditioner< VectorType > &prec, VectorType &x, const VectorType &b, double tol)> LinearSolveFunction
::sunrealtype realtype
STL namespace.
SundialsOperator(void *A_data, SUNATimesFn a_times_fn, SUNContext linsol_ctx)
void vmult(VectorType &dst, const VectorType &src) const
SundialsPreconditioner(void *P_data, SUNPSolveFn p_solve_fn, SUNContext linsol_ctx, double tol)
void vmult(VectorType &dst, const VectorType &src) const
LinearSolveFunction< VectorType > lsolve
LinearSolverContent(std::exception_ptr &pending_exception)
#define AssertSundialsSolver(code)