Initial values of sampled fields are not in output file

[erikvansebille]

Parcels version

v4-alpha

Description

Since #2568 which moved the initial pfile.write() to before the kernels are executed, we are back at the problem that particle Variables that are used to sample fields start with the initial value of these particles (which is generally not the field value); see also example below

Note that this was one of the main reasons to change the kernel loop in v3 (#1402) so that the writing to the particle file is only done after the kernels are executed. But that created other problems that now seem fixed in #2568.

In any case, we need to come up with a solution to the issue that initial values are not written. Options could be

1. again rethink the kernel loop and consider moving the writing until after the kernels have been executed
2. find a way to automatically execute all the sampling kernels before going into the for-loop in the pset.execute()
3. teach users that they first need to execute their sampling kernels before they run a pset.execute (this is how it was like in v2)

Option 3 is the easiest to implement, but will require most effort from users. I have no idea yet how to implement option 2.

Code sample

import numpy as np

import parcels
import parcels.tutorial

# Load the CopernicusMarine data in the Agulhas region from the example_datasets
ds_fields = parcels.tutorial.open_dataset(
    "CopernicusMarine_data_for_Argo_tutorial/data"
)
ds_fields.load()  # load the dataset into memory

# Convert to SGRID-compliant dataset and create FieldSet
fields = {"thetao": ds_fields["thetao"]}
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)

SampleParticle = parcels.Particle.add_variable(
    parcels.Variable("temperature", dtype=np.float32, initial=np.nan)
)

def SampleT(particles, fieldset):
    particles.temperature = fieldset.thetao[
        particles.time, particles.z, particles.lat, particles.lon
    ]

pset = parcels.ParticleSet(
    fieldset=fieldset, pclass=SampleParticle, lon=[32, 33], lat=[-31.5, -30.5]
)

output_file = parcels.ParticleFile("sampletemp.parquet", outputdt=np.timedelta64(1, "h"))

pset.execute(
    SampleT,
    runtime=np.timedelta64(4, "h"),
    dt=np.timedelta64(30, "m"),
    output_file=output_file,
)

df = parcels.read_particlefile(output_file._path)
df

Output (note the NaNs on the first two rows):

[erikvansebille]

I now also added a (failing) unit test in #2622