Merge branch 'optimizing'

Author: sigvaldm

langmuir/finite_length.py

@@ -24,7 +24,7 @@
 from langmuir.geometry import *
 from langmuir.species import *
 from langmuir.misc import *
-from scipy.interpolate import griddata
+from scipy.interpolate import griddata, RectBivariateSpline
 from scipy.constants import value as constants
 from copy import deepcopy
 import numpy as np
@@ -101,7 +101,7 @@ def finite_length_current_density(geometry, species, V=None, eta=None,
         eta_isarray = isinstance(eta, (np.ndarray, list, tuple))
         eta = make_array(eta)
 
-    eta = eta[:, None] # Make eta rows
+    # eta = eta[:, None] # Make eta rows
 
     if not isinstance(geometry, Cylinder):
         raise ValueError('Geometry not supported: {}'.format(geometry))
@@ -124,6 +124,13 @@ def finite_length_current_density(geometry, species, V=None, eta=None,
 
     C, A, alpha, delta = get_lerped_coeffs(lambd_t, eta)
 
+    # Make these column vectors, i.e. one eta per row
+    C = C[:, None]
+    A = A[:, None]
+    alpha = alpha[:, None]
+    delta = delta[:, None]
+    eta = eta[:, None]
+
     if normalization is None: # i0 = i_OML => i = i_OML * g
         geonorm = deepcopy(geometry)
         geonorm.l = 1
@@ -253,47 +260,97 @@ def H(zeta):
 
     return I if eta_isarray else I[0]
 
-def get_lerped_coeffs(lambd, eta):
-    """
-    Fetches and interpolates Marholm-Marchand coefficients.
+# def get_lerped_coeffs(lambd, eta):
+#     """
+#     Fetches and interpolates Marholm-Marchand coefficients.
 
-    Parameters
-    ----------
-    lambd: float
-        Normalized probe length (lambda) to get coefficients for
+#     Parameters
+#     ----------
+#     lambd: float
+#         Normalized probe length (lambda) to get coefficients for
 
-    eta: float
-        Normalized probe voltage (eta) to get coefficients for
+#     eta: float
+#         Normalized probe voltage (eta) to get coefficients for
 
-    Returns
-    -------
-    4-tuple of coefficients (C, A, alpha, delta)
-    """
+#     Returns
+#     -------
+#     4-tuple of coefficients (C, A, alpha, delta)
+#     """
 
-    fname = os.path.join(os.path.dirname(os.path.abspath(__file__)),'params.npz')
-    file = np.load(fname)
+#     fname = os.path.join(os.path.dirname(os.path.abspath(__file__)),'params.npz')
+#     file = np.load(fname)
+
+#     lambds = file['lambds']
+#     etas = file['etas']
+#     Cs = file['Cs']
+#     As = file['As']
+#     alphas = file['alphas']
+#     deltas = file['deltas']
+
+#     # Extrapolate for larger lambda by using the largest available
+#     lambd_coeff = min(lambd, max(lambds))
+
+#     # Tabulated values contains datapoints as described in paper
+#     C     = griddata((lambds, etas), Cs    , (lambd_coeff, eta))
+#     A     = griddata((lambds, etas), As    , (lambd_coeff, eta))
+#     alpha = griddata((lambds, etas), alphas, (lambd_coeff, eta))
+#     delta = griddata((lambds, etas), deltas, (lambd_coeff, eta))
+
+#     # Make repelled species identical to OML through these coefficients
+#     ind = np.where(eta>=0)[0]
+#     C[ind] = 1
+#     A[ind] = 0
+#     alpha[ind] = 1
+#     delta[ind] = 1
+
+#     return C, A, alpha, delta
 
+class lerper(RectBivariateSpline):
+    def __init__(self, coeffname, degree=1):
+
+        assert coeffname in ['C', 'A', 'alpha', 'delta']
+        coeffname += 's'
+
+        fname = os.path.join(os.path.dirname(os.path.abspath(__file__)),'params_structured.npz')
+        file = np.load(fname)
+
+        lambds = file['lambds']
+        etas = file['etas']
+        coeff = file[coeffname]
+
+        super(lerper, self).__init__(lambds, etas, coeff, kx=degree, ky=degree)
+
+def get_max_lambd():
+    fname = os.path.join(os.path.dirname(os.path.abspath(__file__)),'params_structured.npz')
+    file = np.load(fname)
     lambds = file['lambds']
-    etas = file['etas']
-    Cs = file['Cs']
-    As = file['As']
-    alphas = file['alphas']
-    deltas = file['deltas']
+    return max(lambds)
 
-    # Extrapolate for larger lambda by using the largest available
-    lambd_coeff = min(lambd, max(lambds))
+lerp_C = lerper('C')
+lerp_A = lerper('A')
+lerp_alpha = lerper('alpha')
+lerp_delta = lerper('delta')
+max_lambd = get_max_lambd()
 
-    # Tabulated values contains datapoints as described in paper
-    C     = griddata((lambds, etas), Cs    , (lambd_coeff, eta))
-    A     = griddata((lambds, etas), As    , (lambd_coeff, eta))
-    alpha = griddata((lambds, etas), alphas, (lambd_coeff, eta))
-    delta = griddata((lambds, etas), deltas, (lambd_coeff, eta))
+def get_lerped_coeffs_new(lambd, eta):
+    eta = -eta
+
+    # Extrapolate for larger lambda by using the largest available
+    lambd_coeff = min(lambd, max_lambd)
 
     # Make repelled species identical to OML through these coefficients
-    ind = np.where(eta>=0)[0]
-    C[ind] = 1
-    A[ind] = 0
-    alpha[ind] = 1
-    delta[ind] = 1
+    C = np.ones_like(eta)
+    A = np.zeros_like(eta)
+    alpha = np.ones_like(eta)
+    delta = np.ones_like(eta)
+
+    # Tabulated values contains datapoints as described in paper
+    ind = np.where(eta>0)[0]
+    C[ind] = lerp_C(lambd_coeff, eta[ind], grid=False)
+    A[ind] = lerp_A(lambd_coeff, eta[ind], grid=False)
+    alpha[ind] = lerp_alpha(lambd_coeff, eta[ind], grid=False)
+    delta[ind] = lerp_delta(lambd_coeff, eta[ind], grid=False)
 
     return C, A, alpha, delta
+
+get_lerped_coeffs=get_lerped_coeffs_new

langmuir/params_structured.npz

@@ -20,6 +20,7 @@
 """
 
 from scipy.constants import value as constants
+from langmuir.misc import *
 import numpy as np
 
 class Species(object):
@@ -109,6 +110,11 @@ def __init__(self, *args, **kwargs):
         else:
             self.T = kwargs['T']
 
+        if self.T<0:
+            self.T=0
+            logger.warning('Negative temperature interpreted as zero')
+
+
         self.vth = np.sqrt(k*self.T/self.m)
         self.eV  = self.T*k/e
         self.omega_p = np.sqrt(self.q**2*self.n/(eps0*self.m))

langmuir/species.py

@@ -197,9 +197,16 @@ def test_current_interpolated_samples():
     r = 1e-3
     ls = np.array([7, 70, 700], dtype=np.float)*1e-3
     etas = np.array([65, 1, 90], dtype=np.float)
-    Ifs = np.array([-6.2794005496528275e-06,
-                    -1.9760221828553347e-06,
-                    -0.00014456804182897328])
+
+    # Ground truth when using scipy.interpolate.griddata
+    # Ifs = np.array([-6.2794005496528275e-06,
+    #                 -1.9760221828553347e-06,
+    #                 -0.00014456804182897328])
+
+    # Ground truth when using 1st order scipy.interpolate.RectBivariateSpline
+    Ifs = np.array([-6.439906296987808e-06,
+                    -1.9562053783242e-06,
+                    -0.00014420124702985558])
 
     for l, eta, If in zip(ls, etas, Ifs):
         geo = Cylinder(r=r, l=l)
@@ -217,9 +224,16 @@ def test_current_density_interpolated_samples():
     ls = np.array([7, 70, 700], dtype=np.float)*1e-3
     etas = np.array([65, 1, 90], dtype=np.float)
     zs = np.array([0, .4, .29, 1])
-    Ifs = np.array([[-0.00045615, -0.00106614, -0.00099727, -0.00045615],
-                    [-2.92817683e-05, -2.75155255e-05, -2.76705163e-05, -2.92817683e-05],
-                    [-0.00036732, -0.00019135, -0.00019135, -0.00036732]])
+
+    # Ground truth when using scipy.interpolate.griddata
+    # Ifs = np.array([[-0.00045615, -0.00106614, -0.00099727, -0.00045615],
+    #                 [-2.92817683e-05, -2.75155255e-05, -2.76705163e-05, -2.92817683e-05],
+    #                 [-0.00036732, -0.00019135, -0.00019135, -0.00036732]])
+
+    # Ground truth when using 1st order scipy.interpolate.RectBivariateSpline
+    Ifs = np.array([[-0.00054336, -0.00105909, -0.00100473, -0.00054336],
+                    [-2.89516395e-05, -2.72642938e-05, -2.74123613e-05, -2.89516395e-05],
+                    [-0.00036957, -0.00019082, -0.00019082, -0.00036957]])
 
     for l, eta, If in zip(ls, etas, Ifs):
         geo = Cylinder(r=r, l=l)

langmuir/test_finite_length.py

@@ -151,3 +151,7 @@ def test_debye():
     plasma.append(Species('proton',   n=1e11, T=1000))
     assert(debye(plasma[0])==plasma[0].debye)
     assert(debye(plasma) == approx(0.004879671013271479))
+
+def test_negative_temperature(caplog):
+    sp = Species(n=1e11, T=-1)
+    assert(caplog.records[0].levelname == 'WARNING')