OpenDrift / opendrift

Open source framework for ocean trajectory modelling
https://opendrift.github.io
GNU General Public License v2.0
245 stars 120 forks source link

Issue with initial particle positions in backward mode #1010

Open simonweppe opened 1 year ago

simonweppe commented 1 year ago

Hi Opendrift dev team @knutfrode @gauteh , I think we came across a bug that affects the initial (i.e. at t0) and following particle positions in backward simulations.

We realized that by looking at the first timestep of our backward simulations (that use a one-off release) and noting that it wasn't the same as specified in seed_elements() (it looked like a step after the initial position specified). The problem doesn't seem to happen for forward simulations.

I prepared a set of simple simulations with constant currents so the issue can be reproduced. test_initial_position_backwards_runs.zip

Content is summarized below

I'll have a go at stepping through the code to see where the issue might come from but I think you'll likely be much faster than me to pinpoint where it might come from !

FORWARD

BACKWARD

BACKWARD_output_step_equals_simulation_step

simonweppe commented 1 year ago

I think what happens is that the first step of the history array (which should be the initial position), gets overwritten until the first output step because time_ind reamains equal to 0 until first output step

at that line : https://github.com/OpenDrift/opendrift/blob/master/opendrift/models/basemodel.py#L3056

knutfrode commented 1 year ago

Yes, this is probably related to this bug/shortcoming: https://github.com/OpenDrift/opendrift/discussions/854

So OpenDrift does presently not store the actual seeding locations, but rather element positions at time step t0, t1, t2, ... tn (for forwards and opposite for backwards run). Thus when seeding over time in a forward run, particles seeded between t0 and t1 are only recorded at t1, when they have moved a little. For backwards run, even instant seeds seem to have this problem.

In the near future, the actual seeding locations and times of each particle will be saved in the netCDF file as well. In the meantime, a workaround is to obtain the locations from the initial seed calls, as discussed in the above mentioned issue.

simonweppe commented 1 year ago

Thanks @knutfrode - copy that.

However it seems that, when time_step!=time_step_output, not only the initial particle positions are not saved to file (that we can workaround easily) but the following positions saved to file are not the same as when using time_step==time_step_output, when looking at the same times (likely shifted in time by that initial offset before initial output step). Final positions look consistent.

See illustration below. The red and magenta particle clouds at t1 and t2 are not at the same positions in the two scenarios ( time_step!=time_step_output and time_step==time_step_output), even though they are labelled with same timestamp in netcdf files. (it's the ones for the run time_step!=time_step_output that are wrong since current is constant in that simple run).

So it's still quite problematic if we need to looks at successive past positions (e.g. SAR or oil backtrack). Workaround is to save all timesteps (ok for small simulations) or use the opendrift native "object"..but definitely something to be aware of.

BACKWARD

BACKWARD_output_step_equals_simulation_step

knutfrode commented 1 year ago

Yes, you are right, there seem to be an independent error here. I guess this is related to what I had forgotten: when elements are supposed to be released between two actual release time steps, they are still advected with a full timestep, and not with the correct fraction which should be < 1. We are in the middle of a fairly large rewrite to OpenDrift 2.0, so will have a look at this issue after this.

Your above test scripts will be useful in creating a fast and minimalistic unit-test that reproduces this issue, and which can be used for test-driven development.

simonweppe commented 1 year ago

@knutfrode ok copy that. Looking forward to Opendrift 2.0 !