17 #include <deal.II/base/config.h> 19 #include <deal.II/sundials/ida.h> 21 #ifdef DEAL_II_WITH_SUNDIALS 23 # include <deal.II/base/utilities.h> 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> 30 # ifdef DEAL_II_WITH_PETSC 31 # include <deal.II/lac/petsc_block_vector.h> 32 # include <deal.II/lac/petsc_vector.h> 34 # include <deal.II/base/utilities.h> 36 # include <deal.II/sundials/copy.h> 38 # ifdef DEAL_II_SUNDIALS_WITH_IDAS 39 # include <idas/idas_impl.h> 41 # include <ida/ida_impl.h> 47 DEAL_II_NAMESPACE_OPEN
55 template <
typename VectorType>
57 t_dae_residual(realtype tt,
63 IDA<VectorType> &solver = *
static_cast<IDA<VectorType> *
>(user_data);
67 solver.reinit_vector(*src_yy);
70 solver.reinit_vector(*src_yp);
73 solver.reinit_vector(*residual);
78 int err = solver.residual(tt, *src_yy, *src_yp, *residual);
87 template <
typename VectorType>
89 t_dae_lsetup(IDAMem IDA_mem,
101 IDA<VectorType> &solver =
102 *
static_cast<IDA<VectorType> *
>(IDA_mem->ida_user_data);
106 solver.reinit_vector(*src_yy);
109 solver.reinit_vector(*src_yp);
114 int err = solver.setup_jacobian(IDA_mem->ida_tn,
123 template <
typename VectorType>
125 t_dae_solve(IDAMem IDA_mem,
136 IDA<VectorType> &solver =
137 *
static_cast<IDA<VectorType> *
>(IDA_mem->ida_user_data);
141 solver.reinit_vector(*src);
144 solver.reinit_vector(*dst);
148 int err = solver.solve_jacobian_system(*src, *dst);
156 template <
typename VectorType>
166 Utilities::MPI::duplicate_communicator(mpi_comm))
171 template <
typename VectorType>
176 # ifdef DEAL_II_WITH_MPI 188 template <
typename VectorType>
192 unsigned int system_size = solution.size();
196 unsigned int step_number = 0;
204 # ifdef DEAL_II_WITH_MPI 207 const IndexSet is = solution.locally_owned_elements();
208 const size_t local_system_size = is.
n_elements();
210 yy = N_VNew_Parallel(
communicator, local_system_size, system_size);
212 yp = N_VNew_Parallel(
communicator, local_system_size, system_size);
217 N_VNew_Parallel(
communicator, local_system_size, system_size);
224 "Trying to use a serial code with a parallel vector."));
225 yy = N_VNew_Serial(system_size);
226 yp = N_VNew_Serial(system_size);
227 diff_id = N_VNew_Serial(system_size);
240 status = IDASolve(
ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
243 status = IDAGetLastStep(
ida_mem, &h);
247 copy(solution_dot, yp);
250 reset(t, h, solution, solution_dot);
254 output_step(t, solution, solution_dot, step_number);
258 # ifdef DEAL_II_WITH_MPI 261 N_VDestroy_Parallel(yy);
262 N_VDestroy_Parallel(yp);
269 N_VDestroy_Serial(yy);
270 N_VDestroy_Serial(yp);
278 template <
typename VectorType>
281 const double ¤t_time_step,
282 VectorType & solution,
283 VectorType & solution_dot)
285 unsigned int system_size;
297 # ifdef DEAL_II_WITH_MPI 300 N_VDestroy_Parallel(yy);
301 N_VDestroy_Parallel(yp);
308 N_VDestroy_Serial(yy);
309 N_VDestroy_Serial(yp);
317 system_size = solution.size();
318 # ifdef DEAL_II_WITH_MPI 321 const IndexSet is = solution.locally_owned_elements();
322 const size_t local_system_size = is.
n_elements();
324 yy = N_VNew_Parallel(
communicator, local_system_size, system_size);
326 yp = N_VNew_Parallel(
communicator, local_system_size, system_size);
331 N_VNew_Parallel(
communicator, local_system_size, system_size);
336 yy = N_VNew_Serial(system_size);
337 yp = N_VNew_Serial(system_size);
338 diff_id = N_VNew_Serial(system_size);
343 copy(yp, solution_dot);
345 status = IDAInit(
ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
356 status = IDASStolerances(
ida_mem,
362 status = IDASetInitStep(
ida_mem, current_time_step);
365 status = IDASetUserData(
ida_mem, (
void *)
this);
372 VectorType diff_comp_vector(solution);
373 diff_comp_vector = 0.0;
375 for (
auto i = dc.begin(); i != dc.end(); ++i)
376 diff_comp_vector[*i] = 1.0;
378 copy(
diff_id, diff_comp_vector);
397 IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
398 IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
399 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) 400 IDA_mem->ida_setupNonNull =
true;
420 IDACalcIC(
ida_mem, IDA_Y_INIT, current_time + current_time_step);
423 status = IDAGetConsistentIC(
ida_mem, yy, yp);
427 copy(solution_dot, yp);
432 IDACalcIC(
ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
435 status = IDAGetConsistentIC(
ida_mem, yy, yp);
439 copy(solution_dot, yp);
443 template <
typename VectorType>
448 AssertThrow(
false, ExcFunctionNotProvided(
"reinit_vector"));
454 VectorType &) ->
int {
456 AssertThrow(
false, ExcFunctionNotProvided(
"residual"));
463 const double) ->
int {
465 AssertThrow(
false, ExcFunctionNotProvided(
"setup_jacobian"));
471 AssertThrow(
false, ExcFunctionNotProvided(
"solve_jacobian_system"));
478 const unsigned int) {
return; };
481 [](
const double, VectorType &, VectorType &) ->
bool {
return false; };
487 const unsigned int size = v->size();
488 return complete_index_set(size);
495 # ifdef DEAL_II_WITH_MPI 497 # ifdef DEAL_II_WITH_TRILINOS 500 # endif // DEAL_II_WITH_TRILINOS 502 # ifdef DEAL_II_WITH_PETSC 503 # ifndef PETSC_USE_COMPLEX 506 # endif // PETSC_USE_COMPLEX 507 # endif // DEAL_II_WITH_PETSC 509 # endif // DEAL_II_WITH_MPI 513 DEAL_II_NAMESPACE_CLOSE
515 #endif // DEAL_II_WITH_SUNDIALS void set_functions_to_trigger_an_assert()
GrowingVectorMemory< VectorType > mem
std::function< void(VectorType &)> reinit_vector
#define AssertNothrow(cond, exc)
std::function< IndexSet()> differential_components
unsigned int maximum_order
InitialConditionCorrection
unsigned maximum_non_linear_iterations_ic
size_type n_elements() const
void reset(const double &t, const double &h, VectorType &y, VectorType &yp)
InitialConditionCorrection reset_type
#define AssertThrow(cond, exc)
IDA(const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD)
std::function< VectorType &()> get_local_tolerances
#define Assert(cond, exc)
bool ignore_algebraic_terms_for_errors
std::function< int(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
std::function< void(const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> output_step
std::function< int(const VectorType &rhs, VectorType &dst)> solve_jacobian_system
std::function< bool(const double t, VectorType &sol, VectorType &sol_dot)> solver_should_restart
unsigned int maximum_non_linear_iterations
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
double relative_tolerance
double absolute_tolerance
InitialConditionCorrection ic_type
static::ExceptionBase & ExcInternalError()
std::function< int(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian