-
Notifications
You must be signed in to change notification settings - Fork 189
Software design
This page documents the way ESPResSo is designed and how it interfaces with libraries.
ESPResSo uses a MPI callback mechanism based on a visitor pattern. At the start of the program, global variables from the head node are communicated to the worker nodes. The script interface sends and receives data from the core on the head node. When particles are created, the head node distributes them to the corresponding worker nodes based on the currently active domain decomposition.
Functions in the script interface that change the global state of the system in the core or accesses particles are implemented using callbacks. There are also "event" callbacks that get triggered when the global state of the system changes, e.g. the box dimensions in NpT simulations, to synchronize the global state on every MPI rank.
Below is one possible implementation of a callback to calculate the total linear momentum of the system. The momentum calculation is written in the callback calc_linear_momentum_local(). It is registered as a callback and tagged for a reduction operation using the macro REGISTER_CALLBACK_REDUCTION(). The calc_linear_momentum() function is responsible for calling the callbacks on all ranks. This function is exported in the Python interface, where it can be called by the user on the head node.
#include "Particle.hpp"
#include "cells.hpp"
#include "communication.hpp"
#include <utils/Vector.hpp>
/** Callback to calculate the linear momentum of the cell on this node. */
static Utils::Vector3d calc_linear_momentum_local() {
auto const particles = cell_structure.local_particles();
auto const momentum =
std::accumulate(particles.begin(), particles.end(), Utils::Vector3d{},
[](Utils::Vector3d &m, Particle const &p) {
return m + p.p.mass * p.m.v;
});
return momentum;
}
REGISTER_CALLBACK_REDUCTION(calc_linear_momentum_local,
std::plus<Utils::Vector3d>())
/** Calculate the total linear momentum of the system.
* Must be called on the head node.
*/
Utils::Vector3d calc_linear_momentum() {
Utils::Vector3d const result = mpi_call(::Communication::Result::reduction,
std::plus<Utils::Vector3d>(),
calc_linear_momentum_local);
return result;
}There are specialized macros and mpi_call() overloads to accommodate for a wide range of collective operations. See the doxygen documentation of header files MpiCallbacks.hpp and communication.hpp for more details.
Functions called exclusively on the head node should throw an exception. The corresponding import in the Python interface should decorate the symbol with except +, so that the C++ exception gets transformed into a Python exception.
# sd.pxd
cdef extern from "sd.hpp":
void set_sd_viscosity(double eta) except +/* sd.hpp */
void set_sd_viscosity(double eta);
/* sd.cpp */
#include <stdexcept>
#include <string>
double sd_viscosity;
void set_sd_viscosity(double eta) {
if (eta < 0.0)
throw std::runtime_error("invalid viscosity: " + std::to_string(eta));
sd_viscosity = eta;
}Functions called by MPI callbacks cannot throw exceptions. Instead they should use the runtime errors logging mechanism to store a string containing the error message. All Python functions that indirectly call such a callback must add a handle_errors() to halt the flow of the program and display the error messages, if the logger is not empty.
/* p3m.hpp */
void p3m_tune();
/* p3m.cpp */
#include "p3m.hpp"
#include "errorhandling.hpp"
#include "communication.hpp"
void p3m_tune_local() {
int error = sanity_checks();
if (error) {
runtimeErrorMsg() << "sanity checks failed";
return;
}
/* implementation not shown */
}
REGISTER_CALLBACK(p3m_tune_local)
void p3m_tune() {
mpi_call_all(p3m_tune_local);
}# electrostatics.pxd
cdef extern from "p3m.hpp":
void p3m_tune()
# electrostatics.pyx
from . cimport electrostatics # get symbols from the pxd file
from .utils cimport handle_errors
def tune():
p3m_tune()
# convert the logged error messages into a Python exception
handle_errors("P3M tuning failed")Functionality that doesn't depend on the state of the system can be extracted from the EspressoCore target and store in a submodule of src/. This is used for example by the EspressoUtils module, which is a dependency of EspressoCore. Such submodules typically adopt the Pitchfork layout, i.e.
src/
\--- library/
|--- include/
| \--- library/
| \--- feature.hpp
|--- src/
| |--- feature_impl.hpp
| \--- feature_impl.cpp
\--- tests/
\--- feature_test.cppWith this layout, you can pass the compiler flag -Isrc/library/include to make the public header file library/feature.hpp visible to consumers of that library while the implementation details in src/library/src remain private.
CMake-based projects that provide functionality to ESPResSo are included using the FetchContent method. This approach enables pinning the library version and patching the source code:
if(WITH_WALBERLA)
include(FetchContent)
FetchContent_Declare(
walberla
GIT_REPOSITORY https://i10git.cs.fau.de/walberla/walberla
GIT_TAG 64a4d68
PATCH_COMMAND patch -p0 --ignore-whitespace -i
${PROJECT_SOURCE_DIR}/cmake/waLBerlaFunctions.patch
)
FetchContent_Populate(walberla)
set(WALBERLA_BUILD_TESTS off CACHE BOOL "")
set(CMAKE_POSITION_INDEPENDENT_CODE on CACHE BOOL "")
add_subdirectory("${walberla_SOURCE_DIR}" "${walberla_BINARY_DIR}")
set(WALBERLA_LIBS walberla::core walberla::blockforest walberla::boundary walberla::field walberla::lbm)
endif(WITH_WALBERLA)Another method consists in loading the library from the system via $PKG_CONFIG_PATH and $LD_LIBRARY_PATH:
find_package(NumPy REQUIRED)
find_program(IPYTHON_EXECUTABLE NAMES jupyter ipython3 ipython)
pkg_check_modules(SCAFACOS scafacos)For better separation of concerns, the external library should not expose its public interface to the EspressoCore directly, but instead rely on a submodule such that the external library is a private dependency of the submodule, and the submodule is a private dependency of EspressoCore. The submodule should use a C++ design pattern to hide the details of the external library.
Using the Pitchfork layout outlined in section Internal libraries, one could declare an abstract base class in feature.hpp and write a concrete derived class in feature_impl.hpp.
/* ExternalFeature.hpp (external library) */
class ExternalFeature {
public:
ExternalFeature(int argument): m_argument(argument) {}
double energy(int position) const;
private:
int m_argument;
};
/* feature.hpp */
class FeatureBase {
public:
virtual double get_energy(int position) const = 0;
};
FeatureBase *new_feature(int argument);
/* feature_impl.hpp */
#include "feature.hpp"
#include <ExternalFeature.hpp>
class Feature: public FeatureBase {
public:
Feature(int argument): m_feature(ExternalFeature(argument)) {}
double get_energy(int position) const override;
private:
ExternalFeature m_feature;
};
/* feature_impl.cpp */
#include "feature.hpp"
FeatureBase *new_feature(int argument) {
FeatureBase *instance = new Feature(argument);
return instance;
}
double Feature::get_energy(int position) const {
return m_feature.energy(position);
}
/* src/core/analysis.hpp */
#include <library/feature.hpp>
FeatureBase *feature();
/* src/core/analysis.cpp */
#include "analysis.hpp"
#include "communication.hpp"
#include <library/feature.hpp>
#include <stdexcept>
namespace {
FeatureBase *feature_instance = nullptr;
}
FeatureBase *feature() {
if (!feature_instance)
throw std::runtime_error("feature is uninitialized.");
return feature_instance;
}
void init_feature_local(int argument) {
feature_instance = new_feature(argument);
}
REGISTER_CALLBACK(init_feature_local)
void init_feature(int argument) {
Communication::mpiCallbacks().call_all(init_feature_local, argument);
}
double feature_energy(int position) {
return feature()->get_energy(position);
}C++ functions and classes from the ESPResSo core are exported in the Python interface using the src/python/espressomd/*.pxd files. These exported symbols can then be used within cdef functions in src/python/espressomd/*.pyx files.
# interface.pxd
cdef extern from "cuda_init.hpp":
int cuda_get_device()
# interface.pyx
from . cimport interface # get symbols from the pxd file
cdef device(self):
"""Get device id."""
dev = cuda_get_device()
if dev == -1:
raise Exception("cuda device get error")
return devLikewise, struct and class declarations can be copy-pasted in .pxd files to allow returning a copy of a core object in a Cython function.
When exposing C++ classes in Python, it is sometimes necessary to write a bridge between the core implementation and the script interface. This is required when the Python interface needs to access the original core object through a pointer instead of working on a local copy of the core object. This is achieved by the "glue" classes in src/script_interface, which derive from AutoParameters. The class derived from AutoParameters stores a shared pointer to the core object and sets up getters and optionally setters to members of the core object.
A corresponding Cython class needs to be written and derived from the ScriptInterfaceHelper class. The Cython class needs to name the shared object symbol and explicitly list the subset of AutoParameters getters and setters that are visible to the user. Values are automatically converted from the python types (int, float, string, etc.) to the corresponding C++ types (int, double, std::string) and passed between the python interface and core via a boost::variant. The type in the variant must be explicitly given in the setters and getters from the AutoParameters in src/script_interface.
Below is an example with observable classes. The Observable parent class has a shape() method that can be called directly from an instance of the observable. The AutoParameters class has an additional calculate() method that is not meant to be used directly by the user, but can be called from within the class with call_method("calculate"). This approach is used when the result of a function needs to be post-processed, in this case to reshape the flat array into an N-dimensional array.
@script_interface_register
class Observable(ScriptInterfaceHelper):
"""
Base class for all observables.
Methods
-------
shape()
Return the shape of the observable.
"""
_so_name = "Observables::Observable"
_so_bind_methods = ("shape",)
_so_creation_policy = "LOCAL"
def calculate(self):
return np.array(self.call_method("calculate")).reshape(self.shape())
@script_interface_register
class ComVelocity(Observable):
"""Calculates the center of mass velocity for particles with given ids."""
_so_name = "Observables::ComVelocity"To enable checkpointing, Python classes need to define special methods to serialize and deserialize themselves. This is usually achieved by implementing adequate __getstate__(self) and __setstate__(self, params) methods. The result of the __getstate__ can be any data type, usually a list or dict, and will be written to a checkpoint file by the pickle module. When loading the system from a checkpoint file, pickle will read the value and pass it to the class as the params argument to the __setstate__ function to instantiate an object. During this instantiation, the __init__(self) method is not called.