Merge pull request #1298 from knutfrode/dev

Leeway capsizing probability now using tanh probability. New method s…

environment.yml

@@ -30,6 +30,6 @@ dependencies:
   - utm
   - roaring-landmask>=0.7
   - trajan>=0.1.3
-  - adios_db
+  - adios_db<1.2
   - copernicusmarine
   - pip

examples/example_leeway.py

@@ -32,12 +32,15 @@
 o.set_config('environment:fallback:y_sea_water_velocity', 0)
 
 #%%
-# Activating capsizing for high winds
+# Activating capsizing for high winds, with probability per hour given by
+# p(windspeed) = 0.5 + 0.5*tanh((windspeed-wind_threshold)/sigma)
 o.set_config('capsizing', True)
-o.set_config('capsizing:wind_threshold', 15)  # Capsizing can take place for winds above 15 m/s
-o.set_config('capsizing:probability_per_hour', .01)  # 1% probability of capsizing per hour for wind above threshold
+o.set_config('capsizing:wind_threshold', 30)
+o.set_config('capsizing:wind_threshold_sigma', 5)
 o.set_config('capsizing:leeway_fraction', 0.4)  # Reducing leeway coefficients to 40% of original after capsize
 
+o.plot_capsize_probability()
+
 #%%
 # Seed leeway elements at defined position and time
 object_type = 26  # 26 = Life-raft, no ballast

opendrift/models/basemodel/__init__.py

@@ -4287,6 +4287,13 @@ def get_time_array(self):
         time_array_relative = [td * i for i in range(self.steps_output)]
         return time_array, time_array_relative
 
+    def simulation_direction(self):
+        """Return 1 for a forward simulation, and -1 for a backward simulation"""
+        if self.time_step.days < 0:
+            return -1
+        else:
+            return 1
+
     @require_mode(mode=Mode.Result)
     def plot_environment(self, filename=None, ax=None, show=True):
         """Plot mean wind and current velocities of element of last run."""

opendrift/models/leeway.py

@@ -253,21 +253,21 @@ def __init__(self, d=None, *args, **kwargs):
             },
               'capsizing:wind_threshold': {
                 'type': 'float',
-                'default': 15,
+                'default': 30,
                 'min': 0,
                 'max': 50,
                 'description':
-                'Probability of capsizing is activated when wind exceeds this value (m/s)',
+                'Probability of capsizing per hour is: 0.5 + 0.5tanh((windspeed-wind_threshold)/wind_threshold_sigma)',
                 'units': 'm/s',
                 'level': CONFIG_LEVEL_BASIC
             },
-               'capsizing:probability_per_hour': {
+               'capsizing:wind_threshold_sigma': {
                 'type': 'float',
-                'default': 0.01,
+                'default': 5,
                 'min': 0,
-                'max': 1,
+                'max': 20,
                 'description':
-                'The probability of capsizing per hour, when wind exceeds value given by config item capsize:wind_threshold',
+                'Sigma parameter in parameterization of capsize probability',
                 'units': 'm/s',
                 'level': CONFIG_LEVEL_BASIC
             },
@@ -400,12 +400,20 @@ def list_object_categories(self, substr=None):
                     continue
             print('%i %s %s' % (i + 1, objkey, description))
 
-    def prepare_run(self):
-        '''Not allowing capsizing for backwards runs'''
-        if self.time_step.days < 0 and self.get_config('capsizing') is True:
-            raise ValueError('Capsizing is not allowed for backwards runs')
-
-        super(Leeway, self).prepare_run()
+    def plot_capsize_probability(self):
+        U = np.linspace(0, 35, 100)
+        wind_threshold = self.get_config('capsizing:wind_threshold')
+        sigma = self.get_config('capsizing:wind_threshold_sigma')
+        p = self.capsize_probability(U, wind_threshold, sigma)
+        import matplotlib.pyplot as plt
+        plt.plot(U, p)
+        plt.title(f'p(u) = 0.5 + 0.5*tanh((u - {wind_threshold} / {sigma})')
+        plt.xlabel('Wind speed  [m/s]')
+        plt.ylabel('Probability of capsizing per hour')
+        plt.show()
+
+    def capsize_probability(self, wind, threshold, sigma):
+        return .5 + .5*np.tanh((wind-threshold)/sigma)
 
     def update(self):
         """Update positions and properties of leeway particles."""
@@ -418,18 +426,20 @@ def update(self):
         # Capsizing
         if self.get_config('capsizing') is True:
             wind_threshold = self.get_config('capsizing:wind_threshold')
-            hw = np.where(np.logical_and(windspeed>wind_threshold, self.elements.capsized==0))[0]
-            capsizing = np.random.rand(len(hw))<self.get_config('capsizing:probability_per_hour')*\
-                                        self.time_step.total_seconds()/3600  # NB: assuming small timestep
-            capsizing = hw[capsizing]
-            logger.debug(f'Capsizing {len(capsizing)} of {len(hw)} elements where wind is above threshold of {wind_threshold} m/s')
-            # Reducing downwind and crosswind leeway slope and eps (not offset) by given factor
-            cf = self.get_config('capsizing:leeway_fraction')
-            self.elements.capsized[capsizing] = 1
-            self.elements.downwind_slope[capsizing] *= cf
-            self.elements.downwind_eps[capsizing] *= cf
-            self.elements.crosswind_slope[capsizing] *= cf
-            self.elements.crosswind_eps[capsizing] *= cf
+            wind_threshold_sigma = self.get_config('capsizing:wind_threshold_sigma')
+            # For forward run, elements can be capsized, but for backwards run, only capsized elements can be un-capsized
+            if self.simulation_direction() == 1:  # forward run
+                can_be_capsized = np.where(self.elements.capsized==0)[0]
+            else:
+                can_be_capsized = np.where(self.elements.capsized==1)[0]
+            if len(can_be_capsized) > 0:
+                probability = self.capsize_probability(windspeed[can_be_capsized],
+                        wind_threshold, wind_threshold_sigma)*np.abs(self.time_step.total_seconds())/3600
+                # NB: assuming small timestep
+                to_be_capsized = np.where(np.random.rand(len(can_be_capsized)) < probability)[0]
+                to_be_capsized = can_be_capsized[to_be_capsized]
+                logger.warning(f'Capsizing {len(to_be_capsized)} of {len(can_be_capsized)} elements')
+                self.elements.capsized[to_be_capsized] = 1 - self.elements.capsized[to_be_capsized]
 
         # Move particles with the leeway CCC TODO
         downwind_leeway = (
@@ -444,6 +454,9 @@ def update(self):
         costh = np.cos(winddir)
         y_leeway = downwind_leeway * costh + crosswind_leeway * sinth
         x_leeway = -downwind_leeway * sinth + crosswind_leeway * costh
+        capsize_fraction = self.get_config('capsizing:leeway_fraction')  # Reducing leeway for capsized elements
+        x_leeway[self.elements.capsized==1] *= capsize_fraction
+        y_leeway[self.elements.capsized==1] *= capsize_fraction
         self.update_positions(-x_leeway, y_leeway)
 
         # Move particles with ambient current

tests/models/test_leeway.py

@@ -19,7 +19,7 @@
 
 import os
 import time
-from datetime import datetime
+from datetime import datetime, timedelta
 from . import *
 
 from opendrift.readers import reader_global_landmask
@@ -77,3 +77,45 @@ def test_leewayrun(tmpdir, test_data):
             print(line)
     import filecmp
     assert filecmp.cmp(asciif, asciitarget)
+
+def test_capsize():
+    o = Leeway(loglevel=20)
+    o.set_config('environment:constant', {'x_sea_water_velocity': 0, 'y_sea_water_velocity': 0,
+                    'x_wind': 25, 'y_wind': 0, 'land_binary_mask': 0})
+    o.set_config('capsizing', True)
+    o.set_config('capsizing:wind_threshold', 30)
+    o.set_config('capsizing:wind_threshold_sigma', 3)
+    o.set_config('capsizing:leeway_fraction', .4)
+    o.seed_elements(lon=0, lat=60, time=datetime.now(), number=100)
+    o.run(time_step=900, time_step_output=900, duration=timedelta(hours=6))
+    assert o.elements.capsized.max() == 1
+    assert o.elements.capsized.min() == 0
+    assert o.elements.capsized.sum() == 18
+
+    # Backward run, checking that forward capsizing is not happening
+    ob = Leeway(loglevel=20)
+    ob.set_config('capsizing', True)
+    ob.set_config('capsizing:wind_threshold', 30)
+    ob.set_config('capsizing:wind_threshold_sigma', 3)
+    ob.set_config('capsizing:leeway_fraction', .4)
+    ob.set_config('environment:constant', {'x_sea_water_velocity': 0, 'y_sea_water_velocity': 0,
+                    'x_wind': 25, 'y_wind': 0, 'land_binary_mask': 0})
+    ob.seed_elements(lon=0, lat=60, time=datetime.now(), number=100)
+    ob.run(time_step=-900, time_step_output=900, duration=timedelta(hours=6))
+    assert ob.elements.capsized.max() == 0
+    assert ob.elements.capsized.min() == 0
+    assert ob.elements.capsized.sum() == 0
+
+    # Backward run, checking that backward capsizing does happen
+    ob = Leeway(loglevel=20)
+    ob.set_config('capsizing', True)
+    ob.set_config('capsizing:wind_threshold', 30)
+    ob.set_config('capsizing:wind_threshold_sigma', 3)
+    ob.set_config('capsizing:leeway_fraction', .4)
+    ob.set_config('environment:constant', {'x_sea_water_velocity': 0, 'y_sea_water_velocity': 0,
+                    'x_wind': 25, 'y_wind': 0, 'land_binary_mask': 0})
+    ob.seed_elements(lon=0, lat=60, time=datetime.now(), number=100, capsized=1)
+    ob.run(time_step=-900, time_step_output=900, duration=timedelta(hours=6))
+    assert ob.elements.capsized.max() == 1
+    assert ob.elements.capsized.min() == 0
+    assert ob.elements.capsized.sum() == 82