v0.0.5 (#8)

Author: ngam

hypersolver/__init__.py

@@ -12,36 +12,77 @@
     f is speed n moves along x, and
     g lumps sources and sinks
 
-    functionally, n(x; t), f(x), and g(x; t; n)
+    functionally, n(x; t), f(x), and g(n; x)
 
     note, fn is the flux across x.
 
     Usage:
-    >>> from hypersolver import select_solver
-    >>> solver = select_solver(method)
+    >>> from hypersolver import solver
     >>> solver(n0, x, t, f, g, **kwargs)
+    >>> # kwargs include "method", "backend", etc.
 
-    available methods:
+    available `method`s:
+    pde:
         - "lax_friedrichs" (default)
-        - "lax_wendroff" (still unstable, wip)
-        - "method_of_characteristics" (experimental)
+        - "lax_wendroff"
+        - "method_of_characteristics" (broken, experimental)
+    ode:
+        - "rk2"
 
+    available `backend`s:
+        - "numpy" (default)
+        - "jax" (experimental)
+
+    available `solver_type`s:
+        - "unsplit" (default)
+        - "split"
 """
 
+import os
+
+from hypersolver.pde_solver_unsplit import solver_ as solver_upde
+# from hypersolver.pde_solver_split import solver_ as solver_spde
+from hypersolver.ode_solver import solver_ as solver_ode
 
-__version__ = "0.0.4"
 
-__hyper_solvers__ = [
+__version__ = "0.0.5"
+
+__hyper_methods__ = [
     "lax_friedrichs",
     "lax_wendroff",
     "method_of_characteristics",
+    "rk2",
 ]
 
-from hypersolver.step_solver import solver_
+__hyper_solver_types__ = [
+    "unsplit",
+    "split",
+]
 
 
-def solver(*args, method="lax_friedrichs", **kwargs):
+def solver(
+    *args,
+    method="lax_friedrichs",
+    backend=os.environ.get("HS_BACKEND", "numpy"),
+    verbosity=os.environ.get("HS_VERBOSITY", "0"),
+    solver_type="unsplit",
+    **kwargs
+):
     """ wrapper function to select solvers """
-    if method not in __hyper_solvers__:
+
+    os.environ["HS_BACKEND"] = str(backend)
+    os.environ["HS_VERBOSITY"] = str(verbosity)
+
+    if method not in __hyper_methods__ or \
+            solver_type not in __hyper_solver_types__:
         raise ValueError("method not supported")
-    return solver_(*args, method=method, **kwargs)
+
+    if method.startswith("rk"):
+        return solver_ode(*args, method=method, **kwargs)
+    # if method.endswith("_split"):
+    #     return solver_spde(*args, method=method, **kwargs)
+    return solver_upde(
+        *args,
+        method=method,
+        **kwargs
+    )

hypersolver/accurate_derivative.py

@@ -0,0 +1,46 @@
+""" derivative: calculating derivatives central differencing
+
+    available functions:
+        - ord1_acc2: order=1, accuracy=2
+        - ord2_acc2: order=2, accuracy=2
+"""
+
+from hypersolver.util import xnp as np
+
+
+def ord1_acc2(_func, _xvar):
+    """ central differencing: order=1, accuracy=2 """
+
+    _func, _xvar = np.array(_func), np.array(_xvar)
+
+    _derivative = (
+        _func[2:] - _func[:-2]
+    )/(_xvar[2:] - _xvar[:-2])
+
+    return np.pad(
+        _derivative,
+        (1, 1),
+        mode='constant',
+        constant_values=(
+            (_func[1] - _func[0])/(_xvar[1] - _xvar[0]),
+            (_func[-1] - _func[-2])/(_xvar[-1] - _xvar[-2])
+        ),)
+
+
+def ord2_acc2(_func, _xvar):
+    """ central differencing: order=2, accuracy=2 """
+
+    _func, _xvar = np.array(_func), np.array(_xvar)
+
+    _derivative = (
+        _func[2:] - 2.0*_func[1:-1] + _func[:-2]
+    )/((_xvar[2:] - _xvar[:-2])/2.0)**2
+
+    return np.pad(
+        _derivative,
+        (1, 1),
+        mode='constant',
+        constant_values=(
+            (_func[2] - 2.0*_func[1] + _func[0])/(_xvar[1] - _xvar[0])**2,
+            (_func[-1] - 2.0*_func[-2] + _func[-3])/(_xvar[-1] - _xvar[-2])**2
+        ),)

hypersolver/derivative.py

@@ -1,6 +1,24 @@
 """ Lax-Friedrics finite-difference scheme """
 
-from hypersolver.util import prep_next_step
+from hypersolver.util import xnp as np
+from hypersolver.util import time_step_util
+from hypersolver.derivative import ord1_acc2
+
+
+def lx_init(init_vals):
+    """ initialize the array
+
+        pad the array with the prescribed first and last values
+    """
+
+    return np.pad(
+        init_vals,
+        (1, 1),
+        mode="constant",
+        constant_values=(
+            0.5 * (init_vals[1] + init_vals[0]),
+            0.5 * (init_vals[-1] + init_vals[-2])
+        ),)
 
 
 def lx_next(
@@ -16,7 +34,7 @@ def lx_next(
 
         inputs
         ------
-        init_step:  n
+        init_vals:  n
         vars_vals:  x
         flux_term:  f
         sink_term:  g
@@ -28,46 +46,16 @@ def lx_next(
 
         numerics
         --------
-        n(j+1, i) = (
-            0.5 * (n(j,i+1) + n(j+1,i-1)) -
-            time_step / (x(i+1) - x(i-1)) * (
-                n(j,i+1) * f(i+1) -
-                n(j,i-1) * f(i-1)
-            ) +
-            g(j,i) * time_step
+        n(j+1, i) = n(j, i) + Δt (g - Δ(fn)/Δx)(j, i)
 
-        time_step = (
-            stability *
-            (x(i+1) - x(i-1)).min() /
-            (f(i)).max()
-        )
+        Δt ≤ λΔx/f ∀ x
+        Δ(fn)/Δx is first-order derivative with accuracy of 2
+        n(j, i) = (n(j, i-1) + n(j, i+1))/2
     """
 
-    (time_step, next_vals) = prep_next_step(
-        stability, vars_vals, flux_term, init_vals)
-
-    next_vals[1:-1] = (
-        0.5 * (init_vals[2:] + init_vals[:-2]) -
-        1.0 * time_step / (vars_vals[2:] - vars_vals[:-2]) * (
-            init_vals[2:] * flux_term[2:] -
-            init_vals[:-2] * flux_term[:-2]
-        ) + sink_term[1:-1] * time_step
-    )
-
-    next_vals[0] = (
-        0.5 * (init_vals[1] + init_vals[0]) -
-        1.0 * time_step / (vars_vals[1] - vars_vals[0]) * (
-            init_vals[1] * flux_term[1] -
-            init_vals[0] * flux_term[0]
-        ) + sink_term[0] * time_step
-    )
+    time_step = time_step_util(vars_vals, flux_term, stability)
 
-    next_vals[-1] = (
-        0.5 * (init_vals[-1] + init_vals[-2]) -
-        1.0 * time_step / (vars_vals[-1] - vars_vals[-2]) * (
-            init_vals[-1] * flux_term[-1] -
-            init_vals[-2] * flux_term[-2]
-        ) + sink_term[-1] * time_step
-    )
+    _init_vals = lx_init(0.5 * (init_vals[2:] + init_vals[:-2]))
 
-    return next_vals
+    return _init_vals - time_step * (
+        ord1_acc2(init_vals*flux_term, vars_vals) - sink_term)