From 3d6d3a7cfe68c8e1f4e8b2de49b6dd0157280052 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 3 Dec 2015 19:17:31 -0600 Subject: [PATCH 1/3] Try trig monomials as an AB function family --- examples/stability/plot-stability-regions.py | 42 +++++++++++---- leap/multistep/__init__.py | 55 +++++++++++++++++--- 2 files changed, 82 insertions(+), 15 deletions(-) diff --git a/examples/stability/plot-stability-regions.py b/examples/stability/plot-stability-regions.py index d69edc9..e216747 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 e543721..b69511c 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 * np.cos(self.alpha*n*x) + else: + return 1/n * np.sin(self.alpha*n*x) + + def emit_ab_integration(cb, name_gen, function_family, time_values, hist_vars, t_start, t_end): if isinstance(time_values, var): @@ -150,18 +179,32 @@ class AdamsBashforthMethod(Method): + component_id: The right hand side """ - 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 @@ -226,7 +269,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) @@ -303,7 +346,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 -- GitLab From acd3d664d3974e4b285f4312d65bb8d2e5b711fd Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 22 Aug 2017 13:18:42 -0500 Subject: [PATCH 2/3] Fix antiderivative in trig monomial AB --- leap/multistep/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index de2b072..0c3050a 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -98,9 +98,9 @@ class ABTrigMonomialIntegrationFunctionFamily(ABIntegrationFunctionFamily): if func_idx == 1: return x elif func_idx % 2 == 0: - return -1/n * np.cos(self.alpha*n*x) + return -1/(n*self.alpha) * np.cos(self.alpha*n*x) else: - return 1/n * np.sin(self.alpha*n*x) + return 1/(n*self.alpha) * np.sin(self.alpha*n*x) def _emit_func_family_operation(cb, name_gen, -- GitLab From ea9afcd1619719d53765185188f3d5b4aff3482b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 22 Aug 2017 13:19:09 -0500 Subject: [PATCH 3/3] Modify AB test to verify high order for trig AB --- test/test_ab.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/test/test_ab.py b/test/test_ab.py index 7da554e..694d5d2 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) -- GitLab