v0.0.4 (#7)

Author: ngam

hypersolver/__init__.py

@@ -3,7 +3,7 @@
     hypersolver revolves around solving hyperbolic
     partial differential equations (PDEs) of the form
 
-    ∂n/∂t + ∂(fn)/∂x = ∂n/∂t + f ∂n/∂x - n ∂f/∂x = g
+    ∂n/∂t + ∂(fn)/∂x = ∂n/∂t + f ∂n/∂x + n ∂f/∂x = g
 
     where
 
@@ -23,7 +23,8 @@
 
     available methods:
         - "lax_friedrichs" (default)
-        - "lax_wendroff"
+        - "lax_wendroff" (still unstable, wip)
+        - "method_of_characteristics" (experimental)
 
 """
 
@@ -36,11 +37,11 @@
     "method_of_characteristics",
 ]
 
-from hypersolver.basic_solver import solver
+from hypersolver.step_solver import solver_
 
 
-def select_solver(method="lax_friedrichs"):
+def solver(*args, method="lax_friedrichs", **kwargs):
     """ wrapper function to select solvers """
     if method not in __hyper_solvers__:
         raise ValueError("method not supported")
-    return solver(method)
+    return solver_(*args, method=method, **kwargs)

hypersolver/accurate_derivative.py

@@ -3,20 +3,26 @@
 import numpy as np
 from scipy.misc import central_diff_weights
 
+weights = np.zeros((4, 4))
+weights[0, -1:] = np.array([-1/2])
+weights[1, -2:] = np.array([1/12, -2/3])
+weights[2, -3:] = np.array([-1/60, 3/20, -3/4])
+weights[3, -4:] = np.array([1/280, -4/105, 1/5, -4/5])
+
 
 def _derivative(_func, _xvar, _nacc):
     """ 1st derivative following central finite differencing
     """
     derivative = np.zeros_like(_xvar)
     derivative[0] = (_func[1] - _func[0])/(_xvar[1] - _xvar[0])
     derivative[-1] = (_func[-1] - _func[-2])/(_xvar[-1] - _xvar[-2])
-    weights = central_diff_weights(_nacc+1)
+    _weights = central_diff_weights(_nacc+1)
 
     for nac, idx in zip(range(2, _nacc + 1, 2), range(_nacc//2-1, -1, -1)):
         derivative[nac//2:-nac//2] += (
-            weights[idx]*_func[:-nac] /
+            _weights[idx]*_func[:-nac] /
             ((_xvar[nac:] - _xvar[:-nac]) / float(nac)) +
-            weights[-idx-1]*_func[nac:] /
+            _weights[-idx-1]*_func[nac:] /
             ((_xvar[nac:] - _xvar[:-nac]) / float(nac)))
 
     if _nacc > 2:
@@ -26,6 +32,30 @@ def _derivative(_func, _xvar, _nacc):
     return derivative
 
 
+def acc4_derivative(_func, _xvar, _nacc=4):
+    """ 1st derivative, central differencing, accuracy up to 4
+    """
+    derivative = np.zeros_like(_xvar)
+    derivative[0] = (_func[1] - _func[0])/(_xvar[1] - _xvar[0])
+    derivative[-1] = (_func[-1] - _func[-2])/(_xvar[-1] - _xvar[-2])
+
+    for nac, idx in zip(range(2, _nacc + 1, 2), range(_nacc//2)):
+
+        derivative[nac//2] += (
+            weights[nac//2-1, -idx-1]*(_func[0] - _func[nac]) /
+            ((_xvar[nac] - _xvar[0]) / float(nac)))
+
+        derivative[-nac//2-1] += (
+            weights[nac//2-1, -idx-1]*(_func[-nac-1] - _func[-1]) /
+            ((_xvar[-1] - _xvar[-nac-1]) / float(nac)))
+
+        derivative[nac//2+1:-nac//2-1] += (
+            weights[_nacc//2-1, -idx-1]*(_func[1:-nac-1] - _func[nac+1:-1]) /
+            ((_xvar[nac+1:-1] - _xvar[1:-nac-1]) / float(nac)))
+
+    return derivative
+
+
 def acc_derivative(func, xvar, nacc):
     """ 1st derivative of a _func wrt _xvar with _nacc accuracy
 
@@ -40,15 +70,15 @@ def acc_derivative(func, xvar, nacc):
         NOTES:
             minor bug necessitates repeating calculation if nacc > 2
     """
-    if nacc == 0 or nacc % 2 == 1:
+    if nacc <= 0 or nacc % 2 == 1:
         raise ValueError("n must be positive even")
 
+    if nacc < 6:
+        return acc4_derivative(func, xvar, nacc)
+
     axx_derivative = np.zeros_like(xvar)
     axx_derivative = _derivative(func, xvar, 2)
 
-    if nacc == 2:
-        return axx_derivative
-
     for nax in range(4, nacc + 1, 2):
         derivative = _derivative(func, xvar, nax)
         nan_idx = np.isnan(derivative)

hypersolver/lax_wendroff.py

@@ -28,6 +28,16 @@ def lw_next(
 
         numerics
         --------
+        n(j+1, i) = (
+            n(j,i) +
+            time_step * ( -(fn)x + g ) + # first term (order)
+            0.5 * (time_step)**2 * (
+                -fx (-(fn)x + g) - f(-(fn)xx + gx) + gt) # second term (order)
+            )
+        )
+
+        # FIXME: wrong implementation for now
+        # TODO: to be revisited
         n(j+1, i) = (
             n(j,i) -
             time_step / (x(i+1) - x(i-1)) * (
@@ -53,9 +63,9 @@ def lw_next(
         1.0 * time_step / (vars_vals[2:] - vars_vals[:-2]) * (
             init_vals[2:] * flux_term[2:] -
             init_vals[:-2] * flux_term[:-2]
-        ) + (time_step / (vars_vals[2:] - vars_vals[:-2])/2)**2 * (
+        ) + 0.5 * (time_step / ((vars_vals[2:] - vars_vals[:-2])/2.0))**2.0 * (
             init_vals[:-2] * flux_term[:-2] -
-            2 * init_vals[1:-1] * flux_term[1:-1] +
+            2.0 * init_vals[1:-1] * flux_term[1:-1] +
             init_vals[2:] * flux_term[2:]
         ) + sink_term[1:-1] * time_step
     )
@@ -65,7 +75,7 @@ def lw_next(
         1.0 * time_step / (vars_vals[1] - vars_vals[0]) * (
             init_vals[1] * flux_term[1] -
             init_vals[0] * flux_term[0]
-        ) + (time_step / (vars_vals[1] - vars_vals[0])/1)**2 * (
+        ) + 0.5 * ((time_step / (vars_vals[1] - vars_vals[0])/1.0))**2.0 * (
             0.0 * init_vals[0] * flux_term[0] -
             2.0 * init_vals[0] * flux_term[0] +
             1.0 * init_vals[1] * flux_term[1]
@@ -77,7 +87,7 @@ def lw_next(
         1.0 * time_step / (vars_vals[-1] - vars_vals[-2]) * (
             init_vals[-1] * flux_term[-1] -
             init_vals[-2] * flux_term[-2]
-        ) + (time_step / (vars_vals[-1] - vars_vals[-2])/1)**2 * (
+        ) + 0.5 * ((time_step / (vars_vals[-1] - vars_vals[-2])/1.0))**2.0 * (
             1.0 * init_vals[-2] * flux_term[-2] -
             2.0 * init_vals[-1] * flux_term[-1] +
             0.0 * init_vals[-1] * flux_term[-1]

hypersolver/method_of_characteristics.py

@@ -2,17 +2,17 @@
 """
 import numpy as np
 from scipy.integrate import odeint
+from scipy.interpolate import interp1d
 
 from hypersolver.accurate_derivative import acc_derivative
 
 
-def moc(
+def moc_next(
     init_vals,
     vars_vals,
-    time_span,
     flux_term,
     sink_term,
-    **kwargs
+    time_step,
 ):
     """ method of characteristics
 
@@ -29,20 +29,18 @@ def moc(
         ------
         init_step:  n
         vars_vals:  x
-        time_span:
         flux_term:  f
         sink_term:  g
+        time_step:
 
         outputs
         -------
-        next_vars:  x
         next_vals:  n
 
         numerics
         --------
-        use from scipy.integrate.odeint
+        use scipy.integrate.odeint
     """
-    _ = kwargs
 
     # pylint: disable=unused-argument
     # pylint: disable=unused-variable
@@ -66,7 +64,14 @@ def _func(yval, tval):
     yval0 = np.empty((vars_vals.size + init_vals.size))
     yval0[::2] = vars_vals
     yval0[1::2] = init_vals
-    tspan = np.linspace(time_span[0], time_span[-1], 100)
+    tspan = np.linspace(0, time_step, 10)
     results = odeint(_func, yval0, tspan, ml=2, mu=2)
 
-    return (results[:, ::2], results[:, 1::2])
+    fill = interp1d(
+        results[-1, ::2],
+        results[-1, 1::2],
+        fill_value=(0.0, 0.0),
+        bounds_error=False,
+        kind='cubic')
+
+    return fill(vars_vals)

hypersolver/basic_solver.py hypersolver/step_solver.pyhypersolver/basic_solver.py renamed to hypersolver/step_solver.py

@@ -1,23 +1,25 @@
 """ shared solver between schemes
 """
+
 import numpy as np
 
 from hypersolver.util import term_util
 from hypersolver.lax_friedrichs import lx_next
 from hypersolver.lax_wendroff import lw_next
-from hypersolver.method_of_characteristics import moc
+from hypersolver.method_of_characteristics import moc_next
 
 
-def solver(method):  # noqa: C901
+def solver_(*args, **kwargs):
     """ set the solver """
+    method = kwargs.get("method", "lax_friedrichs")
     if method == "lax_friedrichs":
         next_step = lx_next
-    if method == "lax_wendroff":
+    elif method == "lax_wendroff":
         next_step = lw_next
-    if method == "method_of_characteristics":
-        return moc
+    else:
+        next_step = moc_next
 
-    def _solver(  # pylint: disable=too-many-arguments
+    def _solver_(
         init_vals,
         vars_vals,
         time_span,
@@ -39,11 +41,15 @@ def _solver(  # pylint: disable=too-many-arguments
         flux_term = term_util(flux_term, init_vals)
         sink_term = term_util(sink_term, init_vals)
 
-        time_step = (
-            stability_factor *
-            np.diff(vars_vals).min() /
-            np.abs(flux_term).max()
-        )
+        stability_factor, time_step = (stability_factor,
+                                       stability_factor *
+                                       np.diff(vars_vals).min() /
+                                       np.abs(flux_term).max()
+                                       ) if method in [
+            "lax_friedrichs", "lax_wendroff"
+        ] else (
+            np.array((time_span[-1] - time_span[0])/5.0),
+            np.array((time_span[-1] - time_span[0])/5.0))
 
         tidx = np.arange(time_span[0], time_span[-1]+time_step, time_step)
         sols = np.zeros((tidx.size, init_vals.size))
@@ -74,4 +80,4 @@ def _solver(  # pylint: disable=too-many-arguments
 
         return sols
 
-    return _solver
+    return _solver_(*args, **kwargs)