elastodynamicsx.pde.timeschemes

Module contents

The timeschemes module contains tools to construct the weak form of a time-dependent problem. Several implicit and explicit schemes are supported. The preferred way to use this module is however through the elastodynamicsx.solvers.timestepper function.

List of supported schemes:
Explicit schemes:

  • ‘leapfrog’

Implicit schemes:

  • ‘midpoint’

  • ‘linear-acceleration-method’

  • ‘newmark’

  • ‘hht-alpha’

  • ‘generalized-alpha’

elastodynamicsx.pde.timeschemes.timescheme(*args, **kwargs)[source]

Builder method that instanciates the desired time scheme

Parameters:

*args – Passed to the required time scheme

Keyword Arguments:

  • scheme – (mandatory) String identifier of the time scheme

  • **kwargs – Passed to the required material

Returns:

An instance of the desired time scheme

class elastodynamicsx.pde.timeschemes.TimeScheme(dt, out: Vec, **kwargs)[source]

Bases: object

Abstract base class for time schemes as needed by the TimeStepper solvers

Parameters:

  • dt – Time step

  • out – Solution vector

Keyword Arguments:

  • explicit (mandatory) – Whether the scheme is explicit or implicit

  • nbsteps (mandatory) – Number of resolution steps

  • linear_ODE (default=True) – Whether the resulting ODE is linear or not

  • intermediate_dt (default=0.) – Time at which time-dependent terms are evaluated

labels: List[str] = ['supercharge me']
petsc_options_t0: dict = {'ksp_type': 'preonly', 'pc_type': 'lu'}
property explicit: bool
property nbsteps: int
property linear_ODE: bool
property dt

The time step

property intermediate_dt
property out: Vec

The solution vector

b_update_function(b: Vec, t) → None[source]
prepareNextIteration() → None[source]
set_initial_condition(u0, v0) → None[source]
initialStep(t0, callfirsts: List[Callable] = [], callbacks: List[Callable] = [], verbose: int = 0) → None[source]

Specific to the initial value step

class elastodynamicsx.pde.timeschemes.FEniCSxTimeScheme(dt, out: Function, bilinear_form: Form, linear_form: Form, mpc: MultiPointConstraint | None = None, bcs=[], **kwargs)[source]

Bases: TimeScheme

Abstract base class based on FEniCSx’s form language

property out_fenicsx: Function

The solution vector

A() → Mat[source]

The time-independent matrix (bilinear form)

init_b() → Vec[source]

Declares a zero vector compatible with the linear form

b_update_function(b: Vec, t) → None[source]

Updates the b vector (in-place) for a given time t

set_initial_condition(u0, v0) → None[source]

Apply initial conditions

Parameters:

  • u0 – u at t0

  • v0 –

    du/dt at t0

    u0 and v0 can be:

    • Python Callable -> interpolated at nodes

      -> e.g. u0 = lambda x: np.zeros((domain.topology.dim, x.shape[1]), dtype=PETSc.ScalarType)

    • scalar (int, float, complex, PETSc.ScalarType)

      -> e.g. u0 = 0

    • array (list, tuple, np.ndarray)

      -> e.g. u0 = [0,0,0]

    • fem.function.Function or fem.function.Constant

      -> e.g. u0 = fem.Function(V)

class elastodynamicsx.pde.timeschemes.LeapFrog(function_space: FunctionSpace, M_fn: Callable[[TrialFunction, TestFunction], Form], C_fn: None | Callable[[TrialFunction, TestFunction], Form], K_fn: Callable[[TrialFunction, TestFunction], Form], b_fn: None | Callable[[TestFunction], Form], dt, bcs: Tuple[BoundaryConditionBase] | Tuple[()] = (), **kwargs)[source]

Bases: FEniCSxTimeScheme

Implementation of the Leapfrog time-stepping scheme, or Explicit central difference scheme. Leapfrog is a special case of Newmark-beta methods with \(\beta=0\) and \(\gamma=0.5\)

Scheme:
\(u_n\), \(v_n\), \(a_n\) are the (known) displacement, velocity and acceleration at current time step
\(u_{n+1}\) is the unknown: displacement at next time step
\(a_n = (u_{n+1} - 2 u_n + u_{n-1}) / dt^2\)
\(v_n = (u_{n+1} - u_n-1) / (2 dt)\)
Implicit / explicit?

Explicit

Accuracy:

Second-order

Parameters:

  • function_space – The Finite Element functionnal space

  • M_fn –

    The mass form. Usually:

    M_fn = lambda u,v: rho* ufl.dot(u, v) * ufl.dx

  • C_fn (optional, ignored if None) –

    The damping form. E.g. for Rayleigh damping:

    C_fn = lambda u,v: eta_m * M_fn(u,v) + eta_k * K_fn(u,v)

  • K_fn –

    The stiffness form. Usually:

    K_fn = lambda u,v: ufl.inner(sigma(u), epsilon(v)) * ufl.dx

  • b_fn (optional, ignored if None) – Right hand term

  • dt – Time step

  • bcs – The set of boundary conditions

Keyword Arguments:

jit_options – (default=PDECONFIG.default_jit_options) options for the just-in-time compiler

See:

https://en.wikipedia.org/wiki/Leapfrog_integration

labels: List[str] = ['leapfrog', 'central-difference']
property u: Function

The displacement field at current time step

property u_nm1: Function

The displacement field at previous time step

prepareNextIteration() → None[source]

Next-time-step function, to prepare next iteration -> Call it after solving

initialStep(t0, callfirsts: List[Callable] = [], callbacks: List[Callable] = [], verbose: int = 0) → None[source]

Specific to the initial value step

class elastodynamicsx.pde.timeschemes.GalphaNewmarkBeta(function_space: FunctionSpace, M_fn: Callable[[TrialFunction, TestFunction], Form], C_fn: None | Callable[[TrialFunction, TestFunction], Form], K_fn: Callable[[TrialFunction, TestFunction], Form], b_fn: None | Callable[[TestFunction], Form], dt, bcs: Tuple[BoundaryConditionBase] | Tuple[()] = (), **kwargs)[source]

Bases: FEniCSxTimeScheme

Implementation of the ‘g-a-newmark’ (Generalized-alpha Newmark) time-stepping scheme, for \(\beta>0\). The special case \(\beta=0\) is implemented in the LeapFrog class.

Implicit / explicit?

Implicit

Accuracy:

First or second order, depending on parameters

Parameters:

  • function_space – The Finite Element functionnal space

  • M_fn –

    The mass form. Usually:

    M_fn = lambda u,v: rho* ufl.dot(u, v) * ufl.dx

  • C_fn (optional, ignored if None) –

    The damping form. E.g. for Rayleigh damping:

    C_fn = lambda u,v: eta_m * M_fn(u,v) + eta_k * K_fn(u,v)

  • K_fn –

    The stiffness form. Usually:

    K_fn = lambda u,v: ufl.inner(sigma(u), epsilon(v)) * ufl.dx

  • b_fn (optional, ignored if None) – Linear form

  • dt – Time step

  • bcs – The set of boundary conditions

Keyword Arguments:

  • rho_inf (default = 0.75) – Spectral radius in the high frequency limit. Bounds: (1/2,1). Setting rho_inf is the preferred way of defining the scheme.

  • alpha_m (default = (2*rho_inf-1)/(rho_inf+1)) – Unconditionnal stability if \(-1 \leq \alpha_m \leq \alpha_f \leq 0.5\)

  • alpha_f (default = rho_inf/(rho_inf+1)) – Unconditionnal stability if \(-1 \leq \alpha_m \leq \alpha_f \leq 0.5\)

  • gamma (default = 1/2 - alpha_m + alpha_f) – Default value ensures second-order accuracy. Other values give first-order accuracy

  • beta (default = 1/4*(gamma+1/2)**2) – Unconditionnal stability if \(\beta \geq 0.25 + 0.5 (\alpha_f - \alpha_m)\)

  • jit_options (default=PDECONFIG.default_jit_options) – Options for the just-in-time compiler

Reference:

J. Chung and G. M. Hulbert, “A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method,” J. Appl. Mech. 60(2), 371–375 (1993).

Adapted from:

https://comet-fenics.readthedocs.io/en/latest/demo/elastodynamics/demo_elastodynamics.py.html

labels: List[str] = ['g-a-newmark', 'generalized-alpha']
property u: Function
property v: Function
property a: Function
prepareNextIteration() → None[source]

Next-time-step function, to prepare next iteration -> Call it after solving

initialStep(t0, callfirsts: List[Callable] = [], callbacks: List[Callable] = [], verbose: int = 0) → None[source]

Specific to the initial value step

class elastodynamicsx.pde.timeschemes.HilberHughesTaylor(*args, **kwargs)[source]

Bases: GalphaNewmarkBeta

Implementation of the Hilber-Hughes-Taylor or alpha-Newmark time-stepping scheme. HHT is a special case of the Generalized-alpha scheme (\(\alpha_m=0\) and \(\alpha_f=\alpha\)).

(!) The definition of the \(\alpha\) parameter is different from that in the original article. Here: \(\alpha = -\alpha_{HHT}\)

Implicit / explicit?

Implicit

Accuracy:

First or second order, depending on parameters

Parameters:

*args – See GalphaNewmarkBeta

Keyword Arguments:

  • rho_inf (default = 0.9) – Spectral radius in the high frequency limit. Bounds: (1/2,1). Setting rho_inf is the preferred way of defining the scheme.

  • alpha (default = (1-rho_inf)/(1+rho_inf)) – Set to e.g. 0.05 or sqrt(2) for low or moderate dissipation

  • gamma (default = 1/2 + alpha) – Default value ensures second-order accuracy. Other values give first-order accuracy

  • beta (default = 1/4*(gamma+1/2)**2) – Unconditionnal stability if \(\beta \geq 0.25 + 0.5 \alpha\)

Reference:

H. M. Hilber, T. J. R. Hughes, and R. L. Taylor, “Improved Numerical Dissipation for Time Integration Algorithms in Structural Dynamics”, Earthquake Engineering and Structural Dynamics, vol. 5, pp. 283–292, 1977.

labels: List[str] = ['hilber-hughes-taylor', 'hht', 'hht-alpha']
class elastodynamicsx.pde.timeschemes.NewmarkBeta(*args, **kwargs)[source]

Bases: GalphaNewmarkBeta

Implementation of the Newmark or Newmark-beta time-stepping scheme, for \(\beta>0\). The special case \(\beta=0\) is implemented in the LeapFrog class. Newmark-beta is a special case of the Generalized-alpha scheme (\(\alpha_m=\alpha_f=0\)).

Implicit / explicit?

Implicit

Accuracy:

First-order unless \(gamma=1/2\) (second-order)

Parameters:

*args – See GalphaNewmarkBeta

Keyword Arguments:

  • gamma (default = 1/2) – Default value ensures second-order accuracy. Other values give first-order accuracy

  • beta (default = 1/4*(gamma+1/2)**2) – Unconditionnal stability if \(\beta \geq 0.25\)

See:

https://en.wikipedia.org/wiki/Newmark-beta_method

labels: List[str] = ['newmark', 'newmark-beta']
class elastodynamicsx.pde.timeschemes.MidPoint(*args, **kwargs)[source]

Bases: NewmarkBeta

Implementation of the Midpoint time-stepping scheme, or Average Acceleration Method. Midpoint is a special case of Newmark-beta methods with \(\beta=1/4\) and \(\gamma=1/2\) and is unconditionally stable.

Implicit / explicit?

Implicit

Accuracy:

Second-order

Parameters:

*args – See GalphaNewmarkBeta

See:

https://en.wikipedia.org/wiki/Midpoint_method

labels: List[str] = ['midpoint', 'average-acceleration-method', 'aam']
class elastodynamicsx.pde.timeschemes.LinearAccelerationMethod(*args, **kwargs)[source]

Bases: NewmarkBeta

Implementation Linear Acceleration Method time-stepping scheme. It is a special case of Newmark-beta methods with \(\beta=1/6\) and \(\gamma=1/2\) and is unconditionally stable.

Implicit / explicit?

Implicit

Accuracy:

Second-order

Parameters:

*args – See GalphaNewmarkBeta

labels: List[str] = ['linear-acceleration-method', 'lam']

elastodynamicsx.pde.timeschemes.leapfrog

class elastodynamicsx.pde.timeschemes.leapfrog.LeapFrog(function_space: FunctionSpace, M_fn: Callable[[TrialFunction, TestFunction], Form], C_fn: None | Callable[[TrialFunction, TestFunction], Form], K_fn: Callable[[TrialFunction, TestFunction], Form], b_fn: None | Callable[[TestFunction], Form], dt, bcs: Tuple[BoundaryConditionBase] | Tuple[()] = (), **kwargs)[source]

Bases: FEniCSxTimeScheme

Implementation of the Leapfrog time-stepping scheme, or Explicit central difference scheme. Leapfrog is a special case of Newmark-beta methods with \(\beta=0\) and \(\gamma=0.5\)

Scheme:
\(u_n\), \(v_n\), \(a_n\) are the (known) displacement, velocity and acceleration at current time step
\(u_{n+1}\) is the unknown: displacement at next time step
\(a_n = (u_{n+1} - 2 u_n + u_{n-1}) / dt^2\)
\(v_n = (u_{n+1} - u_n-1) / (2 dt)\)
Implicit / explicit?

Explicit

Accuracy:

Second-order

Parameters:

  • function_space – The Finite Element functionnal space

  • M_fn –

    The mass form. Usually:

    M_fn = lambda u,v: rho* ufl.dot(u, v) * ufl.dx

  • C_fn (optional, ignored if None) –

    The damping form. E.g. for Rayleigh damping:

    C_fn = lambda u,v: eta_m * M_fn(u,v) + eta_k * K_fn(u,v)

  • K_fn –

    The stiffness form. Usually:

    K_fn = lambda u,v: ufl.inner(sigma(u), epsilon(v)) * ufl.dx

  • b_fn (optional, ignored if None) – Right hand term

  • dt – Time step

  • bcs – The set of boundary conditions

Keyword Arguments:

jit_options – (default=PDECONFIG.default_jit_options) options for the just-in-time compiler

See:

https://en.wikipedia.org/wiki/Leapfrog_integration

labels: List[str] = ['leapfrog', 'central-difference']
property u: Function

The displacement field at current time step

property u_nm1: Function

The displacement field at previous time step

prepareNextIteration() → None[source]

Next-time-step function, to prepare next iteration -> Call it after solving

initialStep(t0, callfirsts: List[Callable] = [], callbacks: List[Callable] = [], verbose: int = 0) → None[source]

Specific to the initial value step

elastodynamicsx.pde.timeschemes.newmark

class elastodynamicsx.pde.timeschemes.newmark.GalphaNewmarkBeta(function_space: FunctionSpace, M_fn: Callable[[TrialFunction, TestFunction], Form], C_fn: None | Callable[[TrialFunction, TestFunction], Form], K_fn: Callable[[TrialFunction, TestFunction], Form], b_fn: None | Callable[[TestFunction], Form], dt, bcs: Tuple[BoundaryConditionBase] | Tuple[()] = (), **kwargs)[source]

Bases: FEniCSxTimeScheme

Implementation of the ‘g-a-newmark’ (Generalized-alpha Newmark) time-stepping scheme, for \(\beta>0\). The special case \(\beta=0\) is implemented in the LeapFrog class.

Implicit / explicit?

Implicit

Accuracy:

First or second order, depending on parameters

Parameters:

  • function_space – The Finite Element functionnal space

  • M_fn –

    The mass form. Usually:

    M_fn = lambda u,v: rho* ufl.dot(u, v) * ufl.dx

  • C_fn (optional, ignored if None) –

    The damping form. E.g. for Rayleigh damping:

    C_fn = lambda u,v: eta_m * M_fn(u,v) + eta_k * K_fn(u,v)

  • K_fn –

    The stiffness form. Usually:

    K_fn = lambda u,v: ufl.inner(sigma(u), epsilon(v)) * ufl.dx

  • b_fn (optional, ignored if None) – Linear form

  • dt – Time step

  • bcs – The set of boundary conditions

Keyword Arguments:

  • rho_inf (default = 0.75) – Spectral radius in the high frequency limit. Bounds: (1/2,1). Setting rho_inf is the preferred way of defining the scheme.

  • alpha_m (default = (2*rho_inf-1)/(rho_inf+1)) – Unconditionnal stability if \(-1 \leq \alpha_m \leq \alpha_f \leq 0.5\)

  • alpha_f (default = rho_inf/(rho_inf+1)) – Unconditionnal stability if \(-1 \leq \alpha_m \leq \alpha_f \leq 0.5\)

  • gamma (default = 1/2 - alpha_m + alpha_f) – Default value ensures second-order accuracy. Other values give first-order accuracy

  • beta (default = 1/4*(gamma+1/2)**2) – Unconditionnal stability if \(\beta \geq 0.25 + 0.5 (\alpha_f - \alpha_m)\)

  • jit_options (default=PDECONFIG.default_jit_options) – Options for the just-in-time compiler

Reference:

J. Chung and G. M. Hulbert, “A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method,” J. Appl. Mech. 60(2), 371–375 (1993).

Adapted from:

https://comet-fenics.readthedocs.io/en/latest/demo/elastodynamics/demo_elastodynamics.py.html

labels: List[str] = ['g-a-newmark', 'generalized-alpha']
property u: Function
property v: Function
property a: Function
prepareNextIteration() → None[source]

Next-time-step function, to prepare next iteration -> Call it after solving

initialStep(t0, callfirsts: List[Callable] = [], callbacks: List[Callable] = [], verbose: int = 0) → None[source]

Specific to the initial value step

class elastodynamicsx.pde.timeschemes.newmark.HilberHughesTaylor(*args, **kwargs)[source]

Bases: GalphaNewmarkBeta

Implementation of the Hilber-Hughes-Taylor or alpha-Newmark time-stepping scheme. HHT is a special case of the Generalized-alpha scheme (\(\alpha_m=0\) and \(\alpha_f=\alpha\)).

(!) The definition of the \(\alpha\) parameter is different from that in the original article. Here: \(\alpha = -\alpha_{HHT}\)

Implicit / explicit?

Implicit

Accuracy:

First or second order, depending on parameters

Parameters:

*args – See GalphaNewmarkBeta

Keyword Arguments:

  • rho_inf (default = 0.9) – Spectral radius in the high frequency limit. Bounds: (1/2,1). Setting rho_inf is the preferred way of defining the scheme.

  • alpha (default = (1-rho_inf)/(1+rho_inf)) – Set to e.g. 0.05 or sqrt(2) for low or moderate dissipation

  • gamma (default = 1/2 + alpha) – Default value ensures second-order accuracy. Other values give first-order accuracy

  • beta (default = 1/4*(gamma+1/2)**2) – Unconditionnal stability if \(\beta \geq 0.25 + 0.5 \alpha\)

Reference:

H. M. Hilber, T. J. R. Hughes, and R. L. Taylor, “Improved Numerical Dissipation for Time Integration Algorithms in Structural Dynamics”, Earthquake Engineering and Structural Dynamics, vol. 5, pp. 283–292, 1977.

labels: List[str] = ['hilber-hughes-taylor', 'hht', 'hht-alpha']
class elastodynamicsx.pde.timeschemes.newmark.NewmarkBeta(*args, **kwargs)[source]

Bases: GalphaNewmarkBeta

Implementation of the Newmark or Newmark-beta time-stepping scheme, for \(\beta>0\). The special case \(\beta=0\) is implemented in the LeapFrog class. Newmark-beta is a special case of the Generalized-alpha scheme (\(\alpha_m=\alpha_f=0\)).

Implicit / explicit?

Implicit

Accuracy:

First-order unless \(gamma=1/2\) (second-order)

Parameters:

*args – See GalphaNewmarkBeta

Keyword Arguments:

  • gamma (default = 1/2) – Default value ensures second-order accuracy. Other values give first-order accuracy

  • beta (default = 1/4*(gamma+1/2)**2) – Unconditionnal stability if \(\beta \geq 0.25\)

See:

https://en.wikipedia.org/wiki/Newmark-beta_method

labels: List[str] = ['newmark', 'newmark-beta']
class elastodynamicsx.pde.timeschemes.newmark.MidPoint(*args, **kwargs)[source]

Bases: NewmarkBeta

Implementation of the Midpoint time-stepping scheme, or Average Acceleration Method. Midpoint is a special case of Newmark-beta methods with \(\beta=1/4\) and \(\gamma=1/2\) and is unconditionally stable.

Implicit / explicit?

Implicit

Accuracy:

Second-order

Parameters:

*args – See GalphaNewmarkBeta

See:

https://en.wikipedia.org/wiki/Midpoint_method

labels: List[str] = ['midpoint', 'average-acceleration-method', 'aam']
class elastodynamicsx.pde.timeschemes.newmark.LinearAccelerationMethod(*args, **kwargs)[source]

Bases: NewmarkBeta

Implementation Linear Acceleration Method time-stepping scheme. It is a special case of Newmark-beta methods with \(\beta=1/6\) and \(\gamma=1/2\) and is unconditionally stable.

Implicit / explicit?

Implicit

Accuracy:

Second-order

Parameters:

*args – See GalphaNewmarkBeta

labels: List[str] = ['linear-acceleration-method', 'lam']

elastodynamicsx.pde.timeschemes.timescheme

Builder method that instanciates the desired time scheme

param *args:

Passed to the required time scheme

keyword scheme:

(mandatory) String identifier of the time scheme

keyword **kwargs:

Passed to the required material

returns:

An instance of the desired time scheme