Merge pull request #24 from mnlevy1981/timestep_bp_with_fish

Timestep benthic prey with fish

feisty/core/interface.py

@@ -70,7 +70,7 @@ def __init__(
         self._init_benthic_prey_settings(self.settings_dict['benthic_prey'])
 
         fish_ic_data = 1e-5 if fish_ic_data is None else fish_ic_data
-        benthic_prey_ic_data = 1e-5 if benthic_prey_ic_data is None else benthic_prey_ic_data
+        benthic_prey_ic_data = 2e-3 if benthic_prey_ic_data is None else benthic_prey_ic_data
         self._init_biomass(fish_ic_data, benthic_prey_ic_data)
 
         self._init_food_web(self.settings_dict['food_web'])
@@ -152,10 +152,11 @@ def _init_biomass(self, fish_ic_data, benthic_prey_ic_data):
         self.ndx_zoo = np.arange(0, self.n_zoo)
         self.ndx_fish = np.arange(self.n_zoo, self.n_zoo + self.n_fish)
         self.ndx_benthic_prey = np.arange(self.n_zoo + self.n_fish, n)
-
         self.ndx_prognostic = np.concatenate((self.ndx_fish, self.ndx_benthic_prey))
+
         self.prog_ndx_fish = np.arange(0, self.n_fish)
         self.prog_ndx_benthic_prey = np.arange(self.n_fish, self.n_fish + self.n_benthic_prey)
+        self.prog_ndx_prognostic = np.concatenate((self.prog_ndx_fish, self.prog_ndx_benthic_prey))
 
         # TODO: make private
         self.biomass = domain.init_array_2d('group', group_coord, name='biomass')
@@ -186,6 +187,7 @@ def _init_tendency_arrays(self):
             self.zooplankton,
             self.fish,
             self.benthic_prey,
+            self.member_obj_list,
             self.food_web,
         )
 
@@ -278,17 +280,6 @@ def _compute_temperature(self):
                 self.tendency_data.t_frac_pelagic[i, :],
             )
 
-    def _update_benthic_biomass(self):
-        process.compute_benthic_biomass_update(
-            self.tendency_data.benthic_biomass_new,
-            self.tendency_data.consumption_rate_link,
-            self.benthic_prey,
-            self.biomass,
-            self.food_web,
-            self.gcm_state.poc_flux,
-        )
-        self._set_benthic_prey_biomass(self.tendency_data.benthic_biomass_new)
-
     def _compute_pred_encounter_consumption_max(self):
         process.compute_pred_encounter_consumption_max(
             self.tendency_data.encounter_rate_pred,
@@ -416,7 +407,11 @@ def _compute_total_tendency(self):
             self.tendency_data.predation_flux,
             self.tendency_data.fish_catch_rate,
             self.biomass,
-            self.fish,
+            self.tendency_data.consumption_rate_link,
+            self.food_web,
+            self.gcm_state.poc_flux,
+            self.member_obj_list,
+            self.ndx_fish[0],
         )
 
     def compute_tendencies(
@@ -455,9 +450,6 @@ def compute_tendencies(
         self._set_zoo_mortality(zoo_mortality_data)
         self.gcm_state.update(**gcm_state_update_kwargs)
 
-        # advance benthic prey concentrations
-        self._update_benthic_biomass()
-
         # compute temperature terms
         self._compute_t_frac_pelagic()
         self._compute_temperature()
@@ -483,12 +475,13 @@ def compute_tendencies(
         return self.tendency_data.total_tendency
 
 
-def _init_tendency_data(zoo_list, fish_list, benthic_prey_list, food_web):
+def _init_tendency_data(zoo_list, fish_list, benthic_prey_list, member_obj_list, food_web):
     """Return an xarray.Dataset with initialized tendency data arrays."""
 
     fish_names = [f.name for f in fish_list]
     zoo_names = [z.name for z in zoo_list]
     benthic_prey_names = [b.name for b in benthic_prey_list]
+    member_obj_names = [mo.name for mo in member_obj_list]
 
     pred_names = [link.predator.name for link in food_web]
     prey_names = [link.prey.name for link in food_web]
@@ -599,7 +592,7 @@ def _init_tendency_data(zoo_list, fish_list, benthic_prey_list, food_web):
     # TODO: include benthic_prey once timestepping has been moved out of feisty
     ds['total_tendency'] = domain.init_array_2d(
         coord_name='group',
-        coord_values=fish_names,
+        coord_values=member_obj_names,
         name='total_tendency',
         attrs={'long_name': 'Total time tendency'},
     )

feisty/core/process.py

@@ -351,45 +351,6 @@ def compute_natural_mortality(mortality_rate, fish_list, T_habitat, mortality_ty
             raise ValueError(f'unknown mortality type {fish.mortality_type}')
 
 
-def compute_benthic_biomass_update(
-    benthic_biomass_new, consumption_rate_link, benthic_prey_list, biomass, food_web, poc_flux
-):
-    """
-    bio_in = benthic biomass
-    det = poc_flux flux to bottom (g/m2/d)
-    con = biomass specific consumption rate by MD & LD
-    bio = biomass of MD & LD
-    """
-
-    for i, benthic_prey in enumerate(benthic_prey_list):
-        # eaten = consumption * biomass_pred
-        # pred = sum(eaten, 2)
-        biomass_bent = biomass.sel(group=benthic_prey.name)
-        predation = (
-            biomass.isel(group=food_web.prey_ndx_pred[benthic_prey.name])
-            * food_web.get_consumption(consumption_rate_link, prey=benthic_prey.name)
-        ).sum('group')
-
-        # Needs to be in units of per time (g/m2/d) * (g/m2)
-        growth = benthic_prey.benthic_efficiency * poc_flux
-
-        if not benthic_prey.lcarrying_capacity:  # no carrying capacity
-            benthic_biomass_new[i, :] = biomass_bent + growth - predation
-        else:
-            # logistic
-            benthic_biomass_new[i, :] = (
-                biomass_bent
-                + growth * (1.0 - biomass_bent / benthic_prey.carrying_capacity)
-                - predation
-            )
-
-    benthic_biomass_new[:, :] = np.where(
-        benthic_biomass_new < 0.0,
-        constants.eps,
-        benthic_biomass_new,
-    )
-
-
 def compute_energy_avail(energy_avail_rate, ingestion_rate, metabolism_rate, fish_list):
     """Compute energy available for growth (nu)."""
 
@@ -470,7 +431,11 @@ def compute_total_tendency(
     predation_flux,
     fish_catch_rate,
     biomass,
-    fish_list,
+    consumption_rate_link,
+    food_web,
+    poc_flux,
+    member_obj_list,
+    first_fish_id,
 ):
     """
     Compute the total time tendency of fish.
@@ -497,24 +462,45 @@ def compute_total_tendency(
 
     biomass : array_like
 
-    fish_list : list
+    member_obj_list : list
       List of feisty.ecosystem.fish object.
     """
-    for i, fish in enumerate(fish_list):
-        total_tendency[i, :] = (
-            recruitment_flux[i, :]
-            + biomass.sel(group=fish.name)
-            * (
-                (
-                    energy_avail_rate[i, :]
-                    - reproduction_rate[i, :]
-                    - growth_rate[i, :]
-                    - mortality_rate[i, :]
-                    - fish_catch_rate[i, :]
+    for i, member_obj in enumerate(member_obj_list):
+        if type(member_obj) == ecosystem.fish_type:
+            fish_i = i - first_fish_id
+            total_tendency[i, :] = (
+                recruitment_flux[fish_i, :]
+                + biomass[i, :]
+                * (
+                    energy_avail_rate[fish_i, :]
+                    - reproduction_rate[fish_i, :]
+                    - growth_rate[fish_i, :]
+                    - mortality_rate[fish_i, :]
+                    - fish_catch_rate[fish_i, :]
                 )
+                - predation_flux[fish_i, :]
             )
-            - predation_flux[i, :]
-        )
+        elif type(member_obj) == ecosystem.benthic_prey_type:
+            # eaten = consumption * biomass_pred
+            # pred = sum(eaten, 2)
+
+            predation = (
+                biomass.isel(group=food_web.prey_ndx_pred[member_obj.name])
+                * food_web.get_consumption(consumption_rate_link, prey=member_obj.name)
+            ).sum('group')
+
+            # Needs to be in units of per time (g/m2/d) * (g/m2)
+            growth = member_obj.benthic_efficiency * poc_flux
+
+            if not member_obj.lcarrying_capacity:  # no carrying capacity
+                total_tendency[i, :] = growth - predation
+            else:
+                # logistic
+                total_tendency[i, :] = (
+                    growth
+                    * (1.0 - biomass.sel(group=member_obj.name) / member_obj.carrying_capacity)
+                    - predation
+                )
 
 
 def compute_fish_catch(fish_catch_rate, fishing_rate, fish_list):

feisty/driver.py

@@ -161,10 +161,10 @@ def _solve(self, nt, method):
     def _solve_foward_euler(self, nt, state_t):
         """use forward-euler to solve feisty model"""
         for n in range(nt):
-            dfdt = self._compute_tendency(self.time[n], state_t)
-            state_t[self.obj.prog_ndx_fish, :] = state_t[self.obj.prog_ndx_fish, :] + dfdt * self.dt
-            state_t[self.obj.prog_ndx_benthic_prey, :] = self.obj.biomass.isel(
-                group=self.obj.ndx_benthic_prey
+            dsdt = self._compute_tendency(self.time[n], state_t)
+            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)
 

tests/test_init.py

@@ -398,7 +398,7 @@ def test_init_tendency_arrays():
             checked.append(key)
         elif key in group_coord_vars:
             assert da.dims == ('group', 'X')
-            assert da.shape == (n_fish, NX)  # TODO: fix when benthic prey are folded in
+            assert da.shape == (n_zoo + n_fish + n_benthic_prey, NX)
             checked.append(key)
         elif key in benthic_prey_coord_vars:
             assert da.dims == ('benthic_prey', 'X')

tests/test_process.py

@@ -146,15 +146,6 @@ def test_t_frac_pelagic():
         assert (fish.t_frac_pelagic == fish_settings['members'][i]['t_frac_pelagic_static']).all()
 
 
-@pytest.mark.weak
-def test_update_benthic_prey():
-    """test benthic update
-
-    Add regression data
-    """
-    F._update_benthic_biomass()
-
-
 def test_compute_metabolism():
     F.gcm_state.update(
         T_pelagic=10.0,