Add test case for multi-model advection

[VeckoTheGecko]

In exploring #2586, I am interested to see how we handle combining multiple models together. In this case, I'm interested in the time handling with different time origins.

I want to create a test case which explores this so that we can add it to test_advection.py

I think the current abstraction of add_field is insufficient to handle this. How was this done in v3 @erikvansebille ? fset1 + fset2?

[VeckoTheGecko]

I have some ideas on how to proceed with this from a design POV - but I first want to get an idea of how it was done in v3 so that I can flesh out my ideas

[erikvansebille]

In v3, you can simply add Fields with different grids in one FieldSet construction, see e.g. the unit test at

Parcels/tests/test_grids.py

Lines 31 to 79 in e6b6c70

	def test_multi_structured_grids(mode):
	def temp_func(lon, lat):
	return 20 + lat / 1000.0 + 2 * np.sin(lon * 2 * np.pi / 5000.0)
	a = 10000
	b = 10000
	# Grid 0
	xdim_g0 = 201
	ydim_g0 = 201
	# Coordinates of the test fieldset (on A-grid in deg)
	lon_g0 = np.linspace(0, a, xdim_g0, dtype=np.float32)
	lat_g0 = np.linspace(0, b, ydim_g0, dtype=np.float32)
	time_g0 = np.linspace(0.0, 1000.0, 2, dtype=np.float64)
	grid_0 = RectilinearZGrid(lon_g0, lat_g0, time=time_g0)
	# Grid 1
	xdim_g1 = 51
	ydim_g1 = 51
	# Coordinates of the test fieldset (on A-grid in deg)
	lon_g1 = np.linspace(0, a, xdim_g1, dtype=np.float32)
	lat_g1 = np.linspace(0, b, ydim_g1, dtype=np.float32)
	time_g1 = np.linspace(0.0, 1000.0, 2, dtype=np.float64)
	grid_1 = RectilinearZGrid(lon_g1, lat_g1, time=time_g1)
	u_data = np.ones((lon_g0.size, lat_g0.size, time_g0.size), dtype=np.float32)
	u_data = 2 * u_data
	u_field = Field("U", u_data, grid=grid_0, transpose=True)
	temp0_data = np.empty((lon_g0.size, lat_g0.size, time_g0.size), dtype=np.float32)
	for i in range(lon_g0.size):
	for j in range(lat_g0.size):
	temp0_data[i, j, :] = temp_func(lon_g0[i], lat_g0[j])
	temp0_field = Field("temp0", temp0_data, grid=grid_0, transpose=True)
	v_data = np.zeros((lon_g1.size, lat_g1.size, time_g1.size), dtype=np.float32)
	v_field = Field("V", v_data, grid=grid_1, transpose=True)
	temp1_data = np.empty((lon_g1.size, lat_g1.size, time_g1.size), dtype=np.float32)
	for i in range(lon_g1.size):
	for j in range(lat_g1.size):
	temp1_data[i, j, :] = temp_func(lon_g1[i], lat_g1[j])
	temp1_field = Field("temp1", temp1_data, grid=grid_1, transpose=True)
	other_fields = {}
	other_fields["temp0"] = temp0_field
	other_fields["temp1"] = temp1_field
	fieldset = FieldSet(u_field, v_field, fields=other_fields)

An alternative is the [FieldSet.add_field()] method.

Parcels/parcels/fieldset.py

Lines 175 to 212 in e6b6c70

	def add_field(self, field: Field | NestedField, name: str | None = None):
	"""Add a :class:`parcels.field.Field` object to the FieldSet.
	Parameters
	----------
	field : parcels.field.Field
	Field object to be added
	name : str
	Name of the :class:`parcels.field.Field` object to be added. Defaults
	to name in Field object.
	Examples
	--------
	For usage examples see the following tutorials:
	* `Nested Fields <../examples/tutorial_NestedFields.ipynb>`__
	* `Unit converters <../examples/tutorial_unitconverters.ipynb>`__ (Default value = None)
	"""
	if self._completed:
	raise RuntimeError(
	"FieldSet has already been completed. Are you trying to add a Field after you've created the ParticleSet?"
	)
	name = field.name if name is None else name
	if hasattr(self, name): # check if Field with same name already exists when adding new Field
	raise RuntimeError(f"FieldSet already has a Field with name '{name}'")
	if isinstance(field, NestedField):
	setattr(self, name, field)
	for fld in field:
	self.gridset.add_grid(fld)
	fld.fieldset = self
	else:
	setattr(self, name, field)
	self.gridset.add_grid(field)
	field.fieldset = self

An (older) example of this in the wild is https://github.com/Parcels-code/SKIM_dFADtracking/blob/1d028d8ef78ed3d3f0b6d47695e274d672eacd79/fadset_run.py#L138-L140 (note, that was using parcels v2)

Note that in v4, we already have an example that requires combining Fields on different grids into one FieldSet: the Nested Grids tutorial. There we do

for i, ds in enumerate(ds_in):
    fset = parcels.FieldSet.from_sgrid_conventions(ds, mesh="spherical")

    for fld in fset.fields.values():
        fld.name = f"{fld.name}{i}"
        fields.append(fld)
fieldset = parcels.FieldSet(fields)

Does this help?

[VeckoTheGecko]

This is insightful, but I think that this highlights that we need a clearer way for people to add fieldsets together.

Something like:

fset1 = FieldSet.from_copernicusmarine(...)
fset2 = FieldSet.from_icon(...)
fset = fset1 + fset2

Since at the moment it requires a lot of boilerplate.

Perhaps this is worth discussing in person at some time (to flesh out how grids, and fields work when combining multiple models).

[erikvansebille]

Yes, I like that syntax of fset = fset1 + fset2.

Could/should we also incorporate fset = fset + field? and fset += field? As you see in the examples above, This is what the add_field() in v3 was used for. Or would that just complicate the API?