Table of contents


Back to the main page Back to the documentation 1. Overview: FML 2. Particle classes

Overview: FML


The domain we are working with in the library is a qubic box $[0,1)^{\rm NDIM}$ of any dimension(s). It is generally assumed periodic for which the point $1$ wraps around to $0$, but several classes and algorithms lets you specify if this is the case or not. There are no units assumed in the library. Grids have the same number of gridcells in each dimension. The parallelization over MPI is a slab-based domain decomposition along the $x$-axis.

The library is divided into several namespaces. Global variables and functions are placed in the top namespace FML and are defined in Global.

We have just a few global variables (that we could just as well have computed where they are needed). The most important ones are FML::ThisTask and FML::NTasks. This is the ID of a given task and the total number of tasks. Next we have FML::xmin_domain and FML::xmax_domain. This is the $x$-domain of the grid the given task is responsible for. Lastly we have FML::NThreads which is the number of threads we have availiable (this is just OMP_NUM_THREADS). We also have some general type definitions and a few useful functions.

In Global we also have two singletons FML::MPISetup and FML::FFTWSetup that are responsible for initializing MPI and FFTW. If you don't want automatic initialization then you will have to make sure that the variables above is set correctly, see Global.cpp. To turn this off use the defines NO_AUTO_MPI_SETUP and NO_AUTO_FFTW_SETUP.

Error handling. This is hard to do properly so fuck it: if we detect an error we typically call assert_fml which gives some info about where the error happened and what went wrong and then aborts. Several of the individual classes can throw errors instead which can be useful in some cases, but see the individual classes to change this behavior.

Some simple methods to find max, min and sum of a single value from all the tasks. The last method gathers a value to all tasks. These operations are done in-place, i.e. it overwrites the value we pass in a pointer to.

  template<class T>
    void MaxOverTasks(T * value);
  template<class T>
    void MinOverTasks(T * value);
  template<class T>
    void SumOverTasks(T * value);
  template<class T>
    std::vector<T> GatherFromTasks(T * value);

In Global we define a container FML::Vector<T> which is just normal std::vector but with a custom allocator to keep track of the memory use. We currently only use Vector for FFTWGrid, MPIGrid and MPIParticles so this is all it tracks. The define MIN_BYTES_TO_LOG=XXX sets how large (in bytes) an allocation has to be to be logged and MAX_ALLOCATIONS_IN_MEMORY=YYY sets the max number of elements we keep track of before giving up and just stop logging.

To add a label to an already allocated object you can call:

FML::Vector< double > data(1000);
FML::MemoryLog::get()->add_label(data.data(), data.size(), "My vector data");

And the classes that uses this has special methods for this:

// FFTWGrid
FFTWGrid grid(Nmesh);
grid.add_memory_label("label");

// MPIPartices
MPIParticles part;
// ...create the particles
part.add_memory_label("label");

To see what we have in memory at the current time:

FML::MemoryLog::get()->print();  

Example output:

#=====================================================

We are only tracking allocation of standard container
and only allocations larger than 0 bytes

Memory in use:               8.8 MB
Peak memory use:             8.8 MB

We have the following things allocated on task 0: 
Address: 0x106453000 Size: 0.8 (MB) Label: [This is GridA]
Address: 0x106517000 Size: 8   (MB) Label: [This is GridB]

Total memory as function of time:
 Time (sec):             0 Memory:             0 (MB)
 Time (sec):       2.3e-05 Memory:           0.8 (MB)
 Time (sec):      0.001174 Memory:           8.8 (MB)

#=====================================================

Particle Classes


A lot of the algorithms inside the library is templated on particle classes. These will accept any kind of particle as long as it has the two methods

class MyParticle{
  ...

public:
  // Here PositionType can be any floating-point type (float, double, long double)
  PositionType *get_pos();

  // The dimension of the position
  int get_ndim();
}

The get_pos method returns a pointer to the position (assumes NDIM consecutive values in memory) and get_ndim gives the dimension of the position vector. If you have dynamically allocated memory inside the class (i.e. sizeof(MyParticle) does not give the actual size) then you must provide methods for how to append the particle to a buffer and assign the particle from a buffer that we use to do the communication between tasks:

class MyParticle{
  ...

public: 
  ...

  // How many bytes to communicate
  int get_particle_byte_size(){ 
  
  // Append the particle  data to a buffer
  void append_to_buffer(char *data);
  
   // Assign the particle data from a buffer
  void assign_from_buffer(char *data);
}

For example if we have a std::vector<double> stuff; inside the class then the method would look something like this:

std::vector<double> stuff; 

// We need to communiate both the data and its size (plus any other things in the class)
int get_particle_byte_size(){ 
  return stuff.size() * sizeof(double) + sizeof(int) + (...other things we have...); 
}
  
// Append the particle  data to a buffer
void append_to_buffer(char *data){
  // Write the size of stuff
  int nelements = stuff.size();
  std::memcpy(data, &nelements, sizeof(int));
  data += sizeof(int);

  // Write the data in stuff
  std::memcpy(data, stuff.data(), stuff.size() * sizeof(double));
  data += stuff.size() * sizeof(double);

  (...other things we have...)
}
  
 // Assign the particle data from a buffer
void assign_from_buffer(char *data){
  // Read size of stuff
  int nelements;
  std::memcpy(&nelements, data, sizeof(int));
  data += sizeof(int);

  // Allocate stuff
  stuff.resize(nelements);
  
  // Read stuff
  std::memcpy(stuff.data(), data, stuff.size() * sizeof(double));
  data += stuff.size() * sizeof(double);

  (...other things we have...)
}

Some of the algorithms can do more stuff if certain quantities are present in the particles and some algorithms *require* certain additional methods. If you don't have this then you will either get a compile error or a runtime error (if this cannot be deduced at compile time). For example when we do grid assignement then the algorithm checks to see if we have a get_mass method. If we do then it will use this, if not then it will assume that all particles have the same mass. If we want to transform particles from realspace to redshift-space then we *must* have a get_vel method that gives us the velocity of the particle. Here is a list of the methods that some algorithms have as either an optional or a required part:

// (Below the types we used are not fixed and can be any type you want)

// Optional: VELOCITY
double *get_vel();  // Returns the velocity of the particle

// Optional: MASS (Fiducial value 1.0 if not present)
double get_mass();           // Returns the mass of the particle
void set_mass(double _mass); // Set the mass of the particle

// Optional: VOLUME (Fiducial value 1.0 if not present)
double get_volume();             // Returns the volume of the particle
void set_volume(double _volume); // Set the mass of the particle

// Optional: ID
int get_id();         // Returns the ID of the particle
void set_id(int _id); // Set the ID of the particle

//================================================================================
// RAMSES related methods
//================================================================================
// (for example when reading a RAMSES file it checks if these are present and stores the related data if they are)

// Optional: FAMILY (only used for RAMSES files)
char get_family();             // Returns the type of the particle
void set_family(char _family); // Set the type of the particle

// Optional: TAG (only used for RAMSES files)
char get_tag();          // Returns the tag of the particle
void set_tag(char _tag); // Set the tag of the particle

// Optional: LEVEL (only used for RAMSES files)
int get_level();            // Returns the level of the particle
void set_level(int _level); // Set the level of the particle

//================================================================================
// Lagrangian perturbation theory related stuff
//================================================================================

// Optional: Lagrangian Position 
double *get_q();  // Returns the Lagrangian coordinate of the particcle

// Optional: Displacement field
double *get_D_1LPT();  // Returns the 1LPT displacement field of the particcle
double *get_D_2LPT();  // Returns the 2LPT displacement field of the particcle
double *get_D_3LPT();  // Returns the 3LPT displacement field of the particcle

//================================================================================
// Galaxy related methods (required for some galaxy related methods)
/================================================================================

// Optional: galaxy position on the sky and weight
double get_RA();         // Returns the right ascention of the galaxy
void set_RA(double _RA); // Set the right ascention of the galaxy

double get_DEC();          // Returns the declination of the galaxy
void set_DEC(doubel _DEC); // Set the declination of the galaxy

double get_z();        // Returns the redshift of the galaxy
void set_z(double _z); // Set the redshift of the galaxy

double get_weight();             // Returns the weight of the galaxy
void set_weight(double _weight); // Set the weight of the galaxy

//================================================================================
// Halo related method
//================================================================================

// Optional: Friends of Friend ID 
size_t get_fofid();            // Returns the FoFID of the particle (SIZE_MAX if it does not belong to a FoF group)
void set_fofid(size_t _FoFID); // Sets the FoFID of the particle (will be set by FoF finder if this method exists)

To get some info about what methods we have detected in your particle class you can call

FML::PARTICLE::info< MyParticle >();