The Coriolis term¶

The Coriolis deflection of a horizontal velocity (u east, v north).

With the Coriolis parameter f = 2 Omega sin phi, du/dt = f v and dv/dt = -f u: to the right of the motion in the north, to the left in the south, nothing at the equator. Source: Holton & Hakim, An Introduction to Dynamic Meteorology, 2.3. Checked against a reference in numpy, which has its own tests.

Source

Holton & Hakim, An Introduction to Dynamic Meteorology, 2.3.

Proof

curve, invariant, limit.

coriolis.trait

namespace lusine.laws;

(float<m/s2>, float<m/s2>) coriolis(float<m/s> u, float<m/s> v, float<rad> latitude, float<rad/s> omega)
{
  let f = 2.0 * omega * sin(latitude);
  return (f * v, -f * u);
}

A horizontal velocity is turned by f = 2 Omega sin(latitude), to the right in the north and the left in the south.

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
EPS = float(np.finfo(np.float32).eps)

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 EPS = float(np.finfo(np.float32).eps)

The kernel¶

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

In [2]:

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

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

// lusine.laws.coriolis: the Coriolis deflection of a horizontal velocity (u east, v north). With the Coriolis
// parameter f = 2 Omega sin phi, du/dt = f v and dv/dt = -f u: to the right of the motion in the north, to the left
// in the south, nothing at the equator.
// Source: Holton & Hakim, An Introduction to Dynamic Meteorology, 2.3.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.laws;

(float<m/s2>, float<m/s2>) coriolis(float<m/s> u, float<m/s> v, float<rad> latitude, float<rad/s> omega)
{
  let f = 2.0 * omega * sin(latitude);
  return (f * v, -f * u);
}

The program¶

Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.

In [3]:

print(trait_run.law_program('laws/coriolis.trait', 'coriolis'))

print(trait_run.law_program('laws/coriolis.trait', 'coriolis'))

namespace probe;

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

struct item {
  host float<m/s> u;
  host float<m/s> v;
  host float<rad> latitude;
  host float<rad/s> omega;
  float<m/s2> result_0;
  float<m/s2> result_1;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { (items.result_0, items.result_1) = lusine.laws.coriolis(items.u, items.v, items.latitude, items.omega); }

Curve: the deflection grows as sin(latitude)¶

The magnitude for a 10 m/s wind against 2 Omega |sin(phi)| u, with the Earth's sidereal rotation rate, within binary32's rounding of the library's sine.

In [4]:

OMEGA = 7.2921159e-5
lat = np.linspace(-np.pi / 2, np.pi / 2, 37)
du, dv = trait_run.law('laws/coriolis.trait', 'coriolis', u=10.0, v=0.0, latitude=lat, omega=OMEGA)
mag = np.hypot(du, dv)
want = 2 * OMEGA * np.abs(np.sin(lat)) * 10.0
assert np.allclose(mag, want, rtol=8 * EPS, atol=8 * EPS * 2 * OMEGA * 10.0)
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(np.degrees(lat), want, label='2 Omega sin(phi) u')
ax.plot(np.degrees(lat), mag, 'o', ms=3, label='coriolis')
ax.set_xlabel('latitude (deg)')
ax.legend()
plt.show()

OMEGA = 7.2921159e-5 lat = np.linspace(-np.pi / 2, np.pi / 2, 37) du, dv = trait_run.law('laws/coriolis.trait', 'coriolis', u=10.0, v=0.0, latitude=lat, omega=OMEGA) mag = np.hypot(du, dv) want = 2 * OMEGA * np.abs(np.sin(lat)) * 10.0 assert np.allclose(mag, want, rtol=8 * EPS, atol=8 * EPS * 2 * OMEGA * 10.0) fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(np.degrees(lat), want, label='2 Omega sin(phi) u') ax.plot(np.degrees(lat), mag, 'o', ms=3, label='coriolis') ax.set_xlabel('latitude (deg)') ax.legend() plt.show()

Invariant: it turns, it does not speed up¶

The acceleration is perpendicular to the velocity, so it does no work: u du + v dv = 0 for any wind, to the rounding of the two products.

In [5]:

rng = np.random.default_rng(1)
u, v, p = rng.normal(size=100) * 20, rng.normal(size=100) * 20, rng.uniform(-1.5, 1.5, 100)
du, dv = trait_run.law('laws/coriolis.trait', 'coriolis', u=u, v=v, latitude=p, omega=7.2921159e-5)
u, v = np.float32(u).astype(float), np.float32(v).astype(float)
f = 2 * 7.2921159e-5 * np.abs(np.sin(p))
assert np.all(np.abs(u * du + v * dv) <= 4 * EPS * f * (u * u + v * v))
print('no work done')

rng = np.random.default_rng(1) u, v, p = rng.normal(size=100) * 20, rng.normal(size=100) * 20, rng.uniform(-1.5, 1.5, 100) du, dv = trait_run.law('laws/coriolis.trait', 'coriolis', u=u, v=v, latitude=p, omega=7.2921159e-5) u, v = np.float32(u).astype(float), np.float32(v).astype(float) f = 2 * 7.2921159e-5 * np.abs(np.sin(p)) assert np.all(np.abs(u * du + v * dv) <= 4 * EPS * f * (u * u + v * v)) print('no work done')

no work done

Limit: nothing at the equator; to the right in the north¶

At the equator the horizontal term vanishes; an eastward wind turns south in the north, north in the south.

In [6]:

du, dv = trait_run.law('laws/coriolis.trait', 'coriolis', u=10.0, v=0.0, latitude=[0.0, np.pi / 4, -np.pi / 4],
                       omega=7.2921159e-5)
assert dv[0] == 0.0 and dv[1] < 0 < dv[2], dv
print('north:', dv[1], 'south:', dv[2])

du, dv = trait_run.law('laws/coriolis.trait', 'coriolis', u=10.0, v=0.0, latitude=[0.0, np.pi / 4, -np.pi / 4], omega=7.2921159e-5) assert dv[0] == 0.0 and dv[1] < 0 < dv[2], dv print('north:', dv[1], 'south:', dv[2])

north: -0.0010312609374523163 south: 0.0010312609374523163