Merge pull request #1639 from knutfrode/dev

[run-ex] reader_oscillating now takes period as timedelta instead of …

Author: knutfrode

examples/example_river_runoff.py

@@ -24,9 +24,9 @@
 o.set_config('general:coastline_action', 'previous')
 t1 = datetime.now()
 t2 = t1 + timedelta(hours=48)
-reader_x = reader_oscillating.Reader('x_sea_water_velocity', period_seconds=3600*24,
+reader_x = reader_oscillating.Reader('x_sea_water_velocity', period=timedelta(hours=24),
                 amplitude=1, zero_time=t1)
-reader_y = reader_oscillating.Reader('y_sea_water_velocity', period_seconds=3600*72,
+reader_y = reader_oscillating.Reader('y_sea_water_velocity', period=timedelta(hours=72),
                 amplitude=.5, zero_time=t2)
 o.add_reader([reader_x, reader_y])
 number = 25000

examples/example_seafloor_interaction.py

@@ -25,7 +25,7 @@
 # readers
 o = OceanDrift(loglevel=50)
 reader_norkyst = reader_netCDF_CF_generic.Reader(tdf + '14Jan2016_NorKyst_z_3d/NorKyst-800m_ZDEPTHS_his_00_3Dsubset.nc')
-reader_osc = reader_oscillating.Reader('x_sea_water_velocity', amplitude=10, period_seconds=3600)
+reader_osc = reader_oscillating.Reader('x_sea_water_velocity', amplitude=10, period=timedelta(hours=1))
 
 runs = []
 seafloor_actions = ['previous', 'deactivate', 'lift_to_seafloor']

examples/example_sediments_resuspension.py

@@ -13,7 +13,7 @@
 reader_oscx = reader_oscillating.Reader('x_sea_water_velocity',
     amplitude=0.6, zero_time=datetime.now())
 reader_oscy = reader_oscillating.Reader('y_sea_water_velocity',
-    amplitude=.3, period_seconds=3600*5, zero_time=datetime.now())
+    amplitude=.3, period=timedelta(hours=5), zero_time=datetime.now())
 
 o = SedimentDrift(loglevel=50)  # 0 for debug output
 

opendrift/models/basemodel/__init__.py

@@ -23,6 +23,7 @@
 
 import sys
 import os
+import copy
 import types
 from typing import Union, List
 import traceback
@@ -241,8 +242,7 @@ def __init__(self,
         # List to store GeoJSON dicts of seeding commands
         self.seed_geojson = []
 
-        self.required_variables = self.required_variables.copy()  # Avoid modifying the class variable
-
+        self.required_variables = copy.deepcopy(self.required_variables)  # Avoid modifying the class variable
         self.env = Environment(self.required_variables, self._config)
 
         # Make copies of dictionaries so that they are private to each instance
@@ -579,16 +579,26 @@ def post_run(self):
     def update_previous_state(self):
         """Store some properties and variables, for access at next time step"""
 
-        # Retrieve previous from storage
+        if self.newly_seeded_IDs is not None:  # For newly seeded elements, we set previous equal to present
+            newly_seeded_rel_indices = np.where(self.elements.age_seconds==self.time_step.total_seconds())[0]    
+
+        # Environment variables
         if self._environment_previous is not None:
             # Update previous environment
+            if self.newly_seeded_IDs is not None:
+                for var in self._environment_previous:
+                    self._environment_previous[var][newly_seeded_rel_indices] = self.environment[var][newly_seeded_rel_indices]
             self.environment_previous = self._environment_previous.isel(trajectory=self.elements.ID)
             # Store present environment variables
             for var in self._environment_previous:
                 self._environment_previous[var][self.elements.ID] = self.environment[var]
 
+        # Element properties
         if self._elements_previous is not None:
             # Update previous elements
+            if self.newly_seeded_IDs is not None:  # For new elements there is no previous, so using present
+                for var in self._elements_previous:
+                    self._elements_previous[var][newly_seeded_rel_indices] = getattr(self.elements, var)[newly_seeded_rel_indices]
             self.elements_previous = self._elements_previous.isel(trajectory=self.elements.ID)
             # Store present element properties
             for var in self._elements_previous:
@@ -1722,16 +1732,26 @@ def run(self,
                     logger.info(f'Skipping environment variable {vn} because of condition {var["skip_if"]}')
                     self.required_variables.pop(vn)
 
-        ####################################################################################
-        # Evaluate conditionals to determine if previous element properties shall be stored
-        ####################################################################################
+        #########################################################################################
+        # Evaluate conditionals to determine if previous properties or variables shall be stored
+        #########################################################################################
+        # Element properties
         for en, prop in self.elements.variables.copy().items():
             if 'store_previous_if' in prop:
                 store = evaluate_conditional(*prop['store_previous_if'], self)
                 if store is True:
                     logger.info(f'Storing previous values of element property {en} because of condition {prop["store_previous_if"]}')
                     self.elements.variables[en]['store_previous'] = True
                 del self.elements.variables[en]['store_previous_if']  # To avoid writing this to netCDF metadata
+        # Environment variables
+        for en, var in self.required_variables.copy().items():
+            if 'store_previous_if' in var:
+                store = evaluate_conditional(*var['store_previous_if'], self)
+                if store is True:
+                    logger.info(f'Storing previous values of environment variable {en} because of condition {var["store_previous_if"]}')
+                    self.required_variables[en]['store_previous'] = True
+                del self.required_variables[en]['store_previous_if']  # To avoid writing this to netCDF metadata
+
 
         ########################
         # Simulation time step

opendrift/models/oceandrift.py

@@ -70,7 +70,8 @@ class OceanDrift(OpenDriftSimulation):
     required_variables = {
         'x_sea_water_velocity': {'fallback': 0},
         'y_sea_water_velocity': {'fallback': 0},
-        'sea_surface_height': {'fallback': 0},
+        'sea_surface_height': {'fallback': 0,
+            'store_previous_if': ['drift:vertical_advection', 'is', True]},
         'x_wind': {'fallback': 0},
         'y_wind': {'fallback': 0},
         'upward_sea_water_velocity': {'fallback': 0,
@@ -177,6 +178,8 @@ def __init__(self, *args, **kwargs):
     def update(self):
         """Update positions and properties of elements."""
 
+        self.water_column_stretching()
+
         # Simply move particles with ambient current
         self.advect_ocean_current()
 
@@ -223,7 +226,6 @@ def simulate_trajectories(self, outfile, trajectories, number=1,
             if np.isnan(lo):
                 continue
             self.seed_elements(lon=lo, lat=la, time=time, number=number, **kwargs)
-        print(self)
         self.run(outfile=outfile, end_time=start_times[-1]+simulation_duration)
         # Simulate and save to file
 
@@ -287,6 +289,12 @@ def prepare_vertical_mixing(self):
     def prepare_run(self):
         super(OceanDrift, self).prepare_run()
 
+    def water_column_stretching(self):
+        """If sea_water_level changes, adjust z for continuity"""
+
+        return  # To be implemented
+        delta_zeta = self.environment.sea_surface_height - self.environment_previous.sea_surface_height
+
     def vertical_advection(self):
         """Move particles vertically according to vertical ocean current
 

opendrift/readers/reader_oscillating.py

@@ -14,7 +14,7 @@
 #
 # Copyright 2015, Knut-Frode Dagestad, MET Norway
 
-from datetime import datetime
+from datetime import datetime, timedelta
 import numpy as np
 
 from opendrift.readers.basereader import BaseReader, ContinuousReader
@@ -23,13 +23,15 @@
 class Reader(BaseReader, ContinuousReader):
     '''Returning values oscillating in time with given amplitude and period'''
 
-    def __init__(self, variable, amplitude, period_seconds=3600*24,
+    def __init__(self, variable, amplitude, period=timedelta(hours=24), period_seconds=None,
                  phase=0, zero_time=datetime(2017, 1, 1, 0)):
         '''init with '''
 
         self.variables = [variable]
         self.amplitude = amplitude
-        self.period_seconds = period_seconds
+        if period_seconds is not None:
+            raise ValueError('Input parameter "period_seconds" is deprecated, please use "period" (timedelta) instead')
+        self.period_seconds = period.total_seconds()
         self.zero_time = zero_time
         self.proj4 = '+proj=latlong +datum=WGS84'
         self.xmin = -180

tests/models/test_basemodel.py

@@ -10,7 +10,7 @@
 def test_logging(tmpdir, capsys):
     # Accepting small variations in log output,
     # depending on machine, and from which folder test is run
-    accepted = (220, 223)
+    accepted = (221, 224)
 
     # Logging to console
     logfile = None

tests/models/test_environment.py

@@ -1,7 +1,9 @@
 from datetime import datetime, timedelta
 import pytest
+import numpy as np
 from opendrift.models.basemodel.environment import Environment
 from opendrift.models.oceandrift import OceanDrift
+from opendrift.readers import reader_oscillating
 from opendrift.config import Configurable
 
 
@@ -28,6 +30,7 @@ def test_add_readers(test_data_roms):
 def test_previous():
     o = OceanDrift()  # With coastline_action = none
     o.set_config('general:coastline_action', 'none')
+    o.set_config('drift:vertical_advection', False)
     o.set_config('environment:constant:land_binary_mask', 0)
     o.set_config('environment:constant:x_sea_water_velocity', 1)
     o.seed_elements(lon=3, lat=60, time=datetime.now())
@@ -38,6 +41,7 @@ def test_previous():
 
     o = OceanDrift()  # With coastline_action = previous
     o.set_config('general:coastline_action', 'previous')
+    o.set_config('drift:vertical_advection', False)
     o.set_config('environment:constant:land_binary_mask', 0)
     o.set_config('environment:constant:x_sea_water_velocity', 1)
     o.seed_elements(lon=3, lat=60, time=datetime.now())
@@ -47,6 +51,31 @@ def test_previous():
     assert o.elements_previous is not None
     assert o.environment_previous is None
 
+    o = OceanDrift()  # For newly seeded, previous shall equal present
+    o.seed_elements(lon=3, lat=60, time=[datetime.now(), datetime.now()+timedelta(hours=1)], number=10)
+    o.run(steps=1, time_step=timedelta(minutes=30))
+    assert len(o.environment.sea_surface_height) == 5
+    assert len(o.environment_previous.sea_surface_height) == 5
+    assert len(o._environment_previous.sea_surface_height) == 10
+    assert o.environment_previous.sea_surface_height[0] == 0
+    assert np.isnan(o._environment_previous.sea_surface_height[-1])
+
+    o = OceanDrift()  # Check gradient
+    time = datetime.now()
+    reader_tidal = reader_oscillating.Reader('sea_surface_height', amplitude=1,
+                                             period=timedelta(hours=6), phase=0, zero_time=time)
+    o.add_reader(reader_tidal)
+    o.seed_elements(lon=3, lat=60, time=[time, time+timedelta(hours=1)], number=10)
+    o.run(steps=2, time_step=timedelta(minutes=30))
+    assert len(o.environment.sea_surface_height) == 9
+    assert len(o.environment_previous.sea_surface_height) == 9
+    assert len(o._environment_previous.sea_surface_height) == 10
+    assert o.environment.sea_surface_height[0] == pytest.approx(0.2588, 3)
+    assert o.environment_previous.sea_surface_height[0] == 0
+    assert np.isnan(o._environment_previous.sea_surface_height[-1])
+    assert not np.isnan(o._environment_previous.sea_surface_height[-2])
+
+
 def test_skip_env_variable():
     o = OceanDrift()
     o.set_config('drift:vertical_mixing', True)  # Diffusivity shall be included

tests/models/test_openberg.py

@@ -85,7 +85,7 @@ def test_grounding_and_degrounding():
     o = OpenBerg(loglevel=0)
     reader_osc = reader_oscillating.Reader('sea_surface_height',
                                                 amplitude=2,
-                                                period_seconds=3600 * 6,
+                                                period=timedelta(hours=6),
                                                 zero_time=datetime(2020,1,1))
     o.add_reader(reader_osc)
     o.set_config('environment:constant:x_sea_water_velocity', 0)

tests/models/test_stranding.py

@@ -90,7 +90,7 @@ def test_stranding_roms(self):
     def test_stranding_options(self):
         reader_osc = reader_oscillating.Reader(
                 'x_sea_water_velocity', amplitude=.5,
-                period_seconds=3600*6,
+                period=timedelta(hours=6),
                 zero_time=datetime.now())
 
         reader_global = reader_global_landmask.Reader()