Skip to content
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Removed spacing in empty lines
  • Loading branch information
rajdakin committed Oct 13, 2018
commit eb0b0dde891dc9aa63cbfcb467c85b47947ca87f
70 changes: 35 additions & 35 deletions contents/split-operator_method/code/c++/split_op.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ struct Params {
dk = M_PI / xmax;
k.reserve(res);
im_time = im;

for (size_t i = 0; i < res; ++i) {
x.push_back(xmax / res - xmax + i * (2.0 * xmax / res));
if (i < res / 2) {
Expand All @@ -36,7 +36,7 @@ struct Params {
}
}
}

double xmax;
unsigned int res;
double dt;
Expand All @@ -57,11 +57,11 @@ struct Operators {
pe.reserve(size);
ke.reserve(size);
wfc.reserve(size);

for (size_t i = 0; i < size; ++i) {
v.emplace_back(0.5 * pow(par.x[i] - voffset, 2));
wfc.emplace_back(exp(-pow(par.x[i] - wfcoffset, 2) / 2.0));

if (par.im_time) {
ke.emplace_back(exp(-0.5 * par.dt * pow(par.k[i], 2)));
pe.emplace_back(exp(-0.5 * par.dt * v[i]));
Expand All @@ -71,7 +71,7 @@ struct Operators {
}
}
}

size_t size;
vector_complex v;
vector_complex pe;
Expand All @@ -83,75 +83,75 @@ void fft(vector_complex &x, int n, bool inverse) {
complex y[n];
memset(y, 0, sizeof(y));
fftw_plan p;

fftw_complex *in = reinterpret_cast<fftw_complex*>(x.data());
fftw_complex *out = reinterpret_cast<fftw_complex*>(y);
p = fftw_plan_dft_1d(n, in, out,
(inverse ? FFTW_BACKWARD : FFTW_FORWARD), FFTW_ESTIMATE);

fftw_execute(p);
fftw_destroy_plan(p);

for (size_t i = 0; i < n; ++i) {
x[i] = y[i] / sqrt(static_cast<double>(n));
}
}

void split_op(Params &par, Operators &opr) {
double density[opr.size];

for (size_t i = 0; i < par.timesteps; ++i) {
for (size_t j = 0; j < opr.size; ++j) {
opr.wfc[j] *= opr.pe[j];
}

fft(opr.wfc, opr.size, false);

for (size_t j = 0; j < opr.size; ++j) {
opr.wfc[j] *= opr.ke[j];
}

fft(opr.wfc, opr.size, true);

for (size_t j = 0; j < opr.size; ++j) {
opr.wfc[j] *= opr.pe[j];
}

for (size_t j = 0; j < opr.size; ++j) {
density[j] = pow(abs(opr.wfc[j]), 2);
}

if (par.im_time) {
double sum = 0;

for (size_t j = 0; j < opr.size; ++j) {
sum += density[j];
}

sum *= par.dx;

for (size_t j = 0; j < opr.size; ++j) {
opr.wfc[j] /= sqrt(sum);
}
}

// Writing data into a file in the format of:
// index, density, real potential.
std::stringstream filename_stream;
filename_stream << "output" << i << ".dat";

std::ofstream fstream = std::ofstream(filename_stream.str());

if (fstream) {
for (int i = 0; i < opr.size; ++i) {
std::stringstream data_stream;

data_stream << i << "\t" << density[i] << "\t" << real(opr.v[i]) << "\n";

fstream.write(data_stream.str().c_str(), data_stream.str().length());
}
}

fstream.close();
}
}
Expand All @@ -161,41 +161,41 @@ double calculate_energy(Params &par, Operators &opr) {
vector_complex wfc_k(opr.wfc);
vector_complex wfc_c(opr.size);
fft(wfc_k, opr.size, false);

for (size_t i = 0; i < opr.size; ++i) {
wfc_c.push_back(conj(wfc_r[i]));
}

vector_complex energy_k(opr.size);
vector_complex energy_r(opr.size);

for (size_t i = 0; i < opr.size; ++i) {
energy_k.push_back(wfc_k[i] * pow(complex(par.k[i], 0.0), 2));
}

fft(energy_k, opr.size, true);

for (size_t i = 0; i < opr.size; ++i) {
energy_k[i] *= 0.5 * wfc_c[i];
energy_r.push_back(wfc_c[i] * opr.v[i] * wfc_r[i]);
}

double energy_final = 0;

for (size_t i = 0; i < opr.size; ++i) {
energy_final += real(energy_k[i] + energy_r[i]);
}

return energy_final * par.dx;
}

int main() {
Params par = Params(5.0, 256, 0.05, 100, true);
Operators opr = Operators(par, 0.0, -1.0);

split_op(par, opr);

std::cout << "The energy is " << calculate_energy(par, opr) << "\n";

return 0;
}