From 361352e5edcef2fd0e9a9ca7fdc884e903c2c35e Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sat, 4 May 2019 00:43:51 -0500 Subject: [PATCH 1/5] Add examples from paper draft --- examples/adaptive-rk/adaptive-rk.py | 164 +++++ examples/semi-implicit-rk/semi-implicit-rk.py | 121 ++++ .../wave-equation.py | 674 ++++++++++++++++++ 3 files changed, 959 insertions(+) create mode 100755 examples/adaptive-rk/adaptive-rk.py create mode 100755 examples/semi-implicit-rk/semi-implicit-rk.py create mode 100755 examples/variable-coeff-wave-equation/wave-equation.py diff --git a/examples/adaptive-rk/adaptive-rk.py b/examples/adaptive-rk/adaptive-rk.py new file mode 100755 index 0000000..792b108 --- /dev/null +++ b/examples/adaptive-rk/adaptive-rk.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python3 +"""This example demonstrates direct construction of an adaptive RK method using +an order 2/3 method pair. + +Source: + + Bogacki, Przemyslaw; Shampine, Lawrence F. (1989), + "A 3(2) pair of Runge–Kutta formulas", + Applied Mathematics Letters, 2 (4): 321–325, + doi:10.1016/0893-9659(89)90079-7. + +The same functionality is supported using leap.rk.ODE23Method. + +""" + +import logging +import numpy as np + + +logging.basicConfig(level=logging.INFO) + + +class KapsProblem(object): + """ + From Kennedy and Carpenter, Section 7.1 + + y_1' = - (epsilon^{-1} + 2) y_1 + epsilon^{-1} y_2^2 + y_2' = y_1 - y_2 - y_2^2 + + 0 <= t <= 1 + + The initial conditions are + + y_1 = y_2 = 1. + + The exact solution is + + y_1 = exp(-2t) + y_2 = exp(-t). + + The stiff component are the terms multiplied by epsilon^{-1}. + """ + + def __init__(self, epsilon): + self._epsilon_inv = 1 / epsilon + self.t_start = 0 + self.t_end = 1 + + def initial(self): + return np.array([1., 1.]) + + def nonstiff(self, t, y): + y_1 = y[0] + y_2 = y[1] + return np.array([-2 * y_1, y_1 - y_2 - y_2 ** 2]) + + def stiff(self, t, y): + y_1 = y[0] + y_2 = y[1] + return np.array([-self._epsilon_inv * (y_1 - y_2 ** 2), 0]) + + def jacobian(self, t, y): + y_2 = y[1] + return np.array([ + [-self._epsilon_inv - 2, 2 * self._epsilon_inv * y_2], + [1, -1 - 2 * y_2]]) + + def exact(self, t): + return np.array([np.exp(-2 * t), np.exp(-t)]) + + def __call__(self, t, y): + return self.stiff(t, y) + self.nonstiff(t, y) + + +_default_dts = 2 ** -np.array(range(5, 10), dtype=np.float64) + + +def get_convergence_data(method_class, problem, dts=_default_dts): + from pytools.convergence import EOCRecorder + + eocrec = EOCRecorder() + + component_id = "y" + + for dt in dts: + t = problem.t_start + y = problem.initial() + final_t = problem.t_end + + interp = method_class(function_map={ + "f": problem, + }) + interp.set_up(t_start=t, dt_start=dt, context={component_id: y}) + + times = [] + values = [] + for event in interp.run(t_end=final_t): + if isinstance(event, interp.StateComputed): + assert event.component_id == component_id + values.append(event.state_component[0]) + times.append(event.t) + + assert abs(times[-1] - final_t) / final_t < 0.1 + times = np.array(times) + + error = abs(values[-1] - problem.exact(final_t)[0]) + eocrec.add_data_point(dt, error) + + return eocrec + + +def demo_rk_adaptive(): + from dagrt.language import CodeBuilder + from pymbolic import var + + # Set tolerance + tol = 1e-3 + # Bulk declare symbolic names + k1, k2, k3, k4, t, dt, dt_old, y, y_hi, y_lo, f, norm = ( + var(name) for name in + "k1 k2 k3 k4
dt_old y y_hi y_lo f norm_inf".split()) # noqa + + with CodeBuilder() as cb: + # Calculate the RK stage values + cb(k1, f(t, y)) + cb(k2, f(t + 1/2 * dt, y + dt * (1/2 * k1))) + cb(k3, f(t + 3/4 * dt, y + dt * (3/4 * k2))) + cb(k4, f(t + dt, y + dt * (2/9 * k1 + 1/3 * k2 + 4/9 * k3))) + # Compute the low and high order solutions + cb(y_lo, y + dt * (7/24 * k1 + 1/4 * k2 + 1/3 * k3 + 1/8 * k4)) + cb(y_hi, y + dt * (2/9 * k1 + 1/3 * k2 + 4/9 * k3)) + # Save the value of dt + cb(dt_old, dt) + # Update dt based on the error estimate + err_est = norm(y_lo - y_hi) + order = 3 + cb.reset_dep_tracking() + cb(dt, 0.9 * dt * (tol / err_est) ** (1 / order)) + # Adapt the step size + with cb.if_(err_est, "<=", tol): + # Update t and y + cb(y, y_hi) + cb(t, t + dt_old) + cb.yield_state(expression=y, component_id="y", time=t, time_id=None) # noqa + with cb.else_(): + cb.fail_step() + + from dagrt.language import DAGCode + code = DAGCode(phases={ + "primary": cb.as_execution_phase(next_phase="primary") + }, + initial_phase="primary") + + print(code) + + # Generate and run the method. + from dagrt.codegen import PythonCodeGenerator + RKAdaptiveMethod = PythonCodeGenerator("RKAdaptiveMethod").get_class(code) + eocrec = get_convergence_data(RKAdaptiveMethod, problem=KapsProblem(0.001)) + print(eocrec.pretty_print()) + + +if __name__ == "__main__": + demo_rk_adaptive() diff --git a/examples/semi-implicit-rk/semi-implicit-rk.py b/examples/semi-implicit-rk/semi-implicit-rk.py new file mode 100755 index 0000000..e778227 --- /dev/null +++ b/examples/semi-implicit-rk/semi-implicit-rk.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 +"""This example demonstrates direct construction of a semi-implicit RK method. + +Source: + + Kennedy, Christopher A., and Mark H. Carpenter. + "Diagonally implicit runge-kutta + methods for ordinary differential equations. a review." (2016). + +""" + +import numpy as np +import logging + + +logging.basicConfig(level=logging.INFO) + + +class SimpleDecayProblem(object): + + t_start = 0 + t_end = 0.5 + + def initial(self): + return np.array([1]) + + def exact(self, t): + return np.array([np.exp(-10 * t)]) + + def __call__(self, t, y): + return -10 * y + + +_default_dts = 2 ** -np.array(range(5, 10), dtype=np.float64) + + +def get_convergence_data(method_class, problem, dts=_default_dts): + from pytools.convergence import EOCRecorder + + eocrec = EOCRecorder() + + component_id = "y" + + for dt in dts: + t = problem.t_start + y = problem.initial() + final_t = problem.t_end + + interp = method_class(function_map={ + "f": problem, + }) + interp.set_up(t_start=t, dt_start=dt, context={component_id: y}) + + times = [] + values = [] + for event in interp.run(t_end=final_t): + if isinstance(event, interp.StateComputed): + assert event.component_id == component_id + values.append(event.state_component[0]) + times.append(event.t) + + assert abs(times[-1] - final_t) / final_t < 0.1 + times = np.array(times) + + error = abs(values[-1] - problem.exact(final_t)[0]) + eocrec.add_data_point(dt, error) + + return eocrec + + +def demo_rk_implicit(): + from dagrt.language import CodeBuilder + from pymbolic import var + + k1, k2, t, dt, y, f = ( + var(name) for name in + "k1 k2
y f".split()) + + gamma = (2 - 2**0.5) / 2 + + with CodeBuilder() as cb: + cb.assign_implicit_1( + k1, + solve_component=k1, + expression=k1 - f(t + gamma * dt, y + dt * gamma * k1), + guess=f(t, y)) + cb.assign_implicit_1( + k2, + solve_component=k2, + expression=k2 - f(t + dt, y + dt * ((1 - gamma) * k1 + gamma * k2)), # noqa + guess=k1) + cb.reset_dep_tracking() + cb(y, y + dt * ((1 - gamma) * k1 + gamma * k2)) + cb(t, t + dt) + cb.yield_state(y, "y", t, None) + + from dagrt.language import DAGCode + code = DAGCode(phases={ + "primary": cb.as_execution_phase(next_phase="primary") + }, + initial_phase="primary") + + def solver_hook(solve_expr, unknown, solver_id, guess): + from dagrt.expression import match, substitute + pieces = match( + "k - rhs(time, y + dt * (c0 + c1 * k))", + solve_expr, + pre_match={"k": unknown}) + return substitute("-10 * (dt * c0 + y) / (10 * dt * c1 + 1)", pieces) + + from leap.implicit import replace_AssignImplicit + code = replace_AssignImplicit(code, solver_hook) + + from dagrt.codegen import PythonCodeGenerator + IRKMethod = PythonCodeGenerator("IRKMethod").get_class(code) + eocrec = get_convergence_data(IRKMethod, SimpleDecayProblem()) + print(eocrec.pretty_print()) + + +if __name__ == "__main__": + demo_rk_implicit() diff --git a/examples/variable-coeff-wave-equation/wave-equation.py b/examples/variable-coeff-wave-equation/wave-equation.py new file mode 100755 index 0000000..b4880c6 --- /dev/null +++ b/examples/variable-coeff-wave-equation/wave-equation.py @@ -0,0 +1,674 @@ +#!/usr/bin/env python +"""Solves the 1D wave equation + + u_tt = c(x) u_xx + u_t(0) = init' + u(0) = init + +with piecewise constant coefficients c(x) using a multirate multistep method. +""" + + +import logging +import numpy as np +import numpy.linalg as la + +import matplotlib + +from contextlib import contextmanager + + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +PAPER_OUTPUT = False + + +# {{{ matplotlib setup + +if PAPER_OUTPUT: + matplotlib.use("pgf") + + +import matplotlib.pyplot as plt # noqa + + +FONTSIZE = 9 +LINEWIDTH = 0.5 + + +def initialize_matplotlib(): + plt.rc("font", family="serif") + plt.rc("text", usetex=False) + plt.rc("xtick", labelsize=FONTSIZE) + plt.rc("ytick", labelsize=FONTSIZE) + plt.rc("axes", labelsize=FONTSIZE) + plt.rc("axes", titlesize=FONTSIZE) + plt.rc("axes", linewidth=LINEWIDTH) + plt.rc("pgf", rcfonts=True) + plt.rc("lines", linewidth=LINEWIDTH) + + +if PAPER_OUTPUT: + initialize_matplotlib() + +# }}} + + +# {{{ problem discretization + +class VariableCoeffWaveEquationProblem(object): + """Discretizes an instance of the 1D variable coefficient wave equation + problem, providing initial values and right-hand side evaluations. + + This turns the problem instance into the 1D system + + u_t = v + v_t = c(x) * u_xx. + + The discretization uses a centered difference stencil. Interface conditions + are imposed using a second-order biased stencil. + + The system state is separated into components, i.e. regions in which c(x) + is constant. Components are equispaced in the interval (0, 1) and + correspond to the number of coefficients passed to the constructor. + + """ + + def __init__(self, ngridpoints, component_coeffs, decay=0, diffusion=0): + """ + Args: + ngridpoints: Number of mesh points + component_coeffs: Coefficients of components + (=squares of wave speeds) + """ + + ncomponents = len(component_coeffs) + + grid = np.linspace(0, 1, ngridpoints + 1, endpoint=False) + component_times = np.r_[grid[::1+ngridpoints//ncomponents], 1] + + # Left endpoint of grid is fixed at 0. + self.grid = grid[1:] + self.h = grid[1] - grid[0] + self.ncomponents = ncomponents + + # n = Number of simulation unknowns + # + # Interface points (those between adjacent components) are not + # included: the solution value at the interface points is inferred + # using the stencil. + self.n = 2 * (ngridpoints - self.ncomponents + 1) + + component_indices = [] + component_sizes = [] + + for i, (start, stop) in enumerate( + zip(component_times, component_times[1:])): + indices = (start < grid) & (grid < stop) + component_indices.append(indices) + component_sizes.append(2 * sum(indices)) + + assert sum(component_sizes) == self.n + + self.component_indices = component_indices + self.component_sizes = component_sizes + self.coeffs = component_coeffs + + self.decay = decay + self.diffusion = diffusion + + # Initial values + f = np.sin(2 * 6 * np.pi * grid) + f_prime = np.zeros_like(f) + + # Initial conditions + self.u0_by_component = [] + + for index in component_indices: + self.u0_by_component.append(np.hstack([f[index], f_prime[index]])) + + self.u0 = np.hstack(self.u0_by_component) + + @property + def dt_ref(self): + return (1/2) * self.h / min(self.coeffs) + + def split_view(self, comp): + return comp[:len(comp) // 2], comp[len(comp) // 2:] + + def rhs(self, from_component_num, to_component_num, t, **kwargs): + """Evaluate the right-hand side for the given component interaction. + + Args: + from_component_num: Source component + to_component_num: Target component + kwargs: Component values by name. "compN" is the N-th component + """ + component = kwargs["comp%d" % from_component_num] + old_u, old_ut = self.split_view(component) + + result = np.zeros( + self.component_sizes[to_component_num], + dtype=component.dtype) + res_ut, res_utt = self.split_view(result) + + # Contribution from left panel + if from_component_num == to_component_num - 1: + # Interface conditions + r = self.coeffs[from_component_num] / (3 * ( + self.coeffs[from_component_num] + + self.coeffs[to_component_num])) + u_l = r * (4 * old_u[-1] - old_u[-2]) + + res_ut[0] = self.diffusion * u_l / self.h ** 2 + res_utt[0] = self.coeffs[to_component_num] * u_l / self.h ** 2 + + # Contribution from self panel + elif from_component_num == to_component_num: + old_uxx_mid = ( + (old_u[:-2] - 2 * old_u[1:-1] + old_u[2:]) + / self.h ** 2) + old_uxx_l = (-2 * old_u[0] + old_u[1]) / self.h ** 2 + old_uxx_r = (old_u[-2] - 2 * old_u[-1]) / self.h ** 2 + + if from_component_num > 0: + # Left interface condition + r = self.coeffs[from_component_num] / (3 * ( + self.coeffs[from_component_num - 1] + + self.coeffs[from_component_num])) + + old_u_l = r * (4 * old_u[0] - old_u[1]) + old_uxx_l += old_u_l / self.h ** 2 + + if from_component_num < self.ncomponents - 1: + # Right interface condition + r = self.coeffs[from_component_num] / (3 * ( + self.coeffs[from_component_num + 1] + + self.coeffs[from_component_num])) + + old_u_r = r * (4 * old_u[-1] - old_u[-2]) + old_uxx_r += old_u_r / self.h ** 2 + + old_uxx = np.hstack([[old_uxx_l], old_uxx_mid, [old_uxx_r]]) + res_ut[:] = old_ut - self.decay * old_u + self.diffusion * old_uxx + res_utt[:] = self.coeffs[from_component_num] * old_uxx + + # Contribution from right panel + elif from_component_num == to_component_num + 1: + # Interface condition + r = self.coeffs[from_component_num] / (3 * ( + self.coeffs[from_component_num] + + self.coeffs[to_component_num])) + u_r = r * (4 * old_u[0] - old_u[1]) + + res_ut[-1] = self.diffusion * u_r / self.h ** 2 + res_utt[-1] = self.coeffs[to_component_num] * u_r / self.h ** 2 + + else: + raise ValueError("Components not spatially adjacent") + + return result + + def full_solution(self, components): + """Reconstruct a full solution of the problem from a list of solution + components. + + """ + result = [] + + for i in range(self.ncomponents): + u = self.split_view(components[i])[0] + result.append(u) + + if i < self.ncomponents - 1: + u_next, _ = self.split_view(components[i + 1]) + # Reconstruct midpoint from interface condition. + r = 1 / (3 * ( + self.coeffs[i] + self.coeffs[i + 1])) + u_l = self.coeffs[i] * r * (4 * u[-1] - u[-2]) + u_r = self.coeffs[i+1] * r * (4 * u_next[0] - u_next[1]) + result.append([u_l + u_r]) + + result = np.hstack(result) + return result + +# }}} + + +# {{{ multirate method generator + +def make_3_component_multirate_method( + problem, ratios, order=3, code_only=False, return_rhs_map=False): + """Return the object that drives the multirate method for the given + parameters. + + """ + from leap.multistep.multirate import ( + MultiRateMultiStepMethod, MultiRateHistory) + + ncomponents = 3 + assert problem.ncomponents == ncomponents + + code = MultiRateMultiStepMethod( + default_order=3, + system_description=( + ("dt", "comp0", + "=", + MultiRateHistory(ratios[0], "rhs0to0", ("comp0",)), + MultiRateHistory(ratios[0], "rhs1to0", ("comp1",))), + ("dt", "comp1", + "=", + MultiRateHistory(ratios[1], "rhs0to1", ("comp0",)), + MultiRateHistory(ratios[1], "rhs1to1", ("comp1",)), + MultiRateHistory(ratios[1], "rhs2to1", ("comp2",))), + ("dt", "comp2", + "=", + MultiRateHistory(ratios[2], "rhs1to2", ("comp1",)), + MultiRateHistory(ratios[2], "rhs2to2", ("comp2",)))), + static_dt=True).generate() + + from functools import partial + + if code_only and not return_rhs_map: + return code + + rhs_map = {} + for i in range(3): + rhs_map["rhs%dto%d" % (i, i)] = partial(problem.rhs, i, i) + if i > 0: + rhs_map["rhs%dto%d" % (i, i - 1)] = ( + partial(problem.rhs, i, i - 1)) + if i < ncomponents - 1: + rhs_map["rhs%dto%d" % (i, i + 1)] = ( + partial(problem.rhs, i, i + 1)) + + if code_only: + if return_rhs_map: + return code, rhs_map + else: + return code + + from dagrt.codegen import PythonCodeGenerator + MRABMethod = PythonCodeGenerator(class_name='MRABMethod').get_class(code) # noqa + + if return_rhs_map: + return MRABMethod(rhs_map), rhs_map + else: + return MRABMethod(rhs_map) + +# }}} + + +# {{{ example plotter + +def plot_example(ngridpoints): + problem = VariableCoeffWaveEquationProblem( + ngridpoints, component_coeffs=[16, 4, 1]) + stepper = make_3_component_multirate_method(problem, ratios=(1, 1, 1)) + + dt = 0.01 * problem.dt_ref + t_end = 0.5 + + # {{{ step solution + + stepper.set_up( + t_start=0, + dt_start=dt, + context={ + "comp0": problem.u0_by_component[0], + "comp1": problem.u0_by_component[1], + "comp2": problem.u0_by_component[2]}) + + vals_by_component = {} + times_by_component = {} + + for event in stepper.run(t_end=t_end): + if isinstance(event, stepper.StateComputed): + vals_by_component.setdefault( + event.component_id, []).append(event.state_component) + times_by_component.setdefault( + event.component_id, []).append(event.t) + + n = len(vals_by_component["comp0"]) + vals = [] + for i in range(n): + soln = problem.full_solution([ + vals_by_component[comp][i] + for comp in ("comp0", "comp1", "comp2")]) + # Add Dirichlet BCs for proper axis labeling + vals.append(np.r_[0, soln, 0]) + + vals = np.vstack(vals) + + # }}} + + figure = plt.figure(figsize=(4, 2), dpi=300) + axis = figure.add_subplot(111) + + image = axis.imshow( + vals, + cmap="Spectral", + aspect="auto", + origin="lower", + extent=(0, 1, 0, 0.5), + vmin=-1.5, + vmax=1.5) + + from mpl_toolkits.axes_grid1 import make_axes_locatable + divider = make_axes_locatable(axis) + cax = divider.append_axes("right", size="2%", pad=0.1) + + plt.colorbar(image, cax=cax) + + """ + axis.plot_surface( + problem.grid.reshape(-1, 1), + np.array(times_by_component["comp0"]).reshape(1, -1), + vals.T, + rcount=100, + ccount=100, + cmap="Spectral") + """ + + axis.set_xlabel("$x$") + axis.set_ylabel("$t$") + axis.set_aspect(1) + axis.set_xticks([0, 1/3, 2/3, 1]) + axis.set_xticklabels(["0", "1/3", "2/3", "1"]) + axis.set_title("Solution to 1D Wave Equation with Variable Coefficients") + + suffix = "pgf" if PAPER_OUTPUT else "pdf" + filename = f"wave-problem.{suffix}" + + plt.savefig(filename, bbox_inches="tight") + logger.info("wrote to '%s'" % filename) + +# }}} + + +# {{{ stability experiment + +@contextmanager +def timer(name): + from time import time + logging.info("start: {name}".format(name=name)) + start = time() + yield + end = time() + logging.info("finished: {name}: {time} seconds".format( + name=name, time=end - start)) + + +def generate_mrab_step_matrix(ngridpoints, coeffs, substep_ratio): + problem = VariableCoeffWaveEquationProblem(ngridpoints, coeffs) + + code, rhs_map = make_3_component_multirate_method( + problem, substep_ratio, code_only=True, return_rhs_map=True) + from leap.step_matrix import StepMatrixFinder + + finder = StepMatrixFinder( + code, + function_map=rhs_map, + exclude_variables=["

bootstrap_step"]) + + with timer("Constructing MRAB({}) matrix".format(substep_ratio)): + component_sizes = {} + for var in finder.variables: + for i in range(problem.ncomponents): + if f"comp{i}" in var: + component_sizes[var] = problem.component_sizes[i] + break + else: + raise ValueError(f"cannot infer size of variable: {var}") + + import pprint + pprint.pprint(component_sizes) + + mat = finder.get_phase_step_matrix( + "primary", + shapes=component_sizes, + sparse=True) + + substep_ratio_suffix = "-".join(str(r) for r in substep_ratio) + + filename = f"mat{ngridpoints}-{substep_ratio_suffix}.pkl" + + with open(filename, "wb") as outf: + import pickle + pickle.dump(mat, outf) + + logging.info(f"{len(mat.data)} nnz, size {mat.shape}") + + return filename + + +def compute_all_stable_timesteps(filenames, stable_dts_outf): + import re + + rows = [["Intervals", r"Stable $\Delta t$", + r"$\Delta t / \Delta t_{(1,1,1)}$"]] + + first = True + + for fname in filenames: + import pickle + + with timer("loading {}".format(fname)): + with open(fname, "rb") as infile: + mat = pickle.load(infile) + + logging.info("Computing stable timestep") + dt = compute_stable_timestep(mat) + + if first: + dt_1 = dt + + intervals = tuple( + int(i) for i in + re.match(r"mat\d+-(\d+)-(\d+)-(\d+)\.pkl", fname).groups()) + + row = [str(intervals), f"{dt:.2e}"] + + if first: + row.append("---") + else: + ratio = dt / dt_1 + row.append(f"{ratio:.1f}") + + first = False + rows.append(row) + + print(tabulate(rows, col_fmt="cSc"), file=stable_dts_outf) + + +def compute_stable_timestep(step_matrix, tol=0, prec=1e-15): + def as_dense(mat): + from scipy.sparse import coo_matrix + indices = np.asarray(mat.indices) + return coo_matrix( + (mat.data, (indices[:, 0], indices[:, 1])), + shape=mat.shape).todense() + + from leap.step_matrix import fast_evaluator + evaluate_mat = fast_evaluator(step_matrix, sparse=True) + + def spectral_radius(dt): + mat = as_dense(evaluate_mat({"

": dt})) + eigvals = la.eigvals(mat) + radius = np.max(np.abs(eigvals)) + logging.info(f"{dt} -> spectral radius {radius}") + return radius + + def is_stable(dt): + return spectral_radius(dt) <= 1 + + from leap.stability import find_truth_bdry + return find_truth_bdry(is_stable, prec=1e-7, start_magnitude=1e-4) + +# }}} + + +# {{{ accuracy experiment + +def multirate_accuracy_experiment(errors_outf): + results_by_substep_ratio = {} + eocs = [] + + substep_ratios = [(1, 1, 1), (1, 2, 2), (1, 2, 4)] + dts_orig = 2 ** np.arange(-11, -17., -1) + problem = VariableCoeffWaveEquationProblem(100, (16, 4, 1)) + component_names = ("comp0", "comp1", "comp2") + + initial_context = { + name: problem.u0_by_component[i] + for i, name in enumerate(component_names)} + + for substep_ratio in substep_ratios: + dts = dts_orig * max(substep_ratio) + + for dt in dts: + logger.info("running ratio '%s' with dt=%e", substep_ratio, dt) + stepper = make_3_component_multirate_method( + problem, + ratios=substep_ratio) + stepper.set_up(t_start=0, dt_start=dt, context=initial_context) + + state_components = {} + + for event in stepper.run(t_end=0.5): + if isinstance(event, stepper.StateComputed): + state_components[event.component_id] = \ + event.state_component + + results_by_substep_ratio.setdefault(substep_ratio, []).append( + problem.full_solution(list( + state_components[component] + for component in component_names))) + + rows = [] + rows.append( + [r"$\Delta t_\text{fast}$"] + + [f"{substep}" for substep in substep_ratios]) + + from pytools.convergence import EOCRecorder + eocs = {s: EOCRecorder() for s in substep_ratios} + + for i_dt, dt in enumerate(dts_orig[:-1]): + row = [f"{dt:.2e}"] + for substep_ratio in substep_ratios: + ref_result = results_by_substep_ratio[substep_ratio][-1] + result = results_by_substep_ratio[substep_ratio][i_dt] + error = la.norm(result - ref_result, ord=np.inf) + row.append(f"{error:.2e}") + eocs[substep_ratio].add_data_point(dt, error) + rows.append(row) + + if PAPER_OUTPUT: + rows.append([r"\midrow"]) + row = [r"\multicolumn{1}{l}{Order}"] + for s in substep_ratios: + row.append(f"{eocs[s].order_estimate():.2f}") + rows.append(row) + else: + row = ["Order"] + for s in substep_ratios: + row.append(f"{eocs[s].order_estimate():.2f}") + rows.append(row) + + print( + tabulate(rows, col_fmt="S" * (1 + len(substep_ratios))), + file=errors_outf) + +# }}} + + +# {{{ table generation + +def tabulate_latex(rows, col_fmt): + result = [] + result.append(f"\\begin{{tabular}}{{{col_fmt}}}") + result.append(r"\toprule") + result.append( + " & ".join((r"\multicolumn{1}{c}{%s}" % t) for t in rows[0]) + r"\\") + result.append(r"\midrule") + for row in rows[1:]: + result.append(" & ".join(row) + r"\\") + result.append(r"\bottomrule") + result.append(r"\end{tabular}") + return "\n".join(result) + + +def tabulate_ascii(rows, col_fmt): + del col_fmt + from pytools import Table + result = Table() + for row in rows: + result.add_row(row) + return str(result) + +# }}} + + +if PAPER_OUTPUT: + tabulate = tabulate_latex +else: + tabulate = tabulate_ascii + + +@contextmanager +def open_or_stdout(filename): + if not PAPER_OUTPUT: + import sys + yield sys.stdout + else: + with open(filename, "w") as outf: + logger.info("opened '%s' for writing", filename) + yield outf + + +def main(): + if 0: + # Plot an example solution. + plot_example(1000) + + elif 0: + # Generate stability results. + # Generating the step matrices takes a long time. + + step_ratios = ( + (1, 1, 1), + (1, 2, 2), + (1, 2, 4), + (1, 2, 6), + ) + + filenames = [] + + ngridpoints = 100 + + for step_ratio in step_ratios: + filename = "mat%d-%d-%d-%d.pkl" % ((ngridpoints,) + step_ratio) + + import os + if not os.path.exists(filename): + filename = generate_mrab_step_matrix( + ngridpoints, (16, 4, 1), step_ratio) + else: + logger.info("using saved step matrix '%s'" % filename) + + filenames.append(filename) + + with open_or_stdout("mrab-stable-dts.tex") as stable_dts_outf: + compute_all_stable_timesteps(filenames, stable_dts_outf) + + elif 1: + # Accuracy results + with open_or_stdout("mrab-errors.tex") as errors_outf: + multirate_accuracy_experiment(errors_outf) + + +if __name__ == '__main__': + main() + +# vim: foldmethod=marker -- GitLab From 0aa07cc0917ba3c8e454585b3ebcaf3cf9b33b68 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sat, 4 May 2019 00:46:27 -0500 Subject: [PATCH 2/5] Print out code --- examples/semi-implicit-rk/semi-implicit-rk.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/semi-implicit-rk/semi-implicit-rk.py b/examples/semi-implicit-rk/semi-implicit-rk.py index e778227..cfc8374 100755 --- a/examples/semi-implicit-rk/semi-implicit-rk.py +++ b/examples/semi-implicit-rk/semi-implicit-rk.py @@ -110,6 +110,7 @@ def demo_rk_implicit(): from leap.implicit import replace_AssignImplicit code = replace_AssignImplicit(code, solver_hook) + print(code) from dagrt.codegen import PythonCodeGenerator IRKMethod = PythonCodeGenerator("IRKMethod").get_class(code) -- GitLab From 0b596997bce2315a2325f72084d3655c38fcc86e Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sat, 4 May 2019 07:55:10 +0200 Subject: [PATCH 3/5] Require python3 --- examples/variable-coeff-wave-equation/wave-equation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/variable-coeff-wave-equation/wave-equation.py b/examples/variable-coeff-wave-equation/wave-equation.py index b4880c6..f8f2040 100755 --- a/examples/variable-coeff-wave-equation/wave-equation.py +++ b/examples/variable-coeff-wave-equation/wave-equation.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 """Solves the 1D wave equation u_tt = c(x) u_xx -- GitLab From 7b3656baf4d7f3c9ec3c054421482883cd4dcb4a Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 5 May 2019 00:21:06 -0500 Subject: [PATCH 4/5] Update citation --- examples/semi-implicit-rk/semi-implicit-rk.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/semi-implicit-rk/semi-implicit-rk.py b/examples/semi-implicit-rk/semi-implicit-rk.py index cfc8374..81b3a77 100755 --- a/examples/semi-implicit-rk/semi-implicit-rk.py +++ b/examples/semi-implicit-rk/semi-implicit-rk.py @@ -3,9 +3,9 @@ Source: - Kennedy, Christopher A., and Mark H. Carpenter. - "Diagonally implicit runge-kutta - methods for ordinary differential equations. a review." (2016). + Diagonally Implicit Runge–Kutta Methods for Stiff O.D.E.’s + Roger Alexander + SIAM Journal on Numerical Analysis 1977 14:6, 1006-1021 """ -- GitLab From 0f8f92307230830b114b1e9f1c913bff1e1a3fe3 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 6 May 2019 23:58:55 -0500 Subject: [PATCH 5/5] Improve citations and other documentation --- examples/adaptive-rk/adaptive-rk.py | 16 +++++++++---- examples/semi-implicit-rk/semi-implicit-rk.py | 7 +++--- .../wave-equation.py | 23 +++++++++++++++---- 3 files changed, 33 insertions(+), 13 deletions(-) diff --git a/examples/adaptive-rk/adaptive-rk.py b/examples/adaptive-rk/adaptive-rk.py index 792b108..3cba6ca 100755 --- a/examples/adaptive-rk/adaptive-rk.py +++ b/examples/adaptive-rk/adaptive-rk.py @@ -4,10 +4,10 @@ an order 2/3 method pair. Source: - Bogacki, Przemyslaw; Shampine, Lawrence F. (1989), - "A 3(2) pair of Runge–Kutta formulas", - Applied Mathematics Letters, 2 (4): 321–325, - doi:10.1016/0893-9659(89)90079-7. + P. Bogacki and L. F. Shampine. + A 3(2) pair of Runge-Kutta formulas. + Appl. Math. Lett., 2(4):321–325, 1989. + doi: 10.1016/0893-9659(89)90079-7. The same functionality is supported using leap.rk.ODE23Method. @@ -22,7 +22,13 @@ logging.basicConfig(level=logging.INFO) class KapsProblem(object): """ - From Kennedy and Carpenter, Section 7.1 + Source: Section 7.1 of + + Christopher A. Kennedy and Mark H. Carpenter. + Additive Runge-Kutta schemes for convection-diffusion-reaction + equations. + Appl. Numer. Math., 44 (1-2):139–181, 2003. + doi: 10.1016/S0168-9274(02)00138-1. y_1' = - (epsilon^{-1} + 2) y_1 + epsilon^{-1} y_2^2 y_2' = y_1 - y_2 - y_2^2 diff --git a/examples/semi-implicit-rk/semi-implicit-rk.py b/examples/semi-implicit-rk/semi-implicit-rk.py index 81b3a77..2988af9 100755 --- a/examples/semi-implicit-rk/semi-implicit-rk.py +++ b/examples/semi-implicit-rk/semi-implicit-rk.py @@ -3,9 +3,10 @@ Source: - Diagonally Implicit Runge–Kutta Methods for Stiff O.D.E.’s - Roger Alexander - SIAM Journal on Numerical Analysis 1977 14:6, 1006-1021 + Roger Alexander. + Diagonally implicit Runge–Kutta methods for stiff O.D.E.'s. + SIAM Journal on Numerical Analysis, 14(6):1006–1021, 1977. + doi: 10.1137/0714068. """ diff --git a/examples/variable-coeff-wave-equation/wave-equation.py b/examples/variable-coeff-wave-equation/wave-equation.py index f8f2040..547c7c5 100755 --- a/examples/variable-coeff-wave-equation/wave-equation.py +++ b/examples/variable-coeff-wave-equation/wave-equation.py @@ -1,11 +1,12 @@ #!/usr/bin/env python3 """Solves the 1D wave equation - u_tt = c(x) u_xx + u_tt - c(x)^2 u_xx = 0 u_t(0) = init' u(0) = init with piecewise constant coefficients c(x) using a multirate multistep method. + """ @@ -64,15 +65,27 @@ class VariableCoeffWaveEquationProblem(object): This turns the problem instance into the 1D system u_t = v - v_t = c(x) * u_xx. - - The discretization uses a centered difference stencil. Interface conditions - are imposed using a second-order biased stencil. + v_t = c(x)^2 * u_xx. The system state is separated into components, i.e. regions in which c(x) is constant. Components are equispaced in the interval (0, 1) and correspond to the number of coefficients passed to the constructor. + The discretization uses a centered difference stencil. + + The interface conditions to the continuous problem require continuity of u + and c(x)^2 * u_x at the interfaces. These conditions are imposed on the + discretization using a second-order biased stencil by ensuring that the + values at the interface implied by the stencil, when applied from both the + left and the right sides, are the same. For more on interface conditions, + see + + David L. Brown. + A note on the numerical solution of the wave equation with piecewise + smooth coefficients. + Math. Comp., 42(166):369–391, 1984. + doi: 10.2307/2007591. + """ def __init__(self, ngridpoints, component_coeffs, decay=0, diffusion=0): -- GitLab