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