refactored md_vt_code into more functions

ipas/executables/collection_no_db/Ice_Agg_mD.py

@@ -125,13 +125,13 @@ def main():
  # savename (relative path)
  if save:
  filename = (
- "/network/rit/lab/sulialab/share/IPAS/ipas/instance_files/mD_vT_vars_flat"
+ "/network/rit/lab/sulialab/share/IPAS/ipas/instance_files/mD_vT_vars_flat_1"
  )
 
  # parallelize IPAS using dask
  use_dask = True
  if use_dask:
- num_workers = 15
+ num_workers = 50
  agg_as, agg_bs, agg_cs, Aps, Acs, Vps, Ves, Dmaxs = compute(
  phioarr,
  reqarr,

ipas/visualizations/mD_vT_relationships.py

@@ -13,31 +13,9 @@ class Relationships:
 
  ASPECT_RATIOS = [0.01, 0.1, 1.0, 10.0, 50.0]
  RHO_A = 1.395 # air density at -20C kg/m^3
-
  RHO_B = 916.8 # bulk density of ice [kg/m^3]
  GRAVITY = 9.81 # [m/sec^2]
 
- # Best # to Reynolds # power law fit coeffs
- ao = 1.0e-5
- bo = 1.0
- # sfc roughness
- delta_o = 5.83
- Co = 0.6
- C1 = 4 / (delta_o ** 2 * Co ** (1 / 2))
- C2 = delta_o ** 2 / 4
- # effective density coefficients
- k = (
- 0.07
- ) # Exponent in the terminal velocity versusdiameter relationship; aggs at ground H11 in chart
- n = 1.5 # Exponent in effective density relationship
- # Mitchell power laws for best number
- # all monomers and for limited size ranges
- # trying bullet rosette between D 200microns and 1000microns
- alpha = 0.00308
- beta = 2.26
- sigma = 1.57
- gamma = 0.0869
-
  def __init__(
  self, agg_as, agg_bs, agg_cs, phi_idxs, r_idxs, Aps, Acs, Vps, Ves, Dmaxs
  ):
@@ -116,42 +94,8 @@ def mass_spheroid_areas(self):
  m_spheroid[n] = 4 / 3 * np.pi * agg_as[n] ** 2 * agg_cs[n] * rho_i # kg
  return m_spheroid
 
- def b1(self, X):
- return (self.C1 * X ** (1 / 2)) / (
- 2
- * ((1 + self.C1 * X ** (1 / 2)) ** (1 / 2) - 1)
- * (1 + self.C1 * X ** (1 / 2)) ** (1 / 2)
- ) - (
- self.ao
- * self.bo
- * X ** self.bo
- / (self.C2 * ((1 + self.C1 * X ** (1 / 2)) ** (1 / 2) - 1) ** 2)
- )
-
- def a1(self, X):
- return self.C2 * (
- (1 + self.C1 * X ** (1 / 2)) ** (1 / 2) - 1
- ) ** 2 - self.ao * X ** self.bo / (X ** self.b1(X))
-
- def terminal_velocity_Mitchell_2005(self, Ar, X):
-
- m_ellipsoid_vol = self.mass_ellipsoid_volumes() # kg
- m_ellipsoid_vol = self.get_modes(m_ellipsoid_vol)
-
- D = self.Dmaxs[self.phi_idx, self.r_idx, :, :]
- D = self.get_modes(D)
-
- u = self.kinematic_viscosity()
+ def best_number_Mitchell(self, Ar, Ap, D, m):
 
- self.vt = (
- self.a1(X)
- * ((4 * self.GRAVITY * self.k) / (3 * self.RHO_A)) ** self.b1(X)
- * u ** (1 - 2 * self.b1(X))
- * D ** (3 * self.b1(X) - 1)
- * Ar ** ((self.n - 1) * self.b1(X))
- )
-
- def best_number(self, Ar, Ap, D, m):
  rho_p = self.RHO_B * (Ar)
  X = (
  (2 * m / rho_p)
@@ -164,27 +108,14 @@ def best_number(self, Ar, Ap, D, m):
  return X
 
  def best_number_Heymsfield(self, Ar, m):
+
  X = (self.RHO_A / self.eta ** 2) * (
  8 * m * self.GRAVITY / (np.pi * np.sqrt(Ar))
  )
  return X
 
- def best_number_Mitchell(self, D):
- """
- From Mitchell and Heymsfield 2005
- only using coeff from other studies,
- not taking into account IPAS values except Dmax
- """
- X = (
- 2
- * self.alpha
- * self.GRAVITY
- * self.RHO_A
- * D ** (self.beta + 2 - self.sigma)
- ) / (self.gamma * self.eta ** 2)
- return X
+ def reynolds_number_Mitchell(self, X):
 
- def reynolds_number(self, X):
  if X <= 10:
  a = 0.04394
  b = 0.970

ipas/visualizations/m_D_relationships.ipynb

@@ -3,7 +3,7 @@
  {
  "cell_type": "code",
  "execution_count": 1,
- "id": "93281d7e",
+ "id": "ac4fb2e6",
  "metadata": {},
  "outputs": [],
  "source": [
@@ -23,7 +23,7 @@
  {
  "cell_type": "code",
  "execution_count": 2,
- "id": "76f17386",
+ "id": "d4f8c5f1",
  "metadata": {},
  "outputs": [],
  "source": [
@@ -42,7 +42,7 @@
  {
  "cell_type": "code",
  "execution_count": 3,
- "id": "b8226db8",
+ "id": "478d6820",
  "metadata": {},
  "outputs": [
  {
@@ -114,7 +114,7 @@
  },
  {
  "cell_type": "markdown",
- "id": "ac67480e",
+ "id": "fe71632f",
  "metadata": {},
  "source": [
  "### Mass calculated using area vs. volume: <br>\n",
@@ -141,7 +141,7 @@
  {
  "cell_type": "code",
  "execution_count": 195,
- "id": "01921fe1",
+ "id": "f3a56677",
  "metadata": {},
  "outputs": [
  {
@@ -198,7 +198,7 @@
  {
  "cell_type": "code",
  "execution_count": 10,
- "id": "aed02808",
+ "id": "c10ef64d",
  "metadata": {},
  "outputs": [
  {
@@ -244,33 +244,35 @@
  "r_idxs = [0, 1, 2]\n",
  "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14,7), sharey=True, sharex=True)\n",
  "\n",
+ "# flag for plotting vt using Mitchells eqs.; otherwise, use Heymsfield\n",
+ "Mitchell = True\n",
  "# RANDOM ORIENTATION\n",
  "locals().update(result_rand)\n",
  "p = plot.Plots(ax1, df_CPI, agg_as, agg_bs, agg_cs, phi_idxs, r_idxs, Aps, Acs, Vps, Ves, Dmaxs)\n",
  "ylabel = \"$V_t$ [m/s]\"\n",
  "title='Random Orientation\\nCalculated Mass using Ar'\n",
- "p.vt_plot(title, ylabel, mflag = 'area', result_rand=True)\n",
+ "p.vt_plot(title, ylabel, mflag = 'area', Mitchell=Mitchell, result_rand=True)\n",
  "\n",
  "# QUASI-HORIZONTAL ORIENTATION\n",
  "locals().update(result_flat)\n",
  "p = plot.Plots(ax2, df_CPI, agg_as, agg_bs, agg_cs, phi_idxs, r_idxs, Aps, Acs, Vps, Ves, Dmaxs)\n",
  "title = f\"Quasi-Horizontal Orientation\\nCalculated Mass using Ar\"\n",
  "ylabel = \" \"\n",
- "p.vt_plot(title, ylabel, mflag = 'area', result_rand=False)\n",
+ "p.vt_plot(title, ylabel, mflag = 'area',  Mitchell=Mitchell, result_rand=False)\n",
  "\n",
  "# RANDOM ORIENTATION\n",
  "locals().update(result_rand)\n",
  "p = plot.Plots(ax3, df_CPI, agg_as, agg_bs, agg_cs, phi_idxs, r_idxs, Aps, Acs, Vps, Ves, Dmaxs)\n",
  "ylabel = \"$V_t$ [m/s]\"\n",
  "title='\\nCalculated Mass using Vr'\n",
- "p.vt_plot(title, ylabel, mflag = 'vol', result_rand=True)\n",
+ "p.vt_plot(title, ylabel, mflag = 'vol', Mitchell=Mitchell, result_rand=True)\n",
  "\n",
  "# QUASI-HORIZONTAL ORIENTATION\n",
  "locals().update(result_flat)\n",
  "p = plot.Plots(ax4, df_CPI, agg_as, agg_bs, agg_cs, phi_idxs, r_idxs, Aps, Acs, Vps, Ves, Dmaxs)\n",
  "ylabel = \" \"\n",
  "title='\\nCalculated Mass using Vr'\n",
- "p.vt_plot(title, ylabel, mflag = 'vol', result_rand=False)\n",
+ "p.vt_plot(title, ylabel, mflag = 'vol',  Mitchell=Mitchell, result_rand=False)\n",
  "#plt.savefig('../plots/vt.png', bbox_inches='tight')\n",
  "\n",
  "print(\"USING MITCHELL'S METHOD\")"
@@ -279,7 +281,7 @@
  {
  "cell_type": "code",
  "execution_count": 140,
- "id": "05a98c75",
+ "id": "157164da",
  "metadata": {},
  "outputs": [
  {
@@ -342,7 +344,7 @@
  {
  "cell_type": "code",
  "execution_count": 152,
- "id": "439655ce",
+ "id": "dab9a7de",
  "metadata": {},
  "outputs": [
  {
@@ -459,7 +461,7 @@
  {
  "cell_type": "code",
  "execution_count": 122,
- "id": "dbf0ca47",
+ "id": "0a922c69",
  "metadata": {},
  "outputs": [
  {
@@ -502,7 +504,7 @@
  {
  "cell_type": "code",
  "execution_count": null,
- "id": "8ee33d12",
+ "id": "b73add8e",
  "metadata": {},
  "outputs": [],
  "source": []