Reference documentation for deal.II version GIT relicensing-12-g89ba5bf9db 2024-02-24 00:00:03+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 - 2023 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
150
151# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
153# endif
154
156
157 void *P_data;
158 void *A_data;
159
164 std::exception_ptr &pending_exception;
165 };
166 } // namespace internal
167
168 namespace
169 {
170 using namespace SUNDIALS::internal;
171
176 template <typename VectorType>
178 access_content(SUNLinearSolver ls)
179 {
180 Assert(ls->content != nullptr, ExcInternalError());
181 return static_cast<LinearSolverContent<VectorType> *>(ls->content);
182 }
183
184
185
186 SUNLinearSolver_Type
187 arkode_linsol_get_type(SUNLinearSolver)
188 {
189 return SUNLINEARSOLVER_ITERATIVE;
190 }
191
192
193
194 template <typename VectorType>
195 int
196 arkode_linsol_solve(SUNLinearSolver LS,
197 SUNMatrix /*ignored*/,
198 N_Vector x,
199 N_Vector b,
201 {
202 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
203
204 auto *src_b = unwrap_nvector_const<VectorType>(b);
205 auto *dst_x = unwrap_nvector<VectorType>(x);
206
208 content->a_times_fn
210 ,
211 content->linsol_ctx
212# endif
213 );
214
216 content->P_data,
217 content->preconditioner_solve,
219 content->linsol_ctx,
220# endif
221 tol);
222
223 return call_and_possibly_capture_exception(content->lsolve,
224 content->pending_exception,
225 op,
226 preconditioner,
227 *dst_x,
228 *src_b,
229 tol);
230 }
231
232
234 template <typename VectorType>
235 int
236 arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix /*ignored*/)
237 {
238 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
239
240 if (content->preconditioner_setup)
241 return content->preconditioner_setup(content->P_data);
242 return 0;
243 }
244
245
246
247 template <typename VectorType>
248 int
249 arkode_linsol_initialize(SUNLinearSolver)
250 {
251 // this method is currently only provided because SUNDIALS 4.0.0 requires
252 // it - no user-set action is implemented so far
253 return 0;
254 }
255
256
257
258 template <typename VectorType>
259 int
260 arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes)
261 {
262 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
263
264 content->A_data = A_data;
265 content->a_times_fn = ATimes;
266 return 0;
267 }
268
269
270
271 template <typename VectorType>
272 int
273 arkode_linsol_set_preconditioner(SUNLinearSolver LS,
274 void *P_data,
275 PSetupFn p_setup,
276 PSolveFn p_solve)
277 {
278 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
279
280 content->P_data = P_data;
281 content->preconditioner_setup = p_setup;
282 content->preconditioner_solve = p_solve;
283 return 0;
284 }
285 } // namespace
286
287
288
289 template <typename VectorType>
292 std::exception_ptr &pending_exception
294 ,
295 SUNContext linsol_ctx
296# endif
297 )
298 : content(
299 std::make_unique<LinearSolverContent<VectorType>>(pending_exception))
300 {
301# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
302 sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx);
303# else
304 sun_linear_solver = SUNLinSolNewEmpty();
305# endif
306
307 sun_linear_solver->ops->gettype = arkode_linsol_get_type;
308 sun_linear_solver->ops->solve = arkode_linsol_solve<VectorType>;
309 sun_linear_solver->ops->setup = arkode_linsol_setup<VectorType>;
310 sun_linear_solver->ops->initialize = arkode_linsol_initialize<VectorType>;
311 sun_linear_solver->ops->setatimes = arkode_linsol_set_a_times<VectorType>;
312 sun_linear_solver->ops->setpreconditioner =
313 arkode_linsol_set_preconditioner<VectorType>;
314
315 content->lsolve = lsolve;
316# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
317 content->linsol_ctx = linsol_ctx;
318# endif
319 sun_linear_solver->content = content.get();
320 }
321
322
323
324 namespace internal
325 {
326 template <typename VectorType>
328 {
329 SUNLinSolFreeEmpty(sun_linear_solver);
330 }
331
332
333
334 template <typename VectorType>
336 {
337 return sun_linear_solver;
338 }
339 } // namespace internal
340
341
342
343# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
344 template <typename VectorType>
346 ATimesFn a_times_fn,
347 SUNContext linsol_ctx)
348 : A_data(A_data)
349 , a_times_fn(a_times_fn)
350 , linsol_ctx(linsol_ctx)
351
352 {
353 Assert(a_times_fn != nullptr, ExcInternalError());
354 Assert(linsol_ctx != nullptr, ExcInternalError());
355 }
356# else
357 template <typename VectorType>
359 ATimesFn a_times_fn)
360 : A_data(A_data)
361 , a_times_fn(a_times_fn)
362
363 {
364 Assert(a_times_fn != nullptr, ExcInternalError());
365 }
366# endif
367
368
369
370 template <typename VectorType>
371 void
373 const VectorType &src) const
374 {
375 auto sun_dst = internal::make_nvector_view(dst
377 ,
378 linsol_ctx
379# endif
380 );
381 auto sun_src = internal::make_nvector_view(src
383 ,
384 linsol_ctx
385# endif
386 );
387 int status = a_times_fn(A_data, sun_src, sun_dst);
388 (void)status;
389 AssertSundialsSolver(status);
390 }
391
392
393
394# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
395 template <typename VectorType>
397 void *P_data,
398 PSolveFn p_solve_fn,
399 SUNContext linsol_ctx,
400 double tol)
401 : P_data(P_data)
402 , p_solve_fn(p_solve_fn)
403 , linsol_ctx(linsol_ctx)
404 , tol(tol)
405 {}
406# else
407 template <typename VectorType>
409 void *P_data,
410 PSolveFn p_solve_fn,
411 double tol)
412 : P_data(P_data)
413 , p_solve_fn(p_solve_fn)
414 , tol(tol)
415 {}
416# endif
417
418
419
420 template <typename VectorType>
421 void
423 const VectorType &src) const
424 {
425 // apply identity preconditioner if nothing else specified
426 if (!p_solve_fn)
427 {
428 dst = src;
429 return;
430 }
431
432 auto sun_dst = internal::make_nvector_view(dst
434 ,
435 linsol_ctx
436# endif
437 );
438 auto sun_src = internal::make_nvector_view(src
440 ,
441 linsol_ctx
442# endif
443 );
444 // for custom preconditioners no distinction between left and right
445 // preconditioning is made
446 int status =
447 p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
448 (void)status;
449 AssertSundialsSolver(status);
450 }
451
452 template struct SundialsOperator<Vector<double>>;
455 template struct SundialsOperator<
457
460 template struct SundialsPreconditioner<
462 template struct SundialsPreconditioner<
464
467 template class internal::LinearSolverWrapper<
469 template class internal::LinearSolverWrapper<
471
472# ifdef DEAL_II_WITH_MPI
473# ifdef DEAL_II_WITH_TRILINOS
476
479
481 template class internal::LinearSolverWrapper<
483# endif // DEAL_II_WITH_TRILINOS
484
485# ifdef DEAL_II_WITH_PETSC
486# ifndef PETSC_USE_COMPLEX
487
490
493
496# endif // PETSC_USE_COMPLEX
497# endif // DEAL_II_WITH_PETSC
498
499# endif // DEAL_II_WITH_MPI
500} // namespace SUNDIALS
502
503#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:502
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition config.h:403
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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)
std::function< void(SundialsOperator< VectorType > &op, SundialsPreconditioner< VectorType > &prec, VectorType &x, const VectorType &b, double tol)> LinearSolveFunction
::sunrealtype realtype
STL namespace.
void vmult(VectorType &dst, const VectorType &src) const
SundialsOperator(void *A_data, ATimesFn a_times_fn, SUNContext linsol_ctx)
SundialsPreconditioner(void *P_data, PSolveFn 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)