Dissipation leads to inconsistent velocity before and after restart

[sy-cui]

Describe the bug
By using load_state immediately after save_state, one should expect to recover the exact same simulator state. However, if the AnalyticalDamper is introduced into the simulator, a small difference occurs after loading state.

This issue may not be related AnalyticalDamper but

To Reproduce

import elastica as ea
import numpy as np


class RestartTestSimulator(
    ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping
):
    ...
restart_test_simulator = RestartTestSimulator()
cosserat_rod = ea.CosseratRod.straight_rod(
    n_elements=20,
    start=np.zeros(3),
    direction=np.array([1.0, 0.0, 0.0]),
    normal=np.array([0.0, 0.0, 1.0]),
    base_length=1,
    base_radius=0.1,
    density=2e3,
    youngs_modulus=1e4,
    shear_modulus=1e4,
)
restart_test_simulator.append(cosserat_rod)

# Left end fixed
restart_test_simulator.constrain(cosserat_rod).using(
    ea.OneEndFixedBC,
    constrained_position_idx=(0,),
    constrained_director_idx=(0,),
)

# Gravity
restart_test_simulator.add_forcing_to(cosserat_rod).using(
    ea.GravityForces, acc_gravity=np.array([10, 0, 0])
)

# The dissipation here leads to the difference
rod_dt = np.float64(1e-4)
damping_constant = np.float64(1e-3)
restart_test_simulator.dampen(cosserat_rod).using(
    ea.AnalyticalLinearDamper,
    damping_constant=damping_constant,
    time_step=rod_dt,
)
restart_test_simulator.finalize()

timestepper = ea.PositionVerlet()
curr_time = ea.integrate(
    stepper=timestepper, 
    systems=restart_test_simulator, 
    final_time=steps * rod_dt,
    n_steps=steps,
    progress_bar=False
)

velocity_copy = cosserat_rod.velocity_collection.copy()
ea.save_state(restart_test_simulator, "restart_data", float(curr_time))
ea.load_state(restart_test_simulator, "restart_data", False)
print(np.amax(np.abs(cosserat_rod.velocity_collection - velocity_copy)))

9.99949900659658e-08 # This should be zero

Possible cause
The load_state function invokes constraint_values and constrain_rates, which includes dissipation.