Fix small bugs and add CDC documentation (#291)

* added t=0 case to E_Advective_Dispersion

* show CDC coagulant functions

* fix cdc doctest, ut.ceil_nearest, ut.floor_nearest

* allowed 0C temperature, began documenting cdc.py

* Fix doctest and add method specs in cdc

* fix physchem tests and warnings

* increment version number

* deprecate cdc.coag_q_max_est and test cdc.py

* add missing cdc test

aguaclara/core/physchem.py

@@ -113,7 +113,7 @@ def viscosity_dynamic_water(Temperature):
     :return: dynamic viscosity of water
     :rtype: u.kg/(u.m*u.s)
     """
-    ut.check_range([Temperature.magnitude, ">0", "Temperature in Kelvin"])
+    ut.check_range([Temperature.magnitude, ">=0", "Temperature in Kelvin"])
     return 2.414 * (10**-5) * u.kg/(u.m*u.s) * 10**(247.8*u.degK /
                                                     (Temperature - 140*u.degK))
 
@@ -139,7 +139,7 @@ def density_water(Temperature=None, *, temp=None):
                       UserWarning)
         Temperature = temp
 
-    ut.check_range([Temperature.magnitude, ">0", "Temperature in Kelvin"])
+    ut.check_range([Temperature.magnitude, ">=0", "Temperature in Kelvin"])
     rhointerpolated = interpolate.CubicSpline(WATER_DENSITY_TABLE[0],
                                               WATER_DENSITY_TABLE[1])
     Temperature = Temperature.to(u.degK).magnitude
@@ -168,7 +168,7 @@ def viscosity_kinematic_water(Temperature):
     :return: kinematic viscosity of water
     :rtype: u.m**2/u.s
     """
-    ut.check_range([Temperature.magnitude, ">0", "Temperature in Kelvin"])
+    ut.check_range([Temperature.magnitude, ">=0", "Temperature in Kelvin"])
     return (viscosity_dynamic_water(Temperature) / density_water(Temperature))
 
 ####################### Hydraulic Radius #######################
@@ -1680,6 +1680,9 @@ def re_ergun(ApproachVel, DiamMedia, Temperature, Porosity):
     ut.check_range([ApproachVel.magnitude, ">0", "ApproachVel"],
                    [DiamMedia.magnitude, ">0", "DiamMedia"],
                    [Porosity, "0-1", "Porosity"])
+    if Porosity == 1:
+        raise ValueError("Porosity is " + str(Porosity) + " must be great than\
+            or equal to 0 and less than 1")
     return (ApproachVel * DiamMedia /
             (viscosity_kinematic_water(Temperature)
              * (1 - Porosity))).to(u.dimensionless)

aguaclara/core/utility.py

@@ -166,18 +166,26 @@ def floor_nearest(x, array):
         - ``x``: Value to compare
         - ``array (numpy.array)``: Array to search
     """
-    i = np.argmax(array >= x) - 1
-    return array[i]
+    sorted_array = np.sort(array)[::-1]  # [::-1] syntax reverses array
+    if x < sorted_array[-1]:
+        raise ValueError(str(x) + " is smaller than all values in the array.")
+    else:
+        i = np.argmax(sorted_array <= x)
+        return sorted_array[i]
 
 def ceil_nearest(x, array):
-    """Get the nearest element of a NumPy array less than or equal to a value.
+    """Get the nearest element of a NumPy array greater than or equal to a value.
 
     Args:
         - ``x``: Value to compare
         - ``array (numpy.array)``: Array to search
     """
-    i = np.argmax(array >= x)
-    return array[i]
+    sorted_array = np.sort(array)
+    if x > sorted_array[-1]:
+        raise ValueError(str(x) + " is larger than all values in the array.")
+    else:
+        i = np.argmax(sorted_array >= x)
+        return sorted_array[i]
 
 def _minmax(*args, func=np.max):
     """Get the minuimum/maximum value of some Pint quantities with units.

aguaclara/design/cdc.py

@@ -4,15 +4,15 @@
 
 Example:
     >>> from aguaclara.design.cdc import *
-    >>> cdc = CDC(q = 20 * L/s, coag_type = 'pacl')
+    >>> cdc = CDC(q = 20 * u.L/u.s, coag_type = 'pacl')
     >>> cdc.coag_stock_vol
-    208.198 liter
+    <Quantity(2500.0, 'liter')>
 """
 import aguaclara.core.physchem as pc
 import aguaclara.core.utility as ut
 from aguaclara.core.units import u
-import aguaclara.core.constants as con
 from aguaclara.design.component import Component
+import warnings
 
 import numpy as np
 
@@ -23,13 +23,14 @@ class CDC(Component):
         - ``q (float * u.L / u.s)``: Flow rate (required)
     """
     def __init__(self, **kwargs):
-        self.hl = 20 * u.cm,
+        self.hl = 20 * u.cm
         self.coag_type='pacl'
         if self.coag_type.lower() not in ['pacl', 'alum']:
             raise ValueError('coag_type must be either PACl or Alum.')
 
-        self.coag_dose_conc_max=2 * u.g / u.L #What should this default to? -Oliver L., 6 Jun 19
-        self.coag_stock_conc_est=150 * u.g / u.L
+        self.coag_dose_conc_max=100 * u.mg / u.L
+        self.coag_stock_conc=150 * u.g / u.L
+        self.coag_stock_conc_est=150 * u.g / u.L # Deprecated since January 2021
         self.coag_stock_min_est_time=1 * u.day
         self.chem_tank_vol_supplier=[208.198, 450, 600, 750, 1100, 2500] * u.L
         self.chem_tank_dimensions_supplier=[
@@ -48,8 +49,9 @@ def __init__(self, **kwargs):
 
         super().__init__(**kwargs)
 
-    def _alum_nu(self, coag_conc):
-        """Return the dynamic viscosity of water at a given temperature.
+    def alum_nu(self, coag_conc):
+        """
+        Return the dynamic viscosity of water at a given temperature.
 
         If given units, the function will automatically convert to Kelvin.
         If not given units, the function will assume Kelvin.
@@ -61,92 +63,153 @@ def _alum_nu(self, coag_conc):
             (1 + (4.255 * 10 ** -6) * coag_conc.magnitude ** 2.289) * \
             pc.viscosity_kinematic_water(self.temp)
         return alum_nu
-
-    def _pacl_nu(self, coag_conc):
+
+    def _alum_nu(self, coag_conc):
+        """
+        .. deprecated::
+            `_alum_nu` is deprecated; use `alum_nu` instead.
+        """
+        # Deprecated since January 2021
+        warnings.warn('_alum_nu is deprecated; use alum_nu instead.',
+                  UserWarning)
+
+        return self.alum_nu(coag_conc)
+
+    def pacl_nu(self, coag_conc):
         """Return the dynamic viscosity of water at a given temperature.
 
         If given units, the function will automatically convert to Kelvin.
         If not given units, the function will assume Kelvin.
         This function assumes that the temperature dependence can be explained
-        based on the effect on water and that there is no confounding effect from
-        the coagulant.
+        based on the effect on water and that there is no confounding effect 
+        from the coagulant.
         """
         pacl_nu = \
             (1 + (2.383 * 10 ** -5) * (coag_conc).magnitude ** 1.893) * \
             pc.viscosity_kinematic_water(self.temp)
         return pacl_nu
 
+    def _pacl_nu(self, coag_conc):
+        """
+        .. deprecated::
+            `_pacl_nu` is deprecated; use `pacl_nu` instead.
+        """
+        # Deprecated since January 2021
+        warnings.warn('_pacl_nu is deprecated; use pacl_nu instead.',
+                  UserWarning)
+
+        return self.pacl_nu(coag_conc)
+
     def _coag_nu(self, coag_conc, coag_type):
+        """
+        .. deprecated::
+            `_coag_nu` is deprecated; use `coag_nu` instead.
+        """
+        # Deprecated since January 2021
+        warnings.warn('_coag_nu is deprecated; use coag_nu instead.',
+                  UserWarning)
+
+        return self.coag_nu(coag_conc, coag_type)
+
+    def coag_nu(self, coag_conc, coag_type):
         """Return the dynamic viscosity of water at a given temperature.
 
         If given units, the function will automatically convert to Kelvin.
         If not given units, the function will assume Kelvin.
         """
         if coag_type.lower() == 'alum':
-            coag_nu = self._alum_nu(coag_conc)
+            coag_nu = self.alum_nu(coag_conc)
         elif coag_type.lower() == 'pacl':
-            coag_nu = self._pacl_nu(coag_conc)
+            coag_nu = self.pacl_nu(coag_conc)
         return coag_nu
 
     @property
     def coag_q_max_est(self):
+        """The estimated maximum permissible flow rate of the coagulant stock,
+        based on a whole number of sacks of coagulant used.
+
+        .. deprecated::
+            `coag_q_max_est` is deprecated; use `coag_q_max` instead, which is
+            based on the exact user-defined coagulant stock concentration, 
+            `coag_stock_conc`.
+        """
+        # Deprecated since January 2021
+        warnings.warn('coag_q_max_est is deprecated; use coag_q_max instead,\
+            which is based on the exact user-defined coagulant stock \
+            concentration, coag_stock_conc.', UserWarning)
+
         coag_q_max_est = self.q * self.coag_dose_conc_max / \
             self.coag_stock_conc_est
         return coag_q_max_est
 
+    @property
+    def coag_q_max(self):
+        """The maximum permissible flow rate of the coagulant stock."""
+        coag_q_max = self.q * self.coag_dose_conc_max / self.coag_stock_conc
+        return coag_q_max.to(u.L / u.s)
+
     @property
     def coag_stock_vol(self):
+        """The volume of the coagulant stock tank, based on available sizes
+        from the supplier.
+        """
         coag_stock_vol = ut.ceil_nearest(
                 self.coag_stock_min_est_time * self.train_n *
-                    self.coag_q_max_est,
+                    self.coag_q_max,
                 self.chem_tank_vol_supplier
             )
         return coag_stock_vol
 
     @property
     def coag_sack_n(self):
+        """The number of sacks of coagulant needed to make the coagulant stock,
+        rounded to the nearest whole number.
+        """
         coag_sack_n = round(
-                (self.coag_stock_vol * self.coag_stock_conc_est /
+                (self.coag_stock_vol * self.coag_stock_conc /
                 self.coag_sack_mass).to_base_units()
             )
         return coag_sack_n
 
-    @property
-    def coag_stock_conc(self):
-        coag_stock_conc = self.coag_sack_n * self.coag_sack_mass / \
-            self.coag_stock_vol
-        return coag_stock_conc
-
-    @property
-    def coag_q_max(self):
-        coag_q_max = self.q * self.coag_dose_conc_max / self.coag_stock_conc
-        return coag_q_max.to(u.L / u.s)
+    # Commented out January 2021
+    # @property
+    # def coag_stock_conc(self):
+    #     """The concentration of the coagulant stock."""
+    #     coag_stock_conc = self.coag_sack_n * self.coag_sack_mass / \
+    #         self.coag_stock_vol
+    #     return coag_stock_conc
 
     @property
     def coag_stock_time_min(self):
+        """The minimum amount of time that the coagulant stock will last."""
         return self.coag_stock_vol / (self.train_n * self.coag_q_max)
 
     @property
     def coag_stock_nu(self):
-        return self._coag_nu(self.coag_stock_conc, self.coag_type)
+        """The kinematic viscosity of the coagulant stock."""
+        return self.coag_nu(self.coag_stock_conc, self.coag_type)
 #==============================================================================
 # Small-diameter Tube Design
 #==============================================================================
     @property
     def _coag_tube_q_max(self):
         """The maximum permissible flow through a coagulant tube."""
         coag_tube_q_max = ((np.pi * self.coag_tube_id ** 2)/4) * \
-            np.sqrt((2 * self.error_ratio * self.hl * con.GRAVITY)/self.tube_k)
-        return coag_tube_q_max
+            np.sqrt((2 * self.error_ratio * self.hl * u.gravity)/self.tube_k)
+        return coag_tube_q_max.to(u.L / u.s)
 
     @property
     def coag_tubes_active_n(self):
+        """The number of coagulant tubes in use."""
         coag_tubes_active_n = \
             np.ceil((self.coag_q_max / self._coag_tube_q_max).to_base_units())
         return coag_tubes_active_n
 
     @property
     def coag_tubes_n(self):
+        """The number of coagulant tubes in use, plus a spare tube for 
+        maintenance.
+        """
         coag_tubes_n = self.coag_tubes_active_n + 1
         return coag_tubes_n
 
@@ -158,8 +221,9 @@ def coag_tube_operating_q_max(self):
 
     @property
     def coag_tube_l(self):
+        """The length of a coagulant tube."""
         coag_tube_l = (
-                self.hl * con.GRAVITY * np.pi * self.coag_tube_id ** 4 /
+                self.hl * u.gravity * np.pi * self.coag_tube_id ** 4 /
                 (128 * self.coag_stock_nu * self.coag_tube_operating_q_max)
             ) - (
                 8 * self.coag_tube_operating_q_max * self.tube_k /
@@ -169,14 +233,18 @@ def coag_tube_l(self):
 
     @property
     def coag_tank_r(self):
-        index = np.where(self.chem_tank_vol_supplier == self.coag_stock_vol)
-        coag_tank_r = self.chem_tank_dimensions_supplier[0][index] / 2
+        """The radius of the coagulant stock tank, based on available sizes
+        from the supplier."""
+        index = np.argmax(self.chem_tank_vol_supplier == self.coag_stock_vol)
+        coag_tank_r = self.chem_tank_dimensions_supplier[index][0] / 2
         return coag_tank_r
 
     @property
     def coag_tank_h(self):
-        index = np.where(self.chem_tank_vol_supplier == self.coag_stock_vol)
-        coag_tank_h = self.chem_tank_dimensions_supplier[1][index]
+        """The height of the coagulant stock tank, based on available sizes
+        from the supplier."""
+        index = np.argmax(self.chem_tank_vol_supplier == self.coag_stock_vol)
+        coag_tank_h = self.chem_tank_dimensions_supplier[index][1]
         return coag_tank_h
 
     def _DiamTubeAvail(self, en_tube_series = True):

aguaclara/design/pipeline.py

@@ -348,7 +348,7 @@ def ID_SDR_all_available(self, SDR):
     @property
     def headloss(self):
         """Return the total head loss from major and minor losses in a pipe."""
-        return pc.headloss_fric(
+        return pc.headloss_major_pipe(
                 self.q, self.id, self.l, self.nu, self.pipe_rough
             ).to(u.cm)
 
@@ -446,7 +446,7 @@ def _get_id(self, size):
     @property
     def headloss(self):
         """The headloss"""
-        return pc.elbow_minor_loss(self.q, self.id, self.k_minor).to(u.cm)
+        return pc.headloss_minor_elbow(self.q, self.id, self.k_minor).to(u.cm)
 
     def format_print(self):
         """The string representation for an Elbow Fitting."""
@@ -553,11 +553,11 @@ def _set_next(self):
 
     def _headloss_left(self):
         """The headloss of the left outlet"""
-        return pc.elbow_minor_loss(self.q, self.id, self.left_k_minor).to(u.cm)
+        return pc.headloss_minor_elbow(self.q, self.id, self.left_k_minor).to(u.cm)
 
     def _headloss_right(self):
         """The headloss of the right outlet"""
-        return pc.elbow_minor_loss(self.q, self.id, self.right_k_minor).to(u.cm)
+        return pc.headloss_minor_elbow(self.q, self.id, self.right_k_minor).to(u.cm)
 
     @property
     def headloss(self):

aguaclara/research/environmental_processes_analysis.py

@@ -293,7 +293,7 @@ def E_CMFR_N(t, N):
     """
     return (N**N)/special.gamma(N) * (t**(N-1))*np.exp(-N*t)
 
-
+@ut.list_handler()
 def E_Advective_Dispersion(t, Pe):
     """Calculate a dimensionless measure of the output tracer concentration from
     a spike input to reactor with advection and dispersion.
@@ -312,11 +312,10 @@ def E_Advective_Dispersion(t, Pe):
     >>> round(E_Advective_Dispersion(0.5, 5), 7)
     0.4774864
     """
-    # replace any times at zero with a number VERY close to zero to avoid
-    # divide by zero errors
-    if isinstance(t, list):
-        t[t == 0] = 10**(-10)
-    return (Pe/(4*np.pi*t))**(0.5)*np.exp((-Pe*((1-t)**2))/(4*t))
+    if t == 0:
+        return 0
+    else:
+        return (Pe/(4*np.pi*t))**(0.5)*np.exp((-Pe*((1-t)**2))/(4*t))
 
 
 def Tracer_CMFR_N(t_seconds, t_bar, C_bar, N):

docs/source/design/cdc.rst

@@ -0,0 +1,7 @@
+.. _design-cdc:
+
+Chemical Dose Controller
+========================
+
+.. automodule:: aguaclara.design.cdc
+    :members:

docs/source/design/design.rst

@@ -4,6 +4,7 @@ Design
 .. toctree::
     :maxdepth: 2
 
+    cdc
     component
     ent_floc
     ent