Building a Physics Engine From Scratch in C++20: A Step-by-Step AAA Blueprint

There are two ways to build a physics engine. The first is the one you see in most weekend tutorials: start with a circle and a velocity vector, add gravity, then add a second circle, then discover that circles overlap, then add a dot product, then spend five days trying to stop a stack of boxes from vibrating itself into the floor. Zanzlanz's "Coding a Physics Engine from scratch!" is an honest record of exactly that journey. It is valuable as a diary of learning and a testament to persistence, but it is also a catalog of anti-patterns: energy injection from constraint drift, no fixed timestep, no warm-starting, no broadphase, no formal solver, and a heavy reliance on trial and error.

The second way is to build the pipeline as a composition of well-understood operators. This post is that blueprint. It is based on the V Engine's production physics system and on our Level-2 manifold physics paper and Level-3 deterministic accumulator paper. The goal is a AAA-grade rigid-body engine in C++20: stable stacks, deterministic behavior, parallel where safe, and mathematically grounded where it matters.

What this post is and is not

This is a step-by-step implementation blueprint, not a complete source dump. It assumes you are comfortable with C++20, glm for vector math, and the idea that a physics system is one part of a larger Entity Component System. Every stage below has a corresponding implementation in engine/systems/physics/ and a formal justification in our papers. The math is shown in inline diagrams; the pseudocode is shown in plain blocks; the tuning parameters are shown in a reference table.

This post is also an anti-design study. Where a naive approach would fail, we point to the Zanzlanz video as a real example of how that failure looks. The critique is not of the author — the author is learning in public, which is admirable — but of the design trajectory that the video represents. If you are building an engine for a shipped game, you want the trajectory in this post.

The whole pipeline in one sentence

The fixed-timestep update is a composed operator on an extended state (bodies, manifold_cache):

In our Level-2 paper this is written as a single composite map:

If you run that operator s times per frame with timestep Δt/s, you get the frame-level operator. That is the architecture. Everything else is engineering inside the stages.

Step 1: The rigid body as a component, not a class

Before writing integration, decide the data layout. In a modern C++20 engine using an ECS, a body is a set of plain components:

enum class BodyType { Dynamic, Static, Kinematic };

struct RigidBodyComponent {
    BodyType type = BodyType::Dynamic;

    float mass = 1.0f;
    float inv_mass = 1.0f;
    glm::vec3 inertia = {1.0f, 1.0f, 1.0f};
    glm::vec3 inv_inertia = {1.0f, 1.0f, 1.0f};

    glm::vec3 velocity = {0.0f, 0.0f, 0.0f};
    glm::vec3 previous_velocity = {0.0f, 0.0f, 0.0f};
    glm::vec3 angular_velocity = {0.0f, 0.0f, 0.0f};
    glm::vec3 previous_angular_velocity = {0.0f, 0.0f, 0.0f};

    glm::vec3 force = {0.0f, 0.0f, 0.0f};
    glm::vec3 torque = {0.0f, 0.0f, 0.0f};

    float linear_damping = 0.01f;
    float angular_damping = 0.15f;
    float velocity_smoothing = 0.15f;
    float restitution = 0.1f;
    float friction = 0.9f;
    bool use_gravity = true;
    bool is_sleeping = false;
    bool is_grounded = false;

    // Sleep counters and support hints (simplified)
    int slow_frames = 0;
    uint32_t wake_cooldown_frames = 0;
    uint32_t pose_still_frames = 0;
    float last_support_dot_up = -1.0f;

    // Axis constraints (optional)
    bool freeze_position_x = false;
    bool freeze_position_y = false;
    bool freeze_position_z = false;
    bool freeze_rotation_x = false;
    bool freeze_rotation_y = false;
    bool freeze_rotation_z = false;
};

struct TransformComponent {
    glm::vec3 position = {0.0f, 0.0f, 0.0f};
    glm::quat orientation = {1.0f, 0.0f, 0.0f, 0.0f};
    glm::vec3 scale = {1.0f, 1.0f, 1.0f};
    bool matrix_dirty = true;
};

enum class ColliderType { Box, Sphere, Mesh };

struct ColliderComponent {
    ColliderType type = ColliderType::Sphere;
    glm::vec3 size = {0.5f, 0.5f, 0.5f}; // half-extents for box/mesh; x=radius for sphere
    glm::vec3 offset = {0.0f, 0.0f, 0.0f};
    bool is_trigger = false;
};

// For convex triangle-mesh colliders:
struct MeshColliderComponent {
    std::string mesh_path;
    std::vector<glm::vec3> vertices;
    std::vector<uint32_t> indices;
    bool convex = true;
};

The key constraint is that components contain data only. No methods, no inheritance, no hidden logic. The physics system owns the update operator and acts on contiguous views of these components. This is the pattern used in v::engine::systems::physics::PhysicsSystem: ECS snapshots are captured on the main thread, then parallel jobs operate on plain arrays.

entt::registry registry;
entt::entity ball = registry.create();

RigidBodyComponent rb;
rb.type = BodyType::Dynamic;
rb.mass = 1.0f;
rb.inv_mass = 1.0f / rb.mass;
rb.inertia = {0.4f, 0.4f, 0.4f};
rb.inv_inertia = {1.0f/0.4f, 1.0f/0.4f, 1.0f/0.4f};
rb.restitution = 0.1f;
rb.friction = 0.3f;
rb.linear_damping = 0.01f;
rb.angular_damping = 0.15f;

TransformComponent t;
t.position = {0.0f, 5.0f, 0.0f};
t.orientation = {1.0f, 0.0f, 0.0f, 0.0f};
t.scale = {1.0f, 1.0f, 1.0f};

ColliderComponent c;
c.type = ColliderType::Sphere;
c.size = {0.5f, 0.5f, 0.5f}; // x holds the radius for spheres

registry.emplace<RigidBodyComponent>(ball, rb);
registry.emplace<TransformComponent>(ball, t);
registry.emplace<ColliderComponent>(ball, c);

// Static ground:
entt::entity ground = registry.create();
RigidBodyComponent ground_rb;
ground_rb.type = BodyType::Static;
ground_rb.inv_mass = 0.0f;
ColliderComponent ground_c;
ground_c.type = ColliderType::Box;
ground_c.size = {10.0f, 0.5f, 10.0f}; // half-extents

TransformComponent ground_t;
ground_t.position = {0.0f, -0.5f, 0.0f};

registry.emplace<RigidBodyComponent>(ground, ground_rb);
registry.emplace<TransformComponent>(ground, ground_t);
registry.emplace<ColliderComponent>(ground, ground_c);

Step 2: Fixed timestep and symplectic Euler

Games do not run at a stable wall-clock interval. A frame may take 16 ms, 8 ms, or 33 ms. If you integrate with variable dt, your simulation is a function of frame pacing, which is a function of GPU load. The first hard rule of a shippable engine is: integrate at a fixed timestep, subdivided into substeps, and clamp or drop excess time.

Symplectic Euler is the standard game integrator. It is "velocity-first":

void IntegrateBody(body, transform, dt):
    if body.is_sleeping or body.type != Dynamic:
        return

    acceleration = gravity + body.force * body.inv_mass
    body.velocity += acceleration * dt
    body.angular_velocity += body.torque * body.inv_inertia * dt

    body.velocity *= exp(-linear_damping * dt)
    body.angular_velocity *= exp(-angular_damping * dt)

    body.velocity = clamp(body.velocity, -max_linear, max_linear)
    body.angular_velocity = clamp(body.angular_velocity, -max_angular, max_angular)

    transform.position += body.velocity * dt
    q_dot = 0.5f * quat(0, body.angular_velocity) * transform.orientation
    transform.orientation = normalize(transform.orientation + q_dot * dt)

The quaternion derivative uses q_dot = 0.5 * ωq ⊗ q, where ωq = (0, ω). The final normalize is essential: without it the quaternion drifts off the unit sphere and the body begins to stretch. In the V Engine this is done inside physics_operators.hpp and called from the jobified integration path.

Step 3: Damping, clamps, and safety belts

A real-time engine must survive extreme inputs: teleporting an object deep inside another, tiny masses colliding with huge ones, high restitution near a wall. The safety system is simple and must be applied before the solver:

void ApplySafetyBelts(body, dt):
    // Exponential damping removes energy predictably.
    body.velocity *= exp(-c_linear * dt)
    body.angular_velocity *= exp(-c_angular * dt)

    // Clamp velocities to a maximum safe magnitude.
    body.velocity = clamp_magnitude(body.velocity, v_max)
    body.angular_velocity = clamp_magnitude(body.angular_velocity, omega_max)

    // Optionally snap tiny velocities to zero when below a settling threshold.
    if length(body.velocity) < settle_linear and length(angular) < settle_angular:
        body.velocity *= 0.95

These are not just heuristics. In the Level-2 paper they are part of the contraction proof: exponential damping contributes a multiplicative gap factor κ = e−cΔt, and clamping is identified with the projection kernel, which is non-expansive. The math is exactly the code.

Step 4: Broadphase — the uniform spatial grid

Why broadphase?

If every object tested every other object, pair generation would be O(N2). At 1,000 bodies that is a million tests, far too many for a 16 ms frame. Broadphase separates potentially colliding pairs from pairs that are clearly too far apart. The goal is to produce a small candidate set that the expensive narrowphase can actually process.

What is an AABB?

An Axis-Aligned Bounding Box (AABB) is the smallest box with sides parallel to the world axes that fully contains a body. For a sphere the AABB is a cube whose side is the diameter. For a rotated box the AABB is computed from the rotated half-extents, so the AABB is typically larger than the box itself. The broadphase only ever sees AABBs; it does not know the real shape.

Grid cell insertion

Choose a cell size h roughly the size of a typical body. For each body, compute the AABB, divide the min and max corners by h, and insert the body into every integer cell in that range. A large body can appear in many cells; a small body usually appears in one. This is why the data structure is sometimes called a coarse grid: it deliberately over-approximates the set of candidate pairs.

Why are there two cells with the same body in Figure 2b? Because the AABB is larger than the visual circle. The broadphase does not care about the visual shape; it only cares about the conservative AABB. If an object moves slowly, it will usually stay in the same set of cells from frame to frame, which keeps the candidate pair list stable and helps the manifold cache warm-start.

struct CellCoord { int x, int y, int z; };

void BuildBroadphaseGrid(WorldState, float cell_size h):
    for body in WorldState.bodies:
        aabb = ComputeWorldAABB(body)
        i0 = floor(aabb.min.x / h); i1 = floor(aabb.max.x / h)
        j0 = floor(aabb.min.y / h); j1 = floor(aabb.max.y / h)
        k0 = floor(aabb.min.z / h); k1 = floor(aabb.max.z / h)
        for i in [i0, i1]:
            for j in [j0, j1]:
                for k in [k0, k1]:
                    grid[CellCoord(i,j,k)].push_back(body)

void GeneratePairs(WorldState):
    pair_set = empty_set()
    for cell, list in grid:
        for each unordered pair (A,B) in list:
            if CanCollide(A,B):
                pair_set.insert(NormalizePairKey(A,B))
    return sorted(pair_set)  // deterministic order

Determinism requires sorting the final pair list. Hash iteration order can differ between runs. In the V Engine this is handled by PhysicsSystem::generate_pairs() and the SpatialGrid implementation. For strict reproducibility, every iteration order in the engine must be a deterministic function of the input state.

What is a manifold? — a visual explainer

Before narrowphase runs, it helps to understand what it is producing. A manifold is a small data packet that summarizes one collision between two bodies. It is not the raw mesh overlap and it is not a list of every intersecting triangle; it is the minimum information the solver needs to push the two bodies apart without letting them tunnel or explode.

Why a "receipt"? Because once the narrowphase has done its expensive geometry work, the solver should not have to think about triangles again. The manifold records:

• entity_a, entity_b: the two participants. Order is normalized so the pair key is the same whichever body is found first.
• normal: the unit vector pointing from A into B. This is the direction the push must go if the bodies are interpenetrating.
• penetration: how far A has moved into B along the normal. The solver uses this to create the Baumgarte velocity bias.
• contact_point: a single world-space point where the push is applied. This point determines the lever arms r_a and r_b and therefore the torque on each body.
• restitution, friction: material properties. They may be combined from the two bodies (restitution is usually the minimum, friction is the geometric mean in the V Engine).
• up to 4 CachedContact slots: for box-on-box or convex-mesh contacts, one contact point is not enough; the solver needs a small patch so the box does not rock on a single point.

Contact force chains and entangled collisions

A single manifold is local: body A and body B. But bodies are usually connected into chains. A stack of three boxes is not three independent manifolds; it is a network where each manifold pushes on its neighbors and the forces propagate down to the ground. This is the physical meaning of "entangled" collisions: the manifold graph is the medium through which gravity, support, and friction forces travel.

This is why warm-starting and contact islands matter so much. If manifold AB is solved first, it may think it is the only contact and push A down into the ground. The solver must iterate over the whole island many times, passing impulses back and forth, before the network converges. A stack that looks still to the eye is actually a solved equilibrium of many tiny forces.

Step 5: Narrowphase — from pairs to manifolds

Broadphase gives candidate pairs. Narrowphase decides which pairs actually overlap and, if they do, produces a contact manifold: a small set of contact points (usually 1–4) that approximate the local collision geometry.

A manifold stores, at minimum:

struct CollisionManifold {
    entt::entity entity_a;
    entt::entity entity_b;
    glm::vec3 normal;          // points A -> B
    float penetration;         // overlap depth
    glm::vec3 contact_point;   // world-space push point
    float restitution;
    float friction;
};

struct CachedContact {
    glm::vec3 local_point_a;   // on body A, in local space
    glm::vec3 local_point_b;   // on body B, in local space
    glm::vec3 world_point;     // in world space
    glm::vec3 normal;          // points A -> B
    glm::vec3 tangent_1;       // first friction direction
    glm::vec3 tangent_2;       // second friction direction
    glm::vec3 r_a;             // world_point - position_a
    glm::vec3 r_b;             // world_point - position_b
    float penetration;
    float normal_impulse;
    float tangent_impulse_1;
    float tangent_impulse_2;
    float normal_mass;
    float tangent_mass_1;
    float tangent_mass_2;
    float velocity_bias;
    uint8_t lifetime;
    bool matched;
};

constexpr uint8_t MAX_CONTACTS_PER_MANIFOLD = 4;

struct CachedManifold {
    entt::entity entity_a;
    entt::entity entity_b;
    std::array<CachedContact, MAX_CONTACTS_PER_MANIFOLD> contacts;
    uint8_t contact_count = 0; // 1..4
    float friction;
    float restitution;
    uint32_t lifetime;
    bool updated_this_frame;
};

For box-box or box-plane, the contact points should cover the contact patch. The most stable policy is to keep the points with the highest normal impulse (the load-bearing contacts) or the points that form the best support polygon. The V Engine uses ContactReductionPolicy::KeepStable by default, with KeepCorners and KeepDeepest as alternatives.

Narrowphase shapes in the V Engine

The current V Engine supports three collider types: sphere, box, and convex mesh. The actual implementation in engine/systems/physics/physics_job.cpp handles each combination:

• Sphere-sphere: The normal is the center-to-center vector. Penetration is the sum of radii minus the distance. If the centers coincide, it falls back to world-up to avoid a zero normal.
• Sphere-box: Transform the sphere center into the box's local space, clamp to the box half-extents to find the closest point, then compute distance and normal. This also handles the case where the sphere center is inside the box by picking the nearest face.
• Box-box: Uses the Separating Axis Theorem (SAT). The candidate separating axes are the 6 face normals of the two boxes plus the 9 cross products of their edge directions (15 axes total). The V Engine batches these 15 axes with SIMD to improve throughput. For each axis it projects both boxes onto the axis and checks for an overlap gap. If any axis has a gap, the boxes are not overlapping. The smallest overlap across all axes gives the penetration direction.
• Box-mesh and mesh-mesh: Convex mesh colliders are stored in MeshColliderComponent with vertex and index lists. The mesh is treated as a collection of triangles and run through the same SAT loop, with the mesh face normals and edge directions used as candidate axes. This is enabled by the collider_type = "mesh" archetype in levels such as gym_physics.toml.

Contact point and material combination

Once the best axis is found, the engine computes the contact point from the support points of the two shapes along ±normal. For ground-like contacts (normal nearly aligned with gravity-up) the point is biased toward the upper body to avoid underestimating the moment arm and creating artificial edge-balanced equilibria. If one body is a huge static (e.g. the ground plane) and the other is a small dynamic, the contact point is kept near the dynamic body so the large static's corners do not create ghost contacts far from the real support.

Material properties are combined per manifold:

• Restitution:min(body_a.restitution, body_b.restitution), unless a static restitution override is configured.
• Friction:sqrt(body_a.friction * body_b.friction), the geometric mean, which is the standard Coulomb combination for two rough surfaces.

This combination is not a stylistic choice; it is what makes the friction cone consistent between two bodies of very different materials. If one body is ice and the other is rubber, the effective friction must be low.

Step 6: The manifold cache — the secret to stable stacks

Here is the most important architectural decision. Do not throw away your contact data at the end of the frame. Store it. The persistent manifold cache is keyed by the normalized entity pair and holds up to four contacts per pair. Each cached contact holds:

• Local anchors on each body
• World contact point and normal
• Tangent basis (t1, t2)
• Accumulated normal impulse λn
• Accumulated tangential impulses (λt1, λt2)
• Effective mass terms and velocity bias
• Lifetime and matched flag

On the next frame, match newly detected contacts to existing slots by comparing local anchors. If a match exists, reuse the accumulated impulse as the initial guess. This is warm-starting, and it is why stable stacks settle in a few frames instead of wobbling forever.

Impulses are not gameplay state. They are the solution guess of the constraint solver. Reusing them is the difference between a simulation that looks like glue and a simulation that looks like physics.

Step 7: The constraint solver — PGS sequential impulses

With bodies integrated, pairs generated, manifolds detected, and the cache matched, the heart of the engine is the solver. We use Projected Gauss–Seidel (PGS) sequential impulses. The idea is to visit each contact repeatedly and apply a small impulse that pushes the relative velocity toward the desired target.

For each contact direction d (normal or one of the two tangents), compute the relative velocity at the contact point, compare it to a target that includes a bias term, multiply by the effective mass, and add the resulting Δλ to the accumulated impulse.

void SolveContact(body_a, body_b, cached_contact, dt):
    r_a = cached_contact.world_point - transform_a.position
    r_b = cached_contact.world_point - transform_b.position

    // Relative velocity at the contact point.
    v_rel = (body_b.velocity + cross(body_b.angular_velocity, r_b))
          - (body_a.velocity + cross(body_a.angular_velocity, r_a))

    // Normal constraint.
    vn = dot(v_rel, cached_contact.normal)
    d_lambda_n = -(vn + cached_contact.velocity_bias) * cached_contact.normal_mass
    new_lambda_n = max(0.0f, cached_contact.normal_impulse + d_lambda_n)
    d_lambda_n = new_lambda_n - cached_contact.normal_impulse
    cached_contact.normal_impulse = new_lambda_n
    ApplyImpulsePair(body_a, body_b, r_a, r_b, d_lambda_n * cached_contact.normal)

    // Friction constraints (two tangents).
    vt1 = dot(v_rel, cached_contact.tangent_1)
    vt2 = dot(v_rel, cached_contact.tangent_2)
    d_lambda_t1 = -vt1 * cached_contact.tangent_mass_1
    d_lambda_t2 = -vt2 * cached_contact.tangent_mass_2

    new_lt1 = cached_contact.tangent_impulse_1 + d_lambda_t1
    new_lt2 = cached_contact.tangent_impulse_2 + d_lambda_t2

    // Coulomb friction disk: ||λ_t|| ≤ μ * λ_n
    max_friction = cached_contact.friction * cached_contact.normal_impulse
    new_lt1, new_lt2 = ProjectToDisk(new_lt1, new_lt2, max_friction)

    d_lambda_t1 = new_lt1 - cached_contact.tangent_impulse_1
    d_lambda_t2 = new_lt2 - cached_contact.tangent_impulse_2
    cached_contact.tangent_impulse_1 = new_lt1
    cached_contact.tangent_impulse_2 = new_lt2
    ApplyImpulsePair(body_a, body_b, r_a, r_b,
                     d_lambda_t1 * cached_contact.tangent_1
                   + d_lambda_t2 * cached_contact.tangent_2)

The ApplyImpulsePair function is the physical law:

void ApplyImpulsePair(body_a, body_b, r_a, r_b, impulse):
    if body_a is dynamic active:
        body_a.velocity -= impulse * body_a.inv_mass
        body_a.angular_velocity -= body_a.inv_inertia * cross(r_a, impulse)
    if body_b is dynamic active:
        body_b.velocity += impulse * body_b.inv_mass
        body_b.angular_velocity += body_b.inv_inertia * cross(r_b, impulse)

Step 8: Effective mass and the velocity bias

The solver needs two numbers at each contact: the effective massd(d) for direction d, and the velocity bias that tells the contact what velocity it should target.

Effective mass converts a velocity change into an impulse:

In code, you compute the angular Jacobian term r × d, multiply by inverse inertia, and cross back with r. This gives the resistance the body offers to an impulse at that point in that direction.

The velocity bias has two parts:

1. Baumgarte stabilization adds a small extra velocity to push bodies out of penetration over time: bias = β · max(0, penetration − slop) / dt. Larger β resolves overlap faster but can inject energy and jitter.
2. Restitution adds bounce only for fast impacts: if vn < −vbounce, bias += −e · vn. For slow resting contacts, restitution is suppressed.

The combined bias is the target the solver aims at. The normal update becomes Δλn = −(vn + bias) · d(n).

Step 9: Contact islands and the sleep system

Once the solver is stable, the next problem is waste. A stack of a hundred sleeping boxes should cost almost nothing. The engine groups bodies into contact islands: connected components of the contact graph. An island can sleep as a unit if all its members are still and supported.

Sleep criteria in the V Engine are multi-gated:

• Pose is still for a consecutive number of frames.
• Linear and angular velocity are below tuned thresholds.
• The body has a support contact (a contact normal pointing upward-ish).
• No active neighbor in the same island is awake.
• A wake-cooldown hysteresis prevents rapid wake/sleep oscillation.

When a body sleeps, its velocity is zeroed and it is skipped by integration and the solver until an impulse, teleport, or moving overlap wakes it. This is a projection onto a closed set, which the Level-2 proof shows is non-expansive.

Axis constraints

Axis constraints are per-body freeze flags on position and rotation components. In the real RigidBodyComponent these are freeze_position_x/y/z and freeze_rotation_x/y/z. During apply_axis_constraints() the solver simply zeroes the corresponding velocity, force, and torque components. A common pattern is the player-style constraint: freeze X and Z rotation but leave Y (yaw) free; the system then converts any accumulated pitch/roll into a yaw-only quaternion so the character never falls over sideways. These constraints are projections onto an admissible subspace, so they are non-expansive in the contraction sense and help keep the player body readable.

Righting, settling, and the sleep stage

The job graph still has a legacy RightingCorrections stage, but in the current V Engine it is a no-op: all righting is now folded into update_sleep_states() so that the same code evaluates sleep, stabilizes, and aligns. When a box is nearly face-up but tilted a few degrees, the sleep evaluator applies a small righting angular acceleration around the axis that would bring the nearest face normal into alignment with gravity-up. This is clamped and only applied when the body has a contact hint and is not already face-aligned. It is a targeted, non-explosive correction that keeps dense stacks from resting on edges.

Step 10: Determinism and the Level-3 proof sketch

Determinism is not optional for AAA networked gameplay, replay systems, or debugging. It requires three things:

1. Stable iteration order: sort pair lists, sort manifold lists, and never iterate unordered hash maps during simulation.
2. Fixed arithmetic mode: same compiler flags, same floating-point settings, same order of operations.
3. No parallel writes to shared constraint state: the solve stage remains serial, while broadphase and narrowphase can be parallel over independent chunks.

The Level-3 paper formalizes the parallel part. Worker-local manifold lists are merged by ordered concatenation. The per-body impulse accumulator is a sum over the flattened contact list, and because integer/floating-point addition is associative and commutative up to the chosen rounding mode, the total contribution is independent of chunk order. In the Lean 4 module this is Level3PhysicsSystem.totalContribution_bulkMerge_chunks.

For the C++ engine, the practical takeaway is: make the merged contact list deterministic, then make the solver iterate it in that deterministic order. The GPU compute path in physics_system_gpu.cpp and physics_compute_pipelines.cpp follows the same invariant by writing worker-local results and then bulk-appending them.

Level-3 physics in action

The pipeline described in this post is the same one shown in the V Engine's Level-3 physics demo. The video below is a capture of that build from about six months ago. It is a good reference for what this blueprint produces, but it is not the latest internal state: the engine has since gained convex mesh colliders and other gameplay features that are part of the in-progress Level-4 work.

Level-3 physics in action on the M7 demo. Video: p3nGu1nZz, "GPU Physics Level 3 | Custom Vulkan 3D Game Engine in C++20 | V Engine | ECS | M7 Demo." Note: this is the proven Level-3 pipeline, not the still-in-research Level-4 system.

Step 11: GPU compute integration

When the body count is high and the GPU is underutilized, the integration and manifold-cache update stages can be moved to Vulkan compute shaders. The V Engine uses two compute pipelines:

• physics_integrate.comp: integrates dynamic bodies in parallel.
• physics_update_cache.comp: updates cached contact geometry from new transforms.

The CPU still owns the solver, because the PGS iteration needs deterministic ordering and each contact modifies shared body state. The GPU path is a performance additive, not a replacement. The compute shaders are loaded relative to the executable directory and compiled from assets/shaders/ by CMake, not by hand.

Tuning parameters and their meaning

The engine has roughly a dozen knobs that control stability, cost, and look. The table below maps each symbol to a typical range and a rule of thumb.

These constants live in engine/systems/physics/physics_constants.hpp and are loaded from config/physics.toml at runtime. Centralizing them makes the engine tunable without recompilation and makes the tuning visible to QA and designers.

What the numbers look like in practice

The V Engine physics system reports per-stage timings through PhysicsSystem::StageTimings and query/sleep telemetry. The benchmarks in docs/bench/data/physics_system.csv give a concrete sense of the cost:

Those numbers are from the linux-debug preset, so they include debug symbols and assertions. Release builds will be faster. The point is not the absolute magnitude — it is that the cost is measurable per stage, which makes the bottleneck visible. If the solver suddenly dominates, you raise the iteration count. If broadphase dominates, you enlarge the cell size. If narrowphase dominates, you simplify shapes or reduce manifold slots. Without stage timings, you are tuning by feel. The StageTimings struct in physics_system.hpp tracks integrate, grid build, pair generation, narrowphase, manifold cache, island build, warm start, solve, righting, axis constraints, and sleep as separate counters.

What the 7-day video gets wrong: a constructive post-mortem

Zanzlanz's video is a faithful record of learning physics by building it. The author explicitly says most of the time was spent struggling with rigid-body physics. That struggle is instructive, but the design choices that cause it are the exact ones a production engine avoids.

The dot product is indeed the key math tool the author celebrates, and it is also the first math tool of a good engine. The difference is that in a good engine the dot product is embedded inside a pipeline, not used as a replacement for one.

The blueprint as a single implementation skeleton

Putting the pieces together, a minimal C++20 physics system has the following skeleton. This is not the full V Engine source, but it is the shape every correct implementation takes:

void PhysicsSystem::fixed_step(float dt):
    substeps = clamp(floor(dt / fixed_timestep), 1, max_substeps)
    sub_dt = dt / substeps

    for i in 0..substeps-1:
        integrate_parallel(sub_dt)
        build_broadphase_grid()
        pairs = generate_pairs()
        manifolds = narrowphase(pairs)
        update_manifold_cache(manifolds)
        build_contact_islands()
        warm_start()
        solve_constraints_iterative(sub_dt)
        apply_righting_corrections(sub_dt)  // currently a no-op; righting is in sleep
        apply_axis_constraints()
        prune_stale_manifolds(3)
        update_sleep_states(sub_dt)         // righting and sleep happen together here

Each stage is a function with a deterministic contract. The extended state (bodies, cache) is transformed by the composite operator ΦΔt. The Level-2 paper proves contraction under explicit boundedness and non-expansiveness assumptions; the Level-3 paper proves that parallel chunking does not change the final accumulator.

Public-domain historical reference

The idea of integrating physical laws at discrete time steps is older than the games industry. Newton's Philosophiæ Naturalis Principia Mathematica (1687, English translation 1729) laid out the continuous laws we now discretize. The diagram below is from the 1729 edition and shows the same projectile-and-gravity problem that a modern engine solves every frame.

The difference between 1729 and now is not the physics. It is the algorithmic craft that makes the discretization stable, deterministic, and fast enough to run at 60 Hz.

The current V Engine: Level-4 is coming

As of the current build, the V Engine runs the Level-3 pipeline described here and is actively extending toward Level-4. The Level-4 additions include convex dynamic mesh colliders (already present as MeshColliderComponent and used in the gym_physics level), richer sleep/righting integration, and more gameplay-driven queries. Those extensions are still being benchmarked, tuned, and formalized; no Level-4 paper exists yet. When you read the source in engine/systems/physics/, the proven, documented pipeline is the Level-3 one; the Level-4 code paths are the active research frontier.

Caveats and limits

This blueprint covers impulse-based rigid-body dynamics for convex shapes. It does not cover everything an AAA engine might need:

• Deformable bodies, cloth, and soft bodies are a different family of solvers.
• Continuous collision detection (CCD) is needed for very fast thin objects to prevent tunneling; the V Engine currently assumes swept or clamped motion keeps objects within the collision envelope. Fast, thin projectiles can still tunnel in the current build.
• Articulated bodies, motors, and joints require additional constraint types beyond contact.
• Exact friction cones and bilateral constraints may need a non-linear complementarity solver rather than PGS.
• The proofs in the Level-2 and Level-3 papers assume bounded reachable sets and non-expansive projections; they do not prove correctness for unbounded inputs or for every possible GPU memory model.
• Convex mesh colliders are present in the current V Engine but they are still being hardened for Level-4. A convex-mesh-vs-mesh SAT pass is much more expensive than box-vs-box and needs careful culling and LOD to ship in dense scenes.
• GPU compute is treated here as a performance additive. Determinism on GPU is best-effort because different vendors, drivers, and warp schedules can introduce rounding or scheduling differences. For strict replay determinism, keep the solver on the CPU and use the GPU only for integration and cache-update stages that are independent and easy to verify.
• The inline SVG diagrams use CSS animations and SVG SMIL. All modern browsers support SMIL, but some CSPs or accessibility tools may strip it. The diagrams are fallbacks; the prose and pseudocode contain the same information.
• The public-domain images in Figure 1 and Figure 7 are historical references. They illustrate the continuous problem, not the discrete algorithms this post describes.

Cross-references and further reading

Final thought

Building a physics engine from scratch in a weekend is a great way to learn what the dot product is for. Building one that ships in a game is a lesson in architecture: fixed timesteps, persistent caches, projected solvers, and deterministic operators. The V Engine's physics system is the result of iterating on those lessons, and the Level-2 and Level-3 papers are the formal record of why they work. Use this post as the blueprint, the papers as the proof, and the 7-day video as the caution.