Simulation Loop¶

STARK advances a simulation by repeatedly solving one nonlinear optimization problem per time step. This is triggered by the .run() method:

stark::Simulation simulation(settings);

// Build objects, materials, contact, boundary conditions, and output...

simulation.run();

Internally, each call to run_one_time_step() updates time-dependent data, runs Newton’s method via SymX, commits or retries the step, writes frame output, and logs solver statistics.

Optimization Problem¶

STARK does not solve directly for final positions. For deformable point sets, the primary unknowns are the next-step velocities v1. The next positions are reconstructed by first-order time integration:

x1 = x0 + dt * v1;

For example, a deformable object contributes one 3D velocity unknown per point:

soft.v1: 3267

Rigid bodies contribute linear and angular velocity unknowns:

rigid.v1: 3
rigid.w1: 3

All energies are written as functions of the predicted next state. A minimum of the global potential corresponds to a force balance, found by Newton’s method.

Initialization¶

The first time a simulation step is requested, STARK initializes the solver:

check that at least one degree of freedom was registered,
create the SymX Newton solver from the global potential,
compile/load the generated SymX kernels,
print the degrees of freedom,
write frame zero,
run before_simulation callbacks,
check that the initial state is valid.

High-level time loop¶

A simplified version of simulation.run() is:

while (!finished) {
    script.run_a_cycle(simulation.get_time());
    user_callback();              // optional callback passed to run(...)

    bool keep_running = stark.run_one_step();
    if (!keep_running) {
        break;
    }
}

stark.print_summary();

The loop stops when one of the configured limits is reached:

settings.execution.end_simulation_time
settings.execution.end_frame
settings.execution.allowed_execution_time
the duration passed to simulation.run(duration)
a callback requests termination
the solver fails in a non-recoverable way

Newton’s method inside each step¶

The nonlinear solve is handled by the SymX Newton solver. Check the original SymX docs for more details.

The simplified inner loop is:

check initial state validity

for each Newton iteration:
    run before_step callbacks
    run before_energy_evaluation callbacks
    evaluate energy, gradient, and local Hessians
    compute residual
    check residual convergence

    repeat:
        project local Hessians if needed
        assemble/update the global Hessian
        solve H * du = -grad
        check that du is a descent direction
        if not, increase projection and retry

    check step-size convergence

    line search:
        apply user step_cap
        apply max_allowed_step callbacks
        backtrack until intermediate state is valid
        backtrack until Armijo sufficient decrease holds

    run optional user convergence callbacks
    run after_step callbacks

check converged-state validity

Hessian projection¶

STARK commonly runs with positive-definite Hessian projection enabled. This makes the Newton system more robust for simulation problems whose local Hessians may become indefinite. Depending on the selected projection mode, SymX can:

run pure Newton without projection,
project all local Hessians,
project only after a failed linear solve or non-descending direction,
progressively project selected blocks based on the gradient magnitude.

The console field ph reports the percentage of projected local Hessians during a step.

Line search¶

The line search has several layers:

Stage	Meaning	Typical use
`cap`	Clamp the update by `settings.newton.step_cap`	Prevent overly large Newton steps
`max`	Clamp by `max_allowed_step` callbacks	Continuous collision detection and other hard step limits
`inv`	Backtrack while the intermediate state is invalid	Avoid inversions, penetrations, or other invalid states
`bt`	Armijo backtracking	Require sufficient energy decrease

These are the four counters printed as

ls (cap|max|inv|bt):  0| 0| 1| 4

Adaptive time step¶

STARK optionally allows for adaptive time step, which allows for automatically reducing the time step when the problem is challenging and otherwise keeps it at a user specified value.

The time step starts as

dt = settings.simulation.max_time_step_size;

After a successful accepted step, STARK grows it by

dt = min(max_time_step_size, dt * time_step_size_success_multiplier);

If Newton fails because the step is too difficult (e.g. maximum Newton iteration reached), STARK halves dt and retries the same simulation time:

dt *= 0.5;

If dt drops below settings.simulation.time_step_size_lower_bound, the simulation exits.

Callbacks¶

Callbacks inject behavior into the solver loop. Most users only need the scripting interface, but lower-level callbacks are useful when extending STARK.

Callbacks are accessed through simulation.get_callbacks(), which returns a reference to stark::Callbacks. That struct exposes the STARK-level methods directly, and a newton member of type symx::SolverCallbacks that carries the SymX/Newton-level hooks.

STARK-level callbacks¶

These callbacks relate to time stepping.

Callback	When it runs	Purpose
`before_simulation`	once during initialization, after frame zero is written	initialize derived data
`before_time_step`	before Newton starts	reset unknowns, update forces, update time-step data
`on_time_step_accepted`	only after Newton succeeds and the state is accepted	commit the solved state
`after_time_step`	after accepting the step, before advancing output/logging is complete	post-step updates
`write_frame`	whenever a frame should be written	write VTK/output data
`should_continue_execution`	before each attempted step	allow systems/users to stop the simulation

The frame-writing callback is how STARK’s built-in output system writes registered meshes every frame. The frame cadence is controlled by settings.output.fps.

simulation.get_callbacks()->add_before_time_step([&]() {
    // runs once before Newton starts for every time step
    // e.g. update step-lagged quantities
});

SymX/Newton callbacks¶

These callbacks relate to the Newton solves. Check the original SymX docs for more details.

Callback	When it runs	Purpose
`before_step`	at the very start of each Newton iteration	logging, resetting per-iteration state
`before_energy_evaluation`	before each energy or gradient/Hessian evaluation	update collision sets, derived fields, or time-dependent data
`is_initial_state_valid`	before Newton starts	reject invalid starting states
`max_allowed_step`	before line-search backtracking	impose step limits, commonly from CCD
`is_intermediate_state_valid`	during line search	reject invalid intermediate states
`on_intermediate_state_invalid`	after too many invalid-state backtracks	harden constraints/contact or report failure
`on_armijo_fail`	after too many Armijo backtracks	diagnostics or model adjustment
`is_converged`	after a line-search update	add custom convergence criteria
`is_converged_state_valid`	after Newton reports success	verify final constraints/contact tolerances
`after_step`	at the end of each Newton iteration	per-iteration diagnostics, adaptive parameter updates

Notes:

Boolean validity callbacks are combined with logical AND. All registered checks must pass.
max_allowed_step callbacks are combined by taking the smallest returned step fraction.
The default residual used by Newton is the infinity norm of the gradient, unless a custom residual callback is installed.

SymX/Newton callbacks live on simulation.get_callbacks()->newton, which is of type symx::spSolverCallbacks (defined in SymX):

simulation.get_callbacks()->newton->add_before_energy_evaluation([&]() {
    // runs before every energy / gradient / Hessian evaluation inside Newton
    // e.g. update broad-phase collision sets
});

Public scripting and control¶

For ordinary scene control, use the public scripting interface rather than low-level core callbacks.

Time events¶

A time event runs while t0 <= simulation.get_time() < t1:

simulation.add_time_event(0.0, 2.0,
    [&](double t)
    {
        const double z = stark::blend(0.0, 1.0, 0.0, 2.0, t, stark::BlendType::Linear);
        boundary.set_transformation({0.0, 0.0, z}, Eigen::Matrix3d::Identity());
    }
);

This is the most convenient way to animate boundary conditions, gravity, material parameters, friction coefficients, target positions, rigid-body motors, or other time-dependent scene data.

Events with state¶

For more advanced cases, the event callback can receive an EventInfo object:

simulation.add_time_event(1.0, 3.0,
    [&](double t, stark::EventInfo& event)
    {
        if (event.is_first_call()) {
            event.msg = "activated";
        }

        // event.get_n_calls(), event.get_begin_time(), event.pack(...), event.unpack(...)
        // can be used to keep small pieces of state between calls.
    }
);

For full control, simulation.get_script().add_event(...) exposes the event system directly with custom run_when and delete_when predicates.

Simple run callback¶

simulation.run(callback) executes callback before every attempted time step:

simulation.run(
    [&]()
    {
        const Eigen::Vector3d center = magnet.rigidbody.get_translation();

        for (int i = 0; i < object.point_set.size(); ++i) {
            const Eigen::Vector3d x = object.point_set.get_position(i);
            const Eigen::Vector3d u = center - x;
            object.point_set.set_force(i, 0.1 * u.normalized() / u.squaredNorm());
        }
    }
);

Use this for simple procedural control that should run throughout the simulation. For time-bounded behavior, prefer add_time_event.

Manual stepping¶

For interactive applications or external control loops, use:

while (viewer_is_open) {
    simulation.run_one_time_step();
    // update viewer, UI, custom data, etc.
}

This executes one script cycle and then attempts one STARK time step.

Quasistatic simulations¶

STARK supports quasistatic problems, though not as a first-class mode: it still solves for velocity-like unknowns using the same time-step machinery. The trick is to choose parameters so the velocity unknown behaves as a positional increment.

The quasistatic_column_extrusion() example uses the following key changes:

const double duration = 1.0;
const double dt = duration * 0.99999;  // Do not overflow the sim duration in a single step
const bool is_quasistatic = true;

settings.execution.end_simulation_time = duration;
settings.simulation.max_time_step_size = dt;
settings.output.fps = 1.0 / dt;

settings.newton.projection_mode = symx::ProjectionToPD::ProjectedNewton; // Advisable
settings.newton.project_to_pd_use_mirroring = true; // Advisable

stark::Volume::Params material;
material.inertia.quasistatic = true;  // Important!