File kernel/src/simulationTools/EulerMoreauOSI.hpp

Go to the source code of this file

Variables

const unsigned int EULERMOREAUSTEPSINMEMORY = 1
class EulerMoreauOSI : public OneStepIntegrator
#include <EulerMoreauOSI.hpp>

One Step time Integrator for First Order Dynamical Systems.

This integrator is the work horse of the eventcapturing time stepping schemes for first order systems. It is mainly based on some extensions of the Backward Euler and \(\theta-\gamma\) schemes proposed in the pionnering work of J.J. Moreau for the sweeping process

J.J. Moreau. Evolution problem associated with a moving convex set in a Hilbert space. Journal of Differential Equations, 26, pp 347374, 1977.

Variants are now used to integrate LCS, Relay systems, Higher order sweeping process see for instance

Consistency of a time-stepping method for a class of piecewise linear networks

M.K. Camlibel, W.P.M.H. Heemels, and J.M. Schumacher IEEE Transactions on Circuits and Systems I, 2002, 49(3):349357

Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics

V Acary, B Brogliato Springer Verlag 2008

Convergence of time-stepping schemes for passive and extended linear complementarity systems L. Han, A. Tiwari, M.K. Camlibel, and J.-S. Pang SIAM Journal on Numerical Analysis 2009, 47(5):3768-3796

On preserving dissipativity properties of linear complementarity dynamical systems with the &theta-method

Greenhalgh Scott, Acary Vincent, Brogliato Bernard Numer. Math., , 2013.

Main timeintegration schemes are based on the following \(\theta-\gamma\) scheme

\begin{cases} \label{eq:toto1} M x_{k+1} = M x_{k} +h\theta f(x_{k+1},t_{k+1})+h(1-\theta) f(x_k,t_k) + h \gamma r(t_{k+1}) + h(1-\gamma)r(t_k) \\[2mm] y_{k+1} = h(t_{k+1},x_{k+1},\lambda _{k+1}) \\[2mm] r_{k+1} = g(x_{k+1},\lambda_{k+1},t_{k+1})\\[2mm] \mbox{nslaw} ( y_{k+1} , \lambda_{k+1}) \end{cases}

where \(\theta = [0,1]\) and \(\gamma \in [0,1]\). As in Acary & Brogliato 2008, we call the previous problem the ``onestep nonsmooth problem’‘.

Another variant can also be used (FullThetaGamma scheme)

\begin{cases} M x_{k+1} = M x_{k} +h f(x_{k+\theta},t_{k+1}) + h r(t_{k+\gamma}) \\[2mm] y_{k+\gamma} = h(t_{k+\gamma},x_{k+\gamma},\lambda _{k+\gamma}) \\[2mm] r_{k+\gamma} = g(x_{k+\gamma},\lambda_{k+\gamma},t_{k+\gamma})\\[2mm] \mbox{nslaw} ( y_{k+\gamma} , \lambda_{k+\gamma}) \end{cases}

EulerMoreauOSI class is used to define some time-integrators methods for a list of first order dynamical systems. A EulerMoreauOSI instance is defined by the value of theta and possibly gamma and the list of concerned dynamical systems.

Each DynamicalSystem is associated to a SiconosMatrix, named “W”, which is the “iteration” matrix. W matrices are initialized and computed in initializeIterationMatrixW and computeW. Depending on the DS type, they may depend on time t and DS state x.

For first order systems, the implementation uses _r for storing the the input due to the nonsmooth law. This EulerMoreauOSI scheme assumes that the relative degree is zero or one and one level for _r is sufficient

Main functions:

  • computeFreeState(): computes xfree (or vfree), dynamical systems state without taking non-smooth part into account
  • updateState(): computes x (q,v), the complete dynamical systems states.

See User’s guide, for details.

Public Types

enum EulerMoreauOSI_ds_workVector_id

Values:

RESIDU
RESIDU_FREE
FREE
X_PARTIAL_NS_FOR_RELATION
DELTA_X_FOR_RELATION
LOCAL_BUFFER
WORK_LENGTH
enum EulerMoreauOSI_interaction_workBlockVector_id

Values:

XFREE
X_PARTIAL_NS
DELTA_X
BLOCK_WORK_LENGTH
enum EulerMoreauOSI_interaction_workMat_id

Values:

MAT_KHAT
MAT_KTILDE
MAT_WORK_LENGTH
enum EulerMoreauOSI_interaction_workVector_id

Values:

OSNSP_RHS
VEC_X
VEC_Z
H_ALPHA
VEC_RESIDU_Y
G_ALPHA
VEC_RESIDU_R
WORK_INTERACTION_LENGTH

Public Functions

EulerMoreauOSI(double theta)

constructor from theta value only

Parameters
  • theta: value for all DS.

EulerMoreauOSI(double theta, double gamma)

constructor from theta value only

Parameters
  • theta: value for all linked DS.
  • gamma: value for all linked DS.

virtual ~EulerMoreauOSI()

destructor

ACCEPT_STD_VISITORS()

visitors hook

virtual void computeFreeOutput(InteractionsGraph::VDescriptor &vertex_inter, OneStepNSProblem *osnsp)

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

Parameters

virtual void computeFreeState()

Perform the integration of the dynamical systems linked to this integrator without taking into account the nonsmooth input r.

void computeKhat(Interaction &inter, SiconosMatrix &m, VectorOfSMatrices &workM, double h) const
double computeResidu()

Computes the residuFree and residu of all the DynamicalSystems.

Return
the maximum of the 2-norm over all the residu

double computeResiduInput(double time, SP::InteractionsGraph indexSet)

compute the residu of the input of the relation (R or p) This computation depends on the type of OSI

Parameters
  • time: time of computation
  • indexSet: the index set of the interaction that are concerned

double computeResiduOutput(double time, SP::InteractionsGraph indexSet)

compute the residu of the output of the relation (y) This computation depends on the type of OSI

Parameters
  • time: time of computation
  • indexSet: the index set of the interaction that are concerned

void computeW(double time, DynamicalSystem &ds, DynamicalSystemsGraph::VDescriptor &dsv, SiconosMatrix &W)

compute W EulerMoreauOSI matrix at time t

Parameters
  • time: the current time
  • ds: the DynamicalSystem
  • dsv: a descriptor of the ds on the graph (redundant to avoid invocation)
  • W: the matrix to compute

void computeWBoundaryConditions(SP::DynamicalSystem ds)

compute WBoundaryConditionsMap[ds] EulerMoreauOSI matrix at time t

Parameters

void display()

Displays the data of the EulerMoreauOSI’s integrator.

double gamma()

get gamma

Return
a double

const SimpleMatrix getW(SP::DynamicalSystem ds = SP::DynamicalSystem())

get the value of W corresponding to DynamicalSystem ds

Return
SimpleMatrix
Parameters
  • ds: a pointer to DynamicalSystem, optional, default = NULL. get W[0] in that case

const SimpleMatrix getWBoundaryConditions(SP::DynamicalSystem ds = SP::DynamicalSystem())

get the value of WBoundaryConditions corresponding to DynamicalSystem ds

Return
SimpleMatrix
Parameters
  • ds: a pointer to DynamicalSystem, optional, default = NULL. get WBoundaryConditions[0] in that case

void initializeIterationMatrixW(double time, SP::DynamicalSystem ds)

initialize iteration matrix W EulerMoreauOSI matrix at time t

Parameters

void initializeIterationMatrixWBoundaryConditions(SP::DynamicalSystem ds)

initialize iteration matrix WBoundaryConditionsMap[ds] EulerMoreauOSI

Parameters

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

initialization of the EulerMoreauOSI integrator; for linear time invariant systems, we compute time invariant operator (example : W)

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

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

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

void integrate(double &tinit, double &tend, double &tout, int &useless)

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

Parameters

unsigned int numberOfIndexSets() const

get the number of index sets required for the simulation

Return
unsigned int

void prepareNewtonIteration(double time)

computes all the W matrices

Parameters
  • time: current time

void setGamma(double newGamma)

set the value of gamma

Parameters
  • newGamma: a double

void setTheta(double newTheta)

set the value of theta

Parameters
  • newTheta: a double

void setUseGamma(bool b)

set the boolean to indicate that we use gamma

Parameters
  • b: true if gamma has to be used, false otherwise

void setUseGammaForRelation(bool newUseGammaForRelation)

set the boolean to indicate that we use gamma for the relation

Parameters
  • newUseGammaForRelation: a bool

double theta()

get theta

Return
a double

void updateInput(double time)

update the input of the Interaction attached to this Integrator

void updateInput(double time, unsigned int level)

update the input of the Interaction attached to this Integrator

Parameters
  • time: current time
  • level: level of interest for the dynamics

void updateOutput(double time)

update the output of the Interaction attached to this Integrator

void updateOutput(double time, unsigned int level)

update the output of the Interaction attached to this Integrator

Parameters
  • time: current time
  • level: level of interest for the dynamics

virtual void updateState(const unsigned int level)

updates the state of the Dynamical Systems

Parameters
  • level: the level of interest for the dynamics: not used at the time

bool useGamma()

get bool useGamma

Return
a bool

bool useGammaForRelation()

get bool gammaForRelation for the relation

Return
a

SP::SimpleMatrix W(SP::DynamicalSystem ds)

get W corresponding to DynamicalSystem ds

Return
pointer to a SiconosMatrix
Parameters

SP::SiconosMatrix WBoundaryConditions(SP::DynamicalSystem ds)

get WBoundaryConditions corresponding to DynamicalSystem ds

Return
pointer to a SiconosMatrix
Parameters
  • ds: a pointer to DynamicalSystem, optional, default = NULL. get WBoundaryConditions[0] in that case

Protected Functions

EulerMoreauOSI()

Default constructor.

ACCEPT_SERIALIZATION(EulerMoreauOSI)

serialization hooks

Protected Attributes

double _gamma

A gamma parameter for the integration scheme to each DynamicalSystem of the OSI This parameter is used to apply a theta-method to the input $r$.

double _theta

Stl map that associates a theta parameter for the integration scheme to each DynamicalSystem of the OSI.

bool _useGamma

a boolean to know if the parameter must be used or not

bool _useGammaForRelation

a boolean to know if the parameter must be used or not

Friends

friend _NSLEffectOnFreeOutput