Update to Wavespectra 4.3.0

examples/tutorial_3_LUPA.ipynb

@@ -44,6 +44,7 @@
     "import matplotlib.pyplot as plt\n",
     "import xarray as xr\n",
     "from scipy.optimize import brute\n",
+    "from wavespectra.construct.frequency import pierson_moskowitz\n",
     "\n",
     "import wecopttool as wot\n",
     "\n",
@@ -430,9 +431,8 @@
     "\n",
     "def irregular_wave(hs, tp):\n",
     "    fp = 1/tp\n",
-    "    spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, hs)\n",
-    "    efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n",
-    "    return wot.waves.long_crested_wave(efth,nrealizations=nrealizations)\n",
+    "    efth = pierson_moskowitz(freq=freq, hs=hs, fp=fp)\n",
+    "    return wot.waves.long_crested_wave(efth,nrealizations=nrealizations,direction=0)\n",
     "\n",
     "for case, data in wave_cases.items():\n",
     "    waves[case] = irregular_wave(data['Hs'], data['Tp'])"

examples/tutorial_4_Pioneer.ipynb

@@ -30,7 +30,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -41,6 +41,7 @@
     "from scipy.linalg import block_diag\n",
     "import xarray as xr\n",
     "from math import comb\n",
+    "from wavespectra.construct.frequency import pierson_moskowitz\n",
     "\n",
     "import wecopttool as wot\n",
     "\n",
@@ -75,9 +76,11 @@
     "fend = 1.875\n",
     "nfreq_irreg = 150\n",
     "f1_irreg = fend / nfreq_irreg\n",
+    "freq_irreg = wot.frequency(f1_irreg, nfreq_irreg, False) # False -> no zero frequency\n",
     "\n",
     "f1_reg = .325/2\n",
-    "nfreq_reg = 12"
+    "nfreq_reg = 12\n",
+    "freq_reg = wot.frequency(f1_reg, nfreq_reg, False) # False -> no zero frequency"
    ]
   },
   {
@@ -97,8 +100,16 @@
     "nrealizations = 2\n",
     "\n",
     "fp = 1/Tp\n",
-    "spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, Hs)\n",
-    "efth = wot.waves.omnidirectional_spectrum(f1_irreg, nfreq_irreg, spectrum, \"Pierson-Moskowitz\")\n",
+    "efth = pierson_moskowitz(freq=freq_irreg, hs=Hs, fp=fp)\n",
+    "# Here, we can add directionality to the spectrum. This is done automatically by \n",
+    "# wot.waves.long_crested_wave(), but can also be done manually to add multiple directions.\n",
+    "if 'dir' not in efth.dims: \n",
+    "    efth = xr.DataArray(\n",
+    "        data=np.expand_dims(efth, 1),\n",
+    "        dims=[\"freq\", \"dir\"],\n",
+    "        coords=dict(freq=efth['freq'], dir=np.array([0])),\n",
+    "        name=\"efth\",\n",
+    "    )\n",
     "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations=nrealizations)"
    ]
   },
@@ -278,11 +289,9 @@
    "outputs": [],
    "source": [
     "rho = 1025. # kg/m^3\n",
-    "freq_reg = wot.frequency(f1_reg, nfreq_reg, False) # False -> no zero frequency\n",
     "bem_data_reg = wot.run_bem(pnr_fb, freq_reg)\n",
     "omega_reg = bem_data_reg.omega.values\n",
     "\n",
-    "freq_irreg = wot.frequency(f1_irreg, nfreq_irreg, False) # False -> no zero frequency\n",
     "bem_data_irreg = wot.run_bem(pnr_fb, freq_irreg)\n",
     "omega_irreg = bem_data_irreg.omega.values\n",
     "\n",

pyproject.toml

@@ -25,7 +25,7 @@ dependencies = [
     "autograd",
     "capytaine",
     "joblib",
-    "wavespectra>=3.13",
+    "wavespectra>=4.0",
     "netcdf4",
 ]
 

tests/test_waves.py

@@ -214,13 +214,8 @@ def pm_spectrum(self, pm_f1, pm_nfreq, pm_hs):
         Tp = 1.2
         Hs = pm_hs
 
-        def spectrum_func(f):
-            return wot.waves.pierson_moskowitz_spectrum(f, fp=1/Tp, hs=Hs)
-        spectrum_name = f"Pierson-Moskowitz ({Tp}s, {Hs}m)"
-
-        efth_xr = wot.waves.omnidirectional_spectrum(
-            pm_f1, pm_nfreq, spectrum_func, spectrum_name
-        )
+        efth_xr = ws.construction.frequency.pierson_moskowitz(
+            freq=wot.frequency(pm_f1, pm_nfreq, False), hs=Hs, fp=1/Tp)
         return efth_xr
 
     def test_coordinates(self, elevation):
@@ -249,11 +244,6 @@ def test_realizations(self, elevation):
         realization_out = elevation.realization.values
         assert (realization_out == [0,1]).all()
 
-    def test_spectrum(self, pm_spectrum, pm_hs):
-        """Test that the constructed spectrum has the expected Hs."""
-        efth = ws.SpecArray(pm_spectrum)
-        assert np.isclose(pm_hs, efth.hs().values)
-
     def test_time_series(self, pm_spectrum, pm_f1, pm_nfreq):
         """Test that the created time series has the desired spectrum."""
         # create time-series
@@ -355,140 +345,3 @@ def test_float(self,):
         """
         phase = wot.waves.random_phase()
         assert (phase < np.pi) and (phase >= -np.pi)
-
-
-# TODO: Move everything below to wavespectra.construct
-class TestWaveSpectra:
-    def test_omnidirectional_spectrum(self, f1, nfreq, fp, hs):
-        def spectrum_func(
-            f): return wot.waves.pierson_moskowitz_spectrum(f, fp, hs)
-        wave_spec = wot.waves.omnidirectional_spectrum(
-            f1, nfreq, spectrum_func, "Pierson-Moskowitz")
-
-        # the values should be the same as calling the spectrum function
-        freq = wot.frequency(f1, nfreq, False)
-        spec_test = spectrum_func(freq)
-
-        assert np.allclose(spec_test, wave_spec.values.flatten())
-
-    def test_spectrum(self, f1, nfreq, fp, hs, ndir):
-        s_max = 10
-        directions = np.linspace(0, 360, ndir, endpoint=False)
-        dm = directions[np.random.randint(0, ndir)]
-
-        def spectrum_func(f): 
-            return wot.waves.pierson_moskowitz_spectrum(f, fp, hs)
-
-        def spread_func(f, d): 
-            return wot.waves.spread_cos2s(f, d, dm, fp, s_max)
-        spectrum_name, spread_name = "Pierson-Moskowitz", "Cos2s"
-        wave_spec = wot.waves.spectrum(
-            f1, nfreq, directions, spectrum_func, spread_func,
-            spectrum_name, spread_name)
-
-        # integral over all angles should be equal to omnidirectional
-        spec_omni = wot.waves.omnidirectional_spectrum(
-            f1, nfreq, spectrum_func, spectrum_name)
-        spec_omni = spec_omni.values.flatten()
-        ddir = (wave_spec.dir[1] - wave_spec.dir[0]).values
-        integral_d = wave_spec.sum(dim='dir').values * ddir
-
-        # mean direction
-        dfreq = (wave_spec.freq[1] - wave_spec.freq[0]).values
-        integral_f = wave_spec.sum(dim='freq').values * dfreq
-
-        assert wave_spec.shape == (nfreq, ndir)  # shape
-        assert np.allclose(integral_d, spec_omni, rtol=0.01)
-        assert directions[np.argmax(integral_f)] == dm
-
-    def test_pierson_moskowitz_spectrum(self, f1, nfreq, fp, hs):
-        spectrum = wot.waves.pierson_moskowitz_spectrum
-
-        # scalar
-        freq_1 = 0.4
-        spec1 = spectrum(freq_1, fp, hs)
-
-        # vector
-        freqs = wot.frequency(f1, nfreq, False)
-        spec = spectrum(freqs, fp, hs)
-
-        # total elevation variance
-        freqs_int = np.linspace(0, 10, 1000)[1:]
-        total_variance_calc = np.trapz(spectrum(freqs_int, fp, hs), freqs_int)
-        a_param, b_param = wot.waves.pierson_moskowitz_params(fp, hs)
-        total_variance_theory = a_param/(4*b_param)
-
-        assert isinstance(spec1, float)  # scalar
-        assert spec.shape == freqs.shape  # vector shape
-        assert np.isclose(total_variance_calc,
-                          total_variance_theory)  # integral
-
-    def test_jonswap_spectrum(self, f1, nfreq, fp, hs):
-        spectrum = wot.waves.jonswap_spectrum
-
-        # scalar
-        freq_1 = 0.4
-        spec1 = spectrum(freq_1, fp, hs)
-
-        # vector
-        freqs = wot.frequency(f1, nfreq, False)
-        spec = spectrum(freqs, fp, hs)
-
-        # reduces to PM
-        spec_gamma1 = spectrum(freqs, fp, hs, gamma=1.0)
-        spec_pm = wot.waves.pierson_moskowitz_spectrum(freqs, fp, hs)
-
-        assert isinstance(spec1, float)  # scalar
-        assert spec.shape == freqs.shape  # vector shape
-        assert np.allclose(spec_gamma1, spec_pm)  # reduces to PM
-
-    def test_spread_cos2s(self, f1, nfreq, fp, ndir):
-        """Confirm that energy is spread correctly accross wave directions.
-        Integral over all directions of the spread function gives (vector)
-        1.
-        """
-        directions = np.linspace(0, 360, ndir, endpoint=False)
-        wdir_mean = directions[np.random.randint(0, ndir)]
-        freqs = wot.frequency(f1, nfreq, False)
-        s_max = 10
-        spread = wot.waves.spread_cos2s(freq=freqs,
-                                        directions=directions,
-                                        dm=wdir_mean,
-                                        fp=fp,
-                                        s_max=s_max)
-        ddir = directions[1]-directions[0]
-        dfreq = freqs[1] - freqs[0]
-        integral_d = np.sum(spread, axis=1)*ddir
-        integral_f = np.sum(spread, axis=0)*dfreq
-
-        assert directions[np.argmax(integral_f)] == wdir_mean  # mean dir
-        assert np.allclose(integral_d, np.ones(
-            (1, nfreq)), rtol=0.01)  # omnidir
-
-    def test_general_spectrum(self, f1, nfreq):
-        freq = wot.frequency(f1, nfreq, False)
-        a_param = np.random.random()*10
-        b_param = np.random.random()*10
-        spec_f1 = wot.waves.general_spectrum(a_param, b_param, 1.0)
-        spec_a0 = wot.waves.general_spectrum(0, b_param, freq)
-        spec_b0 = wot.waves.general_spectrum(a_param, 0, freq)
-
-        a_vec = np.random.random(freq.shape)*10
-        b_vec = np.random.random(freq.shape)*10
-        spec_vec = wot.waves.general_spectrum(a_vec, b_vec, freq)
-
-        # types and shapes
-        assert isinstance(spec_f1, float)
-        assert spec_a0.shape == spec_b0.shape == freq.shape
-        assert spec_vec.shape == freq.shape
-        # values
-        assert np.isclose(spec_f1, a_param * np.exp(-b_param))
-        assert np.allclose(spec_a0, 0.0)
-        assert np.allclose(spec_b0, a_param * freq**(-5))
-
-    def test_pierson_moskowitz_params(self, fp, hs):
-        params = wot.waves.pierson_moskowitz_params(fp, hs)
-
-        assert len(params) == 2  # returns two floats
-        for iparam in params:
-            assert isinstance(iparam, float)