I have written a program that implements the ADI method and Crank-Nicolson method for solving Schrodinger equations.
The program is working, but it takes a very long time to run. I am looking for tips on how to improve the performance of the program.
#define _USE_MATH_DEFINES
#include <cmath>
#include <vector>
#include <string.h>
#include <array>
#include <iostream>
#include <complex>
#include <stdlib.h>
#include <fstream>
#include <stdio.h>
#include <string>
#include <math.h>
#include <iomanip>
#define PI M_PI
#define h_bar 1.0
#define M 1.0
#define L .5
#define sig .05
#define grid_point 100
#define K1 1.0E6
#define K2 1.0E6
#define omega1 std::sqrt(K1/M)
#define omega2 std::sqrt(K2/M)
typedef std::complex<double> compx;
#define I compx(0.,1.)
#define one compx(1.,0.)
#define two compx(2.,0.)
class wave_function {
public:
wave_function();
wave_function(bool);
double dt = 4.2E-06;
double dx = L / grid_point;
std::vector<std::vector<compx>> value;
void solve_triag();
double potential(int, int);
double real_space(int);
double density_x1[grid_point];
double density_x2[grid_point];
compx sec_space_deriv(char, int, int);
void normalize();
void rho();
void alpha_beta_solver(wave_function &, compx, std::vector<compx> &mid, compx, std::vector<compx> &R, int, char);
compx a = -h_bar / (two*M), b = -h_bar / (two*M);
double r = dt / (dx*dx);
compx A = (I*r*a / two), B = (I*r*b / two), C = I*dt / two;
};
wave_function::wave_function() {
value.resize(grid_point);
for (int l = 0; l < grid_point; l++) {
value[l].resize(grid_point);
for (int m = 0; m < grid_point; m++) {
value[l][m] = compx(0., 0.);
}
}
}
wave_function::wave_function(bool init) {
if (init) {
compx psi00, psi10, psi12;
value.resize(grid_point);
for (int l = 0; l < grid_point; l++) {
value[l].resize(grid_point);
for (int m = 0; m < grid_point; m++) {
/*value[l][m] = exp(I*compx(k1, 0)*compx(real_space(l), 0))*exp(-pow(real_space(l) - real_space(x_01), 2.) / (4.*sig*sig))
*exp(I*compx(k2, 0)*compx(real_space(m), 0))*exp(-pow(real_space(m) - real_space(x_02), 2.) / (4.*sig*sig));*/
psi00 = pow(M*omega1 / (PI*h_bar), .25)*exp(-M*omega1*pow(real_space(l), 2.) / (2.*h_bar))*pow(M*omega2 / (PI*h_bar), .25)*exp(-M*omega2*pow(real_space(m), 2.) / (2.*h_bar));
psi10 = pow(M*omega1 / (PI*h_bar), .25)*pow(2., .5)*pow(M*omega1 / (h_bar), .5)*real_space(l)*exp(-M*omega1*pow(real_space(l), 2.) / (2.*h_bar))*pow(M*omega2 / (PI*h_bar), .25)*exp(-M*omega2*pow(real_space(m), 2.) / (2.*h_bar));
value[l][m] = psi00+psi10;
}
}
normalize();
}
else if (!init) {
wave_function();
}
}
void wave_function::solve_triag() {
std::vector<compx> mid1, mid2, R1, R2;
wave_function tmp;
mid1.resize(grid_point);
R1.resize(grid_point);
for (int x2 = 0; x2 < grid_point; x2++) {
for (int i = 0; i < grid_point; i++) {
mid1[i] = one - two*A;
R1[i] = (one - C*potential(i, x2))*value[i][x2] - B*sec_space_deriv('y', i, x2);
}
for (int x1 = 0; x1 < grid_point; x1++) {
std::cout << x2 << " " << x1 << " calling solver" << std::endl;
alpha_beta_solver(tmp, A, mid1, A, R1, x2, 'x');
}
}
tmp.normalize();
mid2.resize(grid_point);
R2.resize(grid_point);
for (int x1 = 0; x1 < grid_point; x1++) {
for (int i = 0; i < grid_point; i++) {
mid2[i] = one - two*B + C*potential(x1, i);
R2[i] = tmp.value[x1][i] - A*tmp.sec_space_deriv('x', x1, i);
}
for (int x2 = 0; x2 < grid_point; x2++) {
alpha_beta_solver(*this, B, mid2, B, R2, x1, 'y');
}
}
normalize();
}
void wave_function::alpha_beta_solver(wave_function &New, compx b4_mid, std::vector<compx> &mid, compx post_mid, std::vector<compx> &R, int coor, char x_or_y) {
compx x_N, R_N;// new_mid[grid_point], new_R[grid_point];
std::vector<compx> alpha(grid_point - 2);
std::vector<compx> beta1(grid_point - 1);//for x_1
std::vector<compx> beta2(grid_point - 1);//for x_2
std::vector<compx> x_1(grid_point), x_2(grid_point);
alpha[0] = post_mid / mid[0];
beta1[0] = R[0] / mid[0];
beta2[0] = -b4_mid / mid[0];
//Forward run
for (int l = 1; l < grid_point - 2; l++) {
alpha[l] = post_mid / (mid[l] - b4_mid*alpha[l - 1]);
beta1[l] = (R[l] - b4_mid*beta1[l - 1]) / (mid[l] - b4_mid*alpha[l - 1]);
beta2[l] = (-b4_mid*beta2[l - 1]) / (mid[l] - b4_mid*alpha[l - 1]);
}
beta1[grid_point - 2] = (R[grid_point - 2] - b4_mid*beta1[grid_point - 3]) / (mid[grid_point - 2] - b4_mid*alpha[grid_point - 3]);
beta2[grid_point - 2] = (-post_mid - b4_mid*beta2[grid_point - 3]) / (mid[grid_point - 2] - b4_mid*alpha[grid_point - 3]);
//Backward run
x_1[grid_point - 2] = beta1[grid_point - 2];
x_2[grid_point - 2] = beta2[grid_point - 2];
for (int l = grid_point - 3; l >= 0; l--) {
x_1[l] = beta1[l] - alpha[l] * x_1[l + 1];
x_2[l] = beta2[l] - alpha[l] * x_2[l + 1];
}
x_N = (R_N - post_mid*x_1[0] - b4_mid*x_1[grid_point - 2])
/ (mid[grid_point - 2] + post_mid*x_2[0] + b4_mid*x_2[grid_point - 2]);
if (x_or_y == 'x') {
value[grid_point - 1][coor] = x_N;
for (int l = 0; l < grid_point - 1; l++) {
value[l][coor] = x_1[l] + x_2[l] * x_N;
//std::cout << x1 << " " << x2 << " " << l << " " << x_1[l] + x_2[l] * x_N << std::endl;
}
}
else if (x_or_y == 'y') {
value[coor][grid_point - 1] = x_N;
for (int l = 0; l < grid_point - 1; l++) {
value[coor][l] = x_1[l] + x_2[l] * x_N;
}
}
}
double wave_function::potential(int x1, int x2) {
return .5*M*omega1*omega1*(real_space(x1))*(real_space(x1)) + .5*M*omega1*omega1*(real_space(x2))*(real_space(x2)) + 0.0*.5*M*omega2*omega2*pow(real_space(x2) - real_space(x1), 2.);
//square well
/*if (al - abs(real_space(x1) - real_space(x2)) > 0.)
return V_0;
else
return 0.0;*/
//gaussian
//return V_0*exp(-pow(abs(real_space(x1)-real_space(x2)),2.)/(2.*al*al));
}
double wave_function::real_space(int index) {
return (-L / 2.0) + (index*dx);
}
//x stands for x1 or particle 1 and y stands for x2 or particle 2
compx wave_function::sec_space_deriv(char x1_or_x2, int l, int m) {
if (x1_or_x2 == 'x') {
if (l == 0)
return value[l + 1][m] - two * value[l][m] + value[grid_point - 1][m];
else if (l == grid_point - 1)
return value[0][m] - two * value[l][m] + value[l - 1][m];
else
return value[l + 1][m] - two * value[l][m] + value[l - 1][m];
}
else if (x1_or_x2 == 'y') {
if (m == 0)
return value[l][m + 1] - two * value[l][m] + value[l][grid_point - 1];
else if (m == grid_point - 1)
return value[l][0] - two * value[l][m] + value[l][m - 1];
else
return value[l][m + 1] - two * value[l][m] + value[l][m - 1];
}
}
template <typename NEW>
std::string to_string_with_precision(const NEW a_value, const int n = 6)
{
std::ostringstream out;
out << std::fixed << std::setprecision(n) << a_value;
return out.str();
}
//normalization for the function above rho
void wave_function::normalize() {
compx sum = 0;
for (int i1 = 0; i1 < grid_point; i1++) {
for (int i2 = 0; i2 < grid_point; i2++) {
sum += pow(abs(value[i1][i2]), 2.0);
}
}
sum *= dx * dx;
compx amplitude = one / pow(sum, .5);
for (int i1 = 0; i1 < grid_point; i1++) {
for (int i2 = 0; i2 < grid_point; i2++) {
value[i1][i2] *= amplitude;
}
}
}
void wave_function::rho() {
for (int k = 0; k < grid_point; k++) {
density_x1[k] = 0.;
density_x2[k] = 0.;
}
for (int i = 0; i < grid_point; i++) {
for (int j = 0; j < grid_point; j++) {
density_x1[i] += pow(abs(value[i][j]), 2.)*dx;
density_x2[i] += pow(abs(value[j][i]), 2.)*dx;
}
}
}
int main() {
double check = 0.0;
wave_function v(true);
std::ofstream file1; //density of both particles
std::ofstream file2;
std::ofstream file5; //potential
file5.open("potential.dat");
for (int i = 0; i < grid_point; i++)
for (int j = 0; j < grid_point; j++)
file5 << v.real_space(i) << "\t" << v.real_space(j) << "\t" << v.potential(i, j) << std::endl;
int index = 0;
for (double k = v.dt; k <= 200 * v.dt; k += v.dt) {
std::cout << index << std::endl;
v.rho();
if (index % 10 == 0) {
file1.open("data_" + to_string_with_precision(index, 0) + ".dat");
file1 << "Time " << k - v.dt << std::endl
<< "x" << "\t" << "y" << "\t" << "imag" << "\t" << "real" << "\t" << "abs" << std::endl;
file2.open("x_" + std::to_string(index) + ".dat");
file2 << "Time " << k - v.dt << std::endl
<< "Coordinate" << "\t" << "x1" << "\t" << "x2" << std::endl;
for (int i = 0; i < grid_point; i++) {
for (int j = 0; j < grid_point; j++) {
file1 << v.real_space(i) << "\t"
<< v.real_space(j) << "\t"
<< imag(v.value[i][j]) << "\t"
<< real(v.value[i][j]) << "\t"
<< abs(v.value[i][j]) << std::endl;
}
}
for (int k = 0; k < grid_point; k++) {
file2 << v.real_space(k) << "\t"
<< abs(v.density_x1[k]) << "\t"
<< abs(v.density_x2[k]) << std::endl;
}
file2.close();
file1.close();
}
std::cout << "About to call solve_triag() " << std::endl;
v.solve_triag();
index++;
}
getchar();
}
{}
button or Ctrl+K to mark it as a code sample). \$\endgroup\$main()
that generates a big enough input set on its own? That would enable some simple benchmarking and help filter out the good suggestions to put in answers. \$\endgroup\$