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 sympy \
  --output /tmp/sympy_decay

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 and the web playground honor it; native CLI runs currently use --t-end (default 1.0 s) or the [sim] section of a scenario file instead.

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:

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.02 "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 are rotated into 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. The model also exposes input Real aoa_cmd; an aoa_motor state follows that command with der(aoa_motor) = (aoa_cmd - aoa_motor) / aoa_tau. A structural interactive flag selects whether the immersed-boundary mask uses the pre-simulation parameter aoa or the lagged state aoa_motor. 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-solve 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). GPU v1 runs in f32 with events and algebraics frozen at their settled initial values, which is exact for the normal batch run because interactive = false makes the mask depend only on pre-simulation parameters. 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 (~6,500 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
    "Use live AoA motor state for the airfoil mask" annotation(Evaluate = true);
  input Real aoa_cmd(start = 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(start = mc0) "Commanded NACA max camber";
  input Real pc(start = pc0) "Commanded NACA camber position";
  input Real tk(start = 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";
  Real sig[NX, NY] "Solid mask (1 inside the airfoil)";
  // States start at rest (default start = 0): an impulsive wind-tunnel
  // start where the freestream sweeps in through the boundary relaxation.
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;
  for i in 1:NX loop
    for j in 1:NY loop
      if interactive then
        sc[i, j] = ((i - 0.5) * dx - xle) * cos(aoa_motor * pi / 180.0)
          - ((j - 0.5) * dy - Ly / 2.0) * sin(aoa_motor * pi / 180.0);
        nc[i, j] = ((i - 0.5) * dx - xle) * sin(aoa_motor * pi / 180.0)
          + ((j - 0.5) * dy - Ly / 2.0) * cos(aoa_motor * pi / 180.0);
        sig[i, j] =
          0.5 * (1.0 - tanh((abs(nc[i, j]
              - (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)));
      else
        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);
        sig[i, j] =
          0.5 * (1.0 - tanh((abs(nc[i, j]
              - (if sc[i, j] < pc0 then mc0 / pc0 ^ 2 * (2.0 * pc0 * sc[i, j] - sc[i, j] ^ 2)
                 else mc0 / (1.0 - pc0) ^ 2
                   * ((1.0 - 2.0 * pc0) + 2.0 * pc0 * sc[i, j] - sc[i, j] ^ 2)))
            - sqrt((5.0 * tk0 * (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)));
      end if;
    end for;
  end for;
  // Interior: momentum + artificial-compressibility continuity.
  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)
        - sig[i, j] * 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)
        - sig[i, j] * 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:

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:

#![allow(unused)]
fn main() {
let mut model = rumoca_opt::DifferentiableModel::from_dae_default(
    &compiled.dae,
    &rumoca_solver::SimOptions::default(),
)?;
let trainables = rumoca_opt::TrainableSet::by_names(&model, &["W[1,1]", "b[1]"])?;
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,
)?;
}

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 checks for initial parameters, observables, gradients, final parameters, and JAX manual gradients versus jax.value_and_grad.

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 native Solve IR CPU simulator executes the Modelica backpropagation equations, while the JAX side uses jax.value_and_grad for the same RHS.

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 scenario file ([sim] t_end). Browser/live runs do honor it.

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 sympy -o out/        # built-in 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[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

import rumoca as rm

session = rm.Session(roots=["libs/CMM"])
m = session.load("Quadrotor.mo", model="Quadrotor")
m                                  # repr shows a summary
m.parameters["mass"].value         # typed, autocompletes
m.parameters["mass"].kind          # "tunable" or "structural"
m.states["body.v[1]"]              # subscripted names resolve

r = m.simulate(t=(0, 10), dt=0.01)
r.plot("body.v[1]")                # 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):

session, model, config = rm.Session.from_scenario("rumoca-scenario.toml")
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("rumoca-scenario.toml")
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(roots=["libs/CMM"])
run = session.run_scenario(
    "rumoca-scenario.toml",
    overrides={"t_end": 2.0, "mode": "lockstep", "output": "results.csv"},
)

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 drag in [0.0, 0.05, 0.1]:
    r = m.with_params(rotor_drag=drag).simulate(t=(0, 10))
    ...

m.simulate(params={...}, start={...}) applies overrides for a single run. Overriding a structural parameter (one that affects sizing/instantiation) or a parameter another parameter’s value depends on raises StructuralParamError — those need a recompile, not a runtime override, so a sweep is never silently wrong.

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

m6 = m.with_params(n_rotors=6, recompile=True)   # re-instantiate with n_rotors=6
m6.simulate(t=(0, 10))

Live symbolic export

Turn a model into a live object in the symbolic framework of your choice — no file dance:

cm = m.to_casadi()                 # CasadiModel (form="dae")
cm.ode                             # explicit RHS (CasADi's native DAE form)
cm.dae                             # {x, z, p, t, ode, alg} — feed ca.integrator
S  = cm.jacobian("ode", "p")       # exact AD parameter sensitivity

sm = m.to_sympy()                  # SymPy model; sm.solve_explicit()
jm = m.to_jax()                    # JAX ode_fn + diffrax simulate

Pass form="solve" for the scalarized, causalized explicit form (the same source the C/FMI/Rust backends use) — you get a SolveExport exposing the explicit right-hand side xdot = rhs(x, u, p) plus state_names/ parameter_names:

se = m.to_casadi(form="solve")     # SolveExport; se.rhs is a ca.Function
se = m.to_jax(form="solve")        # SolveExport; se.rhs is jit/grad/vmap-able

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 m.codegen("casadi-mx").save_all("out/").

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

session = rm.Session(roots=["libs/Modelica"])
written = session.codegen_file(
    "model.mo",
    "MyModel",
    "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:

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(roots=[r'${CMAKE_CURRENT_SOURCE_DIR}/libs/Modelica']); \
session.codegen_file(r'${CMAKE_CURRENT_SOURCE_DIR}/model.mo', 'MyModel', \
'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}/model.mo"
  VERBATIM
)

Scenario-driven build or test steps use the same pattern:

add_custom_target(run_rumoca_scenario
  COMMAND "${Python3_EXECUTABLE}" -c
          "import rumoca as rm; rm.Session().run_scenario(r'${CMAKE_CURRENT_SOURCE_DIR}/rumoca-scenario.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 CasadiModel. --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.

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.

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.

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).

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 <id>   # repeatable
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 1.0), --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

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       # simulation stop time
atol = 1e-6        # optional absolute solver tolerance
rtol = 1e-6        # optional relative solver tolerance
solver = "auto"    # auto | bdf | esdirk34 | trbdf2 | 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.

[[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"]

[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.

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

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 currently take end time from --t-end (default 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.

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
t_end = 10.0
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.

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, including after the run reports Finished, to restore the vehicle and restart at t = 0. 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

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.

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 sympy -o out -v

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

Targets and Templates

Rumoca can render a compiled model into other languages and ecosystems: symbolic math packages, compiled simulation kernels, FMUs, or Modelica source at any pipeline stage. 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:

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 sympy \
  --output /tmp/sympy_decay

--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.ball_jax.toml — built-in JAX target
• examples/codegen/rumoca-scenario.sympy_decay_sympy.toml — built-in SymPy target
• examples/codegen/rumoca-scenario.sympy_decay_standalone_web.toml — custom web target
• examples/codegen/rumoca-scenario.sympy_decay_custom_casadi.toml — raw 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 serialized IR as its context. The repository example examples/codegen/custom_casadi.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/standalone_web/target.toml renders a standalone HTML page plus companion JavaScript from one model.

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. See Using Modelica Libraries.

“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")). For genuinely stiff smooth systems, try esdirk34 or trbdf2 instead.

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 run used t_end = 1.0 even though the model has experiment(StopTime=...) — Native direct CLI runs take the end time from --t-end (default 1.0); scenario runs use [sim] t_end. The browser examples and playground do honor the annotation.

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.