Reference documentation for deal.II version 9.5.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// Copyright (C) 2021 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#include <deal.II/base/config.h>
18
20
21#ifdef DEAL_II_WITH_SUNDIALS
22
24
28# include <deal.II/lac/vector.h>
29# ifdef DEAL_II_WITH_TRILINOS
32# endif
33# ifdef DEAL_II_WITH_PETSC
36# endif
37
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 arkode_linsol_get_type(SUNLinearSolver)
187 {
188 return SUNLINEARSOLVER_ITERATIVE;
189 }
190
191
192
193 template <typename VectorType>
194 int
195 arkode_linsol_solve(SUNLinearSolver LS,
196 SUNMatrix /*ignored*/,
197 N_Vector x,
198 N_Vector b,
199 realtype tol)
200 {
201 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
202
203 auto *src_b = unwrap_nvector_const<VectorType>(b);
204 auto *dst_x = unwrap_nvector<VectorType>(x);
205
207 content->a_times_fn
209 ,
210 content->linsol_ctx
211# endif
212 );
213
215 content->P_data,
216 content->preconditioner_solve,
218 content->linsol_ctx,
219# endif
220 tol);
221
222 return call_and_possibly_capture_exception(content->lsolve,
223 content->pending_exception,
224 op,
225 preconditioner,
226 *dst_x,
227 *src_b,
228 tol);
230
231
232
233 template <typename VectorType>
234 int
235 arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix /*ignored*/)
236 {
237 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
238
239 if (content->preconditioner_setup)
240 return content->preconditioner_setup(content->P_data);
241 return 0;
242 }
243
244
245
246 template <typename VectorType>
247 int arkode_linsol_initialize(SUNLinearSolver)
248 {
249 // this method is currently only provided because SUNDIALS 4.0.0 requires
250 // it - no user-set action is implemented so far
251 return 0;
252 }
253
254
255
256 template <typename VectorType>
257 int
258 arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes)
259 {
260 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
261
262 content->A_data = A_data;
263 content->a_times_fn = ATimes;
264 return 0;
265 }
266
267
268
269 template <typename VectorType>
270 int
271 arkode_linsol_set_preconditioner(SUNLinearSolver LS,
272 void * P_data,
273 PSetupFn p_setup,
274 PSolveFn p_solve)
275 {
276 LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
277
278 content->P_data = P_data;
279 content->preconditioner_setup = p_setup;
280 content->preconditioner_solve = p_solve;
281 return 0;
282 }
283 } // namespace
284
285
286
287 template <typename VectorType>
290 std::exception_ptr & pending_exception
292 ,
293 SUNContext linsol_ctx
294# endif
295 )
296 : content(
297 std::make_unique<LinearSolverContent<VectorType>>(pending_exception))
298 {
299# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
300 sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx);
301# else
302 sun_linear_solver = SUNLinSolNewEmpty();
303# endif
304
305 sun_linear_solver->ops->gettype = arkode_linsol_get_type;
306 sun_linear_solver->ops->solve = arkode_linsol_solve<VectorType>;
307 sun_linear_solver->ops->setup = arkode_linsol_setup<VectorType>;
308 sun_linear_solver->ops->initialize = arkode_linsol_initialize<VectorType>;
309 sun_linear_solver->ops->setatimes = arkode_linsol_set_a_times<VectorType>;
310 sun_linear_solver->ops->setpreconditioner =
311 arkode_linsol_set_preconditioner<VectorType>;
312
313 content->lsolve = lsolve;
314# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
315 content->linsol_ctx = linsol_ctx;
316# endif
317 sun_linear_solver->content = content.get();
318 }
319
320
321
322 namespace internal
323 {
324 template <typename VectorType>
326 {
327 SUNLinSolFreeEmpty(sun_linear_solver);
328 }
329
330
331
332 template <typename VectorType>
334 {
335 return sun_linear_solver;
336 }
337 } // namespace internal
338
339
340
341# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
342 template <typename VectorType>
344 ATimesFn a_times_fn,
345 SUNContext linsol_ctx)
346 : A_data(A_data)
347 , a_times_fn(a_times_fn)
348 , linsol_ctx(linsol_ctx)
349
350 {
351 Assert(a_times_fn != nullptr, ExcInternalError());
352 Assert(linsol_ctx != nullptr, ExcInternalError());
353 }
354# else
355 template <typename VectorType>
357 ATimesFn a_times_fn)
358 : A_data(A_data)
359 , a_times_fn(a_times_fn)
360
361 {
362 Assert(a_times_fn != nullptr, ExcInternalError());
363 }
364# endif
365
366
367
368 template <typename VectorType>
369 void
371 const VectorType &src) const
372 {
373 auto sun_dst = internal::make_nvector_view(dst
375 ,
376 linsol_ctx
377# endif
378 );
379 auto sun_src = internal::make_nvector_view(src
381 ,
382 linsol_ctx
383# endif
384 );
385 int status = a_times_fn(A_data, sun_src, sun_dst);
386 (void)status;
387 AssertSundialsSolver(status);
388 }
389
390
391
392# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
393 template <typename VectorType>
395 void * P_data,
396 PSolveFn p_solve_fn,
397 SUNContext linsol_ctx,
398 double tol)
399 : P_data(P_data)
400 , p_solve_fn(p_solve_fn)
401 , linsol_ctx(linsol_ctx)
402 , tol(tol)
403 {}
404# else
405 template <typename VectorType>
407 void *P_data,
408 PSolveFn p_solve_fn,
409 double tol)
410 : P_data(P_data)
411 , p_solve_fn(p_solve_fn)
412 , tol(tol)
413 {}
414# endif
415
416
417
418 template <typename VectorType>
419 void
421 const VectorType &src) const
422 {
423 // apply identity preconditioner if nothing else specified
424 if (!p_solve_fn)
425 {
426 dst = src;
427 return;
428 }
429
430 auto sun_dst = internal::make_nvector_view(dst
432 ,
433 linsol_ctx
434# endif
435 );
436 auto sun_src = internal::make_nvector_view(src
438 ,
439 linsol_ctx
440# endif
441 );
442 // for custom preconditioners no distinction between left and right
443 // preconditioning is made
444 int status =
445 p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
446 (void)status;
447 AssertSundialsSolver(status);
448 }
449
450 template struct SundialsOperator<Vector<double>>;
453 template struct SundialsOperator<
455
458 template struct SundialsPreconditioner<
460 template struct SundialsPreconditioner<
462
465 template class internal::LinearSolverWrapper<
467 template class internal::LinearSolverWrapper<
469
470# ifdef DEAL_II_WITH_MPI
471# ifdef DEAL_II_WITH_TRILINOS
474
477
479 template class internal::LinearSolverWrapper<
481# endif // DEAL_II_WITH_TRILINOS
482
483# ifdef DEAL_II_WITH_PETSC
484# ifndef PETSC_USE_COMPLEX
485
488
491
494# endif // PETSC_USE_COMPLEX
495# endif // DEAL_II_WITH_PETSC
496
497# endif // DEAL_II_WITH_MPI
498} // namespace SUNDIALS
500
501#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:472
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition config.h:359
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
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
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)