Time step for symba7 integration
This is a spin off oh #523 (closed)
The question is "how do we "synchornize" the symba7 time step and the fargo DT.
@alainm asked:
@elena @gpichier if I understand correctly, the last thing to do to complete the implementation s to deal with the dime step issue.
Right now, it is decided and continuously updated in real Disk::conditionCFL(real dtLeft)
and is based on the state on the gas at a given time. Gariele mentioned that the Symba7 time step has to be constant over the whole integration but am not sure if that means de whole run (which seems strange: even if symba7 has an internal state, can't we change over a restart ?)
If yes, does it mean we should choose a dt that would be a whole fraction of an evaluation of the expected minimum of conditionCFL over the whole simulation ?
@elena proposed:
I would do the following:
- Set the Symba_dt == disk.physic().timeStep
- divide a timestep in nt equals intervals of lenght dt:
real dtemp = 0.;
real dtLeft = disk.physic().timeStep-dtemp;
real dt = disk.conditionCFL(dtLeft);
int nt = ceil(dtLeft/dt);
- do a symba7 step of size Symba_dt. Keep memory of the planets position and velocities at t-Symba_dt in order to do a keplerian interpolation of the orbits. This allow to have nt values of planets positions and velocities to be used to update the planetary system at each step dt
- do a loop on nt gas steps and update the planetaty spositions and velocities from the interpolated values. Compute the forces from the disk onto the planets at each dt and sum them on a cumulated force
- at time t+Symba_dt do a dissipative step for planet velocities croviding the cumulated forces from the disk
This implies a change in the simulation.cpp: instead of continuously updating dt we use the same dt and loop over nt steps to complete a disk.physic().timestep
The point is that we need to advance the planet(s) at each dt and this is reasonably done via a Keplerian interpolation of the planet orbit over disk.physic().timestep
There is may be an alternative way to have positions and velocities of the planet at each dt , this would require less changes in the code
I did mention all this staff in https://gitlab.oca.eu/DISC/swift/-/issues/1 but we did not continue the discussion, but now it is time to do it