pyyeti.ode.SolveExp2.generator

SolveExp2.generator(nt, F0, d0=None, v0=None, static_ic=False)

Python “generator” version of SolveExp2.tsolve(); interactively solve (or re-solve) one step at a time.

Parameters:

• nt (integer) – Number of time steps

• F0 (1d ndarray) – Initial force vector

• d0 (1d ndarray or None; optional) – Displacement initial conditions; if None, zero ic’s are used.

• v0 (1d ndarray or None; optional) – Velocity initial conditions; if None, zero ic’s are used.

• static_ic (bool; optional) – If True and d0 is None, then d0 is calculated such that static (steady-state) initial conditions are used. Uses the pseudo-inverse in case there are rigid-body modes. static_ic is ignored if d0 is not None.

Returns:

• gen (generator function) – Generator function for solving equations one step at a time

• d, v (2d ndarrays) – The displacement and velocity arrays. Only the first column of d and v are set; other values are all zero.

Notes

To use the generator:

1. Instantiate a SolveExp2 instance:

   ts = SolveExp2(m, b, k, h)
   

2. Retrieve a generator and the arrays for holding the solution:

   gen, d, v = ts.generator(len(time), f0)
   

3. Use gen.send() to send a tuple of the next index and corresponding force vector in a loop. Re-do time-steps as necessary (note that index zero cannot be redone):

   for i in range(1, len(time)):
       # Do whatever to get i'th force
       # - note: d[:, :i] and v[:, :i] are available
       gen.send((i, fi))
   

   There is a second usage of gen.send(): if the index is negative, the force is treated as an addon to forces already included for the i’th step. This is for efficiency and only does the necessary calculations. This feature was originally written for running Henkel-Mar simulations, where interface forces are computed after computing the solution with all the other forces applied. There may be other similar situations. To demonstrate this usage:

   for i in range(1, len(time)):
       # Do whatever to get i'th force
       # - note: d[:, :i] and v[:, :i] are available
       gen.send((i, fi))
       # Do more calculations to compute an addon
       # force. Then, update the i'th solution:
       gen.send((-1, fi_addon))
   

   The class instance will have the attributes _d, _v (same objects as d and v), _a, and _force. d, v and ts._force are updated on every gen.send(). (ts._a is not used until step 4.)

4. Call ts.finalize() to get final solution in standard form:

   sol = ts.finalize()
   

   The internal references _d, _v, _a, and _force are deleted.

The generator solver currently has these limitations:

1. Unlike the normal solver, equations cannot be interspersed. That is, each type of equation (rigid-body, elastic, residual-flexibility) must be contained in a contiguous group (so that self.slices is True).

2. Unlike the normal solver, the pre_eig option is not available.

3. The first time step cannot be redone.

Examples

Set up some equations and solve both the normal way and via the generator:

>>> from pyyeti import ode
>>> import numpy as np
>>> m = np.array([10., 30., 30., 30.])    # diag mass
>>> k = np.array([0., 6.e5, 6.e5, 6.e5])  # diag stiffness
>>> zeta = np.array([0., .05, 1., 2.])    # % damping
>>> b = 2.*zeta*np.sqrt(k/m)*m            # diag damping
>>> h = 0.001                             # time step
>>> t = np.arange(0, .3001, h)            # time vector
>>> c = 2*np.pi
>>> f = np.vstack((3*(1-np.cos(c*2*t)),   # ffn
...                4.5*(np.cos(np.sqrt(k[1]/m[1])*t)),
...                4.5*(np.cos(np.sqrt(k[2]/m[2])*t)),
...                4.5*(np.cos(np.sqrt(k[3]/m[3])*t))))
>>> f *= 1.e4
>>> ts = ode.SolveExp2(m, b, k, h)

Solve the normal way:

>>> sol = ts.tsolve(f, static_ic=1)

Solve via the generator:

>>> nt = f.shape[1]
>>> gen, d, v = ts.generator(nt, f[:, 0], static_ic=1)
>>> for i in range(1, nt):
...     # Could do stuff here using d[:, :i] & v[:, :i] to
...     # get next force
...     gen.send((i, f[:, i]))
>>> sol2 = ts.finalize()

Confirm the solutions are the same:

>>> np.allclose(sol2.a, sol.a)
True
>>> np.allclose(sol2.v, sol.v)
True
>>> np.allclose(sol2.d, sol.d)
True