Add option for cyclic forcing to driver.run()

Also cleaned up some numpy warnings that were coming about when
fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)) was
causing double-precision overflow (and a similar divide-by-0 error in the
matlab-comparison notebook)

Lastly, realized that I needed the nc-time-axis package for making plots now
that we are using cftime objects

docs/matlab-comparison.ipynb

@@ -61,7 +61,10 @@
     "    nx = np.min([len(matlab_vals[\"X\"].data), len(py_vals[\"X\"].data)])\n",
     "    matlab_vals = matlab_vals.isel(time=slice(0, nstep), X=slice(0, nx)).data\n",
     "    py_vals = py_vals.isel(time=slice(0, nstep), X=slice(0, nx)).data\n",
-    "    rel_errs = np.where(matlab_vals != 0, np.abs((py_vals - matlab_vals) / matlab_vals), np.nan)\n",
+    "    rel_err_denom = np.where(\n",
+    "        matlab_vals != 0, matlab_vals, 1\n",
+    "    )  # avoid dividing by 0 in next np.where() statement\n",
+    "    rel_errs = np.where(matlab_vals != 0, np.abs((py_vals - matlab_vals) / rel_err_denom), np.nan)\n",
     "    rel_errs = np.where(np.logical_and(np.isnan(rel_errs), py_vals == 0), 0, rel_errs)\n",
     "    rel_errs = np.where(np.logical_and(np.isnan(rel_errs), py_vals > 0), np.inf, rel_errs)\n",
     "    rel_errs = np.where(np.logical_and(np.isnan(rel_errs), py_vals < 0), -np.inf, rel_errs)\n",
@@ -131,11 +134,16 @@
     "baselines_from_nc = xr.open_dataset(f'../matlab_baselines/{matlab_script}.nc')\n",
     "\n",
     "if matlab_script == \"test_case\":\n",
+    "    cyclic_forcing = True\n",
     "    testcase = feisty.driver.simulate_testcase(\"tanh_shelf\", \"cyclic\")\n",
     "elif matlab_script == \"test_locs3\":\n",
+    "    cyclic_forcing = True\n",
     "    kwargs = {'forcing_yaml': './forcing.yaml', 'forcing_key': matlab_script}\n",
     "    testcase = feisty.driver.simulate_testcase(\n",
-    "        \"from_disk\", \"from_disk\", domain_kwargs=kwargs, forcing_kwargs=kwargs\n",
+    "        \"from_disk\",\n",
+    "        \"from_disk\",\n",
+    "        domain_kwargs=kwargs,\n",
+    "        forcing_kwargs=kwargs,\n",
     "    )\n",
     "else:\n",
     "    raise ValueError(f\"unknown matlab_script '{matlab_script}'\")\n",
@@ -168,14 +176,6 @@
       "\n",
       "Max diff in depth: 0.0\n"
      ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/glade/scratch/mlevy/tmpdir/ipykernel_262885/1751831492.py:6: RuntimeWarning: invalid value encountered in true_divide\n",
-      "  rel_errs = np.where(matlab_vals != 0, np.abs((py_vals - matlab_vals) / matlab_vals), np.nan)\n"
-     ]
     }
    ],
    "source": [
@@ -199,41 +199,20 @@
    "id": "27086302",
    "metadata": {},
    "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/glade/work/mlevy/codes/feisty/feisty/core/process.py:392: RuntimeWarning: overflow encountered in power\n",
-      "  1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))\n",
-      "/glade/work/mlevy/codes/feisty/feisty/core/process.py:392: RuntimeWarning: overflow encountered in power\n",
-      "  1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))\n",
-      "/glade/work/mlevy/codes/feisty/feisty/core/process.py:392: RuntimeWarning: overflow encountered in power\n",
-      "  1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))\n",
-      "/glade/work/mlevy/codes/feisty/feisty/core/process.py:392: RuntimeWarning: overflow encountered in power\n",
-      "  1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))\n",
-      "/glade/work/mlevy/codes/feisty/feisty/core/process.py:392: RuntimeWarning: overflow encountered in power\n",
-      "  1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))\n",
-      "/glade/work/mlevy/codes/feisty/feisty/core/process.py:392: RuntimeWarning: overflow encountered in power\n",
-      "  1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))\n",
-      "/glade/work/mlevy/codes/feisty/feisty/core/process.py:392: RuntimeWarning: overflow encountered in power\n",
-      "  1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))\n",
-      "/glade/work/mlevy/codes/feisty/feisty/core/process.py:392: RuntimeWarning: overflow encountered in power\n",
-      "  1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))\n"
-     ]
-    },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "CPU times: user 20.3 s, sys: 29.7 ms, total: 20.3 s\n",
-      "Wall time: 20.8 s\n"
+      "CPU times: user 43.2 s, sys: 66.5 ms, total: 43.3 s\n",
+      "Wall time: 44.8 s\n"
      ]
     }
    ],
    "source": [
-    "nsteps = 365\n",
+    "nyears = 2\n",
+    "nsteps = nyears * 365\n",
     "make_plot = False\n",
-    "%time testcase.run(nsteps)"
+    "%time testcase.run(nsteps, cyclic_forcing=cyclic_forcing)"
    ]
   },
   {
@@ -267,15 +246,15 @@
      "text": [
       "| group | Matlab Value | Python Value | Rel Err |\n",
       "| --- | --- | --- | --- |\n",
-      "| Sf (t=298, X=0) | 1.2378e-05 | 1.2378e-05 | 2.7371e-16 |\n",
+      "| Sf (t=535, X=0) | 1.7784e-05 | 1.7784e-05 | 5.7154e-16 |\n",
       "| Sp (t=0, X=0) | 9.9989e-06 | 9.9989e-06 | 0.0000e+00 |\n",
-      "| Sd (t=0, X=0) | 9.9989e-06 | 9.9989e-06 | 0.0000e+00 |\n",
-      "| Mf (t=61, X=0) | 4.1794e-05 | 4.1794e-05 | 3.2427e-16 |\n",
+      "| Sd (t=706, X=1) | 8.7526e-06 | 8.7526e-06 | 1.9355e-16 |\n",
+      "| Mf (t=642, X=0) | 5.2255e-04 | 5.2255e-04 | 4.1497e-16 |\n",
       "| Mp (t=4, X=2) | 1.5865e-05 | 1.5865e-05 | 2.1356e-16 |\n",
       "| Md (t=258, X=0) | 1.5325e-04 | 1.5325e-04 | 3.5373e-16 |\n",
       "| Lp (t=187, X=2) | 3.0789e-04 | 3.0789e-04 | 3.5214e-16 |\n",
       "| Ld (t=214, X=1) | 2.0739e-05 | 2.0739e-05 | 1.7971e-15 |\n",
-      "| benthic_prey (t=327, X=1) | 2.9618e-01 | 2.9618e-01 | 1.6868e-15 |\n"
+      "| benthic_prey (t=370, X=1) | 3.3466e-01 | 3.3466e-01 | 1.9905e-15 |\n"
      ]
     }
    ],

docs/requirements-docs.txt

@@ -3,4 +3,5 @@ pydata-sphinx-theme
 nbsphinx
 ipython
 ipykernel
+nc-time-axis
 - r ../requirements.txt

feisty/core/process.py

@@ -388,12 +388,15 @@ def compute_growth(
         death = predation_rate.data[i, :] + mortality_rate.data[i, :] + fish_catch_rate.data[i, :]
         somatic_growth_potential = fish.energy_frac_somatic_growth * energy_avail_rate.data[i, :]
 
-        gg = (somatic_growth_potential - death) / (
-            1.0 - (fish.size_class_bnds_ratio ** (1.0 - (death / somatic_growth_potential)))
+        exp_in_denom = np.log10(fish.size_class_bnds_ratio) * (
+            1.0 - (death / somatic_growth_potential)
         )
+        exp_in_denom = np.where(exp_in_denom < 300, exp_in_denom, 300)
+        gg = (somatic_growth_potential - death) / (1.0 - (10.0 ** exp_in_denom))
         growth_rate.data[i, :] = np.where(
             gg < energy_avail_rate.data[i, :], gg, energy_avail_rate.data[i, :]
         )
+        # print(f"{log_of_exp} -> {gg}")
         lndx = np.isnan(gg) | (gg < 0)
         growth_rate.data[i, lndx] = 0.0
 

feisty/driver.py

@@ -124,8 +124,15 @@ def __init__(
             benthic_prey_ic_data=benthic_prey_ic_data,
         )
 
-    def _forcing_t(self, t):
-        return self.forcing.interp(time=t)
+    def _forcing_t(self, t, ignore_year):
+        if ignore_year:
+            units = 'days since 0001-01-01 00:00:00'
+            interp_time = cftime.num2date(
+                cftime.date2num(t, units) % 365, units, calendar=t.calendar
+            )
+        else:
+            interp_time = t
+        return self.forcing.interp(time=interp_time)
 
     def _init_output_arrays(self, nt):
         self.time = xr.cftime_range(start=self.start_date, periods=nt)
@@ -145,9 +152,9 @@ def ds(self):
         """Data comprising the output from a ``feisty`` simulation."""
         return self._ds
 
-    def _compute_tendency(self, t, state_t):
+    def _compute_tendency(self, t, state_t, cyclic_forcing):
         """Return the feisty time tendency."""
-        gcm_data_t = self._forcing_t(t)
+        gcm_data_t = self._forcing_t(t, cyclic_forcing)
         return self.obj.compute_tendencies(
             state_t.isel(group=self.obj.prog_ndx_fish),
             state_t.isel(group=self.obj.prog_ndx_benthic_prey),
@@ -158,35 +165,35 @@ def _compute_tendency(self, t, state_t):
             poc_flux=gcm_data_t.poc_flux_bottom,
         )
 
-    def _solve(self, nt, method):
+    def _solve(self, nt, method, cyclic_forcing):
         """Call a numerical ODE solver to integrate the feisty model in time."""
 
         state_t = self.obj.get_prognostic().copy()
         self._init_output_arrays(nt)
         if method == 'euler':
-            self._solve_foward_euler(nt, state_t)
+            self._solve_foward_euler(nt, state_t, cyclic_forcing)
 
         elif method in ['Radau', 'RK45']:
             # TODO: make input arguments
-            self._solve_scipy(nt, state_t, method)
+            self._solve_scipy(nt, state_t, method, cyclic_forcing)
         else:
             raise ValueError(f'unknown method: {method}')
 
-    def _solve_foward_euler(self, nt, state_t):
+    def _solve_foward_euler(self, nt, state_t, cyclic_forcing):
         """use forward-euler to solve feisty model"""
         for n in range(nt):
-            dsdt = self._compute_tendency(self.time[n], state_t)
+            dsdt = self._compute_tendency(self.time[n], state_t, cyclic_forcing)
             state_t[self.obj.prog_ndx_prognostic, :] = (
                 state_t[self.obj.prog_ndx_prognostic, :]
                 + dsdt[self.obj.ndx_prognostic, :] * self.dt
             )
             self._post_data(n, state_t)
 
-    def _solve_scipy(self, nt, state_t, method):
+    def _solve_scipy(self, nt, state_t, method, cyclic_forcing):
         """use a SciPy solver to integrate the model equation."""
         raise NotImplementedError('scipy solvers not implemented')
 
-    def run(self, nt, file_out=None, method='euler'):
+    def run(self, nt, file_out=None, method='euler', cyclic_forcing=False):
         """Integrate the FEISTY model.
 
         Parameters
@@ -205,7 +212,7 @@ def run(self, nt, file_out=None, method='euler'):
              Only ``method='euler'`` is supported currently.
 
         """
-        self._solve(nt, method)
+        self._solve(nt, method, cyclic_forcing)
         self._shutdown(file_out)
 
     def _shutdown(self, file_out):

tests/test_driver.py

@@ -96,5 +96,5 @@ def test_cyclic_interpolation():
     testcase1 = feisty.driver.simulate_testcase('tanh_shelf', 'cyclic')
     testcase2 = feisty.driver.simulate_testcase('tanh_shelf', 'cyclic', start_date='0002-01-01')
     testcase1.run(1)
-    testcase2.run(1)
+    testcase2.run(1, cyclic_forcing=True)
     assert (testcase1.ds['biomass'].data == testcase2.ds['biomass'].data).all()