Physics

Simulation — `simulate`

Declare params, states, laws, a scheme and a stop condition. The build writes the C# class that runs it — there is no loop in the source.

A simulate declaration in a .hph file says what a physical system is. The build generates a class whose Run() returns the trajectory and whose Stream() yields it sample by sample. Every quantity is dimension-checked, every unit is converted to SI base measures, and the only loop is the Physis driver.

open Physis.Constants
open Physis.Schemes

simulate Projectile {
  param v₀ : Velocity = 30 m/s
  param α : Angle = 45 deg
  param k : Frequency = 0.1 Hz
  state x : Length[2] = [0 m, 0 m]
  state v : Velocity[2] = v₀·[cos(α), sin(α)]
  law x' = v
  law v' = [0 m/s², −g₀] − k·v
  integrate rk45 tolerance 1e-9
  record every 10 ms
  observe height ≔ x[1]
  observe speed ≔ norm(v)
  stop when x[1] < 0 m ∧ t > 0 s
  stop when t ≥ 60 s
  run Steep with α = 70 deg
}

var projectile = new Projectile();
var run = projectile.Run();                 // or Run(projectile.Defaults with { V0 = 50 })
run.Reason                                  // Stop { Kind = When, Clause = "x[1] < 0 m ∧ t > 0 s", Line = 14 }
run.Time; run.Height; run.Speed             // double[] — one column per observe, SI base measures
run.Final.State.X                           // double[2]
run.ToFrame(); run.ToCsv();
await foreach (var s in projectile.Stream()) { /* s.T, s.Height, s.State */ }
projectile.Steep();                         // a named run

The clauses

clause meaning
param p : Q = value a run input with a constant default; a member of Params
state y : Q = value a state and its initial value (may read params); Length[2] is two components
law y' = e a first-order law; law y'' = e is second order and needs state y' : … = … for the initial velocity
integrate scheme [step Δ] [tolerance ε] [limit n] the scheme; default rk45 tolerance 1e-8
record every Δ samples on a time grid; without it, every accepted step is a sample
observe o ≔ e a recorded column (one number: index with x[i], or take norm(x), dot(a, b))
stop when c a condition that ends the run, located exactly; comparisons joined by ∧ and ∨
run Name with p = value, … a named scenario, generated as a method

In every expression, t is the time, x[i] a component, and y' the derivative of a state: the velocity of a second-order state, or the right side of a first-order law.

Motion from a Lagrangian or a Hamiltonian

simulate CartPendulum {
  param M : Mass = 2 kg
  param m : Mass = 1 kg
  param ℓ : Length = 1 m
  state s : Length = 0 m
  state s' : Velocity = 0 m/s
  state φ : Angle = 0.5 rad
  state φ' : AngularVelocity = 0 rad/s
  lagrangian ≔ 0.5·(M + m)·s'² + m·ℓ·s'·φ'·cos(φ) + 0.5·m·ℓ²·φ'² + m·g₀·ℓ·cos(φ)
  integrate rk4 step 1 ms
  observe energy ≔ hamiltonian
  conserve energy within 1e-6
  stop when t ≥ 5 s
}
  • lagrangian ≔ L moves every state q that has a velocity state q'. The build takes p = ∂L/∂q̇, the mass matrix M = ∂p/∂q̇ and the force b = ∂L/∂q − (∂p/∂q)·q̇ − ∂p/∂t with the laws of Prelude/Calculus.hb, and the run solves M·a = b (directly when M is diagonal, by elimination otherwise). L is an Energy.
  • hamiltonian ≔ H moves every pair declared state p : Momentum conjugate q = …: q' = ∂H/∂p, p' = −∂H/∂q.
  • hamiltonian read in an expression is the energy: H, or Σ q̇·∂L/∂q̇ − L.
  • conserve o within ε ends the run as Violated when the observed column o drifts past ε·max(|o₀|, 1).

Effort and flow: circuits, motors and loudspeakers

A circuit is its energies. An inductor is the kinetic part of a Lagrangian, a capacitor the potential part, a resistor Rayleigh’s dissipation, and a source or a coupling an effort on a coordinate — the same four words for a mass, a spring, a damper and a force.

simulate DcMotor {
  param V : Voltage = 12 V
  param R : Resistance = 1 Ω
  param L : Inductance = 0.5 mH
  param J : MomentOfInertia = 0.01 kg·m²
  param b : RotationalDrag = 0.001 N·m·s
  param k : TorqueConstant = 0.05 N·m/A
  state q : Charge = 0 C
  state q' : Current = 0 A
  state θ : Angle = 0 rad
  state θ' : AngularVelocity = 0 rad/s
  lagrangian ≔ 0.5·L·q'² + 0.5·J·θ'²
  dissipation ≔ 0.5·R·q'² + 0.5·b·θ'²
  effort q ≔ V − k·θ'
  effort θ ≔ k·q'
  integrate rk45 tolerance 1e-9
  observe speed ≔ θ'
  stop when t ≥ 60 s
}
  • dissipation ≔ D is a Power, a function of the velocities; the force on q gains −∂D/∂q̇. It needs a lagrangian.
  • effort q ≔ Q is the generalised force on the coordinate q, an Energy per [q]: a Voltage on a Charge, a Force on a Length, a Torque on an Angle. It adds to the Euler–Lagrange force, or to p' under a hamiltonian.
  • The winding’s k·θ' and the shaft’s k·q' are one constant read both ways — the reciprocity that makes a motor a generator. The run settles at ω = kV/(k² + Rb), which Physis.Electromechanics.motorSpeed also states.

Rigid bodies

simulate FreeBody {
  param I : MomentOfInertia[3] = [1 kg·m², 1 kg·m², 2 kg·m²]
  state ω : AngularVelocity[3] = [1 rad/s, 0 rad/s, 2 rad/s]
  law ω' = cross(I·ω, ω) / I
  integrate rk45 tolerance 1e-11
  observe energy ≔ 0.5·dot(I·ω, ω)
  conserve energy within 1e-8
  stop when t ≥ 10 s
}
  • cross(a, b) takes two quantities of size 3. Euler’s equations in the body frame are one law, I·ω being the component product.
  • An attitude is a Dimensionless[4] quaternion state whose law is ½·q ⊗ (0, ω) written out component by component; a heavy top is a lagrangian in Euler angles (samples Physics/RigidBody.hph).

Fields on grids

simulate HeatedRod {
  param α : Diffusivity = 0.01 m²/s
  param s : HeatingRate = 1 K/s
  grid (x) : Length[51] from 0 m to 1 m
  field u : Temperature = 0 K
  boundary u dirichlet 0 K
  law u' = α·laplacian(u) + s
  integrate backwardEuler step 10 s
  observe middle ≔ u[25]
  stop when maxAbs(u') < 1e-6 K/s
}
  • grid (x) : Q[N] from a to b or grid (x, y) : Q[N, M] from a to b; inside a field’s value x and y are the node’s coordinates.
  • field u : Q = e — one component per node. boundary u dirichlet v, neumann (the default) or periodic.
  • dx(u), dy(u), laplacian(u) are central-difference stencils; sum, mean, maxAbs reduce a field. Two quantities of one size multiply component by component; the inner product is dot(a, b).
  • dxBackward(u), dxForward(u), dyBackward(u), dyForward(u) are one-sided differences. Advection at a positive speed is upwind with dxBackward — monotone and conservative across a jump — where the central dx is second order but ripples. The side is yours to choose; no speed’s sign is read.
  • A steady state is reached by relaxing to rest: stop when maxAbs(u') < … — which is how a Poisson problem is solved.

Events, noise and ensembles

simulate BouncingBall {
  param e : Dimensionless = 0.8
  state y : Length = 10 m
  state y' : Velocity = 0 m/s
  law y'' = −g₀
  on y < 0 m ∧ y' < 0 m/s do y' ≔ −e·y'
  integrate rk45 tolerance 1e-12
  stop when t ≥ 6 s
}

simulate OrnsteinUhlenbeck {
  param γ : Frequency = 1 Hz
  param σ : Velocity = 1 m/s
  state v : Velocity = 0 m/s
  noise w
  law v' = −γ·v + σ·√γ·w
  integrate eulerMaruyama step 10 ms
  stop when t ≥ 10 s
}
  • on c do s ≔ e, … fires when c turns true, located on the step’s interpolant; every reset reads the state from before the event, and the run continues.
  • noise w is white noise of dimension Time^-1/2. A law must be linear in it; the build splits it into drift and diffusion, and a stochastic scheme reads them. integrate … seed n and Run(seed: n) fix the path.
  • Ensemble(runs, parameters, seed) runs many in parallel; Ensemble(sweep) one per parameter set.

Schemes

Physis.Schemes states each integrator as mathematics over the rate f(t, y) (or the acceleration a(t, x)):

scheme kind
euler first order, fixed step
midpoint second order, fixed step
rk4 fourth order, fixed step
rk45 Dormand–Prince 5(4), adaptive
verlet, leapfrog second order, symplectic — every law must be x'' = a(t, x)
backwardEuler implicit, first order, A-stable — solved by Newton
crankNicolson implicit, second order, A-stable — solved by Newton
eulerMaruyama stochastic, for laws with noise

A fixed explicit scheme declares its stability radius (scheme rk4 stable 2.78 { … }); before a run takes its first step, the spectral radius of the rate’s Jacobian times the step is checked against it, and a step past it ends the run as Unstable.

A scheme of your own is a scheme declaration: bind next (and error ≔ … for adaptive <order>), or nextX and nextV for a second scheme.

How a run proceeds

  1. The initial state is the first sample; a clause that already holds ends the run.
  2. A fixed scheme steps by step; an adaptive one retries until its error is within tolerance and scales the next step.
  3. A non-finite state ends the run as Diverged; limit accepted steps end it as Limit. Neither throws.
  4. After each step the first clause that holds is located by bisection over the cubic Hermite interpolant of the step, so a projectile stops at the ground and not one step below it.
  5. Run(step:, tolerance:) overrides the declared settings.

Diagnostics

code meaning
PH0101 a state with no law or two, a law naming no state, a second-order law with no velocity state
PH0102 a second-order scheme with a first-order law, or an acceleration that reads a velocity
PH0103 a fixed-step scheme with no step
PH0104 no stop when
PH0105 a scheme that reads what it may not or binds the wrong outputs
PH0106 a symbolic extent (Length[n])
PH0107 a scheme not in scope
PH0108 a run overriding no param
PH0109 shapes that do not fit
PH0110 a setting of the wrong dimension, or a param default that is not a constant
PH0111 a simulation named like another generated class
PH0112 a derivative the calculus laws do not take
PH0113 a grid operator on something that is not a field, or a field with no grid
PH0114 a law not linear in its noise, or noise under a scheme that is not stochastic
PH0115 a lagrangian/hamiltonian with nothing to move
PH0116 a conserve naming no observed column
PH0117 a dissipation or effort with no lagrangian (or a dissipation beside a hamiltonian), an effort on no coordinate, two efforts on one, or an effort of the wrong size