Skip to content

Orientation time integration algorithms and NewtonIntegrator.densityScaling

I just happened to see (in current master) some suspicious duplicated applications of the (local) density scaling idea:

if (densityScaling) angular_velocity *= state->densityScaling;

in orientation time integration methods introduced (in the present code shape for what concerns the present discussion) in 58738755:

  • NewtonIntegrator::leapfrogAsphericalRotateOmelyan_1998 with here being maybe duplicated there

  • NewtonIntegrator::leapfrogAsphericalRotateCarlos_2023: with here duplicated there

While I'm really not familiar with this NewtonIntegrator.densityScaling idea (and the underlying modifications of State.densityScaling by the GlobalStiffnessTimeStepper) I strongly suspect an error, especially in the latter method that implements a incremental update of angular velocity (as per Eq. (9) of [delValle2024] SPIRAL paper).

Where, because of the incremental update, the application of the densityScaling multiplication shall not apply this way, I presume.

For the Omelyan method, I do not really know at the moment.

PS: And I would tend to think that the situation in NewtonIntegrator::leapfrogAsphericalRotate with two direct computations of old and new values of angular velocities, instead of an incremental update, much better accommodates the two multiplications with state->densityScaling on that line and that other one