Coverage for pyxfem/core/dynamic_solver.py: 21%
90 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1"""Explicit and implicit time integration schemes for dynamics.
3This module provides second-order ordinary differential equation (ODE)
4time integrators suited for dynamic finite element systems:
5- Velocity-Verlet (Störmer-Verlet) explicit integration.
6- Newmark-beta family of integrators using Richardson fixed-point
7 or Newton-Raphson nonlinear solvers.
8"""
9# pylint: disable=invalid-name, too-few-public-methods
10import numpy as np
13class Verlet:
14 """Velocity-Verlet (Störmer-Verlet) explicit time integration scheme.
16 Integrates equations of motion of the form:
17 a = acc(u, v)
19 Parameters
20 ----------
21 acc : callable
22 Function evaluating the acceleration vector `a = acc(u, v)`
23 given displacement `u` and velocity `v`.
24 verbose : bool, optional
25 If True, prints diagnostic information (default is False).
26 """
27 def __init__(self, acc, verbose = False):
28 self.acc = acc
29 self.verbose = verbose
31 def step(self, t, u, v, a, dt):
32 """Advance the solution by a single Velocity-Verlet time step.
34 Parameters
35 ----------
36 t : float
37 Current simulation time.
38 u : NDArray
39 Current displacement vector at time `t`.
40 v : NDArray
41 Current velocity vector at time `t`.
42 a : NDArray
43 Current acceleration vector at time `t`.
44 dt : float
45 Time step increment.
47 Returns
48 -------
49 t_next : float
50 Updated time (`t + dt`).
51 u_next : NDArray
52 Displacement vector at time `t + dt`.
53 v_next : NDArray
54 Velocity vector at time `t + dt`.
55 a_next : NDArray
56 Acceleration vector at time `t + dt`.
57 """
58 u = u+ dt*v + dt*dt/2*a
59 a1 = self.acc(u, v)
60 v = v + 0.5*dt*(a+a1)
61 a = a1
62 t = t+dt
63 return t, u ,v, a
65class Newmark:
66 """Newmark-beta family of implicit and explicit time integration schemes.
68 Solves second-order dynamic equilibrium systems:
69 M(u) * a = F(u)
71 Constructeur :
72 Parameters
73 ----------
74 getF : callable
75 Function returning the force vector `F(u)` for displacement `u`.
76 getM : callable
77 Function returning the mass matrix `M(u)` for displacement `u`.
78 getK : callable or None
79 Function returning the tangent stiffness matrix `K(u)` for displacement `u`.
80 Required when `solver='newton'`.
81 linsolve : callable or None
82 Linear solver function `(A, b) -> x`. Required when `solver='newton'`.
83 gamma : float, optional
84 Newmark parameter governing numerical damping (default is 0.5,
85 ensuring second-order accuracy without numerical dissipation).
86 beta : float, optional
87 Newmark parameter governing acceleration interpolation (default is 0.25,
88 corresponding to average acceleration / unconditional stability).
89 If `beta=0.0` and `gamma=0.5`, the scheme simplifies to explicit Verlet.
90 solver : {'richardson', 'newton'}, optional
91 Nonlinear iteration strategy used for implicit steps (default is 'richardson').
92 verbose : bool, optional
93 If True, prints iteration residuals during time stepping (default is False).
95 Raises
96 ------
97 ValueError
98 If `solver='newton'` and either `getK` or `linsolve` is None.
99 """
101 def __init__(self, getF, getM, getK, linsolve, \
102 gamma = 0.5, beta = 0.25, solver = 'richardson', verbose =False):
104 self.getF = getF
105 self.getM = getM
106 self.getK = getK
107 self.linsolve = linsolve
108 self.acc = lambda u : linsolve(self.getM(u), self.getF(u))
109 self.verbose = verbose
110 self.gamma = gamma
111 self.beta = beta
112 self.eps_abs = 1.e-8
113 self.eps_rel = 1.e-6
114 self.itmax = 10
116 if self.beta == 0. and gamma == 0.5:
117 self.step = Verlet(self.acc).step
118 elif solver == 'richardson':
119 self.step = self.step_richardson
120 else:
121 if self.getK is None :
122 raise ValueError('getK required for newton solver')
123 if self.linsolve is None :
124 raise ValueError('linsolve required for newton solver')
125 self.step = self.step_newton
127 def _up_u(self, u0, v0, a0, a, dt):
128 return u0 + dt*v0 + dt*dt/2*(( 1-2*self.beta)*a0+ 2*self.beta*a)
130 def _up_v(self, v0, a0, a, dt):
131 return v0 + dt*((1-self.gamma)*a0 + self.gamma*a)
133 def _up(self, u0, v0, a0, a, dt):
134 u = self._up_u(u0, v0, a0, a, dt)
135 v = self._up_v(v0, a0, a, dt)
136 a = self.acc(u)
137 return u, v, a
139 def step_richardson(self, t, u0, v0, a0, dt):
140 """Advance the solution by a single step using Richardson fixed-point iterations.
142 Solves the implicit Newmark-beta equations by iteratively updating
143 acceleration, displacement, and velocity until successive acceleration
144 increments satisfy convergence tolerances or `itmax` is reached.
146 Parameters
147 ----------
148 t : float
149 Current simulation time.
150 u0 : NDArray
151 Current displacement vector at time `t`.
152 v0 : NDArray
153 Current velocity vector at time `t`.
154 a0 : NDArray
155 Current acceleration vector at time `t`.
156 dt : float
157 Time step increment.
159 Returns
160 -------
161 t_next : float
162 Updated simulation time (`t + dt`).
163 u : NDArray
164 Displacement vector at time `t + dt`.
165 v : NDArray
166 Velocity vector at time `t + dt`.
167 a : NDArray
168 Acceleration vector at time `t + dt`.
169 """
170 a = a0.copy()
171 u, v, a = self._up(u0, v0, a0, a, dt)
172 nr = np.linalg.norm(a-a0)
173 k = 0
174 if self.verbose :
175 print('newmarkstep - richardson step ', k, nr)
176 nr0= np.linalg.norm(a0)
178 while (nr > self.eps_abs and nr > self.eps_rel*nr0) and k < self.itmax :
179 a0 = a.copy()
180 u, v, a = self._up(u0, v0, a0, a, dt)
181 nr = np.linalg.norm(a-a0)
182 k=k+1
183 if self.verbose :
184 print('newmarkr - richardson step ', k, nr)
185 if k >= self.itmax:
186 print("WARNING : newmark - richardson not conveged after ", k, 'iterations')
187 return t+dt, u ,v, a
189 def step_newton(self, t, u0, v0, a0, dt):
190 """Advance the solution by a single step using Newton-Raphson iterations.
192 Solves the nonlinear dynamic equilibrium equation for displacement `u`
193 at time `t + dt` using the effective dynamic tangent matrix:
194 K_eff = M / (beta * dt^2) + K
195 Once converged, velocity and acceleration are updated via Newmark formulae.
197 Parameters
198 ----------
199 t : float
200 Current simulation time.
201 u0 : NDArray
202 Current displacement vector at time `t`.
203 v0 : NDArray
204 Current velocity vector at time `t`.
205 a0 : NDArray
206 Current acceleration vector at time `t`.
207 dt : float
208 Time step increment.
210 Returns
211 -------
212 t_next : float
213 Updated simulation time (`t + dt`).
214 u : NDArray
215 Displacement vector at time `t + dt`.
216 v : NDArray
217 Velocity vector at time `t + dt`.
218 a : NDArray
219 Acceleration vector at time `t + dt`.
220 """
221 beta = self.beta
222 def res(u):
223 M = self.getM(u)
224 F = self.getF(u)
225 return M@(u-u0 -dt*v0 - 0.5*(1.-2.*beta)*dt*dt*a0)/(beta*dt*dt) - F
226 def dres(u):
227 M = self.getM(u)
228 K = self.getK(u)
229 return M*1./(beta*dt*dt) + K
230 u = self._up_u(u0, v0, a0, a0, dt)
231 r = res(u)
232 nr = np.linalg.norm(r)
233 k = 0
234 if self.verbose :
235 print('newmark - newton step ', k, nr)
236 nr0= np.linalg.norm(a0)
237 while (nr > self.eps_abs and nr > self.eps_rel*nr0) and k < self.itmax :
238 u -= self.linsolve(dres(u), r)
239 r = res(u)
240 nr = np.linalg.norm(r)
241 k += 1
242 if self.verbose :
243 print('newmark - newton step ', k, nr)
244 if k >= self.itmax:
245 print("WARNING newmark - newton not conveged after ", k, 'iterations')
246 a = self.acc(u)
247 v = self._up_v(v0, a0, a, dt)
248 return t+dt, u, v, a