Introduction

Rumoca is a Modelica compiler, simulator, and code generation toolkit written in Rust. It takes equation-based models of physical systems and turns them into simulations you can run from the command line, VS Code, or a browser — or into code for other ecosystems such as Python (SymPy, JAX, CasADi), C, Rust, and FMI.

What is Modelica?

Modelica is an open, equation-based language for modeling physical systems. Instead of writing step-by-step simulation code, you declare the equations that govern your system and let the compiler decide how to solve them.

The example below is a small hot-air-balloon model. It is live: edit the code, press ▶ Simulate to integrate it right here in your browser, or Show DAE to see the equation system the compiler produces. The editor has the same syntax highlighting, completion, and error checking as the Rumoca VS Code extension, powered by the same compiler running in WebAssembly.

model HotAirBalloon "Hot-air balloon warmed by a burner"
  parameter Real m = 320.0 "Balloon, basket, and passenger mass [kg]";
  parameter Real g = 9.81 "Gravity [m/s2]";
  parameter Real Tamb = 293.15 "Ambient air temperature [K]";
  parameter Real tau = 35.0 "Envelope cooling time constant [s]";
  parameter Real burnerHeat = 0.8 "Temperature rise from burner [K/s]";
  parameter Real liftPerKelvin = 130.0 "Buoyant lift per kelvin [N/K]";
  parameter Real drag = 240.0 "Vertical drag [N.s/m]";
  parameter Real targetHeight = 30.0 "Altitude where the burner switches off [m]";
  parameter Real fuelBurnTime = 28.0 "Fuel burn time with burner fully on [s]";
  Real h(start = 2.0) "Altitude [m]";
  Real v(start = 0.0) "Vertical speed [m/s]";
  Real T(start = 320.0) "Envelope air temperature [K]";
  Real fuel(start = 1.0) "Fuel fraction";
  Real fuelPercent "Fuel remaining [%]";
  Real burner "Burner command";
equation
  burner = if fuel > 0.0 then if h < targetHeight then 1.0 else 0.0 else 0.0;
  fuelPercent = if fuel > 0.0 then 100.0 * fuel else 0.0;
  der(fuel) = if burner > 0.0 then -1.0 / fuelBurnTime else 0.0;
  der(T) = burnerHeat * burner - (T - Tamb) / tau;
  der(h) = if h > 0.0 then v else if v > 0.0 then v else 0.0;
  m * der(v) = liftPerKelvin * (T - Tamb) - m * g - drag * v;
  annotation(experiment(StopTime = 70.0, Interval = 0.1, Solver = "rk-like"));
end HotAirBalloon;

// Render the balloon state as a small Three.js scene.
// Only the integrated states (h, v, T, fuel) are read; the algebraics
// `burner` and `fuelPercent` are recomputed from the states below so the
// viz works on every solver path (the GPU integrator returns states only
// and freezes algebraics).
api.plotSeries(['h', 'fuel']);
const { THREE } = await api.loadThree();
const h = api.series('h');
const fuel = api.series('fuel');
// burner = if fuel > 0 and h < targetHeight then 1 else 0 (model eqn);
// targetHeight matches the model parameter.
const targetHeight = 30.0;
const burner = h.map((hi, k) => (fuel[k] > 0.0 && hi < targetHeight) ? 1.0 : 0.0);

container.classList.add('rumoca-live-surface');
const host = document.createElement('div');
host.className = 'rumoca-live-surface-host';
container.appendChild(host);

const renderer = new THREE.WebGLRenderer({ antialias: true, preserveDrawingBuffer: true });
renderer.setPixelRatio(Math.min(window.devicePixelRatio || 1, 2));
renderer.setClearColor(0x8fd4ff, 1);
renderer.outputColorSpace = THREE.SRGBColorSpace;
host.appendChild(renderer.domElement);

const scene = new THREE.Scene();
scene.fog = new THREE.Fog(0x9fd5f1, 16, 42);
const camera = new THREE.PerspectiveCamera(42, 1, 0.1, 80);
camera.position.set(0, 4.2, 13.0);

scene.add(new THREE.HemisphereLight(0xfff2d4, 0x5a8a55, 1.6));
const sunLight = new THREE.DirectionalLight(0xffcc83, 2.4);
sunLight.position.set(5, 8, 4);
scene.add(sunLight);

const sun = new THREE.Mesh(
  new THREE.SphereGeometry(0.55, 32, 16),
  new THREE.MeshBasicMaterial({ color: 0xffe08a, fog: false })
);
sun.position.set(5.8, 5.2, -8);
scene.add(sun);

const ground = new THREE.Mesh(
  new THREE.PlaneGeometry(50, 35),
  new THREE.MeshStandardMaterial({ color: 0x7db35c, roughness: 0.9 })
);
ground.rotation.x = -Math.PI / 2;
ground.position.set(0, -2.4, -7);
scene.add(ground);

for (const [x, z, sx, sy, color] of [
  [-8, -13, 5.4, 1.5, 0x5f9461],
  [2, -14, 7.0, 1.9, 0x6fa76b],
  [10, -12, 4.8, 1.3, 0x8ab86e],
]) {
  const hill = new THREE.Mesh(
    new THREE.SphereGeometry(1, 24, 12),
    new THREE.MeshStandardMaterial({ color, roughness: 0.95, flatShading: true })
  );
  hill.scale.set(sx, sy, 2);
  hill.position.set(x, -2.0, z);
  scene.add(hill);
}

let randomSeed = 9;
function rand() {
  randomSeed = (randomSeed * 1664525 + 1013904223) >>> 0;
  return randomSeed / 4294967296;
}

function cloud(x, y, z, scale) {
  const group = new THREE.Group();
  const material = new THREE.MeshStandardMaterial({
    color: 0xffffff,
    roughness: 0.8,
    transparent: true,
    opacity: 0.82,
  });
  const puffCount = 5 + Math.floor(rand() * 4);
  for (let i = 0; i < puffCount; i++) {
    const puff = new THREE.Mesh(
      new THREE.SphereGeometry(0.55 + rand() * 0.55, 24, 16),
      material
    );
    puff.position.set((rand() - 0.5) * 1.7, (rand() - 0.35) * 0.45, (rand() - 0.5) * 0.35);
    puff.scale.set(1.0 + rand() * 0.9, 0.58 + rand() * 0.35, 0.82 + rand() * 0.45);
    group.add(puff);
  }
  group.position.set(x, y, z);
  group.scale.setScalar(scale);
  group.userData.speed = 0.0014 + rand() * 0.003;
  scene.add(group);
  return group;
}
const clouds = [];
for (let i = 0; i < 7; i++) {
  clouds.push(cloud(-7 + rand() * 14, 3.0 + rand() * 2.5, -4.5 - rand() * 5.5, 0.32 + rand() * 0.42));
}

function envelopeTexture() {
  const canvas = document.createElement('canvas');
  canvas.width = 512;
  canvas.height = 256;
  const ctx = canvas.getContext('2d');
  const colors = ['#2453d6', '#ffe2bd', '#2d60d8', '#ffd8aa'];
  const stripeWidth = canvas.width / 10;
  for (let i = 0; i < 10; i++) {
    const grad = ctx.createLinearGradient(i * stripeWidth, 0, (i + 1) * stripeWidth, 0);
    grad.addColorStop(0, '#173b9d');
    grad.addColorStop(0.16, colors[i % colors.length]);
    grad.addColorStop(0.58, colors[i % colors.length]);
    grad.addColorStop(1, '#173b9d');
    ctx.fillStyle = grad;
    ctx.fillRect(i * stripeWidth, 0, stripeWidth + 1, canvas.height);
  }
  ctx.globalAlpha = 0.28;
  ctx.fillStyle = '#ffffff';
  ctx.fillRect(canvas.width * 0.74, canvas.height * 0.12, 18, canvas.height * 0.52);
  ctx.globalAlpha = 1.0;
  const texture = new THREE.CanvasTexture(canvas);
  texture.colorSpace = THREE.SRGBColorSpace;
  return texture;
}

const balloon = new THREE.Group();
scene.add(balloon);

const envelopeProfile = [
  [0.14, 0.0],
  [0.44, 0.12],
  [0.72, 0.38],
  [1.08, 0.82],
  [1.55, 1.5],
  [2.02, 2.32],
  [2.3, 3.1],
  [2.24, 3.72],
  [1.9, 4.28],
  [1.3, 4.72],
  [0.62, 5.0],
  [0.06, 5.08],
].map(([radius, y]) => new THREE.Vector2(radius, y));
const envelope = new THREE.Mesh(
  new THREE.LatheGeometry(envelopeProfile, 160),
  new THREE.MeshStandardMaterial({ map: envelopeTexture(), roughness: 0.55 })
);
balloon.add(envelope);

const throat = new THREE.Mesh(
  new THREE.TorusGeometry(0.45, 0.08, 12, 40),
  new THREE.MeshStandardMaterial({ color: 0x1743b8, roughness: 0.6 })
);
throat.rotation.x = Math.PI / 2;
throat.position.y = 0.14;
balloon.add(throat);

const basketMaterial = new THREE.MeshStandardMaterial({ color: 0x8a5429, roughness: 0.8 });
const basket = new THREE.Mesh(new THREE.BoxGeometry(0.9, 0.55, 0.65), basketMaterial);
basket.position.y = -0.45;
balloon.add(basket);

const ropeMaterial = new THREE.LineBasicMaterial({ color: 0x463021 });
for (const x of [-0.38, 0.38]) {
  for (const z of [-0.28, 0.28]) {
    balloon.add(new THREE.Line(
      new THREE.BufferGeometry().setFromPoints([
        new THREE.Vector3(x, -0.16, z),
        new THREE.Vector3(x * 1.25, 0.2, z * 1.25),
      ]),
      ropeMaterial
    ));
  }
}

const flameLight = new THREE.PointLight(0xff7a1a, 0, 3.0);
flameLight.position.set(0, 0.38, 0);
balloon.add(flameLight);
const flame = new THREE.Mesh(
  new THREE.ConeGeometry(0.14, 0.45, 24),
  new THREE.MeshBasicMaterial({ color: 0xff7a1a, transparent: true, opacity: 0.85 })
);
flame.rotation.x = Math.PI;
flame.position.y = 0.38;
balloon.add(flame);

function bird(x, y, z, phase) {
  const positions = new Float32Array([
    0, 0, 0.5, 1.2, 0, 0, 0, 0, -0.5,
    0, 0, -0.5, -1.2, 0, 0, 0, 0, 0.5,
  ]);
  const geometry = new THREE.BufferGeometry();
  geometry.setAttribute('position', new THREE.BufferAttribute(positions, 3));
  geometry.computeBoundingSphere();
  const mesh = new THREE.Mesh(
    geometry,
    new THREE.MeshBasicMaterial({ color: 0xffffff, side: THREE.DoubleSide })
  );
  mesh.position.set(x, y, z);
  mesh.scale.setScalar(0.18 + rand() * 0.08);
  mesh.rotation.y = Math.PI * 0.5;
  mesh.userData = {
    homeX: x,
    homeY: y,
    phase,
    speed: 0.011 + rand() * 0.007,
    drift: 0.12 + rand() * 0.18,
    positions,
  };
  scene.add(mesh);
  return mesh;
}
const birds = [];
for (let i = 0; i < 9; i++) {
  birds.push(bird(-4.5 + rand() * 9, 3.9 + rand() * 1.5, -7.5 + rand() * 2.5, rand() * Math.PI * 2));
}

function resize() {
  const width = Math.max(320, Math.floor(host.clientWidth || 720));
  const height = Math.max(320, Math.min(520, Math.floor(width * 0.58)));
  renderer.setSize(width, height, false);
  camera.aspect = width / height;
  camera.updateProjectionMatrix();
}
new ResizeObserver(resize).observe(host);
resize();

api.addAnimation(times, (frame) => {
  const t = times[frame];
  const altitude = Math.max(0, h[frame] || 0);
  const burn = Math.max(0, Math.min(1, burner[frame] || 0));
  const fuelPercent = Math.max(0, Math.min(100, Math.round((fuel[frame] || 0) * 100)));
  balloon.position.set(0.5 * Math.sin(t * 0.16), -1.2 + altitude / 18, 0);
  balloon.rotation.z = burn * 0.05 * Math.sin(t * 0.8);
  envelope.scale.set(1 + burn * 0.02, 1 + burn * 0.03, 1 + burn * 0.02);
  flame.visible = burn > 0.12;
  flame.scale.setScalar(0.8 + burn * 0.8);
  flameLight.intensity = burn > 0.12 ? 2.0 + burn * 5.0 : 0;

  for (let i = 0; i < clouds.length; i++) {
    clouds[i].position.x += clouds[i].userData.speed;
    if (clouds[i].position.x > 7) clouds[i].position.x = -7;
  }
  for (let i = 0; i < birds.length; i++) {
    const bird = birds[i];
    bird.position.z += bird.userData.speed;
    if (bird.position.z > -2.0) bird.position.z = -8.0 - rand() * 1.5;
    bird.position.x = bird.userData.homeX + Math.sin(t * 0.7 + bird.userData.phase) * bird.userData.drift;
    bird.position.y = bird.userData.homeY + Math.sin(t * 1.1 + bird.userData.phase) * 0.08;
    const flap = Math.sin((bird.position.z + t) * 5.2 + bird.userData.phase) * 0.85;
    bird.userData.positions[4] = flap;
    bird.userData.positions[13] = flap;
    bird.geometry.attributes.position.needsUpdate = true;
  }

  camera.lookAt(balloon.position.x, balloon.position.y + 2.2, balloon.position.z);
  renderer.render(scene, camera);
  return `t = ${api.formatTick(times[frame])} s · fuel ${fuelPercent}%`;
}, 8000);

Note what you did not have to write: no integration loop, no state vector bookkeeping, no event logic. der(h) means the time derivative of h, and the compiler transforms the equations into a form a numerical solver can integrate.

What Rumoca Gives You

The Normal Workflow

1. Write or open a Modelica model (.mo file).
2. Configure any external Modelica package roots, such as the Modelica Standard Library (MSL).
3. Run a direct command (rumoca sim model.mo) or a colocated rumoca-scenario.toml scenario (rumoca sim -c rumoca-scenario.toml).
4. Inspect results in the CLI, VS Code, the browser viewer, or generated target output.

Project Status

Rumoca is in active development. It compiles and simulates a growing subset of Modelica, validated continuously against the Modelica Standard Library, but it is not yet a complete replacement for mature tools such as OpenModelica or Dymola. See Language Support Status for an honest description of what works today.

How This Book Is Organized

• Getting Started installs Rumoca and walks you through your first model.
• The Modelica Language explains equation-based modeling and what Rumoca supports.
• Tools covers the CLI, VS Code extension, playground, formatter, and linter.
• Simulation covers direct runs, scenario files, solvers, interactive simulation, and debugging.
• Code Generation covers built-in and custom targets.

Developers who want to understand or modify the compiler itself should read the companion Rumoca Dev Guide book.

Installation

Rumoca is distributed from GitHub Releases: prebuilt binaries, Python wheels, the VS Code extension, and WASM assets. It is not published to crates.io.

Prebuilt Binaries (Recommended)

Linux and macOS:

curl --proto '=https' --tlsv1.2 -LsSf https://raw.githubusercontent.com/cognipilot/rumoca/main/infra/install/install.sh | bash

Install a specific version, and optionally the rumoca-lsp language server:

curl --proto '=https' --tlsv1.2 -LsSf https://raw.githubusercontent.com/cognipilot/rumoca/main/infra/install/install.sh | bash -s -- --version v0.8.0 --with-lsp

Windows PowerShell uses infra/install/install.ps1 from the same directory.

The installer places binaries in ~/.local/bin by default; override with --bin-dir <path>. Check the result with:

rumoca --version

VS Code Extension

Install Rumoca Modelica from the VS Code marketplace. The extension bundles its own rumoca-lsp server, so it works without a separate compiler install. See VS Code Extension.

Python Package

pip install rumoca

The Python package exposes the compiler for scripting and notebook use.

Shell Completions

rumoca completions bash > ~/.local/share/bash-completion/completions/rumoca

Other shells (zsh, fish, …) are supported; run rumoca completions --help.

From Source

Install the Rust toolchain selected by rust-toolchain.toml, then build the workspace:

git clone https://github.com/CogniPilot/rumoca
cd rumoca
cargo build --workspace

For interactive simulation, prefer release builds:

cargo run -p rumoca --release -- --help

The repository includes an xtask developer CLI used by CI and local development (cargo xtask verify quick, cargo xtask vscode test, …). It is documented in the Rumoca Dev Guide book and CONTRIBUTING.md.

Modelica Library Dependencies

The repository examples use pinned Modelica dependencies declared in examples/modelica_dependencies.toml. Fetch them with:

cargo xtask repo modelica-deps ensure

This downloads the Modelica Standard Library (MSL) and the CogniPilot Modelica Models (CMM) into target/. The repository’s committed VS Code settings point at those directories with workspace-relative paths, so the examples work as soon as the download finishes. See Using Modelica Libraries for how library lookup works in general.

Quick Start

This page assumes rumoca is on your PATH (Installation). If you are working from a source checkout, replace rumoca with cargo run -p rumoca --release --.

Simulate a Model Directly

Save this as Ball.mo:

model Ball
  Real x(start=10);
  Real v(start=1);
  parameter Real g = 9.81;
equation
  der(x) = v;
  der(v) = -g;
  when x < 0 then
    reinit(v, -0.8*pre(v));
  end when;
end Ball;

Then run:

rumoca sim Ball.mo --t-end 10

Rumoca compiles the model, simulates to t = 10 s, and writes an HTML report (Ball_results.html by default) with interactive plots of every variable. The model name is inferred from the file; pass --model to pick a specific class.

Useful options:

rumoca sim Ball.mo --model Ball --t-end 10 --solver rk-like --dt 0.01 -o ball.html

Use --source-root for packages that are not in the same source tree:

rumoca sim my_model.mo \
  --model MyPackage.MyModel \
  --source-root target/msl/ModelicaStandardLibrary-4.1.0

Run a Scenario

The preferred repeatable workflow is a colocated rumoca-scenario.toml scenario file that records the model, solver, plots, and viewer settings in one place:

rumoca sim -c examples/simulation/rumoca-scenario.ball.toml

Interactive examples use the same command shape:

rumoca sim -c examples/interactive/quadrotor/rumoca-scenario.acro.toml

Generate a commented starter scenario and validate it without running:

rumoca sim init > rumoca-scenario.toml
rumoca sim check -c rumoca-scenario.toml

See Scenario Files for the full format.

Compile or Generate Code

List the built-in code generation targets:

rumoca targets

Render a target:

rumoca compile examples/models/SympyDecay.mo \
  --model SympyDecay \
  --target c-ode \
  --output /tmp/decay_c_ode

Dump an intermediate representation of the compiler instead:

rumoca compile Ball.mo --emit dae-mo     # the DAE as Modelica source
rumoca compile Ball.mo --emit solve-json # the solver IR as JSON

Next Steps

• Build a model from scratch in Your First Model.
• Learn the scenario format in Scenario Files.
• Explore the full CLI reference.

Your First Model

This tutorial builds a small physical model from scratch and explains each piece of Modelica syntax as it appears. The code blocks are live: you can edit and simulate them directly on this page.

A Cooling Cup of Coffee

Newton’s law of cooling says the temperature T of an object approaches the ambient temperature T_amb at a rate proportional to their difference:

model Coffee "Newton's law of cooling"
  parameter Real T_amb = 20.0 "Ambient temperature [degC]";
  parameter Real tau = 300.0 "Cooling time constant [s]";
  Real T(start = 90.0) "Coffee temperature [degC]";
equation
  der(T) = (T_amb - T) / tau;
  annotation(experiment(StopTime = 1800.0));
end Coffee;

Press ▶ Simulate and watch the temperature decay toward 20 °C.

Reading the model line by line:

• model Coffee "..." — declares a class named Coffee. The string after the name is a description; tools display it but it has no effect on the equations.
• parameter Real T_amb = 20.0 — a parameter is fixed during a simulation but adjustable between runs. Try changing it and re-running.
• Real T(start = 90.0) — a continuous variable. start gives its initial value. Because der(T) appears in an equation, T is a state: the solver integrates it through time.
• der(T) = (T_amb - T) / tau — an equation, not an assignment. You could equally write (T_amb - T) / tau = der(T); the compiler decides how to solve the system.
• annotation(experiment(StopTime = 1800.0)) — the standard Modelica way to store default simulation settings with the model. The live editors on these pages, the web playground, and direct CLI runs honor it; --t-end or the [sim] section of a scenario file overrides it.

Adding a Second State

Models grow by adding variables and equations — one equation per unknown. Here is the same idea applied to a mass hanging on a spring, which needs two states (position and velocity):

model HangingMass
  parameter Real m = 0.5 "Mass [kg]";
  parameter Real k = 20.0 "Spring constant [N/m]";
  parameter Real c = 0.5 "Damping [N.s/m]";
  parameter Real g = 9.81;
  Real x(start = 0.0) "Displacement below natural length [m]";
  Real v(start = 0.0);
  Real f_spring "Spring force [N]";
equation
  f_spring = -k * x;
  der(x) = v;
  m * der(v) = f_spring - c * v + m * g;
  annotation(experiment(StopTime = 5.0));
end HangingMass;

f_spring is an algebraic variable: it has no der(), so the compiler solves it from its equation at every step instead of integrating it. Mixing differential and algebraic equations like this is what makes the system a DAE (differential-algebraic equation system) — press Show DAE to see the sorted system Rumoca produces.

Running It Natively

Save the model as HangingMass.mo and run:

rumoca sim HangingMass.mo --t-end 5

The HTML report HangingMass_results.html plots all variables.

Recording the Run as a Scenario

Once a model has settings worth keeping — solver, plots, output paths — record them in a rumoca-scenario.toml scenario next to the model:

[rumoca]
version = "1"
task = "simulate"

[model]
file = "HangingMass.mo"
name = "HangingMass"

[sim]
t_end = 5.0
solver = "auto"
output = "hanging_mass.html"

rumoca sim -c rumoca-scenario.toml

From here:

• Modeling with Equations explains the language concepts systematically.
• Events and Discrete Behavior adds switching, impacts, and sampled control.
• Scenario Files documents everything a rumoca-scenario.toml can do.

Modeling with Equations

This chapter explains the core ideas of equation-based modeling as Rumoca implements them. It is not a full Modelica tutorial — the Modelica Language Specification is the authority — but it covers what you need to write and read most models.

Equations, Not Assignments

A Modelica equation declares a relationship that must hold at every instant. The compiler — not you — decides which variable each equation solves for, in what order, and which variables become integrator states. This is why you can write m * der(v) = f instead of der(v) := f / m, and why models stay readable as they grow: adding a component adds equations, and the compiler re-sorts the whole system.

A model is balanced when it has exactly one equation per unknown. Rumoca checks this during compilation and reports counts when they disagree.

Variability

Every variable has a variability that tells the compiler how it can change:

States and Initial Values

A continuous variable that appears under der(...) becomes a state. Give states an initial value with the start attribute:

Real x(start = 1.0);

Variables without der() are algebraic: solved from the equation system at each step rather than integrated.

A Worked Example

The Van der Pol oscillator shows a nonlinear two-state system. The mu parameter controls how strongly nonlinear (and numerically stiff) it is — edit it and re-run:

model VanDerPol "Van der Pol oscillator"
  parameter Real mu = 5.0 "Nonlinearity / stiffness";
  Real x(start = 2.0);
  Real y(start = 0.0);
equation
  der(x) = y;
  der(y) = mu * (1 - x^2) * y - x;
  annotation(experiment(StopTime = 30.0));
end VanDerPol;

The same equation style handles coupled systems. Here two masses exchange energy through a weak coupling spring:

model CoupledOscillators
  parameter Real m = 1.0;
  parameter Real k = 10.0 "Outer springs";
  parameter Real kc = 1.0 "Coupling spring";
  Real x1(start = 1.0);
  Real v1(start = 0.0);
  Real x2(start = 0.0);
  Real v2(start = 0.0);
equation
  der(x1) = v1;
  m * der(v1) = -k * x1 - kc * (x1 - x2);
  der(x2) = v2;
  m * der(v2) = -k * x2 - kc * (x2 - x1);
  annotation(experiment(StopTime = 30.0));
end CoupledOscillators;

And nonlinear equations do not need to be mechanical. This predator-prey model has two interacting states and closed orbits:

model LotkaVolterra
  parameter Real alpha = 1.1 "Prey growth";
  parameter Real beta = 0.4 "Predation";
  parameter Real delta = 0.1 "Predator efficiency";
  parameter Real gamma = 0.4 "Predator death";
  Real prey(start = 10.0);
  Real predator(start = 10.0);
equation
  der(prey) = alpha * prey - beta * prey * predator;
  der(predator) = delta * prey * predator - gamma * predator;
  annotation(experiment(StopTime = 50.0));
end LotkaVolterra;

Structure: Models, Packages, Extends

Modelica organizes code in classes:

• model / block / class — components with equations.
• package — a namespace holding other classes, usually one per library.
• record — pure data structures.
• function — algorithmic code callable from equations.
• connector — interface definitions used by connect(...).

Classes compose by instantiation (declaring a component of another class) and inheritance (extends). Modifications customize an instance in place:

Here is a minimal ideal tank and its hydraulic connector. Unconnected ports have zero flow. Connected ports have equal head, and their flows sum to zero:

connector TankPort
  Real head "Hydraulic head";
  flow Real q "Volume flow into the component";
end TankPort;

model Tank
  parameter Real area = 1.0;
  Real h(start = 1.0) "Liquid head";
  TankPort inlet;
  TankPort outlet;
equation
  inlet.head = h;
  outlet.head = h;
  area * der(h) = inlet.q + outlet.q;
end Tank;

The following composition uses that Tank definition. Compile both blocks together. With no external inflow, their initial common head remains constant.

model TwoTanks
  Tank tank1(area = 2.0);
  Tank tank2(area = 0.5, h(start = 1.0));
equation
  connect(tank1.outlet, tank2.inlet);
end TwoTanks;

connect generates equality equations for potential variables and sum-to-zero equations for flow variables, which is how component-based physical modeling works without manually wiring every force and current.

What the Compiler Does With Your Equations

When you press run, Rumoca:

1. parses and resolves the model and everything it references,
2. type-checks and instantiates components with their modifications,
3. flattens the class hierarchy and connect sets into one equation system,
4. sorts and analyzes the system structurally (matching equations to unknowns, finding simultaneous blocks, tearing algebraic loops),
5. hands the prepared system to a numerical solver.

You can watch each stage with rumoca compile --emit <stage> and --inspect structure — see Inspecting and Debugging Models. The full pipeline is documented for contributors in the Rumoca Dev Guide book.

Events and Discrete Behavior

Physical models often mix continuous dynamics with discrete switching: impacts, controllers that sample, valves that open. Modelica expresses these with events, and the solver locates them precisely instead of stepping blindly past them.

When Equations and reinit

A when equation activates at the instant its condition becomes true. The classic example is the bouncing ball, which reverses its velocity at impact:

model Ball
  Real x(start = 10) "Height [m]";
  Real v(start = 1) "Velocity [m/s]";
  parameter Real g = 9.81;
  parameter Real e = 0.8 "Coefficient of restitution";
equation
  der(x) = v;
  der(v) = -g;
  when x < 0 then
    reinit(v, -e * pre(v));
  end when;
  annotation(experiment(StopTime = 10.0));
end Ball;

• when x < 0 then ... end when — fires once at the crossing instant, not continuously while the condition holds.
• pre(v) — the value of v immediately before the event.
• reinit(v, ...) — restarts the state v from a new value.

The solver detects the zero crossing of x, locates the impact time, applies the reinit, and continues. Run the example and zoom in on the velocity trace: each bounce is a clean jump, not a smoothed-over spike.

Sampled (Clocked-Style) Control

sample(start, period) fires periodically, which is the standard way to model a digital controller around a continuous plant:

model SampledControl "Continuous plant with a sampled P controller"
  parameter Real kp = 2.0 "Proportional gain";
  parameter Real Ts = 0.2 "Sample period [s]";
  parameter Real r = 1.0 "Reference";
  Real x(start = 0.0) "Plant state";
  discrete Real u(start = 0.0) "Held control signal";
equation
  der(x) = -x + u;
  when sample(0, Ts) then
    u = kp * (r - pre(x));
  end when;
  annotation(experiment(StopTime = 6.0));
end SampledControl;

u is discrete: it changes only when the when fires and is held constant (zero-order hold) in between. Try shortening Ts toward continuous control, or raising kp until the loop rings.

Hysteresis with elsewhen

elsewhen chains mutually exclusive switching conditions. A thermostat with a hysteresis band:

model Thermostat "Bang-bang temperature control with hysteresis"
  parameter Real T_set = 21.0 "Setpoint [degC]";
  parameter Real band = 1.0 "Hysteresis half-width [degC]";
  parameter Real T_amb = 5.0 "Outside temperature [degC]";
  parameter Real tau = 600.0 "Thermal time constant [s]";
  parameter Real heat = 0.05 "Heater authority [degC/s]";
  Real T(start = 15.0) "Room temperature [degC]";
  Boolean on(start = true) "Heater state";
equation
  der(T) = (T_amb - T) / tau + (if on then heat else 0.0);
  when T > T_set + band then
    on = false;
  elsewhen T < T_set - band then
    on = true;
  end when;
  annotation(experiment(StopTime = 7200.0, Solver = "rk-like"));
end Thermostat;

This model pins Solver = "rk-like" in its experiment annotation: the explicit solver handles its rapid relay switching robustly, while the default implicit solver can currently stall on it (see Solvers and Accuracy).

Conditional Expressions vs Events

An if-expression inside an equation (if on then heat else 0.0) also generates events at its switch points so the integrator never smears across a discontinuity. When a discontinuity is harmless and you want to suppress event handling, Modelica provides noEvent(...); smooth(...) asserts differentiability.

Things to Keep In Mind

• when bodies relate discrete values; use reinit to restart continuous states.
• pre(x) is only meaningful for discrete-valued variables or at event instants.
• Event-heavy models simulate best with explicit solvers today; pin Solver = "rk-like" in the experiment annotation or pass --solver rk-like.

Arrays and Discretized PDEs

Modelica has no built-in partial differential equations, but array variables plus for-equations make method of lines discretizations natural: slice the spatial domain into cells, give each cell a state, and let the compiler unroll the equations.

Cooking a Turkey

A turkey in the oven is (approximately!) a sphere heated from the outside — the 1-D radial heat equation. We split the sphere into N concentric shells, write an energy balance for each shell, and drive the outer surface with oven convection and radiation.

Press ▶ Simulate, then press play on the cross-section: colors show the temperature field conducting inward over four hours of cooking. Drag the slider to scrub through time.

model Turkey "Roasting a turkey: 1-D spherical heat equation, method of lines"
  parameter Integer N = 10 "Number of radial shells";
  final parameter Integer Np1 = N + 1;
  parameter Real M = 5.0 "Turkey mass [kg]";
  parameter Real rho = 1050.0 "Density [kg/m3]";
  parameter Real cp = 3500.0 "Specific heat [J/(kg.K)]";
  parameter Real kc = 0.5 "Thermal conductivity [W/(m.K)]";
  parameter Real T_oven = 450.0 "Oven temperature [K] (~177 degC)";
  parameter Real h = 15.0 "Convective film coefficient [W/(m2.K)]";
  parameter Real epsilon = 0.85 "Surface emissivity";
  parameter Real sigma = 5.67e-8 "Stefan-Boltzmann constant";
  parameter Real pi = 3.14159265359;
  parameter Real R = (3.0 * M / (4.0 * pi * rho)) ^ (1.0 / 3.0) "Radius [m]";
  parameter Real dr = R / N "Shell thickness [m]";
  Real T[N](each start = 277.0) "Shell temperatures [K] (fridge-cold start)";
  Real Q_cond[Np1] "Conductive heat flow across shell interfaces [W]";
  Real Q_surf "Heat into the surface from the oven [W]";
  Real r[Np1] "Interface radii [m]";
  Real A_interface[Np1] "Interface areas [m2]";
  Real V_shell[N] "Shell volumes [m3]";
  Real m_shell[N] "Shell masses [kg]";
equation
  for i in 1:Np1 loop
    r[i] = (i - 1) * dr;
    A_interface[i] = 4.0 * pi * r[i] ^ 2;
  end for;
  for i in 1:N loop
    V_shell[i] = (4.0 / 3.0) * pi * (r[i + 1] ^ 3 - r[i] ^ 3);
    m_shell[i] = rho * V_shell[i];
  end for;
  Q_cond[1] = 0.0 "No flux through the center";
  for i in 2:N loop
    Q_cond[i] = kc * A_interface[i] * (T[i - 1] - T[i]) / dr;
  end for;
  Q_cond[Np1] = 0.0;
  Q_surf = h * A_interface[Np1] * (T_oven - T[N])
    + epsilon * sigma * A_interface[Np1] * (T_oven ^ 4 - T[N] ^ 4);
  for i in 1:N - 1 loop
    m_shell[i] * cp * der(T[i]) = Q_cond[i] - Q_cond[i + 1];
  end for;
  m_shell[N] * cp * der(T[N]) = Q_cond[N] + Q_surf;
  annotation(experiment(StopTime = 14400.0, Interval = 60.0));
end Turkey;

The cross-section animation above is itself an editable JavaScript block — expand it, change the colors or geometry, and re-run ▶ Simulate to see your version. It receives the simulation results and a small helper API (api.arrayField, api.makeCanvas, api.addAnimation, api.addColorbar, api.heatColor, …) from the book’s live harness.

// Draw the turkey cross-section: concentric shells colored by temperature.
const field = api.arrayField();            // T[1..N], sorted by index
const { vMin, vMax } = api.valueRange(field.members);
const T_done = 347;                        // 74 degC: poultry-safe core temp

const size = 300;
const { ctx2d } = api.makeCanvas(size, size);
const n = field.members.length;

api.addAnimation(times, (frame) => {
  ctx2d.clearRect(0, 0, size, size);
  const maxR = size / 2 - 8;
  // Outermost shell first so inner shells paint on top.
  for (let i = n - 1; i >= 0; i--) {
    const T = field.members[i].values[frame];
    ctx2d.beginPath();
    ctx2d.arc(size / 2, size / 2, maxR * ((i + 1) / n), 0, 2 * Math.PI);
    ctx2d.fillStyle = api.heatColor((T - vMin) / (vMax - vMin));
    ctx2d.fill();
  }
  ctx2d.beginPath();
  ctx2d.arc(size / 2, size / 2, maxR, 0, 2 * Math.PI);
  ctx2d.strokeStyle = '#777';
  ctx2d.lineWidth = 2;
  ctx2d.stroke();

  const T_core = field.members[0].values[frame];
  const doneness = T_core >= T_done ? ' — done!' : '';
  return `t = ${api.formatClock(times[frame])} · core `
    + `${(T_core - 273.15).toFixed(1)} degC${doneness}`;
}, 12000);

api.addColorbar(vMin, vMax, api.heatColor);

The poultry-safe core temperature is 347 K (74 °C / 165 °F). Watch T[1] (the center) in the plot: with these parameters a 5 kg bird is not done after four hours at 177 °C — try a hotter oven, a smaller turkey, or a longer StopTime in the experiment annotation.

What to Notice in the Model

• One state per cell. Real T[N](each start = 277.0) declares N states at once; each applies the modification to every element.
• for-equations are unrolled at compile time. The loop range must be known structurally, which is why N is a parameter Integer. After flattening, the compiler sees 5 * N + 4 plain scalar equations — press Show DAE to look at them.
• Energy balances, not finite-difference formulas. Writing m_shell[i] * cp * der(T[i]) = Q_cond[i] - Q_cond[i+1] per shell, with an explicit interface flux array, conserves energy exactly by construction and reads like the physics. (This formulation follows the Dyad turkey demo.)
• Mixed boundary condition. The surface shell receives both convection (h·A·ΔT) and radiation (ε·σ·A·(T⁴_oven − T⁴)) — the T⁴ terms make the system nonlinear, which the solver handles without any special treatment.

2-D Fields: A Vibrating Membrane

The same technique extends to two dimensions with matrix states and nested for-equations. This is the 2-D wave equation on a square membrane with clamped edges, started from a Gaussian pluck in the center. The grid resolution is the N parameter — edit it (try 8 to 20 on CPU, or larger with GPU enabled) and re-run:

model Wave2D "2-D wave equation on a square membrane, method of lines"
  parameter Integer N = 20 "Grid cells per side";
  parameter Real L = 1.0 "Side length [m]";
  parameter Real c = 1.0 "Wave speed [m/s]";
  parameter Real d = 0.05 "Damping [1/s]";
  parameter Real dx = L / (N - 1);
  Real u[N, N] "Displacement";
  Real w[N, N] "Velocity";
initial equation
  for i in 1:N loop
    for j in 1:N loop
      u[i, j] = exp(-200.0 * (((i - 1) * dx - 0.5 * L) ^ 2
                            + ((j - 1) * dx - 0.5 * L) ^ 2));
      w[i, j] = 0.0;
    end for;
  end for;
equation
  for i in 1:N loop
    for j in 1:N loop
      der(u[i, j]) = w[i, j];
    end for;
  end for;
  // Fixed boundary: edges clamped to zero motion.
  for i in 1:N loop
    der(w[i, 1]) = 0.0;
    der(w[i, N]) = 0.0;
  end for;
  for j in 2:N - 1 loop
    der(w[1, j]) = 0.0;
    der(w[N, j]) = 0.0;
  end for;
  // Interior: five-point Laplacian with light damping.
  for i in 2:N - 1 loop
    for j in 2:N - 1 loop
      der(w[i, j]) = c ^ 2 * (u[i + 1, j] + u[i - 1, j] + u[i, j + 1]
                            + u[i, j - 1] - 4.0 * u[i, j]) / dx ^ 2
                   - d * w[i, j];
    end for;
  end for;
  annotation(experiment(StopTime = 2.0, Interval = 0.02, Solver = "rk-like"));
end Wave2D;

The surface below is an editable visualization script, like the turkey’s. api.matrixField() collects the u[i,j] states into a grid; the script draws a smoothed blue Three.js surface. Drag to rotate, shift-drag or right-drag to pan, and scroll to zoom:

// Animate the membrane displacement field u[i,j] as a rotatable surface.
api.hideDefaultPlot();
const field = api.matrixField();
const { vMin, vMax } = api.valueRange(field.members);
const span = Math.max(Math.abs(vMin), Math.abs(vMax)) || 1;
const { THREE } = await api.loadThree();

container.classList.add('rumoca-live-surface');
const host = document.createElement('div');
host.className = 'rumoca-live-surface-host';
container.appendChild(host);

const renderer = new THREE.WebGLRenderer({ antialias: true });
renderer.setPixelRatio(Math.min(window.devicePixelRatio || 1, 2));
renderer.setClearColor(0x071825, 1);
renderer.outputColorSpace = THREE.SRGBColorSpace;
host.appendChild(renderer.domElement);

const scene = new THREE.Scene();
scene.fog = new THREE.Fog(0x071825, 3.5, 8.5);

const camera = new THREE.PerspectiveCamera(42, 1, 0.05, 30);
const target = new THREE.Vector3(0, 0, 0);
const orbit = { theta: -0.75, phi: 0.9, radius: 3.1 };

scene.add(new THREE.HemisphereLight(0x96d7ff, 0x082744, 1.4));
const sun = new THREE.DirectionalLight(0xffffff, 1.9);
sun.position.set(-2.2, 3.5, 2.4);
scene.add(sun);

const positions = new Float32Array(field.rows * field.cols * 3);
const colors = new Float32Array(field.rows * field.cols * 3);
const uvs = new Float32Array(field.rows * field.cols * 2);
const indices = [];
for (let i = 0; i < field.rows; i++) {
  for (let j = 0; j < field.cols; j++) {
    const k = i * field.cols + j;
    uvs[k * 2] = j / Math.max(1, field.cols - 1);
    uvs[k * 2 + 1] = i / Math.max(1, field.rows - 1);
  }
}
for (let i = 0; i < field.rows - 1; i++) {
  for (let j = 0; j < field.cols - 1; j++) {
    const a = i * field.cols + j;
    const b = a + 1;
    const c = a + field.cols;
    const d = c + 1;
    indices.push(a, c, b, b, c, d);
  }
}

const geometry = new THREE.BufferGeometry();
geometry.setAttribute('position', new THREE.BufferAttribute(positions, 3));
geometry.setAttribute('color', new THREE.BufferAttribute(colors, 3));
geometry.setAttribute('uv', new THREE.BufferAttribute(uvs, 2));
geometry.setIndex(indices);

function makeWaterTexture() {
  const canvas = document.createElement('canvas');
  canvas.width = 512;
  canvas.height = 512;
  const ctx = canvas.getContext('2d');
  const gradient = ctx.createLinearGradient(0, 0, 512, 512);
  gradient.addColorStop(0, '#0b4f88');
  gradient.addColorStop(0.45, '#1391b9');
  gradient.addColorStop(1, '#56d6df');
  ctx.fillStyle = gradient;
  ctx.fillRect(0, 0, 512, 512);
  let seed = 7;
  const rand = () => {
    seed = (1664525 * seed + 1013904223) >>> 0;
    return seed / 4294967296;
  };
  ctx.globalAlpha = 0.28;
  for (let y = -80; y < 600; y += 26) {
    ctx.strokeStyle = y % 52 === 0 ? '#d8fbff' : '#74e7f1';
    ctx.lineWidth = y % 52 === 0 ? 2.2 : 1.2;
    ctx.beginPath();
    for (let x = -20; x <= 540; x += 20) {
      const yy = y + Math.sin((x + y) * 0.035) * 8 + Math.sin(x * 0.08) * 3;
      if (x === -20) {
        ctx.moveTo(x, yy);
      } else {
        ctx.lineTo(x, yy);
      }
    }
    ctx.stroke();
  }
  ctx.globalAlpha = 0.16;
  for (let i = 0; i < 180; i++) {
    const x = rand() * 512;
    const y = rand() * 512;
    const r = 0.8 + rand() * 2.4;
    ctx.fillStyle = '#e9ffff';
    ctx.beginPath();
    ctx.ellipse(x, y, r * 2.5, r, rand() * Math.PI, 0, Math.PI * 2);
    ctx.fill();
  }
  const texture = new THREE.CanvasTexture(canvas);
  texture.colorSpace = THREE.SRGBColorSpace;
  texture.wrapS = THREE.RepeatWrapping;
  texture.wrapT = THREE.RepeatWrapping;
  texture.repeat.set(1.4, 1.4);
  return texture;
}

const surface = new THREE.Mesh(
  geometry,
  new THREE.MeshStandardMaterial({
    color: 0xffffff,
    map: makeWaterTexture(),
    vertexColors: true,
    roughness: 0.28,
    metalness: 0.02,
    side: THREE.DoubleSide,
  })
);
scene.add(surface);

function updateCamera() {
  orbit.phi = Math.max(0.22, Math.min(1.38, orbit.phi));
  orbit.radius = Math.max(1.6, Math.min(7.0, orbit.radius));
  camera.position.set(
    target.x + orbit.radius * Math.sin(orbit.phi) * Math.sin(orbit.theta),
    target.y + orbit.radius * Math.cos(orbit.phi),
    target.z + orbit.radius * Math.sin(orbit.phi) * Math.cos(orbit.theta)
  );
  camera.lookAt(target);
}

function render() {
  updateCamera();
  renderer.render(scene, camera);
}

function resize() {
  const width = Math.max(320, Math.floor(host.clientWidth || 640));
  const height = Math.max(300, Math.min(500, Math.floor(width * 0.62)));
  renderer.setSize(width, height, false);
  camera.aspect = width / height;
  camera.updateProjectionMatrix();
  render();
}

function updateSurface(frameIndex) {
  const blue = new THREE.Color();
  let k = 0;
  for (let i = 1; i <= field.rows; i++) {
    for (let j = 1; j <= field.cols; j++) {
      const u = field.at(i, j)[frameIndex];
      const x = ((j - 1) / (field.cols - 1) - 0.5) * 2;
      const z = ((i - 1) / (field.rows - 1) - 0.5) * 2;
      const y = u * 0.85 / span;
      positions[k * 3] = x;
      positions[k * 3 + 1] = y;
      positions[k * 3 + 2] = z;
      const shade = 0.42 + 0.24 * (0.5 + 0.5 * (u / span));
      blue.setHSL(0.53, 0.9, shade);
      colors[k * 3] = blue.r;
      colors[k * 3 + 1] = blue.g;
      colors[k * 3 + 2] = blue.b;
      k++;
    }
  }
  geometry.attributes.position.needsUpdate = true;
  geometry.attributes.color.needsUpdate = true;
  geometry.computeVertexNormals();
  render();
}

let pointer = null;
renderer.domElement.addEventListener('contextmenu', (event) => event.preventDefault());
renderer.domElement.addEventListener('pointerdown', (event) => {
  renderer.domElement.setPointerCapture(event.pointerId);
  pointer = {
    id: event.pointerId,
    x: event.clientX,
    y: event.clientY,
    pan: event.shiftKey || event.button === 1 || event.button === 2,
  };
});
renderer.domElement.addEventListener('pointermove', (event) => {
  if (!pointer || pointer.id !== event.pointerId) {
    return;
  }
  const dx = event.clientX - pointer.x;
  const dy = event.clientY - pointer.y;
  pointer.x = event.clientX;
  pointer.y = event.clientY;
  if (pointer.pan) {
    const direction = new THREE.Vector3();
    camera.getWorldDirection(direction);
    const right = new THREE.Vector3().crossVectors(direction, camera.up).normalize();
    const up = new THREE.Vector3().copy(camera.up).normalize();
    const scale = orbit.radius * 0.0018;
    target.addScaledVector(right, -dx * scale);
    target.addScaledVector(up, dy * scale);
  } else {
    orbit.theta -= dx * 0.008;
    orbit.phi -= dy * 0.008;
  }
  render();
});
for (const type of ['pointerup', 'pointercancel']) {
  renderer.domElement.addEventListener(type, (event) => {
    if (pointer && pointer.id === event.pointerId) {
      pointer = null;
    }
  });
}
renderer.domElement.addEventListener('wheel', (event) => {
  event.preventDefault();
  orbit.radius *= Math.exp(event.deltaY * 0.001);
  render();
}, { passive: false });

new ResizeObserver(resize).observe(host);
resize();
function fixedTime(seconds) {
  return seconds.toFixed(2).padStart(5, '0');
}
api.addAnimation(times, (frame) => {
  updateSurface(frame);
  return `t = ${fixedTime(times[frame])} s`;
}, 8000);

Notice the cost of resolution: each cell adds two states (u and w), so N = 20 is an 800-state system. Rumoca compiles and simulates it fine, but output volume grows as N² — keep Interval coarse enough for the browser. The regular interior loops are also the shape Rumoca preserves as source-proven affine stencils in Solve IR: GPU/codegen targets can emit a single parametric kernel for the repeated grid operation, while CPU and embedded targets still have an exact scalar fallback.

2-D Navier–Stokes: Flow over a NACA 2412 Airfoil

The same machinery scales up to fluid dynamics. This example solves the 2-D incompressible Navier–Stokes equations around a NACA 2412 airfoil at an adjustable angle of attack, using two classic tricks that keep the system a pure ODE — exactly what the method of lines wants:

• Artificial compressibility: instead of the pressure-Poisson algebraic constraint, pressure gets its own fast dynamics der(q) = -cs² · div(V). The continuity error propagates away as an artificial acoustic wave, and no algebraic loop is needed.
• Brinkman penalization: the airfoil is not meshed. Each cell gets a smooth solid fraction sig ∈ [0, 1] (computed per grid column from the standard NACA camber and thickness polynomials) and a drag term -sig·V/tau that drives the velocity to zero, so the flow sees a solid body on a plain Cartesian grid. The mask is a differentiable tanh indicator rather than a hard if/else threshold: it is ≈1 deep inside the body, 0.5 on the geometric contour, and ≈0 in the fluid, with a transition band about one grid cell wide. That smoothness is what makes the shape sensitivities ∂sig/∂(camber, thickness) nonzero in the boundary band — the prerequisite for differentiating the flow with respect to the airfoil shape.

The freestream is horizontal and the airfoil itself pitches: each cell’s coordinates sc (chordwise) and nc (chord-normal) live in the airfoil frame, so the solid mask turns with the angle of attack the way a real wind-tunnel model would. aoa is the pre-simulation angle parameter; initial equations rotate every cell center into the frame at that angle. The model also exposes input Real aoa_cmd = aoa; an aoa_motor state follows that command with der(aoa_motor) = (aoa_cmd - aoa_motor) / aoa_tau, and the frame coordinates follow the exact rotation kinematics der(sc) = -nc·ω, der(nc) = sc·ω at the pitch rate ω = der(aoa_motor). The camber inputs mc, pc, and tk drive their own lagged motor states the same way. A structural interactive flag selects whether the motors follow the command inputs or hold their pre-simulation values. The inputs carry their pre-simulation values as declaration defaults, so a plain batch run needs no input provider. With Interactive off, the AoA slider is a normal parameter tuner and re-runs the simulation from the selected pre-simulation angle. With Interactive on, the same slider feeds aoa_cmd during stepping, so the physical airfoil angle moves through the first-order lag.

This example defaults the GPU checkbox on: the compiler’s experimental wgsl-ode backend lowers the system to WebGPU compute kernels and an in-page RK4 integrator runs them. Interior finite-volume loops are preserved as source-proven affine stencils, so the WebGPU path emits native row-parallel stencil kernels instead of rediscovering grid structure from scalarized equations. If WebGPU is unavailable the run fails with a clear error instead of silently falling back (uncheck GPU for the CPU path). wgsl-ode accepts explicit ODEs only, in f32: no algebraic variables and no events. That is why the mask is written out inside the two momentum equations instead of being a separate sig variable, and why the camber’s branch switch is wrapped in noEvent (the mean line is continuously differentiable there, so no event is needed). The named-input interactive stepping path reads aoa_cmd before each step and uses interactive = true. The run is an impulsive wind-tunnel start: the field begins at rest and the freestream sweeps in from the inlet and far-field boundaries. This is the heaviest example in the book (about 2,700 integrated states on the default grid): expect the first run to take a while.

model AirfoilFlow "2-D flow over a NACA 2412: artificial compressibility + penalization"
  parameter Integer NX = 30 "Cells along the channel";
  parameter Integer NY = 18 "Cells across the channel";
  parameter Real Lx = 4.0 "Domain length [chords]";
  parameter Real Ly = 1.5 "Domain height [chords]";
  parameter Real xle = 1.0 "Leading edge distance from inlet [chords]";
  parameter Real aoa = 8.0 "Initial/pre-simulation angle of attack [deg]";
  parameter Boolean interactive = false
    "Drive the airfoil motors from the live command inputs" annotation(Evaluate = true);
  input Real aoa_cmd = aoa "Commanded angle of attack [deg]";
  parameter Real aoa_tau = 1.0 "First-order AoA motor time constant [s]";
  parameter Real U = 1.0 "Freestream speed (horizontal)";
  parameter Real nu = 0.01 "Kinematic viscosity (Re = U/nu = 100)";
  parameter Real cs = 3.0 "Artificial-compressibility wave speed";
  parameter Real qnu = 0.01 "Pressure-mode damping diffusivity";
  parameter Real tau = 0.02 "Solid penalization time constant [s]";
  parameter Real taub = 0.05 "Boundary relaxation time constant [s]";
  parameter Real mc0 = 0.02 "Initial/pre-simulation NACA max camber";
  parameter Real pc0 = 0.4 "Initial/pre-simulation NACA camber position";
  parameter Real tk0 = 0.12 "Initial/pre-simulation NACA thickness";
  input Real mc = mc0 "Commanded NACA max camber";
  input Real pc = pc0 "Commanded NACA camber position";
  input Real tk = tk0 "Commanded NACA thickness";
  parameter Real shape_tau = 1.0 "First-order airfoil shape actuator time constant [s]";
  parameter Real dx = Lx / NX;
  parameter Real dy = Ly / NY;
  parameter Real pi = 3.14159265359;
  parameter Real epsn = 0.6 * dy "Mask transition width, chord-normal [chords]";
  parameter Real epss = 0.8 * dx "Mask transition width, chordwise [chords]";
  parameter Real tmin = 0.6 * dy "Smooth half-thickness floor: keeps the coarse mask closed";
  Real aoa_motor(start = aoa, fixed = true) "Lagged physical angle of attack [deg]";
  Real mc_motor(start = mc0, fixed = true) "Lagged NACA max camber";
  Real pc_motor(start = pc0, fixed = true) "Lagged NACA camber position";
  Real tk_motor(start = tk0, fixed = true) "Lagged NACA thickness";
  Real u[NX, NY] "x-velocity";
  Real v[NX, NY] "y-velocity";
  Real q[NX, NY] "pressure / rho";
  Real sc[NX, NY] "Chordwise coordinate in the pitched airfoil frame";
  Real nc[NX, NY] "Chord-normal coordinate in the pitched airfoil frame";
  // Flow states start at rest (default start = 0): an impulsive wind-tunnel
  // start where the freestream sweeps in through the boundary relaxation.
initial equation
  // Each cell center rotated into the airfoil frame (pitched about the
  // leading edge) at the initial angle of attack.
  for i in 1:NX loop
    for j in 1:NY loop
      sc[i, j] = ((i - 0.5) * dx - xle) * cos(aoa * pi / 180.0)
        - ((j - 0.5) * dy - Ly / 2.0) * sin(aoa * pi / 180.0);
      nc[i, j] = ((i - 0.5) * dx - xle) * sin(aoa * pi / 180.0)
        + ((j - 0.5) * dy - Ly / 2.0) * cos(aoa * pi / 180.0);
    end for;
  end for;
equation
  der(aoa_motor) =
    if interactive then (aoa_cmd - aoa_motor) / aoa_tau else 0.0;
  der(mc_motor) = if interactive then (mc - mc_motor) / shape_tau else 0.0;
  der(pc_motor) = if interactive then (pc - pc_motor) / shape_tau else 0.0;
  der(tk_motor) = if interactive then (tk - tk_motor) / shape_tau else 0.0;
  // The airfoil frame turns at the pitch rate der(aoa_motor), so each cell's
  // frame coordinates follow the exact rotation kinematics.
  for i in 1:NX loop
    for j in 1:NY loop
      der(sc[i, j]) = -nc[i, j] * pi / 180.0
        * (if interactive then (aoa_cmd - aoa_motor) / aoa_tau else 0.0);
      der(nc[i, j]) = sc[i, j] * pi / 180.0
        * (if interactive then (aoa_cmd - aoa_motor) / aoa_tau else 0.0);
    end for;
  end for;
  // Interior: momentum + artificial-compressibility continuity. The last term
  // of each momentum equation is the Brinkman drag sig * V / tau, with the
  // smooth solid fraction sig written out in place:
  //   sig = 0.5 * (1 - tanh((|nc - camber(sc)| - thickness(sc)) / epsn))
  //         * 0.5 * (1 + tanh(sc / epss)) * 0.5 * (1 + tanh((1 - sc) / epss))
  // The NACA mean line is continuously differentiable where its two branches
  // meet (sc = pc), so the branch switch needs no event: noEvent keeps it a
  // smooth expression.
  for i in 2:NX - 1 loop
    for j in 2:NY - 1 loop
      der(u[i, j]) = -u[i, j] * (u[i + 1, j] - u[i - 1, j]) / (2.0 * dx)
        - v[i, j] * (u[i, j + 1] - u[i, j - 1]) / (2.0 * dy)
        - (q[i + 1, j] - q[i - 1, j]) / (2.0 * dx)
        + nu * ((u[i + 1, j] - 2.0 * u[i, j] + u[i - 1, j]) / dx ^ 2
              + (u[i, j + 1] - 2.0 * u[i, j] + u[i, j - 1]) / dy ^ 2)
        - 0.5 * (1.0 - tanh((abs(nc[i, j]
              - noEvent(if sc[i, j] < pc_motor
                  then mc_motor / pc_motor ^ 2 * (2.0 * pc_motor * sc[i, j] - sc[i, j] ^ 2)
                  else mc_motor / (1.0 - pc_motor) ^ 2
                    * ((1.0 - 2.0 * pc_motor) + 2.0 * pc_motor * sc[i, j] - sc[i, j] ^ 2)))
            - sqrt((5.0 * tk_motor * (0.2969 * sqrt(max(sc[i, j], 0.0)) - 0.1260 * sc[i, j]
                    - 0.3516 * sc[i, j] ^ 2 + 0.2843 * sc[i, j] ^ 3
                    - 0.1036 * sc[i, j] ^ 4)) ^ 2 + tmin ^ 2)) / epsn))
          * (0.5 * (1.0 + tanh(sc[i, j] / epss)))
          * (0.5 * (1.0 + tanh((1.0 - sc[i, j]) / epss)))
          * u[i, j] / tau;
      der(v[i, j]) = -u[i, j] * (v[i + 1, j] - v[i - 1, j]) / (2.0 * dx)
        - v[i, j] * (v[i, j + 1] - v[i, j - 1]) / (2.0 * dy)
        - (q[i, j + 1] - q[i, j - 1]) / (2.0 * dy)
        + nu * ((v[i + 1, j] - 2.0 * v[i, j] + v[i - 1, j]) / dx ^ 2
              + (v[i, j + 1] - 2.0 * v[i, j] + v[i, j - 1]) / dy ^ 2)
        - 0.5 * (1.0 - tanh((abs(nc[i, j]
              - noEvent(if sc[i, j] < pc_motor
                  then mc_motor / pc_motor ^ 2 * (2.0 * pc_motor * sc[i, j] - sc[i, j] ^ 2)
                  else mc_motor / (1.0 - pc_motor) ^ 2
                    * ((1.0 - 2.0 * pc_motor) + 2.0 * pc_motor * sc[i, j] - sc[i, j] ^ 2)))
            - sqrt((5.0 * tk_motor * (0.2969 * sqrt(max(sc[i, j], 0.0)) - 0.1260 * sc[i, j]
                    - 0.3516 * sc[i, j] ^ 2 + 0.2843 * sc[i, j] ^ 3
                    - 0.1036 * sc[i, j] ^ 4)) ^ 2 + tmin ^ 2)) / epsn))
          * (0.5 * (1.0 + tanh(sc[i, j] / epss)))
          * (0.5 * (1.0 + tanh((1.0 - sc[i, j]) / epss)))
          * v[i, j] / tau;
      der(q[i, j]) = -cs ^ 2 * ((u[i + 1, j] - u[i - 1, j]) / (2.0 * dx)
                              + (v[i, j + 1] - v[i, j - 1]) / (2.0 * dy))
        + qnu * ((q[i + 1, j] - 2.0 * q[i, j] + q[i - 1, j]) / dx ^ 2
               + (q[i, j + 1] - 2.0 * q[i, j] + q[i, j - 1]) / dy ^ 2);
    end for;
  end for;
  // Inlet (left): horizontal freestream; pressure zero-gradient.
  for j in 1:NY loop
    der(u[1, j]) = (U - u[1, j]) / taub;
    der(v[1, j]) = (0.0 - v[1, j]) / taub;
    der(q[1, j]) = (q[2, j] - q[1, j]) / taub;
    // Outlet (right): zero-gradient velocities, reference pressure.
    der(u[NX, j]) = (u[NX - 1, j] - u[NX, j]) / taub;
    der(v[NX, j]) = (v[NX - 1, j] - v[NX, j]) / taub;
    der(q[NX, j]) = (0.0 - q[NX, j]) / taub;
  end for;
  // Far field (top/bottom): freestream; pressure zero-gradient.
  for i in 2:NX - 1 loop
    der(u[i, 1]) = (U - u[i, 1]) / taub;
    der(v[i, 1]) = (0.0 - v[i, 1]) / taub;
    der(q[i, 1]) = (q[i, 2] - q[i, 1]) / taub;
    der(u[i, NY]) = (U - u[i, NY]) / taub;
    der(v[i, NY]) = (0.0 - v[i, NY]) / taub;
    der(q[i, NY]) = (q[i, NY - 1] - q[i, NY]) / taub;
  end for;
  // Interval controls the output/readback cadence; __rumoca(Solver(FixedStep))
  // controls the internal fixed-step RK4 step. CFL note: the explicit step must
  // stay below the acoustic/diffusive limit ~ 1 / (cs/h + 2*nu/h^2) with
  // h = min(dx, dy); if you refine NX/NY, drop FixedStep roughly in proportion
  // or the run will diverge.
  annotation(__rumoca(Solver(FixedStep = 0.005)), experiment(StopTime = 30, Interval = 0.1, Solver = "rk-like"));
end AirfoilFlow;

// Field heatmap over the smooth solid mask (gray) with the true NACA 2412
// contour on top. The default field is vorticity, which makes separated shear
// layers and stall visible. The picker can also show pressure q or speed |V|.
// Velocity-direction and streamline overlays are computed from the returned
// u/v states only; they add no solver work or extra readback. Grid size is
// discovered from the result names, so editing NX/NY just works. Only the
// integrated states u/v/q plus the lagged airfoil motor states are read; the
// mask is recomputed below so the overlay follows the same geometry on every
// solver path.
const cell = new Map();
let NX = 0, NY = 0;
names.forEach((n, k) => {
  const m = /^([uvq])\[(\d+),(\d+)\]$/.exec(n);
  if (!m) return;
  cell.set(`${m[1]}:${m[2]},${m[3]}`, data[k]);
  if (m[1] === 'u') {
    NX = Math.max(NX, Number(m[2]));
    NY = Math.max(NY, Number(m[3]));
  }
});
const speed = (i, j, f) => {
  const u = cell.get(`u:${i},${j}`);
  const v = cell.get(`v:${i},${j}`);
  if (!u || !v) return 0;
  return Math.hypot(u[f], v[f]);
};
const press = (i, j, f) => {
  const q = cell.get(`q:${i},${j}`);
  return q ? q[f] : 0;
};
const fieldDx = api.parameter('Lx', 4.0) / NX;
const fieldDy = api.parameter('Ly', 1.5) / NY;
const clippedIndex = (value, hi) => Math.max(1, Math.min(hi, value));
const vorticity = (i, j, f) => {
  const im = clippedIndex(i - 1, NX);
  const ip = clippedIndex(i + 1, NX);
  const jm = clippedIndex(j - 1, NY);
  const jp = clippedIndex(j + 1, NY);
  const dvdx = (cell.get(`v:${ip},${j}`)?.[f] - cell.get(`v:${im},${j}`)?.[f])
    / ((ip - im) * fieldDx);
  const dudy = (cell.get(`u:${i},${jp}`)?.[f] - cell.get(`u:${i},${jm}`)?.[f])
    / ((jp - jm) * fieldDy);
  return Number.isFinite(dvdx) && Number.isFinite(dudy) ? dvdx - dudy : 0;
};
// Velocity arrows use a global speed reference so their visibility does not
// flicker as the colorbar rescales frame-by-frame.
let vMax = 0;
for (let f = 0; f < times.length; f += 5) {
  for (let i = 1; i <= NX; i++) {
    for (let j = 1; j <= NY; j++) {
      vMax = Math.max(vMax, speed(i, j, f));
    }
  }
}
if (vMax <= 0) vMax = 1;
const clamp01 = (x) => Math.max(0, Math.min(1, x));

function framePressureRange(frame) {
  let lo = Infinity, hi = -Infinity;
  for (let i = 1; i <= NX; i++) {
    for (let j = 1; j <= NY; j++) {
      const q = press(i, j, frame);
      if (!Number.isFinite(q)) continue;
      lo = Math.min(lo, q);
      hi = Math.max(hi, q);
    }
  }
  if (!(hi > lo)) return { lo: -1, hi: 1 };
  return { lo, hi };
}

function frameSpeedRange(frame) {
  let hi = 0;
  for (let i = 1; i <= NX; i++) {
    for (let j = 1; j <= NY; j++) hi = Math.max(hi, speed(i, j, frame));
  }
  return { lo: 0, hi: hi > 0 ? hi : 1 };
}

function frameVorticityRange(frame) {
  let span = 0;
  for (let i = 1; i <= NX; i++) {
    for (let j = 1; j <= NY; j++) {
      const w = vorticity(i, j, frame);
      if (Number.isFinite(w)) span = Math.max(span, Math.abs(w));
    }
  }
  span = span > 0 ? span : 1;
  return { lo: -span, hi: span };
}

// Available fields. `norm` maps a cell value into [0,1] for the heat colormap.
const fields = {
  vorticity: {
    label: 'Vorticity',
    range: frameVorticityRange,
    value: vorticity,
  },
  q: {
    label: 'Pressure q',
    range: framePressureRange,
    value: press,
  },
  speed: {
    label: 'Speed |V|',
    range: frameSpeedRange,
    value: speed,
  },
};
let mode = 'vorticity';   // default to the field that shows stall/separation.
let refreshColorbar = () => {};
// Geometry for the overlay — keep in sync with the model parameters.
// The command inputs seed the run; the lagged motor states drive the moving
// contour and mask frame-by-frame.
const mc0 = api.parameter('mc', api.parameter('mc0', 0.02));
const pc0 = api.parameter('pc', api.parameter('pc0', 0.4));
const tk0 = api.parameter('tk', api.parameter('tk0', 0.12));
const geo = {
  Lx: api.parameter('Lx', 4.0),
  Ly: api.parameter('Ly', 1.5),
  xle: api.parameter('xle', 1.0),
  mc: mc0,
  pc: pc0,
  tk: tk0,
};
const aoa = api.parameter('aoa', 8.0);
const aoaSeries = api.series('aoa_motor') || [aoa];
const mcSeries = api.series('mc_motor') || [mc0];
const pcSeries = api.series('pc_motor') || [pc0];
const tkSeries = api.series('tk_motor') || [tk0];
const frameSeriesValue = (series, frame, fallback) => {
  const value = series[Math.min(frame, series.length - 1)];
  return Number.isFinite(value) ? value : fallback;
};
const frameAoa = (frame) => frameSeriesValue(aoaSeries, frame, aoa);
function refreshGeometryParameters(frame) {
  geo.xle = api.parameter('xle', 1.0);
  geo.mc = frameSeriesValue(mcSeries, frame, mc0);
  geo.pc = Math.max(1e-3, Math.min(0.999, frameSeriesValue(pcSeries, frame, pc0)));
  geo.tk = Math.max(1e-6, frameSeriesValue(tkSeries, frame, tk0));
};
const camber = (sc) => sc < geo.pc
  ? geo.mc / geo.pc ** 2 * (2 * geo.pc * sc - sc ** 2)
  : geo.mc / (1 - geo.pc) ** 2 * ((1 - 2 * geo.pc) + 2 * geo.pc * sc - sc ** 2);
const halfThick = (sc) => 5 * geo.tk * (0.2969 * Math.sqrt(sc) - 0.1260 * sc
  - 0.3516 * sc ** 2 + 0.2843 * sc ** 3 - 0.1036 * sc ** 4);

// Smooth solid fraction sig(i, j) in [0, 1], recomputed here exactly as the
// model does so the overlay never depends on the solver returning algebraics.
// The grid spacing and tanh band widths mirror the AirfoilFlow parameters.
const dx = geo.Lx / NX, dy = geo.Ly / NY;
const epsn = api.parameter('epsn', 0.6 * dy);
const epss = api.parameter('epss', 0.8 * dx);
const tmin = api.parameter('tmin', 0.6 * dy);
function maskAt(i, j, angleDeg) {
  const ca = Math.cos(angleDeg * Math.PI / 180);
  const sa = Math.sin(angleDeg * Math.PI / 180);
  const xa = (i - 0.5) * dx - geo.xle;        // chord-frame offsets, pitched
  const ya = (j - 0.5) * dy - geo.Ly / 2;     // about the leading edge
  const sc = xa * ca - ya * sa;               // chordwise coordinate
  const nc = xa * sa + ya * ca;               // chord-normal coordinate
  const s = Math.max(sc, 0);
  const traw = 5 * geo.tk * (0.2969 * Math.sqrt(s) - 0.1260 * sc
    - 0.3516 * sc ** 2 + 0.2843 * sc ** 3 - 0.1036 * sc ** 4);
  const teff = Math.sqrt(traw ** 2 + tmin ** 2);   // softly floored half-thick
  const dthick = Math.abs(nc - camber(sc)) - teff; // signed distance to surface
  return 0.5 * (1 - Math.tanh(dthick / epsn))      // inside the thickness band
    * 0.5 * (1 + Math.tanh(sc / epss))             // past the leading edge
    * 0.5 * (1 + Math.tanh((1 - sc) / epss));      // before the trailing edge
}
const solidCache = new Map();
function solidFracFor(angleDeg) {
  const key = [
    angleDeg.toFixed(3),
    geo.xle.toFixed(4),
    geo.mc.toFixed(4),
    geo.pc.toFixed(4),
    geo.tk.toFixed(4),
  ].join(':');
  const cached = solidCache.get(key);
  if (cached) return cached;

  const solidFrac = new Map();
  for (let i = 1; i <= NX; i++) {
    for (let j = 1; j <= NY; j++) solidFrac.set(`${i},${j}`, maskAt(i, j, angleDeg));
  }
  solidCache.set(key, solidFrac);
  return solidFrac;
}

const W = 600;
const H = Math.round(W * (geo.Ly / geo.Lx));
const { ctx2d } = api.makeCanvas(W, H);
const cw = W / NX;
const ch = H / NY;
const px = (x) => (x / geo.Lx) * W;                    // physical x -> canvas
const py = (y) => H - ((y + geo.Ly / 2) / geo.Ly) * H; // physical y -> canvas

function airfoilFrame(angleDeg) {
  const ca = Math.cos(angleDeg * Math.PI / 180);
  const sa = Math.sin(angleDeg * Math.PI / 180);
  return {
    fx: (sc, h) => geo.xle + sc * ca + h * sa,
    fy: (sc, h) => -sc * sa + h * ca,
  };
}

function drawAirfoil(angleDeg) {
  const { fx, fy } = airfoilFrame(angleDeg);
  ctx2d.beginPath();
  for (let k = 0; k <= 60; k++) {            // upper surface, LE -> TE
    const sc = k / 60;
    const h = camber(sc) + halfThick(sc);
    const fn = k === 0 ? 'moveTo' : 'lineTo';
    ctx2d[fn](px(fx(sc, h)), py(fy(sc, h)));
  }
  for (let k = 60; k >= 0; k--) {            // lower surface, TE -> LE
    const sc = k / 60;
    const h = camber(sc) - halfThick(sc);
    ctx2d.lineTo(px(fx(sc, h)), py(fy(sc, h)));
  }
  ctx2d.closePath();
  ctx2d.fillStyle = '#111';
  ctx2d.fill();
  ctx2d.strokeStyle = '#fff';
  ctx2d.lineWidth = 1;
  ctx2d.stroke();
}

let showVelocityDirections = true;

function drawVelocityDirections(frame, solidFrac) {
  if (!showVelocityDirections || !NX || !NY || !Number.isFinite(vMax)) return;

  const stride = Math.max(1, Math.ceil(Math.max(NX, NY) / 28));
  const len = Math.max(4, 0.65 * stride * Math.min(cw, ch));

  function sampledVelocity(i0, j0) {
    let su = 0, sv = 0, sw = 0;
    for (let di = -1; di <= 1; di++) {
      for (let dj = -1; dj <= 1; dj++) {
        const i = i0 + di, j = j0 + dj;
        if (i < 1 || i > NX || j < 1 || j > NY) continue;

        const m = solidFrac.get(`${i},${j}`) ?? 0;

        const uArr = cell.get(`u:${i},${j}`);
        const vArr = cell.get(`v:${i},${j}`);
        if (!uArr || !vArr) continue;

        const u = uArr[frame];
        const v = vArr[frame];
        if (!Number.isFinite(u) || !Number.isFinite(v)) continue;

        const w = (di === 0 && dj === 0 ? 2 : 1) * Math.max(0.05, 1 - m);
        su += w * u;
        sv += w * v;
        sw += w;
      }
    }
    return sw > 0 ? [su / sw, sv / sw] : [NaN, NaN];
  }

  ctx2d.save();
  ctx2d.lineCap = 'round';
  ctx2d.lineJoin = 'round';

  for (let i = 1; i <= NX; i += stride) {
    for (let j = 1; j <= NY; j += stride) {
      const [u, v] = sampledVelocity(i, j);
      const sp = Math.hypot(u, v);
      if (!Number.isFinite(sp) || sp <= 0) continue;

      const ux = u / sp;
      const uy = v / sp;
      const cx = (i - 0.5) * cw;
      const cy = (NY - j + 0.5) * ch;
      const dxp = ux * len;
      const dyp = -uy * len; // Physical +v points up; canvas +y points down.
      const x1 = cx - 0.5 * dxp;
      const y1 = cy - 0.5 * dyp;
      const x2 = cx + 0.5 * dxp;
      const y2 = cy + 0.5 * dyp;
      const ang = Math.atan2(y2 - y1, x2 - x1);
      const head = Math.max(3, 0.25 * len);
      const alpha = 0.35 + 0.55 * clamp01(sp / (0.45 * vMax));

      ctx2d.strokeStyle = `rgba(0,0,0,${0.55 * alpha})`;
      ctx2d.lineWidth = 3.5;
      ctx2d.beginPath();
      ctx2d.moveTo(x1, y1);
      ctx2d.lineTo(x2, y2);
      ctx2d.lineTo(
        x2 - head * Math.cos(ang - Math.PI / 6),
        y2 - head * Math.sin(ang - Math.PI / 6)
      );
      ctx2d.moveTo(x2, y2);
      ctx2d.lineTo(
        x2 - head * Math.cos(ang + Math.PI / 6),
        y2 - head * Math.sin(ang + Math.PI / 6)
      );
      ctx2d.stroke();

      ctx2d.strokeStyle = `rgba(255,255,255,${alpha})`;
      ctx2d.lineWidth = 1.3;
      ctx2d.beginPath();
      ctx2d.moveTo(x1, y1);
      ctx2d.lineTo(x2, y2);
      ctx2d.lineTo(
        x2 - head * Math.cos(ang - Math.PI / 6),
        y2 - head * Math.sin(ang - Math.PI / 6)
      );
      ctx2d.moveTo(x2, y2);
      ctx2d.lineTo(
        x2 - head * Math.cos(ang + Math.PI / 6),
        y2 - head * Math.sin(ang + Math.PI / 6)
      );
      ctx2d.stroke();
    }
  }

  ctx2d.restore();
}

let showStreamlines = true;

function sampleGridValue(getValue, x, y, frame) {
  const fi = x / dx + 0.5;
  const fj = (y + geo.Ly / 2) / dy + 0.5;
  if (fi < 1 || fi > NX || fj < 1 || fj > NY) return NaN;

  const i0 = Math.max(1, Math.min(NX - 1, Math.floor(fi)));
  const j0 = Math.max(1, Math.min(NY - 1, Math.floor(fj)));
  const tx = clamp01(fi - i0);
  const ty = clamp01(fj - j0);

  const v00 = getValue(i0, j0, frame);
  const v10 = getValue(i0 + 1, j0, frame);
  const v01 = getValue(i0, j0 + 1, frame);
  const v11 = getValue(i0 + 1, j0 + 1, frame);
  if (![v00, v10, v01, v11].every(Number.isFinite)) return NaN;

  const a = v00 * (1 - tx) + v10 * tx;
  const b = v01 * (1 - tx) + v11 * tx;
  return a * (1 - ty) + b * ty;
}

function sampleVelocityAt(x, y, frame) {
  const u = sampleGridValue((i, j, f) => cell.get(`u:${i},${j}`)?.[f], x, y, frame);
  const v = sampleGridValue((i, j, f) => cell.get(`v:${i},${j}`)?.[f], x, y, frame);
  return Number.isFinite(u) && Number.isFinite(v) ? [u, v] : null;
}

function sampleSolidAt(x, y, solidFrac) {
  return sampleGridValue((i, j) => solidFrac.get(`${i},${j}`) ?? 0, x, y, 0);
}

function traceStreamline(seed, frame, solidFrac) {
  const ds = geo.Lx / 120;
  const pts = [];
  let x = seed.x, y = seed.y;

  for (let step = 0; step < 170; step++) {
    if (x < 0 || x > geo.Lx || y < -geo.Ly / 2 || y > geo.Ly / 2) break;
    if (sampleSolidAt(x, y, solidFrac) > 0.85) break;

    const vel = sampleVelocityAt(x, y, frame);
    if (!vel) break;

    const sp = Math.hypot(vel[0], vel[1]);
    if (!Number.isFinite(sp) || sp <= 0) break;

    pts.push([px(x), py(y)]);

    let ux = vel[0] / sp;
    let uy = vel[1] / sp;
    const midVel = sampleVelocityAt(x + 0.5 * ds * ux, y + 0.5 * ds * uy, frame);
    const midSpeed = midVel ? Math.hypot(midVel[0], midVel[1]) : 0;
    if (Number.isFinite(midSpeed) && midSpeed > 0) {
      ux = midVel[0] / midSpeed;
      uy = midVel[1] / midSpeed;
    }

    x += ds * ux;
    y += ds * uy;
  }

  return pts;
}

function streamlineSeeds() {
  const seeds = [];
  const inletX = 0.5 * dx;
  for (let k = 1; k <= 13; k++) {
    seeds.push({
      x: inletX,
      y: -geo.Ly / 2 + (k / 14) * geo.Ly,
    });
  }
  return seeds;
}

function drawStreamlines(frame, solidFrac) {
  if (!showStreamlines || !NX || !NY) return;

  ctx2d.save();
  ctx2d.lineCap = 'round';
  ctx2d.lineJoin = 'round';

  for (const seed of streamlineSeeds()) {
    const pts = traceStreamline(seed, frame, solidFrac);
    if (pts.length < 2) continue;

    ctx2d.beginPath();
    pts.forEach(([x, y], index) => {
      if (index === 0) ctx2d.moveTo(x, y);
      else ctx2d.lineTo(x, y);
    });
    ctx2d.strokeStyle = 'rgba(0,0,0,0.45)';
    ctx2d.lineWidth = 3;
    ctx2d.stroke();

    ctx2d.beginPath();
    pts.forEach(([x, y], index) => {
      if (index === 0) ctx2d.moveTo(x, y);
      else ctx2d.lineTo(x, y);
    });
    ctx2d.strokeStyle = 'rgba(255,255,255,0.72)';
    ctx2d.lineWidth = 1.15;
    ctx2d.stroke();
  }

  ctx2d.restore();
}

let lastFrame = 0;
const anim = api.addAnimation(times, (frame) => {
  lastFrame = frame;
  refreshGeometryParameters(frame);
  const fld = fields[mode];
  const range = fld.range(frame);
  const span = Math.max(1e-12, range.hi - range.lo);
  const angle = frameAoa(frame);
  const solidFrac = solidFracFor(angle);
  refreshColorbar(range);
  for (let i = 1; i <= NX; i++) {
    for (let j = 1; j <= NY; j++) {
      // The mask is smooth (sig in [0,1]); shade the selected field, then
      // overlay gray with opacity = sig so the differentiable boundary band
      // shows as a soft halo rather than a hard on/off cell edge.
      const m = solidFrac.get(`${i},${j}`);
      // j = 1 is the bottom row: flip the y axis for drawing.
      const rx = (i - 1) * cw, ry = (NY - j) * ch;
      ctx2d.fillStyle = api.heatColor(clamp01((fld.value(i, j, frame) - range.lo) / span));
      ctx2d.fillRect(rx, ry, cw + 1, ch + 1);
      if (m > 0.01) {
        ctx2d.fillStyle = `rgba(70,70,70,${Math.min(1, m)})`;
        ctx2d.fillRect(rx, ry, cw + 1, ch + 1);
      }
    }
  }
  try {
    drawStreamlines(frame, solidFrac);
    drawVelocityDirections(frame, solidFrac);
  } catch (e) {
    console.warn('Velocity overlay failed:', e);
    showStreamlines = false;
    showVelocityDirections = false;
  }
  drawAirfoil(angle);
  return `t = ${times[frame].toFixed(1)} s · ${fld.label} `
    + `∈ [${api.formatTick(range.lo)}, ${api.formatTick(range.hi)}]`;
}, 10000);

// Colorbar we can relabel when the field changes (api.addColorbar is static).
const bar = document.createElement('div');
bar.className = 'rumoca-live-radial-colorbar';
const grad = document.createElement('span');
grad.className = 'rumoca-live-radial-gradient';
const stops = [];
for (let i = 0; i <= 10; i++) stops.push(api.heatColor(i / 10));
grad.style.background = `linear-gradient(to right, ${stops.join(', ')})`;
const loEl = document.createElement('span'), hiEl = document.createElement('span');
bar.append(loEl, grad, hiEl);
container.appendChild(bar);
refreshColorbar = (range = fields[mode].range(lastFrame)) => {
  loEl.textContent = api.formatTick(range.lo);
  hiEl.textContent = api.formatTick(range.hi);
};
refreshColorbar();

// Field picker: pressure (default) or speed. Switching repaints the current
// frame and relabels the colorbar; no recompile or re-run needed.
const fieldRow = document.createElement('div');
fieldRow.className = 'rumoca-live-tuner';
const fieldLabel = document.createElement('span');
fieldLabel.textContent = 'Field';
const fieldSel = document.createElement('select');
[
  ['vorticity', 'Vorticity'],
  ['q', 'Pressure'],
  ['speed', 'Speed |V|'],
].forEach(([value, text]) => {
  const opt = document.createElement('option');
  opt.value = value; opt.textContent = text;
  fieldSel.appendChild(opt);
});
fieldSel.value = mode;
fieldSel.addEventListener('change', () => {
  mode = fieldSel.value;
  refreshColorbar();
  anim.redraw(lastFrame);
});
const dirLabel = document.createElement('label');
dirLabel.style.display = 'inline-flex';
dirLabel.style.alignItems = 'center';
dirLabel.style.gap = '0.35rem';
dirLabel.style.marginLeft = '0.75rem';
const dirCheck = document.createElement('input');
dirCheck.type = 'checkbox';
dirCheck.checked = showVelocityDirections;
dirCheck.addEventListener('change', () => {
  showVelocityDirections = dirCheck.checked;
  anim.redraw(lastFrame);
});
dirLabel.append(dirCheck, document.createTextNode('Velocity direction'));

const streamLabel = document.createElement('label');
streamLabel.style.display = 'inline-flex';
streamLabel.style.alignItems = 'center';
streamLabel.style.gap = '0.35rem';
streamLabel.style.marginLeft = '0.75rem';
const streamCheck = document.createElement('input');
streamCheck.type = 'checkbox';
streamCheck.checked = showStreamlines;
streamCheck.addEventListener('change', () => {
  showStreamlines = streamCheck.checked;
  anim.redraw(lastFrame);
});
streamLabel.append(streamCheck, document.createTextNode('Streamlines'));

fieldRow.append(fieldLabel, fieldSel, dirLabel, streamLabel);
container.appendChild(fieldRow);

// Pitch the airfoil. With Interactive off this is a normal pre-run `aoa`
// parameter override. With Interactive on, the same slider drives the named
// model input `aoa_cmd`; the model's `aoa_motor` state follows it with a
// first-order lag.
api.addTuner('aoa', {
  min: -45,
  max: 45,
  step: 1,
  value: aoa,
  label: 'AoA °',
  interactiveInput: 'aoa_cmd',
});

api.addTuner('mc', {
  min: 0,
  max: 0.08,
  step: 0.005,
  value: mc0,
  label: 'Camber',
  interactiveInput: 'mc',
});

api.addTuner('pc', {
  min: 0.2,
  max: 0.8,
  step: 0.05,
  value: pc0,
  label: 'Camber pos',
  interactiveInput: 'pc',
});

api.addTuner('tk', {
  min: 0.06,
  max: 0.24,
  step: 0.01,
  value: tk0,
  label: 'Thickness',
  interactiveInput: 'tk',
});

In the animation, the black shape is the true NACA 2412 contour and the gray haze around it is the smooth solid fraction sig — darker toward the core of the body, fading through the one-cell tanh transition band to the fluid. That soft edge is the differentiable mask the flow actually feels at this resolution. The Field picker chooses what the heatmap shows: it defaults to vorticity, so separated shear layers, vortices, and stall are visible directly. Switch to Pressure to see the stagnation/suction pattern that drives lift, or Speed |V| to see the kinematics. Watch the impulsive start settle: the stagnation point appears at the leading edge, flow accelerates over the upper surface, and the wake trails downstream. Things to try:

• Slide AoA to 0 — the wake straightens and the up/down asymmetry mostly disappears (the residual comes from camber, the 2 in 2412). This is a simulation parameter update; the GPU path refreshes the prepared vectors and reruns without relowering the model.
• Slide AoA negative — the airfoil visibly pitches nose-down and the suction side flips.
• Slide AoA toward 25–45 in Interactive mode — the upper-surface shear layer separates and the vorticity/streamline overlays show the qualitative onset of stall.
• Switch Field between Vorticity, Pressure, and Speed — the same run, recolored instantly with no recompile.

Honest caveats: at this grid and Reynolds number (U/nu = 100), this is a qualitative separated-flow visualization, not an aerodynamic prediction. The half-thickness is softly floored to tmin (a fraction of a cell) so the thin profile stays closed, pressure uses artificial compressibility plus damping, and central differencing limits how low the viscosity may go. Resolving boundary layers at flight Reynolds numbers needs orders of magnitude more cells and a more specialized incompressible-flow discretization.

Scaling the Resolution

Increase the cell counts for finer fields. Each extra cell adds states and equations; the structural analysis and solver scale with the system size. For 1-D problems, tens of cells are usually plenty; for large 2-D/3-D fields you would generate the Modelica programmatically or move to a dedicated PDE solver. GPU-accelerated execution of large discretized fields is on the roadmap — the targets table already includes experimental CUDA backends (rumoca targets), and the same solve-IR pathway is how a WebGPU/WGSL backend would land.

Neural ODEs

A Neural ODE replaces a hand-written right-hand side with a neural network:

der(x) = f_theta(x, time, inputs)

The model is still an ODE. The difference is that f_theta is an MLP, CNN, or other differentiable computation with trained weights theta. In Rumoca, those weights are ordinary Modelica parameters, and the state trajectory is simulated with the same solver/runtime path used by physical models.

The Correll Lab Neural ODE review uses the common spiral demonstration: train a network to describe a two-dimensional rotating trajectory, then integrate the learned vector field. The repository example below follows that shape. It is not a training loop; it is the deploy/simulate side of the workflow, where a trained or initialized network is written as tensor-shaped Modelica arrays.

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.neural_ode_tensor.toml

Run the same scenario from the repository root:

rumoca sim -c examples/simulation/rumoca-scenario.neural_ode_tensor.toml

What the Model Does

The core equations are three matrix-vector layers with smooth activations. This excerpt keeps the activation equations inside their respective for i in 1:nHidden loop sections in the complete NeuralODETensor.mo model:

z1 = W1 * x;
a1[i] = tanh(z1[i] + b1[i]);
z2 = Wmid * a1;
a2[i] = tanh(z2[i] + bmid[i]);
der(x) = W2 * a2 + b2;

x[1] and x[2] are the latent state coordinates. The default parameters seed the first two hidden channels with a damped spiral vector field and add small dense tensor layers around it, so the phase portrait behaves like the blog-post spiral without requiring an external training step.

The useful part for larger models is the array shape:

• W1 * x, Wmid * a1, and W2 * a2 are Modelica matrix-vector products.
• Rumoca lowers these as native tensor operations in Solve IR instead of requiring you to write every scalar multiply-add by hand.
• CPU, browser, and codegen paths still have scalar fallback behavior, but the model source keeps the neural network structure visible.

Parameter Count

For this two-hidden-layer network with nState = 2 and hidden width H:

W1:   H x 2
b1:   H
Wmid: H x H
bmid: H
W2:   2 x H
b2:   2
total = H^2 + 6H + 2

The committed scenario uses H = 32, which is small enough for a quick browser run and gives 1186 trainable parameters. To exercise a roughly 100k-parameter model, edit nHidden in examples/models/NeuralODETensor.mo:

parameter Integer nHidden(min = 2) = 314;

That instantiates 100482 trainable parameters. For this size, prefer the CLI over the browser guide page, because the report and default plots can dominate runtime and memory.

Training and Backpropagation

The examples above are simulation/deployment examples. They keep weights as Modelica parameters and run the learned vector field. Rumoca’s long-term training surface lives one layer above Modelica in the rumoca-opt crate: Modelica defines the differentiable model, while rumoca-opt defines the trainables, objective, gradient method, and optimizer.

The first supported API shape is RHS fitting, which is the core Neural ODE vector-field training problem:

The following code uses compiled from compiling the unchanged NeuralODETensor.mo model and a two-element target_derivatives vector. It selects independent first- and second-layer biases from parameter_slots(). The excerpt is included directly from its executable Rust test:

#![allow(unused)]
fn main() {
    let mut model = rumoca_opt::DifferentiableModel::from_dae_default(
        &compiled.dae,
        &rumoca_solver::SimOptions::default(),
    )
    .expect("prepare the compiled neural model");
    let trainables = rumoca_opt::TrainableSet::by_names(&model, &["b1[1]", "bmid[1]"])
        .expect("select the documented independent network biases");
    let objective = rumoca_opt::RhsMseObjective::new(0.0, target_derivatives);
    let gradient = rumoca_opt::rhs_mse_value_and_gradient(
        &model,
        &objective,
        &trainables,
        rumoca_opt::GradientMode::Auto,
    )
    .expect("differentiate the neural RHS fitting objective");
}

GradientMode::Auto uses reverse-mode VJP for pure ODE models and falls back to forward-mode parameter Jacobians for models with solver algebraics. That keeps backpropagation correct now while preserving a clear extension point for future algebraic-projection reverse mode. A simple optimizer loop is available through rumoca_opt::GradientDescent.

NeuralODEBackprop.mo is also included as a native-equation demonstration: the trainable weights are states and the backprop equations are written directly in Modelica. It proves the tensor math and training dynamics can run inside Rumoca, but it is intentionally not the preferred public API.

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.neural_ode_backprop.toml

Benchmarking Against JAX

The repository includes a reproducible Neural ODE parity benchmark:

python3 -m venv target/jax-bench-venv
target/jax-bench-venv/bin/python -m pip install --upgrade pip
target/jax-bench-venv/bin/python -m pip install jax numpy
target/jax-bench-venv/bin/python examples/benchmarks/neural_ode_jax_parity.py --platform cpu

For CUDA-capable JAX installs, use --platform gpu. The script compares NeuralODEBackprop1k against an equivalent JAX implementation with x64 enabled. It reports:

• Rumoca front-end compile_seconds, Solve IR/runtime prepare_seconds, and hot simulation time.
• JAX JIT compile-plus-first-execute time for the reverse-mode RHS, manual RHS, diagnostics, and full integration.
• A comparison block with Rumoca ready time versus JAX integrated compile time, hot-loop ratio, structural-prepare fraction, and max parity error.
• Accuracy comparisons for the six Rumoca output observables at the initial and final times, plus all 1082 JAX manual gradients versus jax.value_and_grad. The visible-output CSV does not expose internal trainable states or gradients, so the report explicitly leaves their cross-backend parity unavailable.

Use Rumoca ready_seconds versus JAX integrate_compile_seconds when you care about first-run wait time. The benchmark also makes clear which Rumoca path is being measured: the Solve IR CPU simulator executes the Modelica backpropagation equations, while the JAX side uses jax.value_and_grad for the same RHS. Mandatory Linux wheel CI and cargo xtask python test --require-guide also check all eleven output samples against the independent JAX trajectory. Timings remain measurements; numerical assertions determine correctness.

Using Trained Weights

A typical workflow is:

1. Train the Neural ODE in Python/JAX/PyTorch using your preferred optimizer.
2. Export the learned arrays with the same shapes as W1, b1, Wmid, bmid, W2, and b2.
3. Replace the parameter declarations in the Modelica model, or generate a small Modelica package containing those parameter arrays.
4. Simulate, inspect, or generate code from the scenario.

Keep activation functions smooth when possible. tanh is a good default for simulation because it is continuous and differentiable. If your training code uses piecewise activations, check event behavior and solver step size before scaling up.

Predator-Prey Neural ODE

Lotka-Volterra predator-prey dynamics are a more interesting Neural ODE test case because the state has physical meaning, nonlinear interaction, and a phase portrait that should remain coherent over long horizons. The original Neural ODE paper frames these models as continuous-depth dynamics, while recent SciML examples use Lotka-Volterra systems as benchmarks for both full Neural ODEs and hybrid Universal Differential Equations.

This second repository example takes the hybrid view. The state stores normalized prey and predator populations. The neural network reads normalized prey, normalized predator, and a seasonal forcing term, then outputs per-capita growth rates:

features[1] = population[1] / preyEquilibrium - 1.0;
features[2] = population[2] / predatorEquilibrium - 1.0;
features[3] = sin(seasonRate * time);
z1 = W1 * features;
z2 = Wmid * a1;
growth = W2 * a2 + b2;
der(population[1]) = population[1] * growth[1];
der(population[2]) = population[2] * growth[2];

Run it in the guide:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.neural_predator_prey.toml

Or from the repository root:

rumoca sim -c examples/simulation/rumoca-scenario.neural_predator_prey.toml

The default width is small for report size. With hidden width H, this network has:

W1:   H x 3
b1:   H
Wmid: H x H
bmid: H
W2:   2 x H
b2:   2
total = H^2 + 7H + 2

H = 313 instantiates 100162 trainable parameters. The two large matrix products, W1 * features and Wmid * a1, remain tensor-shaped in the model source, so the same example scales from a small guide run to a much larger compiled model.

Larger Latent Neural ODE

The Latent ODE idea is to evolve a hidden continuous state with an ODE and decode that latent trajectory into observed signals. That shape needs more parameters than the two-state examples above: the vector field lives in latent space, and the decoder is another learned map.

NeuralLatentOscillator uses eight latent states, three decoded output channels, time features, a large recurrent tensor block, and a learned decoder:

z1 = W1 * features;
z2 = Wmid * a1;
latentResidual = Wdyn * a2 + bdyn;
observed = Wobs * a2 + bobs;

The baseline dynamics are four lightly damped oscillator pairs. The neural network adds a learned residual to the latent derivative and decodes the hidden trajectory into three observable signals:

der(latent[2 * k - 1]) =
  -damping * latent[2 * k - 1] - freq[k] * latent[2 * k]
  + residualGain * latentResidual[2 * k - 1];
der(latent[2 * k]) =
  freq[k] * latent[2 * k - 1] - damping * latent[2 * k]
  + residualGain * latentResidual[2 * k];

Run it in the guide:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.neural_latent_oscillator.toml

Or from the repository root:

rumoca sim -c examples/simulation/rumoca-scenario.neural_latent_oscillator.toml

With nLatent = 8, nFeature = 10, and nObserved = 3, hidden width H gives:

W1:   H x 10
b1:   H
Wmid: H x H
bmid: H
Wdyn: 8 x H
bdyn: 8
Wobs: 3 x H
bobs: 3
total = H^2 + 23H + 11

The committed default H = 128 instantiates 19339 trainable parameters, so this is already much larger than the spiral and predator-prey guide examples while staying reasonable for the examples smoke gate. For heavier runs, set nHidden = 256 for 71435 trainable parameters, or nHidden = 512 for 273931 trainable parameters while keeping the same source structure.

Using Modelica Libraries

Real models build on libraries — most importantly the Modelica Standard Library (MSL). This chapter explains how Rumoca finds library code.

Source Roots

A source root is a file or directory added to the compiler’s search path. When a model references Modelica.Blocks.Continuous.PID, Rumoca resolves the name through the source roots you configured.

On the command line, pass --source-root (repeatable):

rumoca sim my_model.mo \
  --model MyPackage.MyModel \
  --source-root target/msl/ModelicaStandardLibrary-4.1.0 \
  --source-root helper.mo

In a rumoca-scenario.toml scenario, use the top-level source_roots key with paths relative to the scenario file:

source_roots = ["../modelica_libraries"]

For workspace-wide editor, playground, and docs examples, use rumoca-workspace.toml (see VS Code Extension).

MODELICAPATH

Rumoca also honors the standard MODELICAPATH environment variable. Entries (:-separated) are appended after explicit --source-root flags, so system-wide library installs resolve without per-command flags.

Pinned Dependencies for the Repository Examples

The examples in the Rumoca repository use pinned library versions declared in examples/modelica_dependencies.toml. Fetch them once:

cargo xtask repo modelica-deps ensure

This downloads MSL and the CogniPilot Modelica Models (CMM) into target/. The repository’s committed examples/rumoca-workspace.toml references those directories with workspace-relative paths, so simulation and completion work for every contributor without machine-specific configuration.

Library-Backed Scenarios

The MSL and CMM-backed scenarios live with the other batch runs under examples/simulation/, but they belong conceptually here because the main thing they demonstrate is package/source-root resolution.

KalmanFilterStepTest

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.kalman_filter_step_test.toml

PIDMSL

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.pidmsl.toml

SwitchedRLC_MSL

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.switched_rlc_msl.toml

For a local book build, run cargo xtask repo modelica-deps ensure before building if you want these blocks to compile against the pinned libraries.

Packages and within

Library code is organized as Modelica packages, either as a single .mo file or as a directory tree with package.mo files. The within clause at the top of a library file records which package the file belongs to. When you add a library’s root directory as a source root, all of its classes become resolvable by their fully qualified names.

Practical Tips

• Prefer pinned library versions over whatever happens to be checked out; scenario files plus pinned paths give reproducible runs.
• Large library trees compile slower in the browser playground than natively; use the native CLI or VS Code for MSL-heavy work.
• If a name does not resolve, check Troubleshooting — the usual cause is a missing or wrong source root.

Language Support Status

Rumoca implements a substantial and growing subset of the Modelica language. This page describes what you can rely on today and where the edges are. The ground truth is the continuously-run Modelica Standard Library (MSL) quality gate, which compiles and simulates MSL models on every change and blocks regressions (documented in the Rumoca Dev Guide book).

What Works Well

• Core equation modeling — models, blocks, packages, records, functions, extends, modifications, nested components, connect with potential/flow semantics.
• Continuous dynamics — DAE compilation with structural analysis: matching, sorting (BLT), algebraic-loop tearing, and dummy-derivative index reduction for many higher-index systems.
• Events — when/elsewhen, pre, edge, reinit, sample, if-expressions with event generation, assert, terminate.
• Arrays and regular stencil loops — declarations, slicing, common array builtins, and source-proven affine stencils from for-equation domains. Targets that do not have native stencil kernels still receive exact scalar fallback code.
• Initialization — start attributes and initial equation handling.
• Experiment annotations — StopTime, StartTime, Tolerance, Interval, Solver are parsed and used by scenario and browser runs.

Known Limitations

• High-index DAEs: index reduction handles many mechanical-style systems, but some structurally singular formulations (for example the Cartesian pendulum with an explicit length constraint) are still rejected with a structurally singular system diagnostic. Reformulate with generalized coordinates, or watch the release notes — this area is under active development.
• Event-heavy models with implicit solvers: the default (auto) solver currently selects an implicit method that can stall on rapidly switching models. Workaround: --solver rk-like or annotation(experiment(Solver = "rk-like")).
• Arbitrary large libraries: expect better results with explicit examples and pinned packages than with arbitrary unported libraries; the MSL gate tracks exactly which MSL models compile, simulate, and match reference results.
• Direct CLI runs ignore experiment stop time — pass --t-end or use a batch scenario file ([sim] t_end). Browser/live runs are operator-terminated and do not stop at that horizon.

Checking a Specific Model

The fastest way to find out whether your model is supported is to compile it:

rumoca compile MyModel.mo --emit dae-mo

A clean DAE dump means parsing, resolution, type checking, flattening, and DAE lowering all succeeded. Then rumoca sim (or --inspect structure) exercises the structural preparation and the solver. Diagnostics carry source locations and are designed to name the offending equation or variable.

If you hit something the compiler should support, please file an issue with the model at https://github.com/CogniPilot/rumoca/issues.

Command-Line Interface

The main binary is rumoca. Every subcommand supports --help; this page is a guided reference of the surface you will actually use.

rumoca <COMMAND>

Commands:
  compile      Compile a Modelica file
  sim          Compile and simulate a model/scenario (subcommands: check, init, bench)
  fmt          Format Modelica files
  lint         Lint Modelica files
  completions  Print shell completion scripts
  targets      List built-in code generation targets and their capabilities
  cache        Inspect or prune the shared Rumoca cache

A global --cache-dir <DIR> overrides the on-disk cache root for any subcommand.

rumoca sim

Direct model run:

rumoca sim path/to/Model.mo --model Package.Model --t-end 10

Scenario run:

rumoca sim -c path/to/rumoca-scenario.toml

Subcommands:

Interactive viewer ports and scene files come from the [transport.http] section of the scenario, not CLI flags.

rumoca compile

compile separates three concerns into three flags:

• --emit <stage>-mo|<stage>-json dumps an intermediate representation. Stages are ast, flat, dae, solve (solve has no Modelica form, so only solve-json exists).
• --target <id|dir|file.jinja> runs code generation: a built-in target id (see rumoca targets), a directory containing target.toml, or a raw Jinja template.
• --phase <ast|flat|dae|solve> picks which IR a raw .jinja template receives (default dae). Ignored for built-in and directory targets, which declare their own IR.

rumoca compile Model.mo --emit dae-mo                 # DAE as Modelica, to stdout
rumoca compile Model.mo --emit solve-json -o out.json # solver IR as JSON
rumoca compile Model.mo --target c-ode -o out/      # checked Solve target
rumoca compile Model.mo --target my.jinja --phase flat

compile shares --model, --source-root, --inspect, and --at with sim, and adds -v/--verbose for friendly per-phase progress lines.

rumoca fmt and rumoca lint

See Formatter and Linter.

rumoca targets

Prints the built-in code generation targets with their consumed IR stage, generation mode, deployment class, readiness level, and per-feature support columns, including native/scalar tensor support such as matmul, linsolve, elem, and stencil. See Targets and Templates.

rumoca cache

Compilation artifacts are cached under the platform cache directory.

rumoca cache status          # size and entry counts
rumoca cache prune           # remove oldest cache files until under a size limit

Both accept --root (or the global --cache-dir) to operate on a non-default location.

rumoca completions

rumoca completions bash   # or zsh, fish, ...

Environment

• MODELICAPATH — :-separated library roots, appended after explicit --source-root flags.

Rumoca deliberately has no behavior-changing RUMOCA_* environment variables: every knob is a documented CLI flag or scenario key.

Python API

Rumoca ships a first-class, typed Python API. Create a Session once, compile models through it, and everything — metadata, simulation, code generation, symbolic export — hangs off the returned Model. Returned objects are real typed classes (full autocomplete), never JSON you have to json.loads.

pip install rumoca            # core
pip install rumoca[data]      # + numpy/pandas (result.to_numpy / to_dataframe)
pip install rumoca[plot]      # + matplotlib (result.plot)
pip install rumoca[casadi]    # + casadi   (model.to_casadi)
pip install rumoca[jax]       # + jax      (model.to_jax)
pip install rumoca[all]       # everything, incl. the %%modelica magic

Optional dependencies are genuinely optional: import rumoca never needs them, and a method that does (e.g. result.plot()) raises a message naming the extra to install.

Quick start

Run these examples from the Rumoca repository root. They use the tracked SympyDecay model, whose state is x and whose tunable decay rate is k.

from pathlib import Path
import rumoca as rm

session = rm.Session()
m = session.load("examples/models/SympyDecay.mo", model="SympyDecay")
m                                  # repr shows a summary
m.parameters["k"].value            # 0.5
m.parameters["k"].kind             # "tunable"
m.states["x"]                     # VariableInfo

r = m.simulate(t=(0, 10), dt=0.01)
r.plot("x")                        # needs rumoca[plot]
df = r.to_dataframe()              # needs rumoca[data]

Compile from a string with session.loads(source, model=...). Keep and reuse the Session; it owns source roots and compiler caches.

Or load everything from a rumoca-scenario.toml in one call — its model, source roots, and solver settings (paths resolved relative to the file):

scenario = "examples/simulation/rumoca-scenario.sympy_decay.toml"
session, model, config = rm.Session.from_scenario(scenario)
r = model.simulate(t=(0, 10), config=config)   # config carries solver/dt

For build systems and scenario-driven runs, use the session scenario entry point:

session = rm.Session()
run = session.run_scenario(scenario)
run.status          # "completed"
run.output_paths    # generated reports, CSVs, codegen files, debug logs
run.result          # Result for batch simulation scenarios, otherwise None
run.codegen         # CodegenResult for codegen scenarios, otherwise None

Runtime knobs are scenario overrides, not a separate execution surface:

session = rm.Session()
run = session.run_scenario(
    scenario,
    overrides={"t_end": 2.0, "mode": "lockstep", "output": str(Path("results.csv").resolve())},
)

Relative output paths are resolved beside the scenario file; an absolute path as above writes into the current working directory.

Inspecting a model

Views behave like sequences: len(m.states), m.states[0], m.states["x"], for v in m.states, m.states.names. An unknown name raises a KeyError that suggests the closest match.

Parameter sweeps

Tunable parameters change without recompiling, so a sweep compiles once:

for rate in [0.0, 0.05, 0.1]:
    r = m.with_params(k=rate).simulate(t=(0, 10))
    print(rate, r["x"][-1])

m.simulate(params={...}, start={...}) applies overrides for a single run. Overriding a structural parameter (one that affects sizing/instantiation) raises StructuralParamError; it needs a recompile. Dependent parameter values are recomputed when their tunable inputs change.

A structural parameter can be changed by re-instantiating the model — array dimensions and conditional components re-evaluate:

array_model = session.loads("""
model ArrayDecay
  parameter Integer n = 2 annotation(Evaluate=true);
  parameter Real k = 0.5;
  Real x[n](each start=1, each fixed=true);
equation
  for i in 1:n loop
    der(x[i]) = -k*x[i];
  end for;
end ArrayDecay;
""", model="ArrayDecay")
m3 = array_model.with_params(n=3, recompile=True)
r3 = m3.simulate(t=(0, 1), dt=0.1)

Live symbolic export

CasADi and JAX exports consume the same checked, computable Solve program used by compiled backends. A SolveExport exposes xdot = rhs(t, x, u, p), state_names, input_names, and parameter_names:

casadi_model = m.to_casadi()       # rhs is a differentiable ca.Function
jax_model = m.to_jax()             # rhs supports JIT and differentiation

These render the same tested codegen targets used by m.codegen(target), so the live object and the generated files never drift. For writing files instead, use the same model:

written = m.codegen("casadi-ode").save_all("out/casadi")

For build systems that need a single callable operation, use:

session = rm.Session()
written = session.codegen_file(
    "examples/models/GalecCounter.mo",
    "GalecCounter",
    "galec-production",
    "generated/galec",
)

CMake Integration

CMake needs an executable build step, but downstream projects should not invent a second Rumoca CLI. Use the Python interpreter as the executable and call the stable API:

Save the following as CMakeLists.txt in the repository root, then configure and build it with a Python interpreter that has Rumoca installed. It uses the same tracked GalecCounter example as the Python code above:

cmake_minimum_required(VERSION 3.20)
project(RumocaExamples NONE)
find_package(Python3 REQUIRED COMPONENTS Interpreter)

add_custom_command(
  OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/rumoca-codegen.stamp"
  COMMAND "${Python3_EXECUTABLE}" -c
          "import pathlib, rumoca as rm; \
session = rm.Session(); \
session.codegen_file(r'${CMAKE_CURRENT_SOURCE_DIR}/examples/models/GalecCounter.mo', 'GalecCounter', \
'galec-production', r'${CMAKE_CURRENT_BINARY_DIR}/generated/galec'); \
pathlib.Path(r'${CMAKE_CURRENT_BINARY_DIR}/rumoca-codegen.stamp').write_text('ok\\n')"
  DEPENDS "${CMAKE_CURRENT_SOURCE_DIR}/examples/models/GalecCounter.mo"
  VERBATIM
)
add_custom_target(generate_rumoca ALL
  DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/rumoca-codegen.stamp")

Append this target to run the tracked simulation scenario with cmake --build <build-directory> --target run_rumoca_scenario:

add_custom_target(run_rumoca_scenario
  COMMAND "${Python3_EXECUTABLE}" -c
          "import rumoca as rm; rm.Session().run_scenario(r'${CMAKE_CURRENT_SOURCE_DIR}/examples/simulation/rumoca-scenario.sympy_decay.toml')"
  VERBATIM
)

Reusable CMake modules should be thin wrappers around Session.run_scenario and Session.codegen_file; they should not parse scenario TOML or probe private binding names. Create one session per build step and reuse it for all Rumoca work in that step so compiler caches stay warm.

Jupyter: the %%modelica magic

With rumoca[notebook], write Modelica in a cell and get back a typed object:

%%modelica -m Decay --name m
model Decay Real x(start=1); equation der(x) = -x; end Decay;

%%modelica sim ... returns a Result; %%modelica export casadi ... returns a SolveExport. --name VAR also binds the result into the notebook namespace.

Diagnostics

rm.validate(path) / rm.validate_source(src) return a list of Diagnostic objects and never raise on model errors. Session.load / Session.loads / Model.simulate raise the typed hierarchy — RumocaError and its subclasses ParseError, CompileError, SimulationError, StructuralParamError — each carrying a message that teaches.

diagnostics = rm.validate("examples/models/SympyDecay.mo")
invalid = rm.validate_source("model Broken Real x = ; end Broken;")

The required Linux x86_64 installed-wheel check is cargo xtask python test --require-guide; it also runs in cargo xtask verify full. It builds the current source unless --wheel supplies a wheel, and uses --python to select an interpreter provisioned with the optional dependencies above. The same check executes these Python and notebook examples, both CMake build steps, and the four linked SympyDecay codegen scenarios. Other wheel platforms run the base API checks without claiming this guide qualification.

VS Code Extension

The Rumoca Modelica extension provides language support (diagnostics, completion, hover, semantic highlighting), simulation commands, scenario settings, and result viewers — all inside VS Code.

Install it from the marketplace. It bundles a rumoca-lsp language server, so no separate install is needed. To use your own server build instead (for compiler development), enable rumoca.useSystemServer and put rumoca-lsp on PATH.

Modelica Source Roots

Modelica package paths for a workspace are configured with rumoca-workspace.toml.

Place the file at the workspace root, or in a subdirectory when the roots only apply below that directory:

source_roots = [
  "../target/msl/ModelicaStandardLibrary-4.1.0",
  "../target/cmm/CMM-a642c381",
]

The Rumoca repository examples include examples/rumoca-workspace.toml. Run cargo xtask repo modelica-deps ensure first so those target directories exist.

Settings Panel

Open Rumoca settings from the toolbar or command palette (Open Rumoca Settings). Shared package roots belong in rumoca-workspace.toml; Scenario Source Root Paths edits the active rumoca-scenario.toml for paths that only apply to that one configured run.

The picker stores workspace-relative paths whenever the selected folder is inside the workspace.

Inlay Hints

Array-dimension hints (for example [2x3] after a component name) are always shown. Parameter-name hints, which label call arguments with the callee’s parameter names, are controlled by rumoca.inlayHints.parameterNames and are off by default. Set it to one of:

• none (default): no parameter-name hints.
• userFunctions: hints only for calls to user-defined functions, using each function’s declared parameter names.
• all: hints for user-defined functions and for known builtin functions.

Single-argument calls (such as der(x)) never receive a parameter-name hint, and each hint is anchored at the start of the full argument, including any leading unary operator. Changing the setting takes effect after the language server restarts.

Running Models

Run actions operate on rumoca-scenario.toml scenario files. The extension contributes rumoca-scenario.toml as a TOML language extension, so normal TOML editor features attach while Rumoca still activates for rumoca-scenario.toml and .mo files. Use Create Rumoca Scenario to generate a rumoca-scenario.toml for a model; the runnable source of truth is always the scenario file.

For interactive examples, open the rumoca-scenario.toml scenario, then use Play. For batch simulation, the results panel opens inside VS Code.

Diagnostics

Compiler and simulation diagnostics are reported to the Problems panel when a source span is available. CLI output uses rich terminal diagnostics; VS Code diagnostics use the raw file/range information without terminal wrappers.

Compile and instantiate errors appear in the Problems panel as well, not only in the compile-status lens above a model. When a model fails its full compile because a component’s class cannot be instantiated, the error is anchored at the enclosing model or class name in your open document, and the originating library location is included in the diagnostic message so the detail is not lost. This means the failure you would otherwise see only when simulating is now visible as a squiggle while you edit.

Documentation Access

The user guide and the contributor guide are available directly from the command palette (Open Rumoca User Guide / Open Rumoca Dev Guide Guide), so you do not need to hunt for URLs. The same live-example pages work in any browser.

Web Playground

The browser playground runs the full Rumoca compiler in WebAssembly: a project file tree, Monaco editors with LSP support, simulation with plots, code generation, and package archive loading — no install required.

https://cognipilot.github.io/rumoca/

It is useful for small models, quick experiments, sharing reproductions in bug reports, and demos.

Runnable Blocks in This Book

The runnable code blocks throughout this book (look for the ▶ Simulate button) use the same WASM package as the playground, embedded as focused mini editors:

• the same compiler, solvers, and diagnostics as the native CLI,
• Monaco-based editing with Rumoca’s completion, hover, and error checking,
• inline plots, DAE views, and per-example visualizations.

The first run on a page downloads the WASM compiler; afterwards it is cached by the browser. Models honor their experiment annotation (StopTime, Interval, Tolerance, Solver).

Choosing Plot Channels

After a simulation completes, the results panel plots the channels selected by the active view. Open the panel’s visualization settings (the gear button) to edit a timeseries view. Each view offers two ways to choose channels, and they stay in sync:

• A channel picker: a scrollable checkbox list of every available channel, grouped into States and Outputs (outputs are the non-state algebraic variables). Toggling a checkbox adds or removes that exact channel from the plot.
• A free-text Y field for advanced entry, including the wildcards *states (all states), *outputs (all outputs), and *all (every channel). Wildcards can be mixed with explicit channel names.

The picker appears only once a run has produced channels; before then the Y field still accepts names and wildcards. The same selection is stored per model in the scenario file under [[plot.views]] (see Scenario TOMLs), so a saved view reopens with the same channels on the next run.

Limitations

• Large package trees compile more slowly than native builds.
• Browser storage and worker memory limits matter for full MSL-sized projects.
• Native interactive examples may have more solver/backend options than the browser build.

For larger models or external package development, prefer the native CLI or the VS Code extension.

Formatter and Linter

Rumoca ships a formatter and a linter for Modelica source. Both accept files or directories (defaulting to the current directory) and are also exposed through the VS Code extension and LSP.

rumoca fmt

rumoca fmt                 # format the current directory in place
rumoca fmt src/ Model.mo   # format specific paths
rumoca fmt --check         # report differences without writing (CI-friendly)

Profiles

rumoca fmt --profile canonical

Individual rules can be toggled regardless of profile:

--coverage reports how much of the eligible source trivia the formatter rules cover, without writing changes.

rumoca lint

rumoca lint                       # lint the current directory
rumoca lint --min-level warning   # filter: help | note | warning | error
rumoca lint --warnings-as-errors  # CI gating
rumoca lint --disable-rule naming-convention  # repeat for other rule names
rumoca lint --max-messages 50

Lint diagnostics use the same rich terminal rendering as compiler diagnostics, with source spans. In VS Code they appear in the Problems panel.

Running Simulations

Rumoca offers two ways to run a simulation. They share the same compiler and runtime; the difference is where the configuration lives.

Direct Runs

Point rumoca sim at a .mo file for quick, one-off runs:

rumoca sim Ball.mo --model Ball --t-end 10 --solver rk-like

Everything is a CLI flag: --t-end (default: the model’s experiment(StopTime), else one time unit after experiment(StartTime), where the run starts), --dt, --solver, --source-root, --output. The result is an HTML report with interactive plots of every variable (default <MODEL>_results.html).

Direct runs are great while developing a model. As soon as a run has settings worth repeating, switch to a scenario.

Scenario Runs

A scenario is a rumoca-scenario.toml file colocated with the model it runs:

rumoca sim -c examples/simulation/rumoca-scenario.ball.toml

The scenario records the model file and name, simulation settings, plots, viewer/transport configuration, and source roots — one runnable thing per file. This keeps the CLI, VS Code, and the playground aligned: the play button runs the scenario instead of guessing solver and source roots from a bare .mo file.

See Scenario Files for the format.

Repository Simulation Scenarios

Runnable simulation scenarios live under examples/simulation/. These are the same files opened by the playground and VS Code scenario GUI.

SympyDecay

A one-state exponential decay, useful as a minimal solver and codegen smoke:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.sympy_decay.toml

rumoca sim -c examples/simulation/rumoca-scenario.sympy_decay.toml

Ball

The bouncing ball event model from Events and Discrete Behavior, saved as a reusable scenario:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.ball.toml

rumoca sim -c examples/simulation/rumoca-scenario.ball.toml

SwitchedRLC

A compact circuit-style example with switching behavior:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.switched_rlc.toml

rumoca sim -c examples/simulation/rumoca-scenario.switched_rlc.toml

NeuralODETensor

A tensor-shaped Neural ODE with matrix-vector network layers and a phase portrait plot:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.neural_ode_tensor.toml

rumoca sim -c examples/simulation/rumoca-scenario.neural_ode_tensor.toml

See Neural ODEs for the model structure and parameter-count formula.

NeuralPredatorPrey

A Lotka-Volterra-style Neural ODE whose tensor network predicts per-capita growth rates for prey and predator populations:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.neural_predator_prey.toml

rumoca sim -c examples/simulation/rumoca-scenario.neural_predator_prey.toml

NeuralLatentOscillator

A larger latent Neural ODE with eight hidden states, a learned residual vector field, a learned decoder, and about 19k trainable parameters by default:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.neural_latent_oscillator.toml

rumoca sim -c examples/simulation/rumoca-scenario.neural_latent_oscillator.toml

NeuralODEBackprop

A native backpropagation training demonstration. The model makes weights into states, computes mini-batch vector-field loss, and integrates gradient-descent updates directly:

// rumoca-live-scenario: ../repo-examples/simulation/rumoca-scenario.neural_ode_backprop.toml

rumoca sim -c examples/simulation/rumoca-scenario.neural_ode_backprop.toml

Driving a Session from Another Process

rumoca sim --serve-stdio steps a session under external control, one JSON command and event per line; see Driving a Session over stdio.

What a Run Produces

• Batch runs write an HTML report with time-series plots of all variables, simulation details (solver, tolerances, timing), and any termination message (terminate(...) in the model).
• Interactive runs (scenarios with transports/viewer sections) launch the runtime with a browser viewer and input routing instead — see Interactive Simulation.

Benchmarking

rumoca sim bench measures compile time, preparation time, and hot simulation throughput separately — useful when you care about iteration speed on a large model or are comparing solver settings:

rumoca sim bench Ball.mo --model Ball
rumoca sim bench -c rumoca-scenario.toml

Caching

Compilation artifacts are cached under the platform cache directory, so repeated runs of unchanged models skip recompilation. rumoca cache status shows usage; rumoca cache prune trims it. --cache-dir overrides the location for any command.

Scenario Files (rumoca-scenario.toml)

Rumoca scenarios are plain TOML files and the preferred way to run repeatable simulation and code generation jobs. They follow a filename convention — rumoca-scenario.toml for the default scenario and rumoca-scenario.<profile>.toml for named profiles (such as rumoca-scenario.f16.toml or rumoca-scenario.bench.toml) — and live next to the model they operate on. The filename is the editor/discovery hook; the required [rumoca] marker section is the authoritative declaration.

Each scenario describes one runnable thing. That keeps VS Code, the CLI, and the playground aligned: the play button runs the active scenario instead of guessing from a .mo file.

Getting Started

rumoca sim init > rumoca-scenario.toml    # commented starter template
rumoca sim check -c rumoca-scenario.toml  # validate without running
rumoca sim -c rumoca-scenario.toml        # run

A Minimal Batch Scenario

[rumoca]
version = "1"
task = "simulate"

[model]
file = "../models/Ball.mo"
name = "Ball"

[sim]
solver = "rk-like"
t_end = 10.0
atol = 1e-6
rtol = 1e-6

[[plot.views]]
id = "states_time"
title = "States vs Time"
type = "timeseries"
x = "time"
y = ["x", "v"]

Paths are resolved relative to the rumoca-scenario.toml file.

Section Reference

[rumoca] (required)

The marker section. version = "1" declares the schema version; task = "simulate" runs the model, task = "codegen" renders a target into an output directory.

[model] (required)

[model]
file = "MyVehicle.mo"   # relative to this scenario
name = "MyVehicle"      # top-level class to compile

Use the top-level source_roots key for package dependencies needed by this scenario:

source_roots = ["../modelica_libraries"]

Workspace-wide library paths (MSL, CMM) belong in editor settings, not in every scenario; scenario source_roots are for paths specific to this run.

[sim]

[sim]
dt = 0.01          # simulation timestep [s]
t_end = 10.0       # batch/results-panel output horizon
atol = 1e-6        # optional absolute solver tolerance
rtol = 1e-6        # optional relative solver tolerance
solver = "auto"    # auto | bdf | rk-like
output = "results.html"
mode = "realtime"  # optional pacing, see below

mode selects schedule pacing:

The default is lockstep when external coupling is configured and realtime standalone.

t_end terminates batch/results-panel simulations. Scheduled and browser-live simulations ignore it as a stop condition and extend the solver horizon while they run; stop those runs explicitly with their configured quit signal, the viewer stop control, or an interrupt. Interactive-only scenarios can omit t_end.

[[plot.views]]

Each view adds a plot to the batch report:

[[plot.views]]
id = "states_time"
title = "States vs Time"
type = "timeseries"
x = "time"
y = ["x", "v"]

The y list selects which channels a timeseries view plots. Entries are either exact channel names (any state or output/algebraic variable) or one of the wildcards below, and the two forms can be mixed in a single list:

For example, to plot all states plus one specific output:

[[plot.views]]
id = "states_and_energy"
title = "States and Energy"
type = "timeseries"
x = "time"
y = ["*states", "*outputs", "energy"]

Duplicates are removed and unknown names are ignored when the view is rendered, so combining a wildcard with an explicit name that the wildcard already covers is harmless. The same y selection can also be built interactively from the results panel’s channel picker (see the playground guide); saving that picker writes the selection back to this section.

[transport.*] — external viewer and coupling

HTTP and WebSocket transports serve an external browser viewer surface:

[transport.websocket]
port = 8081

[transport.http]
port  = 8080
scene = "my_scene.js"  # 3D scene, relative to this scenario

UDP is only needed when coupling to an external process:

[transport.udp]
listen = "0.0.0.0:4244"
send   = "127.0.0.1:4242"

[external_interface]
command = "/path/to/external-process"

[schema], [receive], [send] — external interface coupling

These three sections are all-or-nothing. Provide them to couple via FlatBuffers over UDP; omit all three for standalone mode (gamepad/keyboard drive the model inputs directly).

[schema]
bfbs = ["/path/to/your_schema.bfbs"]

[receive]
root_type = "your.namespace.MotorOutput"

[receive.route]
"motors.m0" = { to = "model:omega_m1", scale = 1100.0 }
"armed"     = { to = "local:armed" }

[send]
root_type = "your.namespace.SimInput"

[send.route]
"gyro.x" = { key = "gyro_x" }

[locals] — named persistent simulation state

[locals]
throttle = { type = "float", default = 0.0 }
my_flag  = { type = "bool", default = false }

Types are "bool", "float", or "array" (with element and len). default is optional.

[signals], [input] — input routing

Map keyboard, gamepad, and browser inputs onto model input variables for input-enabled simulations. Input routing is independent from viewer panels or external web presentation; [sim].mode controls the clock and viewer/transport sections control where the run is shown. The worked examples are the best reference:

• examples/interactive/quadrotor/rumoca-scenario.acro.toml
• examples/interactive/rover/rumoca-scenario.toml

Validation

rumoca sim check -c <file> validates structure and paths without running. The VS Code extension surfaces the same validation when editing scenario files.

Solvers and Accuracy

Rumoca ships several integration methods behind one --solver flag (CLI), solver key ([sim] in scenarios), or Solver experiment annotation.

The table below is the whole set. A name that is not in it is an error naming the valid ones — rumoca never falls back to a different solver than the one you asked for, because a run that quietly changed integrator would report timings and trajectories for a method you did not select.

Choosing a Solver

Rules of thumb:

• Start with auto.
• Stiff systems — fast and slow dynamics together (chemical kinetics, stiff mechanical contacts, mu-large Van der Pol) — need an implicit solver; explicit methods crawl with tiny steps.
• Models that switch rapidly (relays, hysteresis controllers) currently run most robustly with rk-like; the implicit path can stall with a step size too small error near dense event cascades. Pin it in the model: annotation(experiment(Solver = "rk-like")).

Tolerances and Output Interval

Implicit solvers control local error against relative/absolute tolerances. Sources, in priority order:

1. CLI/scenario settings (--dt, [sim] dt).
2. The model’s experiment annotation: Tolerance and Interval are used by scenario and browser runs.
3. Runtime defaults.

--dt (or [sim] dt) sets the output interval — and the fixed step for the explicit solver path. If omitted, the runtime chooses automatically.

Experiment Annotations

Place this annotation inside a model, before its closing end, for example in the complete decay example:

annotation(experiment(
  StartTime = 0.0,
  StopTime = 30.0,
  Tolerance = 1e-6,
  Interval = 0.01,
  Solver = "rk-like"
));

StopTime, StartTime, Tolerance, Interval, and Solver (also accepted: Algorithm, __Dymola_Algorithm) are parsed from the annotation. Browser/live runs and the playground honor them fully; native direct CLI runs take end time from --t-end, else StopTime, else 1.0 s, and scenario runs from [sim].

Events and the Solver

All solvers cooperate with the event system: zero crossings from when conditions and if-expressions are located, the state is re-initialized (reinit), and integration restarts cleanly at the event instant. Sampled events (sample(t0, period)) are scheduled exactly, not detected.

Try It

The Van der Pol oscillator becomes stiff as mu grows. Compare solvers by editing the annotation — try mu = 1000 with Solver = "bdf" versus Solver = "rk-like" and watch the run time in the status line:

model VanDerPolStiff
  parameter Real mu = 1000.0 "Try 1, 5, 1000";
  Real x(start = 2.0);
  Real y(start = 0.0);
equation
  der(x) = y;
  der(y) = mu * (1 - x^2) * y - x;
  annotation(experiment(StopTime = 3000.0, Solver = "bdf"));
end VanDerPolStiff;

Interactive Simulation

Interactive simulation is regular task = "simulate" with live inputs enabled. The scenario still owns the model, solver, and output/viewer panels; [sim] chooses the clock policy, [input] maps devices into local state, and [signals.model_inputs] routes that state into Modelica input variables. External browser scenes are just one presentation surface for that same run. Interactive runs are operator-terminated: they continue until you press Q, use the viewer’s Stop control, or interrupt the native process. [sim].t_end is a batch/results-panel horizon and does not stop a live run.

Running the Examples

The interactive examples live under examples/interactive/. They use the same scenario settings surface as batch simulations, but enable input routing and viewer presentation.

Quadrotor SIL

// rumoca-live-scenario: ../repo-examples/interactive/quadrotor/rumoca-scenario.acro.toml

Native run:

cargo xtask repo modelica-deps ensure
cargo run -p rumoca --release -- \
  sim -c examples/interactive/quadrotor/rumoca-scenario.acro.toml

Rover

// rumoca-live-scenario: ../repo-examples/interactive/rover/rumoca-scenario.toml

Native run:

cargo run -p rumoca --release -- \
  sim -c examples/interactive/rover/rumoca-scenario.toml

Fixed-Wing SIL

// rumoca-live-scenario: ../repo-examples/interactive/fixedwing/rumoca-scenario.toml

Native run:

cargo xtask repo modelica-deps ensure
cargo run -p rumoca --release -- \
  sim -c examples/interactive/fixedwing/rumoca-scenario.toml

Reusable Booster Autonomous Landing

// rumoca-live-scenario: ../repo-examples/interactive/reusable_booster/rumoca-scenario.toml

Native run:

cargo xtask repo modelica-deps ensure
cargo run -p rumoca --release -- \
  sim -c examples/interactive/reusable_booster/rumoca-scenario.toml

See Reusable Booster Autonomous Landing for the differential-flatness planner, SE_2(3) controller, gimbal/RCS allocation, and visualization details.

The CLI prints the HTTP and WebSocket endpoints for the viewer. Open the HTTP URL in a browser if it does not open automatically. Use release builds — interactive simulation is real-time work.

Quadrotor Scenario Shape

The quadrotor example is the reference shape:

[rumoca]
version = "1"
task = "simulate"

[model]
file = "QuadrotorSIL.mo"
name = "QuadrotorAcro"

[sim]
dt = 0.01
mode = "realtime"
solver = "rk-like"

[input]
mode = "auto"

Input mappings are scenario-level simulation settings, not viewer settings. The same file then declares local state, keyboard/gamepad bindings, and the bridge into the model:

[locals.throttle]
default = 0.0
type = "float"

[input.gamepad.integrators.throttle]
source = "LeftStickY"
write = "throttle"
deadband = 0.1
rate = 0.7
clamp = [0.0, 1.0]

[input.keyboard.keys.Space]
action = "toggle"
state = "armed"
debounce_ms = 500
precondition = "throttle <= 0.05"

[signals.model_inputs]
stick_throttle = "local:throttle"

[signals.model_inputs.armed]
from = "local:armed"
when_false = 0.0
when_true = 1.0

The external 3D browser view is presentation:

[transport.http]
port = 8080
scene = "quadrotor_scene.js"
asset_dir = "../../assets"

[transport.websocket]
port = 8081

[viewer]
status_title = "Quadrotor"
show_armed = true

Open Settings on the scenario in VS Code, the playground, or the user guide to edit the same TOML through typed controls: input enablement, locals, keys, gamepad axes/buttons/integrators, model-input routes, solver clocking, and viewer panels.

How a Scenario Runs Live

A scenario with only [model] and [sim] is a regular simulation. Live input and live presentation are independent additions:

• [transport.http] + [transport.websocket] — serve the browser viewer (scene = "my_scene.js" selects the 3D scene file). This asks for an external web surface.
• [input], [locals], [signals.model_inputs] — named simulation state and routing from input devices to model input variables. Input alone does not select a special viewer.
• [transport.udp] + [schema]/[receive]/[send] — couple an external process (e.g. a flight controller) over FlatBuffers/UDP.

See Scenario Files for each section.

Pacing

The runtime supports three pacing styles via [sim] mode:

Defaults: lockstep when external coupling is configured, otherwise realtime.

Input Routing

Keyboard, gamepad, and browser inputs go through the generic Rumoca input path. The model side is plain Modelica input variables, so the same model runs in batch (inputs from equations or defaults) and live input runs (inputs from devices) without changing the Modelica source.

Viewer Surfaces

Viewer surfaces are separate from input routing. Standard timeseries/scatter panels use the results viewer. External web scenes use the native realtime scheduled loop and the configured scene script. Both are launched from scenarios, so the same configuration works from the CLI, VS Code, and editor tooling.

Driving a Session over stdio

rumoca sim --serve-stdio runs one simulation session as a child process that an external program steers: set inputs, advance time, read values, repeat. It is the entry point for co-simulation masters, hardware-in-the-loop harnesses, and controllers written in another language. The session uses the real solver (the same one batch sim selects), so events, tolerances, and step control behave as in a batch run. The mode is headless: there is no viewer and no network listener.

rumoca sim --serve-stdio Plant.mo -m Plant --input u=0

--solver, --atol, --rtol, --dt, and --t-end have their batch meaning (--t-end is only the initial horizon; the session extends it on demand). Initialization reads every input before the controller has written one, so an input without a default binding needs --input NAME=VALUE (repeatable); otherwise startup fails with [EX002]. --serve-stdio cannot be combined with --config, --inspect, or --output.

Protocol

The controller writes one JSON command per line to the child’s stdin and reads one JSON event per line from its stdout. Each command is answered by exactly one event, in order. Stdout carries only protocol events; compiler and solver diagnostics go to stderr.

On startup the child writes a hello event with the protocol version and a receipt of the execution engine the session selected:

{"event":"hello","protocol_version":1,"engine":{"engine":"cranelift","refusal":null,"declined":[]}}

engine.engine is interpreter, cranelift, or wasm_program. refusal is null unless compiled execution was refused for a request that allowed it, in which case it is no_continuous_states or external_tables and the engine is interpreter. declined lists each kind of program the selected backend refused to compile, with its reason and a count; those programs run in the interpreter. Compilation is lazy, so declined reflects the session at the time the receipt is read.

Commands

An input takes effect on the next advance. set_inputs applies a whole frame atomically and restarts the integrator history once for the frame, so write coupled signals with one set_inputs rather than several set_input commands. Writing a value that is bit-identical to the current input does not restart the integrator at all.

Events

JSON has no representation for NaN or infinity; a non-finite value is written as null.

Errors

A rejected or failed command is answered with an error event and the session stays usable. code is a stable diagnostic code:

Exit status

Example

-> {"command":"set_input","name":"u","value":1.0}
<- {"event":"ok","time":0.0}
-> {"command":"step","dt":0.001}
<- {"event":"ok","time":0.001}
-> {"command":"state"}
<- {"event":"state","time":0.001,"values":{"u":1.0,"y":0.003996}}
-> {"command":"close"}
<- {"event":"closed","time":0.001}

Closed loop with an external controller

Given Plant.mo:

model Plant
  input Real u;
  Real y(start = 0, fixed = true);
equation
  der(y) = -2 * y + 4 * u;
end Plant;

A proportional controller u = 0.8*(setpoint - y) that lives in the controlling process closes the loop through get, set_input, and step:

import json, subprocess

child = subprocess.Popen(
    ["rumoca", "sim", "--serve-stdio", "Plant.mo", "-m", "Plant", "--input", "u=0"],
    stdin=subprocess.PIPE, stdout=subprocess.PIPE, text=True,
)
assert json.loads(child.stdout.readline())["event"] == "hello"

def call(**command):
    child.stdin.write(json.dumps(command) + "\n")
    child.stdin.flush()
    return json.loads(child.stdout.readline())

setpoint, dt = 1.0, 0.1
for _ in range(10):
    y = call(command="get", name="y")["value"]
    call(command="set_input", name="u", value=0.8 * (setpoint - y))
    call(command="step", dt=dt)
print(call(command="get", name="y")["value"])  # 0.6143..., approaching 8/13
call(command="close")

Reusable Booster Autonomous Landing

This live example starts a reusable orbital booster on a moving drone ship. Press Space to ignite the center engine. The stage performs a six-second ascent burn, coasts through apogee, captures its descending state, and plans a return to the deck. A differential-flatness planner generates the landing reference, an SE_2(3) GPS/IMU error-state Kalman filter estimates the vehicle state, and a geometric controller drives an illustrative hybrid allocation between the center-engine gimbal and RCS.

// rumoca-live-scenario: ../repo-examples/interactive/reusable_booster/rumoca-scenario.toml

The standard live widget opens the Modelica source and TOML scenario in its editor. Open Tune environment and controller over the scene to change wind, gust intensity, deck heave/roll/pitch, or the four controller-gain multipliers while the simulation is running. Those sliders write scenario locals directly to Modelica inputs, so they do not recompile or reset the run. This makes it possible to search for the landing-failure boundary. The editor’s Settings panel remains available for mission constants stored in TOML [parameters].

Press R at any time to restore the vehicle and restart at t = 0. The live run has no preset end time; it continues after touchdown so the landed state and deck motion remain observable until you stop or reset it. The upper-right badge shows mission phase, wind and gust speed, and deck roll/pitch. T switches between realtime and fast pacing, F toggles fullscreen, and Q stops the run.

For a native run:

cargo xtask repo modelica-deps ensure
cargo run -p rumoca --release -- \
  sim -c examples/interactive/reusable_booster/rumoca-scenario.toml

Reference Planning

Position is the translational flat output. For each axis, the Modelica planner constructs a quintic from the position, velocity, and acceleration captured when the coasting stage reaches 8 m/s downward speed. Its terminal state is the predicted deck CG target along the tilted deck normal, including the target point’s translational and rotational velocity and acceleration:

p(t) = c0 + c1 t + c2 t^2 + c3 t^3 + c4 t^4 + c5 t^5
(p, v, a)(0) = (p0, v0, a0),  (p, v, a)(T) = (pf, vf, af)

Analytic derivatives provide reference velocity and acceleration. The vector a_ref + g e3 determines the flatness-derived reference thrust direction and attitude. The amber Three.js line is the complete polynomial reference; its marker advances along the path while the cyan line records the simulated trajectory. The green line is the navigation estimate and the magenta marker is the latest GPS fix consumed by the filter. Green and magenta wireframes show their respective axis-aligned 3-sigma marginal-variance envelopes at 4x visual scale so the sub-metre posterior remains legible beside a 40.9 m stage. They show per-axis marginal variance, not covariance orientation.

The nominal landing-plan duration is 18 seconds. Guidance continues to hold the terminal reference until a safe physical deck contact is detected; reaching the polynomial end does not shut down the engine above the ship.

A 21-point sampled preflight check covers throttle, tilt, positive vertical thrust, deck clearance, and approximate propellant budget. The scene explicitly reports REFERENCE INFEASIBLE when the unconstrained quintic exceeds those modeled limits; this warning is not a constrained trajectory optimizer.

GPS/IMU SE_2(3) Kalman Filter

The navigation state uses the same group representation as the controller:

Xhat = (phat, vhat, Rhat) in SE_2(3)

At 10 Hz, piecewise-constant gyroscope and accelerometer measurements drive LieGroups.SE23.Quat.exp_mixed. The function implements the closed-form mixed invariant propagation described in Purdue’s On Closed-Form Preintegration for a Class of Mixed-Invariant Systems in SE_n(3):

Xhat(k+1) = exp(M dt) Xhat(k) exp(N dt)

The call includes gravity as a world-frame right increment and the nilpotent position-velocity coupling matrix. This is the library implementation from the pinned CogniPilot modelica_models checkout; the example does not duplicate the exponential.

A 5 Hz synthetic differential-GPS position fix performs the Kalman correction. The filter propagates a coupled 9 x 9 position, velocity, and attitude-error covariance. Its inertial transition includes position-velocity coupling and the sensitivity of specific-force integration to attitude error, allowing a position fix to correct velocity and attitude through cross-covariance. The GPS update uses Joseph form for covariance robustness and a multiplicative attitude correction. The model intentionally omits sensor biases, time delay, outlier rejection, Earth rotation, and geodetic effects.

Geometric Tracking

The controller represents estimated pose and velocity together on SE_2(3):

Xhat = (phat, vhat, Rhat)
Xref = (pref, vref, Rref)
xi   = Log(Xhat^-1 Xref)
delta = J_left(xi) K xi

inverse, product, log_map, and left_jacobian come directly from the pinned CogniPilot modelica_models LieGroups.SE23.Quat package. Translation feedback corrects the flatness acceleration, and an SO(3) inner loop aligns the booster axis with the resulting force vector. The inner loop tracks the geodesically rate-limited attitude command with per-axis gain vectors and includes rigid-body gyroscopic and rotating-reference feedforward terms.

The geometric controller runs at 20 Hz with a zero-order hold, while the rigid-body, engine, drag, mass, and contact dynamics remain continuous. This matches a sampled flight-control implementation and avoids reevaluating the Lie-group maps at every adaptive integrator stage.

The control allocator prioritizes center-engine thrust vectoring for pitch/yaw moment. RCS receives only the saturated residual and the requested axial moment. RCS availability fades in between 80 and 120 m engine-plane altitude, its aggregate force/moment has an 80 ms first-order response, and its representative cold-gas mass flow contributes to vehicle mass depletion. Engine and RCS plume geometry, light, and opacity follow actual actuator state. Landing-foot tangential damping is regularized by a Coulomb cone, so friction is capped by normal load and vanishes for an unloaded foot.

Embeddable eFMI Control Law

The continuous vehicle and the CMM-based Lie-group preprocessing remain in the simulation model. The force/moment feedback law behind them is a separate deep module with four vector inputs: translational terms, rotational terms, gain scales, and principal inertia. The 20 Hz simulation controller and the export model call the same functions, so the generated code does not maintain a second copy of the control equations.

The export model is fixed-sample and discrete, which makes it admissible for Rumoca’s galec-production target. Select Generate .alg below to inspect the eFMI Algorithm Code, then Generate C/H to inspect the corresponding C99 Production Code.

// rumoca-live-scenario: ../repo-examples/interactive/reusable_booster/rumoca-scenario.controller-galec-production.toml

The native packaging command emits both a directory-form eFMU and a matching .efmu archive with Algorithm Code and Production Code representations:

cargo xtask repo modelica-deps ensure
cargo run -p rumoca -- \
  compile examples/interactive/reusable_booster/ReusableBoosterLanding.mo \
  --model ReusableBoosterEmbeddedControlLaw \
  --source-root target/cmm/CMM-a642c381 \
  --target galec-production \
  --output examples/interactive/reusable_booster/gen/control_law_efmu

This exported seam starts after the geometric adapter has computed the world-frame translational correction, attitude error, and reference body rate. Exporting the full estimator and Lie-group adapter would additionally require GALEC projection support for the indexed array outputs and matrix operations used by those library functions.

Dryden Wind and Moving Deck

The default mean wind is 5 m/s. Its direction swings by 15 degrees on a 10-second period, while longitudinal, lateral, and vertical deterministic inputs pass through standard Dryden-shaped filters:

Hu(s)   = sigma sqrt(2 Lu / (pi V)) / (1 + Lu s / V)
Hv,w(s) = sigma sqrt(L / (pi V)) (1 + sqrt(3) L s / V) / (1 + L s / V)^2

Band-limited multi-sine inputs replace ideal white noise so a browser run is deterministic and repeatable. The resulting three-axis gust is rotated with the mean-wind direction, added to the mean wind, and subtracted from inertial velocity before the body-axis drag calculation. The numeric indicator, windsock, and deck arrow use this same changing total wind vector.

The drone ship has analytic heave, roll, and pitch kinematics. The planner predicts the terminal CG position along the tilted deck normal and matches that point’s velocity and acceleration, including rotational motion. For each landing foot, the contact model transforms its position into the deck frame and computes

penetration = max(0, -z_foot_deck)
vrel = vfoot - (vdeck + omega_deck x (pfoot - pdeck))

Normal spring/damping acts along the tilted deck normal and tangential damping uses the relative deck-point velocity. Contact is limited to the physical deck bounds. The same Modelica deck position and quaternion drive the Three.js ship, so visible wave motion and collision geometry remain synchronized. The telemetry badge reports deck roll and pitch beside the changing wind and gust values.

Touchdown and Failure State

At landing-phase contact, at least three feet must carry meaningful normal load and the CG projection must remain at least 0.25 m inside their support polygon. The maximum foot speed normal to the deck must be at most 2.0 m/s, tangential speed at most 1.5 m/s, and booster tilt at most 10 deg from the local deck normal. Safe first contact cuts the engine and starts a bounded suspension- settling interval. Stable support during that interval latches LANDED; unsafe kinematics or failure to establish support latches a crash. The scene reports foot count, support margin, and separate normal/tangential speeds.

Vehicle and Contact Model

The visual and physical model use a representative 40.9 m stage and 3.66 m diameter. The plant includes an illustrative landing mass and fixed landing inertia/CG, a generic center-engine propulsion model, propellant depletion, thrust response, anisotropic aerodynamic drag, quaternion 6-DOF dynamics, and four spring-damper landing contacts. The scene includes an interstage, nine-engine mount, deployed grid fins, articulated landing legs, and a scale-matched autonomous drone ship.

After pressing Capture, mouse movement orbits and tilts the scene, the middle or right mouse button pans, and the wheel zooms. Press Escape to release capture.

The geometry, landing mass, inertia, RCS force, throttle floor, gains, and contact coefficients are mutually consistent educational values. They are not intended to reproduce any particular operational vehicle.

Inspecting and Debugging Models

When a model misbehaves — fails to compile, fails to initialize, or produces wrong dynamics — Rumoca gives you structured views into every stage of the pipeline. All of these work with both rumoca compile and rumoca sim.

Dump an Intermediate Representation

--emit prints the model as the compiler sees it after each stage:

rumoca compile Model.mo --emit flat-mo          # to stdout
rumoca compile Model.mo --emit dae-json -o m.json

ast-json is the one stage before name resolution. It dumps the file you pass exactly as the parser read it: --model and --source-root are accepted but ignored, so the dump covers the whole file and stays available while a library the file references is missing or does not resolve. Every later stage needs the model’s whole class closure, so flat-mo and beyond still report unresolved names. For the Modelica Standard Library that means naming the directory that holds Modelica/, ModelicaServices/, and Complex.mo as one --source-root, since Modelica references both of the other top-level packages.

--emit belongs to rumoca compile; a scenario file has no equivalent key and rumoca sim --config does not dump stages. To inspect the model a scenario runs, pass its [model] entries and source_roots to compile:

rumoca compile models/Ball.mo --model Ball --source-root ../modelica_libraries --emit ast-json

Reading flat-mo answers “what did my modifications and connects actually produce?”. Reading dae-mo answers “what equation system is the solver given?” — the live examples in this book expose the same view through their Show DAE button.

The JSON dumps are for debugging and for replay with the same compiler build, not a stable interface. dae-json carries a schema_version that changes as the IR evolves, and its reader accepts only the current version; flat-json has no version. For a model’s inputs and outputs, read the modelDescription.xml of an fmi2 or fmi3 export instead.

In dae-json, a variable’s role is its runtime classification (state, algebraic, output, discrete_value, …), not the declared prefix: a declared output that is a state has role state, and a discrete one has role discrete_value. The declared prefix is declared_causality (none, input, or output), kept for nested components too, while causality is input or output only at the top level. FMI exports mark a nested declaration’s prefix with a per-variable annotation and keep its causality local:

<!-- FMI 2, inside the ScalarVariable after its type element -->
<Annotations><Tool name="rumoca"><DeclaredCausality value="output"/></Tool></Annotations>
<!-- FMI 3, first child of the variable -->
<Annotations><Annotation type="rumoca.declaredCausality">output</Annotation></Annotations>

Structural Analysis

rumoca compile Model.mo --inspect structure

Prints the structural preparation of the system: the matching between equations and unknowns, the block lower-triangular (BLT) ordering, simultaneous (coupled) blocks, and tearing decisions. This is the first place to look when compilation fails with structurally singular system — it names the unmatched equations and unknowns.

Numerical Evaluation at a Point

rumoca sim Model.mo --inspect eval
rumoca sim Model.mo --inspect eval --at "x=1.5,v=0@2.0"

Evaluates all solver values and state derivatives at a point and names any non-finite results — the fastest way to find the division-by-zero or domain error behind a NaN. With no --at, it evaluates at the initial state (which also discovers the state names for you).

The --at syntax is <name=value,...@t>: states by name, unset states keep their initial values, time after @ (default 0).

Jacobian Analysis

rumoca sim Model.mo --inspect jacobian --at "x=1.0@0"

Prints the dense state Jacobian at a point and flags singular columns and zero pivots — useful for diagnosing initialization failures and stiff or degenerate dynamics.

NaN Tracing

When a simulation fails with a non-finite value, rumoca sim automatically re-runs with NaN tracing to locate the offending variables, so the diagnostic names the variable instead of just reporting a solver failure.

Performance

rumoca sim bench Model.mo            # compile / prepare / hot-loop timing
rumoca cache status                  # compilation cache usage

Verbose Compilation

rumoca compile Model.mo --target c-ode -o out -v

-v prints friendly [rumoca] Phase ... progress lines, which localizes slow or failing phases on large models.

A Neural Network in Modelica

A neural network is a composition of ordinary array functions. The shared NeuralNetwork.mo package keeps those functions separate from the plant, policy parameters and learning algorithm. CartPole uses a single neuron; CanyonRun.ContinuousControl.CanyonFly uses two hidden layers and the learned CanyonRun weights.

Rows Are Neurons

dense(W, b, x) calculates y = W*x + b. Each row of W belongs to one output neuron, so its second dimension must match the input length. The array dimensions in each function signature express those contracts directly. The loops make the accumulation easy to inspect.

activate(x) applies tanh element by element. Its outputs lie between −1 and 1. forward composes two dense/tanh hidden layers and a final linear layer; the controller applies its final activation and physical units.

within;
package NeuralNetwork "Small neural-network functions shared by the RL examples"
  function dense "Affine layer: each row of W produces one output"
    input Real W[:, :];
    input Real b[size(W, 1)];
    input Real x[size(W, 2)];
    output Real y[size(W, 1)];
  algorithm
    for i in 1:size(W, 1) loop
      y[i] := b[i];
      for j in 1:size(W, 2) loop
        y[i] := y[i] + W[i, j]*x[j];
      end for;
    end for;
  end dense;

  function activate "Smooth bounded activation, element by element"
    input Real x[:];
    output Real y[size(x, 1)];
  algorithm
    for i in 1:size(x, 1) loop
      y[i] := tanh(x[i]);
    end for;
  end activate;

  function clip "Elementwise observation bound before inference"
    input Real x[:];
    input Real limit;
    output Real y[size(x,1)];
  algorithm
    for i in 1:size(x,1) loop
      y[i] := max(-limit,min(limit,x[i]));
    end for;
  end clip;

  function forward "Two tanh hidden layers and a linear output layer"
    input Real W1[:, :];
    input Real b1[size(W1, 1)];
    input Real W2[:, size(W1, 1)];
    input Real b2[size(W2, 1)];
    input Real W3[:, size(W2, 1)];
    input Real b3[size(W3, 1)];
    input Real x[size(W1, 2)];
    output Real y[size(W3, 1)];
  protected
    Real h1[size(W1, 1)];
    Real h2[size(W2, 1)];
  algorithm
    h1 := activate(dense(W1, b1, x));
    h2 := activate(dense(W2, b2, h1));
    y := dense(W3, b3, h2);
  end forward;

  function normal "Box-Muller sample with a small, explicit integer RNG state"
    input Integer seed;
    output Real value;
    output Integer nextSeed;
  protected
    Integer first;
    Real u1;
    Real u2;
  algorithm
    // Park-Miller products stay exactly representable in a 53-bit number.
    first := mod(16807*seed, 2147483647);
    nextSeed := mod(16807*first, 2147483647);
    u1 := first/2147483647.0;
    u2 := nextSeed/2147483647.0;
    value := sqrt(-2*log(u1))*cos(2*3.141592653589793*u2);
  end normal;

  function perturb "Sample a candidate policy around the best retained weights"
    input Real best[:];
    input Real scale;
    input Integer seed;
    output Real candidate[size(best, 1)];
    output Integer nextSeed;
  protected
    Real noise;
  algorithm
    nextSeed := seed;
    for i in 1:size(best, 1) loop
      (noise, nextSeed) := normal(nextSeed);
      candidate[i] := best[i] + scale*noise;
    end for;
  end perturb;
end NeuralNetwork;

Learning Is Separate from Inference

The Modelica CartPole model calls dense and activate to choose a force. The swing-up learner uses normal to sample candidate gains and refits its search distribution from the two best retained policies. The older upright example uses perturb to propose weights around its best policy. normal uses an explicit Park–Miller state and Box–Muller transform, so learning and resets do not depend on a hidden global generator. Use a seed between 1 and 2147483646.

The Canyon actor uses forward with frozen matrices from its trained policy. The matrices are stored in a separate data package, TrainedCanyonPolicy.mo. They are transposed once during export from the Python trainer’s input-by-output layout to this package’s output-by-input layout. The exported actor does not include the critic or optimizer state.

The current continuous Canyon training route still uses the Python PPO implementation shown in that chapter. Calling a trained actor in Modelica does not mean that PPO’s optimizer has also been ported to Modelica.

CartPole: Learn Swing-Up and Balance

Start with CartPole before Canyon Run. The cart applies one horizontal force to swing a hanging pole up and then balance it. The physics, learner, random-number state and episode resets all run in Modelica. The viewer displays their results; it does not select weights.

Train in the Editor

Press Run below. The initial pole angle is π − 0.01 radians: 0.57° off straight down, with the cart and pole at rest. This small asymmetry is explicit; exactly symmetric down-start performance is not claimed. All seven policy gains begin at zero.

The example starts in Realtime. Choose Fast to learn sooner and switch back to watch the learned policy. It freezes its gains after five simulated seconds of balance. The cart centre must stay within ±2 m. Attempting to leave a rail ends and penalizes that episode, then resets the plant.

// rumoca-live-scenario: ../repo-examples/reinforcement_learning/cartpole/rumoca-scenario.toml

Drag the 3D view to rotate; right-drag to pan; scroll or pinch to zoom. Reset view restores the camera. R restarts training. After learning, Bump left/right requests an actual ±20 N force for 0.3 simulated seconds. A disturbance can still cause failure: the plant resets, retaining the frozen gains. The model bounds the combined policy and bump force to ±30 N.

The diagrams show actual learned gains and observed episode rewards. Fast mode can skip intermediate episodes in the displayed plot; training still evaluates each episode. Realtime pacing and simulation time are separate from rendering.

A Small Policy with an Explicit Physics Prior

A single tanh neuron combines seven clipped features. Three pump the pole’s energy while it hangs down; four emphasize balance near the upright position. The features depend on cart position/velocity, sin/cos of the pole angle and angular velocity. This physical feature map is fixed and documented in the source. All seven gains are learned; no stabilizing gains are preset. This is a small policy-search example, not an end-to-end deep actor.

The learner retains the two best policies across episodes and refits a diagonal Gaussian every 32 evaluations, with a minimum standard deviation of 1. Each candidate gets at most 20 simulated seconds. Its reward favors pole height and sustained balance, penalizes energy error, cart displacement and force, and subtracts 50 on a rail failure. Success requires cos(angle) > 0.98, |angular velocity| < 1 rad/s and |cart position| < 1.2 m for 250 control ticks. The control period is 0.02 s. A one-sample observation/weight delay and the one-candidate distribution-refit delay are explicit in the source.

The shared NeuralNetwork package supplies the dense layer, tanh, clipping and Gaussian samples. PolicySearch.mo owns candidate sampling, elite retention and refitting; CartPoleSwingUp.mo owns the task and force.

One checked native seed-1 run evaluated 84 policies before freezing. A separate 30-second run from the same near-down start, with learning disabled, held the balance criterion from 1.92–30.00 s: 28.08 s. Maximum absolute cart position was 1.498 m, with no rail contacts or resets and unchanged weights. Six separate boundary controls checked both signs, inward noncontacts and outward starts at the limits. These are nominal development results, not held-out, multiple-seed, disturbance or Gymnasium benchmark acceptance.

Native equivalent, from the repository root:

rumoca sim examples/reinforcement_learning/cartpole/CartPoleSwingUp.mo \
  --source-root examples/reinforcement_learning/NeuralNetwork.mo \
  --source-root examples/reinforcement_learning/cartpole/PolicySearch.mo \
  --source-root examples/reinforcement_learning/cartpole/CartPole.mo \
  --model CartPoleSwingUp --solver rk-like

within;
package PolicySearch "Diagonal cross-entropy search over neural policy weights"
  function candidate "Draw one policy; the caller owns the random state"
    input Real mean[:];
    input Real deviation[size(mean,1)];
    input Integer seeds[size(mean,1)];
    output Real weights[size(mean,1)];
    output Integer nextSeeds[size(mean,1)];
  protected
    Real noise;
  algorithm
    for i in 1:size(mean,1) loop
      (noise,nextSeeds[i]) := NeuralNetwork.normal(seeds[i]);
      weights[i] := mean[i]+deviation[i]*noise;
    end for;
  end candidate;

  function retain "Replace only the worst elite when this episode is better"
    input Real elite[:, :];
    input Real scores[size(elite,1)];
    input Real weights[size(elite,2)];
    input Real score;
    output Real updated[size(elite,1),size(elite,2)];
    output Real updatedScores[size(elite,1)];
  protected
    Integer worst;
  algorithm
    updated := elite;
    updatedScores := scores;
    worst := 1;
    for i in 2:size(scores,1) loop
      if scores[i]<scores[worst] then
        worst := i;
      end if;
    end for;
    if score>scores[worst] then
      for j in 1:size(weights,1) loop
        updated[worst,j] := weights[j];
      end for;
      updatedScores[worst] := score;
    end if;
  end retain;

  function fit "Refit each weight's Gaussian from the retained elite policies"
    input Real elite[:, :];
    output Real mean[size(elite,2)];
    output Real deviation[size(elite,2)];
  protected
    Real variance;
  algorithm
    for j in 1:size(elite,2) loop
      mean[j] := 0;
      for i in 1:size(elite,1) loop
        mean[j] := mean[j]+elite[i,j]/size(elite,1);
      end for;
      variance := 0;
      for i in 1:size(elite,1) loop
        variance := variance+(elite[i,j]-mean[j])^2/size(elite,1);
      end for;
      deviation[j] := max(1.0,sqrt(variance));
    end for;
  end fit;

  model Learner "Episode-level search, separate from the fast plant/controller clock"
    parameter Integer randomSeeds[7];
    input Boolean ended;
    input Boolean reached;
    input Real score;
    output Boolean trained(start=false,fixed=true);
    output Integer evaluations(start=0,fixed=true);
    output Integer generation(start=1,fixed=true);
    output Real weights[7](each start=0,each fixed=true);
  protected
    Integer rng[7](start=randomSeeds,each fixed=true);
    Real mean[7](each start=0,each fixed=true);
    Real deviation[7](each start=1,each fixed=true);
    Real elite[2,7](each start=0,each fixed=true);
    Real eliteScores[2](each start=-1e9,each fixed=true);
    Real retained[2,7](each start=0,each fixed=true);
    Real scores[2](each start=-1e9,each fixed=true);
  algorithm
    when ended then
      trained := pre(trained) or reached;
      if not pre(trained) then
        evaluations := pre(evaluations)+1;
        (retained,scores) := retain(pre(elite),pre(eliteScores),pre(weights),score);
        if mod(evaluations,32)==0 then
          (mean,deviation) := fit(retained);
          generation := pre(generation)+1;
        else
          mean := pre(mean);
          deviation := pre(deviation);
          generation := pre(generation);
        end if;
        elite := retained;
        eliteScores := scores;
        if trained then
          weights := pre(weights);
          rng := pre(rng);
        else
          // The first candidate after a refit uses the prior distribution;
          // this explicit event delay keeps search updates separate from sampling.
          (weights,rng) := candidate(pre(mean),pre(deviation),pre(rng));
        end if;
      else
        evaluations := pre(evaluations);
        generation := pre(generation);
        retained := pre(retained);
        scores := pre(scores);
        mean := pre(mean);
        deviation := pre(deviation);
        elite := pre(elite);
        eliteScores := pre(eliteScores);
        weights := pre(weights);
        rng := pre(rng);
      end if;
    end when;
  end Learner;
end PolicySearch;

Compare the Python Upright-Balance Route

This older comparison teaches upright balancing from a small initial tilt; it does not perform the swing-up task above. It uses a four-gain policy and a separate cross-entropy search. The checked replay below belongs to that separate run; it is not presented as new training.

Run the Python Learner

From the repository root, with a native rumoca executable on your path:

python3 examples/reinforcement_learning/cartpole/train.py \
  --rumoca rumoca --output cartpole-run

If you built the CLI locally, use --rumoca target/release/rumoca instead. The output directory must be new. The default run uses seed 1, eight generations of 16 candidate policies, and a 180-second wall-time limit.

The program first measures a zero-force baseline, then prints the best episode length and population mean after each generation. It saves the learned gains, training history and a final evaluation in cartpole-run/result.json. cartpole-run/flight.csv records that final flight’s states and force, with learning disabled.

In a checked native run, the zero-force policy lasted 34 steps. The best policy reached the 500-step limit in generation 2 and again in its final frozen evaluation. The population mean rose from 44.5 steps in generation 1 to 479.9 in generation 8. The run took about 9 seconds on one core; elapsed time depends on the machine and compiler build.

What Is Learned?

The policy is a linear combination of the four normalized states, passed through tanh to bound the force to ±10 N:

scales = [2.4, 2.0, math.radians(12), 2.0]
force = 10 * math.tanh(sum(w*x/s for w, x, s in zip(weights, state, scales)))

The four weights start unknown. The learner uses cross-entropy policy search:

1. Sample 16 sets of gains from a Gaussian distribution.
2. Run one episode for each set and score it by how many steps it lasts.
3. Keep the best quarter and fit the next Gaussian to their gains.
4. Repeat, then evaluate the best policy without changing its gains.

This is reinforcement learning by direct policy search: episode reward selects the gains. It does not use a hand-tuned balancing controller, a value network, physics derivatives or backpropagation. It is smaller than the actor-critic methods used for more complex continuous-control tasks.

Each control step advances the Modelica plant by 0.02 s. An episode ends if the cart leaves ±2.4 m or the pole tilts beyond ±12 degrees; otherwise it is truncated after 500 steps, or 10 s. Every step, including a terminating step, contributes one point. The old episode is scored before resetting the simulator. The policy search retains its gains and random-number state across resets.

Physics and Scope

The model is an upright, frictionless cart with a hinged rigid pole. Its parameters are a 1 kg cart, a 0.1 kg pole, a 0.5 m pivot-to-centre-of-mass distance and gravity of 9.8 m/s². The cart starts at rest, with a 0.03 rad pole tilt. Rumoca integrates the equations; the Python program supplies force and reads the four states.

The equations and failure limits follow the familiar CartPole problem. This example uses a continuous bounded force, a fixed initial state and Rumoca’s solver, so its scores are not Gymnasium CartPole benchmark results. Balancing from this one start also does not prove robustness to disturbances.

Try changing initialAngle in the model, increasing the training population, or evaluating the saved gains from a different starting state. Keep those evaluation results separate from the episodes used to select the gains.

"""Learn four CartPole feedback gains from episode reward, using only stdlib."""
import argparse
import csv
import json
import math
import os
from pathlib import Path
import random
import selectors
import shutil
import subprocess
import tempfile
import time


class Session:
    def __init__(self, binary, scratch, deadline):
        self.deadline = deadline
        self.time = 0.0
        self.buffer = b''
        self.log = (scratch / 'stderr.log').open('wb')
        self.child = subprocess.Popen([binary, 'sim', str(Path(__file__).with_name('CartPole.mo')),
            '-m', 'CartPole', '--solver', 'rk-like', '--serve-stdio',
            '--cache-dir', str(scratch / 'cache')], stdin=subprocess.PIPE,
            stdout=subprocess.PIPE, stderr=self.log)
        self.selector = selectors.DefaultSelector()
        self.selector.register(self.child.stdout, selectors.EVENT_READ)
        try:
            hello = self.read()
            if hello.get('event') != 'hello' or hello.get('protocol_version') != 1:
                raise ValueError('unsupported session handshake')
        except Exception:
            self.close()
            raise

    def read(self):
        while b'\n' not in self.buffer:
            remaining = self.deadline - time.monotonic()
            if remaining <= 0 or not self.selector.select(remaining):
                raise TimeoutError('CartPole session deadline')
            data = os.read(self.child.stdout.fileno(), 65536)
            if not data:
                raise RuntimeError('session ended; inspect stderr.log')
            self.buffer += data
            if len(self.buffer) > 1024*1024:
                raise ValueError('oversized session response')
        line, self.buffer = self.buffer.split(b'\n', 1)
        return json.loads(line)

    def call(self, command, expected='ok'):
        if time.monotonic() >= self.deadline:
            raise TimeoutError('CartPole session deadline')
        self.child.stdin.write((json.dumps(command) + '\n').encode())
        self.child.stdin.flush()
        event = self.read()
        observed = event.get('time')
        if event.get('event') != expected:
            raise ValueError('unexpected session response: ' + str(event))
        if type(observed) not in (float, int) or not math.isfinite(observed):
            raise ValueError('invalid response time')
        return event

    def state(self):
        event = self.call({'command': 'state'}, 'state')
        if abs(event['time'] - self.time) > 1e-9:
            raise ValueError('stale observation')
        values = [event['values'][name] for name in ['x', 'v', 'theta', 'omega']]
        if not all(type(v) in (int, float) and math.isfinite(v) for v in values):
            raise ValueError('invalid physical state')
        return values

    def set_force(self, force):
        if not math.isfinite(force) or abs(force) > 10:
            raise ValueError('force outside policy bounds')
        event = self.call({'command': 'set_inputs', 'inputs': [['force', force]]})
        if abs(event['time'] - self.time) > 1e-9:
            raise ValueError('input update advanced time')

    def reset(self):
        self.set_force(0.0)
        event = self.call({'command': 'reset'})
        if event['time'] != 0:
            raise ValueError('reset did not return time zero')
        self.time = 0.0
        return self.state()

    def step(self, force):
        self.set_force(force)
        target = self.time + .02
        event = self.call({'command': 'advance_to', 'time': target})
        if abs(event['time'] - target) > 1e-9:
            raise ValueError('step time mismatch')
        self.time = event['time']
        return self.state()

    def close(self):
        try:
            if self.child.poll() is None:
                try:
                    self.child.stdin.write(b'{"command":"close"}\n')
                    self.child.stdin.flush()
                    self.child.wait(timeout=3)
                except (BrokenPipeError, subprocess.TimeoutExpired):
                    pass
        finally:
            if self.child.poll() is None:
                self.child.kill()
                self.child.wait(timeout=3)
            self.selector.close()
            self.log.close()
        return self.child.returncode


def force_for(weights, state):
    scales = [2.4, 2.0, math.radians(12), 2.0]
    return 10 * math.tanh(sum(w*x/s for w, x, s in zip(weights, state, scales)))


def episode(session, weights, trace=None):
    state = session.reset()
    for step in range(1, 501):
        force = force_for(weights, state)
        state = session.step(force)
        failed = abs(state[0]) > 2.4 or abs(state[2]) > math.radians(12)
        if trace is not None:
            trace.writerow([session.time, *state, force, failed])
        if failed:
            return step
    return 500


def search(session, rng, generations, population):
    mean, deviation = [0.0]*4, [1.0]*4
    best, best_score, history = mean[:], 0, []
    for generation in range(generations):
        candidates = [[rng.gauss(m, s) for m, s in zip(mean, deviation)]
                      for _ in range(population)]
        scored = [(episode(session, weights), weights) for weights in candidates]
        scored.sort(key=lambda item: item[0], reverse=True)
        if scored[0][0] > best_score:
            best_score, best = scored[0]
        elites = [weights for _, weights in scored[:max(2, population//4)]]
        mean = [sum(w[j] for w in elites)/len(elites) for j in range(4)]
        deviation = [max(.1, math.sqrt(sum((w[j]-mean[j])**2 for w in elites)/len(elites)))
                     for j in range(4)]
        row = {'generation': generation+1, 'best_steps': best_score,
               'population_mean_steps': sum(score for score, _ in scored)/population}
        history.append(row)
        print(json.dumps(row), flush=True)
    return best, history


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--rumoca', default='rumoca')
    parser.add_argument('--output', type=Path, default=Path('cartpole-run'))
    parser.add_argument('--generations', type=int, default=8)
    parser.add_argument('--population', type=int, default=16)
    parser.add_argument('--seed', type=int, default=1)
    parser.add_argument('--wall-seconds', type=int, default=180)
    args = parser.parse_args()
    if args.population < 4 or min(args.generations, args.wall_seconds) < 1:
        parser.error('positive generations/wall budget and population >=4 required')
    args.output.mkdir(exist_ok=False)
    result = {'completed': False, 'seed': args.seed, 'method': 'cross-entropy policy search',
              'scope': 'fixed-start instructional CartPole, not Gymnasium benchmark acceptance'}
    started = time.monotonic()
    with tempfile.TemporaryDirectory(prefix='cartpole-') as scratch_name:
        scratch, session = Path(scratch_name), None
        try:
            session = Session(args.rumoca, scratch, started + args.wall_seconds)
            result['zero_policy_steps'] = episode(session, [0.0]*4)
            weights, history = search(session, random.Random(args.seed), args.generations, args.population)
            with (args.output/'flight.csv').open('w', newline='') as handle:
                writer = csv.writer(handle)
                writer.writerow(['time', 'x', 'v', 'theta', 'omega', 'force', 'failed'])
                result['frozen_steps'] = episode(session, weights, writer)
            result.update(completed=True, weights=weights, training=history)
        except Exception as error:
            result['error'] = repr(error)
            raise
        finally:
            if session is not None:
                result['child_exit_code'] = session.close()
            if (scratch/'stderr.log').exists():
                shutil.copy2(scratch/'stderr.log', args.output/'stderr.log')
            result['elapsed_seconds'] = time.monotonic() - started
            (args.output/'result.json').write_text(json.dumps(result, indent=2)+'\n')
            print(json.dumps(result), flush=True)


if __name__ == '__main__':
    main()

Canyon Run: Continuous Reinforcement Learning

A fixed-wing aircraft learns to fly a winding, 600 m canyon. Rumoca simulates the aircraft and its inner rate controller. A PPO actor and critic learn four simultaneous, continuous commands: forward thrust and body roll, pitch and yaw rates. Start with CartPole for a smaller learning example.

The embedded view keeps the aircraft, canyon, occupancy grid and network display together. Replay trained flight shows the recorded frozen policy. Train from scratch uses the native companion server described below. On a static guide host, the replay works without that server and the training buttons are disabled.

Run the Live Trainer

The current live trainer uses the native CLI on Linux, Python 3.10 or newer, and NumPy. It does not require a new Rumoca compiler build if you already have a native executable. Nix is optional. From the repository root:

python3 -m pip install -r examples/reinforcement_learning/canyon_run/requirements.txt
cargo xtask repo modelica-deps ensure
npm --prefix packages/rumoca-web ci
npm --prefix packages/rumoca-web run build
cargo xtask docs build --book user
python3 examples/reinforcement_learning/canyon_run/serve.py \
  --rumoca rumoca --port 8080

Open http://127.0.0.1:8080/user-guide/simulation/canyon-run.html. If using a locally built executable, pass --rumoca target/release/rumoca. The companion also needs the generated packages/rumoca-web/vendor/three_viewer.js bundle; the repository web build supplies it. --guide-root and --vendor-root can point to an existing built guide and its vendor directory.

Training starts only when requested. One run is allowed at a time. Fast mode uses one CPU core, while Watch training in real time slows the same native loop to its 0.05 s control period. Switch back to fast mode to speed learning up. The default budget is 131,072 transitions with a 900-second wall-time limit, including time spent watching. The native child has a 2 GiB RSS guard and a 4 GiB address-space limit. Results, rollouts, episode records and checkpoints go into a new directory under canyon-runs.

The view uses the same detailed airplane and canyon scene as the earlier example. The occupancy panel shows the actual 8×8 observations. The network panel shows the actor’s actual weights, transposed for the diagram’s destination-row convention, with four continuous outputs. N toggles those actor panel while keeping the occupancy map visible, P toggles paths, T switches training speed and F toggles fullscreen. Buttons expose the same visibility controls. Drag to orbit, right-drag to pan, and scroll or pinch to zoom. Overview camera looks down the course; Top-down camera shows its plan view. Follow camera returns to the aircraft. These camera changes do not affect the simulated flight. The HUD reports a rolling ratio of simulated seconds to wall seconds during live training, including policy-update pauses and episode resets. Replay labels its playback speed separately. The actor panel is reduced in the default view; the episode-return overlay is hidden during replay. Live pose updates are fetched at the control rate and visually interpolated between samples; episode resets snap to the new start. The position readout follows that interpolated pose; other HUD values and occupancy remain sampled simulator data. Decorative rock erosion and animated surfaces are visualization; termination uses the Modelica canyon geometry at the sample times.

Train or Evaluate Without the Viewer

python3 examples/reinforcement_learning/canyon_run/train.py \
  --rumoca rumoca --output canyon-training --steps 131072 --wall 900
python3 examples/reinforcement_learning/canyon_run/evaluate_policy.py \
  --rumoca rumoca --checkpoint canyon-training/checkpoint-131072.npz \
  --output canyon-frozen

Both output directories must be new. Override --source-root if the CMM library is somewhere other than target/cmm/CMM-a642c381. To fly the included trained actor, use --checkpoint examples/reinforcement_learning/canyon_run/actor.npz instead. That actor-only export preserves the trained arrays; it does not contain optimizer or critic state and is not a full training checkpoint.

Each training checkpoint includes the actor, critic, Adam state and random generator state. The native environment is not checkpointed, so resuming the entire simulator state is not supported. Starting another experiment from saved learning state must start a new episode and be recorded as a new run.

Controls and Observations

The aircraft and canyon physics remain those of the earlier guide. The agent now supplies all four commands at once; it does not select bank/climb bins:

Rates use body forward-left-up axes. Positive q is nose-down in those axes; the inner elevator controller converts that sign once. The thrust command scales to throttle using the aircraft mass and maximum propeller thrust. This is a fixed-wing adaptation of continuous thrust/rate control, not a quadrotor’s vertical collective-thrust model.

The actor and critic each read 86 values: normalized position and world velocity, a rotation matrix, body rates, a heading-aligned 8×8 occupancy grid and the previous applied controls. Inputs are clipped to [−5, 5] before the networks. The occupancy grid is computed from true canyon geometry. This is a control-learning example, without camera perception or state estimation.

The policy has two 64-unit tanh layers. It samples four Gaussian latent actions and maps them through tanh to the physical command bounds. The critic has the same hidden widths and one value output. PPO uses the stored latent actions for likelihood ratios; it never substitutes clipped commands or discrete action indices into that ratio. Reported entropy is the latent Gaussian entropy, rather than the entropy after the tanh transform.

The defaults are Adam at 3×10⁻⁴, discount 0.99, GAE parameter 0.95, PPO clip 0.2, 1,024-step rollouts, eight update epochs and 128-row minibatches. A KL check can stop an update early. The core and its finite-difference controls are in ppo.py and test_ppo.py.

Rewards, Termination and Reset

Reward combines forward progress, +0.1 per metre, with lateral-centering, height, command-rate and command-change costs. A crash costs 5 and finishing adds 5. The quadratic height cost coefficient is 0.1; the earlier 0.01 development variant finished but flew too close to the floor.

An episode ends at a sampled wall/floor/rim collision, speed below 4 m/s, or the 600 m finish. It is truncated after 1,200 steps, or 60 simulated seconds. The host captures the old episode’s state, reward and terminal flags before resetting. A true terminal state has zero value bootstrap. A time limit bootstraps from that old final observation, then stops GAE across the reset. The Modelica model does not reset itself underneath the learner.

Measured Development Result

A seed-1 run completed 131,072 transitions and 342 episodes, with 110 finishes. Its last 50 training episodes all finished. The final frozen policy, with learning and action noise disabled, flew 600.60 m in 34.35 simulated seconds. Its minimum sampled altitude was 6.50 m and minimum sampled lateral clearance was 4.95 m. The weights stayed unchanged during that flight.

The same policy also finished six predeclared development checks with ±1 m lateral/height offsets and ±0.05 rad heading offsets. Two initial model-rendering setup failures were retained; those two cases ran successfully after the signed-offset expressions were corrected. No policy tuning occurred between the six cases. These are small perturbations on the training course, not unseen-course, wind, sensor-noise or hardware generalization results.

The included replay records the simulated pose, controls, occupancy observations and actuator commands. demo-flight.json binds it to the source training checkpoint; actor.npz is a separately hashed actor-only export. The replay is labeled as recorded data. Live training comes from the native companion and is labeled separately.

Relationship to the Earlier Example

CanyonRun.Run remains the earlier five-action, fixed-throttle Modelica Q-learner. This chapter’s PPO path uses CanyonRun.ContinuousControl.Run and the Python learner. The old Q-learning results do not measure this new controller. The native PPO demo also does not imply that PPO training has been ported into the WebAssembly learner.

The airborne physics sensitivity examples remain available in CanyonRun.Physics. PPO itself learns from reward and does not use those physics gradients. For the continuous-control motivation, see Kaufmann et al., Nature 2023. The fixed-wing plant, observations, bounds and rewards here are an engineering adaptation, not a reproduction of that quadrotor racing system.

Run the Frozen Actor in Modelica

The trained continuous actor can also run without a Python control loop. CanyonRun.ContinuousControl.CanyonFly samples the same observation, calls the shared NeuralNetwork functions, then holds its four commands for 0.05 s. TrainedCanyonPolicy.mo contains only exported weights. Previous commands are explicit state, so a newly computed action cannot feed back into the observation used to select it.

rumoca sim examples/interactive/canyon_run/CanyonRun.mo \
  --model CanyonRun.ContinuousControl.CanyonFly \
  --source-root examples/reinforcement_learning/NeuralNetwork.mo \
  --source-root examples/reinforcement_learning/canyon_run/TrainedCanyonPolicy.mo \
  --source-root target/cmm/CMM-a642c381 --solver rk-like

A checked native frozen-actor flight reached 600.60 m at 34.35 simulated seconds. Its endpoint differed from the Python-driven flight by less than 1e-11 m. This checks one nominal inference flight, not a new training result or held-out generalization. Continuous PPO training still uses the Python implementation; its optimizer has not been ported to Modelica.

"""Small NumPy PPO core for four continuous, tanh-mapped control channels.

The stored Gaussian latent action defines the score-function likelihood.
The deterministic tanh/physical-unit mapping is part of the environment action;
ratios use the same latent action under old/new policies, never discrete bins.
"""
import math
import numpy as np


class MLP:
    def __init__(self, inputs, outputs, rng, hidden=(64, 64), output_scale=1.0):
        widths = (inputs, *hidden, outputs)
        self.params = {}
        for layer, (left, right) in enumerate(zip(widths, widths[1:])):
            scale = output_scale if layer == len(widths) - 2 else 1.0
            self.params[f'W{layer}'] = rng.normal(size=(left, right)) * scale / math.sqrt(left)
            self.params[f'b{layer}'] = np.zeros(right, dtype=np.float64)

    def forward(self, x, cache=False):
        x = np.asarray(x, dtype=np.float64)
        if x.ndim != 2 or not np.isfinite(x).all():
            raise ValueError('network inputs must be a finite matrix')
        activations = [x]
        layers = len(self.params) // 2
        for layer in range(layers):
            x = x @ self.params[f'W{layer}'] + self.params[f'b{layer}']
            if layer != layers - 1:
                x = np.tanh(x)
            activations.append(x)
        if not np.isfinite(x).all():
            raise FloatingPointError('nonfinite network output')
        return (x, activations) if cache else x

    def backward(self, activations, output_gradient):
        gradient = np.asarray(output_gradient, dtype=np.float64)
        grads = {}
        for layer in reversed(range(len(self.params) // 2)):
            grads[f'W{layer}'] = activations[layer].T @ gradient
            grads[f'b{layer}'] = gradient.sum(axis=0)
            if layer:
                gradient = (gradient @ self.params[f'W{layer}'].T) * (1 - activations[layer] ** 2)
        return grads


def gaussian_log_prob(latent, mean, log_std):
    normalized = (latent - mean) * np.exp(-log_std)
    return -0.5 * np.sum(normalized ** 2 + 2 * log_std + math.log(2 * math.pi), axis=1)


def actor_loss_gradient(actor, log_std, observations, latent, old_log_prob, advantage,
                        clip=0.2, entropy_coefficient=0.001):
    mean, cache = actor.forward(observations, cache=True)
    log_prob = gaussian_log_prob(latent, mean, log_std)
    log_ratio = log_prob - old_log_prob
    bounded_log_ratio = np.clip(log_ratio, -20, 20)
    ratio = np.exp(bounded_log_ratio)
    clipped_ratio = np.clip(ratio, 1 - clip, 1 + clip)
    entropy = np.sum(log_std + 0.5 * (1 + math.log(2 * math.pi)))
    loss = -np.minimum(ratio * advantage, clipped_ratio * advantage).mean() - entropy_coefficient * entropy
    inactive = ((advantage > 0) & (ratio > 1 + clip)) | ((advantage < 0) & (ratio < 1 - clip))
    inactive |= (log_ratio < -20) | (log_ratio > 20)
    d_log_prob = np.where(inactive, 0.0, -advantage * ratio) / len(observations)
    residual = latent - mean
    inv_variance = np.exp(-2 * log_std)
    grads = actor.backward(cache, d_log_prob[:, None] * residual * inv_variance)
    grads['log_std'] = np.sum(d_log_prob[:, None] * (residual ** 2 * inv_variance - 1), axis=0)
    grads['log_std'] -= entropy_coefficient
    metrics = {'actor_loss': float(loss), 'entropy': float(entropy),
               'approx_kl': float(np.mean((ratio - 1) - bounded_log_ratio)),
               'clip_fraction': float(np.mean(np.abs(ratio - 1) > clip))}
    if not math.isfinite(float(loss)) or not all(np.isfinite(g).all() for g in grads.values()):
        raise FloatingPointError('nonfinite PPO actor loss/gradient')
    return float(loss), grads, metrics


def critic_loss_gradient(critic, observations, returns):
    prediction, cache = critic.forward(observations, cache=True)
    error = prediction[:, 0] - returns
    loss = 0.5 * np.mean(error ** 2)
    grads = critic.backward(cache, (error / len(observations))[:, None])
    if not math.isfinite(float(loss)) or not all(np.isfinite(g).all() for g in grads.values()):
        raise FloatingPointError('nonfinite critic loss/gradient')
    return float(loss), grads


def generalized_advantages(reward, value, next_value, terminated, truncated,
                           gamma=0.99, lam=0.95):
    """Bootstrap time limits from the OLD terminal observation; stop GAE at reset."""
    result = np.empty_like(reward, dtype=np.float64)
    running = 0.0
    for index in reversed(range(len(reward))):
        delta = reward[index] + gamma * (not terminated[index]) * next_value[index] - value[index]
        continuation = not (terminated[index] or truncated[index])
        running = delta + gamma * lam * continuation * running
        result[index] = running
    return result, result + value


class Adam:
    def __init__(self, params, learning_rate=3e-4, max_gradient_norm=0.5):
        self.params = params
        self.first = {name: np.zeros_like(value) for name, value in params.items()}
        self.second = {name: np.zeros_like(value) for name, value in params.items()}
        self.steps = 0
        self.learning_rate = learning_rate
        self.max_gradient_norm = max_gradient_norm

    def step(self, grads):
        if grads.keys() != self.params.keys():
            raise ValueError('gradient/parameter roster mismatch')
        norm = math.sqrt(sum(float(np.sum(g * g)) for g in grads.values()))
        if not math.isfinite(norm):
            raise FloatingPointError('nonfinite gradient norm')
        factor = min(1.0, self.max_gradient_norm / max(norm, 1e-12))
        self.steps += 1
        for name, parameter in self.params.items():
            gradient = grads[name] * factor
            self.first[name] *= 0.9
            self.first[name] += 0.1 * gradient
            self.second[name] *= 0.999
            self.second[name] += 0.001 * gradient ** 2
            first = self.first[name] / (1 - 0.9 ** self.steps)
            second = self.second[name] / (1 - 0.999 ** self.steps)
            parameter -= self.learning_rate * first / (np.sqrt(second) + 1e-8)
            if not np.isfinite(parameter).all():
                raise FloatingPointError('nonfinite parameter update')
        return norm


def physical_action(normalized):
    action = np.asarray(normalized, dtype=np.float64)
    if action.shape != (4,) or not np.isfinite(action).all() or np.any(np.abs(action) > 1):
        raise ValueError('action must have four finite normalized channels')
    return np.array([(action[0] + 1) * 4, action[1] * 4, action[2] * 2.5, action[3]], dtype=np.float64)

"""Development PPO training over the released Rumoca CLI, no compiler build."""
import argparse
from datetime import datetime, timezone
import hashlib
import json
import math
import os
from pathlib import Path
import shutil
import statistics
import time
import numpy as np
from environment import ContinuousEnv, MODEL, REPO, configure, sha
from ppo import MLP, Adam, gaussian_log_prob, actor_loss_gradient
from ppo import critic_loss_gradient, generalized_advantages


def preprocess(obs):
    return np.clip(obs, -5, 5)


def episode_summary(rows):
    selected = rows[-50:]
    return {'episodes': len(rows), 'finishes': sum(r['cause'] == 'finish' for r in rows),
            'last50_mean_distance': statistics.mean(r['distance'] for r in selected) if selected else None,
            'last50_mean_return': statistics.mean(r['return'] for r in selected) if selected else None,
            'max_distance': max((r['distance'] for r in rows), default=None),
            'causes': {c: sum(r['cause'] == c for r in rows) for c in ['wall', 'floor', 'rim', 'slow', 'finish', 'timeout']}}


def collect(env, actor, critic, log_std, rng, count, carry, journal, total):
    data = {name: [] for name in ['obs', 'latent', 'logp', 'value', 'reward', 'next_value',
                                  'terminated', 'truncated', 'next_obs', 'applied', 'position']}
    completed = []
    obs, episode_return, episode_length = carry
    for offset in range(count):
        encoded = preprocess(obs)
        mean = actor.forward(encoded[None])[0]
        latent = mean + np.exp(log_std) * rng.normal(size=4)
        normalized = np.tanh(latent)
        logp = gaussian_log_prob(latent[None], mean[None], log_std)[0]
        value = critic.forward(encoded[None])[0, 0]
        next_obs, reward, terminated, truncated, info = env.step(normalized)
        next_value = 0.0 if terminated else critic.forward(preprocess(next_obs)[None])[0, 0]
        values = {'obs': encoded, 'latent': latent, 'logp': logp, 'value': value,
                  'reward': reward, 'next_value': next_value, 'terminated': terminated,
                  'truncated': truncated, 'next_obs': next_obs, 'applied': info['applied'],
                  'position': info['position']}
        for name in data:
            data[name].append(values[name])
        episode_return += reward
        episode_length += 1
        if terminated or truncated:
            _, _, z = info['position']
            cause = ('finish' if info['finished'] else 'timeout' if truncated else
                     'floor' if z <= 0 else 'rim' if z > 30 else 'slow' if info['airspeed'] < 4 else 'wall')
            record = {'global_step': total + offset + 1, 'return': float(episode_return),
                      'length': episode_length, 'distance': info['distance'], 'cause': cause,
                      'terminal_position': info['position'], 'sim_time': info['sim_time']}
            journal.write(json.dumps(record, allow_nan=False) + '\n')
            journal.flush()
            completed.append(record)
            # next_value above was computed BEFORE resetting the native world.
            obs = env.reset()
            episode_return, episode_length = 0.0, 0
        else:
            obs = next_obs
    return {k: np.asarray(v) for k, v in data.items()}, (obs, episode_return, episode_length), completed


def optimize(data, actor, critic, log_std, actor_optimizer, critic_optimizer, rng, epochs, batch_size):
    advantage, returns = generalized_advantages(data['reward'], data['value'], data['next_value'],
                                               data['terminated'], data['truncated'])
    advantage = (advantage - advantage.mean()) / max(float(advantage.std()), 1e-8)
    metrics = {}
    updates = 0
    for epoch in range(epochs):
        order = rng.permutation(len(advantage))
        for start in range(0, len(order), batch_size):
            indices = order[start:start + batch_size]
            _, grads, metrics = actor_loss_gradient(actor, log_std, data['obs'][indices],
                data['latent'][indices], data['logp'][indices], advantage[indices])
            if metrics['approx_kl'] > .03:
                return {**metrics, 'optimizer_steps': updates, 'early_stop_kl': True}
            actor_optimizer.step(grads)
            np.clip(log_std, -3, -.3, out=log_std)
            value_loss, value_grads = critic_loss_gradient(critic, data['obs'][indices], returns[indices])
            critic_optimizer.step(value_grads)
            metrics['critic_loss'] = value_loss
            updates += 1
    return {**metrics, 'optimizer_steps': updates, 'early_stop_kl': False}


def checkpoint(path, actor, critic, log_std, actor_optimizer, critic_optimizer, rng, step):
    arrays = {'log_std': log_std.copy()}
    for prefix, network, optimizer in [('actor', actor, actor_optimizer), ('critic', critic, critic_optimizer)]:
        arrays.update({prefix + '_' + name: value for name, value in network.params.items()})
        arrays.update({prefix + '_adam_m_' + name: value for name, value in optimizer.first.items()})
        arrays.update({prefix + '_adam_v_' + name: value for name, value in optimizer.second.items()})
        arrays[prefix + '_adam_steps'] = np.array(optimizer.steps)
    np.savez_compressed(path, **arrays)
    path.with_suffix('.json').write_text(json.dumps({'global_step': step, 'rng_state': rng.bit_generator.state,
        'checkpoint_sha256': sha(path), 'native_environment_resume': 'NOT_SUPPORTED; any continuation starts a fresh episode and is a new experiment'}, indent=2) + '\n')


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--output', type=Path, required=True)
    parser.add_argument('--steps', type=int, default=32768)
    parser.add_argument('--rollout', type=int, default=1024)
    parser.add_argument('--epochs', type=int, default=8)
    parser.add_argument('--batch', type=int, default=128)
    parser.add_argument('--wall', type=int, default=1800)
    parser.add_argument('--cpu', type=int, default=min(os.sched_getaffinity(0)))
    parser.add_argument('--rumoca', default='rumoca')
    parser.add_argument('--source-root', type=Path, default=REPO / 'target/cmm/CMM-a642c381')
    args = parser.parse_args()
    if min(args.steps, args.rollout, args.epochs, args.batch, args.wall) < 1:
        raise ValueError('budgets must be positive')
    configure(args.rumoca, args.source_root)
    os.sched_setaffinity(0, {args.cpu})
    args.output.mkdir(exist_ok=False)
    for name in ['train.py', 'ppo.py', 'environment.py', 'stdio.py', 'test_ppo.py']:
        shutil.copy2(Path(__file__).with_name(name), args.output / name)
    rng = np.random.default_rng(1)
    actor = MLP(86, 4, rng, output_scale=.01)
    critic = MLP(86, 1, rng)
    actor.params['b2'][0] = np.arctanh(-.5)  # Initial mean = 2m/s² forward thrust.
    log_std = np.full(4, math.log(.15))
    actor_optimizer = Adam({**actor.params, 'log_std': log_std})
    critic_optimizer = Adam(critic.params)
    result = {'completed': False, 'method': 'continuous Gaussian latent PPO, tanh physical controls',
              'seed': 1, 'numpy': np.__version__, 'python': os.sys.version,
              'model_sha256': sha(MODEL),
              'steps_requested': args.steps, 'rollout': args.rollout, 'epochs': args.epochs,
              'batch_size': args.batch, 'wall_seconds': args.wall, 'cpu': args.cpu,
              'gamma': .99, 'gae_lambda': .95, 'clip': .2, 'learning_rate': 3e-4,
              'actor_hidden': [64, 64], 'critic_hidden': [64, 64], 'sample_dt': .05,
              'purpose': 'ordinary continuous-controller development; not held-out or paper reproduction'}
    (args.output / 'START.json').write_text(json.dumps(result, indent=2) + '\n')
    env = None
    started = time.monotonic()
    total, rollout_number, episodes = 0, 0, []
    try:
        env = ContinuousEnv(args.output, args.wall, args.cpu)
        (args.output / 'HELLO.json').write_text(json.dumps(env.hello, indent=2) + '\n')
        carry = (env.reset(), 0.0, 0)
        checkpoint(args.output / 'checkpoint-000000.npz', actor, critic, log_std,
                   actor_optimizer, critic_optimizer, rng, 0)
        with (args.output / 'episodes.jsonl').open('w') as journal:
            while total < args.steps:
                count = min(args.rollout, args.steps - total)
                data, carry, completed = collect(env, actor, critic, log_std, rng, count, carry, journal, total)
                total += count
                rollout_number += 1
                episodes.extend(completed)
                np.savez_compressed(args.output / f'rollout-{rollout_number:04}.npz', **data)
                metrics = optimize(data, actor, critic, log_std, actor_optimizer, critic_optimizer,
                                   rng, args.epochs, args.batch)
                checkpoint(args.output / f'checkpoint-{total:06}.npz', actor, critic, log_std,
                           actor_optimizer, critic_optimizer, rng, total)
                progress = {**episode_summary(episodes), 'global_steps': total,
                            'elapsed_wall': time.monotonic() - started, 'log_std': log_std.tolist(), **metrics}
                (args.output / 'PROGRESS.json').write_text(json.dumps(progress, indent=2) + '\n')
                print(json.dumps(progress), flush=True)
        result.update(completed=True, global_steps=total, **episode_summary(episodes))
    except Exception as error:
        result.update(error=repr(error), global_steps=total, **episode_summary(episodes))
        raise
    finally:
        if env is not None:
            result['native_close'] = env.close()
        result['elapsed_wall'] = time.monotonic() - started
        result['ended_utc'] = datetime.now(timezone.utc).isoformat()
        (args.output / 'RESULT.json').write_text(json.dumps(result, indent=2) + '\n')
        files = [p for p in sorted(args.output.iterdir()) if p.is_file()]
        (args.output / 'SHA256SUMS').write_text(''.join(sha(p) + '  ' + p.name + '\n' for p in files))
        print(json.dumps(result), flush=True)


if __name__ == '__main__':
    main()

Targets and Templates

Rumoca can render a compiled model into symbolic frameworks, compiled simulation kernels, eFMI artifacts, or Modelica source. Code generation is target-directory based: a target is a target.toml manifest plus Jinja templates, and each target declares which compiler IR stage it consumes.

Listing Targets

rumoca targets

Built-in targets include:

Targets without a complete checked artifact and executable or independent validation evidence are intentionally absent. Rumoca does not expose aliases for removed target names or route those names through a weaker IR.

The rumoca targets table also reports a readiness level (0 = experimental … 2 = validated) and per-feature support columns (scalarization, tensor features such as matmul/linear solve/elementwise/stencil kernels, events, AD, …) for each target. Treat the table — not this page — as the current source of truth.

Rendering a Target

rumoca compile examples/models/SympyDecay.mo \
  --model SympyDecay \
  --target c-ode \
  --output /tmp/decay_c_ode

--output may be a file or directory depending on what the target renders.

Codegen Scenarios

Like simulations, generation jobs worth repeating belong in a rumoca-scenario.toml with task = "codegen". Runnable examples live under examples/codegen/ and write into examples/codegen/gen/ (git-ignored):

• examples/codegen/rumoca-scenario.sympy_decay_jax_ode.toml — checked ODE RHS JAX target
• examples/codegen/rumoca-scenario.sympy_decay_c_ode.toml — checked ODE RHS C target
• examples/codegen/rumoca-scenario.sympy_decay_checked_dae_report.toml — custom checked-DAE report target
• examples/codegen/rumoca-scenario.sympy_decay_custom_checked_variables.toml — raw checked-DAE Jinja template

IR Dumps vs Targets

If what you want is to see a compiler stage rather than generate project code, use --emit instead of a target — see Inspecting and Debugging Models.

Custom Targets

When the built-in targets do not fit, write your own. There are two levels:

Raw Jinja Templates

For one-off generation, pass a .jinja file directly. --phase chooses which IR the template receives (default dae):

rumoca compile Model.mo --target my_template.jinja --phase flat -o out.txt

The template gets the canonical stage projection as its context. The repository example examples/codegen/custom_checked_variables.jinja shows this workflow.

To learn the available fields, dump the matching IR as JSON first:

rumoca compile Model.mo --emit flat-json | head -50

Target Directories (target.toml)

For anything reusable, create a directory containing a target.toml manifest and the templates it references, then pass the directory:

rumoca compile Model.mo --target path/to/my_target -o out/

The manifest declares which IR stage the target consumes and which templates render which output files. The target — not individual templates — owns the IR choice, so a bundle stays consistent.

The repository ships a complete worked example: examples/codegen/checked_dae_report/target.toml renders a readable report from the canonical checked DAE projection.

Design Rule: Language Knowledge Lives in Targets

Rumoca’s compiler phases are deliberately target-agnostic: no Rust code special-cases C, CUDA, Python, or MLIR. Everything language-specific belongs in target.toml metadata and templates. If a custom target needs information the IR does not expose, that is a compiler feature request — not something to hack around in a template.

eFMI Algorithm Code Export (GALEC)

Rumoca can project a compiled model into eFMI Algorithm Code — the GALEC (Guarded Algorithmic Language for Embedded Control) .alg representation — and package it as a schema-valid eFMU container.

The three targets

All three GALEC targets consume the dae IR and accept fixed-sample discrete models only — models with no continuous states and no der().

The eFMI container? column describes the CLI packaging step. The GUI’s Generate Code (below) renders the inspectable .alg/.h/.c sources for any target but does not itself build the container.

Exporting from the CLI

rumoca compile Model.mo --target galec -o out/
rumoca compile Model.mo --target galec-production -o out/
rumoca compile Model.mo --target embedded-c-galec -o out/

The galec and galec-production targets write the eFMU container in two forms — a directory and the equivalent .efmu zip:

out/
  Model/                eFMU container, directory form
    __content.xml
    schemas/
    AlgorithmCode/      Model.alg + manifest.xml
    ProductionCode/     galec-production only: Model.c/.h + manifest.xml
  Model.efmu            eFMU container, zip form (same content)

The embedded-c-galec target instead writes plain out/Model.h and out/Model.c — no manifest and no container.

Code generation in the GUI

The scenario config editor (the VS Code custom editor and its web build) exposes a Code Generation task. Pick a GALEC target from the codegen target dropdown (galec, galec-production, or embedded-c-galec), choose an output location, and select Generate Code.

The GUI renders the inspectable GALEC sources in stages: first the .alg Algorithm Code text, then generated C header/source from the current .alg editor contents. This browser render is identity-free — it does not mint the eFMU container (the __content.xml registry, per-representation manifest.xml, and checksum web). For galec/galec-production, those container artifacts are produced by the CLI packaging step (rumoca compile … --target galec-production). Use the GUI to inspect and edit the target-language code, and use the CLI when you need the full eFMU package.

Try GALEC Production Code in the guide

The example below is a fixed-sample discrete counter, which is the subset the current GALEC projection accepts. Select Generate .alg to project the Modelica source into GALEC Algorithm Code. The .alg artifact opens in the same Monaco editor surface as the Modelica input, with GALEC syntax highlighting and the GALEC language service diagnostics/hover/definition hooks active. Then select Generate C/H in the .alg panel to generate C header and source from the current GALEC editor text.

// rumoca-live-scenario: ../repo-examples/codegen/rumoca-scenario.galec_counter_production.toml

Native run:

cargo run -p rumoca -- \
  compile examples/models/GalecCounter.mo \
  --model GalecCounter \
  --target galec-production \
  --output examples/codegen/gen/galec_counter_production

See also

• Targets and Templates — the full target list and the live rumoca targets readiness table.
• The authoritative contract for what each target emits and which conformance rung it claims is SPEC_0034; this page does not restate its rules.

Troubleshooting and FAQ

Compilation

“class not found” / unresolved names — The model references a package Rumoca cannot see. Add the library with --source-root (CLI), top-level source_roots (scenario), or rumoca-workspace.toml (workspace). Check MODELICAPATH if you rely on it. For the Modelica Standard Library, name the directory that contains Modelica/, ModelicaServices/, and Complex.mo, not the Modelica/ package itself: Modelica references the other two top-level packages. See Using Modelica Libraries. The parser’s view of a file, rumoca compile Model.mo --emit ast-json, does not resolve names and so stays available while you sort out the library paths.

“Duplicate class ‘X’ … with non-identical definition” (EM001) — Two files in the same source root define the same class name. This commonly happens when an old copy of a model sits next to a new one in the same directory; direct rumoca sim file.mo runs include the file’s directory as a source root.

“structurally singular system: N matched out of M equations” — The equation system cannot be matched one-to-one with its unknowns. The diagnostic names the unmatched equations and unknowns; run rumoca compile --inspect structure for the full matching. Common causes:

• Genuinely unbalanced models (forgotten equation, extra variable).
• High-index DAE formulations that current index reduction does not yet handle, such as a Cartesian pendulum with an explicit constraint — see Language Support Status. Reformulating with generalized coordinates usually fixes it.

Model is balanced but produces wrong dynamics — Dump the system the solver actually integrates with rumoca compile --emit dae-mo and compare it to your intent. Please file an issue with a minimal model if the lowering looks wrong.

Simulation

“Step size is too small at time = …” — The implicit solver stalled, most often near a dense cascade of state events (rapid relay switching). Use the explicit solver: --solver rk-like or annotation(experiment(Solver = "rk-like")).

NaN/Inf failures — rumoca sim automatically re-runs with NaN tracing and names the offending variables. To investigate further, evaluate the model at a chosen point: rumoca sim Model.mo --inspect eval --at "x=...@t". Typical causes: division by a variable that crosses zero, sqrt/log of a negative value, or missing initial values.

Simulation stops early with a message — The model called terminate(...); the message is recorded in the report. assert failures likewise carry their message and source location.

My batch run ignored experiment(StopTime=...): native direct CLI runs start at the model’s experiment(StartTime) (else 0) and take the end time from --t-end, else the model’s experiment(StopTime), else one time unit after the start; rumoca sim bench does the same. Batch scenario runs use [sim] t_end. Browser and scheduled interactive runs are user-terminated and intentionally do not stop at either value.

Results

The HTML report did not open — It is written next to where you ran the command (<MODEL>_results.html by default, -o to override). Open it in any browser.

Too many variables in the report — Add [[plot.views]] sections to the scenario to define focused plots.

VS Code

No diagnostics / completion — Check that the extension is active for the file (it activates on .mo and rumoca-scenario.toml). If you enabled rumoca.useSystemServer, ensure rumoca-lsp is on PATH; otherwise the bundled server is used.

Library completion missing — Add shared library roots to rumoca-workspace.toml. For the repository’s own examples, run cargo xtask repo modelica-deps ensure first.

Runnable Blocks in This Book

The ▶ Simulate button reports the WASM package is missing — Live examples need the WASM package deployed next to the book. They work on the published site; for a local build, use cargo xtask docs serve; it builds the missing local WASM package for live examples before serving the books.

The editor has no syntax highlighting — Monaco loads from a CDN; when offline, the examples fall back to plain text editors but still simulate.

Performance

Compilation feels slow on repeat runs — Check the cache: rumoca cache status. Direct file runs are cached by content; use rumoca sim bench to separate compile time from simulation throughput.

Browser runs are slower than native — Expected, especially for large package trees; use the native CLI for MSL-heavy work.

Reporting Bugs

File issues at https://github.com/CogniPilot/rumoca/issues with a minimal model. The Web Playground is a convenient way to confirm a reproduction without local setup.