The kick-drift: velocity Verlet's first half¶

A half-kick, then the drift.

Requires

the Position x (m) and Velocity v (m/s), scalars or vectors; the Acceleration a(x) (m/s^2) at x, the game's force kernels' over the mass; the step dt (s).

Provides

the position after the step x' and the half-step velocity v_half.

Law

v_half = v + a(x) dt/2; x' = x + v_half dt. The first of velocity Verlet's three passes: the game's force kernels then give a(x') at the new position, and integrator.kick closes the step (a kernel never calls a kernel, so the force between them is the schedule's).

Source

Swope, Andersen, Berens & Wilson, J. Chem. Phys. 76, 637 (1982); Hairer, Lubich & Wanner, Geometric Numerical Integration, I.1.4.

Proof

limit, invariant.

integrator.kick_drift(position, velocity, acceleration, dt)

A half-kick from the acceleration at the position, then the drift with that half-step velocity. It is written once over a vector of any length: a binding over single values takes kick_drift@1, the scalar form, and one over 3-vectors kick_drift@3 (as standalone functions, kick_drift_1 to kick_drift_3). The method's order and its orbit are proved on the kick's page, where whole steps run.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.kernels._si import G0
EPS = float(np.finfo(np.float32).eps)   # Trait's float is binary32

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run from lusine.kernels._si import G0 EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32

The kernel¶

In [2]:

print(trait_run.source('integrator/kick.trait'))

print(trait_run.source('integrator/kick.trait'))

// lusine.integrator.kick and kick_drift: velocity Verlet's halves. kick_drift takes the first half kick and the drift,
// kick the second half kick with the acceleration at the new position. Both are written once over a vector of any
// length, componentwise: a binding's first vector argument fixes the length, a single value the
// length 1, so one pair of halves steps a rod, a plane and a space.
// Checked against references in numpy, which have their own tests.

namespace lusine.integrator;

vector<float<m/s>, _> kick(vector<float<m/s>, _> velocity, vector<float<m/s2>, _> acceleration, float<s> dt)
{
  return velocity + acceleration * (0.5 * dt);
}

(vector<float<m>, _>, vector<float<m/s>, _>) kick_drift(vector<float<m>, _> position, vector<float<m/s>, _> velocity,
                                                    vector<float<m/s2>, _> acceleration, float<s> dt)
{
  let half = velocity + acceleration * (0.5 * dt);
  return (position + half * dt, half);
}

The program¶

One step from states the host sets, a body each: a component, or a case.

In [3]:

STEP = '''namespace step;

lusine.version version = { { number = 1 } };

struct body { host float<m> position; host float<m/s> velocity; host float<m/s2> acceleration; host float<s> dt;
              float<m> moved; float<m/s> half; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);

void main()
{
  (bodies.moved, bodies.half) = lusine.integrator.kick_drift(position = bodies.position, velocity = bodies.velocity,
                                                             acceleration = bodies.acceleration, dt = bodies.dt);
}
'''

def kick_drift(position, velocity, acceleration, dt):
    n = len(position)
    run = trait_run.Program(['integrator/kick.trait'], STEP, {'step.bodies': [n]})
    for name, value in (('position', position), ('velocity', velocity), ('acceleration', acceleration), ('dt', dt)):
        run.store('step.bodies', name, np.broadcast_to(np.asarray(value, dtype=float), (n,)))
    moved, half = run.steps('step.main', 1, [('step.bodies', 'moved'), ('step.bodies', 'half')])
    return moved[0], half[0]

STEP = '''namespace step; lusine.version version = { { number = 1 } }; struct body { host float<m> position; host float<m/s> velocity; host float<m/s2> acceleration; host float<s> dt; float<m> moved; float<m/s> half; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); void main() { (bodies.moved, bodies.half) = lusine.integrator.kick_drift(position = bodies.position, velocity = bodies.velocity, acceleration = bodies.acceleration, dt = bodies.dt); } ''' def kick_drift(position, velocity, acceleration, dt): n = len(position) run = trait_run.Program(['integrator/kick.trait'], STEP, {'step.bodies': [n]}) for name, value in (('position', position), ('velocity', velocity), ('acceleration', acceleration), ('dt', dt)): run.store('step.bodies', name, np.broadcast_to(np.asarray(value, dtype=float), (n,))) moved, half = run.steps('step.main', 1, [('step.bodies', 'moved'), ('step.bodies', 'half')]) return moved[0], half[0]

Limit: no force, a drift; a constant force, a parabola¶

With no acceleration the step is the drift x + v dt; from rest under gravity it lands on -g dt²/2, for every step length.

In [4]:

x, v = np.array([0.3, -1.2]), np.array([2.0, 0.5])
x1, half = kick_drift(x, v, np.zeros(2), 0.01)
assert np.allclose(x1, x + v * 0.01, rtol=4 * EPS, atol=0) and np.array_equal(half, v.astype(np.float32))
dts = np.array([0.001, 0.01, 0.1, 0.5])
drop, _ = kick_drift(np.zeros(4), np.zeros(4), np.full(4, -G0), dts)
assert np.allclose(drop, -0.5 * G0 * dts ** 2, rtol=4 * EPS, atol=0)
fig, ax = plt.subplots(figsize=(7, 3))
ax.loglog(dts, -drop, 'o', label='kick_drift from rest')
ax.loglog(dts, 0.5 * G0 * dts ** 2, '--', label='g dt^2 / 2')
ax.set_xlabel('dt (s)')
ax.set_ylabel('fall (m)')
ax.legend()
plt.show()

x, v = np.array([0.3, -1.2]), np.array([2.0, 0.5]) x1, half = kick_drift(x, v, np.zeros(2), 0.01) assert np.allclose(x1, x + v * 0.01, rtol=4 * EPS, atol=0) and np.array_equal(half, v.astype(np.float32)) dts = np.array([0.001, 0.01, 0.1, 0.5]) drop, _ = kick_drift(np.zeros(4), np.zeros(4), np.full(4, -G0), dts) assert np.allclose(drop, -0.5 * G0 * dts ** 2, rtol=4 * EPS, atol=0) fig, ax = plt.subplots(figsize=(7, 3)) ax.loglog(dts, -drop, 'o', label='kick_drift from rest') ax.loglog(dts, 0.5 * G0 * dts ** 2, '--', label='g dt^2 / 2') ax.set_xlabel('dt (s)') ax.set_ylabel('fall (m)') ax.legend() plt.show()

Invariant: the mirror¶

The mirrored state under the mirrored force steps to the mirrored state, bit for bit: nothing in the step prefers a direction.

In [5]:

a_ = kick_drift([0.3, -0.7], [1.1, 0.2], [-0.5, 0.9], 0.013)
b_ = kick_drift([-0.3, 0.7], [-1.1, -0.2], [0.5, -0.9], 0.013)
assert np.array_equal(a_[0], -b_[0]) and np.array_equal(a_[1], -b_[1])
print('mirrored exactly')

a_ = kick_drift([0.3, -0.7], [1.1, 0.2], [-0.5, 0.9], 0.013) b_ = kick_drift([-0.3, 0.7], [-1.1, -0.2], [0.5, -0.9], 0.013) assert np.array_equal(a_[0], -b_[0]) and np.array_equal(a_[1], -b_[1]) print('mirrored exactly')

mirrored exactly

Limit: in three dimensions¶

A body's position, velocity and acceleration 3-vectors, so the program binds the instance at 3 (kick_drift@3): from rest under a constant force the drift reaches a t² / 2 in each component, and over random states every component is the scalar form's on that component alone, bit for bit.

In [6]:

SPACE = STEP.replace('namespace step;', 'namespace space;').replace(
    '''struct body { host float<m> position; host float<m/s> velocity; host float<m/s2> acceleration; host float<s> dt;
              float<m> moved; float<m/s> half; };''',
    '''struct body { host vector<float<m>, 3> position; host vector<float<m/s>, 3> velocity;
              host vector<float<m/s2>, 3> acceleration; host float<s> dt;
              vector<float<m>, 3> moved; vector<float<m/s>, 3> half; };''')


def kick_drift3(position, velocity, acceleration, dt):
    n = len(position)
    run = trait_run.Program(['integrator/kick.trait'], SPACE, {'space.bodies': [n]})
    for name, value in (('position', position), ('velocity', velocity), ('acceleration', acceleration), ('dt', dt)):
        run.store('space.bodies', name, np.asarray(value, dtype=float))
    moved, half = run.steps('space.main', 1, [('space.bodies', 'moved'), ('space.bodies', 'half')])
    return moved[0], half[0]

a = np.array([0.5, -1.0, -9.81])
x1, half = kick_drift3(np.zeros((1, 3)), np.zeros((1, 3)), [a], [0.1])
assert np.allclose(x1[0], 0.5 * a * 0.01, rtol=4 * EPS, atol=0) and np.allclose(half[0], 0.05 * a, rtol=4 * EPS, atol=0)
rng = np.random.default_rng(90)
p, v, acc = (rng.uniform(-50, 50, (64, 3)) for _ in range(3))
dt = rng.uniform(1e-3, 0.1, 64)
x3, h3 = kick_drift3(p, v, acc, dt)
for c in range(3):
    x1c, h1c = kick_drift(p[:, c], v[:, c], acc[:, c], dt)
    assert np.array_equal(x3[:, c], x1c) and np.array_equal(h3[:, c], h1c)
print(x1[0], "; 64 random 3-vector states, every component the scalar form's bit for bit")

SPACE = STEP.replace('namespace step;', 'namespace space;').replace( '''struct body { host float<m> position; host float<m/s> velocity; host float<m/s2> acceleration; host float<s> dt; float<m> moved; float<m/s> half; };''', '''struct body { host vector<float<m>, 3> position; host vector<float<m/s>, 3> velocity; host vector<float<m/s2>, 3> acceleration; host float<s> dt; vector<float<m>, 3> moved; vector<float<m/s>, 3> half; };''') def kick_drift3(position, velocity, acceleration, dt): n = len(position) run = trait_run.Program(['integrator/kick.trait'], SPACE, {'space.bodies': [n]}) for name, value in (('position', position), ('velocity', velocity), ('acceleration', acceleration), ('dt', dt)): run.store('space.bodies', name, np.asarray(value, dtype=float)) moved, half = run.steps('space.main', 1, [('space.bodies', 'moved'), ('space.bodies', 'half')]) return moved[0], half[0] a = np.array([0.5, -1.0, -9.81]) x1, half = kick_drift3(np.zeros((1, 3)), np.zeros((1, 3)), [a], [0.1]) assert np.allclose(x1[0], 0.5 * a * 0.01, rtol=4 * EPS, atol=0) and np.allclose(half[0], 0.05 * a, rtol=4 * EPS, atol=0) rng = np.random.default_rng(90) p, v, acc = (rng.uniform(-50, 50, (64, 3)) for _ in range(3)) dt = rng.uniform(1e-3, 0.1, 64) x3, h3 = kick_drift3(p, v, acc, dt) for c in range(3): x1c, h1c = kick_drift(p[:, c], v[:, c], acc[:, c], dt) assert np.array_equal(x3[:, c], x1c) and np.array_equal(h3[:, c], h1c) print(x1[0], "; 64 random 3-vector states, every component the scalar form's bit for bit")

[ 0.0025  -0.005   -0.04905] ; 64 random 3-vector states, every component the scalar form's bit for bit