Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
kinsol.cc
1 //-----------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2018 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 
19 #include <deal.II/sundials/kinsol.h>
20 
21 #ifdef DEAL_II_WITH_SUNDIALS
22 
23 # include <deal.II/base/utilities.h>
24 
25 # include <deal.II/lac/block_vector.h>
26 # ifdef DEAL_II_WITH_TRILINOS
27 # include <deal.II/lac/trilinos_parallel_block_vector.h>
28 # include <deal.II/lac/trilinos_vector.h>
29 # endif
30 # ifdef DEAL_II_WITH_PETSC
31 # include <deal.II/lac/petsc_block_vector.h>
32 # include <deal.II/lac/petsc_vector.h>
33 # endif
34 # include <deal.II/base/utilities.h>
35 
36 # include <deal.II/sundials/copy.h>
37 
38 # include <sundials/sundials_config.h>
39 # if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0)
40 # include <kinsol/kinsol_direct.h>
41 # include <sunlinsol/sunlinsol_dense.h>
42 # include <sunmatrix/sunmatrix_dense.h>
43 # else
44 # include <kinsol/kinsol_dense.h>
45 # endif
46 
47 # include <iomanip>
48 # include <iostream>
49 
50 DEAL_II_NAMESPACE_OPEN
51 
52 namespace SUNDIALS
53 {
54  using namespace internal;
55 
56  namespace
57  {
58  template <typename VectorType>
59  int
60  t_kinsol_function(N_Vector yy, N_Vector FF, void *user_data)
61  {
62  KINSOL<VectorType> &solver =
63  *static_cast<KINSOL<VectorType> *>(user_data);
64  GrowingVectorMemory<VectorType> mem;
65 
66  typename VectorMemory<VectorType>::Pointer src_yy(mem);
67  solver.reinit_vector(*src_yy);
68 
69  typename VectorMemory<VectorType>::Pointer dst_FF(mem);
70  solver.reinit_vector(*dst_FF);
71 
72  copy(*src_yy, yy);
73 
74  int err = 0;
75  if (solver.residual)
76  err = solver.residual(*src_yy, *dst_FF);
77  else if (solver.iteration_function)
78  err = solver.iteration_function(*src_yy, *dst_FF);
79  else
80  Assert(false, ExcInternalError());
81 
82  copy(FF, *dst_FF);
83 
84  return err;
85  }
86 
87 
88 
89  template <typename VectorType>
90  int
91  t_kinsol_setup_jacobian(KINMem kinsol_mem)
92  {
93  KINSOL<VectorType> &solver =
94  *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
95  GrowingVectorMemory<VectorType> mem;
96 
97  typename VectorMemory<VectorType>::Pointer src_ycur(mem);
98  solver.reinit_vector(*src_ycur);
99 
100  typename VectorMemory<VectorType>::Pointer src_fcur(mem);
101  solver.reinit_vector(*src_fcur);
102 
103  copy(*src_ycur, kinsol_mem->kin_uu);
104  copy(*src_fcur, kinsol_mem->kin_fval);
105 
106  int err = solver.setup_jacobian(*src_ycur, *src_fcur);
107  return err;
108  }
109 
110 
111 
112  template <typename VectorType>
113  int
114  t_kinsol_solve_jacobian(KINMem kinsol_mem,
115  N_Vector x,
116  N_Vector b,
117  realtype *sJpnorm,
118  realtype *sFdotJp)
119  {
120  KINSOL<VectorType> &solver =
121  *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
122  GrowingVectorMemory<VectorType> mem;
123 
124  typename VectorMemory<VectorType>::Pointer src_ycur(mem);
125  solver.reinit_vector(*src_ycur);
126 
127  typename VectorMemory<VectorType>::Pointer src_fcur(mem);
128  solver.reinit_vector(*src_fcur);
129 
130  copy(*src_ycur, kinsol_mem->kin_uu);
131  copy(*src_fcur, kinsol_mem->kin_fval);
132 
133  typename VectorMemory<VectorType>::Pointer src(mem);
134  solver.reinit_vector(*src);
135 
136  typename VectorMemory<VectorType>::Pointer dst(mem);
137  solver.reinit_vector(*dst);
138 
139  copy(*src, b);
140 
141  int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src, *dst);
142  copy(x, *dst);
143 
144  *sJpnorm = N_VWL2Norm(b, kinsol_mem->kin_fscale);
145  N_VProd(b, kinsol_mem->kin_fscale, b);
146  N_VProd(b, kinsol_mem->kin_fscale, b);
147  *sFdotJp = N_VDotProd(kinsol_mem->kin_fval, b);
148 
149  return err;
150  }
151  } // namespace
152 
153  template <typename VectorType>
154  KINSOL<VectorType>::KINSOL(const AdditionalData &data,
155  const MPI_Comm mpi_comm)
156  : data(data)
157  , kinsol_mem(nullptr)
158  , solution(nullptr)
159  , u_scale(nullptr)
160  , f_scale(nullptr)
161  , communicator(is_serial_vector<VectorType>::value ?
162  MPI_COMM_SELF :
163  Utilities::MPI::duplicate_communicator(mpi_comm))
164  {
165  set_functions_to_trigger_an_assert();
166  }
167 
168 
169 
170  template <typename VectorType>
171  KINSOL<VectorType>::~KINSOL()
172  {
173  if (kinsol_mem)
174  KINFree(&kinsol_mem);
175 # ifdef DEAL_II_WITH_MPI
176  if (is_serial_vector<VectorType>::value == false)
177  {
178  const int ierr = MPI_Comm_free(&communicator);
179  (void)ierr;
180  AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
181  }
182 # endif
183  }
184 
185 
186 
187  template <typename VectorType>
188  unsigned int
189  KINSOL<VectorType>::solve(VectorType &initial_guess_and_solution)
190  {
191  unsigned int system_size = initial_guess_and_solution.size();
192 
193  // The solution is stored in
194  // solution. Here we take only a
195  // view of it.
196 # ifdef DEAL_II_WITH_MPI
197  if (is_serial_vector<VectorType>::value == false)
198  {
199  const IndexSet is = initial_guess_and_solution.locally_owned_elements();
200  const unsigned int local_system_size = is.n_elements();
201 
202  solution =
203  N_VNew_Parallel(communicator, local_system_size, system_size);
204 
205  u_scale = N_VNew_Parallel(communicator, local_system_size, system_size);
206  N_VConst_Parallel(1.e0, u_scale);
207 
208  f_scale = N_VNew_Parallel(communicator, local_system_size, system_size);
209  N_VConst_Parallel(1.e0, f_scale);
210  }
211  else
212 # endif
213  {
214  Assert(is_serial_vector<VectorType>::value,
215  ExcInternalError(
216  "Trying to use a serial code with a parallel vector."));
217  solution = N_VNew_Serial(system_size);
218  u_scale = N_VNew_Serial(system_size);
219  N_VConst_Serial(1.e0, u_scale);
220  f_scale = N_VNew_Serial(system_size);
221  N_VConst_Serial(1.e0, f_scale);
222  }
223 
224  if (get_solution_scaling)
225  copy(u_scale, get_solution_scaling());
226 
227  if (get_function_scaling)
228  copy(f_scale, get_function_scaling());
229 
230  copy(solution, initial_guess_and_solution);
231 
232  if (kinsol_mem)
233  KINFree(&kinsol_mem);
234 
235  kinsol_mem = KINCreate();
236 
237  int status = KINInit(kinsol_mem, t_kinsol_function<VectorType>, solution);
238  (void)status;
239  AssertKINSOL(status);
240 
241  status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
242  AssertKINSOL(status);
243 
244  status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
245  AssertKINSOL(status);
246 
247  status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
248  AssertKINSOL(status);
249 
250  status = KINSetScaledStepTol(kinsol_mem, data.step_tolerance);
251  AssertKINSOL(status);
252 
253  status = KINSetMaxSetupCalls(kinsol_mem, data.maximum_setup_calls);
254  AssertKINSOL(status);
255 
256  status = KINSetNoInitSetup(kinsol_mem, data.no_init_setup);
257  AssertKINSOL(status);
258 
259  status = KINSetMaxNewtonStep(kinsol_mem, data.maximum_newton_step);
260  AssertKINSOL(status);
261 
262  status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
263  AssertKINSOL(status);
264 
265  status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
266  AssertKINSOL(status);
267 
268  status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
269  AssertKINSOL(status);
270 
271 # if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0)
272  SUNMatrix J = nullptr;
273  SUNLinearSolver LS = nullptr;
274 # endif
275 
276  if (solve_jacobian_system)
277  {
278  auto KIN_mem = static_cast<KINMem>(kinsol_mem);
279  KIN_mem->kin_lsolve = t_kinsol_solve_jacobian<VectorType>;
280  if (setup_jacobian)
281  {
282  KIN_mem->kin_lsetup = t_kinsol_setup_jacobian<VectorType>;
283 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
284  KIN_mem->kin_setupNonNull = true;
285 # endif
286  }
287  }
288  else
289  {
290 # if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0)
291  J = SUNDenseMatrix(system_size, system_size);
292  LS = SUNDenseLinearSolver(u_scale, J);
293  status = KINDlsSetLinearSolver(kinsol_mem, LS, J);
294 # else
295  status = KINDense(kinsol_mem, system_size);
296 # endif
297  AssertKINSOL(status);
298  }
299 
300  if (data.strategy == AdditionalData::newton ||
301  data.strategy == AdditionalData::linesearch)
302  Assert(residual, ExcFunctionNotProvided("residual"));
303 
304  if (data.strategy == AdditionalData::fixed_point ||
305  data.strategy == AdditionalData::picard)
306  Assert(iteration_function, ExcFunctionNotProvided("iteration_function"));
307 
308  // call to KINSol
309  status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
310  AssertKINSOL(status);
311 
312  copy(initial_guess_and_solution, solution);
313 
314  // Free the vectors which are no longer used.
315 # ifdef DEAL_II_WITH_MPI
316  if (is_serial_vector<VectorType>::value == false)
317  {
318  N_VDestroy_Parallel(solution);
319  N_VDestroy_Parallel(u_scale);
320  N_VDestroy_Parallel(f_scale);
321  }
322  else
323 # endif
324  {
325  N_VDestroy_Serial(solution);
326  N_VDestroy_Serial(u_scale);
327  N_VDestroy_Serial(f_scale);
328  }
329 
330  long nniters;
331  status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
332  AssertKINSOL(status);
333 
334 # if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0)
335  SUNMatDestroy(J);
336  SUNLinSolFree(LS);
337 # endif
338  KINFree(&kinsol_mem);
339 
340  return static_cast<unsigned int>(nniters);
341  }
342 
343  template <typename VectorType>
344  void
345  KINSOL<VectorType>::set_functions_to_trigger_an_assert()
346  {
347  reinit_vector = [](VectorType &) {
348  AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
349  };
350  }
351 
352  template class KINSOL<Vector<double>>;
353  template class KINSOL<BlockVector<double>>;
354 
355 # ifdef DEAL_II_WITH_MPI
356 
357 # ifdef DEAL_II_WITH_TRILINOS
358  template class KINSOL<TrilinosWrappers::MPI::Vector>;
359  template class KINSOL<TrilinosWrappers::MPI::BlockVector>;
360 # endif
361 
362 # ifdef DEAL_II_WITH_PETSC
363 # ifndef PETSC_USE_COMPLEX
364  template class KINSOL<PETScWrappers::MPI::Vector>;
365  template class KINSOL<PETScWrappers::MPI::BlockVector>;
366 # endif
367 # endif
368 
369 # endif
370 
371 } // namespace SUNDIALS
372 
373 DEAL_II_NAMESPACE_CLOSE
374 
375 #endif
IndexSet
Definition: index_set.h:74
SUNDIALS::KINSOL
Definition: kinsol.h:198
SUNDIALS
Definition: arkode.h:59
SUNDIALS::KINSOL::AdditionalData
Definition: kinsol.h:204
GrowingVectorMemory
Definition: vector_memory.h:320
IndexSet::n_elements
size_type n_elements() const
Definition: index_set.h:1732
SUNDIALS::KINSOL::set_functions_to_trigger_an_assert
void set_functions_to_trigger_an_assert()
Definition: kinsol.cc:345
is_serial_vector
Definition: vector_type_traits.h:49
VectorMemory::Pointer
Definition: vector_memory.h:192
SUNDIALS::KINSOL::~KINSOL
~KINSOL()
Definition: kinsol.cc:171
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
SUNDIALS::KINSOL::solve
unsigned int solve(VectorType &initial_guess_and_solution)
Definition: kinsol.cc:189
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1407
StandardExceptions::ExcMPI
Definition: exceptions.h:1156
SUNDIALS::KINSOL::KINSOL
KINSOL(const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD)
Definition: kinsol.cc:154
AssertNothrow
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1471
internal
Definition: aligned_vector.h:357
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
Utilities
Definition: cuda.h:31
Generated by   doxygen 1.8.16