diff --git a/.gitignore b/.gitignore index bbaaa01..36588df 100644 --- a/.gitignore +++ b/.gitignore @@ -30,3 +30,9 @@ *.vrb *.xdy *.tdo +*.loa + +*.bak +*.brf +*.pre +LWB/propagators/*.pdf diff --git a/LWB/propagators/Makefile b/LWB/propagators/Makefile new file mode 100644 index 0000000..1d18ca7 --- /dev/null +++ b/LWB/propagators/Makefile @@ -0,0 +1,22 @@ +PDFLATEX = pdflatex +BIBTEX = bibtex +SECTIONS = $(wildcard sections/*.tex) + +all: main.pdf + +.PHONY: main.pdf +main.pdf: main.tex header.tex $(SECTIONS) uml.png # figs + $(PDFLATEX) main.tex + $(BIBTEX) main.aux + $(PDFLATEX) main.tex + $(PDFLATEX) main.tex + +uml.png: uml.dot + dot -Tpng $< -o $@ + +# figs: +# $(MAKE) -C figures + +clean: + $(MAKE) clean -C figures + rm -f *.aux *.bbl *.blg *.log *.out *.toc diff --git a/LWB/propagators/declaration-orig.pdf b/LWB/propagators/declaration-orig.pdf new file mode 100644 index 0000000..4bec829 Binary files /dev/null and b/LWB/propagators/declaration-orig.pdf differ diff --git a/LWB/propagators/figures/coefficient_analysis.png b/LWB/propagators/figures/coefficient_analysis.png new file mode 100644 index 0000000..ef1fc5b Binary files /dev/null and b/LWB/propagators/figures/coefficient_analysis.png differ diff --git a/LWB/propagators/figures/coefs.png b/LWB/propagators/figures/coefs.png new file mode 100644 index 0000000..32eaf0d Binary files /dev/null and b/LWB/propagators/figures/coefs.png differ diff --git a/LWB/propagators/figures/convergence_harmonic_1D_lt.png b/LWB/propagators/figures/convergence_harmonic_1D_lt.png new file mode 100644 index 0000000..7ae4234 Binary files /dev/null and b/LWB/propagators/figures/convergence_harmonic_1D_lt.png differ diff --git a/LWB/propagators/figures/convergence_harmonic_1D_y4.png b/LWB/propagators/figures/convergence_harmonic_1D_y4.png new file mode 100644 index 0000000..ac6f643 Binary files /dev/null and b/LWB/propagators/figures/convergence_harmonic_1D_y4.png differ diff --git a/LWB/propagators/figures/convergence_legend.png b/LWB/propagators/figures/convergence_legend.png new file mode 100644 index 0000000..3861e0b Binary files /dev/null and b/LWB/propagators/figures/convergence_legend.png differ diff --git a/LWB/propagators/figures/convergence_morse_1D_lt.png b/LWB/propagators/figures/convergence_morse_1D_lt.png new file mode 100644 index 0000000..7401c24 Binary files /dev/null and b/LWB/propagators/figures/convergence_morse_1D_lt.png differ diff --git a/LWB/propagators/figures/convergence_morse_1D_y4.png b/LWB/propagators/figures/convergence_morse_1D_y4.png new file mode 100644 index 0000000..0da4544 Binary files /dev/null and b/LWB/propagators/figures/convergence_morse_1D_y4.png differ diff --git a/LWB/propagators/figures/convergence_torsional_1D_lt.png b/LWB/propagators/figures/convergence_torsional_1D_lt.png new file mode 100644 index 0000000..8ad2741 Binary files /dev/null and b/LWB/propagators/figures/convergence_torsional_1D_lt.png differ diff --git a/LWB/propagators/figures/convergence_torsional_1D_y4.png b/LWB/propagators/figures/convergence_torsional_1D_y4.png new file mode 100644 index 0000000..e927a1c Binary files /dev/null and b/LWB/propagators/figures/convergence_torsional_1D_y4.png differ diff --git a/LWB/propagators/figures/dim.png b/LWB/propagators/figures/dim.png new file mode 100644 index 0000000..ca37dc9 Binary files /dev/null and b/LWB/propagators/figures/dim.png differ diff --git a/LWB/propagators/figures/dimension_analysis.png b/LWB/propagators/figures/dimension_analysis.png new file mode 100644 index 0000000..ca37dc9 Binary files /dev/null and b/LWB/propagators/figures/dimension_analysis.png differ diff --git a/LWB/propagators/figures/error_analysis_lt.png b/LWB/propagators/figures/error_analysis_lt.png new file mode 100644 index 0000000..015196f Binary files /dev/null and b/LWB/propagators/figures/error_analysis_lt.png differ diff --git a/LWB/propagators/figures/error_analysis_y4.png b/LWB/propagators/figures/error_analysis_y4.png new file mode 100644 index 0000000..fbde3ee Binary files /dev/null and b/LWB/propagators/figures/error_analysis_y4.png differ diff --git a/LWB/propagators/figures/eth-logo.pdf b/LWB/propagators/figures/eth-logo.pdf new file mode 100644 index 0000000..34d6635 Binary files /dev/null and b/LWB/propagators/figures/eth-logo.pdf differ diff --git a/LWB/propagators/figures/harmonic_1D_Hagedorn_drift.png b/LWB/propagators/figures/harmonic_1D_Hagedorn_drift.png new file mode 100644 index 0000000..37f96be Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_Hagedorn_drift.png differ diff --git a/LWB/propagators/figures/harmonic_1D_Hagedorn_energies.png b/LWB/propagators/figures/harmonic_1D_Hagedorn_energies.png new file mode 100644 index 0000000..1c118a3 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_Hagedorn_energies.png differ diff --git a/LWB/propagators/figures/harmonic_1D_MG4_drift.png b/LWB/propagators/figures/harmonic_1D_MG4_drift.png new file mode 100644 index 0000000..16033cb Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_MG4_drift.png differ diff --git a/LWB/propagators/figures/harmonic_1D_MG4_energies.png b/LWB/propagators/figures/harmonic_1D_MG4_energies.png new file mode 100644 index 0000000..9ccf9d0 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_MG4_energies.png differ diff --git a/LWB/propagators/figures/harmonic_1D_McL42_drift.png b/LWB/propagators/figures/harmonic_1D_McL42_drift.png new file mode 100644 index 0000000..64c01ce Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_McL42_drift.png differ diff --git a/LWB/propagators/figures/harmonic_1D_McL42_energies.png b/LWB/propagators/figures/harmonic_1D_McL42_energies.png new file mode 100644 index 0000000..9ccf9d0 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_McL42_energies.png differ diff --git a/LWB/propagators/figures/harmonic_1D_McL84_drift.png b/LWB/propagators/figures/harmonic_1D_McL84_drift.png new file mode 100644 index 0000000..423dba3 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_McL84_drift.png differ diff --git a/LWB/propagators/figures/harmonic_1D_McL84_energies.png b/LWB/propagators/figures/harmonic_1D_McL84_energies.png new file mode 100644 index 0000000..9ccf9d0 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_McL84_energies.png differ diff --git a/LWB/propagators/figures/harmonic_1D_Pre764_drift.png b/LWB/propagators/figures/harmonic_1D_Pre764_drift.png new file mode 100644 index 0000000..1e2e937 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_Pre764_drift.png differ diff --git a/LWB/propagators/figures/harmonic_1D_Pre764_energies.png b/LWB/propagators/figures/harmonic_1D_Pre764_energies.png new file mode 100644 index 0000000..6529f02 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_Pre764_energies.png differ diff --git a/LWB/propagators/figures/harmonic_1D_Semiclassical_drift.png b/LWB/propagators/figures/harmonic_1D_Semiclassical_drift.png new file mode 100644 index 0000000..b8a8679 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_Semiclassical_drift.png differ diff --git a/LWB/propagators/figures/harmonic_1D_Semiclassical_energies.png b/LWB/propagators/figures/harmonic_1D_Semiclassical_energies.png new file mode 100644 index 0000000..9ccf9d0 Binary files /dev/null and b/LWB/propagators/figures/harmonic_1D_Semiclassical_energies.png differ diff --git a/LWB/propagators/figures/makeplots.py b/LWB/propagators/figures/makeplots.py new file mode 100644 index 0000000..2e5e166 --- /dev/null +++ b/LWB/propagators/figures/makeplots.py @@ -0,0 +1,278 @@ +#!/usr/bin/env python + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.ticker import MaxNLocator + +####################################### +# Plot 1: Dimension analysis +####################################### + +D = [ 2,3,4,5 ] +time_D = [ 0.08,5.12,24.06,778.75 ] +time_D = [ 0.0448, 0.3493, 24.0235, 2105.5222 ] + +fig = plt.figure() +ax = fig.gca() + +ax.semilogy(D, time_D, 'k-o') + +# plt.margins(x=0.1, y=0.1) +ax.set_xlim([1.5,5.5]) +ax.set_ylim([1e-2,1e4]) +ax.xaxis.set_major_locator(MaxNLocator(integer=True)) +# ax.set_xlim([0,1]) +# ax.set_ylim(view[2:]) +# ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="y") +ax.set_title("Compute Time vs. Dimension of Wave Packet") +ax.grid(True) +ax.set_xlabel(r"Dimension $D$") +ax.set_ylabel(r"Compute time $t$ [s]") +fig.savefig("dimension_analysis.png") +plt.show() +plt.close(fig) + + + +####################################### +# Plot 2: Coefficient Analysis +####################################### + +size_coef = np.array([ 1,2,4,7,8,10,15,18,34 ]) +time_python = np.array([47.155, 52.934, 65.351, 84.075, 90.260, 102.531, 133.415, 151.438, 249.749]) +time_cpp = np.array([1.3905, 1.3817, 1.4386, 1.5122, 1.4568, 1.6049, 1.7763, 1.8147, 2.2168]) +speedup = time_python / time_cpp + +print(speedup) + +fig, ax = plt.subplots(2,1) + +ax[0].loglog(size_coef, time_python, 'b-o', label=r"Python") +ax[0].loglog(size_coef, time_cpp, 'g-o', label=r"C++") +ax[1].loglog(size_coef, speedup, 'r-o') + +# ax.set_xlim([0,1]) +# ax.set_ylim(view[2:]) +# ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="y") +ax[0].set_xlim([.8,40]) +ax[1].set_xlim([.8,40]) +ax[1].set_ylim([25,160]) +ax[1].set_yscale('log', basey=2) +ax[0].set_title("Cost of IntSplit splitting coefficients") +ax[0].legend(loc='upper left') +ax[0].grid(True) +ax[1].grid(True) +ax[0].set_ylabel(r"Compute Time $t$ [s]") +ax[1].set_ylabel(r"Speedup") +ax[1].set_xlabel(r"Size of coefficient array") +fig.savefig("coefficient_analysis.png") +plt.show() +plt.close(fig) + + + +####################################### +# Plot 3: Error analysis +####################################### + +n_steps = [ 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.007812, 0.003906, 0.001953, 0.000977, 0.000488, 0.000244 ] +n_steps_morse = [ 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.007812, 0.003906, 0.001953 ] + +# harmonic LT +error_harmonic_1D_lt_hagedorn = [ 9.545477e-07, 3.818193e-06, 1.527278e-05, 6.109122e-05, 2.443662e-04, 9.774862e-04, 3.910282e-03, 1.564622e-02, 6.264607e-02, 2.502071e-01, 9.070318e-01] +error_harmonic_1D_lt_mg4 = [ 3.968109e-12, 3.105810e-12, 2.779794e-12, 8.872414e-13, 3.629015e-12, 2.571130e-11, 3.863204e-10, 6.157092e-09, 9.848602e-08, 1.575815e-06, 2.521706e-05] +error_harmonic_1D_lt_mcl42 = [ 4.170521e-12, 2.910990e-12, 2.896429e-12, 7.860211e-13, 3.391014e-12, 2.569097e-11, 3.861969e-10, 6.156957e-09, 9.848581e-08, 1.575815e-06, 2.521706e-05] +error_harmonic_1D_lt_mcl84 = [ 4.424280e-12, 7.103594e-13, 4.890004e-12, 3.017013e-13, 4.049679e-12, 2.928067e-11, 4.425939e-10, 7.059074e-09, 1.129233e-07, 1.806804e-06, 2.891260e-05] +error_harmonic_1D_lt_pre764 = [ 3.154724e-12, 1.388741e-12, 2.953003e-12, 1.909897e-12, 4.346275e-12, 4.220303e-11, 6.432055e-10, 1.007762e-08, 1.554063e-07, 2.301564e-06, 3.102992e-05] +error_harmonic_1D_lt_semiclassical = [ 6.570944e-13, 2.810143e-12, 8.654158e-13, 2.010771e-12, 3.424029e-12, 2.499841e-11, 3.718457e-10, 5.920513e-09, 9.470839e-08, 1.515352e-06, 2.424855e-05] + +# torsional LT +error_torsional_1D_lt_hagedorn = [ 4.413159e-08, 1.766168e-07, 7.065577e-07, 2.826307e-06, 1.130505e-05, 4.521597e-05, 1.807949e-04, 7.220749e-04, 2.870631e-03, 1.120047e-02, 4.033724e-02] +error_torsional_1D_lt_mg4 = [ 3.389891e-11, 3.717412e-11, 3.493770e-11, 3.225172e-11, 3.261378e-11, 3.397044e-11, 5.169600e-11, 3.423753e-10, 5.017209e-09, 7.981489e-08, 1.276773e-06] +error_torsional_1D_lt_mcl42 = [ 3.379270e-11, 3.702897e-11, 3.479613e-11, 3.327578e-11, 4.706530e-11, 1.434777e-10, 5.594621e-10, 2.195010e-09, 9.091761e-09, 7.849469e-08, 1.247707e-06] +error_torsional_1D_lt_mcl84 = [ 2.446365e-11, 3.216573e-11, 3.219421e-11, 3.325376e-11, 3.309425e-11, 3.384717e-11, 5.432106e-11, 3.926380e-10, 5.815999e-09, 9.259400e-08, 1.481307e-06] +error_torsional_1D_lt_pre764 = [ 3.414883e-11, 3.064713e-11, 3.417420e-11, 3.024380e-11, 3.211200e-11, 3.464738e-11, 6.443673e-11, 5.365296e-10, 7.648273e-09, 1.077038e-07, 1.305210e-06] +error_torsional_1D_lt_semiclassical = [ 3.279210e-11, 3.760648e-11, 3.617805e-11, 6.278279e-11, 1.945545e-10, 7.432790e-10, 2.947945e-09, 1.180447e-08, 4.801044e-08, 2.158989e-07, 1.608904e-06] + +# morse LT +error_morse_1D_lt_hagedorn = [ 1.325977e-08, 5.307810e-08, 2.123536e-07, 8.494546e-07, 3.397868e-06, 1.359167e-05, 5.436919e-05, 2.175166e-04] +error_morse_1D_lt_mg4 = [ 9.032965e-06, 1.806527e-05, 3.612797e-05, 7.224567e-05, 1.444503e-04, 2.887368e-04, 5.768199e-04, 1.151037e-03] +error_morse_1D_lt_mcl42 = [ 9.032949e-06, 1.806521e-05, 3.612772e-05, 7.224467e-05, 1.444463e-04, 2.887210e-04, 5.767573e-04, 1.150793e-03] +error_morse_1D_lt_mcl84 = [ 1.208038e-05, 2.416038e-05, 4.831927e-05, 9.663264e-05, 1.932418e-04, 3.863905e-04, 7.724157e-04, 1.543424e-03] +error_morse_1D_lt_pre764 = [ 1.357298e-05, 2.716454e-05, 5.440344e-05, 1.091042e-04, 2.193969e-04, 4.435388e-04, 9.059879e-04, 1.887033e-03] +error_morse_1D_lt_semiclassical = [ 1.068625e-05, 2.137201e-05, 4.274209e-05, 8.547655e-05, 1.709227e-04, 3.417245e-04, 6.829718e-04, 1.364083e-03] + +######################################### + +# harmonic Y4 +error_harmonic_1D_y4_hagedorn = [ 9.545477e-07, 3.818193e-06, 1.527278e-05, 6.109122e-05, 2.443662e-04, 9.774862e-04, 3.910282e-03, 1.564622e-02, 6.264607e-02, 2.502071e-01, 9.070318e-01] +error_harmonic_1D_y4_mg4 = [ 3.968109e-12, 3.105810e-12, 2.779794e-12, 8.872414e-13, 3.629015e-12, 2.571130e-11, 3.863204e-10, 6.157092e-09, 9.848602e-08, 1.575815e-06, 2.521706e-05] +error_harmonic_1D_y4_mcl42 = [ 4.170521e-12, 2.910990e-12, 2.896429e-12, 7.860211e-13, 3.391014e-12, 2.569097e-11, 3.861969e-10, 6.156957e-09, 9.848581e-08, 1.575815e-06, 2.521706e-05] +error_harmonic_1D_y4_mcl84 = [ 4.424280e-12, 7.103594e-13, 4.890004e-12, 3.017013e-13, 4.049679e-12, 2.928067e-11, 4.425939e-10, 7.059074e-09, 1.129233e-07, 1.806804e-06, 2.891260e-05] +error_harmonic_1D_y4_pre764 = [ 3.154724e-12, 1.388741e-12, 2.953003e-12, 1.909897e-12, 4.346275e-12, 4.220303e-11, 6.432055e-10, 1.007762e-08, 1.554063e-07, 2.301564e-06, 3.102992e-05] +error_harmonic_1D_y4_semiclassical = [ 6.570944e-13, 2.810143e-12, 8.654158e-13, 2.010771e-12, 3.424029e-12, 2.499841e-11, 3.718457e-10, 5.920513e-09, 9.470839e-08, 1.515352e-06, 2.424855e-05] + +# torsional Y4 +error_torsional_1D_y4_hagedorn = [ 4.413159e-08, 1.766168e-07, 7.065577e-07, 2.826307e-06, 1.130505e-05, 4.521597e-05, 1.807949e-04, 7.220749e-04, 2.870631e-03, 1.120047e-02, 4.033724e-02] +error_torsional_1D_y4_mg4 = [ 3.389891e-11, 3.717412e-11, 3.493770e-11, 3.225172e-11, 3.261378e-11, 3.397044e-11, 5.169600e-11, 3.423753e-10, 5.017209e-09, 7.981489e-08, 1.276773e-06] +error_torsional_1D_y4_mcl42 = [ 3.379270e-11, 3.702897e-11, 3.479613e-11, 3.327578e-11, 4.706530e-11, 1.434777e-10, 5.594621e-10, 2.195010e-09, 9.091761e-09, 7.849469e-08, 1.247707e-06] +error_torsional_1D_y4_mcl84 = [ 2.446365e-11, 3.216573e-11, 3.219421e-11, 3.325376e-11, 3.309425e-11, 3.384717e-11, 5.432106e-11, 3.926380e-10, 5.815999e-09, 9.259400e-08, 1.481307e-06] +error_torsional_1D_y4_pre764 = [ 3.414883e-11, 3.064713e-11, 3.417420e-11, 3.024380e-11, 3.211200e-11, 3.464738e-11, 6.443673e-11, 5.365296e-10, 7.648273e-09, 1.077038e-07, 1.305210e-06] +error_torsional_1D_y4_semiclassical = [ 3.279210e-11, 3.760648e-11, 3.617805e-11, 6.278279e-11, 1.945545e-10, 7.432790e-10, 2.947945e-09, 1.180447e-08, 4.801044e-08, 2.158989e-07, 1.608904e-06] + +# morse Y4 +error_morse_1D_y4_hagedorn = [ 1.325977e-08, 5.307810e-08, 2.123536e-07, 8.494546e-07, 3.397868e-06, 1.359167e-05, 5.436919e-05, 2.175166e-04] +error_morse_1D_y4_mg4 = [ 1.460253e-11, 1.472816e-11, 1.457033e-11, 1.520796e-11, 1.527166e-11, 3.794749e-11, 5.579321e-10, 8.913417e-09] +error_morse_1D_y4_mcl42 = [ 2.962581e-10, 1.182166e-09, 4.726521e-09, 1.890463e-08, 7.561734e-08, 3.024688e-07, 1.209888e-06, 4.839784e-06] +error_morse_1D_y4_mcl84 = [ 1.517644e-11, 1.480756e-11, 1.431413e-11, 1.569338e-11, 1.420991e-11, 2.342160e-11, 2.857094e-10, 4.524417e-09] +error_morse_1D_y4_pre764 = [ 1.579363e-11, 1.498456e-11, 1.553254e-11, 1.435766e-11, 1.516359e-11, 2.999104e-11, 4.127790e-10, 6.227530e-09] +error_morse_1D_y4_semiclassical = [ 1.125289e-09, 4.496678e-09, 1.798135e-08, 7.192046e-08, 2.876774e-07, 1.150712e-06, 4.602971e-06, 1.841389e-05] + +######################################### + +# torsional 2D + + +fig = plt.figure() +ax = fig.gca() + +ax.loglog(n_steps, error_harmonic_1D_lt_hagedorn, 'r-o', label=r"Hagedorn") +ax.loglog(n_steps, error_harmonic_1D_lt_semiclassical, 'm-o', label=r"Semiclassical") +ax.loglog(n_steps, error_harmonic_1D_lt_mg4, 'b-o', label=r"Magnus (4)") +ax.loglog(n_steps, error_harmonic_1D_lt_mcl42, 'g-o', label=r"McLachlan (4,2)") +ax.loglog(n_steps, error_harmonic_1D_lt_mcl84, 'k-o', label=r"McLachlan (8,4)") +ax.loglog(n_steps, error_harmonic_1D_lt_pre764, 'c-o', label=r"Processing (7,6,4)") + +# lgd = ax.legend(bbox_to_anchor=(1.02, 1), loc=2, borderaxespad=0.) +ax.set_xlim([5e-1,1e-4]) +# ax.set_ylim(view[2:]) +# ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="y") +ax.set_title(r"Convergence for Harmonic 1D potential (LT splitting)") +ax.grid(True) +ax.set_xlabel(r"Size of time step $\Delta t$") +ax.set_ylabel(r"$L_2$ error $\frac{\| u (x) - u_* (x) \|}{\| u_* (x) \|}$") +# fig.savefig("convergence_harmonic_1D_lt.png",bbox_extra_artists=(lgd,), bbox_inches='tight') +fig.savefig("convergence_harmonic_1D_lt.png") +plt.show() +plt.close(fig) + + +fig = plt.figure() +ax = fig.gca() + +ax.loglog(n_steps, error_torsional_1D_lt_hagedorn, 'r-o', label=r"Hagedorn") +ax.loglog(n_steps, error_torsional_1D_lt_semiclassical, 'm-o', label=r"Semiclassical") +ax.loglog(n_steps, error_torsional_1D_lt_mg4, 'b-o', label=r"Magnus (4)") +ax.loglog(n_steps, error_torsional_1D_lt_mcl42, 'g-o', label=r"McLachlan (4,2)") +ax.loglog(n_steps, error_torsional_1D_lt_mcl84, 'k-o', label=r"McLachlan (8,4)") +ax.loglog(n_steps, error_torsional_1D_lt_pre764, 'c-o', label=r"Processing (7,6,4)") + +# lgd = ax.legend(bbox_to_anchor=(1.02, 1), loc=2, borderaxespad=0.) +ax.set_xlim([5e-1,1e-4]) +# ax.set_ylim(view[2:]) +# ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="y") +ax.set_title(r"Convergence for Torsional 1D potential (LT splitting)") +ax.grid(True) +ax.set_xlabel(r"Size of time step $\Delta t$") +ax.set_ylabel(r"$L_2$ error $\frac{\| u (x) - u_* (x) \|}{\| u_* (x) \|}$") +# fig.savefig("convergence_torsional_1D_lt.png",bbox_extra_artists=(lgd,), bbox_inches='tight') +fig.savefig("convergence_torsional_1D_lt.png") +plt.show() +plt.close(fig) + + +fig = plt.figure() +ax = fig.gca() + +ax.loglog(n_steps_morse, error_morse_1D_lt_hagedorn, 'r-o', label=r"Hagedorn") +ax.loglog(n_steps_morse, error_morse_1D_lt_semiclassical, 'm-o', label=r"Semiclassical") +ax.loglog(n_steps_morse, error_morse_1D_lt_mg4, 'b-o', label=r"Magnus (4)") +ax.loglog(n_steps_morse, error_morse_1D_lt_mcl42, 'g-o', label=r"McLachlan (4,2)") +ax.loglog(n_steps_morse, error_morse_1D_lt_mcl84, 'k-o', label=r"McLachlan (8,4)") +ax.loglog(n_steps_morse, error_morse_1D_lt_pre764, 'c-o', label=r"Processing (7,6,4)") + +# lgd = ax.legend(bbox_to_anchor=(1.02, 1), loc=2, borderaxespad=0.) +ax.set_xlim([5e-1,1e-3]) +# ax.set_ylim(view[2:]) +# ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="y") +ax.set_title(r"Convergence for Morse 1D potential (LT splitting)") +ax.grid(True) +ax.set_xlabel(r"Size of time step $\Delta t$") +ax.set_ylabel(r"$L_2$ error $\frac{\| u (x) - u_* (x) \|}{\| u_* (x) \|}$") +# fig.savefig("convergence_morse_1D_lt.png",bbox_extra_artists=(lgd,), bbox_inches='tight') +fig.savefig("convergence_morse_1D_lt.png") +plt.show() +plt.close(fig) + + + +fig = plt.figure() +ax = fig.gca() + +ax.loglog(n_steps, error_harmonic_1D_y4_hagedorn, 'r-o', label=r"Hagedorn") +ax.loglog(n_steps, error_harmonic_1D_y4_semiclassical, 'm-o', label=r"Semiclassical") +ax.loglog(n_steps, error_harmonic_1D_y4_mg4, 'b-o', label=r"Magnus (4)") +ax.loglog(n_steps, error_harmonic_1D_y4_mcl42, 'g-o', label=r"McLachlan (4,2)") +ax.loglog(n_steps, error_harmonic_1D_y4_mcl84, 'k-o', label=r"McLachlan (8,4)") +ax.loglog(n_steps, error_harmonic_1D_y4_pre764, 'c-o', label=r"Processing (7,6,4)") + +# lgd = ax.legend(bbox_to_anchor=(1.02, 1), loc=2, borderaxespad=0.) +ax.set_xlim([5e-1,1e-4]) +# ax.set_ylim(view[2:]) +# ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="y") +ax.set_title(r"Convergence for Harmonic 1D potential (Y4 splitting)") +ax.grid(True) +ax.set_xlabel(r"Size of time step $\Delta t$") +ax.set_ylabel(r"$L_2$ error $\frac{\| u (x) - u_* (x) \|}{\| u_* (x) \|}$") +# fig.savefig("convergence_harmonic_1D_y4.png",bbox_extra_artists=(lgd,), bbox_inches='tight') +fig.savefig("convergence_harmonic_1D_y4.png") +plt.show() +plt.close(fig) + + +fig = plt.figure() +ax = fig.gca() + +ax.loglog(n_steps, error_torsional_1D_y4_hagedorn, 'r-o', label=r"Hagedorn") +ax.loglog(n_steps, error_torsional_1D_y4_semiclassical, 'm-o', label=r"Semiclassical") +ax.loglog(n_steps, error_torsional_1D_y4_mg4, 'b-o', label=r"Magnus (4)") +ax.loglog(n_steps, error_torsional_1D_y4_mcl42, 'g-o', label=r"McLachlan (4,2)") +ax.loglog(n_steps, error_torsional_1D_y4_mcl84, 'k-o', label=r"McLachlan (8,4)") +ax.loglog(n_steps, error_torsional_1D_y4_pre764, 'c-o', label=r"Processing (7,6,4)") + +# lgd = ax.legend(bbox_to_anchor=(1.02, 1), loc=2, borderaxespad=0.) +ax.set_xlim([5e-1,1e-4]) +# ax.set_ylim(view[2:]) +# ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="y") +ax.set_title(r"Convergence for Torsional 1D potential (Y4 splitting)") +ax.grid(True) +ax.set_xlabel(r"Size of time step $\Delta t$") +ax.set_ylabel(r"$L_2$ error $\frac{\| u (x) - u_* (x) \|}{\| u_* (x) \|}$") +# fig.savefig("convergence_torsional_1D_y4.png",bbox_extra_artists=(lgd,), bbox_inches='tight') +fig.savefig("convergence_torsional_1D_y4.png") +plt.show() +plt.close(fig) + + +fig = plt.figure() +ax = fig.gca() + +ax.loglog(n_steps_morse, error_morse_1D_y4_hagedorn, 'r-o', label=r"Hagedorn") +ax.loglog(n_steps_morse, error_morse_1D_y4_semiclassical, 'm-o', label=r"Semiclassical") +ax.loglog(n_steps_morse, error_morse_1D_y4_mg4, 'b-o', label=r"Magnus (4)") +ax.loglog(n_steps_morse, error_morse_1D_y4_mcl42, 'g-o', label=r"McLachlan (4,2)") +ax.loglog(n_steps_morse, error_morse_1D_y4_mcl84, 'k-o', label=r"McLachlan (8,4)") +ax.loglog(n_steps_morse, error_morse_1D_y4_pre764, 'c-o', label=r"Processing (7,6,4)") + +# lgd = ax.legend(bbox_to_anchor=(1.02, 1), loc=2, borderaxespad=0.) +ax.set_xlim([5e-1,1e-3]) +# ax.set_ylim(view[2:]) +# ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="y") +ax.set_title(r"Convergence for Morse 1D potential (Y4 splitting)") +ax.grid(True) +ax.set_xlabel(r"Size of time step $\Delta t$") +ax.set_ylabel(r"$L_2$ error $\frac{\| u (x) - u_* (x) \|}{\| u_* (x) \|}$") +# fig.savefig("convergence_morse_1D_y4.png",bbox_extra_artists=(lgd,), bbox_inches='tight') +fig.savefig("convergence_morse_1D_y4.png") +plt.show() +plt.close(fig) diff --git a/LWB/propagators/figures/morse_1D_Hagedorn_drift.png b/LWB/propagators/figures/morse_1D_Hagedorn_drift.png new file mode 100644 index 0000000..eb14dbc Binary files /dev/null and b/LWB/propagators/figures/morse_1D_Hagedorn_drift.png differ diff --git a/LWB/propagators/figures/morse_1D_Hagedorn_energies.png b/LWB/propagators/figures/morse_1D_Hagedorn_energies.png new file mode 100644 index 0000000..2a555ff Binary files /dev/null and b/LWB/propagators/figures/morse_1D_Hagedorn_energies.png differ diff --git a/LWB/propagators/figures/morse_1D_MG4_drift.png b/LWB/propagators/figures/morse_1D_MG4_drift.png new file mode 100644 index 0000000..977fa1d Binary files /dev/null and b/LWB/propagators/figures/morse_1D_MG4_drift.png differ diff --git a/LWB/propagators/figures/morse_1D_MG4_energies.png b/LWB/propagators/figures/morse_1D_MG4_energies.png new file mode 100644 index 0000000..01e9283 Binary files /dev/null and b/LWB/propagators/figures/morse_1D_MG4_energies.png differ diff --git a/LWB/propagators/figures/morse_1D_McL42_drift.png b/LWB/propagators/figures/morse_1D_McL42_drift.png new file mode 100644 index 0000000..be702cf Binary files /dev/null and b/LWB/propagators/figures/morse_1D_McL42_drift.png differ diff --git a/LWB/propagators/figures/morse_1D_McL42_energies.png b/LWB/propagators/figures/morse_1D_McL42_energies.png new file mode 100644 index 0000000..01e9283 Binary files /dev/null and b/LWB/propagators/figures/morse_1D_McL42_energies.png differ diff --git a/LWB/propagators/figures/morse_1D_McL84_drift.png b/LWB/propagators/figures/morse_1D_McL84_drift.png new file mode 100644 index 0000000..977fa1d Binary files /dev/null and b/LWB/propagators/figures/morse_1D_McL84_drift.png differ diff --git a/LWB/propagators/figures/morse_1D_McL84_energies.png b/LWB/propagators/figures/morse_1D_McL84_energies.png new file mode 100644 index 0000000..01e9283 Binary files /dev/null and b/LWB/propagators/figures/morse_1D_McL84_energies.png differ diff --git a/LWB/propagators/figures/morse_1D_Pre764_drift.png b/LWB/propagators/figures/morse_1D_Pre764_drift.png new file mode 100644 index 0000000..2417f12 Binary files /dev/null and b/LWB/propagators/figures/morse_1D_Pre764_drift.png differ diff --git a/LWB/propagators/figures/morse_1D_Pre764_energies.png b/LWB/propagators/figures/morse_1D_Pre764_energies.png new file mode 100644 index 0000000..c071786 Binary files /dev/null and b/LWB/propagators/figures/morse_1D_Pre764_energies.png differ diff --git a/LWB/propagators/figures/morse_1D_Semiclassical_drift.png b/LWB/propagators/figures/morse_1D_Semiclassical_drift.png new file mode 100644 index 0000000..fe51568 Binary files /dev/null and b/LWB/propagators/figures/morse_1D_Semiclassical_drift.png differ diff --git a/LWB/propagators/figures/morse_1D_Semiclassical_energies.png b/LWB/propagators/figures/morse_1D_Semiclassical_energies.png new file mode 100644 index 0000000..01e9283 Binary files /dev/null and b/LWB/propagators/figures/morse_1D_Semiclassical_energies.png differ diff --git a/LWB/propagators/figures/old/coefficient_analysis.png b/LWB/propagators/figures/old/coefficient_analysis.png new file mode 100644 index 0000000..3cbc428 Binary files /dev/null and b/LWB/propagators/figures/old/coefficient_analysis.png differ diff --git a/LWB/propagators/figures/old/dimension_analysis.png b/LWB/propagators/figures/old/dimension_analysis.png new file mode 100644 index 0000000..a17a7b9 Binary files /dev/null and b/LWB/propagators/figures/old/dimension_analysis.png differ diff --git a/LWB/propagators/figures/old/error_analysis.png b/LWB/propagators/figures/old/error_analysis.png new file mode 100644 index 0000000..06299b7 Binary files /dev/null and b/LWB/propagators/figures/old/error_analysis.png differ diff --git a/LWB/propagators/figures/old/error_analysis_lt.png b/LWB/propagators/figures/old/error_analysis_lt.png new file mode 100644 index 0000000..d3334ed Binary files /dev/null and b/LWB/propagators/figures/old/error_analysis_lt.png differ diff --git a/LWB/propagators/figures/old/error_analysis_y4.png b/LWB/propagators/figures/old/error_analysis_y4.png new file mode 100644 index 0000000..dee0616 Binary files /dev/null and b/LWB/propagators/figures/old/error_analysis_y4.png differ diff --git a/LWB/propagators/figures/old/harmonic_drift_Hagedorn.png b/LWB/propagators/figures/old/harmonic_drift_Hagedorn.png new file mode 100644 index 0000000..240b256 Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_drift_Hagedorn.png differ diff --git a/LWB/propagators/figures/old/harmonic_drift_MG4.png b/LWB/propagators/figures/old/harmonic_drift_MG4.png new file mode 100644 index 0000000..de8d027 Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_drift_MG4.png differ diff --git a/LWB/propagators/figures/old/harmonic_drift_McL42.png b/LWB/propagators/figures/old/harmonic_drift_McL42.png new file mode 100644 index 0000000..300356e Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_drift_McL42.png differ diff --git a/LWB/propagators/figures/old/harmonic_drift_McL84.png b/LWB/propagators/figures/old/harmonic_drift_McL84.png new file mode 100644 index 0000000..7a9e7f4 Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_drift_McL84.png differ diff --git a/LWB/propagators/figures/old/harmonic_drift_Pre764.png b/LWB/propagators/figures/old/harmonic_drift_Pre764.png new file mode 100644 index 0000000..724ad1b Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_drift_Pre764.png differ diff --git a/LWB/propagators/figures/old/harmonic_drift_Semiclassical.png b/LWB/propagators/figures/old/harmonic_drift_Semiclassical.png new file mode 100644 index 0000000..58b2e1d Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_drift_Semiclassical.png differ diff --git a/LWB/propagators/figures/old/harmonic_energies_Hagedorn.png b/LWB/propagators/figures/old/harmonic_energies_Hagedorn.png new file mode 100644 index 0000000..7f7f9fa Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_energies_Hagedorn.png differ diff --git a/LWB/propagators/figures/old/harmonic_energies_MG4.png b/LWB/propagators/figures/old/harmonic_energies_MG4.png new file mode 100644 index 0000000..86a00d3 Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_energies_MG4.png differ diff --git a/LWB/propagators/figures/old/harmonic_energies_McL42.png b/LWB/propagators/figures/old/harmonic_energies_McL42.png new file mode 100644 index 0000000..5397c62 Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_energies_McL42.png differ diff --git a/LWB/propagators/figures/old/harmonic_energies_McL84.png b/LWB/propagators/figures/old/harmonic_energies_McL84.png new file mode 100644 index 0000000..a4201b5 Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_energies_McL84.png differ diff --git a/LWB/propagators/figures/old/harmonic_energies_Pre764.png b/LWB/propagators/figures/old/harmonic_energies_Pre764.png new file mode 100644 index 0000000..b87826d Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_energies_Pre764.png differ diff --git a/LWB/propagators/figures/old/harmonic_energies_Semiclassical.png b/LWB/propagators/figures/old/harmonic_energies_Semiclassical.png new file mode 100644 index 0000000..06d19ef Binary files /dev/null and b/LWB/propagators/figures/old/harmonic_energies_Semiclassical.png differ diff --git a/LWB/propagators/figures/old/henon_drift_Hagedorn.png b/LWB/propagators/figures/old/henon_drift_Hagedorn.png new file mode 100644 index 0000000..cc67a6b Binary files /dev/null and b/LWB/propagators/figures/old/henon_drift_Hagedorn.png differ diff --git a/LWB/propagators/figures/old/henon_drift_MG4.png b/LWB/propagators/figures/old/henon_drift_MG4.png new file mode 100644 index 0000000..9b3660c Binary files /dev/null and b/LWB/propagators/figures/old/henon_drift_MG4.png differ diff --git a/LWB/propagators/figures/old/henon_drift_McL42.png b/LWB/propagators/figures/old/henon_drift_McL42.png new file mode 100644 index 0000000..911bac8 Binary files /dev/null and b/LWB/propagators/figures/old/henon_drift_McL42.png differ diff --git a/LWB/propagators/figures/old/henon_drift_McL84.png b/LWB/propagators/figures/old/henon_drift_McL84.png new file mode 100644 index 0000000..1f35ad6 Binary files /dev/null and b/LWB/propagators/figures/old/henon_drift_McL84.png differ diff --git a/LWB/propagators/figures/old/henon_drift_Pre764.png b/LWB/propagators/figures/old/henon_drift_Pre764.png new file mode 100644 index 0000000..1d9cf8f Binary files /dev/null and b/LWB/propagators/figures/old/henon_drift_Pre764.png differ diff --git a/LWB/propagators/figures/old/henon_drift_Semiclassical.png b/LWB/propagators/figures/old/henon_drift_Semiclassical.png new file mode 100644 index 0000000..5c0f4ff Binary files /dev/null and b/LWB/propagators/figures/old/henon_drift_Semiclassical.png differ diff --git a/LWB/propagators/figures/old/henon_energies_Hagedorn.png b/LWB/propagators/figures/old/henon_energies_Hagedorn.png new file mode 100644 index 0000000..e8a7c25 Binary files /dev/null and b/LWB/propagators/figures/old/henon_energies_Hagedorn.png differ diff --git a/LWB/propagators/figures/old/henon_energies_MG4.png b/LWB/propagators/figures/old/henon_energies_MG4.png new file mode 100644 index 0000000..bef478c Binary files /dev/null and b/LWB/propagators/figures/old/henon_energies_MG4.png differ diff --git a/LWB/propagators/figures/old/henon_energies_McL42.png b/LWB/propagators/figures/old/henon_energies_McL42.png new file mode 100644 index 0000000..285a9fa Binary files /dev/null and b/LWB/propagators/figures/old/henon_energies_McL42.png differ diff --git a/LWB/propagators/figures/old/henon_energies_McL84.png b/LWB/propagators/figures/old/henon_energies_McL84.png new file mode 100644 index 0000000..224b9cd Binary files /dev/null and b/LWB/propagators/figures/old/henon_energies_McL84.png differ diff --git a/LWB/propagators/figures/old/henon_energies_Pre764.png b/LWB/propagators/figures/old/henon_energies_Pre764.png new file mode 100644 index 0000000..62776eb Binary files /dev/null and b/LWB/propagators/figures/old/henon_energies_Pre764.png differ diff --git a/LWB/propagators/figures/old/henon_energies_Semiclassical.png b/LWB/propagators/figures/old/henon_energies_Semiclassical.png new file mode 100644 index 0000000..0c1f8ac Binary files /dev/null and b/LWB/propagators/figures/old/henon_energies_Semiclassical.png differ diff --git a/LWB/propagators/figures/rubbish/drift_Hagedorn.png b/LWB/propagators/figures/rubbish/drift_Hagedorn.png new file mode 100644 index 0000000..b485d75 Binary files /dev/null and b/LWB/propagators/figures/rubbish/drift_Hagedorn.png differ diff --git a/LWB/propagators/figures/rubbish/drift_MG4.png b/LWB/propagators/figures/rubbish/drift_MG4.png new file mode 100644 index 0000000..afb544f Binary files /dev/null and b/LWB/propagators/figures/rubbish/drift_MG4.png differ diff --git a/LWB/propagators/figures/rubbish/drift_McL42.png b/LWB/propagators/figures/rubbish/drift_McL42.png new file mode 100644 index 0000000..ed16bbe Binary files /dev/null and b/LWB/propagators/figures/rubbish/drift_McL42.png differ diff --git a/LWB/propagators/figures/rubbish/drift_McL84.png b/LWB/propagators/figures/rubbish/drift_McL84.png new file mode 100644 index 0000000..08f5fb5 Binary files /dev/null and b/LWB/propagators/figures/rubbish/drift_McL84.png differ diff --git a/LWB/propagators/figures/rubbish/drift_Pre764.png b/LWB/propagators/figures/rubbish/drift_Pre764.png new file mode 100644 index 0000000..1921736 Binary files /dev/null and b/LWB/propagators/figures/rubbish/drift_Pre764.png differ diff --git a/LWB/propagators/figures/rubbish/drift_Semiclassical.png b/LWB/propagators/figures/rubbish/drift_Semiclassical.png new file mode 100644 index 0000000..4c07a83 Binary files /dev/null and b/LWB/propagators/figures/rubbish/drift_Semiclassical.png differ diff --git a/LWB/propagators/figures/rubbish/energy_Hagedorn.png b/LWB/propagators/figures/rubbish/energy_Hagedorn.png new file mode 100644 index 0000000..96f7e76 Binary files /dev/null and b/LWB/propagators/figures/rubbish/energy_Hagedorn.png differ diff --git a/LWB/propagators/figures/rubbish/energy_MG4.png b/LWB/propagators/figures/rubbish/energy_MG4.png new file mode 100644 index 0000000..b91f673 Binary files /dev/null and b/LWB/propagators/figures/rubbish/energy_MG4.png differ diff --git a/LWB/propagators/figures/rubbish/energy_McL42.png b/LWB/propagators/figures/rubbish/energy_McL42.png new file mode 100644 index 0000000..b91f673 Binary files /dev/null and b/LWB/propagators/figures/rubbish/energy_McL42.png differ diff --git a/LWB/propagators/figures/rubbish/energy_McL84.png b/LWB/propagators/figures/rubbish/energy_McL84.png new file mode 100644 index 0000000..e731a68 Binary files /dev/null and b/LWB/propagators/figures/rubbish/energy_McL84.png differ diff --git a/LWB/propagators/figures/rubbish/energy_Pre764.png b/LWB/propagators/figures/rubbish/energy_Pre764.png new file mode 100644 index 0000000..2dd9907 Binary files /dev/null and b/LWB/propagators/figures/rubbish/energy_Pre764.png differ diff --git a/LWB/propagators/figures/rubbish/energy_Semiclassical.png b/LWB/propagators/figures/rubbish/energy_Semiclassical.png new file mode 100644 index 0000000..b91f673 Binary files /dev/null and b/LWB/propagators/figures/rubbish/energy_Semiclassical.png differ diff --git a/LWB/propagators/figures/torsional_1D_Hagedorn_drift.png b/LWB/propagators/figures/torsional_1D_Hagedorn_drift.png new file mode 100644 index 0000000..6761603 Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_Hagedorn_drift.png differ diff --git a/LWB/propagators/figures/torsional_1D_Hagedorn_energies.png b/LWB/propagators/figures/torsional_1D_Hagedorn_energies.png new file mode 100644 index 0000000..4356132 Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_Hagedorn_energies.png differ diff --git a/LWB/propagators/figures/torsional_1D_MG4_drift.png b/LWB/propagators/figures/torsional_1D_MG4_drift.png new file mode 100644 index 0000000..ea23180 Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_MG4_drift.png differ diff --git a/LWB/propagators/figures/torsional_1D_MG4_energies.png b/LWB/propagators/figures/torsional_1D_MG4_energies.png new file mode 100644 index 0000000..e7e169e Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_MG4_energies.png differ diff --git a/LWB/propagators/figures/torsional_1D_McL42_drift.png b/LWB/propagators/figures/torsional_1D_McL42_drift.png new file mode 100644 index 0000000..16d239a Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_McL42_drift.png differ diff --git a/LWB/propagators/figures/torsional_1D_McL42_energies.png b/LWB/propagators/figures/torsional_1D_McL42_energies.png new file mode 100644 index 0000000..e7e169e Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_McL42_energies.png differ diff --git a/LWB/propagators/figures/torsional_1D_McL84_drift.png b/LWB/propagators/figures/torsional_1D_McL84_drift.png new file mode 100644 index 0000000..5e49298 Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_McL84_drift.png differ diff --git a/LWB/propagators/figures/torsional_1D_McL84_energies.png b/LWB/propagators/figures/torsional_1D_McL84_energies.png new file mode 100644 index 0000000..e7e169e Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_McL84_energies.png differ diff --git a/LWB/propagators/figures/torsional_1D_Pre764_drift.png b/LWB/propagators/figures/torsional_1D_Pre764_drift.png new file mode 100644 index 0000000..3c88234 Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_Pre764_drift.png differ diff --git a/LWB/propagators/figures/torsional_1D_Pre764_energies.png b/LWB/propagators/figures/torsional_1D_Pre764_energies.png new file mode 100644 index 0000000..ccd7427 Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_Pre764_energies.png differ diff --git a/LWB/propagators/figures/torsional_1D_Semiclassical_drift.png b/LWB/propagators/figures/torsional_1D_Semiclassical_drift.png new file mode 100644 index 0000000..4840f3a Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_Semiclassical_drift.png differ diff --git a/LWB/propagators/figures/torsional_1D_Semiclassical_energies.png b/LWB/propagators/figures/torsional_1D_Semiclassical_energies.png new file mode 100644 index 0000000..e7e169e Binary files /dev/null and b/LWB/propagators/figures/torsional_1D_Semiclassical_energies.png differ diff --git a/LWB/propagators/figures/torsional_2D_Semiclassical_drift.png b/LWB/propagators/figures/torsional_2D_Semiclassical_drift.png new file mode 100644 index 0000000..88b107b Binary files /dev/null and b/LWB/propagators/figures/torsional_2D_Semiclassical_drift.png differ diff --git a/LWB/propagators/figures/torsional_2D_Semiclassical_energies.png b/LWB/propagators/figures/torsional_2D_Semiclassical_energies.png new file mode 100644 index 0000000..00c3dff Binary files /dev/null and b/LWB/propagators/figures/torsional_2D_Semiclassical_energies.png differ diff --git a/LWB/propagators/figures/torsional_Hagedorn_drift.png b/LWB/propagators/figures/torsional_Hagedorn_drift.png new file mode 100644 index 0000000..32569db Binary files /dev/null and b/LWB/propagators/figures/torsional_Hagedorn_drift.png differ diff --git a/LWB/propagators/figures/torsional_Hagedorn_energies.png b/LWB/propagators/figures/torsional_Hagedorn_energies.png new file mode 100644 index 0000000..a2b5546 Binary files /dev/null and b/LWB/propagators/figures/torsional_Hagedorn_energies.png differ diff --git a/LWB/propagators/figures/torsional_MG4_drift.png b/LWB/propagators/figures/torsional_MG4_drift.png new file mode 100644 index 0000000..cd44cab Binary files /dev/null and b/LWB/propagators/figures/torsional_MG4_drift.png differ diff --git a/LWB/propagators/figures/torsional_MG4_energies.png b/LWB/propagators/figures/torsional_MG4_energies.png new file mode 100644 index 0000000..00c3dff Binary files /dev/null and b/LWB/propagators/figures/torsional_MG4_energies.png differ diff --git a/LWB/propagators/figures/torsional_McL42_drift.png b/LWB/propagators/figures/torsional_McL42_drift.png new file mode 100644 index 0000000..e199600 Binary files /dev/null and b/LWB/propagators/figures/torsional_McL42_drift.png differ diff --git a/LWB/propagators/figures/torsional_McL42_energies.png b/LWB/propagators/figures/torsional_McL42_energies.png new file mode 100644 index 0000000..00c3dff Binary files /dev/null and b/LWB/propagators/figures/torsional_McL42_energies.png differ diff --git a/LWB/propagators/figures/torsional_McL84_drift.png b/LWB/propagators/figures/torsional_McL84_drift.png new file mode 100644 index 0000000..29ae311 Binary files /dev/null and b/LWB/propagators/figures/torsional_McL84_drift.png differ diff --git a/LWB/propagators/figures/torsional_McL84_energies.png b/LWB/propagators/figures/torsional_McL84_energies.png new file mode 100644 index 0000000..00c3dff Binary files /dev/null and b/LWB/propagators/figures/torsional_McL84_energies.png differ diff --git a/LWB/propagators/figures/torsional_Pre764_drift.png b/LWB/propagators/figures/torsional_Pre764_drift.png new file mode 100644 index 0000000..2b5293b Binary files /dev/null and b/LWB/propagators/figures/torsional_Pre764_drift.png differ diff --git a/LWB/propagators/figures/torsional_Pre764_energies.png b/LWB/propagators/figures/torsional_Pre764_energies.png new file mode 100644 index 0000000..cd37fc0 Binary files /dev/null and b/LWB/propagators/figures/torsional_Pre764_energies.png differ diff --git a/LWB/propagators/figures/torsional_Semiclassical_drift.png b/LWB/propagators/figures/torsional_Semiclassical_drift.png new file mode 100644 index 0000000..88b107b Binary files /dev/null and b/LWB/propagators/figures/torsional_Semiclassical_drift.png differ diff --git a/LWB/propagators/figures/torsional_Semiclassical_energies.png b/LWB/propagators/figures/torsional_Semiclassical_energies.png new file mode 100644 index 0000000..00c3dff Binary files /dev/null and b/LWB/propagators/figures/torsional_Semiclassical_energies.png differ diff --git a/LWB/propagators/figures/uml.png b/LWB/propagators/figures/uml.png new file mode 100644 index 0000000..ef02792 Binary files /dev/null and b/LWB/propagators/figures/uml.png differ diff --git a/LWB/propagators/header.tex b/LWB/propagators/header.tex new file mode 100644 index 0000000..668fd13 --- /dev/null +++ b/LWB/propagators/header.tex @@ -0,0 +1,209 @@ +\documentclass[a4paper,10pt]{scrartcl} + +\usepackage{standalone} + +% Load ``float'' before ``hyperref`` before ''algorithm`` +% Note: ''algorithm`` would load ''float`` by itself +\usepackage{float} +\usepackage[pagebackref,hyperindex=true]{hyperref} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Math +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{amsfonts} +\usepackage{amsopn} +\usepackage{braket} +\usepackage{bbm} +\usepackage{dsfont} +% \usepackage{mathabx} + + +% Various new commands that ease typesetting math even further +% \newcommand{\assign}{\ensuremath{\coloneq}} +% \newcommand{\rassign}{\ensuremath{\eqcolon}} +% TODO: remove +% \newcommand{\assign}{\ensuremath{:=}} +% \newcommand{\rassign}{\ensuremath{=:}} +% \newcommand{\seteq}{\ensuremath{\overset{!}{=}}} +% \newcommand{\of}[1]{\ensuremath{\left( #1 \right)}} +% \newcommand{\ofs}[1]{\ensuremath{\left( #1 \right)}} +% \newcommand{\norm}[1]{\ensuremath{\| #1 \|}} +% \newcommand{\tmop}[1]{\ensuremath{\operatorname{#1}}} +% \newcommand{\id}{\ensuremath{\mathds{1}}} +% % \newcommand{\id}{\ensuremath{I}} +% \newcommand{\kron}[1]{\ensuremath{\delta_{#1}}} +% \newcommand{\conj}[1]{\ensuremath{\overline{#1}}} +% \renewcommand{\vec}[1]{\ensuremath{\underline{#1}}} +% \newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}} +% \newcommand{\inv}{\ensuremath{{}^{-1}}} +% \newcommand{\T}{\ensuremath{{}^{\textnormal{T}}}} +% \renewcommand{\H}{\ensuremath{{}^{\textnormal{H}}}} +% \newcommand{\Tinv}{\ensuremath{{}^{\textnormal{-T}}}} +% \newcommand{\Hinv}{\ensuremath{{}^{\textnormal{-H}}}} +% \newcommand{\tr}{\ensuremath{\textnormal{Tr}}} +% \newcommand{\ft}[1]{\ensuremath{\mathcal{F}\left(#1\right)}} +% \newcommand{\ift}[1]{\ensuremath{\mathcal{F}^{-1}\left(#1\right)}} +% \newcommand{\fft}[1]{\ensuremath{\mathtt{FFT}\left(#1\right)}} +% \newcommand{\ifft}[1]{\ensuremath{\mathtt{IFFT}\left(#1\right)}} +% \newcommand{\dotp}[2]{\ensuremath{\left\langle #1 , #2 \right\rangle}} +% \newcommand{\bigO}[1]{\ensuremath{\mathcal{O}\left( #1 \right)}} +% \newcommand{\laplace}{\ensuremath{\operatorname{\Delta}}} +% \newcommand{\diff}[3][]{\frac{\mathrm{d}^{#1}#2}{\mathrm{d}#3^{#1}}} +% \newcommand{\pdiff}[3][]{\frac{\partial^{#1}#2}{\partial #3^{#1}}} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Graphics +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\usepackage{graphicx} +\usepackage{subfig} +\usepackage{asymptote} +\usepackage{tikz} +\usetikzlibrary{topaths,calc} +\usetikzlibrary{positioning} +\usetikzlibrary{arrows,shapes,backgrounds} +\definecolor{gray}{gray}{0.55} + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Math +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\usepackage{algorithm} +% \usepackage{algorithmic} +\usepackage{color} +\usepackage{relsize} +\usepackage[scaled=0.8]{beramono} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Code +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\usepackage{listings} +% C++ Code macro +\newcommand{\cpp}[1] { + \lstset{language=c++, + basicstyle=\smaller, + basewidth=0.58em, + columns=fixed, + tabsize=2, + fontadjust=true, + frame=l, + xleftmargin=4.2pt, + numbers=left, + stepnumber=2, + breaklines=true, + breakindent=0pt, + prebreak=\mbox{\tiny$\searrow$}, + postbreak=\mbox{{\color{gray}$\cdots$}}, + numberstyle=\color{gray}, + commentstyle=\color{gray}, + stringstyle=\textit, + showstringspaces=false, + } + \lstinputlisting{#1} +} + + \lstset{language=c++, + basicstyle=\smaller, + basewidth=0.58em, + columns=fixed, + tabsize=2, + fontadjust=true, + breaklines=true, + breakindent=0pt, + prebreak=\mbox{\tiny$\searrow$}, + postbreak=\mbox{{\color{gray}$\cdots$}}, + numberstyle=\color{gray}, + commentstyle=\color{gray}, + stringstyle=\textit, + showstringspaces=false, + basicstyle=\footnotesize + } + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% More Packages +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\usepackage[utf8x]{inputenc} +\usepackage{cancel} +\usepackage{url} +\usepackage{placeins} +\usepackage{microtype} +\usepackage{amsthm} +\newtheorem{theorem}{Theorem} +\newtheorem{definition}{Definition} + +% Links in pdf +\usepackage{color} +\newcommand{\blue}{ \color{blue} } +\definecolor{linkcol}{rgb}{0,0,0.4} +\definecolor{citecol}{rgb}{0.5,0,0} + +\hypersetup +{ +colorlinks=true, +linkcolor=linkcol, +citecolor=citecol, +urlcolor=linkcol +} + +% nicer backref links +\renewcommand*{\backref}[1]{} +\renewcommand*{\backrefalt}[4]{% +\ifcase #1 % +(Not cited.)% +\or +(Cited on page~#2.)% +\else +(Cited on pages~#2.)% +\fi} +\renewcommand*{\backrefsep}{, } +\renewcommand*{\backreftwosep}{ and~} +\renewcommand*{\backreflastsep}{ and~} + +% TODO: set some indent +% \parindent 0cm + +\newcommand{\clearemptydoublepage}{\newpage{\pagestyle{empty}\cleardoublepage}} + + +% Mycommands + +\usepackage{bm} + +% \newcommand*{\op}[1]{\hat{\mathcal{#1}}} +\newcommand*{\C}{\ensuremath{\mathbb{C}}} +\newcommand*{\Dt}{\ensuremath{\Delta t}} +\newcommand*{\K}{\ensuremath{\mathfrak{K}}} +\newcommand*{\R}{\ensuremath{\mathbb{R}}} +\newcommand*{\bmat}[1]{\bm{#1}} +\newcommand*{\bvec}[1]{\bm{#1}} +\newcommand*{\conj}[1]{\overline{#1}} +\newcommand*{\del}{\ensuremath{\partial}} +\newcommand*{\dif}{\ensuremath{d}} +\newcommand*{\dt}{\ensuremath{\delta t}} +\newcommand*{\eps}{\ensuremath{\varepsilon}} +\newcommand*{\im}{\ensuremath{\imath}} +\newcommand*{\opH}{\mathcal{H}} +\newcommand*{\opT}{\mathcal{T}} +\newcommand*{\opU}{\mathcal{U}} +\newcommand*{\opV}{\mathcal{V}} +\newcommand*{\opW}{\mathcal{W}} +\newcommand*{\opA}{\mathcal{A}} +\newcommand*{\opB}{\mathcal{B}} +\newcommand*{\opX}{\mathcal{X}} +\newcommand*{\upic}{\ensuremath{u[\Pi,\bm{c}]}} +\newcommand*{\proc}[1]{\textsc{#1}} + +\usepackage{algorithm} +\usepackage{verbatim} +\usepackage{physics} +\usepackage{amssymb} +\usepackage{bm} +% \usepackage[noend]{algpseudocode} +\usepackage{algpseudocode} + + +\usepackage{tikz} +\usetikzlibrary{arrows,chains,matrix,positioning,scopes} +\usepackage{multicol} diff --git a/LWB/propagators/main.tex b/LWB/propagators/main.tex new file mode 100644 index 0000000..cac4275 --- /dev/null +++ b/LWB/propagators/main.tex @@ -0,0 +1,33 @@ +\input{./header.tex} + +\begin{document} + +\input{sections/title.tex} + +\clearemptydoublepage + +\mbox{} \vfill \input{sections/abstract.tex} \vfill +\tableofcontents + +\clearpage + +% \listoffigures +\listofalgorithms + + +\clearemptydoublepage \input{sections/introduction.tex} +\clearpage \input{sections/operatorsplitting.tex} +\clearpage \input{sections/buildingblocks.tex} +\clearpage \input{sections/propagators.tex} +\clearpage \input{sections/implementation.tex} +\clearpage \input{sections/results.tex} +\clearpage \input{sections/conclusion.tex} + +\clearemptydoublepage \input{sections/appendix.tex} + +\clearemptydoublepage + +\bibliographystyle{plain} +\bibliography{references,mt,bt,chem,wp,own,propagators} + +\end{document} diff --git a/LWB/propagators/sections/abstract.tex b/LWB/propagators/sections/abstract.tex new file mode 100644 index 0000000..1934905 --- /dev/null +++ b/LWB/propagators/sections/abstract.tex @@ -0,0 +1,11 @@ +\begin{abstract} + This article explains the central design decisions and most important code details related to the efficient implementation of time propagation for Hagedorn wave packets in the WaveBlocks framework. + A short introduction to operator splitting and a summary of some of the most important time propagation schemes for Hagedorn wave packets are also given. + Finally, some benchmark results are presented which analyze energy conservation and compare different propagators, splitting schemes, step sizes and dimensionalities. +\end{abstract} +% +\vspace{5mm} +% +\begin{center} + The code to the project can be found on GitHub at \cite{libwaveblocks}. +\end{center} diff --git a/LWB/propagators/sections/appendix.tex b/LWB/propagators/sections/appendix.tex new file mode 100644 index 0000000..daed1b3 --- /dev/null +++ b/LWB/propagators/sections/appendix.tex @@ -0,0 +1,105 @@ +\appendix + +\section{Energy Evolution and Drift Plots} +% +The following images show the energy and energy drift plots for all the propagators of this report. + +\begin{figure}[ht] + \centering + \begin{minipage}[c]{\textwidth} + \begin{center} + \large Hagedorn Propagator \\[1mm] + \normalsize Energy Evolution and Drift + \vspace{4mm} + \end{center} + \end{minipage} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_Hagedorn_energies.png} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_Hagedorn_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/torsional_1D_Hagedorn_energies.png} + \includegraphics[width=.45\textwidth]{figures/torsional_1D_Hagedorn_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/morse_1D_Hagedorn_energies.png} + \includegraphics[width=.45\textwidth]{figures/morse_1D_Hagedorn_drift.png} + \caption{Energy evolution and drift for a 1D wave packet propagated with the Hagedorn propagator in a harmonic potential (top), a torsional potential (middle) and a Morse potential (bottom). + (Parameters: $N=1$, $D=1$ $|\K|=16$ $\eps=0.01$ (Morse $0.0484$), $T=10$ (Morse $T=50$), $\Dt=0.001$ (Morse $\Dt=0.005$), Hagedorn propagator)} + \label{fig:energy_Hagedorn} +\end{figure} +% +\begin{figure}[ht] + \centering + \begin{minipage}[c]{\textwidth} + \begin{center} + \large MG4 Propagator \\[1mm] + \normalsize Energy Evolution and Drift + \vspace{4mm} + \end{center} + \end{minipage} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_MG4_energies.png} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_MG4_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/torsional_1D_MG4_energies.png} + \includegraphics[width=.45\textwidth]{figures/torsional_1D_MG4_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/morse_1D_MG4_energies.png} + \includegraphics[width=.45\textwidth]{figures/morse_1D_MG4_drift.png} + \caption{Energy evolution and drift for a 1D wave packet propagated with the MG4 propagator in a harmonic potential (top), a torsional potential (middle) and a Morse potential (bottom). + (Parameters: $N=1$, $D=1$ $|\K|=16$ $\eps=0.01$ (Morse $0.0484$), $T=10$ (Morse $T=50$), $\Dt=0.001$ (Morse $\Dt=0.005$), MG4 propagator with \emph{Y4} splitting for \proc{IntSplit})} + \label{fig:energy_MG4} +\end{figure} +% +\begin{figure}[ht] + \centering + \begin{minipage}[c]{\textwidth} + \begin{center} + \large McL42 Propagator \\[1mm] + \normalsize Energy Evolution and Drift + \vspace{4mm} + \end{center} + \end{minipage} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_McL42_energies.png} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_McL42_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/torsional_1D_McL42_energies.png} + \includegraphics[width=.45\textwidth]{figures/torsional_1D_McL42_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/morse_1D_McL42_energies.png} + \includegraphics[width=.45\textwidth]{figures/morse_1D_McL42_drift.png} + \caption{Energy evolution and drift for a 1D wave packet propagated with the McL42 propagator in a harmonic potential (top), a torsional potential (middle) and a Morse potential (bottom). + (Parameters: $N=1$, $D=1$ $|\K|=16$ $\eps=0.01$ (Morse $0.0484$), $T=10$ (Morse $T=50$), $\Dt=0.001$ (Morse $\Dt=0.005$), McL42 propagator with \emph{Y4} splitting for \proc{IntSplit})} + \label{fig:energy_McL42} +\end{figure} +% +\begin{figure}[ht] + \centering + \begin{minipage}[c]{\textwidth} + \begin{center} + \large McL84 Propagator \\[1mm] + \normalsize Energy Evolution and Drift + \vspace{4mm} + \end{center} + \end{minipage} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_McL84_energies.png} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_McL84_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/torsional_1D_McL84_energies.png} + \includegraphics[width=.45\textwidth]{figures/torsional_1D_McL84_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/morse_1D_McL84_energies.png} + \includegraphics[width=.45\textwidth]{figures/morse_1D_McL84_drift.png} + \caption{Energy evolution and drift for a 1D wave packet propagated with the McL84 propagator in a harmonic potential (top), a torsional potential (middle) and a Morse potential (bottom). + (Parameters: $N=1$, $D=1$ $|\K|=16$ $\eps=0.01$ (Morse $0.0484$), $T=10$ (Morse $T=50$), $\Dt=0.001$ (Morse $\Dt=0.005$), McL84 propagator with \emph{Y4} splitting for \proc{IntSplit})} + \label{fig:energy_McL84} +\end{figure} +% +\begin{figure}[ht] + \centering + \begin{minipage}[c]{\textwidth} + \begin{center} + \large Pre764 Propagator \\[1mm] + \normalsize Energy Evolution and Drift + \vspace{4mm} + \end{center} + \end{minipage} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_Pre764_energies.png} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_Pre764_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/torsional_1D_Pre764_energies.png} + \includegraphics[width=.45\textwidth]{figures/torsional_1D_Pre764_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/morse_1D_Pre764_energies.png} + \includegraphics[width=.45\textwidth]{figures/morse_1D_Pre764_drift.png} + \caption{Energy evolution and drift for a 1D wave packet propagated with the Pre764 propagator in a harmonic potential (top), a torsional potential (middle) and a Morse potential (bottom). + (Parameters: $N=1$, $D=1$ $|\K|=16$ $\eps=0.01$ (Morse $0.0484$), $T=10$ (Morse $T=50$), $\Dt=0.001$ (Morse $\Dt=0.005$), Pre764 propagator with \emph{Y4} splitting for \proc{IntSplit})} + \label{fig:energy_Pre764} +\end{figure} diff --git a/LWB/propagators/sections/buildingblocks.tex b/LWB/propagators/sections/buildingblocks.tex new file mode 100644 index 0000000..2b47706 --- /dev/null +++ b/LWB/propagators/sections/buildingblocks.tex @@ -0,0 +1,174 @@ +\section{Common building blocks for time evolution schemes} +\label{sec:buildingblocks} +% +In order to be able to create an efficient, still flexible, implementation of various time propagators, it is helpful to break them down in common, optimized building blocks in terms of which all (or most) propagators can be expressed. +The present section will briefly introduce the most important, generic building blocks that were found necessary in order to build the propagators in section~\ref{sec:propagators}. +For each of them, it will be assumed that the following wave packet variables are implicitly accessible, without passing them as an argument to the procedure: the parameter pack $\bvec{\Pi} = (\bvec{q},\bvec{p},\bvec{Q},\bvec{P})$, the phase factor $S$, the coefficients $\{ c_k \}_{k \in \K}$, the potential function $V(\bvec{q})$ (including Jacobian and Hessian thereof) and the potential remainder $W(\bvec{q})$. +\par\medskip +% +Many of the building blocks depend on the number of energy levels and on whether the wave packet is homogeneous or inhomogeneous. +If any of these cases require a special implementation, it is denoted in the signature of the corresponding algorithm. + + +\subsection{Evolve} +\label{subsec:evolve} +The convenience function \proc{Evolve} (algorithm~\ref{alg:Evolve}) is a wrapper that will carry out the \proc{PrePropagate}, \proc{Propagate} and \proc{PostPropagate} functions. +The concrete implementation of the methods \proc{PrePropagate}, \proc{Propagate} and \proc{PostPropagate} is delegated to the derived classes (see section~\ref{sec:implementation} for more details), +but will generally be a combination of the basic building blocks that follow in this section. +% +\begin{algorithm}[ht] + \caption{Evolve the wave packet for a time period $T$} + \label{alg:Evolve} + \begin{algorithmic} + \State + \Procedure{Evolve}{$T, \Dt$} + \State + \State $M_{tot} = T/\Dt$ + \State \Call{PrePropagate}{$\Dt$} + \For{$m = 1, \dots, M_{tot}$} + \State \Call{Propagate}{$\Dt$} + \EndFor + \State \Call{PostPropagate}{$\Dt$} + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} + + +\subsection{StepT, StepU and IntSplit} +\label{subsec:tuintsplit} +% +%%% StepT & StepU +The time stepping with operators $\opT$ and $\opU = U(\bvec{x})$ follows directly from the propositions in \cite{FGL_semiclassical_dynamics} and is outlined in the algorithms~\ref{alg:StepT} and~\ref{alg:StepU} respectively. +Both functions take a variable $\xi$ as an argument which describes the size of the time step. \\ +In the inhomogeneous implementation of the algorithms, the subscript $n$ is used to denote the relative parameters at energy level $n$. +\par\medskip +% +\begin{algorithm}[ht] + \caption{Propagate with Kinetic Energy Operator $\opT$} + \label{alg:StepT} + \begin{algorithmic} + \State + \Procedure{StepT[homogeneous]}{$\xi$} + \State + \State $\bvec{q} = \bvec{q} + \xi \bvec{p}$ + \Comment update $\bvec{q}$ + \State $\bmat{Q} = \bmat{Q} + \xi \bmat{P}$ + \Comment update $\bmat{Q}$ + \State $S = S + \frac{\xi}{2} \bvec{p}^T \bvec{p}$ + \Comment update $S$ + \State + \EndProcedure + \\\hrulefill + \State + \Procedure{StepT[inhomogeneous]}{$\xi$} + \State + \For{$n=1,...,N$} + \Comment loop over all energy levels + \State $\bvec{q}_n = \bvec{q}_n + \xi \bvec{p}_n$ + \Comment update $\bvec{q}$ for level $n$ + \State $\bmat{Q}_n = \bmat{Q}_n + \xi \bmat{P}_n$ + \Comment update $\bmat{Q}$ for level $n$ + \State $S_n = S_n + \frac{\xi}{2} \bvec{p}_n^T \bvec{p}_n$ + \Comment update $S$ for level $n$ + \EndFor + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} +% +\begin{algorithm}[ht] + \caption{Propagate with (Quadratic) Potential Energy Operator $\opU$} + \label{alg:StepU} + \begin{algorithmic} + \State + \Procedure{StepU[homogeneous]}{$\xi$} + \State + \State $\bvec{p} = \bvec{p} - \xi \nabla V (\bvec{q})$ + \Comment update $\bvec{p}$ + \State $\bvec{P} = \bvec{P} - \xi \nabla^2 V (\bvec{q}) \bvec{Q}$ + \Comment update $\bvec{P}$ + \State $S = S - \xi V (\bvec{q})$ + \Comment update $S$ + \State + \EndProcedure + \\\hrulefill + \State + \Procedure{StepU[inhomogeneous]}{$\xi$} + \State + \For{$n=1,...,N$} + \Comment loop over all energy levels + \State $\bvec{p}_n = \bvec{p}_n - \xi \nabla V (\bvec{q}_n)$ + \Comment update $\bvec{p}$ for level $n$ + \State $\bvec{P}_n = \bvec{P}_n - \xi \nabla^2 V (\bvec{q}_n) \bvec{Q}_n$ + \Comment update $\bvec{P}$ for level $n$ + \State $S_n = S_n - \xi V (\bvec{q}_n)$ + \Comment update $S$ for level $n$ + \EndFor + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} +% +%%% IntSplit +In the practical implementations of the propagators that will be presented in the next section, propagation with $\opT$ and $\opU$ is usually replaced by a series of smaller, alternating steps with the two operators. +In order to facilitate the use of such an alternating pattern, the helper function \proc{IntSplit} (algorithm~\ref{alg:intsplit}) is introduced: +it splits the time step $\Dt$ into $M_{split}$ smaller time steps of size $\dt$ and uses the pair of weight lists $\{ w_T, w_U \}$ to alternatingly apply operators $\opT$ and $\opU$ on each of the sub-intervals. +More detailed information on the individual parameters to the \proc{IntSplit} function and the practical implementation can be found in section \ref{sec:implementation}. +Information about the various coefficient sets $\{w_T, w_U\}$ can be found in the Python source code of the WaveBlocks project \cite{libwaveblocks}. +In this report, the most used splittings are the \emph{LT} splitting (size 1), \emph{Y4} splitting (size 4) and \emph{KL10} splitting (size 34). +% +\begin{algorithm}[ht] + \caption{Split a time interval and alternatingly apply $\opT$ and $\opU$} + \label{alg:intsplit} + \begin{algorithmic} + \State + \Procedure{IntSplit}{$\Dt, M_{split}, \{w_T,w_U\}$} + \State + \State $\dt = \Dt/M_{split}$ + \Comment step size for splitting + \For{$n = 1,\dots,M_{split}$} + \Comment split interval in $M_{split}$ sub-steps + \For{$\alpha$ in $w_T$ and $\beta$ in $w_U$} + \Comment alternatingly apply $\opT$ and $\opU$ + \State \Call{StepT}{$\alpha \cdot \dt$} + \State \Call{StepU}{$\beta \cdot \dt$} + \EndFor + \EndFor + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} + + +\subsection{StepW and building of the interaction matrix} +% +As pointed out in the previous section, the time propagation for non-quadratic potentials $W(\bvec{x})$ (see algorithm~\ref{alg:StepW}) can be achieved by updating only the coefficients $\{ c_k \}_{k \in \K}$. +The update rule is +% +\begin{align} + \bvec{c}(t + \dt) = \exp \left( - \frac{\im t}{\eps^2} \bmat{F} \right) \bvec{c}(t) +\end{align} +% +where the matrix $\bmat{F} = \left[ f_{k,l} \right]_{k,l \in \K}$ has entries +% +\begin{align} + f_{k,l} = \matrixel{\varphi_k}{W}{\varphi_l} + = \int_{\R^D} \conj{\varphi_k(\bvec{x})} W(\bvec{x}) \varphi_l(\bvec{x}) \; \dif \bvec{x} \; . +\end{align} +% +The calculation of such interaction matrices $\bmat{F}$ makes heavy use of the work on inner products that has been presented in \cite{LWB_innerproducts}, another student project in the WaveBlocks framework. +For the further treatment in this report, it is assumed that an efficient implementation of \proc{BuildF} is readily available. +% +\begin{algorithm}[ht] + \caption{Propagate with (Non-Quadratic) Potential Energy Operator $\opW$} + \label{alg:StepW} + \begin{algorithmic} + \State + \Procedure{StepW}{$\xi$} + \State $\bm{\Sigma} = - \im \frac{\xi}{\eps^2} \cdot$ \Call{BuildF}{$\Pi$} + \State $\bm{c} = \exp{\bm{\Sigma}} \bm{c}$ + \EndProcedure + \State + \end{algorithmic} +\end{algorithm} diff --git a/LWB/propagators/sections/conclusion.tex b/LWB/propagators/sections/conclusion.tex new file mode 100644 index 0000000..f0f8371 --- /dev/null +++ b/LWB/propagators/sections/conclusion.tex @@ -0,0 +1,7 @@ +\section{Conclusion} +% +In this project, a series of time propagators for Hagedorn wave packets was studied, analyzed and implemented. +The C++ code, which is the main outcome of the conduced work, is based on a clean software design that makes it very easy to add new Propagators and integrate them into the existing framework. +At the same time, as shown by the benchmark results, the performance of the resulting code is by orders of magnitude better than that of the corresponding Python code. +In particular, the additional cost for using higher-order splitting coefficients is very contained thanks to a templated implementation of the \proc{IntSplit} method which is resolved at compile time and avoids the usage of any loops. +Given these performance improvements, more complex quantum mechanical simulations are now possible with the WaveBlocks framework. diff --git a/LWB/propagators/sections/implementation.tex b/LWB/propagators/sections/implementation.tex new file mode 100644 index 0000000..793cc08 --- /dev/null +++ b/LWB/propagators/sections/implementation.tex @@ -0,0 +1,223 @@ +\section{Implementation in C++} +\label{sec:implementation} +% +This section will explain the most important design decisions and technical details on which the implementation of the presented algorithms is based. +While the first priority was clearly to make the code efficient and reduce computation time, another aim was also to create simple and readable code interfaces that make it easy to use and re-use components and extend the existing functionality. \\ +% +Even though these two requirements often seem to mutually exclude each other, in this case it was possible to meet both targets. +In fact, it is quite remarkable how similar the C++ implementations of the individual propagators look to their pseudo-algorithms as presented in section \ref{sec:propagators}. +\par\medskip +% +Also, it should be mentioned here that some of the most important bottlenecks of the code were not strictly related to the implementation of the propagators, such as for example the building of the $\bmat{F}$-matrix through numerical quadrature of $\matrixel{\varphi_k}{W}{\varphi_l}$. +In fact, the current implementation of the \proc{build\_matrix} routine for computing the inner product requires to always compute a full $N \times N$ sub-matrix, even if just one single entry is needed, which is a very severe inefficiency if the number of levels $N>1$. +\par\medskip +% +However, the aim of the work presented here was to provide a fast implementation of the propagators by re-using existing functions, therefore possible performance enhancements through changing of the existing code was out of the scope of this project and could be targeted in future work. +Consequently, providing an efficient implementation for time propagation mainly translates into finding an efficient way for calling existing functions with minimal overhead and by making use of compiler optimizations where possible. +\par\medskip +% +First, subsection \ref{subsec:poly} outlines how C++ features were used in order to achieve class inheritance and polymorphism without any additional runtime cost. +Secondly, in subsection \ref{subsec:callback}, it will be discussed how the \proc{Evolve} method (see algorithm \ref{alg:Evolve}) was implemented in order to provide a minimal working interface to the user while still preserving all the flexibility that may be required for intermediate measurements or other interaction. +Finally, subsection \ref{subsec:intsplit} dives into the technical details of achieving inlined, highly optimized function calls for the \proc{IntSplit} method (see algorithm \ref{alg:intsplit}) that is used intensively in all the presented propagators. +\par\medskip +% +Some details about the code were not found worthy to be mentioned here, +but more information may be found in the Doxygen documentation. +\par\medskip +% +The Hagedorn Propagator was already implemented in C++ as part of another project on the WaveBlocks library \cite{libwaveblocks}. +Many components of that code were re-used for the present implementation, but they were fundamentally re-organized into a new structure. + + +\subsection{Static Polymorphism through CRTP and enable\_if} +\label{subsec:poly} +% +In the current software design, all propagators from section \ref{sec:propagators} inherit from a common base class \proc{Propagator}, see the UML diagram \ref{fig:uml}. +The strategy for the implementation was the following: +Provide an efficient and generic implementation of all the building blocks from section \ref{sec:buildingblocks} in a common base class (\proc{Propagator}) and then let the derived class decide which combination of building blocks to use for its own, specific implementation of the propagation methods. +This idea is essentially an application of the Strategy Pattern as introduced by Gamma et al. in \cite{Gamma1995}, one of the most important references for software design patterns. +\par\medskip +% +A popular C++ design pattern for achieving static polymorphism is the \emph{Curiously Recurring Template Pattern}, \emph{CRTP} for short, which is described in more detail in \cite{C_CRTP} and in earlier work on the WaveBlocks project \cite{LWB_wavepackets}. +\par\medskip +% +While CRTP is an excellent candidate for providing static polymorphism for derived classes, there was also the need for some polymorphism in the Propagator base class itself, since some of the building blocks need to have different implementations based on the dimensionality or homogeneity of the wave packet. +Providing a partially specialized class for every possible combination of parameters is often tedious and unnecessarily bloats the code. +Therefore, the C++11 \proc{enable\_if} keyword was used, which is a handy C++ meta-function from the \proc{type\_traits} header. +It allows to define multiple versions of a function and conditionally (de)activate them based on a compile-time boolean value. +If the boolean is true, the function is enabled. Otherwise it is omitted. +In the code, such boolean values are indicating whether the number of levels $N$ is equal or greater to one or to distinguish between functions for homogeneous or inhomogeneous wave functions. +The presented use of \proc{enable\_if} effectively provides a form of compile time polymorphism. +% +\begin{figure}[ht] + \centering + \includegraphics[height=10cm]{figures/uml.png} + \caption{Simplified UML diagram for the propagator implementation. + The derived propagator provides implementations for \proc{PrePropagate}, \proc{Propagate} and \proc{PostPropagate} which are called from the public \proc{Evolve} function. + Some functions and function parameters are omitted for better readability of the graph.} + \label{fig:uml} +\end{figure} + + + +\subsection{Evolve and callback function} +\label{subsec:callback} +% +The simple \proc{Evolve} method as described in algorithm \ref{alg:Evolve} in the section about propagator building blocks provides a handy wrapper for calling the \proc{PreProcessing}, \proc{Processing} and \proc{PostProcessing} methods in sequence and dividing the propagation time into smaller steps. \\ +% +However, one may still want to interact with the wave packet during the simulation, for example for writing the wave packet norm to file every 100 steps, throwing errors when energy conservation is violated or to carry out any other task that requires intermediate simulation results. +In order not to lose any flexibility in this regard, an optional callback-function argument was added to the \proc{Evolve} function. +The callback function is passed as a lambda function taking two arguments, the current simulation time \proc{t} and the index of the current time step \proc{m}. +If no callback function is passed, it defaults to an empty function. +\par\medskip +% +Furthermore, since the wave packet is passed to the \proc{Propagator} class \emph{by reference}, it is of course available in the callback function as well and does not need to be passed as an argument explicitly (provided that it is captured in the angular brackets \proc{[]} of the lambda function). \\ +\emph{Disclaimer} if the value of the wave packet is (intentionally or unintentionally) changed in the callback function, the further results of the simulation will of course also be affected. +Also, care must be taken when measuring energies of packets that are being propagated with propagators such as \emph{Pre764} which apply a pre-processor to the packet before starting the simulation. +In order to get the ``physically correct'' wave packet, the corresponding post-processor needs to be applied before making any meaningful measurement. +\par\medskip +% +The callback function is called once prior to simulation start and then in every time step until termination. +\par\medskip +% +Finally, a few helper functions were introduced in order to facilitate pretty printing of console output and some information about the currently propagated wave packet is printed at the beginning of the \proc{Evolve} function. +These console outputs can (and should) of course be removed for optimal performance, but in most cases it was found that some feedback about the current status of the simulation is highly valuable. + + +\subsection{IntSplit using alternating templates} +\label{subsec:intsplit} +% +The \proc{IntSplit} function is merely a helper function that alternatingly applies \proc{StepT} and \proc{StepU}. +Steps with the operators $\opT$ and $\opU$ are usually much cheaper than steps with $\opW$ (which requires numerical integration for approximating the inner product $\matrixel{\varphi_k}{W}{\varphi_l}$), but are in many cases found to be one of the main sources of error. +Therefore the idea of sub-dividing the steps into many smaller intervals - an approach introduced in the semiclassical splitting \ref{sub:semiclassical_propagator} - is applied in many other propagators as well. +Consequently, it is desirable to have a generic and efficient implementation available. +\par\medskip +% +The current \proc{IntSplit} algorithm described in \ref{alg:intsplit} takes three arguments: +the size of the time interval to be covered $\Dt$, the number $M_{split}$ of sub-intervals of size $\dt = \Dt/M_{split}$ in which the initial interval should be divided, and a pair of coefficient lists $\{ w_T, w_U \}$ denoting the weights with which the single steps of size $\dt$ should be scaled. +In the simplest case, the lists $w_T$ and $w_U$ are both only of length one (Lie-Trotter), but higher order splitting strategies can have much larger coefficient sets. +In general, it must hold that $\text{length}(w_T) = \text{length}(w_U)$ (for TU..TU schemes) or $\text{length}(w_T) = \text{length}(w_U)+1$ (for TU..UT schemes). +\par\medskip +% +The obvious implementation of using a dynamic container (like a C++ vector) for the coefficients and a simple \emph{for loop} over all its entries brings a two (minimal) drawbacks: First, the overhead of the loop itself and second, the fact that the number of coefficients is not known at compile time and thus no optimizations can be done. +Given that the functions \proc{StepT} and \proc{StepU} themselves are very cheap (just simple updates of parameters $\bvec{\Pi}$ and $S$), it is preferable to remove any kind of additional cost, even if only marginal. +\par\medskip +% +By making use of C++ arrays and templates, one can remove the issues mentioned above. \\ +% +Firstly, by replacing C++ vectors (dynamically allocated) through C++ arrays (statically allocated), the container size is known at compile time. +For convenience, a new struct \proc{SplitCoefs} was introduced that holds a pair of such arrays. +Secondly, and with the information about container size readily available at compile time, one can then implement two sets of templates \proc{TU} and \proc{UT} that execute a step with operator $\opT$ or $\opU$ respectively and then mutually invoke each other, until both arrays of coefficients are used up. +As a termination condition, when all coefficients have been used, a partially specialized instance of the templates will be called and interrupt the mutual invocation. \\ +% +One technical subtlety of the presented approach is, that the C++ language does not allow partially specialized function templates. This can, however, be overcome quite easily by declaring \proc{TU} and \proc{UT} as classes with a static function. +In contrast to function templates, class templates can be partially specialized. +By using static member functions, the actual class objects do not need to be created and all the function calls will be resolved at compile time\footnote{One limitation is the maximal template recursion depth. The default value on gcc compiler is 900 and can be increased if required. Typical template nesting depths for coefficient pairs should lie significantly below 100.}. +All templates, no loops, and the resulting sequential code can then be further optimized by the compiler. +\par\medskip +% +Some minimal code snippets explaining the working principle of \proc{IntSplit} are shown in Figure~\ref{fig:intsplit}. +% +In order to prevent the inappropriate use of the \proc{IntSplit} function through incompatible coefficient lists $w_T$ and $w_U$, their sizes are checked in the code at compile time via static assertions. +\par\medskip +% +The fact that $\opW$ is significantly more costly to compute than its counterparts $\opT$ and $\opU$ also justifies the approach of using templates that are hard-coded to call functions \proc{StepT} and \proc{StepU}, as these are the only logical candidates for \proc{IntSplit}. + +\begin{figure}[ht] + \begin{minipage}[c]{\textwidth} + \begin{minipage}[l][10cm][t]{.4\textwidth} + \begin{minipage}[t]{\textwidth} +\begin{lstlisting} +// Specialized Template +// N_T == S_T +// N_U == S_U +template +struct TU { + static void split( + array& w_T, + array& w_U, + float dt) + { /* terminate */ } +}; +\end{lstlisting} + \end{minipage} \\ + \vspace{5mm} + \begin{minipage}[t]{\textwidth} +\begin{lstlisting} +template +struct TU { + static void split( + array& w_T, + array& w_U, + float dt) + /* T */ + stepT(w_T(N_T)*dt); + /* UT... */ + UT + ::split(w_U,w_T,dt); + } +}; +\end{lstlisting} + \end{minipage} + \end{minipage} + \begin{minipage}[l][10cm][t]{.18\textwidth} + \vspace{25mm} + \resizebox{\textwidth}{!}{ +\begin{tikzpicture}[->,>=stealth',auto,node distance=3cm, + thick,main node/.style={circle,draw,font=\sffamily\Large\bfseries}] + \node (TU) at (0,0) {}; + \node (UT) at (2,0) {}; + \node (TUexit) at (0,2) {}; + \node (UTexit) at (2,2) {}; + \node (TUend) at (0,4) {}; + \node (UTend) at (2,4) {}; + \path[every node/.style={font=\sffamily\small}] + (TU) edge[bend right] node {} (UT) + (UT) edge[bend right] node {} (TU) + (TUexit) edge node {} (UTend) + (UTexit) edge node {} (TUend); +\end{tikzpicture} + } + \end{minipage} + \begin{minipage}[r][10cm][t]{.4\textwidth} + \begin{minipage}[t]{\textwidth} +\begin{lstlisting} +// Specialized Template +// N_U == S_U +// N_T == S_T +template +struct UT { + static void split( + array& w_U, + array& w_T, + float dt) + { /* terminate */ } +}; +\end{lstlisting} + \end{minipage} \\ + \vspace{5mm} + \begin{minipage}[t]{\textwidth} +\begin{lstlisting} +template +struct UT { + static void split( + array& w_U, + array& w_T, + float dt) + /* U */ + stepU(w_U(N_U)*dt); + /* TU... */ + TU + ::split(w_T,w_U,dt); + } +}; +\end{lstlisting} + \end{minipage} + \end{minipage} + \end{minipage} + \caption{Simplified C++ snippets for outlining the template invocation sequence in the \proc{IntSplit} function. The template parameters \proc{S\_T} and \proc{S\_U} denote the total size of the coefficient arrays \proc{w\_T} and \proc{w\_U} while the \proc{N\_T} and \proc{N\_U} are the indices of the current coefficient to be applied. The \proc{TU::split} and \proc{UT::split} methods first apply their own step and then invoke each other by gradually incrementing the indices \proc{N\_T} and \proc{N\_U}. \\ + When the conditions \proc{N\_T == S\_T} and \proc{N\_U == S\_U} are satisfied, the partially specialized templates become active and the recursion is terminated.} + \label{fig:intsplit} +\end{figure} diff --git a/LWB/propagators/sections/introduction.tex b/LWB/propagators/sections/introduction.tex new file mode 100644 index 0000000..7531dfb --- /dev/null +++ b/LWB/propagators/sections/introduction.tex @@ -0,0 +1,16 @@ +\section*{Introduction} +% +With the growing understanding for quantum mechanics and dynamics, one can attempt to solve increasingly complex systems and computational problems. +Without doubt, time propagation is at the core of simulating quantum systems and getting a better understanding for them. +\par\medskip +% +The aim of this project was to give a summary of some state of the art quantum time propagators and provide an efficient C++ implementation for evolving Hagedorn wave packets in time. +The produced code is an addition to the C++ WaveBlocks framework whose functionality has already been extended in the context of preceding student projects. +The entire WaveBlocks project, which is partially based on an existing Python project with the same name, can be found at \cite{libwaveblocks}. +\par\medskip +% +Section \ref{sec:operatorsplitting} repeats some of the most important mathematical ideas underlying the time propagation of wave packets and introduces the notations that will be used throughout the rest of the document. +Based on these mathematical observations, section \ref{sec:buildingblocks} prepares a toolbox of essential building blocks that will help in the construction of propagators as well as for the final understanding of the source code. +In section \ref{sec:propagators}, a series of numerical propagators for quantum wave packets are presented and expressed in terms of the components provided earlier. +Some detailed insights into the most important technical ideas of the C++ implementation are given in section \ref{sec:implementation}. +Finally, section \ref{sec:results} presents the results of a series of numerical experiments to show the correctness of the implemented algorithms and analyze some performance aspects like the scaling with step size, dimensionality of the wave packet and order of the splitting coefficients. diff --git a/LWB/propagators/sections/operatorsplitting.tex b/LWB/propagators/sections/operatorsplitting.tex new file mode 100644 index 0000000..d262305 --- /dev/null +++ b/LWB/propagators/sections/operatorsplitting.tex @@ -0,0 +1,71 @@ +\section{Operator Splitting for Time Propagation} +\label{sec:operatorsplitting} +% +In this report, we will use the variable $D$ for describing the number of dimensions and $N$ for the number of energy levels in the system. +The semiclassical time-dependent Schrödinger equation is of the form +% +\begin{align} + \label{math:tdse} + \im \eps \; \del_t \psi (\bvec{x},t) = \opH (\eps) \; \psi (\bvec{x},t) +\end{align} +% +with Hamiltonian operator $\opH$, wave function $\psi(\bvec{x},t)$ depending on the spatial variables $\bvec{x} = (x_1,\dots,x_D) \in \R^D$ and the time variable $t \in \R$. +The parameter $\eps$, also known as the semiclassical parameter, determines the degree of influence of quantum dynamics, approaching classical dynamics in the limit $\eps \rightarrow 0$ and full quantum mechanics for $\eps = 1$. \\ +Moreover, $\eps$ plays an important role in the stability of numerical schemes. +As underlined in \cite{GH_convsemiclassical}, a small parameter $\eps$ can often impose severe constraints on the step size of splitting methods and significantly increase the error. For example, the error was shown to be proportional to $\eps^{-2}$ for Lie-Trotter splitting, Strang splitting and other related methods, meaning that smaller $\eps$ would lead to a larger error. +\par\medskip +% +The solution to the above Schrödinger equation \ref{math:tdse} can be approximated as a finite linear combination of Hagedorn functions +\begin{align} + \label{math:hagedornwp} + \psi(\bvec{x},t) \approx u(\bvec{x},t) + = e^{\im S(t)/\eps^2} \sum_{k \in \K} c_k(t) \varphi_k^\eps [ \bvec{\Pi} (t) ] (\bvec{x}) +\end{align} +where $\K$ represents a finite multi-index set, $c_k \in \C$ are coefficients that depend on time only, $S(t)$ is a phase factor and $\bvec{q},\bvec{p} \in \R^D$ and $\bmat{Q},\bmat{P} \in \C^{D \times D}$ are some parameters. +For convenience, the notation for the parameter pack $\bvec{\Pi} (t) = (\bvec{q}(t),\bvec{p}(t),\bvec{Q}(t),\bvec{P}(t))$ is introduced. +It is worth noting that $\varphi_k^\eps [ \bvec{\Pi} (t) ]$ depends on time $t$ only implicitly through the time dependency of $\bvec{q}(t),\bvec{p}(t),\bvec{Q}(t),\bvec{P}(t)$. +A more detailed analysis of Hagedorn functions is given in \cite{FGL_semiclassical_dynamics} and a comprehensive overview over the relevant mathematical details can be found in \cite{B_master_thesis}. +\par\medskip +% +The Schrödinger equation \ref{math:tdse} can be solved using the ansatz +\begin{align} + \psi (\bvec{x},t+dt) = \exp \left( - \frac{\im}{\eps^2} \opH t \right) \; \psi (\bvec{x},t) \; . +\end{align} +% +However, the matrix exponential is very expensive to compute. +Instead, in the search for a numerical solution, one can exploit the structure of the Hamiltonian operator $\opH$, by splitting it into its kinetic part $\opT$ and potential part $\opV$. +Furthermore, for the rest of the report we will denote by $\opU$ and $\opW$ the quadratic and non-quadratic component of $\opV$ respectively and it will be used that +% +\begin{align} + \opV &= V(\bvec{x}) \\ + \opU &= U(\bvec{q},\bvec{x}) := V(\bvec{q}) + \nabla V(\bvec{q}) (\bvec{x}-\bvec{q}) + + \frac{1}{2} (\bvec{x}-\bvec{q})^T \nabla^2 V(\bvec{q}) (\bvec{x}-\bvec{q}) \\ + \opW &= W(\bvec{q},\bvec{x}) := V(\bvec{x}) - U(\bvec{q},\bvec{x}) \; . +\end{align} +% +In other words, $\opU = U(\bvec{q},\bvec{x})$ is the second order Taylor expansion of the potential energy $\opV = V(\bvec{x})$ around $\bvec{q}$ (the coordinate $\bvec{q}$ will be omitted for simplicity of notation) and $\opW = W(\bvec{q},\bvec{x})$ is the corresponding remainder. +\par\medskip +% +By further recalling the form of the kinetic operator $\opT$ (for a particle of mass $m=1$) +\begin{align} + \opT &= - \sum_{j=1}^D \frac{\eps^4}{2} \frac{\del^2}{\del x_j^2} +\end{align} +% +one can then split the Hamiltonian operator $\opH$ into the components +% +\begin{align} + \opH = \opT + \opV = \opT + \opU + \opW \; . +\end{align} +% +This splitting is so central to the construction of the time propagators in the following sections that it is worthwhile repeating the most important findings from \cite{FGL_semiclassical_dynamics}: +\begin{itemize} + \item The free particle Schrödinger equation ($\opV=0$) can be solved exactly and the wave packet remains in Hagedorn wave packet form (equation~\ref{math:hagedornwp}). + For time propagation, only the parameters $\bvec{\Pi},S$ need to be updated, the coefficients $\{c_k\}_{k \in \K}$ remain unchanged. + \item The Schrödinger equation~\ref{math:tdse} can be solved exactly in a pure quadratic potential, i.e. in a potential $\opV=U(\bvec{x})$ with $\opT=0$. + Again, time propagation only affects the parameters $\bvec{\Pi},S$, not the coefficients $\{c_k\}_{k \in \K}$. + \item In the case of an arbitrary potential $\opV = W(\bvec{x})$ that is not quadratic, a variational approach can be used. + A set of Galerkin functions can be propagated by adapting the coefficients $\{c_k\}_{k \in \K}$ without changing the parameters $\bvec{\Pi},S$. +\end{itemize} +% +The following sections will heavily make use of these properties in order to build composition methods for time integration of quantum wave packets. +A more in depth analysis of composition methods can be found in \cite{GeoNumInt}. diff --git a/LWB/propagators/sections/propagators.tex b/LWB/propagators/sections/propagators.tex new file mode 100644 index 0000000..3acb32d --- /dev/null +++ b/LWB/propagators/sections/propagators.tex @@ -0,0 +1,301 @@ +\section{Propagators} +\label{sec:propagators} + +This section gives a short overview of the numerical properties of the propagators that were implemented in the context of time propagation. +Where appropriate, references are given to a more detailed mathematical analysis of the propagators. +\par\medskip +The methods are called in the sequence determined by the \proc{Evolve} function (\ref{alg:Evolve}) and generally are expressed in terms of the basic building blocks that were presented in the previous section. +In order to specify to which propagator each of the implementations belongs, the C++ inspired notation \proc{PropagatorName.Propagate} is adapted also for pseudo algorithms. +If the methods \proc{PrePropagate} and \proc{PostPropagate} are omitted, it means that they are empty. \\ +Whenever convenient, the direct succession of functions \proc{StepT} and \proc{StepU} was replaced by a call to \proc{IntSplit} which effectively corresponds to a series of $M_{split}$ alternating calls to \proc{StepT} and \proc{StepU} (see explanation in \ref{subsec:tuintsplit}) for better accuracy. + +\subsection{Hagedorn Propagator} +\label{sub:hagedorn_propagator} +% +The Hagedorn propagator, shown in algorithm \ref{alg:hagedorn}, is one of the simplest propagators that can be built by exploiting the numerical properties of the Hamiltonian splitting $\opH = \opT + \opU + \opW$ discussed earlier. +As pointed out in \cite{FGL_semiclassical_dynamics}, this simple time stepping scheme has many beneficial properties like the preservation of the $L^2$ norm of the wave packet, time reversibility and stability in the classical limit $\eps \rightarrow 0$. +Also, in the limit of $\K$ approaching the full basis set, the variational approximation used for the propagation with the non-quadratic part $\opW$ becomes exact. +\begin{algorithm}[ht] + \caption{Single time step with Hagedorn propagator} + \label{alg:hagedorn} + \begin{algorithmic} + \State + \Procedure{Hagedorn.Propagate}{$\Dt$} + \State + \State \Call{StepT}{$\frac{\Dt}{2}$} + \Comment{Step of size $\Dt/2$ with $\opT$} + \State \Call{StepU}{$\Dt$} + \Comment{Step of size $\Dt$ with $\opU$} + \State \Call{StepW}{$\Dt$} + \Comment{Step of size $\Dt$ with $\opW$} + \State \Call{StepT}{$\frac{\Dt}{2}$} + \Comment{Step of size $\Dt/2$ with $\opT$} + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} + + +\subsection{Semiclassical Propagator} +\label{sub:semiclassical_propagator} +% +The central idea of the semiclassical splitting, as introduced in \cite{GH_convsemiclassical}, +is to split the propagation with operators $\opT$ and $\opU$ into many smaller, alternating steps, thereby reducing the dominating error\footnote{for small $\eps$, +the main source of error lies in the updating of $\bvec{\Pi}$ and $S$} in the update of $\bvec{\Pi}$ and $S$. +While there is some additional computational cost caused by a higher number of updates for the parameters $\bvec{\Pi}$ and $S$, the extra effort is usually negligible compared to the propagation with $\opW$ which requires numerical evaluation of multi-dimensional integrals. +\par\medskip +% +In addition, due to the numerical properties of the semiclassical splitting, +it even allows to take larger time steps $\Dt$ than conventional +splitting methods like the YL-splitting (see \cite{GH_convsemiclassical}), while maintaining the same error. +\par\medskip +% +Finally and most importantly, the error is no longer proportional to $1/\eps^2$ but instead +scales linearly in the semiclassical parameter $\eps$, +meaning that a smaller $\eps$ will now reduce the error instead of increasing it. +The error scales with $\eps (\Delta t)^2$ for the semiclassical splitting using the Y-splitting (see \cite{GH_convsemiclassical}), +but the dependency on the time step $\Delta t$ can be improved even further by using different +splittings which effectively corresponds to higher order coefficient pairs $w_T$ and $w_U$. +\par\medskip +% +The steps for the semiclassical propagator are shown in algorithm \ref{alg:semiclassical} and +further details can be found in the original paper \cite{GH_convsemiclassical}. +% +\begin{algorithm}[ht] + \caption{Single time step with Semiclassical propagator} + \label{alg:semiclassical} + \begin{algorithmic} + \State + \Procedure{Semiclassical.Propagate}{$\Dt$} + \State + \State $M_{split} := \lceil 1 + \frac{\sqrt{\Dt}}{\eps^{3/4}} \rceil$ + \Comment{Divide $\Dt$ into smaller steps} + \State + \State \Call{IntSplit}{$\frac{\Dt}{2}, M_{split}, \{ w_T, w_U \}$} + \Comment{$M_{split}$ split steps with $\opT+\opU$} + \State \Call{StepW}{$\Dt$} + \Comment{Single step with $\opW$} + \State \Call{IntSplit}{$\frac{\Dt}{2}, M_{split}, \{ w_T, w_U \}$} + \Comment{$M_{split}$ split steps with $\opT+\opU$} + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} + + +\subsection{Magnus Propagator} +\label{sub:magnus_propagator} +% +As was noted by Magnus in \cite{Magnus1954}, the solution to a differential equation of the form +% +\begin{align} + y'(t) = a(t) y(t) \qquad t \ge 0 +\end{align} +% +with initial value $y(0) = y_0$ can be written as +% +\begin{align} + \label{math:magnussolution} + y(t) = e^{\sigma (t)} y_0 +\end{align} +% +where $\sigma (t)$ is an infinite sum of iterated integrals and commutators, also known as the Magnus series. +The idea behind the Magnus approximation was extensively studied using the Baker-Campbell-Hausdorff (BCH) formula and rooted trees techniques, more details can be found in \cite{Blanes2006}, \cite{Blanes2000}, \cite{Iserles1999}. +\par\medskip +% +In order to approximate the solution from equation \ref{math:magnussolution}, one can take only a finite number of terms from the infinite series whereby a truncation error is committed +(additional error sources in the Magnus method are the discretization of integrals and the approximation of matrix exponentials). +% +The Magnus Propagator has several beneficial numerical properties that are described in \cite{Iserles1999}. +In particular, for solutions that evolve within a Lie group, the same holds for the approximate numerical solution calculated through a truncated Magnus series. +It was also shown in the same paper that the Magnus series can compete with - and in fact may even outplay - classical schemes like Runge-Kutta or Gauss-Legendre. +Although the method is not a symplectic scheme in the usual sense, \cite{Iserles1999} has shown that in practical applications it conserves the Hamiltonian energy just as well as symplectic integrators. +\par\medskip +Moreover, the numerical stability and good performance of the Magnus propagator are not limited to problems where the solution evolves within a Lie Group, but also apply to various problems where this is not the case. +Algorithm \ref{alg:magnus} shows the \emph{MG4} method as presented in \cite{Iserles1999} and implemented in C++ in the scope of this project. + +\begin{algorithm}[ht] + \caption{Single time step with Magnus propagator} + \label{alg:magnus} + \begin{algorithmic} + \State + \Procedure{Magnus.Propagate}{$\Dt$} + \State + \State $h_1 = \frac{3-\sqrt{3}}{6} \Dt$ $h_2 = \frac{2\sqrt{3}}{6} \Dt$ + \Comment Gauss-Legendre coefficients on $[0,\Dt]$ + \State $M_{k} = 1 + \lfloor \sqrt{h_{k} \eps^{-3/8}} \rfloor, \quad k=1,2$ + \Comment number of time steps for splitting + \State + \State \Call{IntSplit}{$h_1, M_1, \{w_T, w_U\}$} + \Comment advance till $\frac{3-\sqrt{3}}{6} \Dt$ + \State $\bmat{A}_1 = - \frac{\im}{\eps^2} \cdot$ \Call{BuildF}{\mbox{}} + \Comment temporarily store interaction matrix + \State \Call{IntSplit}{$h_2, M_2, \{w_T, w_U\}$} + \Comment advance till $\frac{3+\sqrt{3}}{6} \Dt$ + \State $\bmat{A}_2 = - \frac{\im}{\eps^2} \cdot$ \Call{BuildF}{\mbox{}} + \Comment temporarily store interaction matrix + \State $\bmat{\Sigma} = \frac{1}{2} \Dt (\bmat{A}_1 + \bmat{A}_2) + \frac{\sqrt{3}}{12} (\Dt)^2 (\bmat{A}_2 \cdot \bmat{A}_1 - \bmat{A}_1 \cdot \bmat{A}_2)$ + \Comment compute $\sigma (t)$ + \State $\bvec{c} = \exp \left( \bmat{\Sigma} \right) \bvec{c}$ + \Comment update coefficients + \State \Call{IntSplit}{$h_1, M_1, \{w_T, w_U\}$} + \Comment advance till $\frac{6}{6} \Dt$ + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} + +\subsection{McLachlan Propagators} +\label{sub:mcl_propagator} +% +McLachlan (see \cite{McLachlan1995}) has investigated various symplectic schemes for computing the effect of operators of the form $\opX = \opA + \epsilon \opB$ where $\opA$ describes an exactly solvable system and $\opB$ is a perturbation. +The factor $\epsilon$ is a small parameter that indicates the limited influence of $\opB$, not to be confused with the semiclassical parameter $\eps$. +\par\medskip +The \emph{McL} integration schemes are characterized by pairs of weighing coefficients $a_i$ and $b_i$ that allow to represent the exponential $e^{\opX}$ as a product +% +\begin{align} + e^\opX = \prod_i e^{b_i t \epsilon \opB} e^{a_i t \opA} \; . +\end{align} +% +McLachlan also goes through the process of deriving the coefficients for the following schemes +\begin{align*} + &\text{\emph{McL(2s,2)}} &&\text{with error of order } \epsilon (\Dt)^{2s} + \epsilon^2 (\Dt)^2 \\ + &\text{\emph{McL(6,4)}} &&\text{with error of order } \epsilon (\Dt)^{6} + \epsilon^2 (\Dt)^4 \\ + &\text{\emph{McL(8,4)}} &&\text{with error of order } \epsilon (\Dt)^{8} + \epsilon^2 (\Dt)^4 \; . +\end{align*} +% +Note that the error depends on two small parameters, the time step $\Dt$ and the parameter $\epsilon$ that quantifies the influence of $\opB$. +\par\medskip +% +In our case, $\opX = \opH = \opT + \opU + \opW$, the computationally expensive operator $\opW$ takes the role of $\epsilon \opB$, while the remaining Hamiltonian $\opT + \opU$ corresponds to $\opA$. +Our preferred method will therefore be of the form $\opA \opB \opA$ because such a scheme minimizes the number of steps with $\opB = \opW$. +\par\medskip +% +The source code of the project contains the C++ implementation of the propagators \emph{McL(4,2)} and \emph{McL(8,4)} that require only two respectively five evaluations of $\opW$ per step. + +\begin{algorithm}[ht] + \caption{Single time step with McL42 propagator} + \label{alg:mcl} + \begin{algorithmic} + \State + \Procedure{McL42.Propagate}{$\Dt$} + \State + \State $M = \lceil \left( \Dt^4 / \eps^3 \right)^{\frac{1}{8}} \rceil$ + \Comment compute number of sub-steps + \State \Call{IntSplit}{$A_0 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State \Call{StepW}{$B_0 \Dt$} + \Comment $\opW = B$ + \State \Call{IntSplit}{$A_1 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State \Call{StepW}{$B_1 \Dt$} + \Comment $\opW = B$ + \State \Call{IntSplit}{$A_2 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State + \EndProcedure + \\\hrulefill + \State + \Procedure{McL84.Propagate}{$\Dt$} + \State + \State $M = \lceil \left( \Dt^4 / \eps^3 \right)^{\frac{1}{8}} \rceil$ + \Comment compute number of sub-steps + \State \Call{IntSplit}{$A_0 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State \Call{StepW}{$B_0 \Dt$} + \Comment $\opW = B$ + \State \Call{IntSplit}{$A_1 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State \Call{StepW}{$B_1 \Dt$} + \Comment $\opW = B$ + \State \Call{IntSplit}{$A_2 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State \Call{StepW}{$B_2 \Dt$} + \Comment $\opW = B$ + \State \Call{IntSplit}{$A_3 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State \Call{StepW}{$B_3 \Dt$} + \Comment $\opW = B$ + \State \Call{IntSplit}{$A_4 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State \Call{StepW}{$B_4 \Dt$} + \Comment $\opW = B$ + \State \Call{IntSplit}{$A_5 \Dt, M, \{w_T, w_U\}$} + \Comment $\opT + \opU = A$ + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} + + +\subsection{Processing Propagators} +\label{sub:pre764_propagator} +% +The idea of processing propagators is to carry out the time evolution with a Hamiltonian that is slightly perturbed from the original one, which can be achieved by applying a pre and post processing step. +The exponential of the Hamiltonian can be re-formulated as +% +\begin{align} + e^{- \Dt \opH (\Dt)} = e^P e^{- \Dt K} e^{-P} +\end{align} +% +where the pre and post processing steps are applied via the multiplication with matrices $e^P$ and $e^{-P}$ (also referred to as \emph{processors}). +While the processors only need to be applied once at the very beginning and end of the propagation (in order to return to the original Hamiltonian), the multiplication with the \emph{kernel} $e^{-\Dt K}$ is repeated $M_{tot}$ times, once for every time step. \\ +As a consequence, one wants to choose the processor $e^P$ in such a way that the evaluation of the kernel $e^{-\Dt K}$ is as simple as possible in order to minimize the work associated with the computation of the kernel. +\par\medskip +% +Unlike most commonly used integration schemes, the family of processing propagators are symplectic integrators. +They are not suited for every kind of Hamiltonian, and the work of Blanes, Casas and Ros in \cite{Blanes1999} gives a method for finding the necessary conditions that need to be satisfied in order for processing methods to be applicable. +\par\medskip +% +Algorithm \ref{alg:pre764} presents the \emph{Pre764} processing method as derived in \cite{Blanes1999}. +The coefficient arrays $\alpha$ and $\beta$ (both of length $k=4$) are used for propagating the wave packet with \proc{StepW} and \proc{IntSplit} in the \proc{Propagate} method, while the coefficient arrays $Y$ and $Z$ (both of length $v=6$) are used in \proc{PrePropagate} and \proc{PostPropagate} to transform the wave packet with the processor. +% +\begin{algorithm}[ht] + \caption{Single time step with Pre764 propagator} + \label{alg:pre764} + \begin{algorithmic} + \State + \Procedure{Pre764.PrePropagate}{$\Dt$} + \State + \State $M_{pre} = 1+ \left\lfloor \sqrt{\Dt \eps^{-\frac{3}{4}}} \right\rfloor$ + \Comment compute number of time steps + \For{$j=1,...,v$} \Comment $v=6$ for Pre(7,6,4) + \State \Call{IntSplit}{$-Z_j \Dt, M_{pre}, \{w_T, w_U\}$} + \Comment $M_{pre}$ alternating steps with $\opT$ and $\opU$ + \State \Call{StepW}{$\upic, -Y_j \Dt$} + \Comment single step with $\opW$ + \EndFor + \State + \EndProcedure + \\\hrulefill + \State + \Procedure{Pre764.Propagate}{$\Dt$} + \State + \State $M = 1+ \left\lfloor \sqrt{\Dt \eps^{-\frac{3}{4}}} \right\rfloor$ + \Comment compute number of time steps + \For{$j=1,...,k$} \Comment $k=4$ for Pre(7,6,4) + \State \Call{StepW}{$\alpha_j \Dt$} + \Comment single step with $\opW$ + \State \Call{IntSplit}{$\beta_j \Dt, M, \{w_T, w_U\}$} + \Comment $M$ alternating steps with $\opT$ and $\opU$ + \EndFor + \State + \EndProcedure + \\\hrulefill + \State + \Procedure{Pre764.PostPropagate}{$\Dt$} + \State + \State $M_{post} = 1+ \left\lfloor \sqrt{\Dt \eps^{-\frac{3}{4}}} \right\rfloor$ + \Comment compute number of time steps + \For{$j=v,...,1$} \Comment $v=6$ for Pre(7,6,4) + \State \Call{StepW}{$Y_j \Dt$} + \Comment single step with $\opW$ + \State \Call{IntSplit}{$Z_j \Dt, M_{post}, \{w_T, w_U\}$} + \Comment $M_{post}$ alternating steps with $\opT$ and $\opU$ + \EndFor + \State + \EndProcedure + \end{algorithmic} +\end{algorithm} + + diff --git a/LWB/propagators/sections/results.tex b/LWB/propagators/sections/results.tex new file mode 100644 index 0000000..1ae1a3c --- /dev/null +++ b/LWB/propagators/sections/results.tex @@ -0,0 +1,223 @@ +\section{Results} +\label{sec:results} +% +The newly developed code was used to carry out some numerical experiments which are presented in this section. +Aspects that need to be checked are the correctness of the propagators from a quantum mechanical point of view (subsection~\ref{subsec:physical}), convergence of the methods (subsection~\ref{subsec:convergence}) and benchmark of the compute time (subsection~\ref{subsec:benchmark}). +\par\medskip +% +The number of steps $M_{split}$ for the \proc{IntSplit} method can be calculated based on an analysis of the numerical properties of the methods, such that the stability of the propagators is guaranteed. +However, it was found that such an adaptive number of splitting intervals (usually depending on $\Dt$ itself) made the numerical experiments hardly comparable (both in terms of reunite and precision). +Therefore, for all the presented results, $M_{split}=4$ was used for every call to \proc{IntSplit} in order to create fair conditions for comparing the different propagators. +\par\medskip +% +The potentials used for these numerical experiments are the harmonic, torsional and Morse potential that are now briefly introduced. +The semiclassical parameter $\eps = 0.01$ was used if not specified otherwise. + + +\paragraph{Harmonic Potential} in $D$ dimensions +\begin{align} + \label{math:harmonic} + V(\bvec{x}) = \frac{1}{2} \sum_{i=1}^D x_i^2 +\end{align} +% +Initial values: $\bvec{q} = (1,\dots,1)^T$, $\bvec{q} = (1,\dots,1,-1,\dots,-1)$, $\bmat{Q}$ and $\bmat{P} = \im \bmat{Q}^{-1}$. \\ +\par\medskip + + +\paragraph{Torsional Potential} +% +\begin{align} + \label{math:torsional} + V(\bvec{x}) = \sum_{i=1}^D \left( 1 - \cos x_i \right) +\end{align} +% +Initial values: $\bmat{Q} = \bmat{I}_N$, $\bmat{P} = \im \bmat{I}_N$, $\bvec{q} = (1,0,\dots,0)$, $\bvec{p} = (0,\dots,0)$ (see \cite{FGL_semiclassical_dynamics}). + +\paragraph{Morse Potential} +% +\begin{align} + \label{math:morse} + V(\bvec{x}) = C e^{-2a(x-x_0)} - 2 e^{-a(x-x_0)} \; . +\end{align} +% +Constants (see \cite{Unpublished}): $C = 0.004164$, $a=0.896696$, $x_0 = 5.542567$. \\ +Initial values: $q_0 = 6.0426$, $p_0 = -0.1100$, $Q_0 = 3.4957$, $P_0 = 0.2861 \im$. +The Morse potential was used with $\eps = 0.0484$ for all simulations. + + +\subsection{Physical Correctness} +\label{subsec:physical} +% +In order to analyze the physical correctness of the implemented propagators, a wave packet was propagated in the various potentials introduced above and the energy conservation was observed. +\par\medskip +% +Figure \ref{fig:energy_Semiclassical} shows the time evolution of energy and the energy drift for the semiclassical propagator. The corresponding energy plots for the remaining propagators can be found in the appendix to the report. +The equivalent conservation plot for the two dimensional case is shown in figure \ref{fig:energy_Semiclassical_2D}. +As is clearly shown by the graphs, the propagators conserve the total energy (except for very small oscillations). +% +\begin{figure}[ht] + \centering + \begin{minipage}[c]{\textwidth} + \begin{center} + \large Semiclassical Propagator \\[1mm] + \normalsize Energy Evolution and Drift + \vspace{4mm} + \end{center} + \end{minipage} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_Semiclassical_energies.png} + \includegraphics[width=.45\textwidth]{figures/harmonic_1D_Semiclassical_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/torsional_1D_Semiclassical_energies.png} + \includegraphics[width=.45\textwidth]{figures/torsional_1D_Semiclassical_drift.png} \\ + \includegraphics[width=.45\textwidth]{figures/morse_1D_Semiclassical_energies.png} + \includegraphics[width=.45\textwidth]{figures/morse_1D_Semiclassical_drift.png} + \caption{Energy evolution and drift for a 1D wave packet propagated with the Semiclassical propagator in a harmonic potential (top), a torsional potential (middle) and a Morse potential (bottom). + (Parameters: $N=1$, $D=1$ $|\K|=16$ $\eps=0.01$ (Morse $0.0484$), $T=10$ (Morse $T=50$), $\Dt=0.001$ (Morse $\Dt=0.005$), Semiclassical propagator with \emph{Y4} splitting for \proc{IntSplit})} + \label{fig:energy_Semiclassical} +\end{figure} +% +\begin{figure}[ht] + \centering + \includegraphics[width=.8\textwidth]{figures/torsional_2D_Semiclassical_energies.png} + \includegraphics[width=.8\textwidth]{figures/torsional_2D_Semiclassical_drift.png} + \caption{Energy evolution and drift for a wave packet propagated with the Semiclassical propagator in a 2D torsional potential. + (Parameters: $N=1$, $D=2$ $|\K|=16$ $\eps=0.01$, $T=4$, $\Dt=0.005$, Semiclassical propagator with \emph{KL10} splitting for \proc{IntSplit})} + \label{fig:energy_Semiclassical_2D} +\end{figure} + + +\subsection{Effect of step size} +\label{subsec:convergence} +% +A convergence analysis was carried out in which the error of each propagator was recorded while reducing the step size $\Dt$. +A \emph{McL84} propagator with \emph{KL10} splitting coefficients and step size $\Dt = 0.0001$ (Morse: $\Dt = 0.0005$) was used to propagate the wave packet to a final time $T = 10$ (Morse: $T=50$) and create a reference solution. +\par\medskip +% +The error between two wave packets was computed by evaluating them on a grid on $[-5,5]$ (Morse: $[0,10]$) with 1000 grid points and taking the $L_2$ norm of the differences. +Figure \ref{fig:error_analysis} shows how the error changes for different step sizes $\Dt$. +All propagators converge towards the reference solution, quickly reaching the limit imposed by the machine precision. +The simplistic Hagedorn propagator shows the largest error and a convergence order of two, similarly to the \emph{Semiclassical} and \emph{McL42} propagators. +The more advanced propagators \emph{MG4}, \emph{McL84} and \emph{Pre764} even show a convergence order of four. +Interestingly enough, if only a \emph{LT} splitting is used in the \proc{IntSplit} method, then the simplistic \emph{Hagedorn} propagator outperforms all the others. However, changing the splitting scheme to \emph{Y4} again gives the expected behavior. +% +\begin{figure}[ht] + \centering + \includegraphics[width=.2\textwidth]{figures/convergence_legend.png} \\ + \includegraphics[width=.45\textwidth]{figures/convergence_harmonic_1D_lt.png} + \includegraphics[width=.45\textwidth]{figures/convergence_harmonic_1D_y4.png} \\ + \includegraphics[width=.45\textwidth]{figures/convergence_torsional_1D_lt.png} + \includegraphics[width=.45\textwidth]{figures/convergence_torsional_1D_y4.png} \\ + \includegraphics[width=.45\textwidth]{figures/convergence_morse_1D_lt.png} + \includegraphics[width=.45\textwidth]{figures/convergence_morse_1D_y4.png} + \caption{Convergence of the various propagators towards the reference solution for \emph{LT} splitting (left) and \emph{Y4} splitting (right). The potentials are the harmonic potential (top), the torsional potential (middle) and the Morse potential (bottom). The $L_2$ norm was measured by projecting the wave function on a grid with $1000$ nodes in the range $[-5,5]$ (Morse $[0,10]$) and taking the differences of the current solution and the reference solution at $T=10$ (Morse $T=50$).} + \label{fig:error_analysis} +\end{figure} + +\subsection{Benchmark} +\label{subsec:benchmark} +% +In order to benchmark the code, three timing experiments were carried out: a simple comparison of run-times for different propagators, an investigation of the computational cost associated with the usage of high order splitting schemes, and an analysis of the compute time as a function of the dimension $D$. +In all cases, the $1D$ torsional potential introduced at the beginning of this section was used. +\par\medskip +% +The run-time analysis was carried out on a Linux (Kernel 4.8.7, Fedora 24 Workstation) Quad-Core machine with an Intel Core i5-3210M Processor and +4GB of RAM. The C++ code was compiled with the GNU compiler, version 6.2.1, and using the \emph{-Ofast} optimization flag. + + + +\subsubsection{Comparison of propagators} +% +A quick comparison of timings for propagating a wave packet in a 1D torsional potential with different propagators is given in Table~\ref{tab:comparison}. +Surprisingly, the Semiclassical propagator seems to take the same time as the Hagedorn propagator, despite its more complex structure. +The Magnus Propagator \emph{MG4} is cheaper than expected, as it involves two evaluations of $\matrixel{\varphi_k}{W}{\varphi_l}$ per time step but is only slightly slower than the Semiclassical operator with one evaluation of the inner product. +The slowest propagator is \emph{McL84} which involves five propagations with $\opW$ per step. +% +\begin{table}[ht] + \centering + \begin{tabular}{|l | r |} + \hline + \multicolumn{1}{|c}{\textbf{Propagator}} & + \multicolumn{1}{|c|}{\textbf{Timing [s]}} \\ + \hline + Hagedorn & 27.83 \\ + Semiclassical & 27.66 \\ + MG4 & 37.85 \\ + McL42 & 49.20 \\ + McL84 & 107.62 \\ + Pre764 & 70.93 \\ + \hline + \end{tabular} + \caption{Run times for different propagators. The standard deviation of these measurements was below $1\%$. + (Parameters: torsional potential, $N=1$, $D=1$, $|\K|=16$ $\eps=0.01$, $T=10$, $\Dt=0.0001$, Semiclassical propagator with \emph{Y4} splitting for \proc{IntSplit}) } + \label{tab:comparison} +\end{table} + + +\subsubsection{Splitting Schemes} +% +As mentioned previously, the splitting coefficients $\{ w_T, w_U \}$ that are used as weights on the time step $\dt$ in the \proc{IntSplit} method can have vastly different complexity, ranging from the \emph{Lie-Trotter} coefficients (one coefficient for propagation with $\opT$, one for propagation with $\opU$) up to the \emph{KL10} coefficients (34 coefficients for each operator). \\ +Higher order schemes are usually preferred in terms of accuracy, but they come at the price of longer computation time. +\par\medskip +% +Therefore, a numerical experiment was carried out in order to analyze how much computational time is required for using different sizes of coefficient pairs $\{ w_T, w_U \}$ in the \proc{IntSplit} method. +The same simulations were carried out for Python and C++ code. +The results are shown in figure \ref{fig:benchmarksplit}. +\par\medskip +% +The Python simulations were run using the Pade matrix exponential, Gauss-Hermite Quadrature, and a \proc{HyperbolicCutShape} as basis shape, while +C++ used the more expensive \proc{HyperCubicShape} as in all other benchmarks. +More information about the different basis shapes can be found in \cite{B_master_thesis} +\par\medskip +% +Looking at the Python timings only, the measurements suggest that in the case of the \emph{KL10} coefficient set, at least 80\% of the total run-time is spend in the \proc{IntSplit} function (since larger coefficient sets $\{ w_T, w_U \}$ only affect the number of steps with operators $\opT$ and $\opU$, but not with operator $\opW$). +\par\medskip +% +When bringing the timings in context with the C++ code, a comparison of absolute run-times is of course not fair since C++ is intrinsically faster than a scripting language like Python and the C++ version has been optimized for speed in many different ways. +However, a very significant result in the context of the work on Propagators is how the speedup increases for larger coefficient pairs. +Using the \emph{KL10} coefficient set instead of a simple \emph{LT} set, the Python code takes more than five times longer. +The same comparison on C++ code shows that in the case of the largest coefficient set \emph{KL10} it doesn't even take twice as long as for the basic \emph{LT} version. +This is a very encouraging result as it means that high order splitting coefficients come at almost no extra cost and the \proc{IntSplit} method was implemented efficiently. +\par\medskip +% +\begin{figure}[ht] + \centering + \begin{tabular}{|l | r | r | r | r |} + \hline + \multicolumn{1}{|c}{\textbf{Splitting}} & + \multicolumn{1}{|c}{\textbf{No. of coefs}} & + \multicolumn{1}{|c}{\textbf{Python timing [s]}} & + \multicolumn{1}{|c}{\textbf{C++ timing [s]}} & + \multicolumn{1}{|c|}{\textbf{Speedup}} \\ + \hline + LT &1 &47.155 &1.3905 &\textbf{33.91} \\ + S2 &2 &52.934 &1.3817 &\textbf{38.31} \\ + Y4 &4 &65.351 &1.4386 &\textbf{45.42} \\ + PRKS6 &7 &84.075 &1.5122 &\textbf{55.59} \\ + Y61 &8 &90.260 &1.4568 &\textbf{61.96} \\ + KL6 &10 &102.531 &1.6049 &\textbf{63.89} \\ + BM63 &15 &133.415 &1.7763 &\textbf{75.11} \\ + KL8 &18 &151.438 &1.8147 &\textbf{83.45} \\ + KL10 &34 &249.749 &2.2168 &\textbf{112.66} \\ + \hline + \end{tabular} + \begin{center} + \includegraphics[width=.8\textwidth]{figures/coefs.png} + \end{center} + \caption{Comparison of computation times for different splitting coefficients $\{ w_T, w_U \}$ for Python code vs. C++ code. The absolute timings are shown in the top graph, the speedup on the bottom. It is remarkable that the initial speedup factor is about 30, but grows for larger coefficient sets. The speedup factor for the largest tested coefficient pair \emph{KL10} is over 110. + All the timings in the table were measured by taking the average of 10 independent runs. The standard deviation of these measurements was below $1\%$ for the Python timings and below $0.1\%$ for the C++ timings. (Parameters: $1D$ torsional potential, $N=1$, $|\K|=16$ $\eps=0.01$, $T=10$, $\Dt=0.001$, Semiclassical propagator)} + \label{fig:benchmarksplit} +\end{figure} + + +\subsubsection{Scaling with dimensionality D} +% +The question that was addressed in this benchmark is how the compute time of the code scales with increasing dimension $D$ of the wave packet. +Figure \ref{fig:dimension_analysis} shows a clear exponential scaling with the dimension $D$. +With a run-time of just about 35 minutes for 100 time steps with the Semiclassical Propagator, basis size of four and \emph{LT} splitting in $D=5$ dimensions, the computation is still quite feasible. +% +\begin{figure}[ht] + \centering + \includegraphics[width=.8\textwidth]{figures/dim.png} + \caption{Scaling of the computation time with wave packet dimension $D$. + (Parameters: $D$-dimensional torsional potential, $N=1$, $|\K|=4$ $\eps=0.01$, $T=1$, $\Dt=0.01$, Semiclassical propagator with \emph{LT} splitting for \proc{IntSplit})} + \label{fig:dimension_analysis} +\end{figure} diff --git a/LWB/propagators/sections/title.tex b/LWB/propagators/sections/title.tex new file mode 100644 index 0000000..88506ac --- /dev/null +++ b/LWB/propagators/sections/title.tex @@ -0,0 +1,52 @@ +\begin{titlepage} + \begin{center} + + \hfill + + \vspace{2cm} + + {\Large \textsc{Time Propagation\\ of Hagedorn Wave Packets\\[5pt]}} + {\large \textsc{an efficient implementation in C++}} + + \vspace{1cm} + + {\large{Bachelor Thesis}} + + \vspace{1cm} + + {\emph{written by}} \\ + { + Matthias Untergassmair + } + + \vspace{1cm} + + {\emph{advised by}} \\ + { + Raoul Bourquin + } + + \vspace{1cm} + + {\emph{supervised by}} \\ + { + Dr. Vasile Gr\u{a}dinaru\\ + Prof. Dr. Ralf Hiptmair + } + + \vfill + + Seminar for Applied Mathematics\\ + ETH Zurich + + \vspace{.5cm} + + \emph{{Fall semester 2016}} + + \vspace{.5cm} + + \includegraphics[width=0.5\linewidth]{figures/eth-logo.pdf} + + \end{center} + +\end{titlepage} diff --git a/LWB/propagators/uml.dot b/LWB/propagators/uml.dot new file mode 100644 index 0000000..979d92d --- /dev/null +++ b/LWB/propagators/uml.dot @@ -0,0 +1,100 @@ +digraph Propagators { + + graph [dpi = 300]; + + fontname = "Bitstream Vera Sans" + fontsize = 48 + + node [ + fontname = "Bitstream Vera Sans" + fontsize = 8 + shape = "record" + ] + + edge [ + fontname = "Bitstream Vera Sans" + fontsize = 8 + ] + + PropagatorBase [ + label = "{Base Propagator|+ evolve(T,Dt,callback) \l|# virtual getName() \l# virtual propagate(Dt) \l# virtual pre_propagate(Dt) \l# virtual post_propagate(Dt) \l|# stepT(Dt) \l# stepU(Dt) \l# stepW(Dt) \l# buildF() \l# intSplit(Dt,M,\{w\}) \l- TU::split(Dt,\{w\}) \l- UT::split(Dt,\{w\}) \l|# Packet_t wpacket \l# Potential_t V \l# Matrix F \l}" + width = "2" + ] + + PropagatorDerived [ + label = "{Derived Propagator|# override getName() \l# override propagate(Dt) \l# override pre_propagate(Dt) \l# override post_propagate(Dt) \l}" + width = "2" + ] + + edge [ + arrowtail = "empty" + ] + + PropagatorBase -> PropagatorDerived [dir="back"] + + /* + subgraph clusterVectorialWavepackets { + label = "Vectorial Wavepackets" + + HomogeneousHaWp__Component [ + label = "{Component|- shape \l- coefficients \l}" + group = groupHomogeneousHaWp + ] + + HomogeneousHaWp [ + label = "{HomogeneousHaWp|- components \l- eps \l- parameters \l+ evaluate(grid) \l}" + group = groupHomogeneousHaWp + ] + + + InhomogeneousHaWp__Component [ + label = "{Component|- parameters \l- shape\l- coefficients \l}" + group = groupInhomogeneousHaWp + ] + + InhomogeneousHaWp [ + label = "{InhomogeneousHaWp|- components \l- eps \l+ evaluate(grid) \l}" + group = groupInhomogeneousHaWp + ] + + + HaWpGradient__Component [ + label = "{Component|- coefficients \l}" + group = groupHaWpGradient + ] + + HaWpGradient [ + label = "{HaWpGradient|- components \l- eps \l- parameters \l- shape \l+ evaluate(grid) \l}" + group = groupHaWpGradient + ] + + edge [ + arrowtail = "none" + headlabel = "N" + taillabel = "1" + ] + + HomogeneousHaWp -> HomogeneousHaWp__Component [dir="back"] + InhomogeneousHaWp -> InhomogeneousHaWp__Component [dir="back"] + + edge [ + arrowtail = "none" + headlabel = "D" + taillabel = "1" + ] + + HaWpGradient -> HaWpGradient__Component [dir="back"] + } + */ + + /* + edge [ + arrowtail = "empty" + ] + + AbstractScalarHaWpBasis -> HaWpGradient [dir="back"] + AbstractScalarHaWp -> HomogeneousHaWp__Component [dir="back"] + AbstractScalarHaWp -> InhomogeneousHaWp__Component [dir="back"] + AbstractScalarHaWp -> HaWpGradient__Component [dir="back"] + */ +} diff --git a/LWB/propagators/uml.png b/LWB/propagators/uml.png new file mode 100644 index 0000000..ef02792 Binary files /dev/null and b/LWB/propagators/uml.png differ diff --git a/LWB/serialization/data_layout.pdf b/LWB/serialization/data_layout.pdf new file mode 100644 index 0000000..2421d87 Binary files /dev/null and b/LWB/serialization/data_layout.pdf differ diff --git a/LWB/serialization/selection.pdf b/LWB/serialization/selection.pdf new file mode 100644 index 0000000..ae02b91 Binary files /dev/null and b/LWB/serialization/selection.pdf differ diff --git a/LWB/serialization/thesis.pdf b/LWB/serialization/thesis.pdf new file mode 100644 index 0000000..2d6cbd6 Binary files /dev/null and b/LWB/serialization/thesis.pdf differ diff --git a/LWB/serialization/writing_dataset.pdf b/LWB/serialization/writing_dataset.pdf new file mode 100644 index 0000000..632a461 Binary files /dev/null and b/LWB/serialization/writing_dataset.pdf differ diff --git a/references/propagators.bib b/references/propagators.bib new file mode 100644 index 0000000..c0a4a53 --- /dev/null +++ b/references/propagators.bib @@ -0,0 +1,130 @@ +@article{Iserles1999, + author="A. Iserles and A. Marthinsen and S. P. N{\o}rsett", + title="On the Implementation of the Method of Magnus Series for Linear Differential Equations", + journal="BIT Numerical Mathematics", + year="1999", + volume="39", + number="2", + pages="281--304", + issn="1572-9125", + doi="10.1023/A:1022393913721", + url="http://dx.doi.org/10.1023/A:1022393913721" +} + +@article{Magnus1954, + author = {Wilhelm Magnus}, + title = {On the exponential solution of differential equations for a linear operator}, + journal = {Communications on Pure and Applied Mathematics}, + volume = {7}, + number = {4}, + publisher = {Wiley Subscription Services, Inc., A Wiley Company}, + issn = {1097-0312}, + url = {http://dx.doi.org/10.1002/cpa.3160070404}, + doi = {10.1002/cpa.3160070404}, + pages = {649--673}, + year = {1954} +} + +@article{Blanes2006, + title = "Fourth- and sixth-order commutator-free Magnus integrators for linear and non-linear dynamical systems ", + journal = "Applied Numerical Mathematics ", + volume = "56", + number = "12", + pages = "1519 - 1537", + year = "2006", + note = "", + issn = "0168-9274", + doi = "http://dx.doi.org/10.1016/j.apnum.2005.11.004", + url = "//www.sciencedirect.com/science/article/pii/S0168927405002163", + author = "S. Blanes and P.C. Moan", + keywords = "Time-dependent differential equations", + keywords = "Initial value problems", + keywords = "Numerical methods", + keywords = "Magnus expansion " +} + +@Article{Blanes2000, + author="S. Blanes and F. Casas and J. Ros", + title="Improved High Order Integrators Based on the Magnus Expansion", + journal="BIT Numerical Mathematics", + year="2000", + volume="40", + number="3", + pages="434--450", + issn="1572-9125", + doi="10.1023/A:1022311628317", + url="http://dx.doi.org/10.1023/A:1022311628317" +} + +@article{Blanes1999, + author = {S. Blanes and F. Casas and J. Ros}, + title = {Symplectic Integration with Processing: A General Study}, + journal = {SIAM Journal on Scientific Computing}, + volume = {21}, + number = {2}, + pages = {711-727}, + year = {1999}, + doi = {10.1137/S1064827598332497}, + URL = { http://dx.doi.org/10.1137/S1064827598332497 }, + eprint = { http://dx.doi.org/10.1137/S1064827598332497 } +} + +@Article{McLachlan1995, + author="Robert I. McLachlan", + title="Composition methods in the presence of small parameters", + journal="BIT Numerical Mathematics", + year="1995", + volume="35", + number="2", + pages="258--268", + issn="1572-9125", + doi="10.1007/BF01737165", + url="http://dx.doi.org/10.1007/BF01737165" +} + + +@book{Gamma1995, + author = {Erich Gamma and Richard Helm and Ralph Johnson and John Vlissides}, + title = {Design Patterns: Elements of Reusable Object-oriented Software}, + year = {1995}, + isbn = {0-201-63361-2}, + publisher = {Addison-Wesley Longman Publishing Co., Inc.}, + address = {Boston, MA, USA}, +} + +@Article{LWB_wavepackets, + author="Lionel Miserez", + title="Porting WaveBlocksND Matrix Potential Functionality to C++", + journal="Libwaveblocks", + year="2015", + url="https://github.com/raoulbq/TeX" +} + + +@Article{LWB_innerproducts, + author="Benedek Vartok", + title="Implementation of WaveBlocksND’s Quadrature and Inner Products in C++", + journal="Libwaveblocks", + year="2015", + url="https://github.com/raoulbq/TeX" +} + + +@book{GeoNumInt, + author = {Ernst Hairer and Christian Lubich and Gerhard Wanner}, + title = {Geometric Numerical Integration}, + year = {2002}, + isbn = {3-540-30663-3}, + publisher = {Springer}, + address = {Berlin, Germany}, +} + +@article{Unpublished, + title = "Raising the Order of Convergence in the Semiclassical Splitting", + author = "Sergio Blanes and Raoul Bourquin and Vasile Gradinaru", + year = {2017}, + isbn = {}, + publisher = {}, + address = {}, + journal = "" +}