Repository: ross-rotordynamics/ross
Pull request #617 (2020-07-06 23:35:14)

Reduce attributes
This pull request is related to issue #616 

@rodrigomoliveira1 , I tried to reduce as many attributes as possible. Some could be removed, but are used many times in the code and calculating them several times would increase the time it takes to run the code.

Do you think we should modify anything else?


Patch:
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index e3bb5265..8371e3cc 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -32,8 +32,6 @@ class FluidFlow:
         Number of points along the Z direction (direction of flow).
     ntheta: int
         Number of points along the direction theta. NOTE: ntheta must be odd.
-    nradius: int
-        Number of points along the direction r.
     length: float
         Length in the Z direction (m).
 
@@ -91,13 +89,11 @@ class FluidFlow:
     dtheta: float
         Range size in the theta direction.
     ntotal: int
-        Number of nodes in the grid., ntheta, nradius, n_interv_z, n_interv_theta,
+        Number of nodes in the grid., ntheta, n_interv_z, n_interv_theta,
     n_interv_z: int
         Number of intervals on Z.
     n_interv_theta: int
       Number of intervals on theta.
-    n_interv_radius: int
-        Number of intervals on r.
     p_mat_analytical : array of shape (nz, ntheta)
         The analytical pressure matrix.
     p_mat_numerical : array of shape (nz, ntheta)
@@ -114,8 +110,6 @@ class FluidFlow:
         The external radius in each position of the grid.
     ri : array of shape (nz, ntheta)
         The internal radius in each position of the grid.
-    z_list : array of shape (1, nz)
-        z along the object. It goes from 0 to lb.
     xre : array of shape (nz, ntheta)
         x value of the external radius.
     xri : array of shape (nz, ntheta)
@@ -124,6 +118,11 @@ class FluidFlow:
         y value of the external radius.
     yri : array of shape (nz, ntheta)
         y value of the internal radius.
+    z_list : array of shape (1, nz)
+        z along the object. It goes from 0 to lb.
+    gama : array of shape (nz, ntheta)
+        Points along the object in the tangential direction.
+        It ranges from 0 to 2 pi, starting at the largest spacing between rotor and stator.
     t : float
         Time.
     xp : float
@@ -156,7 +155,6 @@ class FluidFlow:
     >>> import numpy as np
     >>> nz = 8
     >>> ntheta = 64
-    >>> nradius = 8
     >>> length = 0.01
     >>> omega = 100.*2*np.pi/60
     >>> p_in = 0.
@@ -167,7 +165,7 @@ class FluidFlow:
     >>> density = 860.
     >>> eccentricity = 0.001
     >>> attitude_angle = np.pi
-    >>> my_fluid_flow = flow.FluidFlow(nz, ntheta, nradius, length,
+    >>> my_fluid_flow = flow.FluidFlow(nz, ntheta, length,
     ...                                omega, p_in, p_out, radius_rotor,
     ...                                radius_stator, viscosity, density,
     ...                                attitude_angle=attitude_angle, eccentricity=eccentricity,
@@ -187,7 +185,6 @@ class FluidFlow:
         self,
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
@@ -206,10 +203,8 @@ class FluidFlow:
             sys.exit("Either load or eccentricity must be given.")
         self.nz = nz
         self.ntheta = ntheta
-        self.nradius = nradius
         self.n_interv_z = nz - 1
         self.n_interv_theta = ntheta - 1
-        self.n_interv_radius = nradius - 1
         self.length = length
         self.ltheta = 2.0 * np.pi
         self.dz = length / self.n_interv_z
@@ -269,21 +264,18 @@ class FluidFlow:
         self.yi = self.eccentricity * np.sin(3 * np.pi / 2 + self.attitude_angle)
         self.re = np.zeros([self.nz, self.ntheta])
         self.ri = np.zeros([self.nz, self.ntheta])
-        self.z_list = np.zeros(self.nz)
         self.xre = np.zeros([self.nz, self.ntheta])
         self.xri = np.zeros([self.nz, self.ntheta])
         self.yre = np.zeros([self.nz, self.ntheta])
         self.yri = np.zeros([self.nz, self.ntheta])
+        self.z_list = np.zeros(self.nz)
+        self.gama = np.zeros([self.nz, self.ntheta])
+        self.t = 0
+        self.xp = 0
+        self.yp = 0
         self.p_mat_analytical = np.zeros([self.nz, self.ntheta])
-        self.c1 = np.zeros([self.nz, self.ntheta])
-        self.c2 = np.zeros([self.nz, self.ntheta])
-        self.c0w = np.zeros([self.nz, self.ntheta])
-        self.M = np.zeros([self.ntotal, self.ntotal])
-        self.f = np.zeros([self.ntotal, 1])
-        self.P = np.zeros([self.ntotal, 1])
         self.p_mat_numerical = np.zeros([self.nz, self.ntheta])
-        self.gama = np.zeros([self.nz, self.ntheta])
-        self.calculate_coefficients()
+        self.geometry_description()
         self.analytical_pressure_matrix_available = False
         self.numerical_pressure_matrix_available = False
         self.calculate_pressure_matrix_numerical()
@@ -319,19 +311,19 @@ class FluidFlow:
                 for i in range(0, self.nz):
                     for j in range(0, self.ntheta):
                         # fmt: off
-                        self.p_mat_analytical[i][j] = (
+                        self.p_mat_analytical[i, j] = (
                                 ((-3 * self.viscosity * self.omega) / self.radial_clearance ** 2) *
                                 ((i * self.dz - (self.length / 2)) ** 2 - (self.length ** 2) / 4) *
                                 (self.eccentricity_ratio * np.sin(j * self.dtheta)) /
                                 (1 + self.eccentricity_ratio * np.cos(j * self.dtheta)) ** 3)
                         # fmt: on
-                        if self.p_mat_analytical[i][j] < 0:
-                            self.p_mat_analytical[i][j] = 0
+                        if self.p_mat_analytical[i, j] < 0:
+                            self.p_mat_analytical[i, j] = 0
             elif method == 1:
                 for i in range(0, self.nz):
                     for j in range(0, self.ntheta):
                         # fmt: off
-                        self.p_mat_analytical[i][j] = (3 * self.viscosity / ((self.radial_clearance ** 2) *
+                        self.p_mat_analytical[i, j] = (3 * self.viscosity / ((self.radial_clearance ** 2) *
                                                                              (1. + self.eccentricity_ratio * np.cos(
                                                                                  j * self.dtheta)) ** 3)) * \
                                                       (-self.eccentricity_ratio * self.omega * np.sin(
@@ -339,18 +331,18 @@ class FluidFlow:
                                                       (((i * self.dz - (self.length / 2)) ** 2) - (
                                                               self.length ** 2) / 4)
                         # fmt: on
-                        if self.p_mat_analytical[i][j] < 0:
-                            self.p_mat_analytical[i][j] = 0
+                        if self.p_mat_analytical[i, j] < 0:
+                            self.p_mat_analytical[i, j] = 0
         elif self.bearing_type == "long_bearing" or force_type == "long":
             if method == 0:
                 for i in range(0, self.nz):
                     for j in range(0, self.ntheta):
-                        self.p_mat_analytical[i][j] = (
+                        self.p_mat_analytical[i, j] = (
                             (
                                 6
                                 * self.viscosity
                                 * self.omega
-                                * (self.ri[i][j] / self.radial_clearance) ** 2
+                                * (self.radius_rotor / self.radial_clearance) ** 2
                                 * self.eccentricity_ratio
                                 * np.sin(self.dtheta * j)
                                 * (
@@ -368,8 +360,8 @@ class FluidFlow:
                             )
                             + self.p_in
                         )
-                        if self.p_mat_analytical[i][j] < 0:
-                            self.p_mat_analytical[i][j] = 0
+                        if self.p_mat_analytical[i, j] < 0:
+                            self.p_mat_analytical[i, j] = 0
         elif self.bearing_type == "medium_size":
             raise ValueError(
                 "The pressure matrix for a bearing that is neither short or long can only be calculated "
@@ -379,158 +371,209 @@ class FluidFlow:
         self.analytical_pressure_matrix_available = True
         return self.p_mat_analytical
 
-    def calculate_coefficients(self, direction=None):
-        """This function calculates the constants that form the Poisson equation
-        of the discrete pressure (central differences in the second
-        derivatives). It is executed when the class is instantiated.
+    def geometry_description(self):
+        """This function calculates the geometry description.
+        It is executed when the class is instantiated.
         Examples
         --------
         >>> my_fluid_flow = fluid_flow_example()
-        >>> my_fluid_flow.calculate_coefficients()
-        >>> my_fluid_flow.c0w # doctest: +ELLIPSIS
-        array([[...
+        >>> my_fluid_flow.geometry_description()
         """
         for i in range(0, self.nz):
             zno = i * self.dz
             self.z_list[i] = zno
-            eccentricity_error = False
             for j in range(0, self.ntheta):
                 # fmt: off
-                self.gama[i][j] = j * self.dtheta + np.pi / 2 + self.attitude_angle
-                [radius_external, self.xre[i][j], self.yre[i][j]] = \
-                    external_radius_function(self.gama[i][j], self.radius_stator)
-                [radius_internal, self.xri[i][j], self.yri[i][j]] = \
-                    internal_radius_function(self.gama[i][j], self.attitude_angle, self.radius_rotor,
+                self.gama[i, j] = j * self.dtheta + np.pi / 2 + self.attitude_angle
+                [radius_external, self.xre[i, j], self.yre[i, j]] = \
+                    external_radius_function(self.gama[i, j], self.radius_stator)
+                [radius_internal, self.xri[i, j], self.yri[i, j]] = \
+                    internal_radius_function(self.gama[i, j], self.attitude_angle, self.radius_rotor,
                                              self.eccentricity)
-                self.re[i][j] = radius_external
-                self.ri[i][j] = radius_internal
+                self.re[i, j] = radius_external
+                self.ri[i, j] = radius_internal
+                # fmt: on
+
+    def calculate_coefficients(self, direction=None):
+        """This function calculates the constants that form the Poisson equation
+        of the discrete pressure (central differences in the second
+        derivatives).
+        Parameters
+        ----------
+        direction: str
+            If defined, it calculates the model based on the disturbance in the chosen direction: 'x' or 'y'.
+        Returns
+        --------
+        c1, c2, c0w: matrix of float
+            Constants that form the Poisson equation.
+        >>> my_fluid_flow = fluid_flow_example()
+        >>> my_fluid_flow.calculate_coefficients()# doctest: +ELLIPSIS
+        (array([[...
+        """
+        c1 = np.zeros([self.nz, self.ntheta])
+        c2 = np.zeros([self.nz, self.ntheta])
+        c0w = np.zeros([self.nz, self.ntheta])
+
+        for i in range(0, self.nz):
+            eccentricity_error = False
+            for j in range(0, self.ntheta):
+                # fmt: off
 
                 w = self.omega * self.radius_rotor
 
-                k = (self.re[i][j] ** 2 * (np.log(self.re[i][j]) - 1 / 2) - self.ri[i][j] ** 2 *
-                     (np.log(self.ri[i][j]) - 1 / 2)) / (self.ri[i][j] ** 2 - self.re[i][j] ** 2)
-
-                self.c1[i][j] = (1 / (4 * self.viscosity)) * ((self.re[i][j] ** 2 * np.log(self.re[i][j]) -
-                                                               self.ri[i][j] ** 2 * np.log(self.ri[i][j]) +
-                                                               (self.re[i][j] ** 2 - self.ri[i][j] ** 2) *
-                                                               (k - 1)) - 2 * self.re[i][j] ** 2 * (
-                                                                      (np.log(self.re[i][j]) + k - 1 / 2) * np.log(
-                                                                       self.re[i][j] / self.ri[i][j])))
-
-                self.c2[i][j] = (- self.ri[i][j] ** 2) / (8 * self.viscosity) * \
-                                ((self.re[i][j] ** 2 - self.ri[i][j] ** 2 -
-                                  (self.re[i][j] ** 4 - self.ri[i][j] ** 4) /
-                                  (2 * self.ri[i][j] ** 2)) +
-                                 ((self.re[i][j] ** 2 - self.ri[i][j] ** 2) /
-                                  (self.ri[i][j] ** 2 *
-                                   np.log(self.re[i][j] / self.ri[i][j]))) *
-                                 (self.re[i][j] ** 2 * np.log(self.re[i][j] / self.ri[i][j]) -
-                                  (self.re[i][j] ** 2 - self.ri[i][j] ** 2) / 2))
-
-                self.c0w[i][j] = (- w * self.ri[i][j] *
-                                  (np.log(self.re[i][j] / self.ri[i][j]) *
-                                   (1 + (self.ri[i][j] ** 2) / (self.re[i][j] ** 2 - self.ri[i][j] ** 2)) - 1 / 2))
+                k = (self.re[i, j] ** 2 * (np.log(self.re[i, j]) - 1 / 2) - self.ri[i, j] ** 2 *
+                     (np.log(self.ri[i, j]) - 1 / 2)) / (self.ri[i, j] ** 2 - self.re[i, j] ** 2)
+
+                c1[i, j] = (1 / (4 * self.viscosity)) * ((self.re[i, j] ** 2 * np.log(self.re[i, j]) -
+                                                               self.ri[i, j] ** 2 * np.log(self.ri[i, j]) +
+                                                               (self.re[i, j] ** 2 - self.ri[i, j] ** 2) *
+                                                               (k - 1)) - 2 * self.re[i, j] ** 2 * (
+                                                                      (np.log(self.re[i, j]) + k - 1 / 2) * np.log(
+                                                                       self.re[i, j] / self.ri[i, j])))
+
+                c2[i, j] = (- self.ri[i, j] ** 2) / (8 * self.viscosity) * \
+                                ((self.re[i, j] ** 2 - self.ri[i, j] ** 2 -
+                                  (self.re[i, j] ** 4 - self.ri[i, j] ** 4) /
+                                  (2 * self.ri[i, j] ** 2)) +
+                                 ((self.re[i, j] ** 2 - self.ri[i, j] ** 2) /
+                                  (self.ri[i, j] ** 2 *
+                                   np.log(self.re[i, j] / self.ri[i, j]))) *
+                                 (self.re[i, j] ** 2 * np.log(self.re[i, j] / self.ri[i, j]) -
+                                  (self.re[i, j] ** 2 - self.ri[i, j] ** 2) / 2))
+
+                c0w[i, j] = (- w * self.ri[i, j] *
+                                  (np.log(self.re[i, j] / self.ri[i, j]) *
+                                   (1 + (self.ri[i, j] ** 2) / (self.re[i, j] ** 2 - self.ri[i, j] ** 2)) - 1 / 2))
                 if direction == "x":
-                    a = self.omegap * (self.xp) * np.cos(self.omegap * self.t)
-                    self.c0w[i][j] += self.ri[i][j] * a * np.sin(self.gama[i][j])
+                    a = self.omegap * self.xp * np.cos(self.omegap * self.t)
+                    c0w[i, j] += self.ri[i, j] * a * np.sin(self.gama[i, j])
                 elif direction == "y":
-                    b = self.omegap * (self.yp) * np.cos(self.omegap * self.t)
-                    self.c0w[i][j] -= self.ri[i][j] * b * np.cos(self.gama[i][j])
+                    b = self.omegap * self.yp * np.cos(self.omegap * self.t)
+                    c0w[i, j] -= self.ri[i, j] * b * np.cos(self.gama[i, j])
                 else:
-                    self.c0w[i][j] += 0
+                    c0w[i, j] += 0
                 # fmt: on
                 if not eccentricity_error:
-                    if abs(self.xri[i][j]) > abs(self.xre[i][j]) or abs(
-                        self.yri[i][j]
-                    ) > abs(self.yre[i][j]):
+                    if abs(self.xri[i, j]) > abs(self.xre[i, j]) or abs(
+                        self.yri[i, j]
+                    ) > abs(self.yre[i, j]):
                         eccentricity_error = True
             if eccentricity_error:
                 raise ValueError(
                     "Error: The given parameters create a rotor that is not inside the stator. "
                     "Check parameters and fix accordingly."
                 )
+        return c1, c2, c0w
 
-    def mounting_matrix(self):
+    def mounting_matrix(self, c1, c2, c0w):
         """This function assembles the matrix M and the independent vector f.
+        Parameters
+        ----------
+        c1, c2, c0w: matrix of float
+            Constants that form the Poisson equation.
+        Returns
+        --------
+        M: matrix of float
+            Matrix composed of coefficients that multiply the pressures at each point in the discrete domain.
+        f: array of float
+            Pressure independent terms.
         Examples
         --------
         >>> my_fluid_flow = fluid_flow_example()
-        >>> my_fluid_flow.mounting_matrix()
-        >>> my_fluid_flow.M # doctest: +ELLIPSIS
-        array([[...
+        >>> c1, c2, c0w = my_fluid_flow.calculate_coefficients()
+        >>> my_fluid_flow.mounting_matrix(c1, c2, c0w)# doctest: +ELLIPSIS
+        (array([[...
         """
         # fmt: off
+        M = np.zeros([self.ntotal, self.ntotal])
+        f = np.zeros([self.ntotal, 1])
         count = 0
         for x in range(self.ntheta):
-            self.M[count][count] = 1
-            self.f[count][0] = self.p_in
+            M[count, count] = 1
+            f[count, 0] = self.p_in
             count = count + self.nz - 1
-            self.M[count][count] = 1
-            self.f[count][0] = self.p_out
+            M[count, count] = 1
+            f[count, 0] = self.p_out
             count = count + 1
         count = 0
         for x in range(self.nz - 2):
-            self.M[self.ntotal - self.nz + 1 + count][1 + count] = 1
-            self.M[self.ntotal - self.nz + 1 + count][self.ntotal - self.nz + 1 + count] = -1
+            M[self.ntotal - self.nz + 1 + count, 1 + count] = 1
+            M[self.ntotal - self.nz + 1 + count, self.ntotal - self.nz + 1 + count] = -1
             count = count + 1
         count = 1
         j = 0
         for i in range(1, self.nz - 1):
-            a = (1 / self.dtheta ** 2) * (self.c1[i][self.ntheta - 1])
-            self.M[count][self.ntotal - 2 * self.nz + count] = a
-            b = (1 / self.dz ** 2) * (self.c2[i - 1, j])
-            self.M[count][count - 1] = b
-            c = -((1 / self.dtheta ** 2) * ((self.c1[i][j]) + self.c1[i][self.ntheta - 1])
-                  + (1 / self.dz ** 2) * (self.c2[i][j] + self.c2[i - 1][j]))
-            self.M[count, count] = c
-            d = (1 / self.dz ** 2) * (self.c2[i][j])
-            self.M[count][count + 1] = d
-            e = (1 / self.dtheta ** 2) * (self.c1[i][j])
-            self.M[count][count + self.nz] = e
+            a = (1 / self.dtheta ** 2) * (c1[i, self.ntheta - 1])
+            M[count, self.ntotal - 2 * self.nz + count] = a
+            b = (1 / self.dz ** 2) * (c2[i - 1, j])
+            M[count, count - 1] = b
+            c = -((1 / self.dtheta ** 2) * ((c1[i, j]) + c1[i, self.ntheta - 1])
+                  + (1 / self.dz ** 2) * (c2[i, j] + c2[i - 1, j]))
+            M[count, count] = c
+            d = (1 / self.dz ** 2) * (c2[i, j])
+            M[count, count + 1] = d
+            e = (1 / self.dtheta ** 2) * (c1[i, j])
+            M[count, count + self.nz] = e
             count = count + 1
         count = self.nz + 1
         for j in range(1, self.ntheta - 1):
             for i in range(1, self.nz - 1):
-                a = (1 / self.dtheta ** 2) * (self.c1[i, j - 1])
-                self.M[count][count - self.nz] = a
-                b = (1 / self.dz ** 2) * (self.c2[i - 1][j])
-                self.M[count][count - 1] = b
-                c = -((1 / self.dtheta ** 2) * ((self.c1[i][j]) + self.c1[i][j - 1])
-                      + (1 / self.dz ** 2) * (self.c2[i][j] + self.c2[i - 1][j]))
-                self.M[count, count] = c
-                d = (1 / self.dz ** 2) * (self.c2[i][j])
-                self.M[count][count + 1] = d
-                e = (1 / self.dtheta ** 2) * (self.c1[i][j])
-                self.M[count][count + self.nz] = e
+                a = (1 / self.dtheta ** 2) * (c1[i, j - 1])
+                M[count, count - self.nz] = a
+                b = (1 / self.dz ** 2) * (c2[i - 1, j])
+                M[count, count - 1] = b
+                c = -((1 / self.dtheta ** 2) * ((c1[i, j]) + c1[i, j - 1])
+                      + (1 / self.dz ** 2) * (c2[i, j] + c2[i - 1, j]))
+                M[count, count] = c
+                d = (1 / self.dz ** 2) * (c2[i, j])
+                M[count, count + 1] = d
+                e = (1 / self.dtheta ** 2) * (c1[i, j])
+                M[count, count + self.nz] = e
                 count = count + 1
             count = count + 2
         count = 1
         for j in range(self.ntheta - 1):
             for i in range(1, self.nz - 1):
                 if j == 0:
-                    self.f[count][0] = (self.c0w[i][j] - self.c0w[i][self.ntheta - 1]) / self.dtheta
+                    f[count, 0] = (c0w[i, j] - c0w[i, self.ntheta - 1]) / self.dtheta
                 else:
-                    self.f[count][0] = (self.c0w[i, j] - self.c0w[i, j - 1]) / self.dtheta
+                    f[count, 0] = (c0w[i, j] - c0w[i, j - 1]) / self.dtheta
                 count = count + 1
             count = count + 2
         # fmt: on
+        return M, f
 
-    def resolves_matrix(self):
+    def resolves_matrix(self, M, f):
         """This function resolves the linear system [M]{P} = {f}.
+         Parameters
+        ----------
+        M: matrix of float
+            Matrix composed of coefficients that multiply the pressures at each point in the discrete domain.
+        f: array of float
+            Pressure independent terms.
+        Returns
+        --------
+        P: array of floats
+            Pressure. Unknowns of the finite difference system.
         Examples
         --------
         >>> my_fluid_flow = fluid_flow_example()
-        >>> my_fluid_flow.mounting_matrix()
-        >>> my_fluid_flow.resolves_matrix()
-        >>> my_fluid_flow.P # doctest: +ELLIPSIS
+        >>> c1, c2, c0w = my_fluid_flow.calculate_coefficients()
+        >>> M, f = my_fluid_flow.mounting_matrix(c1, c2, c0w)
+        >>> my_fluid_flow.resolves_matrix(M, f)# doctest: +ELLIPSIS
         array([[...
         """
-        sparse_matrix = sp.sparse.csc_matrix(self.M)
-        self.P = sp.sparse.linalg.spsolve(sparse_matrix, self.f)
-        self.P.shape = (self.P.size, 1)
+        sparse_matrix = sp.sparse.csc_matrix(M)
+        P = sp.sparse.linalg.spsolve(sparse_matrix, f)
+        P.shape = (P.size, 1)
+        return P
 
-    def calculate_pressure_matrix_numerical(self):
+    def calculate_pressure_matrix_numerical(self, direction=None):
         """This function calculates the pressure matrix numerically.
+        Parameters
+        ----------
+        direction: str
+            If defined, it calculates the model based on the disturbance in the chosen direction: 'x' or 'y'.
         Returns
         -------
         p_mat_numerical: matrix of float
@@ -541,15 +584,21 @@ class FluidFlow:
         >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
         array([[...
         """
-        self.mounting_matrix()
-        self.resolves_matrix()
+        if direction == "x":
+            c1, c2, c0w = self.calculate_coefficients(direction="x")
+        elif direction == "y":
+            c1, c2, c0w = self.calculate_coefficients(direction="y")
+        else:
+            c1, c2, c0w = self.calculate_coefficients()
+        M, f = self.mounting_matrix(c1, c2, c0w)
+        P = self.resolves_matrix(M, f)
         for i in range(self.nz):
             for j in range(self.ntheta):
                 k = j * self.nz + i
-                if self.P[k] < 0:
-                    self.p_mat_numerical[i][j] = 0
+                if P[k] < 0:
+                    self.p_mat_numerical[i, j] = 0
                 else:
-                    self.p_mat_numerical[i][j] = self.P[k]
+                    self.p_mat_numerical[i, j] = P[k]
         self.numerical_pressure_matrix_available = True
         return self.p_mat_numerical
 
@@ -575,7 +624,6 @@ def fluid_flow_example():
     my_pressure_matrix = FluidFlow(
         nz=8,
         ntheta=32,
-        nradius=8,
         length=0.04,
         omega=100.0 * 2 * np.pi / 60,
         p_in=0.0,
@@ -611,7 +659,6 @@ def fluid_flow_example2():
     """
     nz = 8
     ntheta = 16
-    nradius = 8
     length = 0.03
     omega = 157.1
     p_in = 0.0
@@ -624,7 +671,6 @@ def fluid_flow_example2():
     return FluidFlow(
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
diff --git a/ross/fluid_flow/fluid_flow_coefficients.py b/ross/fluid_flow/fluid_flow_coefficients.py
index e3c20087..9813904c 100644
--- a/ross/fluid_flow/fluid_flow_coefficients.py
+++ b/ross/fluid_flow/fluid_flow_coefficients.py
@@ -197,8 +197,8 @@ def calculate_stiffness_and_damping_coefficients(fluid_flow_object):
             * fluid_flow_object.xp
             * np.cos(fluid_flow_object.omegap * fluid_flow_object.t)
         )
-        fluid_flow_object.calculate_coefficients(direction="x")
-        fluid_flow_object.calculate_pressure_matrix_numerical()
+        fluid_flow_object.geometry_description()
+        fluid_flow_object.calculate_pressure_matrix_numerical(direction="x")
         [
             radial_force[i],
             tangential_force[i],
@@ -217,8 +217,8 @@ def calculate_stiffness_and_damping_coefficients(fluid_flow_object):
             * fluid_flow_object.yp
             * np.cos(fluid_flow_object.omegap * fluid_flow_object.t)
         )
-        fluid_flow_object.calculate_coefficients(direction="y")
-        fluid_flow_object.calculate_pressure_matrix_numerical()
+        fluid_flow_object.geometry_description()
+        fluid_flow_object.calculate_pressure_matrix_numerical(direction="y")
         [
             radial_force[i],
             tangential_force[i],
@@ -226,7 +226,7 @@ def calculate_stiffness_and_damping_coefficients(fluid_flow_object):
             force_yy[i],
         ] = calculate_oil_film_force(fluid_flow_object, force_type="numerical")
         move_rotor_center(fluid_flow_object, 0, -delta_y)
-        fluid_flow_object.calculate_coefficients()
+        fluid_flow_object.geometry_description()
         fluid_flow_object.calculate_pressure_matrix_numerical()
 
         X1[i] = [1, dx[i], xdot[i]]
diff --git a/ross/fluid_flow/fluid_flow_graphics.py b/ross/fluid_flow/fluid_flow_graphics.py
index 69d69be7..d951d9ab 100644
--- a/ross/fluid_flow/fluid_flow_graphics.py
+++ b/ross/fluid_flow/fluid_flow_graphics.py
@@ -394,11 +394,7 @@ def plot_pressure_theta_cylindrical(
         else:
             p_mat = fluid_flow_object.p_mat_numerical
 
-    r = np.linspace(
-        fluid_flow_object.radius_rotor,
-        fluid_flow_object.radius_stator,
-        fluid_flow_object.nradius,
-    )
+    r = np.linspace(fluid_flow_object.radius_rotor, fluid_flow_object.radius_stator)
     theta = np.linspace(
         0.0, 2.0 * np.pi + fluid_flow_object.dtheta / 2, fluid_flow_object.ntheta
     )
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 74adcbcf..61ca10a1 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -1043,7 +1043,7 @@ class Rotor(object):
         >>> rotor = rotor_example()
         >>> speed_range = rotor._clustering_points(num_modes=12, num_points=5)
         >>> speed_range.shape
-        (60,)
+        (61,)
         """
         critical_speeds = self.run_critical_speed(num_modes=num_modes, rtol=rtol)
         omega = critical_speeds.wd
@@ -1067,6 +1067,7 @@ class Rotor(object):
 
         omega = omega.reshape((len(omega), 1))
         speed_range = np.sort(np.ravel(np.concatenate((omega / a, omega * a))))
+        speed_range = np.insert(speed_range, 0, 0)
 
         return speed_range
 
@@ -1353,7 +1354,7 @@ class Rotor(object):
 
         >>> response = rotor.run_freq_response(cluster_points=True, num_points=5)
         >>> response.speed_range.shape
-        (60,)
+        (61,)
 
         # plot frequency response function:
         >>> fig = response.plot(inp=13, out=13)
@@ -1468,7 +1469,7 @@ class Rotor(object):
         ...     force=force, cluster_points=True, num_modes=12, num_points=5
         ... )
         >>> response.speed_range.shape
-        (60,)
+        (61,)
         """
         if speed_range is None:
             if cluster_points:
@@ -1633,7 +1634,7 @@ class Rotor(object):
         ...     node=3, magnitude=0.01, phase=0.0, cluster_points=True, num_points=5
         ... )
         >>> response2.speed_range.shape
-        (60,)
+        (61,)
 
         plot unbalance response:
         >>> fig = response.plot(dof=13)


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index 2515cd4d..2951e5fb 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -9,6 +9,7 @@ from numpy.testing import assert_allclose
 
 from ross.bearing_seal_element import (BallBearingElement, BearingElement,
                                        BearingElement6DoF,
+                                       MagneticBearingElement,
                                        RollerBearingElement)
 
 # fmt: on
@@ -292,6 +293,41 @@ def test_roller_bearing_element():
     assert_allclose(rollerbearing.G(), G)
 
 
+def test_magnetic_bearing_element():
+    n = 0
+    g0 = 1e-3
+    i0 = 1.0
+    ag = 1e-4
+    nw = 200
+    alpha = 0.392
+    kp_pid = 1.0
+    kd_pid = 1.0
+    k_amp = 1.0
+    k_sense = 1.0
+    tag = "magneticbearing"
+    mbearing = MagneticBearingElement(
+        n=n,
+        g0=g0,
+        i0=i0,
+        ag=ag,
+        nw=nw,
+        alpha=alpha,
+        kp_pid=kp_pid,
+        kd_pid=kd_pid,
+        k_amp=k_amp,
+        k_sense=k_sense,
+    )
+    M = np.array([[0.0, 0.0], [0.0, 0.0]])
+    K = np.array([[-4640.62337718, 0.0], [0.0, -4640.62337718]])
+    C = np.array([[4.64526865, 0.0], [0.0, 4.64526865]])
+    G = np.array([[0.0, 0.0], [0.0, 0.0]])
+
+    assert_allclose(mbearing.M(), M)
+    assert_allclose(mbearing.K(0), K)
+    assert_allclose(mbearing.C(0), C)
+    assert_allclose(mbearing.G(), G)
+
+
 @pytest.fixture
 def bearing_6dof():
     bearing_6dof = BearingElement6DoF(
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
index 5ea73416..2a419115 100644
--- a/ross/tests/test_fluid_flow.py
+++ b/ross/tests/test_fluid_flow.py
@@ -23,7 +23,6 @@ from ross.fluid_flow.fluid_flow_graphics import (
 def fluid_flow_short_eccentricity():
     nz = 8
     ntheta = 32
-    nradius = 8
     omega = 100.0 * 2 * np.pi / 60
     p_in = 0.0
     p_out = 0.0
@@ -36,7 +35,6 @@ def fluid_flow_short_eccentricity():
     return flow.FluidFlow(
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
@@ -53,7 +51,6 @@ def fluid_flow_short_eccentricity():
 def fluid_flow_short_friswell(set_load=True):
     nz = 8
     ntheta = 32
-    nradius = 8
     length = 0.03
     omega = 157.1
     p_in = 0.0
@@ -68,7 +65,6 @@ def fluid_flow_short_friswell(set_load=True):
         return flow.FluidFlow(
             nz,
             ntheta,
-            nradius,
             length,
             omega,
             p_in,
@@ -84,7 +80,6 @@ def fluid_flow_short_friswell(set_load=True):
         return flow.FluidFlow(
             nz,
             ntheta,
-            nradius,
             length,
             omega,
             p_in,
@@ -185,7 +180,6 @@ def test_damping_matrix():
 def fluid_flow_short_numerical():
     nz = 8
     ntheta = 32
-    nradius = 8
     length = 0.01
     omega = 100.0 * 2 * np.pi / 60
     p_in = 0.0
@@ -199,7 +193,6 @@ def fluid_flow_short_numerical():
     return flow.FluidFlow(
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
@@ -217,7 +210,6 @@ def fluid_flow_short_numerical():
 def fluid_flow_long_numerical():
     nz = 8
     ntheta = 32
-    nradius = 8
     omega = 100.0 * 2 * np.pi / 60
     p_in = 0.0
     p_out = 0.0
@@ -231,7 +223,6 @@ def fluid_flow_long_numerical():
     return flow.FluidFlow(
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,