mj_step is short. In MuJoCo 3.14.0 it checks positions and velocities, calls mj_forward, checks accelerations, and hands the result to the integrator (src/engine/engine_forward.c, lines 2012 to 2058). Everything that computes physics lives in mj_forward, which runs three stages in a fixed order: position, velocity, acceleration. Sensors and energy follow each stage, and the control callback runs before actuation. Each stage reads what the earlier ones wrote and writes a known set of mjData fields. This lesson walks the stages in the source, then measures which fields each one writes by filling mjData with NaN and running the stages one at a time.
[!established] A map of the source (3.14.0, tag commit 9ecbb9d) All paths are under
src/engine/.
Where What engine_forward.cmj_step(2012),mj_forwardSkip(1947), the stagesmj_fwdPosition(133),mj_fwdVelocity(183),mj_fwdActuation(365),mj_fwdAcceleration(1006),mj_fwdConstraint(1314), the integratorsmj_Euler(1597),mj_RungeKutta(1624),mj_implicit(1915), andmj_step1/mj_step2(2061, 2094)engine_core_smooth.cmj_kinematics(240),mj_comPos(250),mj_tendon(931),mj_transmission(1269),mj_makeM(1970),mj_factorM(1979),mj_comVel(2206),mj_rne(2355),mj_rnePostConstraint(2488)engine_collision_driver.cmj_collision(606), with broad and narrow phaseengine_core_constraint.cmj_makeConstraint(2919),mj_projectConstraint(3275),mj_referenceConstraint(3443)engine_passive.cmj_passive(1006): springs, dampers, fluid, gravity compensation, adhesion intoqfrc_passiveengine_solver.cmj_solPGS(748),mj_solCG(3024),mj_solNewton(3036)engine_sensor.cmj_sensorPos(1465),mj_sensorVel(1516),mj_sensorAcc(1570)Line numbers are where each function’s definition starts.
The lab simulates pick_place for 0.5 s, so the cube rests on the table and the gripper holds its pose: 12 contacts and 49 constraint rows. Fill mjData with NaN overwrites every derived array. Then press the stages in order. Each button runs one public stage function and lists the fields it filled; after the last one, the lab compares qacc with what mj_forward computed from the same state.
```lab pipeline {“dock”: true, “title”: “mj_forward, one stage at a time”}
## Mathematics
### The stages and what they write
`mj_fwdPosition` runs kinematics, then builds and factors the mass matrix, then detects collisions, then builds the constraint rows. `mj_fwdVelocity` adds everything that depends on velocity: passive forces, the bias forces $c(q,\dot q)$ from recursive Newton-Euler without acceleration, and the constraint reference accelerations. `mj_fwdActuation` turns controls into forces. `mj_fwdAcceleration` forms the unconstrained acceleration,
$$\ddot q_\text{smooth} = M^{-1}\big(\tau_\text{passive} - c + \tau_\text{applied} + \tau_\text{actuator}\big),$$
and `mj_fwdConstraint` warmstarts and runs the solver to get the constrained $\ddot q$. Measured on `pick_place`, with every non-input float array set to NaN first:
| Stage | Fields it fills |
|---|---|
| `mj_fwdPosition` | `xpos`, `xquat`, `xmat`, `xipos`, `ximat`, `xanchor`, `xaxis`, `geom_xpos`, `geom_xmat`, `site_xpos`, `site_xmat`, `cam_*`, `light_*`, `subtree_com`, `cdof`, `cinert`, `crb`, `M`, `qLD`, `qLDiagInv`, `ten_length`, `ten_J`, `actuator_length`, `actuator_moment`, `efc_J`, `efc_pos`, `efc_margin`, `efc_frictionloss`, `efc_diagA`, `efc_R`, `efc_D`, `efc_KBIP` |
| `mj_fwdVelocity` | `cvel`, `cdof_dot`, `ten_velocity`, `actuator_velocity`, `qfrc_spring`, `qfrc_damper`, `qfrc_gravcomp`, `qfrc_fluid`, `qfrc_adhesion`, `qfrc_passive`, `qfrc_bias`, `efc_vel`, `efc_aref` |
| `mj_fwdActuation` | `actuator_force`, `qfrc_actuator` |
| `mj_fwdAcceleration` | `qfrc_smooth`, `qacc_smooth` |
| `mj_fwdConstraint` | `efc_b`, `efc_force`, `qfrc_constraint`, `qacc`, and the per-island copies `iacc`, `iefc_force`, ... |
| `mj_sensorPos`, `mj_sensorVel`, `mj_sensorAcc` | `sensordata`, each stage its own sensors: here 8 of its 9 values after the position stage and the last after the acceleration stage |
Run in this order, the stages reproduce `mj_forward`'s `qacc` exactly (difference 0.0). Some fields were left as NaN by every stage, and each has a reason. `cacc`, `cfrc_int`, `cfrc_ext`, `subtree_linvel` and `subtree_angmom` are computed lazily: `mj_forward` clears the flags `flg_rnepost` and `flg_subtreevel`, and a sensor that needs those quantities computes them on demand (`engine_sensor.c`, around lines 843, 962 and 1600). `qDeriv`, `qLU`, `qH` and `qHDiagInv` belong to the integrator: `mj_implicit` fills them after `mj_forward`. `qfrc_inverse` belongs to `mj_inverse`, and `energy` stays empty unless the energy flag is on.
### After the step, mjData describes the state before it
`mj_step` integrates after `mj_forward` and does not run the position stage again. So when it returns, `qpos` and `qvel` are new, but every derived quantity, from `xpos` to `sensordata`, still describes the state the step started from. Measured with the arm's joints moving at 1 rad/s: after one step, the stored body positions differ by 0.000 mm from those before the step and by 3.631 mm from those of the new `qpos`. The sensors likewise match the old state exactly and differ from the new one by up to 0.0242. This is why Lesson 12.1 calls `mj_forward` at the end of an environment's `step`, and why `mujoco.rollout` pairs each step's sensors with the state before it (Lesson 18.2). `mj_step1` and `mj_step2` split the same work at the control callback, so that a controller can read fresh sensors and set `ctrl` in between.
### Skipping stages
`mj_forwardSkip(m, d, skipstage, skipsensor)` reuses what is already in `mjData` for the stages below `skipstage`. After a change of `qvel` only, everything the position stage computes is still valid, so `mjSTAGE_POS` skips it. Over 2000 random velocity perturbations on `pick_place`, the result equals `mj_forward`'s to the last bit (difference 0.0), at 26.3 µs against 38.2 µs. Finite-difference derivatives (Lesson 20.4) use exactly this trick: perturbations of velocity and control skip the position stage.
> [!recommendation] Reading the source
> Start at `mj_step` and follow one quantity: find every function that writes it (`grep` the field name in `src/engine`), then confirm the order by instrumenting a run, as part (1) does with NaN. Read the version you run: `mujoco.__version__` names the tag, and the line numbers here are those of 3.14.0.
## Implementation
```io
INPUT: pick_place.xml after 0.5 s of simulation, MuJoCo 3.14.0
PROCESS: fill derived mjData fields with NaN and run the stage functions one at a time; compare with mj_forward; the state stored after mj_step; mj_forwardSkip against mj_forward
OUTPUT: printed tables
```python file=examples/l20_1_pipeline.py “"”Lesson 20.1: inside mj_step, stage by stage (Level 20).
INPUT pick_place.xml (arm, gripper with a tendon and an equality, cube, sensors) simulated for 0.5 s from its keyframe so that contacts are active; MuJoCo 3.14.0, whose mj_step is in src/engine/engine_forward.c PROCESS (1) every float array of mjData that is not an input is filled with NaN; the public stage functions then run one at a time in the order mj_forward calls them, and each stage’s newly written fields are listed; the resulting qacc is compared with mj_forward’s; (2) the fields no stage wrote, and why; (3) after mj_step: how far the positions and sensors stored in mjData are from the new state; (4) mj_forwardSkip: after a change of velocity only, the position stage can be skipped; equality with mj_forward and the time saved OUTPUT printed tables
Run: python examples/l20_1_pipeline.py “””
import time
import mujoco import numpy as np
from mjcourse import model_path
INPUTS = {“qpos”, “qvel”, “act”, “ctrl”, “qfrc_applied”, “xfrc_applied”, “mocap_pos”, “mocap_quat”, “userdata”, “qacc_warmstart”, “plugin_state”, “history”, “eq_active”} STAGES = [(“mj_fwdPosition”, mujoco.mj_fwdPosition), (“mj_sensorPos”, mujoco.mj_sensorPos), (“mj_fwdVelocity”, mujoco.mj_fwdVelocity), (“mj_sensorVel”, mujoco.mj_sensorVel), (“mj_fwdActuation”, mujoco.mj_fwdActuation), (“mj_fwdAcceleration”, mujoco.mj_fwdAcceleration), (“mj_fwdConstraint”, mujoco.mj_fwdConstraint), (“mj_sensorAcc”, mujoco.mj_sensorAcc)]
def float_fields(data: mujoco.MjData) -> dict[str, np.ndarray]: “"”Every non-empty float64 array of mjData that is not an input (fresh views: the arena moves).””” out = {} for name in dir(data): if name.startswith(“_”) or name in INPUTS: continue try: value = getattr(data, name) except Exception: continue if isinstance(value, np.ndarray) and value.dtype == np.float64 and value.size: out[name] = value return out
def scene() -> tuple[mujoco.MjModel, mujoco.MjData]: model = mujoco.MjModel.from_xml_path(str(model_path(“pick_place”))) data = mujoco.MjData(model) mujoco.mj_resetDataKeyframe(model, data, 0) while data.time < 0.5: mujoco.mj_step(model, data) return model, data
if name == “main”: model, data = scene() mujoco.mj_forward(model, data) reference_qacc = data.qacc.copy() print(f”MuJoCo {mujoco.version}; pick_place after 0.5 s: nv {model.nv}, ncon {data.ncon}, nefc {data.nefc}, “ f”solver {mujoco.mjtSolver(model.opt.solver).name}, integrator {mujoco.mjtIntegrator(model.opt.integrator).name}”) print(“(1) fields each stage writes (filled with NaN beforehand)”) for name, value in float_fields(data).items(): value[…] = np.nan written = set() for stage, fn in STAGES: fn(model, data) now = {name for name, value in float_fields(data).items() if np.isfinite(value).all()} new = sorted(now - written) written |= now partial = [f”{n} ({int(np.isfinite(v).sum())} of {v.size})” for n, v in float_fields(data).items() if n not in now and np.isfinite(v).any()] print(f” {stage:<19} {‘, ‘.join(new) if new else ‘(nothing complete)’}” + (f”; partly: {‘, ‘.join(partial)}” if partial else “”)) print(f” qacc from the stages against mj_forward: max difference {np.abs(data.qacc - reference_qacc).max():.1e}”) print(“(2) fields no stage wrote”) untouched = sorted(name for name, value in float_fields(data).items() if not np.isfinite(value).all()) print(“ “ + “, “.join(untouched)) print(“(3) after one mj_step, with the arm’s joints moving at 1 rad/s”) model, data = scene() data.qvel[:7] = 1.0 mujoco.mj_forward(model, data) before_xpos, before_sensor = data.xpos.copy(), data.sensordata.copy() mujoco.mj_step(model, data) fresh = mujoco.MjData(model) fresh.qpos[:], fresh.qvel[:], fresh.act[:], fresh.ctrl[:] = data.qpos, data.qvel, data.act, data.ctrl mujoco.mj_forward(model, fresh) print(f” body positions in mjData against those of the new qpos: max difference “ f”{1000 * np.abs(data.xpos - fresh.xpos).max():.3f} mm; against those before the step: “ f”{1000 * np.abs(data.xpos - before_xpos).max():.3f} mm”) print(f” sensordata in mjData against the new state’s: max difference {np.abs(data.sensordata - fresh.sensordata).max():.2e}; “ f”against the state before the step: {np.abs(data.sensordata - before_sensor).max():.2e}”) print(“(4) mj_forwardSkip after a change of qvel only”) model, data = scene() mujoco.mj_forward(model, data) work = mujoco.MjData(model) work.qpos[:], work.qvel[:], work.act[:], work.ctrl[:] = data.qpos, data.qvel, data.act, data.ctrl mujoco.mj_forward(model, work) rng = np.random.default_rng(0) errors, t_full, t_skip, n = [], 0.0, 0.0, 2000 for _ in range(n): dv = rng.normal(0.0, 0.01, model.nv) work.qvel[:] = data.qvel + dv t0 = time.perf_counter() mujoco.mj_forwardSkip(model, work, mujoco.mjtStage.mjSTAGE_POS, 0) t_skip += time.perf_counter() - t0 skipped = work.qacc.copy() t0 = time.perf_counter() mujoco.mj_forward(model, work) t_full += time.perf_counter() - t0 errors.append(np.abs(work.qacc - skipped).max()) print(f” {n} random velocity perturbations: max |qacc(skip) - qacc(full)| {max(errors):.1e}; “ f”mj_forward {1e6 * t_full / n:.1f} us, mj_forwardSkip(mjSTAGE_POS) {1e6 * t_skip / n:.1f} us”)
Output:
```text
MuJoCo 3.14.0; pick_place after 0.5 s: nv 27, ncon 12, nefc 49, solver mjSOL_NEWTON, integrator mjINT_IMPLICITFAST
(1) fields each stage writes (filled with NaN beforehand)
mj_fwdPosition M, actuator_length, actuator_moment, cam_xmat, cam_xpos, cdof, cinert, crb, efc_D, efc_J, efc_KBIP, efc_R, efc_diagA, efc_frictionloss, efc_margin, efc_pos, geom_xmat, geom_xpos, iefc_D, iefc_R, iefc_frictionloss, light_xdir, light_xpos, qLD, qLDiagInv, site_xmat, site_xpos, subtree_com, ten_J, ten_length, xanchor, xaxis, ximat, xipos, xmat, xpos, xquat
mj_sensorPos (nothing complete); partly: sensordata (8 of 9)
mj_fwdVelocity actuator_velocity, cdof_dot, cvel, efc_aref, efc_vel, qfrc_adhesion, qfrc_bias, qfrc_damper, qfrc_fluid, qfrc_gravcomp, qfrc_passive, qfrc_spring, ten_velocity; partly: sensordata (8 of 9)
mj_sensorVel (nothing complete); partly: sensordata (8 of 9)
mj_fwdActuation actuator_force, qfrc_actuator; partly: sensordata (8 of 9)
mj_fwdAcceleration qacc_smooth, qfrc_smooth; partly: sensordata (8 of 9)
mj_fwdConstraint efc_b, efc_force, iacc, iacc_smooth, iefc_aref, iefc_force, ifrc_constraint, ifrc_smooth, qacc, qfrc_constraint; partly: sensordata (8 of 9)
mj_sensorAcc sensordata
qacc from the stages against mj_forward: max difference 0.0e+00
(2) fields no stage wrote
cacc, cfrc_ext, cfrc_int, efm_ak, efm_as, efm_ck, efm_diag, efm_fluid, efm_sdiag, efm_ts, energy, qDeriv, qH, qHDiagInv, qLU, qfrc_inverse, solver_fwdinv, subtree_angmom, subtree_linvel, wrap_xpos
(3) after one mj_step, with the arm's joints moving at 1 rad/s
body positions in mjData against those of the new qpos: max difference 3.631 mm; against those before the step: 0.000 mm
sensordata in mjData against the new state's: max difference 2.42e-02; against the state before the step: 0.00e+00
(4) mj_forwardSkip after a change of qvel only
2000 random velocity perturbations: max |qacc(skip) - qacc(full)| 0.0e+00; mj_forward 38.2 us, mj_forwardSkip(mjSTAGE_POS) 26.3 us
Positions or sensors read after mj_step are one step old. They are, by design. Call mj_forward (or at least mj_kinematics) before reading them, or read them before stepping.
A quantity is zero or stale although its sensor works. It is computed lazily. Call mj_rnePostConstraint for cacc and cfrc_*, or mj_subtreeVel for subtree_linvel and subtree_angmom.
Stages called by hand give a different result from mj_forward. The order or a stage is missing. mj_forwardSkip also runs the control callback, energy, and, for the discrete integrator, constraint regularization between the stages.
A field holds data from the previous step after mj_forwardSkip. You skipped a stage whose inputs changed. Skip the position stage only when qpos and the model are unchanged.
Run the stage experiment on push.xml after 1 s, and on a model with a fluid medium (set opt.density and opt.viscosity). Which fields appear that did not on pick_place, and in which stage?
Checkpoint 20: pick qfrc_passive or efc_force, find every function in src/engine that writes it, and confirm the order by instrumenting a run. Hint: qfrc_passive is assembled in mj_passive from qfrc_spring, qfrc_damper, qfrc_fluid, qfrc_adhesion and qfrc_gravcomp (engine_passive.c, lines 1018 to 1078). efc_force is first set by the warmstart (engine_forward.c, around lines 1094 to 1149) and then by the solver. Document one surprise.
[!research] Why the pipeline’s shape matters for learning Differentiable and batched simulators reuse this pipeline’s structure. MJX and MuJoCo Warp (Lesson 20.3) reimplement these stages for accelerators, and the finite-difference derivatives of Lesson 20.4 rely on
mj_forwardSkip. When a learning pipeline reads observations frommjDataright aftermj_step, it reads the previous state’s kinematics. That one-step lag is easy to miss and shows up as a fixed latency in a trained policy (inference; Lessons 12.1, 18.2 and 19.1 measure related effects).
{"id": "20.1-check", "title": "Knowledge check", "questions": [
{"kind": "mcq", "q": "Which stage writes qfrc_bias?",
"options": ["mj_fwdPosition", "mj_fwdVelocity", "mj_fwdActuation", "mj_fwdConstraint"],
"answer": 1,
"explain": "<p>mj_fwdVelocity runs RNE without acceleration (mj_rne with flg_acc 0) to get c(q, qdot).</p>"},
{"kind": "mcq", "q": "After mj_step returns, data.xpos describes:",
"options": ["the new qpos", "the state before the step", "an interpolation of both", "nothing until mj_forward"],
"answer": 1,
"explain": "<p>Measured: 0.000 mm from the state before the step, 3.631 mm from the new one.</p>"},
{"kind": "mcq", "q": "Why are cacc and cfrc_int still NaN after every stage?",
"options": ["They are inputs", "They are computed lazily, when a sensor needs them, via mj_rnePostConstraint", "They are unused in 3.14", "The integrator writes them"],
"answer": 1,
"explain": "<p>mj_forward clears flg_rnepost; the acceleration sensors call mj_rnePostConstraint only if needed.</p>"},
{"kind": "mcq", "q": "When is mj_forwardSkip(m, d, mjSTAGE_POS, 0) safe?",
"options": ["Always", "When only qvel, ctrl or applied forces changed since the last full forward", "When qpos changed", "Only with RK4"],
"answer": 1,
"explain": "<p>Everything the position stage computes depends on qpos and the model only.</p>"}
]}
Lesson 20.2 opens the last stage: the constraint solver.