diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index fc73fe04319659a6de5e1f0abbaafb4230904232..0787df3b6a76b13e6bbffd05d58f58daa6a181ec 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -84,6 +84,9 @@ def _emit_func_family_operation(cb, name_gen, array = var("array") linear_solve = var("linear_solve") + svd = var("svd") + matmul = var("matmul") + transpose = var("transpose") # use: # Vandermonde^T * ab_coeffs = integrate(t_start, t_end, monomials) @@ -92,7 +95,7 @@ def _emit_func_family_operation(cb, name_gen, cb(vdmt, array(nfunctions*hist_len)) coeff_rhs = var(name_gen("coeff_rhs")) - cb(coeff_rhs, array(hist_len)) + cb(coeff_rhs, array(nfunctions)) j = var(name_gen("vdm_j")) @@ -104,7 +107,42 @@ def _emit_func_family_operation(cb, name_gen, cb(coeff_rhs[i], rhs_func(i)) ab_coeffs = var(name_gen("ab_coeffs")) - cb(ab_coeffs, linear_solve(vdmt, coeff_rhs, nfunctions, 1)) + + if hist_len == nfunctions: + cb(ab_coeffs, linear_solve(vdmt, coeff_rhs, nfunctions, 1)) + else: + # Least squares with SVD builtin + u = var(name_gen("u")) + ut = var(name_gen("ut")) + intermed = var(name_gen("intermed")) + ainv = var(name_gen("ainv")) + sigma = var(name_gen("sigma")) + sig_array = var(name_gen("sig_array")) + v = var(name_gen("v")) + vt = var(name_gen("vt")) + + cb(ainv, array(nfunctions*hist_len)) + cb(intermed, array(nfunctions*hist_len)) + + cb((u, sigma, vt), svd(vdmt, hist_len)) + + cb(ut, transpose(u, nfunctions)) + cb(v, transpose(vt, hist_len)) + + # Make singular value array + cb(sig_array, array(nfunctions*nfunctions)) + + for j in range(len(function_family)*len(function_family)): + cb(sig_array[j], 0) + + cb.fence() + + for i in range(len(function_family)): + cb(sig_array[i*(nfunctions+1)], sigma[i]**-1) + + cb(intermed, matmul(v, sig_array, nfunctions, nfunctions)) + cb(ainv, matmul(intermed, ut, nfunctions, nfunctions)) + cb(ab_coeffs, matmul(ainv, coeff_rhs, nfunctions, 1)) return _linear_comb( [ab_coeffs[ii] for ii in range(hist_len)], @@ -127,7 +165,14 @@ def _emit_func_family_operation(cb, name_gen, coeff_rhs[i] = rhs_func(i) - ab_coeffs = la.solve(vdm_t, coeff_rhs) + if hist_len == nfunctions: + ab_coeffs = la.solve(vdm_t, coeff_rhs) + else: + # SVD-based least squares solve + u, sigma, v = la.svd(vdm_t, full_matrices=False) + ainv = np.dot(v.transpose(), np.dot(la.inv(np.diag(sigma)), + u.transpose())) + ab_coeffs = np.dot(ainv, coeff_rhs) return _linear_comb(ab_coeffs, hist_vars) @@ -164,18 +209,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 +302,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 +378,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/leap/multistep/multirate/__init__.py b/leap/multistep/multirate/__init__.py index cbd79da7f38db86f7931edf5ae44736280a007ce..855bf573ff550b491caabf903e6de6b932a35e6a 100644 --- a/leap/multistep/multirate/__init__.py +++ b/leap/multistep/multirate/__init__.py @@ -66,7 +66,8 @@ class MultiRateHistory(Record): .. automethod:: __init__ """ def __init__(self, interval, func_name, arguments, order=None, - rhs_policy=rhs_policy.late, invalidate_computed_state=False): + rhs_policy=rhs_policy.late, invalidate_computed_state=False, + hist_length=None): """ :arg interval: An integer indicating the interval (relative to the smallest available timestep) at which this history is to be @@ -81,6 +82,9 @@ class MultiRateHistory(Record): :arg invalidate_dependent_state: Whether evaluating this right-hand side should force a recomputation of any state that depended upon now-superseded state. + :arg hist_length: history length. If greater than order, we use a + least-squares solve rather than a linear solve to obtain the AB + coefficients for this history """ super(MultiRateHistory, self).__init__( interval=interval, @@ -88,11 +92,12 @@ class MultiRateHistory(Record): arguments=arguments, order=order, rhs_policy=rhs_policy, - invalidate_computed_state=invalidate_computed_state) + invalidate_computed_state=invalidate_computed_state, + hist_length=hist_length) @property def history_length(self): - return self.order + return self.hist_length class RHS(MultiRateHistory): @@ -303,17 +308,21 @@ class MultiRateMultiStepMethod(Method): # }}} - # {{{ plug default order into rhss + # {{{ plug default order, history length into rhss new_rhss = [] for component_rhss in rhss: new_component_rhss = [] for rhs in component_rhss: order = rhs.order + hist_length = rhs.hist_length if order is None: order = default_order + if hist_length is None: + hist_length = order - new_component_rhss.append(rhs.copy(order=order)) + new_component_rhss.append(rhs.copy(order=order, + hist_length=hist_length)) new_rhss.append(tuple(new_component_rhss)) @@ -336,6 +345,10 @@ class MultiRateMultiStepMethod(Method): for component_rhss in self.rhss for rhs in component_rhss) + self.max_hist_length = max(rhs.hist_length + for component_rhss in self.rhss + for rhs in component_rhss) + # {{{ process intervals intervals = sorted(rhs.interval @@ -587,7 +600,7 @@ class MultiRateMultiStepMethod(Method): """Initialize the stepper with an RK method. Return the code that computes the startup history.""" - bootstrap_steps = self.max_order - 1 + bootstrap_steps = self.max_hist_length - 1 final_iglobal_substep = bootstrap_steps * self.nsubsteps @@ -682,12 +695,12 @@ class MultiRateMultiStepMethod(Method): for irhs, rhs in enumerate(component_rhss): if (substeps_from_start % rhs.interval == 0 and (substeps_from_start // rhs.interval - < rhs.order)): + < rhs.history_length)): intervals_from_start = ( substeps_from_start // rhs.interval) - i = rhs.order - 1 - intervals_from_start + i = rhs.history_length - 1 - intervals_from_start assert i >= 0 with cb.if_(self.bootstrap_step, "==", test_step): @@ -740,7 +753,7 @@ class MultiRateMultiStepMethod(Method): key = comp_name, irhs temp_hist_substeps[key] = list(range( - -rhs.interval*(rhs.order-1), 1, rhs.interval)) + -rhs.interval*(rhs.history_length-1), 1, rhs.interval)) if self.static_dt: temp_time_vars[key] = list( @@ -1026,13 +1039,14 @@ class MultiRateMultiStepMethod(Method): if not self.static_dt: for time_var, time_expr in zip( self.time_vars[key], - temp_time_vars[comp_name, irhs][-rhs.order:]): + temp_time_vars[comp_name, + irhs][-rhs.history_length:]): cb(time_var, time_expr) cb.fence() for hist_var, hist_expr in zip( self.history_vars[key], - temp_hist_vars[comp_name, irhs][-rhs.order:]): + temp_hist_vars[comp_name, irhs][-rhs.history_length:]): cb(hist_var, hist_expr) cb.fence() @@ -1126,7 +1140,9 @@ class TwoRateAdamsBashforthMethod(MultiRateMultiStepMethod): def __init__(self, method, order, step_ratio, slow_state_filter_name=None, fast_state_filter_name=None, - static_dt=False): + static_dt=False, + hist_length_slow=None, + hist_length_fast=None): from warnings import warn warn("TwoRateAdamsBashforthMethod is a compatibility shim that should no " "longer be used. Use the fully general " @@ -1158,21 +1174,29 @@ class TwoRateAdamsBashforthMethod(MultiRateMultiStepMethod): else: s2f_policy = rhs_policy.late + if hist_length_slow is None: + hist_length_slow = order + + if hist_length_fast is None: + hist_length_slow = order + super(TwoRateAdamsBashforthMethod, self).__init__( order, ( ( "dt", "fast", "=", - MultiRateHistory(1, "f2f", ("fast", "slow",)), + MultiRateHistory(1, "f2f", ("fast", "slow",), + hist_length=hist_length_fast), MultiRateHistory(s2f_interval, "s2f", - ("fast", "slow",), rhs_policy=s2f_policy), + ("fast", "slow",), rhs_policy=s2f_policy, + hist_length=hist_length_fast), ), ( "dt", "slow", "=", MultiRateHistory(step_ratio, "f2s", ("fast", "slow",), - rhs_policy=f2s_policy), + rhs_policy=f2s_policy, hist_length=hist_length_slow), MultiRateHistory(step_ratio, "s2s", ("fast", "slow",), - rhs_policy=s2s_policy), + rhs_policy=s2s_policy, hist_length=hist_length_slow), ),), state_filter_names={ diff --git a/test/test_ab.py b/test/test_ab.py index 7da554e72847b65b37b561a78f845140fc887a23..f3db445d74c791cbc95bd093adaf2add73fd003b 100755 --- a/test/test_ab.py +++ b/test/test_ab.py @@ -33,15 +33,21 @@ from utils import ( # noqa 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", order, hist_length=order+1, + static_dt=static_dt), order) + for order in [1, 3, 5] + for static_dt in [True, False] + ]) +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) diff --git a/test/test_codegen_fortran.py b/test/test_codegen_fortran.py index d04a51709e581b15a6b723f4658f0053a9843fa7..781b723e539d6a5feca050d479b4b674af0aab32 100755 --- a/test/test_codegen_fortran.py +++ b/test/test_codegen_fortran.py @@ -403,14 +403,15 @@ def test_adaptive_rk_codegen_error(): ]) -@pytest.mark.parametrize("min_order", [2, 3, 4, 5]) -def test_singlerate_squarewave(min_order): +@pytest.mark.parametrize(("min_order", "hist_length"), [(5, 5), (4, 4), (4, 5), + (3, 3), (3, 4), (2, 2), ]) +def test_singlerate_squarewave(min_order, hist_length): from leap.multistep import AdamsBashforthMethod component_id = 'y' rhs_function = 'y' - stepper = AdamsBashforthMethod("y", min_order) + stepper = AdamsBashforthMethod("y", min_order, hist_length=hist_length) from dagrt.function_registry import ( base_function_registry, register_ode_rhs) @@ -438,6 +439,8 @@ def test_singlerate_squarewave(min_order): """) code_str = codegen(code) + with open("abmethod_test.f90", "wt") as outf: + outf.write(code_str) run_fortran([ ("abmethod.f90", code_str), @@ -448,9 +451,11 @@ def test_singlerate_squarewave(min_order): @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) -@pytest.mark.parametrize("min_order", [4, 3, 2]) -def test_multirate_squarewave(min_order, method_name): - stepper = TwoRateAdamsBashforthMethod(method_name, min_order, 4) +@pytest.mark.parametrize(("min_order", "hist_length"), [(4, 4), (4, 5), + (3, 3), (3, 4), (2, 2), ]) +def test_multirate_squarewave(min_order, hist_length, method_name): + stepper = TwoRateAdamsBashforthMethod(method_name, min_order, 4, + hist_length_slow=hist_length, hist_length_fast=hist_length) code = stepper.generate() diff --git a/test/test_multirate.py b/test/test_multirate.py index 07af8e752277f67239d14ac2c05f5308093d1a28..c701620cc86a5aac76ee5c09aa03ff79221bcd94 100644 --- a/test/test_multirate.py +++ b/test/test_multirate.py @@ -46,10 +46,11 @@ from utils import ( # noqa class MultirateTimestepperAccuracyChecker(object): """Check that the multirate timestepper has the advertised accuracy.""" - def __init__(self, method, order, step_ratio, static_dt, ode, method_impl, - display_dag=False, display_solution=False): + def __init__(self, method, order, hist_length, step_ratio, static_dt, ode, + method_impl, display_dag=False, display_solution=False): self.method = method self.order = order + self.hist_length = hist_length self.step_ratio = step_ratio self.static_dt = static_dt self.ode = ode @@ -61,7 +62,9 @@ class MultirateTimestepperAccuracyChecker(object): def get_code(self): method = TwoRateAdamsBashforthMethod( self.method, self.order, self.step_ratio, - static_dt=self.static_dt) + static_dt=self.static_dt, + hist_length_slow=self.hist_length, + hist_length_fast=self.hist_length) return method.generate() @@ -166,7 +169,8 @@ class MultirateTimestepperAccuracyChecker(object): @pytest.mark.slowtest -@pytest.mark.parametrize("order", [1, 3]) +@pytest.mark.parametrize(("order", "hist_length"), [(1, 1), + (3, 3), (3, 4), ]) @pytest.mark.parametrize("system", [ #"Basic", "Full", @@ -175,7 +179,8 @@ class MultirateTimestepperAccuracyChecker(object): ]) @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) @pytest.mark.parametrize("static_dt", [True, False]) -def test_multirate_accuracy(method_name, order, system, static_dt, step_ratio=2): +def test_multirate_accuracy(method_name, order, hist_length, + system, static_dt, step_ratio=2): """Check that the multirate timestepper has the advertised accuracy""" import multirate_test_systems @@ -187,12 +192,12 @@ def test_multirate_accuracy(method_name, order, system, static_dt, step_ratio=2) print("------------------------------------------------------") MultirateTimestepperAccuracyChecker( - method_name, order, step_ratio, static_dt=static_dt, + method_name, order, hist_length, step_ratio, static_dt=static_dt, ode=system(), method_impl=pmi_cg)() -def test_single_rate_identical(order=3): +def test_single_rate_identical(order=3, hist_length=3): from leap.multistep import AdamsBashforthMethod from dagrt.exec_numpy import NumpyInterpreter @@ -204,7 +209,8 @@ def test_single_rate_identical(order=3): # {{{ single rate - single_rate_method = AdamsBashforthMethod("y", order=order) + single_rate_method = AdamsBashforthMethod("y", order=order, + hist_length=hist_length) single_rate_code = single_rate_method.generate() def single_rate_rhs(t, y): @@ -244,12 +250,13 @@ def test_single_rate_identical(order=3): ( ( 'dt', 'fast', '=', - MRHistory(1, "f", ("fast", "slow",)), + MRHistory(1, "f", ("fast", "slow",), + hist_length=hist_length), ), ( 'dt', 'slow', '=', MRHistory(1, "s", ("fast", "slow",), - rhs_policy=rhs_policy.late), + rhs_policy=rhs_policy.late, hist_length=hist_length), ),) ) multi_rate_code = multi_rate_method.generate()