Commit 642f8bae authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

Improvements

parent 9ef7059c
......@@ -175,7 +175,7 @@ void IO::initNCFile(const shared_ptr<GridManager> &gridman, const shared_ptr<ato
auto v = g.defineVar<float>("adf_intensities",
vd({defocusDim, xDim, yDim, fpDim, sliceDim, plasmonDim, angleDim}));
g.defineVar<float>("center_of_mass", vd({defocusDim, xDim, yDim, fpDim, sliceDim, coordinateDim}));
g.defineVar<float>("center_of_mass", vd({defocusDim, xDim, yDim, fpDim, sliceDim, plasmonDim, coordinateDim}));
g.defineVar<double>("adf_probe_x_grid", vd({xDim}));
g.defineVar<double>("adf_probe_y_grid", vd({yDim}));
g.defineVar<double>("adf_detector_grid", vd({angleDim}));
......@@ -251,7 +251,7 @@ void IO::initBuffers(const shared_ptr<GridManager> &gridman) {
_adf_intensities_buffer.resize(gridman->adfDetectorGrid().size() *
gridman->adfSliceCoords().size() *
gridman->energyLossGrid().size(), 0);
_com_buffer.resize(2 * gridman->adfSliceCoords().size(), 0);
_com_buffer.resize(2 * gridman->adfSliceCoords().size() * gridman->energyLossGrid().size(), 0);
}
if(p.cbed()) {
......@@ -520,6 +520,7 @@ IO::writeCOM(unsigned int idefocus, unsigned int iconf, const shared_ptr<GridMan
unsigned int idefocus_store = p.adfAverageDefoci() ? 0 : idefocus;
unsigned int iconf_store = p.adfAverageConfigurations() ? 0 : iconf;
unsigned long nplasmons = gridman->energyLossGrid().size();
float w = 1.0;
......@@ -531,8 +532,8 @@ IO::writeCOM(unsigned int idefocus, unsigned int iconf, const shared_ptr<GridMan
if(iconf > iconf_store || idefocus > idefocus_store) {
v.get(vs({idefocus_store, point.adf_row, point.adf_col, iconf_store, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, 2}),
v.get(vs({idefocus_store, point.adf_row, point.adf_col, iconf_store, 0, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nplasmons, 2}),
_com_buffer);
for(size_t j = 0; j < _com_buffer.size(); ++j)
......@@ -545,8 +546,8 @@ IO::writeCOM(unsigned int idefocus, unsigned int iconf, const shared_ptr<GridMan
}
v.put(vs({idefocus_store, point.adf_row, point.adf_col, iconf_store, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, 2}),
v.put(vs({idefocus_store, point.adf_row, point.adf_col, iconf_store, 0, 0, 0}),
vs({1, 1, 1, 1, n_slices_total, nplasmons, 2}),
_com_buffer);
} catch(NcException &e) {
......
......@@ -21,6 +21,7 @@
#include "Simulation.hpp"
#include <thread>
#include "../utilities/initialize.hpp"
#include "../libatomic/Scattering.hpp"
#include "../utilities/timings.hpp"
#include "../libatomic/Plasmons.hpp"
......@@ -50,17 +51,6 @@ void Simulation::init() {
_fpman = std::make_shared<FPConfManager>(p.cell());
_io = std::make_shared<IO>();
if(mpi_env.isMpi())
output::print("[%s, %d] Hello from %s (rank %d/%d), spawning %d threads.\n",
mpi_env.name(),
mpi_env.rank(),
mpi_env.name(),
mpi_env.rank(),
mpi_env.size(),
p.numberOfThreads());
else
output::print("This is a single-node calculation using %d threads.\n", p.numberOfThreads());
// Read and communicate the mean square displacements of the frozen phonon
// approximation.
{
......@@ -87,11 +77,28 @@ void Simulation::init() {
_gridman->generateGrids();
if(mpi_env.isMaster())
printStatus(_gridman);
if(mpi_env.isMpi())
output::print("[%s, %d] Hello from %s (rank %d/%d), spawning %d threads.\n",
mpi_env.name(),
mpi_env.rank(),
mpi_env.name(),
mpi_env.rank(),
mpi_env.size(),
p.numberOfThreads());
else
output::print("This is a single-node calculation using %d threads.\n", p.numberOfThreads());
{
// plasmon stuff. Needs to come after crystal and Gridmanager initialization!
if(p.plasmonsEnabled()) {
printMaster("Initializing plasmon scattering functions...\n");
atomic::Plasmons &plasmons = atomic::Plasmons::getInstance();
plasmons.populateCache(_gridman, p.cell()->electronDensity(), p.cell()->atomicDensity());
plasmons.populateCache(_gridman, p.cell()->electronDensity());
}
}
......@@ -117,15 +124,6 @@ void Simulation::init() {
}
_io->initBuffers(_gridman);
printMaster(output::fmt("Calculating on a real and k space lattice of %dx%d points.\n",
_gridman->samplingX(),
_gridman->samplingY()));
printMaster(output::fmt("Calculating %d defoci, %d FP confs, %dx%d scan points.\n",
p.numberOfDefoci(),
p.numberOfConfigurations(),
_gridman->adfXGrid().size(),
_gridman->adfYGrid().size()));
printMaster(output::fmt("Initialization took %s.\n", output::humantime(algorithms::getTimeSince(init_start))));
initBuffers();
......@@ -290,8 +288,10 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
vector<float> detector_intensities;
if(point.adf) {
adf_intensities.reserve(_gridman->adfSliceCoords().size() * _gridman->adfDetectorGrid().size() * _gridman->energyLossGrid().size());
center_of_mass.reserve(_gridman->adfSliceCoords().size() * 2);
adf_intensities.reserve(_gridman->adfSliceCoords().size() *
_gridman->adfDetectorGrid().size() *
_gridman->energyLossGrid().size());
center_of_mass.reserve(_gridman->adfSliceCoords().size() * _gridman->energyLossGrid().size() * 2);
detector_intensities.resize(_gridman->adfDetectorGrid().size());
}
......@@ -321,7 +321,7 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
if(p.plasmonsEnabled()) {
// scattering
for(unsigned long energy_i = _gridman->energyLossGrid().size()-1ul; energy_i > 0; energy_i--) {
for(unsigned long energy_i = _gridman->energyLossGrid().size() - 1ul; energy_i > 0; energy_i--) {
waves[energy_i] *= sqrt(plasmons.beta0());
......@@ -353,7 +353,7 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
if(p.plasmonsEnabled()) {
double total_intensity = 0;
double desired_intensity = 1-waves[0].integratedIntensity();
double desired_intensity = 1 - waves[0].integratedIntensity();
for(unsigned long ei = 1; ei < _gridman->energyLossGrid().size(); ei++)
total_intensity += waves[ei].integratedIntensity();
......@@ -388,22 +388,24 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
adf_intensities.insert(adf_intensities.end(), detector_intensities.begin(), detector_intensities.end());
}
// calculate center of mass
double total_intensity = waves[0].integratedIntensity(), cx = 0, cy = 0;
for(unsigned int ix = 0; ix < _gridman->samplingX(); ix++) {
for(unsigned int iy = 0; iy < _gridman->samplingY(); iy++) {
auto intens = (float) pow(abs(waves[0](ix, iy)), 2);
cx += intens * _gridman->kx(ix);
cy += intens * _gridman->ky(iy);
// calculate center of mass
double total_intensity = waves[ip].integratedIntensity(), cx = 0, cy = 0;
for(unsigned int ix = 0; ix < _gridman->samplingX(); ix++) {
for(unsigned int iy = 0; iy < _gridman->samplingY(); iy++) {
auto intens = (float) pow(abs(waves[ip](ix, iy)), 2);
cx += intens * _gridman->kx(ix);
cy += intens * _gridman->ky(iy);
}
}
cx /= total_intensity;
cy /= total_intensity;
center_of_mass.push_back((float) (cx * p.wavelength() * 1e3));
center_of_mass.push_back((float) (cy * p.wavelength() * 1e3));
}
cx /= total_intensity;
cy /= total_intensity;
center_of_mass.push_back((float) (cx * p.wavelength() * 1e3));
center_of_mass.push_back((float) (cy * p.wavelength() * 1e3));
}
......@@ -729,7 +731,9 @@ void Simulation::initBuffers() {
_gridman->cbedSliceCoords().size() *
_gridman->energyLossGrid().size()) : 0;
_com_buffer = std::make_shared<buf_type>(prms.adf() ? (2 * _gridman->adfSliceCoords().size()) : 0,
_com_buffer = std::make_shared<buf_type>(prms.adf() ? (2 *
_gridman->adfSliceCoords().size() *
_gridman->energyLossGrid().size()) : 0,
prms.workPackageSize());
_adf_intensity_buffer = std::make_shared<buf_type>(number_intensities_per_pixel, prms.workPackageSize());
_cbed_intensity_buffer = std::make_shared<buf_type>(number_cbed_per_pixel, prms.workPackageSize());
......
......@@ -43,7 +43,7 @@ namespace stemsalabim { namespace atomic {
return instance;
}
void populateCache(const std::shared_ptr<GridManager> &gridman, double electron_density, double atom_density) {
void populateCache(const std::shared_ptr<GridManager> &gridman, double electron_density) {
if(_cache_populated)
return;
......@@ -80,7 +80,7 @@ namespace stemsalabim { namespace atomic {
// integrate over energy
std::function<double(double)> func = [&](double x) {
return plasmonDoubleDifferentialCrossSection(x, theta, atom_density, electron_density);
return plasmonDoubleDifferentialCrossSection(x, theta, electron_density);
};
F.params = &func;
......@@ -142,12 +142,12 @@ namespace stemsalabim { namespace atomic {
* @return double differential cross section in angstroms^2/eV
*/
double
plasmonDoubleDifferentialCrossSection(double E, double theta, double atom_density, double electron_density) {
plasmonDoubleDifferentialCrossSection(double E, double theta, double electron_density) {
// calculate fermi energy in eV
// this factor is pow(3*pi^2, 2./3.)
constexpr double numerical_factor = 9.570780000627304;
// this factor is hbar^2 / (2 * m_e * e) * 10^20 (eV A^2)
// this factor is hbar^2 / (2 * m_e * e) * 10^20 (eV nm^2)
constexpr double numerical_factor_2 = 0.03809982080226884;
constexpr double fermi_factor = numerical_factor_2 * numerical_factor;
......@@ -162,9 +162,6 @@ namespace stemsalabim { namespace atomic {
1e-3;
double factor = 2 * _incident_energy * 1000 * (theta * theta + characteristic_angle * characteristic_angle);
double g = 3. / 5. * fermi_energy / _plasmon_energy;
// plasmon dispersion
double Ep = _plasmon_energy +
3. / 5. * fermi_energy / _plasmon_energy * 2 * _incident_energy * 1000 *
......
......@@ -26,17 +26,18 @@
#include "utilities/Wave.hpp"
#include "classes/Simulation.hpp"
#include "utilities/initialize.hpp"
#include "config.h"
using namespace std;
using namespace stemsalabim;
namespace op = stemsalabim::output;
const char *getTmpDir (void) {
const char *getTmpDir(void) {
char *tmpdir;
if ((tmpdir = getenv ("TEMP")) != NULL) return tmpdir;
if ((tmpdir = getenv ("TMP")) != NULL) return tmpdir;
if ((tmpdir = getenv ("TMPDIR")) != NULL) return tmpdir;
if((tmpdir = getenv("TEMP")) != NULL) return tmpdir;
if((tmpdir = getenv("TMP")) != NULL) return tmpdir;
if((tmpdir = getenv("TMPDIR")) != NULL) return tmpdir;
return "/tmp";
}
......@@ -124,44 +125,7 @@ int main(int argc, const char **argv) {
// -------------- checks
double kx_max = 1.0 / (3.0 * max(1. / gridman->scaleX(), 1. / gridman->scaleY())) * 1000;
size_t mem_per_proc = 0, result_mem_per_proc = 0;
size_t fs = sizeof(float);
// transmission functions
mem_per_proc += gridman->samplingX() * gridman->samplingY() * 2 * fs * gridman->slices().size();
// wave functions for each thread during pixel calculation
mem_per_proc += gridman->samplingX() * gridman->samplingY() * 2 * fs * (
gridman->energyLossGrid().size() + 2) * p.numberOfThreads();
// results
unsigned long number_intensities_per_pixel = p.adf() ? (gridman->adfDetectorGrid().size() *
gridman->adfSliceCoords().size() *
gridman->energyLossGrid().size()) : 0;
unsigned long number_cbed_per_pixel = p.cbed() ? (gridman->storedCbedSizeX() *
gridman->storedCbedSizeY() *
gridman->cbedSliceCoords().size() *
gridman->energyLossGrid().size()) : 0;
// buffers at least!
result_mem_per_proc += p.workPackageSize() * (number_intensities_per_pixel + number_cbed_per_pixel) * fs;
op::nakedprint("General information\n");
op::nakedprint("-------------------------------\n");
op::nakedprint(" Grid size: %d x %d\n", gridman->samplingX(), gridman->samplingY());
op::nakedprint(" Scan points: %d x %d\n", gridman->adfXGrid().size(), gridman->adfYGrid().size());
op::nakedprint(" Energy bins: %d\n", gridman->energyLossGrid().size());
op::nakedprint(" Max angle %.1f mrad\n", kx_max*p.wavelength());
op::nakedprint(" Angle resolution: %.3f mrad\n", gridman->kx(1)*p.wavelength()*1000);
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(" Total mem/proc: %.1f GB\n", (result_mem_per_proc+mem_per_proc)/(float)(1024*1024*1024));
printStatus(gridman);
// -------------- stop checks
......
......@@ -26,6 +26,7 @@
#include <thread>
#include "mpi.hpp"
#include "config.h"
#include "TaskQueue.hpp"
#include "../classes/Params.hpp"
#include "../classes/IO.hpp"
......@@ -34,7 +35,10 @@
namespace stemsalabim {
void initialize_input_file(int argc, const char **argv, bool silent = false) {
using namespace std;
namespace op = stemsalabim::output;
inline void initialize_input_file(int argc, const char **argv, bool silent = false) {
auto &mpi_env = mpi::Environment::getInstance();
Params &p = Params::getInstance();
......@@ -101,5 +105,142 @@ namespace stemsalabim {
}
}
inline void printStatus(const std::shared_ptr<GridManager> & gridman) {
Params &p = Params::getInstance();
double kx_max = 1.0 / (3.0 * max(1. / gridman->scaleX(), 1. / gridman->scaleY())) * 1000;
size_t mem_per_proc = 0, result_mem_per_proc = 0;
size_t fs = sizeof(float);
// transmission functions
mem_per_proc += gridman->samplingX() * gridman->samplingY() * 2 * fs * gridman->slices().size();
// wave functions for each thread during pixel calculation
mem_per_proc += gridman->samplingX() *
gridman->samplingY() *
2 *
fs *
(gridman->energyLossGrid().size() + 2) *
p.numberOfThreads();
// results
unsigned long number_intensities_per_pixel = p.adf() ? (gridman->adfDetectorGrid().size() *
gridman->adfSliceCoords().size() *
gridman->energyLossGrid().size()) : 0;
unsigned long number_cbed_per_pixel = p.cbed() ? (gridman->storedCbedSizeX() *
gridman->storedCbedSizeY() *
gridman->cbedSliceCoords().size() *
gridman->energyLossGrid().size()) : 0;
// buffers at least!
result_mem_per_proc += p.workPackageSize() * (number_intensities_per_pixel + number_cbed_per_pixel) * fs;
// output file size
// adf
unsigned long num_adf = gridman->adfDetectorGrid().size() *
gridman->adfSliceCoords().size() *
gridman->adfXGrid().size() *
gridman->adfYGrid().size() *
gridman->energyLossGrid().size() *
fs;
if(!p.adfAverageConfigurations())
num_adf *= p.numberOfConfigurations();
if(!p.adfAverageDefoci())
num_adf *= p.numberOfDefoci();
// cbed
unsigned long num_cbed = gridman->storedCbedSizeX() *
gridman->storedCbedSizeY() *
gridman->cbedSliceCoords().size() *
gridman->cbedXGrid().size() *
gridman->cbedYGrid().size() *
gridman->energyLossGrid().size() *
fs;
if(!p.cbedAverageConfigurations())
num_cbed *= p.numberOfConfigurations();
if(!p.cbedAverageDefoci())
num_cbed *= p.numberOfDefoci();
// com
// cbed
unsigned long num_com = 2 *
gridman->adfSliceCoords().size() *
gridman->adfXGrid().size() *
gridman->adfYGrid().size() *
gridman->energyLossGrid().size() *
fs;
if(!p.adfAverageConfigurations())
num_com *= p.numberOfConfigurations();
if(!p.adfAverageDefoci())
num_com *= p.numberOfDefoci();
if(!p.cbed())
num_cbed = 0;
if(!p.adf()) {
num_com = 0;
num_adf = 0;
}
// elements
string el_str;
unsigned int counter;
for(auto & e: p.cell()->elements()) {
el_str += " " + e->symbol();
counter = 0;
for(auto & a: p.cell()->getAtoms()) {
if(a->getElement()->symbol() == e->symbol()) {
counter++;
}
}
el_str += output::fmt(" (%d)", counter);
}
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);
op::nakedprint(" build type %s\n", BUILD_TYPE);
op::nakedprint(" build date %s\n", BUILD_DATE);
op::nakedprint("\nGeneral information\n");
op::nakedprint("-----------------------------------\n");
op::nakedprint(" Grid size: %d x %d\n", gridman->samplingX(), gridman->samplingY());
op::nakedprint(" Scan points: %d x %d\n", gridman->adfXGrid().size(), gridman->adfYGrid().size());
op::nakedprint(" Energy bins: %d\n", gridman->energyLossGrid().size());
op::nakedprint(" Max angle %.1f mrad\n", kx_max * p.wavelength());
op::nakedprint(" Angle resolution: %.3f mrad\n", gridman->kx(1) * p.wavelength() * 1000);
op::nakedprint(" Num FP configs: %d\n", p.numberOfConfigurations());
op::nakedprint(" Num Defoci: %d\n", p.numberOfDefoci());
op::nakedprint("\nCrystal\n");
op::nakedprint("-----------------------------------\n");
op::nakedprint(" Number of atoms: %d\n", p.cell()->numberOfAtoms());
op::nakedprint(" Number of slices: %d\n", gridman->numberOfSlices());
op::nakedprint(" Cell size: %.2f x %.2f x %.2f nm\n", p.cell()->sizeX(), p.cell()->sizeY(), p.cell()->sizeZ());
op::nakedprint(" Elements %s\n", el_str);
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(" Total mem/proc: %.1f GB\n",
(result_mem_per_proc + mem_per_proc) / (float) (1024 * 1024 * 1024));
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("\n\n");
}
}
#endif // INITIALIZE_HPP_
\ No newline at end of file
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