Verified order reduction

Author: brownbaerchen

pySDC/playgrounds/time_dep_BCs/run_time_dep_heat.py

@@ -8,7 +8,7 @@
 from pySDC.implementations.hooks.plotting import PlotPostStep
 
 
-def run_heat(dt=1e-1, Tend=4, plotting=False):
+def run_heat(dt=1e-1, Tend=4, kmax=5, ft=np.pi, plotting=False):
     level_params = {}
     level_params['dt'] = dt
 
@@ -17,13 +17,13 @@ def run_heat(dt=1e-1, Tend=4, plotting=False):
     sweeper_params['num_nodes'] = 3
     sweeper_params['QI'] = 'LU'
 
-    problem_params = {}
+    problem_params = {'ft': ft, 'spectral_space': False}
 
     step_params = {}
-    step_params['maxiter'] = 4
+    step_params['maxiter'] = kmax
 
     controller_params = {}
-    controller_params['logger_level'] = 15
+    controller_params['logger_level'] = 30
     if plotting:
         controller_params['hook_class'] = PlotPostStep
 
@@ -44,9 +44,57 @@ def run_heat(dt=1e-1, Tend=4, plotting=False):
 
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
-    return stats, controller
+    return uend, stats, controller
+
+
+def compute_errors(Tend, N, ft, kmax):
+    solutions = []
+    dts = []
+
+    for n in range(N):
+        dt = Tend / (2**n)
+        u, _, _ = run_heat(dt=dt, Tend=Tend, ft=ft, kmax=kmax, plotting=False)
+
+        solutions.append(u)
+        dts.append(dt)
+
+    errors = [abs(solutions[i] - solutions[-1]) / abs(solutions[-1]) for i in range(N - 1)]
+    _dts = [dts[i] for i in range(N - 1)]
+
+    return _dts, errors
+
+
+def compute_order(dts, errors):
+    orders = []
+    for i in range(len(errors) - 1):
+        if errors[i + 1] < 1e-12:
+            break
+        orders.append(np.log(errors[i] / errors[i + 1]) / np.log(dts[i] / dts[i + 1]))
+    return np.median(orders)
+
+
+def plot_order():  # pragma: no cover
+    fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(9, 4))
+
+    for k in range(1, 10):
+        dts, errors = compute_errors(4, 8, 0, k)
+        order = compute_order(dts, errors)
+        axs[0].loglog(dts, errors, label=f'{k=}, {order=:.1f}')
+
+        dts, errors = compute_errors(4, 8, 1, k)
+        order = compute_order(dts, errors)
+        axs[1].loglog(dts, errors, label=f'{k=}, {order=:.1f}')
+
+    axs[0].set_ylabel(r'relative global error')
+    axs[0].set_xlabel(r'$\Delta t$')
+    axs[1].set_xlabel(r'$\Delta t$')
+    axs[0].legend(frameon=False)
+    axs[1].legend(frameon=False)
+    axs[0].set_title('Constant in time BCs')
+    axs[1].set_title('Time-dependent BCs')
+    fig.tight_layout()
 
 
 if __name__ == '__main__':
-    run_heat(plotting=True)
+    plot_order()
     plt.show()

pySDC/playgrounds/time_dep_BCs/tests.py

@@ -19,7 +19,7 @@ def test_time_dep_heat_eq(a, b, t):
     expect_boundary = np.empty((2, 2))
     expect_boundary = problem.put_time_dep_BCs_in_rhs(expect_boundary, t)
 
-    # we use T_n(1) = 1 and T_n(-1) = (1)^n, to compute the values at the boundaries from the spectral representation
+    # we use T_n(1) = 1 and T_n(-1) = (-1)^n, to compute the values at the boundaries from the spectral representation
     # see Wikipedia for more details: https://en.wikipedia.org/wiki/Chebyshev_polynomials#Roots_and_extrema
     right_boundary = u.sum()
     expect_right_boundary = expect_boundary[0, -2]