Reference documentation for deal.II version Git d228448 2017-10-17 21:17:08 +0200
Classes | Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
SUNDIALS::IDA< VectorType > Class Template Reference

#include <deal.II/sundials/ida.h>

Classes

class  AdditionalData
 

Public Member Functions

 IDA (const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD)
 
 ~IDA ()
 
unsigned int solve_dae (VectorType &solution, VectorType &solution_dot)
 
void reset (const double &t, const double &h, VectorType &y, VectorType &yp)
 
 DeclException1 (ExcIDAError, int,<< "One of the SUNDIALS IDA internal functions "<< " returned a negative error code: "<< arg1<< ". Please consult SUNDIALS manual.")
 

Public Attributes

std::function< void(VectorType &)> reinit_vector
 
std::function< int(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
 
std::function< int(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
 
std::function< int(const VectorType &rhs, VectorType &dst)> solve_jacobian_system
 
std::function< void(const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> output_step
 
std::function< bool(const double t, VectorType &sol, VectorType &sol_dot)> solver_should_restart
 
std::function< IndexSet()> differential_components
 
std::function< VectorType &()> get_local_tolerances
 

Private Member Functions

 DeclException1 (ExcFunctionNotProvided, std::string,<< "Please provide an implementation for the function \""<< arg1<< "\"")
 
void set_functions_to_trigger_an_assert ()
 

Private Attributes

AdditionalData data
 
void * ida_mem
 
N_Vector yy
 
N_Vector yp
 
N_Vector abs_tolls
 
N_Vector diff_id
 
MPI_Comm communicator
 
GrowingVectorMemory< VectorType > mem
 

Detailed Description

template<typename VectorType = Vector<double>>
class SUNDIALS::IDA< VectorType >

Interface to SUNDIALS Implicit Differential-Algebraic (IDA) solver.

The class IDA is a wrapper to SUNDIALS Implicit Differential-Algebraic solver which is a general purpose solver for systems of Differential-Algebraic Equations (DAEs).

The user has to provide the implementation of the following std::functions:

Optionally, also the following functions could be rewritten. By default they do nothing, or are not required. If you call the constructor in a way that requires a not-implemented function, an Assertion will be thrown.

To output steps, connect a function to the signal

Citing from the SUNDIALS documentation:

Consider a system of Differential-Algebraic Equations written in the general form

\[ \begin{cases} F(t,y,\dot y) = 0\, , \\ y(t_0) = y_0\, , \\ \dot y (t_0) = \dot y_0\, . \end{cases} \]

where \(y,\dot y\) are vectors in \(\R^n\), \(t\) is often the time (but can also be a parametric quantity), and \(F:\R\times\R^n\times\R^n\rightarrow\R^n\). Such problem is solved using Newton iteration augmented with a line search global strategy. The integration method used in IDA is the variable-order, variable-coefficient BDF (Backward Differentiation Formula), in fixed-leading-coefficient. The method order ranges from 1 to 5, with the BDF of order \(q\) given by the multistep formula

\[ \sum_{i=0}^q \alpha_{n,i}\,y_{n-i}=h_n\,\dot y_n\, , \label{eq:bdf} \]

where \(y_n\) and \(\dot y_n\) are the computed approximations of \(y(t_n)\) and \(\dot y(t_n)\), respectively, and the step size is \(h_n=t_n-t_{n-1}\). The coefficients \(\alpha_{n,i}\) are uniquely determined by the order \(q\), and the history of the step sizes. The application of the BDF method to the DAE system results in a nonlinear algebraic system to be solved at each time step:

\[ G(y_n)\equiv F\left(t_n,y_n,\dfrac{1}{h_n}\sum_{i=0}^q \alpha_{n,i}\,y_{n-i}\right)=0\, . \]

The Newton method leads to a linear system of the form

\[ J[y_{n(m+1)}-y_{n(m)}]=-G(y_{n(m)})\, , \]

where \(y_{n(m)}\) is the \(m\)-th approximation to \(y_n\), and \(J\) is the approximation of the system Jacobian

\[ J=\dfrac{\partial G}{\partial y} = \dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot y}\, , \]

and \(\alpha = \alpha_{n,0}/h_n\). It is worth metioning that the scalar \(\alpha\) changes whenever the step size or method order changes.

To provide a simple example, consider the following harmonic oscillator problem:

\[ \begin{split} u'' & = -k^2 u \\ u (0) & = 0 \\ u'(0) & = k \end{split} \]

We write it in terms of a first order ode:

\[ \begin{matrix} y_0' & -y_1 & = 0 \\ y_1' & + k^2 y_0 & = 0 \end{matrix} \]

That is \(F(y', y, t) = y' + A y = 0 \) where

\[ \begin{matrix} 0 & -1 \\ k^2 &0 \end{matrix} \]

and \(y(0)=(0, k)\), \(y'(0) = (k, 0)\).

The exact solution is \(y_0(t) = \sin(k t)\), \(y_1(t) = y_0'(t) = k \cos(k t)\), \(y_1'(t) = -k^2 \sin(k t)\).

The Jacobian to assemble is the following: \(J = \alpha I + A\).

This is achieved by the following snippet of code:

typedef Vector<double> VectorType;
VectorType y(2);
VectorType y_dot(2);
A(0,1) = -1.0;
A(1,0) = kappa*kappa;
double kappa = 1.0;
IDA time_stepper;
time_stepper.reinit_vector = [&] (VectorType&v)
{
v.reinit(2);
};
time_stepper.residual = [&](const double t,
const VectorType &y,
const VectorType &y_dot,
VectorType &res) ->int
{
res = y_dot;
A.vmult_add(res, y);
return 0;
};
time_stepper.setup_jacobian = [&](const double ,
const VectorType &,
const VectorType &,
const double alpha) ->int
{
J = A;
J(0,0) = alpha;
J(1,1) = alpha;
Jinv.invert(J);
return 0;
};
time_stepper.solve_jacobian_system = [&](const VectorType &src,
VectorType &dst) ->int
{
Jinv.vmult(dst,src);
return 0;
};
y[1] = kappa;
y_dot[0] = kappa;
time_stepper.solve_dae(y,y_dot);
Author
Luca Heltai, Alberto Sartori, 2017.

Definition at line 231 of file ida.h.

Constructor & Destructor Documentation

template<typename VectorType >
SUNDIALS::IDA< VectorType >::IDA ( const AdditionalData data = AdditionalData(),
const MPI_Comm  mpi_comm = MPI_COMM_WORLD 
)

Constructor. It is possible to fine tune the SUNDIALS IDA solver by passing an AdditionalData() object that sets all of the solver parameters.

IDA is a Differential Algebraic solver. As such, it requires initial conditions also for the first order derivatives. If you do not provide consistent initial conditions, (i.e., conditions for which F(y_dot(0), y(0), 0) = 0), you can ask SUNDIALS to compute initial conditions for you by using the ic_type parameter at construction time.

You have three options

  • none: do not try to make initial conditions consistent.
  • use_y_diff: compute the algebraic components of y and differential components of y_dot, given the differential components of y. This option requires that the user specifies differential and algebraic components in the function get_differential_components.
  • use_y_dot: compute all components of y, given y_dot.

By default, this class assumes that all components are differential, and that you want to solve a standard ode. In this case, the initial component type is set to use_y_diff, so that the y_dot at time t=initial_time is computed by solving the nonlinear problem \(F(y_dot, y(t0), t0) = 0\) in the variable y_dot.

Notice that a Newton solver is used for this computation. The Newton solver parameters can be tweaked by acting on ic_alpha and ic_max_iter.

If you reset the solver at some point, you may want to select a different computation for the initial conditions after reset. Say, for example, that you have refined a grid, and after transferring the solution to the new grid, the initial conditions are no longer consistent. Then you can choose how these are made consistent, using the same three options that you used for the initial conditions in reset_type.

The MPI communicator is simply ignored in the serial case.

Parameters
dataIDA configuration data
mpi_commMPI communicator

Definition at line 141 of file ida.cc.

template<typename VectorType >
SUNDIALS::IDA< VectorType >::~IDA ( )

Destructor.

Definition at line 151 of file ida.cc.

Member Function Documentation

template<typename VectorType >
unsigned int SUNDIALS::IDA< VectorType >::solve_dae ( VectorType &  solution,
VectorType &  solution_dot 
)

Integrate differential-algebraic equations. This function returns the final number of computed steps.

Definition at line 163 of file ida.cc.

template<typename VectorType >
void SUNDIALS::IDA< VectorType >::reset ( const double &  t,
const double &  h,
VectorType &  y,
VectorType &  yp 
)

Clear internal memory and start with clean objects. This function is called when the simulation start and when the user returns true to a call to solver_should_restart().

By default solver_should_restart() returns false. If the user needs to implement, for example, local adaptivity in space, he or she may assign a different function to solver_should_restart() that performs all mesh changes, transfers the solution and the solution dot to the new mesh, and returns true.

During reset(), both y and yp are checked for consistency, and according to what was specified as ic_type (if t==initial_time) or reset_type (if t>initial_time), yp, y, or both are modified to obtain a consistent set of initial data.

Parameters
[in]tThe new starting time
[in]hThe new (tentative) starting time step
[in,out]yThe new (tentative) initial solution
[in,out]ypThe new (tentative) initial solution_dot

Definition at line 264 of file ida.cc.

template<typename VectorType = Vector<double>>
SUNDIALS::IDA< VectorType >::DeclException1 ( ExcIDAError  ,
int  ,
<< "One of the SUNDIALS IDA< VectorType > internal functions "<< " returned a negative error code: "<< arg1<< ". Please consult SUNDIALS manual."   
)

Handle IDA exceptions.

template<typename VectorType = Vector<double>>
SUNDIALS::IDA< VectorType >::DeclException1 ( ExcFunctionNotProvided  ,
std::string  ,
<< "Please provide an implementation for the function \""<< arg1<< "\""   
)
private

Throw an exception when a function with the given name is not implemented.

template<typename VectorType >
void SUNDIALS::IDA< VectorType >::set_functions_to_trigger_an_assert ( )
private

This function is executed at construction time to set the std::function above to trigger an assert if they are not implemented.

Definition at line 430 of file ida.cc.

Member Data Documentation

template<typename VectorType = Vector<double>>
std::function<void(VectorType &)> SUNDIALS::IDA< VectorType >::reinit_vector

Reinit vector to have the right size, MPI communicator, etc.

Definition at line 606 of file ida.h.

template<typename VectorType = Vector<double>>
std::function<int(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> SUNDIALS::IDA< VectorType >::residual

Compute residual. Return \(F(t, y, \dot y)\).

This function should return:

  • 0: Success
  • >0: Recoverable error (IDAReinit will be called if this happens, and then last function will be attempted again
  • <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 621 of file ida.h.

template<typename VectorType = Vector<double>>
std::function<int(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> SUNDIALS::IDA< VectorType >::setup_jacobian

Compute Jacobian. This function is called by IDA any time a Jacobian update is required. The user should compute the Jacobian (or update all the variables that allow the application of the Jacobian). This function is called by IDA once, before any call to solve_jacobian_system().

The Jacobian \(J\) should be a (possibly inexact) computation of

\[ J=\dfrac{\partial G}{\partial y} = \dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot y}. \]

If the user uses a matrix based computation of the Jacobian, than this is the right place where an assembly routine shoulde be called to assemble both a matrix and a preconditioner for the Jacobian system. Subsequent calls (possibly more than one) to solve_jacobian_system() can assume that this function has been called at least once.

Notice that no assumption is made by this interface on what the user should do in this function. IDA only assumes that after a call to setup_jacobian() it is possible to call solve_jacobian_system(), to obtain a solution \(x\) to the system \(J x = b\).

This function should return:

  • 0: Success
  • >0: Recoverable error (IDAReinit will be called if this happens, and then last function will be attempted again
  • <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 656 of file ida.h.

template<typename VectorType = Vector<double>>
std::function<int(const VectorType &rhs, VectorType &dst)> SUNDIALS::IDA< VectorType >::solve_jacobian_system

Solve the Jacobian linear system. This function will be called by IDA (possibly several times) after setup_jacobian() has been called at least once. IDA tries to do its best to call setup_jacobian() the minimum amount of times. If convergence can be achieved without updating the Jacobian, then IDA does not call setup_jacobian() again. If, on the contrary, internal IDA convergence tests fail, then IDA calls again setup_jacobian() with updated vectors and coefficents so that successive calls to solve_jacobian_systems() lead to better convergence in the Newton process.

The jacobian \(J\) should be (an approximation of) the system Jacobian

\[ J=\dfrac{\partial G}{\partial y} = \dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot y}. \]

A call to this function should store in dst the result of \(J^{-1}\) applied to src, i.e., J*dst = src. It is the users responsability to set up proper solvers and preconditioners inside this function.

This function should return:

  • 0: Success
  • >0: Recoverable error (IDAReinit will be called if this happens, and then last function will be attempted again
  • <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 686 of file ida.h.

template<typename VectorType = Vector<double>>
std::function<void (const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> SUNDIALS::IDA< VectorType >::output_step

Process solution. This function is called by IDA at fixed time steps, every output_period seconds, and it is passed a polynomial interpolation of the solution and of its time derivative, computed using the current BDF order and the (internally stored) previously computed solution steps.

Notice that it is well possible that internally IDA computes a time step which is much larger than the output_period step, and therefore calls this function consecutively several times by simply performing all intermediate interpolations. There is no relationship between how many times this function is called and how many time steps have actually been computed.

Definition at line 705 of file ida.h.

template<typename VectorType = Vector<double>>
std::function<bool (const double t, VectorType &sol, VectorType &sol_dot)> SUNDIALS::IDA< VectorType >::solver_should_restart

Evaluate wether the solver should be restarted (for example because the number of degrees of freedom has changed).

This function is supposed to perform all operations that are necessary in sol and sol_dot to make sure that the resulting vectors are consistent, and of the correct final size.

For example, one may decide that a local refinement is necessary at time t. This function should then return true, and change the dimension of both sol and sol_dot to reflect the new dimension. Since IDA does not know about the new dimension, an internal reset is necessary.

The default implementation simply returns false, i.e., no restart is performed during the evolution.

Definition at line 725 of file ida.h.

template<typename VectorType = Vector<double>>
std::function<IndexSet()> SUNDIALS::IDA< VectorType >::differential_components

Return an index set containing the differential components. Implementation of this function is optional. The default is to return a complete index set. If your equation is also algebraic (i.e., it contains algebraic constraints, or lagrange multipliers), you should overwrite this function in order to return only the differential components of your system.

Definition at line 735 of file ida.h.

template<typename VectorType = Vector<double>>
std::function<VectorType&()> SUNDIALS::IDA< VectorType >::get_local_tolerances

Return a vector whose components are the weights used by IDA to compute the vector norm. The implementation of this function is optional. If the user does not provide an implementation, the weights are assumed to be all ones.

Definition at line 743 of file ida.h.

template<typename VectorType = Vector<double>>
AdditionalData SUNDIALS::IDA< VectorType >::data
private

IDA configuration data.

Definition at line 771 of file ida.h.

template<typename VectorType = Vector<double>>
void* SUNDIALS::IDA< VectorType >::ida_mem
private

IDA memory object.

Definition at line 776 of file ida.h.

template<typename VectorType = Vector<double>>
N_Vector SUNDIALS::IDA< VectorType >::yy
private

IDA solution vector.

Definition at line 781 of file ida.h.

template<typename VectorType = Vector<double>>
N_Vector SUNDIALS::IDA< VectorType >::yp
private

IDA solution derivative vector.

Definition at line 786 of file ida.h.

template<typename VectorType = Vector<double>>
N_Vector SUNDIALS::IDA< VectorType >::abs_tolls
private

IDA absolute tolerances vector.

Definition at line 791 of file ida.h.

template<typename VectorType = Vector<double>>
N_Vector SUNDIALS::IDA< VectorType >::diff_id
private

IDA differential components vector.

Definition at line 796 of file ida.h.

template<typename VectorType = Vector<double>>
MPI_Comm SUNDIALS::IDA< VectorType >::communicator
private

MPI communicator. SUNDIALS solver runs happily in parallel. Note that if the library is compiled without MPI support, MPI_Comm is typedefed as int.

Definition at line 803 of file ida.h.

template<typename VectorType = Vector<double>>
GrowingVectorMemory<VectorType> SUNDIALS::IDA< VectorType >::mem
private

Memory pool of vectors.

Definition at line 808 of file ida.h.


The documentation for this class was generated from the following files: