From bdab0f0b9df2c3a880a41f27395205de9e9edf94 Mon Sep 17 00:00:00 2001 From: Knut-Frode Dagestad Date: Tue, 14 May 2024 15:24:41 +0200 Subject: [PATCH 1/3] Leeway capsizing probability now using tanh probability. New method simulation_direction() is 1 for forward runs, and -1 for backwards runs --- examples/example_leeway.py | 9 ++-- opendrift/models/basemodel/__init__.py | 7 +++ opendrift/models/leeway.py | 61 ++++++++++++++++---------- tests/models/test_leeway.py | 38 +++++++++++++++- 4 files changed, 87 insertions(+), 28 deletions(-) diff --git a/examples/example_leeway.py b/examples/example_leeway.py index c547b1a3e..5acd7398b 100755 --- a/examples/example_leeway.py +++ b/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 diff --git a/opendrift/models/basemodel/__init__.py b/opendrift/models/basemodel/__init__.py index c526e7809..7ad7a2ee7 100644 --- a/opendrift/models/basemodel/__init__.py +++ b/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.""" diff --git a/opendrift/models/leeway.py b/opendrift/models/leeway.py index b333efce0..e3001a034 100644 --- a/opendrift/models/leeway.py +++ b/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)) 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 diff --git a/tests/models/test_leeway.py b/tests/models/test_leeway.py index 0ecf79bb6..57efa9b24 100644 --- a/tests/models/test_leeway.py +++ b/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,39 @@ 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': 15, 'y_wind': 0, 'land_binary_mask': 0}) + o.set_config('capsizing', True) + o.set_config('capsizing:wind_threshold', 15) + o.set_config('capsizing:wind_threshold_sigma', 5) + o.set_config('capsizing:leeway_fraction', .4) + o.seed_elements(lon=0, lat=60, time=datetime.now(), number=10) + o.run(time_step=900, time_step_output=900, duration=timedelta(hours=2)) + assert o.elements.capsized.max() == 1 + assert o.elements.capsized.min() == 0 + assert o.elements.capsized.sum() == 6 + + # Backward run, checking that forward capsizing is not happening + ob = Leeway(loglevel=20) + ob.set_config('capsizing', True) + ob.set_config('environment:constant', {'x_sea_water_velocity': 0, 'y_sea_water_velocity': 0, + 'x_wind': 15, 'y_wind': 0, 'land_binary_mask': 0}) + ob.seed_elements(lon=0, lat=60, time=datetime.now(), number=10) + ob.run(time_step=-900, time_step_output=900, duration=timedelta(hours=2)) + 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('environment:constant', {'x_sea_water_velocity': 0, 'y_sea_water_velocity': 0, + 'x_wind': 15, 'y_wind': 0, 'land_binary_mask': 0}) + ob.seed_elements(lon=0, lat=60, time=datetime.now(), number=10, capsized=1) + ob.run(time_step=-900, time_step_output=900, duration=timedelta(hours=2)) + assert ob.elements.capsized.max() == 1 + assert ob.elements.capsized.min() == 0 + assert ob.elements.capsized.sum() == 7 From 68be30dfe177bec426800279e342c67f225b88eb Mon Sep 17 00:00:00 2001 From: Knut-Frode Dagestad Date: Tue, 14 May 2024 15:43:21 +0200 Subject: [PATCH 2/3] Requiring adios_db<1.2 until OpenOil is adapted --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 2234dc927..fde49ab62 100644 --- a/environment.yml +++ b/environment.yml @@ -30,6 +30,6 @@ dependencies: - utm - roaring-landmask>=0.7 - trajan>=0.1.3 - - adios_db + - adios_db<1.2 - copernicusmarine - pip From c3890602d5988085a38b61ebe58c261a56d87db0 Mon Sep 17 00:00:00 2001 From: Knut-Frode Dagestad Date: Tue, 14 May 2024 16:03:57 +0200 Subject: [PATCH 3/3] Updated capsizing test numbers --- tests/models/test_leeway.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/tests/models/test_leeway.py b/tests/models/test_leeway.py index 57efa9b24..1ab9149cf 100644 --- a/tests/models/test_leeway.py +++ b/tests/models/test_leeway.py @@ -81,24 +81,27 @@ def test_leewayrun(tmpdir, test_data): def test_capsize(): o = Leeway(loglevel=20) o.set_config('environment:constant', {'x_sea_water_velocity': 0, 'y_sea_water_velocity': 0, - 'x_wind': 15, 'y_wind': 0, 'land_binary_mask': 0}) + 'x_wind': 25, 'y_wind': 0, 'land_binary_mask': 0}) o.set_config('capsizing', True) - o.set_config('capsizing:wind_threshold', 15) - o.set_config('capsizing:wind_threshold_sigma', 5) + 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=10) - o.run(time_step=900, time_step_output=900, duration=timedelta(hours=2)) + 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() == 6 + 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': 15, 'y_wind': 0, 'land_binary_mask': 0}) - ob.seed_elements(lon=0, lat=60, time=datetime.now(), number=10) - ob.run(time_step=-900, time_step_output=900, duration=timedelta(hours=2)) + '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 @@ -106,10 +109,13 @@ def test_capsize(): # 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': 15, 'y_wind': 0, 'land_binary_mask': 0}) - ob.seed_elements(lon=0, lat=60, time=datetime.now(), number=10, capsized=1) - ob.run(time_step=-900, time_step_output=900, duration=timedelta(hours=2)) + '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() == 7 + assert ob.elements.capsized.sum() == 82