File kernel/src/simulationTools/LsodarOSI.hpp

Contents

File kernel/src/simulationTools/LsodarOSI.hpp#

Go to the source code of this file

LsodarOSI solver (from odepack)

class LsodarOSI : public OneStepIntegrator
#include <LsodarOSI.hpp>

LsodarOSI solver (odepack)

Many parameters are required as input/output for LSODAR. See the documentation of this function in externals/odepack/opkdmain.f to have a full description of these parameters. Most of them are read-only parameters (ie can not be set by user).

Except:

  • jt: Jacobian type indicator (1 means a user-supplied full Jacobian, 2 means an internally generated full Jacobian). Default = 2.

  • itol, rtol and atol

ITOL = an indicator for the type of error control.

RTOL = a relative error tolerance parameter, either a scalar or array of length NEQ.

ATOL = an absolute error tolerance parameter, either a scalar or an array of length NEQ. Input only.

Public Types

enum LsodarOSI_ds_workVector_id#

Values:

enumerator FREE#
enumerator WORK_LENGTH#
enum LsodarOSI_interaction_workVector_id#

Values:

enumerator OSNSP_RHS#
enumerator WORK_INTERACTION_LENGTH#
enum LsodarOSI_interaction_workBlockVector_id#

Values:

enumerator xfree#
enumerator BLOCK_WORK_LENGTH#

Public Functions

LsodarOSI()

Default and only constructor.

~LsodarOSI() noexcept = default

destructor

inline const std::vector<int> intData() const
Returns:

int parameters for lsodar

inline int intData(unsigned int i) const
Parameters:

i – index number (starting from 0)

Returns:

int parameter number i

inline void setIntData(unsigned int i, int newValue)

set _intData[i]

Parameters:

  • i – index number (starting from 0)

  • newValue – the new value

inline const std::vector<double> &getRtol() const
Returns:

relative tolerance parameter for lsodar

inline const std::vector<double> &getAtol() const
Returns:

absolute tolerance parameter for lsodar

inline int getMaxNstep() const
Returns:

the maximum number of steps for one call

inline const std::vector<double> &getRwork() const
Returns:

real work vector parameter for lsodar

inline const std::vector<int> &getIwork() const
Returns:

iwork

inline const std::vector<int> &getJroot() const
Returns:

root information

inline void setJT(int newJT)

set Jt value, Jacobian type indicator.

Excerpts from the lsodar documentation. 1 means a user-supplied full (neq by neq) jacobian. 2 means an internally generated (difference quotient) full jacobian (using neq extra calls to f per df/dy value). 4 means a user-supplied banded jacobian. 5 means an internally generated banded jacobian (using ml+mu+1 extra calls to f per df/dy evaluation). if jt = 1 or 4, the user must supply a subroutine jac (the name is arbitrary) as described above under jac. if jt = 2 or 5, a dummy argument can be used.

Parameters:

newJT – new value for the jt parameter.

void setTol(int newItol, std::vector<double> &&newRtol, std::vector<double> &&newAtol)

set itol, rtol and atol (tolerance parameters for lsodar)

Parameters:

  • newItol – itol value

  • newRtol – rtol value

  • newAtol – atol value

void setTol(int newItol, double newRtol, double newAtol)

set itol, rtol and atol (scalar tolerance parameters for lsodar)

Parameters:

  • newItol – itol value

  • newRtol – rtol value

  • newAtol – atol value

void setMaxNstep(int maxNumberSteps)

set the maximum number of steps for one call of Lsodar

Parameters:

maxNumberSteps – the maximum number of steps

void setMinMaxStepSizes(double minStep, double maxStep)

set the minimum and maximum step sizes

Parameters:

  • minStep – minimum step size

  • maxStep – maximum step size

void setMaxOrder(int maxorderNonStiff, int maxorderStiff)

set maximum method order

Parameters:

  • maxorderNonStiff – maximum order for nonstiff methods

  • maxorderStiff – maximum order for stiff methods

void updateData()

update doubleData and iwork memory size, when changes occur in _intData.

void fillXWork(int *size, double *array)

fill xWork with a double

Parameters:

  • size – size of x array

  • array – x array of double

void computeRhs(double t)

compute rhs(t) for all dynamical systems in the set

Parameters:

t – current time of simulation

void computeJacobianRhs(double t, DynamicalSystemsGraph &DSG0)

compute jacobian of the rhs at time t for all dynamical systems in the set

Parameters:
void f(int *sizeOfX, double *time, double *x, double *xdot)#
void g(int *nEq, double *time, double *x, int *ng, double *gOut)#
void jacobianfx(int*, double*, double*, int*, int*, double*, int*)#
virtual void initialize() override

initialization of the integrator

virtual void initializeWorkVectorsForDS(double t, SP::DynamicalSystem ds) override

initialization of the work vectors and matrices (properties) related to one dynamical system on the graph and needed by the osi

Parameters:

  • t – time of initialization

  • ds – the dynamical system

virtual void initializeWorkVectorsForInteraction(Interaction &inter, InteractionProperties &interProp, DynamicalSystemsGraph &DSG) override

initialization of the work vectors and matrices (properties) related to one interaction on the graph and needed by the osi

Parameters:

  • inter – the interaction

  • interProp – the properties on the graph

  • DSG – the dynamical systems graph

inline virtual unsigned int numberOfIndexSets() const override

get the number of index sets required for the simulation

Returns:

unsigned int

virtual void integrate(double &tinit, double &tend, double &tout, int &ioparam) override

integrate the system, between tinit and tend (->iout=true), with possible stop at tout (->iout=false)

Parameters:

  • tinit – initial time

  • tend – end time

  • tout – real end time

  • ioparam – in-out parameter, input: 1 for first call, else 2. Output: 2 if no root was found, else 3.

virtual void updateState(const unsigned int level) override

update the state of the DynamicalSystems attached to this Integrator

Parameters:

level – level of interest for the dynamics

inline virtual void prepareNewtonIteration(double time) override#
virtual void computeFreeOutput(InteractionsGraph::VDescriptor &vertex_descr, OneStepNSProblem *osnsp) override

integrates the Interaction linked to this integrator, without taking non-smooth effects into account

Parameters:

  • vertex_descr – descriptor vertex of the interaction graph

  • osnsp – pointer to OneStepNSProblem

inline virtual SiconosVector &osnsp_rhs(InteractionsGraph::VDescriptor &vertex_inter, InteractionsGraph &indexSet) override

return the workVector corresponding to the right hand side of the OneStepNonsmooth problem

virtual void display() override

print the data to the screen

ACCEPT_STD_VISITORS()#

Public Static Functions

static inline int count_rhs_call()

Return current number of rhs call (for all lsodar-like OSIs!)

Returns:

int

static inline int count_steps()

Return the number of lsodar steps already done (for all lsodar-like OSIs!)

Returns:

int

Public Static Attributes

static int count_NST

Lsodar counter : Number of steps taken for the problem so far.

static int count_NFE

Number of RHS evaluations for the problem so far.

Private Functions

ACCEPT_SERIALIZATION(LsodarOSI)#

Private Members

std::vector<int> _intData = {0, 0, 1, 1, 0, 0, 2}#

neq, ng, itol, itask, istate, iopt, lrw, liw, jt See opkdmain.f and lsodar routine for details on those variables.

std::vector<double> rtol = {RTOL_DEFAULT}#

relative tolerance

std::vector<double> atol = {ATOL_DEFAULT}#

absolute tolerance

std::vector<double> rwork = {}#

real work array

std::vector<int> iwork = {}#

integer work array

std::vector<int> jroot = {}#

integer array used for output of root information

SP::BlockVector _xWork = {nullptr}#

temporary vector to save x values

SP::SiconosVector _xtmp = {nullptr}#

Private Static Attributes

static constexpr double ATOL_DEFAULT = 100 * siconos::internal::MACHINE_PREC#
static constexpr double RTOL_DEFAULT = 10 * siconos::internal::MACHINE_PREC#

Friends

friend struct _NSLEffectOnFreeOutput
Contents