|
1 | 1 | /** |
2 | 2 | * @file Schrodinger2D.cpp |
3 | | - * @brief Solves a 2D Time Dependant Schrodinger's Equation using mimetic methods (MOLE). |
| 3 | + * @brief Solves a 2D Time Dependant Schrodinger's Equation using mimetic |
| 4 | + * methods (MOLE). |
4 | 5 | * |
5 | 6 | * Can be visualized as a surface using gnuplot by running the command: |
6 | | - * |
| 7 | + * |
7 | 8 | * ./Schrodinger2D | sed 1d | gnuplot -p -e "splot '<cat'" |
8 | 9 | */ |
9 | 10 |
|
10 | | -#include "mole.h" |
11 | | -#include <iostream> |
| 11 | +#include "mole.h" |
12 | 12 | #include <cmath> |
13 | 13 | #include <iomanip> |
| 14 | +#include <iostream> |
| 15 | + |
| 16 | +using namespace arma; |
| 17 | + |
| 18 | +int main() { |
| 19 | + int p = 105; // Number of time steps for the simulation |
| 20 | + double Lxy = 1.0; |
| 21 | + int k = 2; // Order of accuracy |
| 22 | + int m = 50; // Grid points in x |
| 23 | + int n = 50; // Grid points in y |
| 24 | + int nx = 2; // Energy level in x |
| 25 | + int ny = 2; // Energy level in y |
| 26 | + double dx = Lxy / m; // Step in x |
| 27 | + double dy = Lxy / n; // Step in y |
| 28 | + double dt = dx; // Time step size |
| 29 | + |
| 30 | + // Define staggered grids |
| 31 | + vec xgrid = |
| 32 | + join_vert(vec({0}), linspace(dx / 2, Lxy - dx / 2, m), vec({Lxy})); |
| 33 | + vec ygrid = |
| 34 | + join_vert(vec({0}), linspace(dy / 2, Lxy - dy / 2, n), vec({Lxy})); |
| 35 | + |
| 36 | + // Initialize 2D staggered grid |
| 37 | + mat X, Y; |
| 38 | + Utils utils; |
| 39 | + utils.meshgrid(xgrid, ygrid, X, Y); |
| 40 | + |
| 41 | + // Initialize Laplacian operator with Robin BC |
| 42 | + Laplacian L(k, m, n, dx, dy); |
| 43 | + RobinBC BC(k, m, 1, n, 1, 1, 0); |
| 44 | + L = L + BC; |
| 45 | + |
| 46 | + // Hamiltonian definition |
| 47 | + auto H = [&](const vec &x) { |
| 48 | + vec result = 0.5 * (L * x); |
| 49 | + return result; |
| 50 | + }; |
| 51 | + |
| 52 | + // Define wave numbers |
| 53 | + auto kx = [&](int nx) { return nx * M_PI / Lxy; }; |
| 54 | + auto ky = [&](int ny) { return ny * M_PI / Lxy; }; |
| 55 | + |
| 56 | + double A = 2 / Lxy; |
| 57 | + |
| 58 | + // Initialize the wavefunction psi_old |
| 59 | + mat Psi_grid(m + 2, n + 2, fill::zeros); |
| 60 | + |
| 61 | + // Manually pad Psi_grid |
| 62 | + for (int i = 1; i <= m; i++) { |
| 63 | + for (int j = 1; j <= n; j++) { |
| 64 | + Psi_grid(i, j) = A * sin(kx(nx) * X(i, j)) * sin(ky(ny) * Y(i, j)); |
| 65 | + } |
| 66 | + } |
| 67 | + |
| 68 | + // Convert to column vector for compatibility |
| 69 | + vec psi_old = vectorise(Psi_grid); // Convert to column vector |
| 70 | + // Create interpolators |
| 71 | + Interpol I(m, n, 0.5, 0.5); |
| 72 | + Interpol I2(true, m, n, 0.5, 0.5); |
| 73 | + |
| 74 | + I *= dt; // Scale by dt |
| 75 | + I2 *= 0.5 * dt; |
| 76 | + |
| 77 | + vec v_old(2 * m * n + m + n, fill::zeros); // Initialize v_old |
| 78 | + |
| 79 | + // Initialize Psi_re with zeros |
| 80 | + mat Psi_re = zeros<mat>(m + 2, n + 2); |
| 81 | + |
| 82 | + try { |
| 83 | + // Time-stepping loop (Position Verlet) |
| 84 | + for (int t = 0; t <= p; ++t) { |
| 85 | + // Position Verlet algorithm: Update psi_old based on v_old |
| 86 | + psi_old = psi_old + I2 * v_old; |
| 87 | + |
| 88 | + // Calculate v_new using the Hamiltonian and psi_old |
| 89 | + vec v_new = v_old + I * H(psi_old); |
| 90 | + |
| 91 | + // Update psi_old based on v_new |
| 92 | + psi_old = psi_old + I2 * v_new; |
| 93 | + |
| 94 | + // Update Psi_re for output |
| 95 | + Psi_re = reshape(psi_old, m + 2, n + 2); |
| 96 | + |
| 97 | + if (t == p) { |
| 98 | + std::cout << "Final Time Step " << t << ": X, Y, Psi" << std::endl; |
| 99 | + for (size_t i = 0; i < Psi_re.n_rows; ++i) { |
| 100 | + for (size_t j = 0; j < Psi_re.n_cols; ++j) { |
| 101 | + std::cout << std::fixed << std::setprecision(5) << X(i, j) << ", " |
| 102 | + << Y(i, j) << ", " << Psi_re(i, j) << std::endl; |
| 103 | + } |
| 104 | + } |
| 105 | + } |
| 106 | + |
| 107 | + // Update variables for the next time step |
| 108 | + v_old = v_new; |
| 109 | + } |
| 110 | + } catch (const std::exception &e) { |
| 111 | + std::cerr << "Error: " << e.what() << std::endl; |
| 112 | + return -1; |
| 113 | + } |
14 | 114 |
|
15 | | -using namespace arma; |
16 | | - |
17 | | -int main() { |
18 | | - int p = 105; // Number of time steps for the simulation |
19 | | - double Lxy = 1.0; |
20 | | - int k = 2; // Order of accuracy |
21 | | - int m = 50; // Grid points in x |
22 | | - int n = 50; // Grid points in y |
23 | | - int nx = 2; // Energy level in x |
24 | | - int ny = 2; // Energy level in y |
25 | | - double dx = Lxy / m; // Step in x |
26 | | - double dy = Lxy / n; // Step in y |
27 | | - double dt = dx; // Time step size |
28 | | - |
29 | | - // Define staggered grids |
30 | | - vec xgrid = join_vert(vec({0}), linspace(dx / 2, Lxy - dx / 2, m), vec({Lxy})); |
31 | | - vec ygrid = join_vert(vec({0}), linspace(dy / 2, Lxy - dy / 2, n), vec({Lxy})); |
32 | | - |
33 | | - // Initialize 2D staggered grid |
34 | | - mat X, Y; |
35 | | - Utils utils; |
36 | | - utils.meshgrid(xgrid, ygrid, X, Y); |
37 | | - |
38 | | - // Initialize Laplacian operator with Robin BC |
39 | | - Laplacian L(k, m, n, dx, dy); |
40 | | - RobinBC BC(k, m, 1, n, 1, 1, 0); |
41 | | - L = L + BC; |
42 | | - |
43 | | - // Ensure the Laplacian is square |
44 | | - int total_size = (m + 2) * (n + 2); // Total size for the grid including boundaries |
45 | | - |
46 | | - // Hamiltonian definition |
47 | | - auto H = [&](const vec &x) { |
48 | | - vec result = 0.5 * (L * x); |
49 | | - return result; |
50 | | - }; |
51 | | - |
52 | | - // Define wave numbers |
53 | | - auto kx = [&](int nx) { return nx * M_PI / Lxy; }; |
54 | | - auto ky = [&](int ny) { return ny * M_PI / Lxy; }; |
55 | | - |
56 | | - double A = 2 / Lxy; |
57 | | - |
58 | | - // Initialize the wavefunction psi_old |
59 | | - mat Psi_grid(m + 2, n + 2, fill::zeros); |
60 | | - |
61 | | - // Manually pad Psi_grid |
62 | | - for (int i = 1; i <= m; i++) { |
63 | | - for (int j = 1; j <= n; j++) { |
64 | | - Psi_grid(i, j) = A * sin(kx(nx) * X(i, j)) * sin(ky(ny) * Y(i, j)); |
65 | | - } |
66 | | - } |
67 | | - |
68 | | - // Convert to column vector for compatibility |
69 | | - vec psi_old = vectorise(Psi_grid); // Convert to column vector |
70 | | - // Create interpolators |
71 | | - Interpol I(m, n, 0.5, 0.5); |
72 | | - Interpol I2(true, m, n, 0.5, 0.5); |
73 | | - |
74 | | - I *= dt; // Scale by dt |
75 | | - I2 *= 0.5 * dt; |
76 | | - |
77 | | - vec v_old(2*m*n+m+n, fill::zeros); // Initialize v_old |
78 | | - |
79 | | - // Initialize Psi_re with zeros |
80 | | - mat Psi_re = zeros<mat>(m + 2, n + 2); |
81 | | - |
82 | | - try { |
83 | | - // Time-stepping loop (Position Verlet) |
84 | | - for (int t = 0; t <= p; ++t) { |
85 | | - // Position Verlet algorithm: Update psi_old based on v_old |
86 | | - psi_old = psi_old + I2*v_old; |
87 | | - |
88 | | - // Calculate v_new using the Hamiltonian and psi_old |
89 | | - vec v_new = v_old + I*H(psi_old); |
90 | | - |
91 | | - // Update psi_old based on v_new |
92 | | - psi_old = psi_old + I2*v_new; |
93 | | - |
94 | | - // Update Psi_re for output |
95 | | - Psi_re = reshape(psi_old, m + 2, n + 2); |
96 | | - |
97 | | - if (t == p) { |
98 | | - std::cout << "Final Time Step " << t << ": X, Y, Psi" << std::endl; |
99 | | - for (size_t i = 0; i < Psi_re.n_rows; ++i) { |
100 | | - for (size_t j = 0; j < Psi_re.n_cols; ++j) { |
101 | | - std::cout << std::fixed << std::setprecision(5) |
102 | | - << X(i, j) << ", " << Y(i, j) << ", " << Psi_re(i, j) << std::endl; |
103 | | - } |
104 | | - } |
105 | | - } |
106 | | - |
107 | | - // Update variables for the next time step |
108 | | - v_old = v_new; |
109 | | - } |
110 | | - } catch (const std::exception &e) { |
111 | | - std::cerr << "Error: " << e.what() << std::endl; |
112 | | - return -1; |
113 | | - } |
114 | | - |
115 | | - std::cout << "Simulation complete." << std::endl; |
116 | | - return 0; |
117 | | -} |
| 115 | + std::cout << "Simulation complete." << std::endl; |
| 116 | + return 0; |
| 117 | +} |
0 commit comments