Commit 8b394515 authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

single plasmon scattering

parent 5f9eafdd
......@@ -21,6 +21,7 @@
#include "IO.hpp"
#include <sys/stat.h>
#include <malloc.h>
#include <regex>
......@@ -927,8 +928,16 @@ void IO::dumpWave(string filename, const Wave &wave) {
fclose(f);
}
bool exists(const std::string& name) {
struct stat buffer;
return (stat (name.c_str(), &buffer) == 0);
}
std::shared_ptr<stemsalabim::atomic::Cell> IO::initCrystalFromXYZFile(const std::string &filename) {
if(!exists(filename))
output::error("Crystal file %s can't be found!\n", filename);
atomic::ElementProvider elp = atomic::ElementProvider::getInstance();
// some useful variable declarations
......
......@@ -278,7 +278,9 @@ void Simulation::run() {
}
void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
seconds Simulation::calculatePixel(ScanPoint &point, double defocus) {
auto start_time = high_resolution_clock::now();
Params &p = Params::getInstance();
atomic::Plasmons &plasmons = atomic::Plasmons::getInstance();
......@@ -320,19 +322,31 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
// zero loss
if(p.plasmonsEnabled()) {
// scattering
for(unsigned long energy_i = _gridman->energyLossGrid().size() - 1ul; energy_i > 0; energy_i--) {
// for now, single scattering only...
for(unsigned int energy_i = _gridman->energyLossGrid().size() - 1u; energy_i > 0; energy_i--) {
waves[energy_i] *= sqrt(plasmons.beta0());
// plural scattering includes single scattering with energy_j = 0
for(unsigned int energy_j = 0; energy_j < (p.plasmonsPluralScattering() ? energy_i : 1u); energy_j++) {
tmpwave = waves[energy_j];
tmpwave *= plasmons.scatteringFunction(energy_i - energy_j - 1);
tmpwave *= sqrt(plasmons.beta1());
waves[energy_i] += tmpwave;
}
tmpwave = waves[0];
tmpwave *= plasmons.scatteringFunction(energy_i - 1);
tmpwave *= sqrt(plasmons.beta1());
waves[energy_i] += tmpwave;
}
//
// // scattering
// for(unsigned int energy_i = _gridman->energyLossGrid().size() - 1u; energy_i > 0; energy_i--) {
//
// waves[energy_i] *= sqrt(plasmons.beta0());
//
// // plural scattering includes single scattering with energy_j = 0
// for(unsigned int energy_j = 0; energy_j < (p.plasmonsPluralScattering() ? energy_i : 1u); energy_j++) {
// tmpwave = waves[energy_j];
// tmpwave *= plasmons.scatteringFunction(energy_i - energy_j - 1);
// tmpwave *= sqrt(plasmons.beta1());
// waves[energy_i] += tmpwave;
// }
// }
//
waves[0] *= sqrt(plasmons.beta0());
......@@ -479,6 +493,9 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
if(point.cbed)
point.storeCBEDIntensities(cbed_intensities, _cbed_intensity_buffer);
auto finish_time = high_resolution_clock::now();
return duration_cast<seconds>(finish_time - start_time);
}
void Simulation::printLine(const SimulationState &st, const string &line) {
......@@ -548,6 +565,9 @@ void Simulation::multisliceMaster(const SimulationState &st, bool do_mpi_queue)
TaskQueue scan_work;
scan_work.append(work_packages);
atomic_ulong px_counter = 0;
atomic_ullong px_calc_seconds = 0;
// keep the master's worker threads busy as well. We will fetch from the task queue
// one by one AFTER the previous task is finished, as the queue is worked on in parallel
// via MPI. It is very important not to distribute all tasks here.
......@@ -556,7 +576,9 @@ void Simulation::multisliceMaster(const SimulationState &st, bool do_mpi_queue)
threads.emplace_back([&] {
unsigned int pix_index = 0;
while(scan_work.pop(pix_index)) {
calculatePixel(work_packages[pix_index], st.defocus());
px_counter++;
auto duration = calculatePixel(work_packages[pix_index], st.defocus());
px_calc_seconds += duration.count() / px_counter;
scan_work.finish(pix_index);
}
});
......@@ -608,18 +630,17 @@ void Simulation::multisliceMaster(const SimulationState &st, bool do_mpi_queue)
} else {
microseconds io_duration(0);
// do some more work on the main thread.
unsigned int pix_index = 0;
unsigned int last_prcnt = 0;
while(scan_work.pop(pix_index)) {
calculatePixel(work_packages[pix_index], st.defocus());
px_counter++;
auto duration = calculatePixel(work_packages[pix_index], st.defocus());
px_calc_seconds += duration.count() / px_counter;
scan_work.finish(pix_index);
// store finished intensities
storeResultData(st.idefocus(), st.iconf(), scan_work.getDirtyAndClean(), work_packages);
float p = scan_work.progress();
......@@ -637,7 +658,7 @@ void Simulation::multisliceMaster(const SimulationState &st, bool do_mpi_queue)
storeResultData(st.idefocus(), st.iconf(), scan_work.getDirtyAndClean(), work_packages);
printLine(st, output::fmt("Calculation took %s.\n", output::humantime(algorithms::getTimeSince(t1))));
printLine(st, output::fmt("Calculation took %s\n", output::humantime(algorithms::getTimeSince(t1))));
}
// join the threads
......
......@@ -85,9 +85,6 @@ namespace stemsalabim {
/// all IO functions, i.e., NC results file management.
std::shared_ptr<IO> _io;
/// A timer to prevent sending HTTP requests too often.
std::chrono::high_resolution_clock::time_point _http_request_timer{std::chrono::high_resolution_clock::now()};
/// The propagator can be represented as two 1D functions, one for each
/// spatial direction.
Wave _propagator;
......@@ -103,7 +100,7 @@ namespace stemsalabim {
* and a specific defocus value. Returns a vector with the intensities, that is aligned by slices * angles, i.e.,
* having M angles and N slices, [0, M[ are the angular intensity of slice 0, [M, 2M[ the second slice etc.
*/
void calculatePixel(ScanPoint &point, const double defocus);
std::chrono::seconds calculatePixel(ScanPoint &point, double defocus);
/*!
* The part of the multi-slice simulation that is carried out by the MPI master, i.e., rank 0.
......
......@@ -71,6 +71,9 @@ namespace stemsalabim { namespace atomic {
_cache[eii].init(lx, ly);
_cache[eii].setIsKSpace(true);
// if(eii != 16)
// continue;
for(unsigned int ix = 0; ix < lx; ix++) {
for(unsigned int iy = 0; iy < ly; iy++) {
double theta = sqrt(pow(gridman->kx(ix), 2) + pow(gridman->ky(iy), 2)) * p.wavelength();
......@@ -95,13 +98,12 @@ namespace stemsalabim { namespace atomic {
w,
&result,
&error);
//result = func((gridman->energyLossGrid()[eii] + gridman->energyLossGrid()[eii+1])*0.5);
_cache[eii](ix, iy) = (float) sqrt(result);
}
}
}
//_cache[eii](0, 0) = _cache[eii](0, 1);
//_cache[eii].normalize();
_cache[eii].backwardFFT();
}
......@@ -122,6 +124,10 @@ namespace stemsalabim { namespace atomic {
return _beta_1;
}
double beta2() const {
return _beta_2;
}
private:
explicit Plasmons() {
......@@ -132,6 +138,7 @@ namespace stemsalabim { namespace atomic {
_incident_energy = p.beamEnergy();
_beta_0 = exp(-p.sliceThickness() / p.plasmonMeanFreePath());
_beta_1 = p.sliceThickness() / p.plasmonMeanFreePath() * _beta_0;
_beta_2 = p.sliceThickness() / p.plasmonMeanFreePath() * 0.5 * _beta_1;
}
/*!
......@@ -167,10 +174,13 @@ namespace stemsalabim { namespace atomic {
3. / 5. * fermi_energy / _plasmon_energy * 2 * _incident_energy * 1000 *
(theta * theta + characteristic_angle * characteristic_angle);
// double g = 3. / 5. * fermi_energy / _plasmon_energy;
// double differential cross section of plasmon scattering in angstroms^2/eV
double B = sin(theta) / (theta * theta + characteristic_angle * characteristic_angle);
B *= E * _plasmon_fwhm * Ep * Ep / (pow(E * E - Ep * Ep, 2.) + pow(E * _plasmon_fwhm, 2.));
//B /= 1000;
// double B = sin(theta) / (theta * theta + characteristic_angle * characteristic_angle);
// B *= E * _plasmon_fwhm * _plasmon_energy * _plasmon_energy / (pow(E * E - _plasmon_energy * _plasmon_energy - 4*_plasmon_energy*_incident_energy*1000*g*(theta * theta + characteristic_angle * characteristic_angle), 2.) + pow(E * _plasmon_fwhm, 2.));
return B;
}
......@@ -180,6 +190,7 @@ namespace stemsalabim { namespace atomic {
double _plasmon_fwhm;
double _beta_0;
double _beta_1;
double _beta_2;
bool _cache_populated{false};
std::vector<Wave> _cache;
......
......@@ -205,6 +205,8 @@ namespace stemsalabim {
el_str += output::fmt(" (%d)", counter);
}
float kb_to_GB = powf(1024.f, 3);
op::nakedprint("\n\n");
op::nakedprint("%s v%s.%s.%s\n", PKG_NAME, PKG_VERSION_MAJOR, PKG_VERSION_MINOR, PKG_VERSION_PATCH);
op::nakedprint(" git commit %s\n", GIT_COMMIT_HASH);
......@@ -230,17 +232,17 @@ namespace stemsalabim {
op::nakedprint("\nMemory estimates (lower bounds!!)\n");
op::nakedprint("-----------------------------------\n");
op::nakedprint(" Working mem/proc: %.1f GB\n", mem_per_proc / (float) (1024 * 1024 * 1024));
op::nakedprint(" Results mem/proc: %.1f GB\n", result_mem_per_proc / (float) (1024 * 1024 * 1024));
op::nakedprint(" Working mem/proc: %.1f GB\n", mem_per_proc / kb_to_GB);
op::nakedprint(" Results mem/proc: %.1f GB\n", result_mem_per_proc / kb_to_GB);
op::nakedprint(" Total mem/proc: %.1f GB\n",
(result_mem_per_proc + mem_per_proc) / (float) (1024 * 1024 * 1024));
(result_mem_per_proc + mem_per_proc) / kb_to_GB);
op::nakedprint("\nOutput file size\n");
op::nakedprint("-----------------------------------\n");
op::nakedprint(" ADF storage: %.1f GB\n", num_adf / (float) (1024 * 1024 * 1024));
op::nakedprint(" Center of mass: %.1f GB\n", num_com / (float) (1024 * 1024 * 1024));
op::nakedprint(" CBED storage: %.1f GB\n", num_cbed / (float) (1024 * 1024 * 1024));
op::nakedprint(" Total size: %.1f GB\n", (num_cbed + num_adf + num_com) / (float) (1024 * 1024 * 1024));
op::nakedprint(" ADF storage: %.1f GB\n", num_adf / kb_to_GB);
op::nakedprint(" Center of mass: %.1f GB\n", num_com / kb_to_GB);
op::nakedprint(" CBED storage: %.1f GB\n", num_cbed / kb_to_GB);
op::nakedprint(" Total size: %.1f GB\n", (num_cbed + num_adf + num_com) / kb_to_GB);
op::nakedprint("\n\n");
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment