fixed multiple errors in the LFOM code. (#138)

* fixed multiple errors in the code. Confirmed that it produces reasonable designs over a range of flow rates.
Added all inputs as keyword inputs.

* Fix unit and test discrepancies

* Increment version number

aguaclara/design/lfom.py

@@ -5,33 +5,30 @@
 import aguaclara.core.utility as ut
 import aguaclara.core.drills as drills
 from aguaclara.core.units import unit_registry as u
-
-# from onshapepy.part import Part
-
 import numpy as np
 import math
 
 
 class LFOM:
 
-    safety_factor = 1.5
-    sdr = 26
-    drill_bits = drills.DRILL_BITS_D_IMPERIAL
-    s_orfice = 1*u.cm
-    # cad = Part(
-    #     "https://cad.onshape.com/documents/e1798ab5f546e1414e86992d/w/104d463fef6c6a71c703abe6/e/890edb42c7884277d8d8711d"
-    # )
+    def __init__(self, q=20*u.L/u.s, hl=20*u.cm, safety_factor=1.5, sdr=26,
+                 drill_bits=drills.DRILL_BITS_D_IMPERIAL, s_orifice=0.5*u.cm):
 
-    def __init__(self, q=20*u.L/u.s, hl=20*u.cm):
         self.q = q
         self.hl = hl
+        self.safety_factor = safety_factor
+        self.sdr = sdr
+        self.drill_bits = drill_bits
+        self.s_orifice = s_orifice
 
-    def width_stout(self, z):
+    def stout_w_per_flow(self, z):
         """Return the width of a Stout weir at elevation z. More info
-        here. <https://confluence.cornell.edu/display/AGUACLARA/LFOM+sutro+weir+research>
+        here. <https://confluence.cornell.edu/display/AGUACLARA/
+        LFOM+sutro+weir+research>
         """
-        w_per_flow = 2 / ((2 * pc.gravity * z) ** (1 / 2) * con.VC_ORIFICE_RATIO * np.pi * self.hl)
-        return w_per_flow*self.q
+        w_per_flow = 2 / ((2 * pc.gravity * z) ** (1 / 2) *
+                          con.VC_ORIFICE_RATIO * np.pi * self.hl)
+        return w_per_flow.to_base_units()
 
     @property
     def n_rows(self):
@@ -44,12 +41,11 @@ def n_rows(self):
         of the orifices. This spacing function also sets the lower depth on
         the high flow rate LFOM with no accurate flows below a depth equal
         to the first row height.
-
         But it might be better to always set then number of rows to 10.
         The challenge is to figure out a reasonable system of constraints that
         reliably returns a valid solution.
         """
-        N_estimated = (self.hl * np.pi / (2 * self.width_stout(self.hl) * self.q))
+        N_estimated = (self.hl * np.pi / (2 * self.stout_w_per_flow(self.hl) * self.q)).to(u.dimensionless)
         variablerow = min(10, max(4, math.trunc(N_estimated.magnitude)))
         return variablerow
 
@@ -59,19 +55,19 @@ def b_rows(self):
         Message how long it took to load everything (minus packages)"""
         return self.hl / self.n_rows
 
-
     @property
     def vel_critical(self):
         """The average vertical velocity of the water inside the LFOM pipe
         at the very bottom of the bottom row of orifices The speed of
         falling water is 0.841 m/s for all linear flow orifice meters of
         height 20 cm, independent of total plant flow rate. """
-        return 4 / (3 * math.pi) * (2 * pc.gravity * self.hl) ** (1 / 2)
+        return (4 / (3 * math.pi) * (2 * pc.gravity * self.hl) ** (1 / 2)).to(u.m/u.s)
 
     @property
     def area_pipe_min(self):
-        """The minimum cross-sectional area of the LFOM pipe that assures a safety factor."""
-        return self.safety_factor * self.q / self.vel_critical
+        """The minimum cross-sectional area of the LFOM pipe that assures
+        a safety factor."""
+        return (self.safety_factor * self.q / self.vel_critical).to(u.cm**2)
 
     @property
     def nom_diam_pipe(self):
@@ -89,7 +85,7 @@ def area_top_orifice(self):
         # Calculate the center of the top row:
         z = self.hl - 0.5 * self.b_rows
         # Multiply the stout weir width by the height of one row.
-        return self.width_stout(z) * self.b_rows
+        return self.stout_w_per_flow(z) * self.q * self.b_rows
 
     @property
     def d_orifice_max(self):
@@ -98,7 +94,8 @@ def d_orifice_max(self):
 
     @property
     def orifice_diameter(self):
-        """The actual orifice diameter. We don't let the diameter extend beyond its row space. """
+        """The actual orifice diameter. We don't let the diameter extend
+        beyond its row space. """
         maxdrill = min(self.b_rows, self.d_orifice_max)
         return ut.floor_nearest(maxdrill, self.drill_bits)
 
@@ -114,14 +111,14 @@ def n_orifices_per_row_max(self):
         structural integrity of the pipe.
         """
         c = math.pi * pipe.ID_SDR(self.nom_diam_pipe, self.sdr)
-        b = self.orifice_diameter + self.s_orfice
+        b = self.orifice_diameter + self.s_orifice
 
         return math.floor(c/b)
 
     @property
     def flow_ramp(self):
         """An equally spaced array representing flow at each row."""
-        return np.linspace(self.q / self.n_rows, self.q, self.n_rows)*self.q.units
+        return np.linspace(1 / self.n_rows, 1, self.n_rows)*self.q
 
     @property
     def height_orifices(self):
@@ -130,58 +127,63 @@ def height_orifices(self):
         point of the LFOM so that the flow goes to zero when the water height
         is at zero.
         """
-        return np.arange((self.orifice_diameter * 0.5), self.hl, self.q)
+
+        return (np.linspace(0, self.n_rows-1, self.n_rows))*self.b_rows + 0.5 * self.orifice_diameter
 
     def flow_actual(self, Row_Index_Submerged, N_LFOM_Orifices):
         """Calculates the flow for a given number of submerged rows of orifices
-        harray is the distance from the water level to the center of the orifices
-        when the water is at the max level.
-
+        harray is the distance from the water level to the center of the
+        orifices when the water is at the max level.
         Parameters
         ----------
         Row_Index_Submerged: int
-            The index of the submerged row. All rows below and including this index are submerged.
+          The index of the submerged row. All rows below and including this
+          index are submerged.
         N_LFOM_Orifices: [int]
-            The number of orifices at each row.
-
+          The number of orifices at each row.
         Returns
         --------
         The flow through all of the orifices that are submerged.
-
         """
-        harray = (np.linspace(self.b_rows, self.hl, self.n_rows))*self.hl.units - 0.5 * self.orifice_diameter
-        FLOW_new = 0
+
+        flow = 0
         for i in range(Row_Index_Submerged + 1):
-            FLOW_new = FLOW_new + (N_LFOM_Orifices[i] * (
-                pc.flow_orifice_vert(self.orifice_diameter, harray[Row_Index_Submerged - i],
-                                     con.VC_ORIFICE_RATIO)))
-        return FLOW_new
+            flow = flow + (N_LFOM_Orifices[i] * (
+               pc.flow_orifice_vert(self.orifice_diameter,
+                                    self.b_rows*(Row_Index_Submerged + 1)
+                                    - self.height_orifices[i],
+                                    con.VC_ORIFICE_RATIO)))
+        return flow
 
     @property
     def n_orifices_per_row(self):
-        """Calculate number of orifices at each level given a diameter.
+        """Calculate number of orifices at each level given an orifice
+        diameter.
         """
-        # H is distance from the elevation between two rows of orifices down to the bottom of the orifices
-        H = self.b_rows/2 + 0.5*self.orifice_diameter
+        # H is distance from the bottom of the next row of orifices to the
+        # center of the current row of orifices
+        H = self.b_rows - 0.5*self.orifice_diameter
         flow_per_orifice = pc.flow_orifice_vert(self.orifice_diameter, H, con.VC_ORIFICE_RATIO)
         n = np.zeros(self.n_rows)
         for i in range(self.n_rows):
-            # calculate the ideal number of orifices at the current row without constraining to an integer
+            # calculate the ideal number of orifices at the current row without
+            # constraining to an integer
             flow_needed = self.flow_ramp[i] - self.flow_actual(i, n)
             n_orifices_real = (flow_needed / flow_per_orifice).to(u.dimensionless)
-            # constrain number of orifices to be less than the max  per row and greater or equal to 0
+            # constrain number of orifices to be less than the max per row and
+            # greater or equal to 0
             n[i] = min((max(0, round(n_orifices_real))), self.n_orifices_per_row_max)
         return n
 
     @property
     def error_per_row(self):
-        """This function calculates the error of the design based on the differences between the predicted flow rate
+        """This function calculates the error of the design based on the
+        differences between the predicted flow rate
         and the actual flow rate through the LFOM."""
-        FLOW_lfom_error = []
-        for i in range(self.n_rows-1):
+        FLOW_lfom_error = np.zeros(self.n_rows)
+        for i in range(self.n_rows):
             actual_flow = self.flow_actual(i, self.n_orifices_per_row)
-            row_error = (actual_flow - self.flow_ramp[i]) / self.q
-            FLOW_lfom_error.append(row_error.to(u.dimensionless))
+            FLOW_lfom_error[i] = (((actual_flow - self.flow_ramp[i]) / self.flow_ramp[i]).to(u.dimensionless)).magnitude
         return FLOW_lfom_error
 
     # def draw(self):

setup.py

@@ -1,7 +1,7 @@
 from setuptools import setup, find_packages
 
 setup(name='aguaclara',
-      version='0.0.18',
+      version='0.0.19',
       description='Open source functions for AguaClara water treatment research and plant design.',
       url='https://github.com/AguaClara/aguaclara',
       author='AguaClara at Cornell',

tests/design/test_lfom.py

@@ -12,7 +12,7 @@ def lfom():
 
 
 def test_lfom(lfom):
-    assert lfom.width_stout(10*u.cm).to(u.m) == 0.002405154520928927*u.m
+    assert lfom.stout_w_per_flow(10 * u.cm) == 2.4051545209289267 * u.s / u.m ** 2
     assert lfom.n_rows == 10
     assert lfom.b_rows.to(u.m) == 0.03*u.m
     assert lfom.vel_critical.to(u.m/u.s) == 1.0294963872061624 * u.m/u.s
@@ -22,8 +22,8 @@ def test_lfom(lfom):
     assert lfom.d_orifice_max.to(u.m) == 0.00737693510978969 * u.m
     assert lfom.orifice_diameter.to(u.m) == 0.00635*u.m
     assert lfom.drillbit_area.to(u.m**2) == 3.1669217443593606e-05*u.m**2
-    assert lfom.n_orifices_per_row_max == 10
+    assert lfom.n_orifices_per_row_max == 15
     assert lfom.flow_ramp[5] == 0.6 * u.L/u.s
     assert lfom.flow_actual(2, [4, 3, 2]) == 0.0001962543726011661 * u.m**3/u.s
-    assert lfom.n_orifices_per_row[0] == 8
-    assert -0.001 < (np.average(lfom.error_per_row) - 0.005194582036259183) < 0.001
+    assert lfom.n_orifices_per_row[0] == 7
+    assert -0.01 < (np.average(lfom.error_per_row) - 0.005194582036259183) < 0.01