Gravity and linear drag¶

Gravity and a linear (Stokes) drag, the acceleration a body takes per unit mass, the force law an integrator steps (lusine.integrator.kick_drift and kick).

Zero at the terminal velocity g / b. A body whose drag is stiff (b dt above 1) wants the implicit step, integrator.drag_gravity, a law of its own. Written once over a vector of any length, componentwise: a binding's first vector argument fixes the length, a single value the length 1. Source: Kundu & Cohen, Fluid Mechanics, 9.12; Taylor, Classical Mechanics, 2.2. Checked against a reference in numpy, which has its own tests.

Source

Kundu & Cohen, Fluid Mechanics, 9.12 (Stokes drag); Taylor, Classical Mechanics, 2.2 (linear air resistance).

Proof

curve, limit, invariant.

drag.trait

namespace lusine.laws;

vector<float<m/s2>, _> drag(vector<float<m/s>, _> velocity, vector<float<m/s2>, _> gravity, float<s^-1> drag_rate)
{
  return gravity - drag_rate * velocity;
}

Uses, without proving it here: integrator.kick_drift, integrator.kick.

The acceleration a body takes under gravity and a drag opposing its motion. It is written once over a vector of any length: a binding over single values takes its scalar form, drag@1, and one over a plane's velocities drag@2 (as standalone functions, drag_1 to drag_3). Stepped by the library's velocity Verlet (lusine.integrator.kick_drift, this force, kick), it falls to its terminal velocity.

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

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

The kernel¶

The law as the library states it, and as this page runs it.

In [2]:

print(trait_run.source('laws/drag.trait'))

print(trait_run.source('laws/drag.trait'))

// lusine.laws.drag: gravity and a linear (Stokes) drag, the acceleration a body takes per unit mass, the force law an
// integrator steps (lusine.integrator.kick_drift and kick). Zero at the terminal velocity g / b. A body whose drag is
// stiff (b dt above 1) wants the implicit step, integrator.drag_gravity, a law of its own. Written once over a vector
// of any length, componentwise: a binding's first vector argument fixes the length, a single value
// the length 1.
// Source: Kundu & Cohen, Fluid Mechanics, 9.12; Taylor, Classical Mechanics, 2.2.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.laws;

vector<float<m/s2>, _> drag(vector<float<m/s>, _> velocity, vector<float<m/s2>, _> gravity, float<s^-1> drag_rate)
{
  return gravity - drag_rate * velocity;
}

The program¶

A body on a one-item lattice, the host's gravity, drag rate and step, and the procedures that bind the library's kernels: start takes the acceleration from the body's velocity, run once by name before the first step, and main, a step, one velocity Verlet step (the host runs any procedure of the program by name, and counts only main's runs as steps).

In [3]:

FALL = '''namespace fall;

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

struct body { float<m> position; float<m/s> velocity; float<m/s2> acceleration; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<m/s2> gravity; float<s^-1> drag_rate; float<s> dt; } settings;

void start() { bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate); }

void main()
{
  (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity,
                                                                    bodies.acceleration, settings.dt);
  bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''

def fall(gravity, drag_rate, dt, steps):
    run = trait_run.Program(['laws/drag.trait', 'integrator/kick.trait'], FALL, {'fall.bodies': [1]})
    for name, value in (('gravity', gravity), ('drag_rate', drag_rate), ('dt', dt)):
        run.store('fall.settings', name, value)
    run.step('fall.start')
    return run.steps('fall.main', steps, [('fall.bodies', 'velocity')])[0].ravel()

FALL = '''namespace fall; lusine.version version = { { number = 1 } }; struct body { float<m> position; float<m/s> velocity; float<m/s2> acceleration; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { float<m/s2> gravity; float<s^-1> drag_rate; float<s> dt; } settings; void start() { bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate); } void main() { (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity, bodies.acceleration, settings.dt); bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' def fall(gravity, drag_rate, dt, steps): run = trait_run.Program(['laws/drag.trait', 'integrator/kick.trait'], FALL, {'fall.bodies': [1]}) for name, value in (('gravity', gravity), ('drag_rate', drag_rate), ('dt', dt)): run.store('fall.settings', name, value) run.step('fall.start') return run.steps('fall.main', steps, [('fall.bodies', 'velocity')])[0].ravel()

Curve: the fall to the terminal velocity¶

From rest, the speed follows v(t) = (g / b)(1 - e^(-b t)): within 0.2% of the terminal speed at a hundredth of a second, the error falling as the step does.

In [4]:

g, b = 9.81, 0.5
ts = 0.01 * np.arange(1, 1001)
want = g / b * (1 - np.exp(-b * ts))
coarse, fine = fall(g, b, 0.01, 1000), fall(g, b, 0.005, 2000)[1::2]
e1, e2 = np.max(np.abs(coarse - want)), np.max(np.abs(fine - want))
print(f'worst {e1:.4f} m/s at dt 0.01, {e2:.4f} at 0.005, of a terminal {g / b:.2f}')
assert e1 < 0.002 * g / b and e1 / e2 > 1.8
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(ts, coarse, label='stepped by the emitted code')
ax.plot(ts, want, '--', label='(g / b)(1 - exp(-b t))')
ax.set_xlabel('t (s)')
ax.set_ylabel('v (m/s)')
ax.legend()
plt.show()

g, b = 9.81, 0.5 ts = 0.01 * np.arange(1, 1001) want = g / b * (1 - np.exp(-b * ts)) coarse, fine = fall(g, b, 0.01, 1000), fall(g, b, 0.005, 2000)[1::2] e1, e2 = np.max(np.abs(coarse - want)), np.max(np.abs(fine - want)) print(f'worst {e1:.4f} m/s at dt 0.01, {e2:.4f} at 0.005, of a terminal {g / b:.2f}') assert e1 < 0.002 * g / b and e1 / e2 > 1.8 fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(ts, coarse, label='stepped by the emitted code') ax.plot(ts, want, '--', label='(g / b)(1 - exp(-b t))') ax.set_xlabel('t (s)') ax.set_ylabel('v (m/s)') ax.legend() plt.show()

worst 0.0180 m/s at dt 0.01, 0.0090 at 0.005, of a terminal 19.62

The law at given speeds¶

The kernel alone, over a lattice of bodies whose speeds the host sets.

In [5]:

PROBE = '''namespace probe;

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

struct body { host float<m/s> velocity; float<m/s2> acceleration; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<m/s2> gravity; float<s^-1> drag_rate; } settings;

void main() { bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate); }
'''

def accelerations(velocities, gravity, drag_rate):
    run = trait_run.Program(['laws/drag.trait'], PROBE, {'probe.bodies': [len(velocities)]})
    run.store('probe.bodies', 'velocity', velocities)
    run.store('probe.settings', 'gravity', gravity)
    run.store('probe.settings', 'drag_rate', drag_rate)
    return run.steps('probe.main', 1, [('probe.bodies', 'acceleration')])[0][0]

PROBE = '''namespace probe; lusine.version version = { { number = 1 } }; struct body { host float<m/s> velocity; float<m/s2> acceleration; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { float<m/s2> gravity; float<s^-1> drag_rate; } settings; void main() { bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate); } ''' def accelerations(velocities, gravity, drag_rate): run = trait_run.Program(['laws/drag.trait'], PROBE, {'probe.bodies': [len(velocities)]}) run.store('probe.bodies', 'velocity', velocities) run.store('probe.settings', 'gravity', gravity) run.store('probe.settings', 'drag_rate', drag_rate) return run.steps('probe.main', 1, [('probe.bodies', 'acceleration')])[0][0]

Limit: no drag is free fall¶

With b = 0 the acceleration is gravity's, whatever the speed.

In [6]:

a = accelerations([-30.0, 0.0, 12.5, 400.0], 9.81, 0.0)
assert np.all(a == np.float32(9.81)), a
print('free fall:', a)

a = accelerations([-30.0, 0.0, 12.5, 400.0], 9.81, 0.0) assert np.all(a == np.float32(9.81)), a print('free fall:', a)

free fall: [9.81000042 9.81000042 9.81000042 9.81000042]

Invariant: the terminal velocity is a rest¶

At v = g / b the acceleration is zero; slower the body speeds up, faster it slows.

In [7]:

g, b = 9.81, 0.5
a = accelerations([0.5 * g / b, g / b, 2 * g / b], g, b)
assert a[0] > 0 and abs(a[1]) < 1e-5 and a[2] < 0, a
print('at rest at g / b:', a)

g, b = 9.81, 0.5 a = accelerations([0.5 * g / b, g / b, 2 * g / b], g, b) assert a[0] > 0 and abs(a[1]) < 1e-5 and a[2] < 0, a print('at rest at g / b:', a)

at rest at g / b: [ 4.90500021  0.         -9.81000042]

Curve: a throw in a plane¶

A body thrown level at 12 m/s, its position, velocity, acceleration and gravity 2-vectors, so the program binds the instances at 2 (drag@2, kick_drift@2, kick@2): its level speed decays as u0 e^(−b t) and its fall approaches the terminal −g / b as (−g / b)(1 − e^(−b t)), each axis its own law, within 0.2% of the terminal speed at a hundredth of a second.

In [8]:

PLANE = '''namespace plane;

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

struct body { vector<float<m>, 2> position; vector<float<m/s>, 2> velocity; vector<float<m/s2>, 2> acceleration; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { vector<float<m/s2>, 2> gravity; float<s^-1> drag_rate; float<s> dt; } settings;

void start() { bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate); }

void main()
{
  (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity,
                                                                    bodies.acceleration, settings.dt);
  bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''
g, b, u0 = 9.81, 0.5, 12.0
run = trait_run.Program(['laws/drag.trait', 'integrator/kick.trait'], PLANE, {'plane.bodies': [1]})
run.store('plane.bodies', 'velocity', [[u0, 0.0]])
run.store('plane.settings', 'gravity', [0.0, -g])
run.store('plane.settings', 'drag_rate', b)
run.store('plane.settings', 'dt', 0.01)
run.step('plane.start')
v = run.steps('plane.main', 1000, [('plane.bodies', 'velocity')])[0][:, 0]
ts = 0.01 * np.arange(1, 1001)
level, fall_ = u0 * np.exp(-b * ts), -g / b * (1 - np.exp(-b * ts))
assert np.max(np.abs(v[:, 0] - level)) < 0.002 * g / b and np.max(np.abs(v[:, 1] - fall_)) < 0.002 * g / b
print(f'after 10 s: level {v[-1, 0]:.4f} m/s (u0 e^(-bt) = {level[-1]:.4f}), falling {v[-1, 1]:.4f} (= {fall_[-1]:.4f})')
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(ts, v[:, 0], label='level, stepped')
ax.plot(ts, v[:, 1], label='falling, stepped')
ax.plot(ts, level, 'k--', lw=0.8)
ax.plot(ts, fall_, 'k--', lw=0.8, label='the laws')
ax.set_xlabel('t (s)')
ax.set_ylabel('v (m/s)')
ax.legend()
plt.show()

PLANE = '''namespace plane; lusine.version version = { { number = 1 } }; struct body { vector<float<m>, 2> position; vector<float<m/s>, 2> velocity; vector<float<m/s2>, 2> acceleration; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { vector<float<m/s2>, 2> gravity; float<s^-1> drag_rate; float<s> dt; } settings; void start() { bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate); } void main() { (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity, bodies.acceleration, settings.dt); bodies.acceleration = lusine.laws.drag(bodies.velocity, settings.gravity, settings.drag_rate); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' g, b, u0 = 9.81, 0.5, 12.0 run = trait_run.Program(['laws/drag.trait', 'integrator/kick.trait'], PLANE, {'plane.bodies': [1]}) run.store('plane.bodies', 'velocity', [[u0, 0.0]]) run.store('plane.settings', 'gravity', [0.0, -g]) run.store('plane.settings', 'drag_rate', b) run.store('plane.settings', 'dt', 0.01) run.step('plane.start') v = run.steps('plane.main', 1000, [('plane.bodies', 'velocity')])[0][:, 0] ts = 0.01 * np.arange(1, 1001) level, fall_ = u0 * np.exp(-b * ts), -g / b * (1 - np.exp(-b * ts)) assert np.max(np.abs(v[:, 0] - level)) < 0.002 * g / b and np.max(np.abs(v[:, 1] - fall_)) < 0.002 * g / b print(f'after 10 s: level {v[-1, 0]:.4f} m/s (u0 e^(-bt) = {level[-1]:.4f}), falling {v[-1, 1]:.4f} (= {fall_[-1]:.4f})') fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(ts, v[:, 0], label='level, stepped') ax.plot(ts, v[:, 1], label='falling, stepped') ax.plot(ts, level, 'k--', lw=0.8) ax.plot(ts, fall_, 'k--', lw=0.8, label='the laws') ax.set_xlabel('t (s)') ax.set_ylabel('v (m/s)') ax.legend() plt.show()

after 10 s: level 0.0798 m/s (u0 e^(-bt) = 0.0809), falling -19.4894 (= -19.4878)