diff --git a/leap/multistep/multirate/__init__.py b/leap/multistep/multirate/__init__.py index cbd79da7f38db86f7931edf5ae44736280a007ce..578184281cf499ac5c5325061873087276b36c1e 100644 --- a/leap/multistep/multirate/__init__.py +++ b/leap/multistep/multirate/__init__.py @@ -146,6 +146,10 @@ def _topologically_sort_comp_names_and_rhss(component_names, rhss): # }}} +class InconsistentHistoryError: + pass + + # {{{ method class MultiRateMultiStepMethod(Method): @@ -166,7 +170,10 @@ class MultiRateMultiStepMethod(Method): def __init__(self, default_order, system_description, state_filter_names=None, component_arg_names=None, - static_dt=False): + static_dt=False, + hist_consistency_threshold=None, + early_hist_consistency_threshold=None): + """ :arg default_order: The order to be used for right-hand sides where no differing order is specified. @@ -363,6 +370,12 @@ class MultiRateMultiStepMethod(Method): # }}} self.static_dt = static_dt + self.hist_consistency_threshold = hist_consistency_threshold + if isinstance(early_hist_consistency_threshold, str): + from dagrt.expression import parse + early_hist_consistency_threshold = parse( + early_hist_consistency_threshold) + self.early_hist_consistency_threshold = early_hist_consistency_threshold if not self.static_dt: self.time_vars = {} @@ -594,7 +607,7 @@ class MultiRateMultiStepMethod(Method): from pytools import UniqueNameGenerator name_gen = UniqueNameGenerator() - for isubstep in range(self.nsubsteps + 1): + for isubstep in range(self.nsubsteps): name_prefix = 'substep' + str(isubstep) current_rhss = {} @@ -699,19 +712,20 @@ class MultiRateMultiStepMethod(Method): # }}} - if isubstep == self.nsubsteps: - cb.fence() - cb(self.bootstrap_step, self.bootstrap_step + 1) - break + cb.fence() if isubstep == 0: with cb.if_(self.bootstrap_step, "==", bootstrap_steps): cb.state_transition("primary") + cb.exit_step() cb.fence() self.emit_small_rk_step(cb, name_prefix, name_gen, current_rhss) + cb.fence() + cb(self.bootstrap_step, self.bootstrap_step + 1) + return cb # }}} @@ -968,9 +982,70 @@ class MultiRateMultiStepMethod(Method): # }}} + def norm(expr): + return var('norm_2')(expr) + + def check_history_consistency(): + # At the start of a macrostep, ensure that the last computed + # RHS history corresponds to the current state + for comp_idx, (comp_name, component_rhss) in enumerate( + zip(self.component_names, self.rhss)): + for irhs, rhs in enumerate(component_rhss): + t_expr = self.t + kwargs = dict( + (self.comp_name_to_kwarg_name[arg_comp_name], + get_state(arg_comp_name, 0)) + for arg_comp_name in rhs.arguments) + test_rhs_var = var( + name_gen( + "test_rhs_{comp_name}_rhs{irhs}_0" + .format(comp_name=comp_name, irhs=irhs))) + + cb(test_rhs_var, var(rhs.func_name)(t=t_expr, **kwargs)) + # Compare this computed RHS with the 0th history point using + # built-in norm. + + zeroth_hist = temp_hist_vars[comp_name, irhs][-1] + rel_rhs_error = ( + norm(test_rhs_var - zeroth_hist) + / + norm(test_rhs_var)) + + cb("rel_rhs_error", rel_rhs_error) + + # cb((), "print(rel_rhs_error)") + + if rhs.rhs_policy == rhs_policy.early: + # Check for scheme-order accuracy + if self.early_hist_consistency_threshold is not None: + with cb.if_("rel_rhs_error", ">=", + self.early_hist_consistency_threshold): + cb((), "print(rel_rhs_error)") + cb.raise_(InconsistentHistoryError, + "MRAB: top-of-history for RHS '%s' is not " + "consistent with current state" + % rhs.func_name) + else: + cb.raise_(InconsistentHistoryError, + "MRAB: RHS '%s' has early policy " + "and requires relaxed threshold input" + % rhs.func_name) + + else: + # Check for floating-point accuracy + with cb.if_("rel_rhs_error", ">=", + self.hist_consistency_threshold): + cb.raise_(InconsistentHistoryError, + "MRAB: top-of-history for RHS '%s' is not " + "consistent with current state" % rhs.func_name) + # {{{ run_substep_loop def run_substep_loop(): + # Check last history value from previous macrostep + if self.hist_consistency_threshold is not None: + check_history_consistency() + for isubstep in range(self.nsubsteps+1): for comp_idx, (comp_name, component_rhss) in enumerate( zip(self.component_names, self.rhss)): @@ -1126,7 +1201,8 @@ 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_consistency_threshold=None, + early_hist_consistency_threshold=None): from warnings import warn warn("TwoRateAdamsBashforthMethod is a compatibility shim that should no " "longer be used. Use the fully general " @@ -1183,7 +1259,10 @@ class TwoRateAdamsBashforthMethod(MultiRateMultiStepMethod): # This is a hack to avoid having to change the 2RAB test # cases, which use these arguments component_arg_names=("f", "s"), - static_dt=static_dt) + + static_dt=static_dt, + hist_consistency_threshold=hist_consistency_threshold, + early_hist_consistency_threshold=early_hist_consistency_threshold) # }}} @@ -1295,5 +1374,4 @@ class TextualSchemeExplainer(SchemeExplainerBase): # }}} - # vim: foldmethod=marker diff --git a/requirements.txt b/requirements.txt index f8db6acc7c93693eabacbc94d95fcc215257a361..e22f333cbd530c0259a7df4fcb84ac4496c7482c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ git+git://github.com/inducer/pytools git+git://github.com/inducer/pymbolic -git+https://gitlab.tiker.net/inducer/dagrt.git +git+https://gitlab.tiker.net/inducer/dagrt.git \ No newline at end of file diff --git a/test/test_codegen_fortran.py b/test/test_codegen_fortran.py index d04a51709e581b15a6b723f4658f0053a9843fa7..f884b1b76f6be03aac61ec3d7e9258d541e06747 100755 --- a/test/test_codegen_fortran.py +++ b/test/test_codegen_fortran.py @@ -223,7 +223,23 @@ def test_multirate_codegen(min_order, method_name): stepper = TwoRateAdamsBashforthMethod( method_name, min_order, 4, slow_state_filter_name="slow_filt", - fast_state_filter_name="fast_filt") + fast_state_filter_name="fast_filt", + # should pass with either, let's alternate by order + # static_dt=True is 'more finnicky', so use that at min_order=5. + static_dt=True if min_order % 2 == 1 else False, + hist_consistency_threshold=1e-8, + early_hist_consistency_threshold="1.25e3 *
**%d" % min_order) + + # Early consistency threshold checked for convergence + # with timestep change - C. Mikida, 2/6/18 (commit hash 2e6ca077) + + # With method 5-Fqs (limiting case), the following maximum relative + # errors were observed: + # for dt = 0.0384: 1.03E-04 + # for dt = 0.0128: 5.43E-08 + + # Reported relative errors motivate constant factor of 1.25e3 for early + # consistency threshold code = stepper.generate() @@ -304,18 +320,12 @@ def test_multirate_codegen(min_order, method_name): with open("abmethod.f90", "wt") as outf: outf.write(code_str) - if (method_name.startswith("S") - and "f" in method_name): - fac = 12 - elif (method_name.startswith("S") - or method_name.startswith("Fs") - or method_name == "F"): - fac = 9 - else: - fac = 5 + fac = 130 + if min_order == 5 and method_name in ["Srs", "Ss"]: + pytest.xfail("Srs,Ss do not achieve fifth order convergence") num_trips_one = 10*fac - num_trips_two = 100*fac + num_trips_two = 30*fac run_fortran([ ("abmethod.f90", code_str), @@ -450,7 +460,21 @@ 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) + stepper = TwoRateAdamsBashforthMethod(method_name, min_order, 4, + hist_consistency_threshold=1e-8, + early_hist_consistency_threshold="3.0e3 *
**%d" % min_order) + + # Early consistency threshold checked for convergence + # with timestep change - C. Mikida, 2/6/18 (commit hash 2e6ca077) + + # With method 4-Ss (limiting case), the following maximum relative + # errors were observed: + # for dt = 0.009167: 5.55E-08 + # for dt = 0.00611: 1.59E-08 + # Corresponding EOC: 3.08 + + # Reported relative errors motivate constant factor of 3.0e3 for early + # consistency threshold code = stepper.generate() @@ -515,6 +539,8 @@ def test_multirate_squarewave(min_order, method_name): if min_order == 2: fac = 200 + elif min_order == 3: + fac = 100 else: fac = 12 num_trips_one = 100*fac diff --git a/test/test_mrab.f90 b/test/test_mrab.f90 index 59d1992927b4b5adc7803928bef67fb8dd535efd..763e4d078f60cb5090663e6af0213060ce92adba 100644 --- a/test/test_mrab.f90 +++ b/test/test_mrab.f90 @@ -15,7 +15,7 @@ program test_mrabmethod integer run_count real*8 t_fin - parameter (run_count=2, t_fin=10d-0) + parameter (run_count=2, t_fin=50d-0) real*8, dimension(run_count):: dt_values, error_slow, error_fast diff --git a/test/test_multirate.py b/test/test_multirate.py index 07af8e752277f67239d14ac2c05f5308093d1a28..ac0e7e4e1de45d0381e8519615ef75f206c64a38 100644 --- a/test/test_multirate.py +++ b/test/test_multirate.py @@ -58,10 +58,24 @@ class MultirateTimestepperAccuracyChecker(object): self.display_solution = display_solution @memoize_method - def get_code(self): + def get_code(self, dt): method = TwoRateAdamsBashforthMethod( self.method, self.order, self.step_ratio, - static_dt=self.static_dt) + static_dt=self.static_dt, + hist_consistency_threshold=1e-8, + early_hist_consistency_threshold=dt**self.order) + + # Early consistency threshold checked for convergence + # with timestep change - C. Mikida, 2/6/18 (commit hash 2e6ca077) + + # With method 4-Ss (limiting case), the following maximum relative + # errors were observed: + # for dt = 0.015625: 3.11E-09 + # for dt = 0.0078125: 1.04e-10 + # Corresponding EOC: 4.90 + + # Reported relative errors show that no constant factor is needed on + # early consistency threshold return method.generate() @@ -80,8 +94,8 @@ class MultirateTimestepperAccuracyChecker(object): self.ode.f2f_rhs, self.ode.s2f_rhs, self.ode.f2s_rhs, self.ode.s2s_rhs)} - print(self.get_code()) - method = self.method_impl(self.get_code(), function_map=function_map) + print(self.get_code(dt)) + method = self.method_impl(self.get_code(dt), function_map=function_map) t = self.ode.t_start y = self.ode.initial_values @@ -250,8 +264,10 @@ def test_single_rate_identical(order=3): 'dt', 'slow', '=', MRHistory(1, "s", ("fast", "slow",), rhs_policy=rhs_policy.late), - ),) - ) + ),), + hist_consistency_threshold=1e-8, + early_hist_consistency_threshold=dt**order) + multi_rate_code = multi_rate_method.generate() def rhs_fast(t, fast, slow): @@ -301,6 +317,7 @@ def test_2rab_scheme_explainers(method_name, order=3, step_ratio=3, explainer=TextualSchemeExplainer()): method = TwoRateAdamsBashforthMethod( method_name, order=order, step_ratio=step_ratio) + method.generate(explainer=explainer) print(explainer) @@ -355,7 +372,8 @@ def test_mrab_with_derived_state_scheme_explainers(order=3, step_ratio=3, def test_dot(order=3, step_ratio=3, method_name="F", show=False): method = TwoRateAdamsBashforthMethod( - method_name, order=order, step_ratio=step_ratio) + method_name, order=order, step_ratio=step_ratio, + hist_consistency_threshold=1e-8) code = method.generate() from dagrt.language import get_dot_dependency_graph