Reference documentation for deal.II version Git 191d06ed00 2021-05-11 21:15:49 -0400
\(\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\}}\)
sunlinsol_wrapper.cc
Go to the documentation of this file.
1 //-----------------------------------------------------------
2 //
3 // Copyright (C) 2021 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 #include <deal.II/base/config.h>
17 
19 
20 #ifdef DEAL_II_WITH_SUNDIALS
21 # if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
22 
23 # include <deal.II/base/exceptions.h>
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 
40 # if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
42 # endif
43 
45 
46 
47 
48 namespace SUNDIALS
49 {
51  int,
52  << "One of the SUNDIALS linear solver internal"
53  << " functions returned a negative error code: " << arg1
54  << ". Please consult SUNDIALS manual.");
55 
56 # define AssertSundialsSolver(code) \
57  Assert(code >= 0, ExcSundialsSolverError(code))
58 
59  namespace internal
60  {
64  template <typename VectorType>
66  {
67  ATimesFn a_times_fn;
70 
72 
73  void *P_data;
74  void *A_data;
75  };
76  } // namespace internal
77 
78  namespace
79  {
80  using namespace SUNDIALS::internal;
81 
86  template <typename VectorType>
88  access_content(SUNLinearSolver ls)
89  {
90  Assert(ls->content != nullptr, ExcInternalError());
91  return static_cast<LinearSolverContent<VectorType> *>(ls->content);
92  }
93 
94 
95 
96  SUNLinearSolver_Type arkode_linsol_get_type(SUNLinearSolver)
97  {
98  return SUNLINEARSOLVER_ITERATIVE;
99  }
100 
101 
102 
103  template <typename VectorType>
104  int
105  arkode_linsol_solve(SUNLinearSolver LS,
106  SUNMatrix /*ignored*/,
107  N_Vector x,
108  N_Vector b,
109  realtype tol)
110  {
111  auto content = access_content<VectorType>(LS);
112 
113  auto *src_b = unwrap_nvector_const<VectorType>(b);
114  auto *dst_x = unwrap_nvector<VectorType>(x);
115 
116  SundialsOperator<VectorType> op(content->A_data, content->a_times_fn);
117 
118  SundialsPreconditioner<VectorType> preconditioner(
119  content->P_data, content->preconditioner_solve, tol);
120 
121  return content->lsolve(op, preconditioner, *dst_x, *src_b, tol);
122  }
123 
124 
125 
126  template <typename VectorType>
127  int
128  arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix /*ignored*/)
129  {
130  auto content = access_content<VectorType>(LS);
131  if (content->preconditioner_setup)
132  return content->preconditioner_setup(content->P_data);
133  return 0;
134  }
135 
136 
137 
138  template <typename VectorType>
139  int arkode_linsol_initialize(SUNLinearSolver)
140  {
141  // this method is currently only provided because SUNDIALS 4.0.0 requires
142  // it - no user-set action is implemented so far
143  return 0;
144  }
145 
146 
147 
148  template <typename VectorType>
149  int
150  arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes)
151  {
152  auto content = access_content<VectorType>(LS);
153  content->A_data = A_data;
154  content->a_times_fn = ATimes;
155  return 0;
156  }
157 
158 
159 
160  template <typename VectorType>
161  int
162  arkode_linsol_set_preconditioner(SUNLinearSolver LS,
163  void * P_data,
164  PSetupFn p_setup,
165  PSolveFn p_solve)
166  {
167  auto content = access_content<VectorType>(LS);
168  content->P_data = P_data;
169  content->preconditioner_setup = p_setup;
170  content->preconditioner_solve = p_solve;
171  return 0;
172  }
173  } // namespace
174 
175 
176 
177  template <typename VectorType>
180  : content(std::make_unique<LinearSolverContent<VectorType>>())
181  {
183  sun_linear_solver->ops->gettype = arkode_linsol_get_type;
184  sun_linear_solver->ops->solve = arkode_linsol_solve<VectorType>;
185  sun_linear_solver->ops->setup = arkode_linsol_setup<VectorType>;
186  sun_linear_solver->ops->initialize = arkode_linsol_initialize<VectorType>;
187  sun_linear_solver->ops->setatimes = arkode_linsol_set_a_times<VectorType>;
188  sun_linear_solver->ops->setpreconditioner =
189  arkode_linsol_set_preconditioner<VectorType>;
190 
191  content->lsolve = lsolve;
192  sun_linear_solver->content = content.get();
193  }
194 
195 
196 
197  namespace internal
198  {
199  template <typename VectorType>
201  {
203  }
204 
205 
206 
207  template <typename VectorType>
209  {
210  return sun_linear_solver;
211  }
212  } // namespace internal
213 
214 
215 
216  template <typename VectorType>
218  ATimesFn a_times_fn)
219  : A_data(A_data)
220  , a_times_fn(a_times_fn)
221 
222  {
223  Assert(a_times_fn != nullptr, ExcInternalError());
224  }
225 
226 
227 
228  template <typename VectorType>
229  void
231  const VectorType &src) const
232  {
233  auto sun_dst = internal::make_nvector_view(dst);
234  auto sun_src = internal::make_nvector_view(src);
235  int status = a_times_fn(A_data, sun_src, sun_dst);
236  (void)status;
237  AssertSundialsSolver(status);
238  }
239 
240 
241 
242  template <typename VectorType>
244  void * P_data,
245  PSolveFn p_solve_fn,
246  double tol)
247  : P_data(P_data)
248  , p_solve_fn(p_solve_fn)
249  , tol(tol)
250  {}
251 
252 
253 
254  template <typename VectorType>
255  void
257  const VectorType &src) const
258  {
259  // apply identity preconditioner if nothing else specified
260  if (!p_solve_fn)
261  {
262  dst = src;
263  return;
264  }
265 
266  auto sun_dst = internal::make_nvector_view(dst);
267  auto sun_src = internal::make_nvector_view(src);
268  // for custom preconditioners no distinction between left and right
269  // preconditioning is made
270  int status =
271  p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
272  (void)status;
273  AssertSundialsSolver(status);
274  }
275 
276  template struct SundialsOperator<Vector<double>>;
277  template struct SundialsOperator<BlockVector<double>>;
279  template struct SundialsOperator<
281 
282  template struct SundialsPreconditioner<Vector<double>>;
284  template struct SundialsPreconditioner<
286  template struct SundialsPreconditioner<
287  LinearAlgebra::distributed::BlockVector<double>>;
288 
291  template class internal::LinearSolverWrapper<
292  LinearAlgebra::distributed::Vector<double>>;
293  template class internal::LinearSolverWrapper<
294  LinearAlgebra::distributed::BlockVector<double>>;
295 
296 # ifdef DEAL_II_WITH_MPI
297 # ifdef DEAL_II_WITH_TRILINOS
300 
303 
305  template class internal::LinearSolverWrapper<
307 # endif // DEAL_II_WITH_TRILINOS
308 
309 # ifdef DEAL_II_WITH_PETSC
310 # ifndef PETSC_USE_COMPLEX
311 
314 
317 
320 # endif // PETSC_USE_COMPLEX
321 # endif // DEAL_II_WITH_PETSC
322 
323 # endif // DEAL_II_WITH_MPI
324 } // namespace SUNDIALS
326 
327 # endif
328 #endif
LinearSolverWrapper(LinearSolveFunction< VectorType > lsolve)
#define AssertSundialsSolver(code)
LinearSolveFunction< VectorType > lsolve
void SUNLinSolFreeEmpty(SUNLinearSolver solver)
STL namespace.
static ::ExceptionBase & ExcSundialsSolverError(int arg1)
void vmult(VectorType &dst, const VectorType &src) const
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
void vmult(VectorType &dst, const VectorType &src) const
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SUNLinearSolver SUNLinSolNewEmpty()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
SundialsOperator(void *A_data, ATimesFn a_times_fn)
DeclException1(ExcARKodeError, int,<< "One of the SUNDIALS ARKode internal functions "<< " returned a negative error code: "<< arg1<< ". Please consult SUNDIALS manual.")
std::function< int(SundialsOperator< VectorType > &op, SundialsPreconditioner< VectorType > &prec, VectorType &x, const VectorType &b, double tol)> LinearSolveFunction
std::unique_ptr< LinearSolverContent< VectorType > > content
SundialsPreconditioner(void *P_data, PSolveFn p_solve_fn, double tol)
static ::ExceptionBase & ExcInternalError()