Reference documentation for deal.II version 9.0.0
ida.cc
1 //-----------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 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 at
12 // the top level of the deal.II distribution.
13 //
14 //-----------------------------------------------------------
15 
16 
17 #include <deal.II/sundials/ida.h>
18 #include <deal.II/base/config.h>
19 
20 #ifdef DEAL_II_WITH_SUNDIALS
21 
22 #include <deal.II/base/utilities.h>
23 #include <deal.II/lac/block_vector.h>
24 #ifdef DEAL_II_WITH_TRILINOS
25 #include <deal.II/lac/trilinos_parallel_block_vector.h>
26 #include <deal.II/lac/trilinos_vector.h>
27 #endif
28 #ifdef DEAL_II_WITH_PETSC
29 #include <deal.II/lac/petsc_parallel_block_vector.h>
30 #include <deal.II/lac/petsc_parallel_vector.h>
31 #endif
32 #include <deal.II/base/utilities.h>
33 #include <deal.II/sundials/copy.h>
34 
35 #ifdef DEAL_II_SUNDIALS_WITH_IDAS
36 #include <idas/idas_impl.h>
37 #else
38 #include <ida/ida_impl.h>
39 #endif
40 
41 #include <iostream>
42 #include <iomanip>
43 
44 DEAL_II_NAMESPACE_OPEN
45 
46 namespace SUNDIALS
47 {
48  using namespace internal;
49 
50  namespace
51  {
52  template<typename VectorType>
53  int t_dae_residual(realtype tt, N_Vector yy, N_Vector yp,
54  N_Vector rr, void *user_data)
55  {
56  IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
57  GrowingVectorMemory<VectorType> mem;
58 
59  typename VectorMemory<VectorType>::Pointer src_yy(mem);
60  solver.reinit_vector(*src_yy);
61 
62  typename VectorMemory<VectorType>::Pointer src_yp(mem);
63  solver.reinit_vector(*src_yp);
64 
65  typename VectorMemory<VectorType>::Pointer residual(mem);
66  solver.reinit_vector(*residual);
67 
68  copy(*src_yy, yy);
69  copy(*src_yp, yp);
70 
71  int err = solver.residual(tt, *src_yy, *src_yp, *residual);
72 
73  copy(rr, *residual);
74 
75  return err;
76  }
77 
78 
79 
80  template<typename VectorType>
81  int t_dae_lsetup(IDAMem IDA_mem,
82  N_Vector yy,
83  N_Vector yp,
84  N_Vector resp,
85  N_Vector tmp1,
86  N_Vector tmp2,
87  N_Vector tmp3)
88  {
89  (void) tmp1;
90  (void) tmp2;
91  (void) tmp3;
92  (void) resp;
93  IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
94  GrowingVectorMemory<VectorType> mem;
95 
96  typename VectorMemory<VectorType>::Pointer src_yy(mem);
97  solver.reinit_vector(*src_yy);
98 
99  typename VectorMemory<VectorType>::Pointer src_yp(mem);
100  solver.reinit_vector(*src_yp);
101 
102  copy(*src_yy, yy);
103  copy(*src_yp, yp);
104 
105  int err = solver.setup_jacobian(IDA_mem->ida_tn,
106  *src_yy,
107  *src_yp,
108  IDA_mem->ida_cj);
109 
110  return err;
111  }
112 
113 
114  template<typename VectorType>
115  int t_dae_solve(IDAMem IDA_mem,
116  N_Vector b,
117  N_Vector weight,
118  N_Vector yy,
119  N_Vector yp,
120  N_Vector resp)
121  {
122  (void) weight;
123  (void) yy;
124  (void) yp;
125  (void) resp;
126  IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
127  GrowingVectorMemory<VectorType> mem;
128 
129  typename VectorMemory<VectorType>::Pointer src(mem);
130  solver.reinit_vector(*src);
131 
132  typename VectorMemory<VectorType>::Pointer dst(mem);
133  solver.reinit_vector(*dst);
134 
135  copy(*src, b);
136 
137  int err = solver.solve_jacobian_system(*src,*dst);
138  copy(b, *dst);
139 
140  return err;
141  }
142 
143  }
144 
145  template <typename VectorType>
146  IDA<VectorType>::IDA(const AdditionalData &data,
147  const MPI_Comm mpi_comm) :
148  data(data),
149  ida_mem(nullptr),
150  yy(nullptr),
151  yp(nullptr),
152  abs_tolls(nullptr),
153  diff_id(nullptr),
154  communicator(is_serial_vector<VectorType>::value ?
155  MPI_COMM_SELF :
156  Utilities::MPI::duplicate_communicator(mpi_comm))
157  {
158  set_functions_to_trigger_an_assert();
159  }
160 
161  template <typename VectorType>
162  IDA<VectorType>::~IDA()
163  {
164  if (ida_mem)
165  IDAFree(&ida_mem);
166 #ifdef DEAL_II_WITH_MPI
167  if (is_serial_vector<VectorType>::value == false)
168  {
169  const int ierr = MPI_Comm_free(&communicator);
170  (void)ierr;
171  AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
172  }
173 #endif
174  }
175 
176 
177 
178  template <typename VectorType>
179  unsigned int IDA<VectorType>::solve_dae(VectorType &solution,
180  VectorType &solution_dot)
181  {
182 
183  unsigned int system_size = solution.size();
184 
185  double t = data.initial_time;
186  double h = data.initial_step_size;
187  unsigned int step_number = 0;
188 
189  int status;
190  (void)status;
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 = solution.locally_owned_elements();
199  const size_t local_system_size = is.n_elements();
200 
201  yy = N_VNew_Parallel(communicator,
202  local_system_size,
203  system_size);
204 
205  yp = N_VNew_Parallel(communicator,
206  local_system_size,
207  system_size);
208 
209  diff_id = N_VNew_Parallel(communicator,
210  local_system_size,
211  system_size);
212 
213  abs_tolls = N_VNew_Parallel(communicator,
214  local_system_size,
215  system_size);
216  }
217  else
218 #endif
219  {
220  Assert(is_serial_vector<VectorType>::value,
221  ExcInternalError("Trying to use a serial code with a parallel vector."));
222  yy = N_VNew_Serial(system_size);
223  yp = N_VNew_Serial(system_size);
224  diff_id = N_VNew_Serial(system_size);
225  abs_tolls = N_VNew_Serial(system_size);
226  }
227  reset(data.initial_time,
228  data.initial_step_size,
229  solution,
230  solution_dot);
231 
232  double next_time = data.initial_time;
233 
234  output_step( 0, solution, solution_dot, 0);
235 
236  while (t<data.final_time)
237  {
238 
239  next_time += data.output_period;
240 
241  status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
242  AssertIDA(status);
243 
244  status = IDAGetLastStep(ida_mem, &h);
245  AssertIDA(status);
246 
247  copy(solution, yy);
248  copy(solution_dot, yp);
249 
250  while (solver_should_restart(t, solution, solution_dot))
251  reset(t, h, solution, solution_dot);
252 
253  step_number++;
254 
255  output_step(t, solution, solution_dot, step_number);
256  }
257 
258  // Free the vectors which are no longer used.
259 #ifdef DEAL_II_WITH_MPI
260  if (is_serial_vector<VectorType>::value == false)
261  {
262  N_VDestroy_Parallel(yy);
263  N_VDestroy_Parallel(yp);
264  N_VDestroy_Parallel(abs_tolls);
265  N_VDestroy_Parallel(diff_id);
266  }
267  else
268 #endif
269  {
270  N_VDestroy_Serial(yy);
271  N_VDestroy_Serial(yp);
272  N_VDestroy_Serial(abs_tolls);
273  N_VDestroy_Serial(diff_id);
274  }
275 
276  return step_number;
277  }
278 
279  template <typename VectorType>
280  void IDA<VectorType>::reset(const double &current_time,
281  const double &current_time_step,
282  VectorType &solution,
283  VectorType &solution_dot)
284  {
285 
286  unsigned int system_size;
287  bool first_step = (current_time == data.initial_time);
288 
289  if (ida_mem)
290  IDAFree(&ida_mem);
291 
292  ida_mem = IDACreate();
293 
294 
295  // Free the vectors which are no longer used.
296  if (yy)
297  {
298 #ifdef DEAL_II_WITH_MPI
299  if (is_serial_vector<VectorType>::value == false)
300  {
301  N_VDestroy_Parallel(yy);
302  N_VDestroy_Parallel(yp);
303  N_VDestroy_Parallel(abs_tolls);
304  N_VDestroy_Parallel(diff_id);
305  }
306  else
307 #endif
308  {
309  N_VDestroy_Serial(yy);
310  N_VDestroy_Serial(yp);
311  N_VDestroy_Serial(abs_tolls);
312  N_VDestroy_Serial(diff_id);
313  }
314  }
315 
316  int status;
317  (void)status;
318  system_size = solution.size();
319 #ifdef DEAL_II_WITH_MPI
320  if (is_serial_vector<VectorType>::value == false)
321  {
322  const IndexSet is = solution.locally_owned_elements();
323  const size_t local_system_size = is.n_elements();
324 
325  yy = N_VNew_Parallel(communicator,
326  local_system_size,
327  system_size);
328 
329  yp = N_VNew_Parallel(communicator,
330  local_system_size,
331  system_size);
332 
333  diff_id = N_VNew_Parallel(communicator,
334  local_system_size,
335  system_size);
336 
337  abs_tolls = N_VNew_Parallel(communicator,
338  local_system_size,
339  system_size);
340 
341  }
342  else
343 #endif
344  {
345  yy = N_VNew_Serial(system_size);
346  yp = N_VNew_Serial(system_size);
347  diff_id = N_VNew_Serial(system_size);
348  abs_tolls = N_VNew_Serial(system_size);
349  }
350 
351  copy(yy, solution);
352  copy(yp, solution_dot);
353 
354  status = IDAInit(ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
355  AssertIDA(status);
356 
357  if (get_local_tolerances)
358  {
359  copy(abs_tolls, get_local_tolerances());
360  status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tolls);
361  AssertIDA(status);
362  }
363  else
364  {
365  status = IDASStolerances(ida_mem, data.relative_tolerance, data.absolute_tolerance);
366  AssertIDA(status);
367  }
368 
369  status = IDASetInitStep(ida_mem, current_time_step);
370  AssertIDA(status);
371 
372  status = IDASetUserData(ida_mem, (void *) this);
373  AssertIDA(status);
374 
375  if (data.ic_type == AdditionalData::use_y_diff ||
376  data.reset_type == AdditionalData::use_y_diff ||
377  data.ignore_algebraic_terms_for_errors)
378  {
379  VectorType diff_comp_vector(solution);
380  diff_comp_vector = 0.0;
381  auto dc = differential_components();
382  for (auto i = dc.begin(); i != dc.end(); ++i)
383  diff_comp_vector[*i] = 1.0;
384 
385  copy(diff_id, diff_comp_vector);
386  status = IDASetId(ida_mem, diff_id);
387  AssertIDA(status);
388  }
389 
390  status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
391  AssertIDA(status);
392 
393 // status = IDASetMaxNumSteps(ida_mem, max_steps);
394  status = IDASetStopTime(ida_mem, data.final_time);
395  AssertIDA(status);
396 
397  status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
398  AssertIDA(status);
399 
400  // Initialize solver
401  IDAMem IDA_mem;
402  IDA_mem = (IDAMem) ida_mem;
403 
404  IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
405  IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
406 #if DEAL_II_SUNDIALS_VERSION_LT(3,0,0)
407  IDA_mem->ida_setupNonNull = true;
408 #endif
409 
410  status = IDASetMaxOrd(ida_mem, data.maximum_order);
411  AssertIDA(status);
412 
413  typename AdditionalData::InitialConditionCorrection type;
414  if (first_step)
415  type = data.ic_type;
416  else
417  type = data.reset_type;
418 
419  status = IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
420  AssertIDA(status);
421 
422  if (type == AdditionalData::use_y_dot)
423  {
424  // (re)initialization of the vectors
425  status = IDACalcIC(ida_mem, IDA_Y_INIT, current_time+current_time_step);
426  AssertIDA(status);
427 
428  status = IDAGetConsistentIC(ida_mem, yy, yp);
429  AssertIDA(status);
430 
431  copy(solution, yy);
432  copy(solution_dot, yp);
433  }
434  else if (type == AdditionalData::use_y_diff)
435  {
436  status = IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time+current_time_step);
437  AssertIDA(status);
438 
439  status = IDAGetConsistentIC(ida_mem, yy, yp);
440  AssertIDA(status);
441 
442  copy(solution, yy);
443  copy(solution_dot, yp);
444  }
445  }
446 
447  template<typename VectorType>
448  void IDA<VectorType>::set_functions_to_trigger_an_assert()
449  {
450 
451  reinit_vector = [](VectorType &)
452  {
453  AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
454  };
455 
456  residual = [](const double,
457  const VectorType &,
458  const VectorType &,
459  VectorType &) ->int
460  {
461  int ret=0;
462  AssertThrow(false, ExcFunctionNotProvided("residual"));
463  return ret;
464  };
465 
466  setup_jacobian = [](const double,
467  const VectorType &,
468  const VectorType &,
469  const double) ->int
470  {
471  int ret=0;
472  AssertThrow(false, ExcFunctionNotProvided("setup_jacobian"));
473  return ret;
474  };
475 
476  solve_jacobian_system = [](const VectorType &,
477  VectorType &) ->int
478  {
479  int ret=0;
480  AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
481  return ret;
482  };
483 
484  output_step = [](const double,
485  const VectorType &,
486  const VectorType &,
487  const unsigned int)
488  {
489  return;
490  };
491 
492  solver_should_restart = [](const double,
493  VectorType &,
494  VectorType &) ->bool
495  {
496  return false;
497  };
498 
499  differential_components = [&]() ->IndexSet
500  {
501  GrowingVectorMemory<VectorType> mem;
502  typename VectorMemory<VectorType>::Pointer v(mem);
503  reinit_vector(*v);
504  const unsigned int size = v->size();
505  return complete_index_set(size);
506  };
507  }
508 
509  template class IDA<Vector<double> >;
510  template class IDA<BlockVector<double> >;
511 
512 #ifdef DEAL_II_WITH_MPI
513 
514 #ifdef DEAL_II_WITH_TRILINOS
515  template class IDA<TrilinosWrappers::MPI::Vector>;
516  template class IDA<TrilinosWrappers::MPI::BlockVector>;
517 #endif // DEAL_II_WITH_TRILINOS
518 
519 #ifdef DEAL_II_WITH_PETSC
520 #ifndef PETSC_USE_COMPLEX
521  template class IDA<PETScWrappers::MPI::Vector>;
522  template class IDA<PETScWrappers::MPI::BlockVector>;
523 #endif // PETSC_USE_COMPLEX
524 #endif // DEAL_II_WITH_PETSC
525 
526 #endif //DEAL_II_WITH_MPI
527 
528 }
529 
530 DEAL_II_NAMESPACE_CLOSE
531 
532 #endif // DEAL_II_WITH_SUNDIALS
SUNDIALS::IDA::set_functions_to_trigger_an_assert
void set_functions_to_trigger_an_assert()
Definition: ida.cc:448
AssertNothrow
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1183
SUNDIALS::IDA::reset
void reset(const double &t, const double &h, VectorType &y, VectorType &yp)
Definition: ida.cc:280
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1221
is_serial_vector
Definition: vector_type_traits.h:47
SUNDIALS::IDA::IDA
IDA(const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD)
Definition: ida.cc:146
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1142
Utilities
Definition: cuda.h:31
SUNDIALS::IDA::solve_dae
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
Definition: ida.cc:179
IndexSet::n_elements
size_type n_elements() const
Definition: index_set.h:1717
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.14