"""Explicit and implicit time integration schemes for dynamics.
This module provides second-order ordinary differential equation (ODE)
time integrators suited for dynamic finite element systems:
- Velocity-Verlet (Störmer-Verlet) explicit integration.
- Newmark-beta family of integrators using Richardson fixed-point
or Newton-Raphson nonlinear solvers.
"""
# pylint: disable=invalid-name, too-few-public-methods
import numpy as np
[docs]
class Verlet:
"""Velocity-Verlet (Störmer-Verlet) explicit time integration scheme.
Integrates equations of motion of the form:
a = acc(u, v)
Parameters
----------
acc : callable
Function evaluating the acceleration vector `a = acc(u, v)`
given displacement `u` and velocity `v`.
verbose : bool, optional
If True, prints diagnostic information (default is False).
"""
def __init__(self, acc, verbose = False):
self.acc = acc
self.verbose = verbose
[docs]
def step(self, t, u, v, a, dt):
"""Advance the solution by a single Velocity-Verlet time step.
Parameters
----------
t : float
Current simulation time.
u : NDArray
Current displacement vector at time `t`.
v : NDArray
Current velocity vector at time `t`.
a : NDArray
Current acceleration vector at time `t`.
dt : float
Time step increment.
Returns
-------
t_next : float
Updated time (`t + dt`).
u_next : NDArray
Displacement vector at time `t + dt`.
v_next : NDArray
Velocity vector at time `t + dt`.
a_next : NDArray
Acceleration vector at time `t + dt`.
"""
u = u+ dt*v + dt*dt/2*a
a1 = self.acc(u, v)
v = v + 0.5*dt*(a+a1)
a = a1
t = t+dt
return t, u ,v, a
[docs]
class Newmark:
"""Newmark-beta family of implicit and explicit time integration schemes.
Solves second-order dynamic equilibrium systems:
M(u) * a = F(u)
Constructeur :
Parameters
----------
getF : callable
Function returning the force vector `F(u)` for displacement `u`.
getM : callable
Function returning the mass matrix `M(u)` for displacement `u`.
getK : callable or None
Function returning the tangent stiffness matrix `K(u)` for displacement `u`.
Required when `solver='newton'`.
linsolve : callable or None
Linear solver function `(A, b) -> x`. Required when `solver='newton'`.
gamma : float, optional
Newmark parameter governing numerical damping (default is 0.5,
ensuring second-order accuracy without numerical dissipation).
beta : float, optional
Newmark parameter governing acceleration interpolation (default is 0.25,
corresponding to average acceleration / unconditional stability).
If `beta=0.0` and `gamma=0.5`, the scheme simplifies to explicit Verlet.
solver : {'richardson', 'newton'}, optional
Nonlinear iteration strategy used for implicit steps (default is 'richardson').
verbose : bool, optional
If True, prints iteration residuals during time stepping (default is False).
Raises
------
ValueError
If `solver='newton'` and either `getK` or `linsolve` is None.
"""
def __init__(self, getF, getM, getK, linsolve, \
gamma = 0.5, beta = 0.25, solver = 'richardson', verbose =False):
self.getF = getF
self.getM = getM
self.getK = getK
self.linsolve = linsolve
self.acc = lambda u : linsolve(self.getM(u), self.getF(u))
self.verbose = verbose
self.gamma = gamma
self.beta = beta
self.eps_abs = 1.e-8
self.eps_rel = 1.e-6
self.itmax = 10
if self.beta == 0. and gamma == 0.5:
self.step = Verlet(self.acc).step
elif solver == 'richardson':
self.step = self.step_richardson
else:
if self.getK is None :
raise ValueError('getK required for newton solver')
if self.linsolve is None :
raise ValueError('linsolve required for newton solver')
self.step = self.step_newton
def _up_u(self, u0, v0, a0, a, dt):
return u0 + dt*v0 + dt*dt/2*(( 1-2*self.beta)*a0+ 2*self.beta*a)
def _up_v(self, v0, a0, a, dt):
return v0 + dt*((1-self.gamma)*a0 + self.gamma*a)
def _up(self, u0, v0, a0, a, dt):
u = self._up_u(u0, v0, a0, a, dt)
v = self._up_v(v0, a0, a, dt)
a = self.acc(u)
return u, v, a
[docs]
def step_richardson(self, t, u0, v0, a0, dt):
"""Advance the solution by a single step using Richardson fixed-point iterations.
Solves the implicit Newmark-beta equations by iteratively updating
acceleration, displacement, and velocity until successive acceleration
increments satisfy convergence tolerances or `itmax` is reached.
Parameters
----------
t : float
Current simulation time.
u0 : NDArray
Current displacement vector at time `t`.
v0 : NDArray
Current velocity vector at time `t`.
a0 : NDArray
Current acceleration vector at time `t`.
dt : float
Time step increment.
Returns
-------
t_next : float
Updated simulation time (`t + dt`).
u : NDArray
Displacement vector at time `t + dt`.
v : NDArray
Velocity vector at time `t + dt`.
a : NDArray
Acceleration vector at time `t + dt`.
"""
a = a0.copy()
u, v, a = self._up(u0, v0, a0, a, dt)
nr = np.linalg.norm(a-a0)
k = 0
if self.verbose :
print('newmarkstep - richardson step ', k, nr)
nr0= np.linalg.norm(a0)
while (nr > self.eps_abs and nr > self.eps_rel*nr0) and k < self.itmax :
a0 = a.copy()
u, v, a = self._up(u0, v0, a0, a, dt)
nr = np.linalg.norm(a-a0)
k=k+1
if self.verbose :
print('newmarkr - richardson step ', k, nr)
if k >= self.itmax:
print("WARNING : newmark - richardson not conveged after ", k, 'iterations')
return t+dt, u ,v, a
[docs]
def step_newton(self, t, u0, v0, a0, dt):
"""Advance the solution by a single step using Newton-Raphson iterations.
Solves the nonlinear dynamic equilibrium equation for displacement `u`
at time `t + dt` using the effective dynamic tangent matrix:
K_eff = M / (beta * dt^2) + K
Once converged, velocity and acceleration are updated via Newmark formulae.
Parameters
----------
t : float
Current simulation time.
u0 : NDArray
Current displacement vector at time `t`.
v0 : NDArray
Current velocity vector at time `t`.
a0 : NDArray
Current acceleration vector at time `t`.
dt : float
Time step increment.
Returns
-------
t_next : float
Updated simulation time (`t + dt`).
u : NDArray
Displacement vector at time `t + dt`.
v : NDArray
Velocity vector at time `t + dt`.
a : NDArray
Acceleration vector at time `t + dt`.
"""
beta = self.beta
def res(u):
M = self.getM(u)
F = self.getF(u)
return M@(u-u0 -dt*v0 - 0.5*(1.-2.*beta)*dt*dt*a0)/(beta*dt*dt) - F
def dres(u):
M = self.getM(u)
K = self.getK(u)
return M*1./(beta*dt*dt) + K
u = self._up_u(u0, v0, a0, a0, dt)
r = res(u)
nr = np.linalg.norm(r)
k = 0
if self.verbose :
print('newmark - newton step ', k, nr)
nr0= np.linalg.norm(a0)
while (nr > self.eps_abs and nr > self.eps_rel*nr0) and k < self.itmax :
u -= self.linsolve(dres(u), r)
r = res(u)
nr = np.linalg.norm(r)
k += 1
if self.verbose :
print('newmark - newton step ', k, nr)
if k >= self.itmax:
print("WARNING newmark - newton not conveged after ", k, 'iterations')
a = self.acc(u)
v = self._up_v(v0, a0, a, dt)
return t+dt, u, v, a