diff --git a/examples/stability/plot-stability-regions.py b/examples/stability/plot-stability-regions.py index d69edc9d698c20fb4e43b6e008c5c49d2ce1a744..e216747c3c2a043cf62278b3a4bd729488bc0a78 100644 --- a/examples/stability/plot-stability-regions.py +++ b/examples/stability/plot-stability-regions.py @@ -13,24 +13,46 @@ def plot_stability_region(code, parallel=None, scale_factor=None, **kwargs): fill(points.real, points.imag, **kwargs) -def main(save_pdfs=False): +def main(save_pdfs=True): import matplotlib.pyplot as pt - pt.rc("font", size=20) + pt.rc("font", size=15) #title("Stability Region") pt.xlabel(r"Re $\lambda$ / RHS calls") pt.ylabel(r"Im $\lambda$ / RHS calls") pt.grid() - import leap.rk as rk - import leap.multistep as multistep + import leap.rk as rk # noqa + import leap.multistep as multistep # noqa for label, method, factor in [ #("ode23", rk.ODE23Method("y", use_high_order=True), 1), #("ab2", multistep.AdamsBashforthMethod("y", 2), 1), - ("ab3", multistep.AdamsBashforthMethod("y", 3), 1), - #("ab4", multistep.AdamsBashforthMethod("y", 4), 1), - ("lserk", rk.LSRK4Method("y"), 1/5), - ("rk4", rk.RK4Method("y"), 1/4), + #("ab2", multistep.AdamsBashforthMethod("y", 2, static_dt=True), 1), + ("ab3", multistep.AdamsBashforthMethod("y", 3, static_dt=True), 1), + ("ab4-trig-0.7", multistep.AdamsBashforthMethod( + "y", multistep.ABTrigMonomialIntegrationFunctionFamily(4, 0.7), + static_dt=True), 1), + ("ab4-trig-0.9", multistep.AdamsBashforthMethod( + "y", multistep.ABTrigMonomialIntegrationFunctionFamily(4, 0.9), + static_dt=True), 1), + ("ab4-trig-1.1", multistep.AdamsBashforthMethod( + "y", multistep.ABTrigMonomialIntegrationFunctionFamily(4, 1.1), + static_dt=True), 1), + ("ab4-trig-1.3", multistep.AdamsBashforthMethod( + "y", multistep.ABTrigMonomialIntegrationFunctionFamily(4, 1.3), + static_dt=True), 1), + ("ab4-trig-1.5", multistep.AdamsBashforthMethod( + "y", multistep.ABTrigMonomialIntegrationFunctionFamily(4, 1.5), + static_dt=True), 1), + ("ab4-trig-2", multistep.AdamsBashforthMethod( + "y", multistep.ABTrigMonomialIntegrationFunctionFamily(4, 2), + static_dt=True), 1), + #("ab3-4", multistep.AdamsBashforthMethod("y", 3, hist_length=4), 1), + #("ab3-4", multistep.AdamsBashforthMethod("y", 3, hist_length=5), 1), + #("ab3-4", multistep.AdamsBashforthMethod("y", 3, hist_length=6), 1), + #("ab4", multistep.AdamsBashforthMethod("y", 4, static_dt=True), 1), + #("lserk", rk.LSRK4Method("y"), 1/5), + #("rk4", rk.RK4Method("y"), 1/4), ]: code = method.generate() @@ -39,12 +61,14 @@ def main(save_pdfs=False): pt.legend(labelspacing=0.1, borderpad=0.3, loc="best") if save_pdfs: + pt.tight_layout() pt.savefig("stab-regions.pdf") - xmin, xmax = pt.xlim() + xmin, xmax = pt.ylim() xsize = xmax-xmin pt.gca().set_aspect("equal") pt.ylim([-xsize/2*0.75, xsize/2*0.75]) + pt.tight_layout() pt.savefig("stab-regions-eq-aspect.pdf") else: diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index fc73fe04319659a6de5e1f0abbaafb4230904232..0c3050ae3d87ac3fae6475c97cdef1d279b39763 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -74,6 +74,35 @@ class ABMonomialIntegrationFunctionFamily(ABIntegrationFunctionFamily): return 1/(func_idx+1) * x**(func_idx+1) +class ABTrigMonomialIntegrationFunctionFamily(ABIntegrationFunctionFamily): + # FIXME: Doesn't yet work for on-the-fly coefficients + + def __init__(self, order, alpha): + self.order = order + self.alpha = alpha + + def __len__(self): + return self.order + + def evaluate(self, func_idx, x): + func_idx += 1 + n = func_idx // 2 + if func_idx % 2 == 0: + return np.sin(self.alpha*n*x) + else: + return np.cos(self.alpha*n*x) + + def antiderivative(self, func_idx, x): + func_idx += 1 + n = func_idx // 2 + if func_idx == 1: + return x + elif func_idx % 2 == 0: + return -1/(n*self.alpha) * np.cos(self.alpha*n*x) + else: + return 1/(n*self.alpha) * np.sin(self.alpha*n*x) + + def _emit_func_family_operation(cb, name_gen, function_family, time_values, hist_vars, rhs_func): if isinstance(time_values, var): @@ -164,18 +193,32 @@ class AdamsBashforthMethod(Method): .. automethod:: generate """ - def __init__(self, component_id, order, state_filter_name=None, - hist_length=None, static_dt=False): + def __init__(self, component_id, function_family=None, state_filter_name=None, + hist_length=None, static_dt=False, order=None): """ + :arg function_family: Accepts an instance of + :class:`ABIntegrationFunctionFamily` + or an integer, in which case the classical monomial function family + with the order given by the integer is used. :arg static_dt: If *True*, changing the timestep during time integration is not allowed. """ + if function_family is not None and order is not None: + raise ValueError("may not specify both function_family and order") + + if function_family is None: + function_family = order + del order + + if isinstance(function_family, int): + function_family = ABMonomialIntegrationFunctionFamily(function_family) + super(AdamsBashforthMethod, self).__init__() - self.order = order + self.function_family = function_family if hist_length is None: - hist_length = order + hist_length = len(function_family) self.hist_length = hist_length self.static_dt = static_dt @@ -243,7 +286,7 @@ class AdamsBashforthMethod(Method): ab_sum = emit_ab_integration( cb_primary, name_gen, - ABMonomialIntegrationFunctionFamily(self.order), + self.function_family, time_hist, history, 0, t_end) @@ -319,7 +362,7 @@ class AdamsBashforthMethod(Method): cb(self.time_history[i], self.t) from leap.rk import ORDER_TO_RK_METHOD - rk_method = ORDER_TO_RK_METHOD[self.order] + rk_method = ORDER_TO_RK_METHOD[self.function_family.order] rk_tableau = tuple(zip(rk_method.c, rk_method.a_explicit)) rk_coeffs = rk_method.output_coeffs diff --git a/test/test_ab.py b/test/test_ab.py index 7da554e72847b65b37b561a78f845140fc887a23..694d5d223005d25f967334f9e513d539ffb0773f 100755 --- a/test/test_ab.py +++ b/test/test_ab.py @@ -27,21 +27,28 @@ THE SOFTWARE. import sys import pytest from leap.multistep import AdamsBashforthMethod +import leap.multistep as multistep from utils import ( # noqa python_method_impl_interpreter as pmi_int, python_method_impl_codegen as pmi_cg) -@pytest.mark.parametrize("order", [1, 2, 3, 4, 5]) -@pytest.mark.parametrize("static_dt", [True, False]) -def test_ab_accuracy(python_method_impl, order, static_dt, show_dag=False, - plot_solution=False): +@pytest.mark.parametrize(("method", "expected_order"), [ + (AdamsBashforthMethod("y", order, static_dt=static_dt), order) + for order in [1, 3, 5] + for static_dt in [True, False] + ] + [ + (AdamsBashforthMethod("y", + multistep.ABTrigMonomialIntegrationFunctionFamily(4, 0.1), + static_dt=True), 3) + ]) +def test_ab_accuracy(python_method_impl, method, expected_order, + show_dag=False, plot_solution=False): from utils import check_simple_convergence - method = AdamsBashforthMethod("y", order, static_dt=static_dt) #method = AdamsBashforthMethod("y", order, hist_length=order+1) check_simple_convergence(method=method, method_impl=python_method_impl, - expected_order=order, show_dag=show_dag, + expected_order=expected_order, show_dag=show_dag, plot_solution=plot_solution)