File kernel/src/simulationTools/EulerMoreauOSI.hpp

File kernel/src/simulationTools/EulerMoreauOSI.hpp#

class EulerMoreauOSI : public OneStepIntegrator

#include <EulerMoreauOSI.hpp>

One Step time Integrator for First Order Dynamical Systems.

This integrator is the work horse of the event—capturing 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 347—374, 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):349—357

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 time—integration schemes are based on the following $ \theta-\gamma $ scheme

\[\begin{split} \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} \end{split}\]

where $ \theta = [0,1] $ and $ \gamma \in [0,1] $. As in Acary & Brogliato 2008, we call the previous problem the `one–step nonsmooth problem’.

Another variant can also be used (FullThetaGamma scheme)

\[\begin{split} \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} \end{split}\]

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:

enumerator RESIDU#

enumerator RESIDU_FREE#

enumerator FREE#

enumerator X_PARTIAL_NS_FOR_RELATION#

enumerator DELTA_X_FOR_RELATION#

enumerator LOCAL_BUFFER#

enumerator WORK_LENGTH#

enum EulerMoreauOSI_interaction_workVector_id#

Values:

enumerator OSNSP_RHS#

enumerator VEC_X#

enumerator H_ALPHA#

enumerator VEC_RESIDU_Y#

enumerator VEC_RESIDU_R#

enumerator YOLD#

enumerator LAMBDAOLD#

enumerator WORK_INTERACTION_LENGTH#

enum EulerMoreauOSI_interaction_workBlockVector_id#

Values:

enumerator XFREE#

enumerator X_PARTIAL_NS#

enumerator DELTA_X#

enumerator G_ALPHA#

enumerator BLOCK_WORK_LENGTH#

enum EulerMoreauOSI_interaction_workMat_id#

Values:

enumerator MAT_KHAT#

enumerator MAT_KTILDE#

enumerator MAT_WORK_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.

inline virtual ~EulerMoreauOSI(): destructor

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

get the value of W corresponding to DynamicalSystem ds

Parameters:: ds – a pointer to DynamicalSystem, optional, default = nullptr. get W[0] in that case
Returns:: SimpleMatrix

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

get W corresponding to DynamicalSystem ds

Parameters:: ds – a pointer to DynamicalSystem
Returns:: pointer to a SiconosMatrix

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

get the value of WBoundaryConditions corresponding to DynamicalSystem ds

Parameters:: ds – a pointer to DynamicalSystem, optional, default = nullptr. get WBoundaryConditions[0] in that case
Returns:: SimpleMatrix

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

get WBoundaryConditions corresponding to DynamicalSystem ds

Parameters:: ds – a pointer to DynamicalSystem, optional, default = nullptr. get WBoundaryConditions[0] in that case
Returns:: pointer to a SiconosMatrix

inline double theta()

get theta

Returns:: a double

inline void setTheta(double newTheta)

set the value of theta

Parameters:: newTheta – a double

inline double gamma()

get gamma

Returns:: a double

inline void setGamma(double newGamma)

set the value of gamma

Parameters:: newGamma – a double

inline bool useGamma()

get bool useGamma

Returns:: a bool

inline void setUseGamma(bool b)

set the boolean to indicate that we use gamma

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

inline bool useGammaForRelation()

get bool gammaForRelation for the relation

Returns:: a

inline void setUseGammaForRelation(bool newUseGammaForRelation)

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

Parameters:: newUseGammaForRelation – a bool

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

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

initialize iteration matrix W EulerMoreauOSI matrix at time t

Parameters:

time – the time (double)
ds – a pointer to DynamicalSystem

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 computeKhat(Interaction &inter, SiconosMatrix &m, VectorOfSMatrices &workM, double h) const#

void computeWBoundaryConditions(SP::DynamicalSystem ds)

compute WBoundaryConditionsMap[ds] EulerMoreauOSI matrix at time t

Parameters:: ds – a pointer to DynamicalSystem

void initializeIterationMatrixWBoundaryConditions(SP::DynamicalSystem ds)

initialize iteration matrix WBoundaryConditionsMap[ds] EulerMoreauOSI

Parameters:: ds – a pointer to DynamicalSystem

virtual double computeResidu() override

Computes the residuFree and residu of all the DynamicalSystems.

Returns:: the maximum of the 2-norm over all the residu

virtual void computeFreeState() override: Perform the integration of the dynamical systems linked to this integrator without taking into account the nonsmooth input r.

virtual void updateOutput(double time) override: update the output of the Interaction attached to this Integrator

virtual void updateInput(double time) override: update the input of the Interaction attached to this Integrator

virtual void updateOutput(double time, unsigned int level) override

update the output of the Interaction attached to this Integrator

Parameters:

time – current time
level – level of interest for the dynamics

virtual void updateInput(double time, unsigned int level) override

update the input of the Interaction attached to this Integrator

Parameters:

time – current time
level – level of interest for the dynamics

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

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

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

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

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

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

Parameters:

vertex_inter – 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 prepareNewtonIteration(double time) override

computes all the W matrices

Parameters:: time – current time

virtual void integrate(double &tinit, double &tend, double &tout, int &useless) 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
useless – flag (for EulerMoreauOSI, used in LsodarOSI)

virtual void updateState(const unsigned int level) override

updates the state of the Dynamical Systems

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

virtual void display() override: Displays the data of the EulerMoreauOSI’s integrator.

ACCEPT_STD_VISITORS()#

Protected Functions

ACCEPT_SERIALIZATION(EulerMoreauOSI)#

inline EulerMoreauOSI()#: Default constructor.

Protected Attributes

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

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$.

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 struct _NSLEffectOnFreeOutput

File kernel/src/simulationTools/EulerMoreauOSI.hpp

Contents

File kernel/src/simulationTools/EulerMoreauOSI.hpp#