Merge pull request #1772 from knutfrode/dev

knutfrode · web-flow · commit 33b871f0e967 · 2026-02-15T22:34:09.000+01:00
Major change: horizontal_diffusivity is not anymore a config setting …
diff --git a/examples/example_biodegradation.py b/examples/example_biodegradation.py
@@ -15,7 +15,8 @@
 # No motion is needed for this test
 o.set_config('environment:constant', {k: 0 for k in
              ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']})
-o.set_config('drift', {'current_uncertainty': 0, 'wind_uncertainty': 0, 'horizontal_diffusivity': 10})
+o.set_config('environment:constant:horizontal_diffusivity', 10)
+o.set_config('drift', {'current_uncertainty': 0, 'wind_uncertainty': 0})
 
 #%%
 # Seeding some particles
diff --git a/examples/example_depth.py b/examples/example_depth.py
@@ -17,7 +17,7 @@
 
 #%%
 # Adding some diffusion
-o.set_config('drift:horizontal_diffusivity', 10)  # m2/s
+o.set_config('environment:constant:horizontal_diffusivity', 10)  # m2/s
 
 #%%
 # Seed 1000 elements at random depths
diff --git a/examples/example_horizontal_diffusion.py b/examples/example_horizontal_diffusion.py
@@ -40,7 +40,7 @@
 o2.add_reader([reader_norkyst, reader_arome])
 #o2.set_config('drift:current_uncertainty', .2) # Difference from first run
 #o2.set_config('drift:wind_uncertainty', 1)     # Difference from first run
-o2.set_config('drift:horizontal_diffusivity', 10)     # Difference from first run
+o2.set_config('environment:constant:horizontal_diffusivity', 10)     # Difference from first run
 o2.seed_elements(lon, lat, radius=500, number=2000, time=time)
 o2.run(duration=timedelta(hours=24))
 
@@ -50,7 +50,7 @@
 o3.add_reader([reader_norkyst, reader_arome])
 #o3.set_config('drift:current_uncertainty', .2) # Difference from first run
 #o3.set_config('drift:wind_uncertainty', 1)     # Difference from first run
-o3.set_config('drift:horizontal_diffusivity', 10)     # Difference from first run
+o3.set_config('environment:constant:horizontal_diffusivity', 10)     # Difference from first run
 o3.seed_elements(lon, lat, radius=500, number=2000, time=time)
 o3.run(duration=timedelta(hours=24), time_step=300, time_step_output=3600)
 
diff --git a/examples/example_huge_output.py b/examples/example_huge_output.py
@@ -13,7 +13,7 @@
 #%%
 # First make a simulation with two seedings, marked by *origin_marker*
 o = OceanDrift(loglevel=10)
-o.set_config('drift:horizontal_diffusivity', 10)
+o.set_config('environment:constant:horizontal_diffusivity', 10)
 t1 = datetime.now()
 t2 = t1 + timedelta(hours=6)
 number = 10000
diff --git a/examples/example_lake.py b/examples/example_lake.py
@@ -24,7 +24,7 @@
 o.add_reader(reader_lakes)
 o.set_config('general:use_auto_landmask', False)  # To use custom landmask instead
 o.set_config('environment:constant:x_sea_water_velocity', 1)
-o.set_config('drift:horizontal_diffusivity', 100)
+o.set_config('environment:constant:horizontal_diffusivity', 100)
 o.seed_elements(lon=48.819, lat=44.959, radius=5000, number=100, time=datetime.now())
 
 o.run(steps=10)
diff --git a/examples/example_long_multiprocessing.py b/examples/example_long_multiprocessing.py
@@ -46,7 +46,7 @@ def RunOceanDrift(pool_number):
     reader_topaz4 = reader_netCDF_CF_generic.Reader("https://thredds.met.no/thredds/dodsC/topaz/dataset-topaz4-arc-myoceanv2-be")
     o.add_reader(reader_topaz4, variables=['x_sea_water_velocity', 'y_sea_water_velocity','sea_water_temperature','sea_water_salinity','sea_floor_depth_below_sea_level'])
 
-    o.set_config('drift:horizontal_diffusivity', 50)
+    o.set_config('environment:constant:horizontal_diffusivity', 50)
 
     time = datetime(2023,1,1)
 
diff --git a/examples/example_long_radionuclides.py b/examples/example_long_radionuclides.py
@@ -27,7 +27,7 @@
 o.set_config('vertical_mixing:diffusivitymodel','environment')  # apply vertical diffusivity from ocean model
 # Vertical mixing requires fast time step
 o.set_config('vertical_mixing:timestep', 600.) # seconds
-o.set_config('drift:horizontal_diffusivity', 10)
+o.set_config('environment:constant:horizontal_diffusivity', 10)
 
 #%%
 o.set_config('radionuclide:particle_diameter',5.e-6)  # m
diff --git a/examples/example_oil_ice.py b/examples/example_oil_ice.py
@@ -20,7 +20,7 @@
 o.set_config('processes:dispersion',  False)
 o.set_config('processes:evaporation',  False)
 o.set_config('processes:emulsification',  False)
-o.set_config('drift:horizontal_diffusivity', 10)
+o.set_config('environment:constant:horizontal_diffusivity', 10)
 o.set_config('drift:truncate_ocean_model_below_m', 3)
 
 #%%
diff --git a/examples/example_oilspill_seafloor_biodegradation.py b/examples/example_oilspill_seafloor_biodegradation.py
@@ -18,7 +18,8 @@
     'https://tds.hycom.org/thredds/dodsC/FMRC_ESPC-D-V02_uv3z/FMRC_ESPC-D-V02_uv3z_best.ncd',
     'https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/ncep_global/NCEP_Global_Atmospheric_Model_best.ncd'])
 o.set_config('environment:constant:ocean_mixed_layer_thickness', 20)
-o.set_config('drift', {'current_uncertainty': 0, 'wind_uncertainty': 0, 'horizontal_diffusivity': 20})
+o.set_config('environment:constant:horizontal_diffusivity', 20)
+o.set_config('drift', {'current_uncertainty': 0, 'wind_uncertainty': 0})
 
 #%%
 # Configuration
diff --git a/examples/example_openberg.py b/examples/example_openberg.py
@@ -31,7 +31,7 @@
 # Simulating drift for 48 hours
 o = OpenBerg()
 o.set_config('drift:vertical_profile', False)
-o.set_config('drift:horizontal_diffusivity', 100)
+o.set_config('environment:constant:horizontal_diffusivity', 100)
 o.add_readers_from_list(forcing)
 o.seed_elements(**icebergs)
 o.run(duration=timedelta(hours=48))
diff --git a/examples/example_river_runoff.py b/examples/example_river_runoff.py
@@ -20,7 +20,7 @@
 #%%
 # First make a simulation with two seedings, marked by *origin_marker*
 o = OceanDrift(loglevel=20)
-o.set_config('drift:horizontal_diffusivity', 300)
+o.set_config('environment:constant:horizontal_diffusivity', 300)
 o.set_config('general:coastline_action', 'previous')
 o.set_config('general:coastline_approximation_precision', None)  # Faster to go back to previous
 t1 = datetime.now()
diff --git a/examples/example_seafloor_interaction.py b/examples/example_seafloor_interaction.py
@@ -37,7 +37,7 @@
 
     o.add_reader([reader_osc, reader_norkyst])
 
-    o.set_config('drift:horizontal_diffusivity', 0)
+    o.set_config('environment:constant:horizontal_diffusivity', 0)
     o.set_config('environment:constant:y_sea_water_velocity', 0)
     o.set_config('environment:constant:land_binary_mask', 0)
     o.set_config('general:use_auto_landmask', False)
diff --git a/examples/example_wind_drift_coefficient_from_trajectory.py b/examples/example_wind_drift_coefficient_from_trajectory.py
@@ -30,7 +30,7 @@
 
 #%%
 # Adding some horizontal diffusivity as "noise"
-ot.set_config('drift:horizontal_diffusivity', 10)
+ot.set_config('environment:constant:horizontal_diffusivity', 10)
 
 #%%
 # Using a wind_drift_factor of 0.33 i.e. drift is current + 3.3% of wind speed
diff --git a/opendrift/models/basemodel/__init__.py b/opendrift/models/basemodel/__init__.py
@@ -687,7 +687,8 @@ def interact_with_coastline(self, final=False):
                                      self.time,
                                      self.elements.lon,
                                      self.elements.lat,
-                                     self.elements.z)
+                                     self.elements.z,
+                                     element_ID = self.elements.ID)
             self.environment.land_binary_mask = en.land_binary_mask
 
         if i == 'stranding':  # Deactivate elements on land, but not in air
@@ -977,7 +978,8 @@ def closest_ocean_points(self, lon, lat):
                                         lon=lon,
                                         lat=lat,
                                         z=0 * lon,
-                                        time=land_reader.start_time)[0]['land_binary_mask']
+                                        time=land_reader.start_time,
+                                        element_ID=self.elements.ID)[0]['land_binary_mask']
             if land.max() == 0:
                 logger.info('All points are in ocean')
                 return lon, lat, None
@@ -1003,7 +1005,8 @@ def closest_ocean_points(self, lon, lat):
                                             lon=longrid,
                                             lat=latgrid,
                                             z=0 * longrid,
-                                            time=land_reader.start_time)[0]['land_binary_mask']
+                                            time=land_reader.start_time,
+                                            element_ID=self.elements.ID)[0]['land_binary_mask']
             if landgrid.size == 0:
                 # Need to catch this before trying .min() on it...
                 logger.warning('Land grid has zero size, cannot move elements.')
@@ -1660,8 +1663,13 @@ def seed_from_ladim(self, ladimfile, roms):
 
     def horizontal_diffusion(self):
         """Move elements with random walk according to given horizontal diffuivity."""
-        D = self.get_config('drift:horizontal_diffusivity')
-        if D == 0:
+        if 'horizontal_diffusivity' not in self.required_variables:
+            logger.debug('No horizontal diffusion')
+            return
+        D = self.environment.horizontal_diffusivity
+        if len(D) == 0:
+            return  # no active elements
+        if D.max() == 0:
             logger.debug('Horizontal diffusivity is 0, no random walk.')
             return
         if self.num_elements_active() == 0:
@@ -1673,9 +1681,12 @@ def horizontal_diffusion(self):
         y_vel = self.elements.moving * np.sqrt(2 * D / dt) * np.random.normal(
             scale=1, size=self.num_elements_active())
         speed = np.sqrt(x_vel * x_vel + y_vel * y_vel)
+        if D.min() == D.max():
+            dstring = f'{D.min()}'
+        else:
+            dstring = f'{D.min()} - {D.max()}'
         logger.debug(
-            'Moving elements according to horizontal diffusivity of %s, with speeds between %s and %s m/s'
-            % (D, speed.min(), speed.max()))
+            f'Moving elements according to horizontal diffusivity of {dstring}, with speeds between {speed.min()} and {speed.max()} m/s')
         self.update_positions(x_vel, y_vel)
 
     def deactivate_elements(self, indices, reason='deactivated'):
@@ -2147,7 +2158,8 @@ def run(self,
                                          self.elements.lat,
                                          self.elements.z,
                                          self.required_profiles,
-                                         self.profiles_depth)
+                                         self.profiles_depth,
+                                         element_ID = self.elements.ID)
 
                 self.calculate_missing_environment_variables()
 
diff --git a/opendrift/models/basemodel/environment.py b/opendrift/models/basemodel/environment.py
@@ -504,7 +504,8 @@ def get_environment(self,
                         lat,
                         z,
                         profiles=None,
-                        profiles_depth=None):
+                        profiles_depth=None,
+                        element_ID=None):
         '''Retrieve environmental variables at requested positions.
 
         Args:
@@ -523,6 +524,8 @@ def get_environment(self,
 
             profiles_depth: depth of profiles in meters, as a positive number
 
+            element_ID: ID of actual elements, for the readers that might provide data based on this ID (exception)
+
         Updates:
             Buffer (raw data blocks) for each reader stored for performance:
                 [readers].var_block_before (last before requested time)
@@ -613,9 +616,12 @@ def get_environment(self,
                             )
                             return self.get_environment(variables, time, lon,
                                                         lat, z, profiles,
-                                                        profiles_depth)
+                                                        profiles_depth, element_ID=element_ID)
                     continue
                 # Fetch given variables at given positions from current reader
+                if reader._element_ID is not None:
+                    logger.debug(f'Setting _element_ID for reader {reader_name}')
+                    reader._element_ID = element_ID
                 try:
                     logger.debug('Data needed for %i elements' %
                                  len(missing_indices))
@@ -632,7 +638,8 @@ def get_environment(self,
                             variable_group, profiles = profiles_from_reader,
                             profiles_depth = profiles_depth, time = time,
                             lon=lon[missing_indices], lat=lat[missing_indices],
-                            z=z[missing_indices], rotate_to_proj=self.proj_latlon)
+                            z=z[missing_indices], rotate_to_proj=self.proj_latlon,
+                            element_ID=element_ID)
 
                 except NotCoveredError as e:
                     logger.info(e)
@@ -645,7 +652,7 @@ def get_environment(self,
                             )
                             return self.get_environment(variables, time, lon,
                                                         lat, z, profiles,
-                                                        profiles_depth)
+                                                        profiles_depth, element_ID=element_ID)
                     continue
 
                 except Exception as e:  # Unknown error
@@ -677,7 +684,7 @@ def get_environment(self,
                             )
                             return self.get_environment(variables, time, lon,
                                                         lat, z, profiles,
-                                                        profiles_depth)
+                                                        profiles_depth, element_ID=element_ID)
                     continue
 
                 # Copy retrieved variables to env array, and mask nan-values
@@ -764,7 +771,7 @@ def get_environment(self,
                                 'Missing variables: calling get_environment recursively'
                             )
                             return self.get_environment(variables, time, lon,
-                                                        lat, z, profiles)
+                                                        lat, z, profiles, profiles_depth, element_ID=element_ID)
 
         logger.debug('---------------------------------------')
         logger.debug('Finished processing all variable groups')
diff --git a/opendrift/models/chemicaldrift.py b/opendrift/models/chemicaldrift.py
@@ -99,6 +99,7 @@ class ChemicalDrift(OceanDrift):
         'land_binary_mask': {'fallback': None},
         'sea_floor_depth_below_sea_level': {'fallback': 10000},
         'ocean_vertical_diffusivity': {'fallback': 0.0001, 'profiles': True},
+        'horizontal_diffusivity': {'fallback': 0, 'important': False},
         'sea_water_temperature': {'fallback': 10, 'profiles': True},
         'sea_water_salinity': {'fallback': 34, 'profiles': True},
         'upward_sea_water_velocity': {'fallback': 0},
diff --git a/opendrift/models/oceandrift.py b/opendrift/models/oceandrift.py
@@ -79,6 +79,7 @@ class OceanDrift(OpenDriftSimulation):
         'ocean_vertical_diffusivity': {'fallback': 0,
              'skip_if': ['drift:vertical_mixing', 'is', False],
              'profiles': True},
+        'horizontal_diffusivity': {'fallback': 0},
         'sea_surface_wave_significant_height': {'fallback': 0},
         'sea_surface_wave_stokes_drift_x_velocity': {'fallback': 0,
             'skip_if': ['drift:stokes_drift', 'is', False]},
diff --git a/opendrift/models/openberg.py b/opendrift/models/openberg.py
@@ -303,6 +303,10 @@ class OpenBerg(OpenDriftSimulation):
         "sea_surface_y_slope": {"fallback": 0, 'important': False},
         "x_wind": {"fallback": None, "important": True},
         "y_wind": {"fallback": None, "important": True},
+        # Since OpenBerg model is deterministic for given iceberg size,
+        # (in contrast to the Leeway model), we use a default diffusivity
+        # to yield some variability.
+        "horizontal_diffusivity": {"fallback": 0, "important": False},
         "sea_surface_wave_significant_height": {"fallback": 0},
         "sea_surface_wave_from_direction": {"fallback": 0},
         "sea_surface_wave_stokes_drift_x_velocity": {"fallback": 0, 'important': False},
@@ -419,11 +423,6 @@ def __init__(self, *args, **kwargs):
             },  
         })
 
-        # Since OpenBerg model is deterministic for given iceberg size,
-        # (in contrast to the Leeway model), we use a default diffusivity
-        # to yield some variability.
-        self._set_config_default('drift:horizontal_diffusivity', 100)
-
 
     def advect_iceberg(self):
         sail = self.elements.sail
diff --git a/opendrift/models/openoil/openoil.py b/opendrift/models/openoil/openoil.py
@@ -282,6 +282,9 @@ class OpenOil(OceanDrift):
         'sea_floor_depth_below_sea_level': {
             'fallback': 10000
         },
+        'horizontal_diffusivity': {
+            'fallback': 0, 'important': False,
+        },
         'ocean_vertical_diffusivity': {
             'fallback': 0.02,
             'important': False,
diff --git a/opendrift/models/physics_methods.py b/opendrift/models/physics_methods.py
@@ -597,7 +597,7 @@ def advect_ocean_current(self, factor=1):
             mid_env, profiles, missing = self.env.get_environment(
                 ['x_sea_water_velocity', 'y_sea_water_velocity'],
                 self.time + self.time_step/2,
-                mid_lon, mid_lat, self.elements.z, profiles=None)
+                mid_lon, mid_lat, self.elements.z, profiles=None, element_ID=self.elements.ID)
             if self.get_config('drift:advection_scheme') == 'runge-kutta4':
                 logger.debug('Runge-kutta 4th order...')
                 x_vel2 = mid_env['x_sea_water_velocity']
@@ -612,7 +612,7 @@ def advect_ocean_current(self, factor=1):
                 env2, profiles, missing = self.env.get_environment(
                     ['x_sea_water_velocity', 'y_sea_water_velocity'],
                     self.time + self.time_step/2,
-                    lon2, lat2, self.elements.z, profiles=None)
+                    lon2, lat2, self.elements.z, profiles=None, element_ID=self.elements.ID)
                 # Third step
                 x_vel3 = env2['x_sea_water_velocity']
                 y_vel3 = env2['y_sea_water_velocity']
@@ -626,7 +626,7 @@ def advect_ocean_current(self, factor=1):
                 env3, profiles, missing = self.env.get_environment(
                     ['x_sea_water_velocity', 'y_sea_water_velocity'],
                     self.time + self.time_step,
-                    lon3, lat3, self.elements.z, profiles=None)
+                    lon3, lat3, self.elements.z, profiles=None, element_ID=self.elements.ID)
                 # Fourth step
                 x_vel4 = env3['x_sea_water_velocity']
                 y_vel4 = env3['y_sea_water_velocity']
diff --git a/opendrift/models/shipdrift.py b/opendrift/models/shipdrift.py
@@ -92,6 +92,10 @@ class ShipDrift(OpenDriftSimulation):
         'land_binary_mask': {'fallback': None},
         'x_sea_water_velocity': {'fallback': None},
         'y_sea_water_velocity': {'fallback': None},
+        # Since the ShipDrift model is deterministic for given ship size,
+        # (in contrast to the Leeway model), we use a default diffusivity
+        # to yield some variability.
+        'horizontal_diffusivity': {'fallback': 100, 'important': False},
         'sea_surface_wave_stokes_drift_x_velocity': {'fallback': 0},
         'sea_surface_wave_stokes_drift_y_velocity': {'fallback': 0},
         'sea_surface_wave_significant_height': {'fallback': 0},
@@ -148,11 +152,6 @@ def __init__(self, *args, **kwargs):
             'description': 'If ships are oriented to the left or right of the downwind direction,'
                 'or whether this is unknown. Left/right means that wind will hit ship from backboard/steerboard'}})
 
-        # Since the ShipDrift model is deterministic for given ship size,
-        # (in contrast to the Leeway model), we use a default diffusivity
-        # to yield some variability.
-        self._set_config_default('drift:horizontal_diffusivity', 100)
-
         self._set_config_default('drift:max_speed', 2)
 
     def seed_elements(self, *args, **kwargs):
diff --git a/opendrift/readers/basereader/__init__.py b/opendrift/readers/basereader/__init__.py
@@ -118,6 +118,8 @@ def __init__(self):
 
         self.is_lazy = False  # Generally False
 
+        self._element_ID = None  # Generally readers should not know ID of elements
+
         # Set projection for coordinate transformations
         if self.proj is None:
             if self.proj4 is not None:
diff --git a/opendrift/readers/basereader/variables.py b/opendrift/readers/basereader/variables.py
diff --git a/opendrift/readers/reader_constant.py b/opendrift/readers/reader_constant.py
diff --git a/tests/models/test_basemodel.py b/tests/models/test_basemodel.py
diff --git a/tests/models/test_chemicaldrift.py b/tests/models/test_chemicaldrift.py
diff --git a/tests/models/test_models.py b/tests/models/test_models.py
diff --git a/tests/readers/test_xarray.py b/tests/readers/test_xarray.py
