Reference documentation for deal.II version 9.6.0
\(\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
41
42
43
44namespace SUNDIALS
45{
47 int,
48 << "One of the SUNDIALS linear solver internal"
49 << " functions returned a negative error code: " << arg1
50 << ". Please consult SUNDIALS manual.");
51
52# define AssertSundialsSolver(code) \
53 Assert(code >= 0, ExcSundialsSolverError(code))
54
55 namespace
56 {
69 template <typename F, typename... Args>
70 int
71 call_and_possibly_capture_exception(const F &f,
72 std::exception_ptr &eptr,
73 Args &&...args)
74 {
75 // See whether there is already something in the exception pointer
76 // variable. This can only happen if we had previously had
77 // a recoverable exception, and the underlying library actually
78 // did recover successfully. In that case, we can abandon the
79 // exception previously thrown. If eptr contains anything other,
80 // then we really don't know how that could have happened, and
81 // should probably bail out:
82 if (eptr)
83 {
84 try
85 {
86 std::rethrow_exception(eptr);
87 }
88 catch (const RecoverableUserCallbackError &)
89 {
90 // ok, ignore, but reset the pointer
91 eptr = nullptr;
92 }
93 catch (...)
94 {
95 // uh oh:
97 }
98 }
99
100 // Call the function and if that succeeds, return zero:
101 try
102 {
103 f(std::forward<Args>(args)...);
104 eptr = nullptr;
105 return 0;
106 }
107 // If the call failed with a recoverable error, then
108 // ignore the exception for now (but store a pointer to it)
109 // and return a positive return value (+1). If the underlying
110 // implementation manages to recover
111 catch (const RecoverableUserCallbackError &)
112 {
113 eptr = std::current_exception();
114 return 1;
115 }
116 // For any other exception, capture the exception and
117 // return -1:
118 catch (const std::exception &)
119 {
120 eptr = std::current_exception();
121 return -1;
122 }
123 }
124 } // namespace
125
126
127 namespace internal
128 {
132 template <typename VectorType>
134 {
136 : a_times_fn(nullptr)
137 , preconditioner_setup(nullptr)
138 , preconditioner_solve(nullptr)
140 , linsol_ctx(nullptr)
141# endif
142 , P_data(nullptr)
143 , A_data(nullptr)
145 {}
146
147# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
151
153# else
154 ATimesFn a_times_fn;
155 PSetupFn preconditioner_setup;
156 PSolveFn preconditioner_solve;
157# endif
158
160
161 void *P_data;
162 void *A_data;
163
168 std::exception_ptr &pending_exception;
169 };
170 } // namespace internal
171
172 namespace
173 {
174 using namespace SUNDIALS::internal;
175
180 template <typename VectorType>
182 access_content(SUNLinearSolver ls)
183 {
184 Assert(ls->content != nullptr, ExcInternalError());
185 return static_cast<LinearSolverContent<VectorType> *>(ls->content);
186 }
187
188
189
190 SUNLinearSolver_Type
191 arkode_linsol_get_type(SUNLinearSolver)
192 {
193 return SUNLINEARSOLVER_ITERATIVE;
194 }
195
196
197
198 template <typename VectorType>
199 int
200 arkode_linsol_solve(SUNLinearSolver LS,
201 SUNMatrix /*ignored*/,
202 N_Vector x,
203 N_Vector b,
205 {
206 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
207
208 auto *src_b = unwrap_nvector_const<VectorType>(b);
209 auto *dst_x = unwrap_nvector<VectorType>(x);
210
212 content->a_times_fn
214 ,
215 content->linsol_ctx
216# endif
217 );
218
220 content->P_data,
221 content->preconditioner_solve,
223 content->linsol_ctx,
224# endif
225 tol);
226
228 content->pending_exception,
229 op,
230 preconditioner,
231 *dst_x,
232 *src_b,
233 tol);
234 }
235
236
237
238 template <typename VectorType>
239 int
240 arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix /*ignored*/)
241 {
242 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
243
244 if (content->preconditioner_setup)
245 return content->preconditioner_setup(content->P_data);
246 return 0;
248
249
250
251 template <typename VectorType>
252 int
253 arkode_linsol_initialize(SUNLinearSolver)
254 {
255 // this method is currently only provided because SUNDIALS 4.0.0 requires
256 // it - no user-set action is implemented so far
257 return 0;
258 }
259
260
261 template <typename VectorType>
262 int
263 arkode_linsol_set_a_times(SUNLinearSolver LS,
264 void *A_data,
266 SUNATimesFn ATimes
267# else
268 ATimesFn ATimes
269# endif
270 )
271 {
272 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
273
274 content->A_data = A_data;
275 content->a_times_fn = ATimes;
276 return 0;
277 }
278
279
280
281 template <typename VectorType>
282 int
283 arkode_linsol_set_preconditioner(SUNLinearSolver LS,
284 void *P_data,
286 SUNPSetupFn p_setup,
287 SUNPSolveFn p_solve
288# else
289 PSetupFn p_setup,
290 PSolveFn p_solve
291# endif
292 )
293 {
294 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
295
296 content->P_data = P_data;
297 content->preconditioner_setup = p_setup;
298 content->preconditioner_solve = p_solve;
299 return 0;
300 }
301 } // namespace
302
303
304
305 template <typename VectorType>
308 std::exception_ptr &pending_exception
310 ,
311 SUNContext linsol_ctx
312# endif
313 )
314 : content(
315 std::make_unique<LinearSolverContent<VectorType>>(pending_exception))
316 {
317# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
318 sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx);
319# else
320 sun_linear_solver = SUNLinSolNewEmpty();
321# endif
322
323 sun_linear_solver->ops->gettype = arkode_linsol_get_type;
324 sun_linear_solver->ops->solve = arkode_linsol_solve<VectorType>;
325 sun_linear_solver->ops->setup = arkode_linsol_setup<VectorType>;
326 sun_linear_solver->ops->initialize = arkode_linsol_initialize<VectorType>;
327 sun_linear_solver->ops->setatimes = arkode_linsol_set_a_times<VectorType>;
328 sun_linear_solver->ops->setpreconditioner =
329 arkode_linsol_set_preconditioner<VectorType>;
330
331 content->lsolve = lsolve;
332# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
333 content->linsol_ctx = linsol_ctx;
334# endif
335 sun_linear_solver->content = content.get();
336 }
337
338
339
340 namespace internal
341 {
342 template <typename VectorType>
344 {
345 SUNLinSolFreeEmpty(sun_linear_solver);
346 }
347
348
349
350 template <typename VectorType>
352 {
353 return sun_linear_solver;
354 }
355 } // namespace internal
356
357
358
359# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
360 template <typename VectorType>
362 SUNATimesFn a_times_fn,
363 SUNContext linsol_ctx)
364 : A_data(A_data)
365 , a_times_fn(a_times_fn)
366 , linsol_ctx(linsol_ctx)
367
368 {
369 Assert(a_times_fn != nullptr, ExcInternalError());
370 Assert(linsol_ctx != nullptr, ExcInternalError());
371 }
372# else
373 template <typename VectorType>
375 ATimesFn a_times_fn)
376 : A_data(A_data)
377 , a_times_fn(a_times_fn)
378
379 {
380 Assert(a_times_fn != nullptr, ExcInternalError());
381 }
382# endif
383
384
385
386 template <typename VectorType>
387 void
389 const VectorType &src) const
390 {
391 auto sun_dst = internal::make_nvector_view(dst
393 ,
394 linsol_ctx
395# endif
396 );
397 auto sun_src = internal::make_nvector_view(src
399 ,
400 linsol_ctx
401# endif
402 );
403 int status = a_times_fn(A_data, sun_src, sun_dst);
404 (void)status;
405 AssertSundialsSolver(status);
406 }
407
408
409
410# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
411 template <typename VectorType>
413 void *P_data,
414 SUNPSolveFn p_solve_fn,
415 SUNContext linsol_ctx,
416 double tol)
417 : P_data(P_data)
418 , p_solve_fn(p_solve_fn)
419 , linsol_ctx(linsol_ctx)
420 , tol(tol)
421 {}
422# else
423 template <typename VectorType>
425 void *P_data,
426 PSolveFn p_solve_fn,
427 double tol)
428 : P_data(P_data)
429 , p_solve_fn(p_solve_fn)
430 , tol(tol)
431 {}
432# endif
433
434
435
436 template <typename VectorType>
437 void
439 const VectorType &src) const
440 {
441 // apply identity preconditioner if nothing else specified
442 if (!p_solve_fn)
443 {
444 dst = src;
445 return;
446 }
447
448 auto sun_dst = internal::make_nvector_view(dst
450 ,
451 linsol_ctx
452# endif
453 );
454 auto sun_src = internal::make_nvector_view(src
456 ,
457 linsol_ctx
458# endif
459 );
460 // for custom preconditioners no distinction between left and right
461 // preconditioning is made
462 int status =
463 p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
464 (void)status;
465 AssertSundialsSolver(status);
466 }
467
468 template struct SundialsOperator<Vector<double>>;
471 template struct SundialsOperator<
473
476 template struct SundialsPreconditioner<
478 template struct SundialsPreconditioner<
480
483 template class internal::LinearSolverWrapper<
485 template class internal::LinearSolverWrapper<
487
488# ifdef DEAL_II_WITH_MPI
489# ifdef DEAL_II_WITH_TRILINOS
492
495
497 template class internal::LinearSolverWrapper<
499# endif // DEAL_II_WITH_TRILINOS
500
501# ifdef DEAL_II_WITH_PETSC
502# ifndef PETSC_USE_COMPLEX
503
506
509
512# endif // PETSC_USE_COMPLEX
513# endif // DEAL_II_WITH_PETSC
514
515# endif // DEAL_II_WITH_MPI
516} // namespace SUNDIALS
518
519#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:503
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition config.h:404
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
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:45
VectorType * unwrap_nvector(N_Vector v)
const VectorType * unwrap_nvector_const(N_Vector v)
::sunrealtype realtype
std::function< void(SundialsOperator< VectorType > &op, SundialsPreconditioner< VectorType > &prec, VectorType &x, const VectorType &b, double tol)> LinearSolveFunction
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)