Commit 1b026e05 authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

still doing the refactoring. THIS IS NOT A USABLE COMMIT!

parent 898310b0
......@@ -18,7 +18,9 @@ set(PACKAGE_AUTHOR_EMAIL "jan.oliver.oelerich@physik.uni-marburg.de")
project(STEMsalabim CXX)
cmake_minimum_required(VERSION 3.3)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake")
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
string(TIMESTAMP DATE "%Y-%m-%dT%H:%M:%S")
......
......@@ -64,8 +64,6 @@ add_executable(ssb-mkin ssb-mkin.cpp)
# linking ################################################
target_compile_features(stemsalabim PRIVATE cxx_range_for)
target_link_libraries(stemsalabim_lib m ${LIBS})
target_link_libraries(stemsalabim stemsalabim_lib)
target_link_libraries(ssb-chk stemsalabim_lib)
......
......@@ -32,13 +32,13 @@ using namespace std;
using namespace stemsalabim;
void stemsalabim::Crystal::init(const std::string &crystal_file_content) {
void stemsalabim::Crystal::initFromXYZ(const std::string &field_file_content) {
// this order is very important!
atomic::ElementProvider p = atomic::ElementProvider::getInstance();
// some useful variable declarations
istringstream iss(crystal_file_content);
istringstream iss(field_file_content);
string line = "";
vector<string> tokens;
double x, y, z, msd;
......
......@@ -56,7 +56,7 @@ namespace stemsalabim {
* initialized.
* @param field_file_content string content of the crystal XYZ file.
*/
void init(const std::string &field_file_content);
void initFromXYZ(const std::string &field_file_content);
void readFieldFile(const std::string &field_file_content);
......
......@@ -28,27 +28,30 @@ using namespace std;
using namespace stemsalabim;
void FPConfManager::init(const std::vector<double> &random_numbers) {
void FPConfManager::init() {
// atomic displacements are cached for all configurations
bool fixed_slicing = Params::getInstance().isFixedSlicing();
Params &p = Params::getInstance();
bool fixed_slicing = p.isFixedSlicing();
normal_distribution<double> distribution(0.0, 1.0);
auto rgen = p.getRandomGenerator();
size_t cnt = 0;
_displacements.resize(_total_configurations);
assert(random_numbers.size() >= _total_configurations * _crystal->numberOfAtoms() * 3);
for(size_t j = 0; j < _total_configurations; j++) {
_displacements[j].resize(_crystal->numberOfAtoms());
for(size_t i = 0; i < _crystal->numberOfAtoms(); i++) {
get<0>(_displacements[j][i]) = sqrt(_crystal->getAtoms()[i]->getMSD()) * random_numbers[cnt++];
get<1>(_displacements[j][i]) = sqrt(_crystal->getAtoms()[i]->getMSD()) * random_numbers[cnt++];
get<0>(_displacements[j][i]) = sqrt(_crystal->getAtoms()[i]->getMSD()) * distribution(rgen);
get<1>(_displacements[j][i]) = sqrt(_crystal->getAtoms()[i]->getMSD()) * distribution(rgen);
if(fixed_slicing) {
get<2>(_displacements[j][i]) = 0;
cnt++;
} else {
get<2>(_displacements[j][i]) = sqrt(_crystal->getAtoms()[i]->getMSD()) * random_numbers[cnt++];
get<2>(_displacements[j][i]) = sqrt(_crystal->getAtoms()[i]->getMSD()) * distribution(rgen);
}
}
}
......
......@@ -52,7 +52,7 @@ namespace stemsalabim {
* square displacement parameters.
* @param random_numbers A vector of sufficiently many random numbers to generate displacements.
*/
void init(const std::vector<double> &random_numbers);
void init();
void initFromMSDs(const std::vector<float> &msds);
......
......@@ -254,8 +254,6 @@ void IO::writeParams(const shared_ptr<GridManager> &gridman) {
pg.att("program_arguments", p.cliArguments());
pg.att("config_file_contents", prms_string);
pg.att("work_package_size", p.workPackageSize());
pg.att("num_threads", p.numberOfThreads());
}
{
......
......@@ -26,135 +26,13 @@
#include <tuple>
#include "../utilities/netcdf_utils.hpp"
#ifdef __GNUC__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
#endif
#include "../3rdparty/args/args.hxx"
#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif
#include "../3rdparty/sole/sole.hpp"
#include "../utilities/output.hpp"
#include "../utilities/algorithms.hpp"
using namespace libconfig;
using namespace args;
using namespace std;
using namespace stemsalabim;
void Params::initFromCLI(int argc, const char **argv) {
ArgumentParser parser("STEMsalabim usage instructions.", "Please refer to the User documentation.");
HelpFlag help(parser, "help", "Display this help menu", {'h', "help"});
Flag skip_simulation(parser, "skip-simulation", "Skip simulation, only write crystal", {"skip-simulation"});
ValueFlag<string> params(parser, "params", "The parameter file path.", {'p', "params"});
ValueFlag<unsigned int> num_threads(parser, "num-threads", "The number of threads per MPI proc.", {"num-threads"});
ValueFlag<unsigned int> package_size(parser,
"package-size",
"Number of tasks in one MPI work package.",
{"package-size"});
ValueFlag<float> defocus(parser, "defocus", "The probe defocus.", {'d', "defocus"});
ValueFlag<unsigned int> num_configurations(parser,
"num-configurations",
"The number of Frozen Phonon configurations.",
{"num-configurations"});
ValueFlag<string> output_file(parser, "output-file", "The output file path.", {'o', "output-file"});
ValueFlag<string> tmp_dir(parser, "tmp-dir", "The path to the temporary directory.", {"tmp-dir"});
ValueFlag<string> crystal_file(parser, "crystal-file", "The crystal file path.", {'c', "crystal-file"});
ValueFlag<string> title(parser, "title", "The simulation title.", {'t', "title"});
ValueFlag<string> diag_dir(parser,
"diagnostics-dir",
"The directory to place diagnostics output files.",
{"diagnostics-dir"});
Flag stored_potentials(parser, "stored-potentials", "Work with stored potentials", {"stored-potentials"});
try {
parser.ParseCLI(argc, argv);
} catch(Help &) {
cout << parser;
exit(0);
} catch(ParseError &e) {
cerr << e.what() << endl;
cerr << parser;
exit(1);
}
if(argc > 1) {
vector<string> all_args;
all_args.assign(argv + 1, argv + argc);
_command_line_arguments = algorithms::join(all_args, " ");
}
if(skip_simulation) {
_cli_options.push_back(params.Name());
_skip_simulation = true;
}
if(stored_potentials) {
_cli_options.push_back(params.Name());
_stored_potentials = true;
}
if(params) {
_cli_options.push_back(params.Name());
_param_file = get(params);
} else {
output::error("a parameter file MUST be specified!\n");
}
if(num_threads) {
_cli_options.push_back(num_threads.Name());
_num_threads = get(num_threads);
}
if(package_size) {
_cli_options.push_back(package_size.Name());
_work_package_size = get(package_size);
}
if(num_configurations) {
_cli_options.push_back(num_configurations.Name());
_number_configurations = get(num_configurations);
}
if(output_file) {
_cli_options.push_back(output_file.Name());
_output_file = get(output_file);
}
if(tmp_dir) {
_cli_options.push_back(tmp_dir.Name());
_tmp_dir = get(tmp_dir);
}
if(diag_dir) {
_cli_options.push_back(diag_dir.Name());
_diagnostics_dir = get(diag_dir);
}
if(crystal_file) {
_cli_options.push_back(crystal_file.Name());
_crystal_file = get(crystal_file);
}
if(title) {
_cli_options.push_back(title.Name());
_title = get(title);
}
if(defocus) {
_cli_options.push_back(defocus.Name());
_probe_defocus = get(defocus);
}
_uuid = sole::uuid4().str();
}
void Params::readParamsFromString(const string &prms) {
......@@ -171,7 +49,6 @@ void Params::readParamsFromString(const string &prms) {
//Setting &root = _cfg.getRoot();
_cfg.lookupValue("application.random_seed", _random_seed);
_cfg.lookupValue("application.skip_simulation", _skip_simulation);
string prec_str = "single";
_cfg.lookupValue("application.precision", prec_str);
......@@ -312,8 +189,9 @@ void Params::broadcast(int rank) {
mpi_env.broadcast(_init_from_ncfile, rank);
mpi_env.broadcast(_num_threads, rank);
mpi_env.broadcast(_work_package_size, rank);
mpi_env.broadcast(_random_seed, rank);
mpi_env.broadcast(_skip_simulation, rank);
mpi_env.broadcast(_title, rank);
......@@ -398,7 +276,6 @@ void Params::readParamsFromNCFile() {
{
auto gg = g.group("application");
_random_seed = gg.att<unsigned int>("random_seed");
_skip_simulation = gg.is("skip_simulation");
}
{
......
......@@ -324,31 +324,6 @@ namespace stemsalabim {
return _normalize_always;
}
/*!
* Return whether diagnostic output should be printed.
* @return true if diagnostics
*/
bool printDiagnostics() const {
return !_diagnostics_dir.empty();
}
/*!
* Get directory of diagnostics output.
* @return directory as string
*/
std::string diagnosticsDir() const {
return _diagnostics_dir;
}
/*!
* Return whether the multi-slice simulation is skipped, i.e., a dry-run is performed
* writing only crystal coordinates to the output.
* @return true if skip simulation
*/
bool skipSimulation() const {
return _skip_simulation;
}
/*!
* Return the package of each thread of slave MPI processes,
* i.e., how many pixels this thread calculates before reporting the results back to the master.
......@@ -570,8 +545,6 @@ namespace stemsalabim {
std::string _uuid{""};
std::string _diagnostics_dir{""};
bool _skip_simulation{false};
double _sample_density{36.0};
bool _output_compress{false};
......
......@@ -64,7 +64,7 @@ void Simulation::init() {
crystal_file_content = buffer.str();
}
mpi_env.broadcast(crystal_file_content, mpi::Environment::MASTER);
_crystal->init(crystal_file_content);
_crystal->initFromXYZ(crystal_file_content);
}
// read and communicate the field file content, when it is given as a parameter
......@@ -103,7 +103,7 @@ void Simulation::init() {
}
mpi_env.broadcast(numbers, mpi::Environment::MASTER);
_fpman->init(numbers);
//_fpman->init(numbers);
} else {
vector<float> msds;
if(mpi_env.isMaster()) {
......@@ -142,9 +142,6 @@ void Simulation::init() {
_gridman->adfXGrid().size(),
_gridman->adfYGrid().size()));
if(mpi_env.isMpi())
PRINT_DIAGNOSTICS("starting initialization of FFTW wisdom\n");
// every MPI processor figures out the FFT plans itself. This should somehow be
// improved at some point
//Wave::initPlans(_gridman->samplingX(), _gridman->samplingY());
......@@ -152,9 +149,6 @@ void Simulation::init() {
// make the atomic potentials
//Scattering::init(_crystal);
if(mpi_env.isMpi())
PRINT_DIAGNOSTICS("finished initialization\n");
printMaster(output::fmt("Initialization took %s.\n", output::humantime(algorithms::getTimeSince(init_start))));
initBuffers();
......@@ -222,18 +216,10 @@ void Simulation::run() {
if(mpi_env.isMaster())
_io->writeCrystal(_fpman, _crystal);
// skip multi-slice simulation when skip_simulation = true in the configuration file.
if(prms.skipSimulation()) {
continue;
}
// generate transmission functions in parallel using both threads and the main thread.
// This is executed on each MPI processor.
t1 = high_resolution_clock::now();
if(mpi_env.isMpi())
PRINT_DIAGNOSTICS(output::fmt("starting slice work\n"));
TaskQueue slice_work;
slice_work.append(_crystal->slices());
......@@ -255,9 +241,6 @@ void Simulation::run() {
slice_work.finish(slic_index);
}
if(mpi_env.isMpi())
PRINT_DIAGNOSTICS("waiting for slice work to finish work\n");
slice_work.waitUntilFinished();
for(auto & t : threads)
t.join();
......@@ -693,8 +676,6 @@ void Simulation::multisliceWorker(const SimulationState &st) {
}
PRINT_DIAGNOSTICS("Exiting loop");
printLine(st,
output::fmt("Finished working for %s and calculated %d pixels.\n",
output::humantime(algorithms::getTimeSince(t1)),
......
......@@ -27,6 +27,10 @@
using namespace std;
using namespace stemsalabim;
void Params::initFromCLI(int argc, const char **argv) {
}
int main(int argc, const char **argv) {
Params &p = Params::getInstance();
......
......@@ -21,94 +21,124 @@
#include <iostream>
#include <thread>
#include <algorithm>
#ifdef __GNUC__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
#endif
#include "3rdparty/args/args.hxx"
#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif
#include "classes/Crystal.hpp"
#include "classes/Simulation.hpp"
using namespace std;
using namespace stemsalabim;
#include "3rdparty/sole/sole.hpp"
#include "utilities/timings.hpp"
using namespace std;
using namespace chrono;
using namespace stemsalabim;
void Params::initFromCLI(int argc, const char **argv) {
using namespace args;
ArgumentParser parser("STEMsalabim usage instructions.", "Please refer to the User documentation.");
HelpFlag help(parser, "help", "Display this help menu", {'h', "help"});
ValueFlag<string> params(parser, "params", "The parameter file path.", {'p', "params"});
ValueFlag<unsigned int> num_threads(parser, "num-threads", "The number of threads per MPI proc.", {"num-threads"});
ValueFlag<float> defocus(parser, "defocus", "The probe defocus.", {'d', "defocus"});
ValueFlag<string> output_file(parser, "output-file", "The output file path.", {'o', "output-file"});
ValueFlag<string> crystal_file(parser, "crystal-file", "The crystal file path.", {'c', "crystal-file"});
Flag stored_potentials(parser, "stored-potentials", "Work with stored potentials", {"stored-potentials"});
try {
parser.ParseCLI(argc, argv);
} catch(Help &) {
cout << parser;
exit(0);
} catch(ParseError &e) {
cerr << e.what() << endl;
cerr << parser;
exit(1);
}
if(argc > 1) {
vector<string> all_args;
all_args.assign(argv + 1, argv + argc);
_command_line_arguments = algorithms::join(all_args, " ");
}
if(stored_potentials) {
_cli_options.push_back(params.Name());
_stored_potentials = true;
}
if(params) {
_cli_options.push_back(params.Name());
_param_file = get(params);
} else {
output::error("a parameter file MUST be specified!\n");
}
if(num_threads) {
_cli_options.push_back(num_threads.Name());
_num_threads = get(num_threads);
}
if(output_file) {
_cli_options.push_back(output_file.Name());
_output_file = get(output_file);
}
if(crystal_file) {
_cli_options.push_back(crystal_file.Name());
_crystal_file = get(crystal_file);
}
if(defocus) {
_cli_options.push_back(defocus.Name());
_probe_defocus = get(defocus);
}
_uuid = sole::uuid4().str();
}
int main(int argc, const char **argv) {
auto &mpi_env = mpi::Environment::getInstance();
// Initialization of the input NC file is done only by the master
// MPI proc, in case the tool is run differently.
if(mpi_env.isSlave()) {
return 0;
}
Params &p = Params::getInstance();
p.initFromCLI(argc, argv);
std::ifstream t(p.paramsFileName());
string params_file_content = std::string((std::istreambuf_iterator<char>(t)), std::istreambuf_iterator<char>());
p.readParamsFromString(params_file_content);
high_resolution_clock::time_point init_start = high_resolution_clock::now();
p.readParamsFromString(algorithms::readTextFile(p.paramsFileName()));
auto _crystal = std::make_shared<Crystal>();
auto _gridman = std::make_shared<GridManager>(_crystal);
auto _fpman = std::make_shared<FPConfManager>(_crystal);
auto _io = std::make_shared<IO>();
// read and communicate the content of the crystal file and initialize crystal
// object.
{
ifstream t(p.crystalFilename());
stringstream buffer;
buffer << t.rdbuf();
_crystal->init(buffer.str());
}
// generate the required random numbers for the frozen phonon stuff.
{
const size_t number_fp_random_numbers = p.numberOfConfigurations() *
p.numberOfDefoci() *
_crystal->numberOfAtoms() *
3;
vector<double> numbers;
numbers.resize(number_fp_random_numbers);
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
normal_distribution<double> distribution(0.0, 1.0);
auto rgen = p.getRandomGenerator();
for(unsigned long i = 0; i < number_fp_random_numbers; ++i)
numbers[i] = distribution(rgen);
#pragma GCC diagnostic pop
_fpman->init(numbers);
}
_crystal->initFromXYZ(algorithms::readTextFile(p.crystalFilename()));
_fpman->init();
_gridman->generateGrids();
// initialize output file if this is master process and write some simulation information
// to it.
auto _io = std::make_shared<IO>();
// As IO should be handled by the master only, we need to initialize
// the instance only when this is the master process.
// The crystal instance will be initialized here and then broadcast to
// tge slaves.
_io->initNcFile(_gridman, _crystal);
_io->writeParams(_gridman);
_io->writeGrids(_gridman);
// in case the slicing is fixed, this can be done outside of the main loops
if(p.isFixedSlicing()) {
_fpman->assignAtomsToSlices();
for(shared_ptr<Slice> &slice: _crystal->slices())
Simulation::printMaster(output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z() / 10,
(slice->z() + slice->thickness()) / 10,
slice->numberOfAtoms()));
}
SimulationStateManager state_manager(_gridman);
vector<vector<vector<float>>> vz(_crystal->slices().size());
......@@ -121,28 +151,19 @@ int main(int argc, const char **argv) {
_fpman->setConfiguration(st.iteration());
// if a new defocus was started, do some output stuff.
if(!st.iconf()) {
Simulation::printMaster(st, output::fmt("Starting Defocus %.2fnm.\n", st.defocus() / 10));
}
Simulation::printMaster(st, "Starting configuration. \n");
// in each configuration, assign atoms to the corresponding slices.
if(!p.isFixedSlicing()) {
_fpman->assignAtomsToSlices();
for(shared_ptr<Slice> &slice: _crystal->slices())
Simulation::printMaster(st,
output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z() / 10,
(slice->z() + slice->thickness()) / 10,
slice->numberOfAtoms()));
}
_fpman->assignAtomsToSlices();
for(shared_ptr<Slice> &slice: _crystal->slices())
Simulation::printMaster(st,
output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z() / 10,
(slice->z() + slice->thickness()) / 10,
slice->numberOfAtoms()));
// in each frozen phonon iteration, first write the crystal to the output
// file.
_io->writeCrystal(_fpman, _crystal);
if(p.storedPotentials()) {
......@@ -150,29 +171,18 @@ int main(int argc, const char **argv) {
slice_work.append(_crystal->slices());
vector<thread> threads;
for(unsigned int nt = 1; nt < p.numberOfThreads(); nt++) {
for(unsigned int nt = 0; nt < p.numberOfThreads(); nt++) {
threads.emplace_back([&] {
unsigned int slic_index;
while(slice_work.pop(slic_index)) {
vz[slic_index] = _crystal->slices()[slic_index]->calculatePotential(_fpman, _gridman);
slice_work.finish(slic_index);
}
});
}
// also do work on the main thread.
unsigned int slic_index;
while(slice_work.pop(slic_index)) {
vz[slic_index] = _crystal->slices()[slic_index]->calculatePotential(_fpman, _gridman);
slice_work.finish(slic_index);
}
if(mpi_env.isMpi())
PRINT_DIAGNOSTICS("waiting for slice work to finish work\n");
slice_work.waitUntilFinished();
for(auto & t : threads)
for(auto &t : threads)
t.join();
_io->writePotentials(_fpman, vz);
......
......@@ -20,6 +20,20 @@
*/
#include <iostream>
#ifdef __GNUC__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
#endif
#include "3rdparty/args/args.hxx"
#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif
#include "3rdparty/sole/sole.hpp"
#include "utilities/Wave.hpp"
#include "classes/Crystal.hpp"
#include "classes/Simulation.hpp"
......@@ -27,6 +41,72 @@
using namespace std;
using namespace stemsalabim;
void Params::initFromCLI(int argc, const char **argv) {
using namespace args;
ArgumentParser parser("STEMsalabim usage instructions.", "Please refer to the User documentation.");
HelpFlag help(parser, "help", "Display this help menu", {'h', "help"});
ValueFlag<string> params(parser, "params", "The parameter file path.", {'p', "params"});
ValueFlag<unsigned int> num_threads(parser, "num-threads", "The number of threads per MPI proc.", {"num-threads"});
ValueFlag<unsigned int> pkgsize(parser, "package-size", "Number of tasks in a MPI work package.", {"package-size"});
ValueFlag<string> output_file(parser, "output-file", "The output file path.", {'o', "output-file"});
ValueFlag<string> tmp_dir(parser, "tmp-dir", "The path to the temporary directory.", {"tmp-dir"});
Flag stored_potentials(parser, "stored-potentials", "Work with stored potentials", {"stored-potentials"});
try {
parser.ParseCLI(argc, argv);
} catch(Help &) {
cout << parser;
exit(0);
} catch(ParseError &e) {
cerr << e.what() << endl;
cerr << parser;
exit(1);
}
if(argc > 1) {
vector<string> all_args;
all_args.assign(argv + 1, argv + argc);
_command_line_arguments = algorithms::join(all_args, " ");
}
if(stored_potentials) {
_cli_options.push_back(params.Name());
_stored_potentials = true;
}
if(params) {
_cli_options.push_back(params.Name());
_param_file = get(params);
} else {
output::error("a parameter file MUST be specified!\n");
}
if(num_threads) {
_cli_options.push_back(num_threads.Name());
_num_threads = get(num_threads);
}
if(pkgsize) {
_cli_options.push_back(pkgsize.Name());
_work_package_size = get(pkgsize);
}
if(output_file) {
_cli_options.push_back(output_file.Name());
_output_file = get(output_file);
}
if(tmp_dir) {
_cli_options.push_back(tmp_dir.Name());
_tmp_dir = get(tmp_dir);
}
_uuid = sole::uuid4().str();
}