Developer's Guide

Links

ChaNGa Coding Guidelines

Table of Contents

Tutorials
General notes

Tutorials

Add runtime parameters (as in a .param file)

By convention, (Hungarian Notation) the parameter name should begin with it's type, so start with a d for doubles, i for ints, etc. Let's add a double dDustParam1

In Main::Main(CkArgMsg* m) in ParallelGravity.cpp add:
```
param.dDustParam1 = 0.0;  // Default value
prmAddParam(prm,"dDustParam1",paramDouble,&param.dDustParam1,
            sizeof(double), "dust1", "<Param 1 for dust> = 0.0");
```
Here "dust1" is the short name for command line flags, "< Param 1 for dust>" is the description and " = 0.0" is the default value. This parses the command line arguments and .param file for dDustParam1. Note that there are a bunch of statements like this, remember where you place this.
Now to declare this new param. In parameters.h in the parameters struct declaration (typedef struct parameters {) add a declaration for dDustParam1:
```
    double dDustParam1;
```
Try to place it more or less in the same order you placed in in step (1)
Make the parameter communicatable (make it PUPable). In the same parameters struct declaration there is a method:
```
    inline void operator|(PUP::er &p, Parameters &param) {
```
In this method, add the following line (in the same order again): This param can be accessed in most places by param.dDustParam1. To set the value of the this runtime to say 1.5: In a .param file add dDustParam1 = 1.5. As command line argument -dust1 1.5

Add a runtime parameters class

If you are adding several parameters, it may be better to make class that will handle your runtime parameters. This helps maintain a clean, easy to read code. Here we present an example of creating dustygas parameters. Note, this is also implemented for e.g. star formation in starform.h in the class Stfm.

Let's say we have some classes, methods, etc... implemented in the header and source files dustygas.h and dustygas.cpp and I want to have some runtime parameters dDustParam1, dDustParam2, etc...

Make a class to store dustygas parameters. For example:

/* DECLARATION in dustygas.h */
class DustParam {
public:
    double dDustParam1;
    double dDustParam2;
    
    void AddParams(PRM prm);
    inline void pup(PUP::er &p);
};
// Make these attributes communicatable over the charm++ framework
inline void DustParam::pup(PUP::er &p) {
    p|dDustParam1;
    p|dDustParam2;
}

Now we have to create a method to handle parsing/storing the dust parameters.

/* DEFINITION in dustygas.cpp */
void DustParam::AddParams(PRM prm){
    dDustParam1 = 0.0;
    // "dDustParam1" is the .param file flag name, "dust1" is the command-line
    // short hand.
    prmAddParam(prm, "dDustParam1", paramDouble, &dDustParam1, 
                sizeof(double), "dust1", "<Description of DustParam1> = 0.0");
    dDustParam2 = 0.0;
    prmAddParam(prm, "dDustParam2", paramDouble, &dDustParam2, 
                sizeof(double), "dust2", "<Description of DustParam2> = 0.0");
}

Now this has to be added to the parameters struct. In parameters.h add a DustParam object and make it PUPable.

typedef struct parameters {
    ...
    DustParam dustparams;
    ...

...
inline void operator|(PUP::er &p, Parameters &param) {
    ...
    p|param.dustparams;
    ...

To make ChaNGa parse the new parameters at runtime, in ParallelGravity.cpp in Main add:

Main::Main(CkArgMsg* m) {
    ...
    param.dustparams.AddParams(prm);
    ...

Add compile flags

Goal: add a compiler flag to toggle #defining DUSTYGAS. At compile time, we'll have the option to do:

./configure --enable-dustygas=yes

Then the macro DUSTYGAS will be defined. Any code such as this will be included:

#ifdef DUSTYGAS
//do stuff...
#endif

Steps

Make sure you have autoconf installed.
In configure.ac add this:

# DUSTYGAS
ARG_ENABLE([dustygas], [Enable dusty gas], [FLAG_DUSTYGAS], [-DDUSTYGAS], [no])

The arguments to ARG_ENABLE are: the name of the option, a description, the Makefile variable, the variable's value if it's enabled, and the default state. Make sure to place this in the right section (where all the other such definitions are placed).

Run autoconf to generate the new configure script by calling autoconf with no arguments.

autoconf

In the section of Makefile.in indicated by Helper Variables, edit the defines assigment to include @FLAG_DUSTYGAS@:

defines := $(strip ... @FLAG_DUSTYGAS@ ...)

These will add the DUSTYGAS flag to the automatically generated Makefile

Good to go! Now re-configure/compile ChaNGa:

./configure --enable-dustygas=yes
make

Add SPH quantities

Let's say we want to add an attribute x to SPH (ie gas) particles. All this is handled in GravityParticle.h This requires modifying the extraSPHData, GravityParticle, and ExternalSmoothParticle classes

Steps

In class extraSPHData, add attributes (with leading underscore) and make them private and make a public method to access the attribute.
Example:
```
private:
    double _x;
public:
    inline double& x() {return _x;}
```
To allow attributes to be communicated through the charm++ framework, make them PUPable. In the class' pup method [void pup(PUP::er &p)], add a line:
```
p | _x;
```
Add an inline method to the GravityParticle class for the attribute to make it accessible from a gravity particle.
For example:
```
inline double& x() { IMAGAS; return (((extraSPHData*)extraData)->x());}
```
IMAGAS is a debugging macro. IF it's defined, then an error is raised if you try to access this method for non-gas particles
If you want the attribute to be shared w/ neighbors, it'll need to be added to the ExternalSmoothParticle class. This class is initialized with a gravity particle to share data across processors to other particles. To make an attribute shareable:
1. Add it as a public attribute
2. Add a line to the constructor to set its value to p->x(). Do this within if(TYPETest(p, TYPE_GAS)){:
```
x = p->x();
```
  Now this attribute can be accessed from other processors
3. Access data from other processors on this one, add a line to the inline getParticle( ) method which stores the attribute in a gravity particle:
```
tmp->x() = x;
```
4. Then make sure x is PUPable by adding to void pup(PUP::er &p) { The following:
```
p | x;
```

Pseudocode

/* Add these lines in --GravityParticle.h-- */

class extraSPHData 
{
 private:
    // Add this line:
    double _x;
    ...
public:
    ...
    // Add this line:
    inline double& x() {return _x;}
    ...
    void pup(PUP::er &p) {
        ...
        // Add this line:
        p | _x;
        ...
    }
}
.
.
.
class GravityParticle : public ExternalGravityParticle {
...
    // Access SPH quantities
	/// @brief Get quantities needed for SPH smooths.
	...
	// Add this line:
    inline double& x() { IMAGAS; return (((extraSPHData*)extraData)->x());}
...
}
.
.
.
class ExternalSmoothParticle {
public:
    // Add this line:
    double x;
    ...
    ExternalSmoothParticle(GravityParticle *p) 
      {
      ...
      if(TYPETest(p, TYPE_GAS)) {
        ...
        // Add this line:
        x = p->x();
        ...}
      }
    ...
    inline void getParticle(GravityParticle *tmp) const {
    ...
        if(TYPETest(tmp, TYPE_GAS)) {
            ...
            // Add this line:
            tmp->x() = x;
            ...
        }
    ...
    void pup(PUP::er &p) {
        ...
        // Add this line:
        p | x;
        ...
    }

Write out an auxiliary array

Let's follow the example of writing out the sound speed, cs.

Create a class which handles getting our values to output. These classes are defined and declared in InOutput.h. The classes are children of OutputParams. See CsOutputParams as an example. Note that currently array reading isn't quite handled correctly in CsOutputParams. If you want to read in an array that is specific to a particle type (e.g., gas) you should add a test for the particle type to setDValue:
```
void DustArrayOutputParams::setDValue(GravityParticle *p, double val) {
    if (p->isGas()) {
        p->dustFracDot() = val;
    }
```
Create the Charm++ interface (needed to handle parallelization). In ParallelGravity.ci add a PUPable statement, eg: PUPable CsOutputParams;
Snapshots are output by Main::writeOutput( ) in ParallelGravity.cpp. To have this write out the array, 3 lines need to be added. This is a long function, finding where to place these lines is hard. So just follow the example of cs and create the following 3 lines:
1. Create an object which can write out cs
```
CsOutputParams pCSOut(achFile, param.iBinaryOut, dOutTime);
```
2. Do binary output (N-chilada) if required
```
if(param.bDoCSound)
    outputBinary(pCSOut, param.bParaWrite, CkCallbackResumeThread());
```
3. Otherwise we'll do ASCII (tipsy) output
```
if(param.bDoCSound)
    treeProxy[0].outputASCII(pCSOut, param.bParaWrite,
                  CkCallbackResumeThread());
```
To write out the array when nSteps=0 for debugging, in ParallelGravity.cpp in the end of doSimulation() you can add a write-out statement after:

writeOutput(0);
      if(param.bDoGas) {
	  ckout << "Outputting gas properties ...";
	  if(printBinaryAcc)
	      CkAssert(0);
	  else {
              ...

Read in an auxiliary tipsy array

Let's follow the example of reading in the sound speed.

Make sure you've created the output class to handle io for this array (see directions in Write out an array.
For this example, the class is CsOutputParams
Reading tispy arrays is handled by treeProxy.readTipsyArray( ) (e.g. Main::initCooling() in Sph.C). Read the array by passing an output params object to readTipsyArray:
```
CsOutputParams pCsOut(param.achInFile, 0, dTime);
treeProxy.readTipsyArray(pCsOut, CkCallbackResumeThread());
```
Here we've set _iBinaryOut=0 (we're reading ASCII, not binary). This will read in the dust fraction array corresponding to param.achInFile

If you want you can check the file exists:

string csfilename=std::string(param.achInFile) + "." + pCsOut.sTipsyExt;
if(arrayFileExists(csfilename, nTotalParticles)){
....
}

NOTE: The TipsyArray reading needs to be called before any force calculation, because force calculations will change the particle order. In general this means that any routine to read in arrays must be called from Main::setupICs() in ParallelGravity.cpp.

SPH quantity time derivatives.

Here are a few notes on implementing time derivatives for SPH quantites. This assumes you are calculating a quantity which requires a loop over neighbor particles and assumes it can be done at the same time as the momentum equation loop. Let's say we have a quantity x and we have a functional form for its time derivative xDot. Here I assume the time derivative can get calculated alongside the SPH pressure terms (i.e. the momentum equation)

Add xDot and x attributes to SPH particles. (see howto:add-SPH-attributes)
The calculation of the time derivatives is handled in Sph.C in the function PressureSmoothParams::fcnSmooth( ), which also uses SphPressureTerms.h.
The time derivative needs to be zeroed at the beginning of each time step. This is handled in PressureSmoothParams::initSmoothParticle( ) which initializes smooth particle quantities whenever fcnSmooth( ) is called. Within the if statement if (p->rung >= activeRung) { add:
```
p->xDot() = 0.0;
```
Handle sharing across processors, e.g. if we have particle p on proc 1 and we want to interactive with proc 2 particles and update p->xDot() on proc 2.
Here's how this is done:
- proc 2 receives a copy of p and sets pcopy->xDot() to zero (to avoid double counting). In PressureSmoothParams::initSmoothCache( ), add:
```
if (p->rung >= activeRung) {
    ...
    p->xDot( ) =0.0;
    ...
    }
```
- p interacts with particles on proc 2 by adding terms to pcopy->xdot()
- pcopy is sent back to proc 1 and combined with p. Here we do something like p->xDot() += pcopy->xDot(). In PressureSmoothParams::combSmoothCache( ) add:
```
if (p1->rung >= activeRung) {
    ...
    p1->xDot() += p2->xDot;
    ...
}
```

Integrate SPH quantities

This assumes the quantity has been added to the SPH particles (see Add sph quantites) and that we can calculate the time derivative (see SPH quantity time derivatives).

Let's call our quantity x. Integration is handled by kick-drift. If we know the x time derivative xDot, we'll have to implement 3-steps for the integration. I'm assuming that the time step condition is larger than the courant time.
Otherwise we might need a higher order integrator of some kind. We'll implement this out of order a bit.

We'll start implementing the kick close. This gets things ready for the the next kick (ie the next timestep). In TreePiece::kick( ), in TreePiece.cpp, after if(bClosing) { the kick closing begins. Here we need to update x and xPred:
```
p->x() += p->xDot()*dDelta[p->rung];
p->xPred() = p->x();
```
x has been integrated a time step now.
We now implement open kick. After the if(bClosing){...} statement is an else { .... Store current x value in xPred and update x:
```
p->xPred() = p->x();
p->x() += p->xDot()*dDelta[p->rung]
```
Now we'll implement the drifts. Drifts are handled in TreePiece::drift( ) in TreePiece.cpp. Here we just update xPred:
```
p->xPred() += p->xDot()*dDelta
```

Use the smooth framework for neighbors

Operations involving interactions with the nearest neighbors of a particle (e.g. SPH are done using the Smooth framework. To implement a new smooth operation, first create a new class derived from the SmoothParams abstract class defined in smoothparams.h. The key method of this class is fcnSmooth(), which will be called for each particle to be smoothed, and is passed a list of neighbors. The neighbor information is contained in an array of pqSmoothNodes which contain the distance squared (fKey), the vector displacement of the particle (dx, which may be different than px - qx because of periodic boundary conditions), and a pointer to the neighbor particle data.

The next important method is isSmoothActive() which should return "1" if the particle is to be operated on by fcnSmooth() and "0" otherwise.

The methods initSmoothPartice() and initTreeParticle() can be used to initialized particles which, respectively, are to be operated on, or in the search tree, but not to be operated on. The postTreeParicle() can be used to process the particles after all smooth operations are complete.

The smooth framework handles neighbors from other processors via a software cache. By default this cache is read-only, but associative/communicative operations can be done on the neighbor data. Any modification of neighbor data by a remote processor can be combined with the data on the home processor using the combSmoothCache() method. When particle data is first copied to a remote processor, it is initialized by the initSmoothCache() method. This method can be use to avoid double counting by initializing any particle attributes used as accumulators.

Any parameters needed for the smooth operation can be included as attributes in the SmoothParams derived class, and set using the class constructor. Like all attributes that are communicated across processors, any such attributes need to be included in the pup() method of the class using the p|variable idiom. The pup() method also needs to call the method of the base class, SmoothParams::pup(). The derived class also needs to have a CkMigrateMessage constructor which calls the base constructor. The class will also need to have the declaration

PUPable_decl(XXXSmoothParams);

and a declaration needs to be added to ParallelGravity.ci:

PUPable XXXSmoothParams;

See the SmoothParams class for reference.

To use the class, construct an instance of the class, and pass a reference to the instance as the first argument to either treeProxy.startSmooth() (to smooth over the nSmooth nearest neighbors) or treeProxy.startReSmooth() (to perform a smooth over all particles within a distance specified by a particle's fBall attribute. treeProxy.startSmooth() sets the particle's fBall attribute based on the distance to the nSmooth'th neighbor.

Loop over all particle data

The particle data is contained in the TreePiece objects in the myParticles attribute. myNumParticles attribute is the number of particles on a given TreePiece. The myParticles array is actually myNumParticles+2 in size: the first and last particles in the array indicate domain boundaries. Therefore, a loop over all particles on a TreePiece will be:

    for(int i = 0; i < myNumParticles; ++i) {
        GravityParticle *p = &myParticles[i+1];
        /// manipulate p
    }

Getting all TreePieces to perform this loop requires creating a remotely callable method (called an "entry method" in Charm++) in the TreePiece class. As well as adding the method to the TreePiece class declaration in ParallelGravity.h, an entry statement needs to be added within the

    array [1D] TreePiece {

declaration in ParallelGravity.ci. This method can only have a return type of void. This method can then be called on all TreePieces from within a Main method using the syntax:

    treeProxy.newMethod(arg1, arg2);

where newMethod() is the new method which takes arguments arg1 and arg2.

The remote entry method call executes asynchronously. I.e., the method call in Main will return before the method has finished on all the TreePieces. If you need to detect when the loop has finished on all TreePieces, a Callback object needs to be created and passed to the entry method. The most straightforward Callback to use is CkCallbackResumeThread(), which will suspend execution of the calling method until the callback has been executed. For example, if we have the following in the TreePiece declaration:

    newMethod(const CkCallback& cb) {
           // Do loop over particles
           contribute(cb);
           }

This can be called from a Main method with:

    treeProxy.newMethod(CkCallbackResumeThread());

which will not return until all TreePieces have performed the loop and called contribute(). The contribute call can also be used to do a parallel reduction. See the charm++ manual for how to use reductions.

Add a source file

Let's say we've created a new source file srcfile.cpp, possibly with a header srcfile.h. To get these to compile and be linked to the rest of the source we need to add them to our objects list and our source list in the Makefile. Open Makefile.in and append to srcfile.o to the end of the objects and srcfile.c to the end of the source files:

OBJECTS = Reductions.o DataManager.o TreePiece.o ...
          ...
          ... srcfile.o \
          @OBJECTS_COOLING@

SRC = Reductions.cpp DataManager.cpp ...
        ...
        ... srcfile.cpp

Running configure will create an updated Makefile. Your source files will now be included. You may want to run make clean before re-compiling.

General notes

Vectors

As a convenience, ChaNGa uses the Vector3D class (implemented in utility) to handle vector quantities such as acceleration, velocity, and position. These are implemented with operator overloading to allowing some nice operations, such as adding vectors or multiplying by scalars.

Implementing new modules

Modules in ChaNGa are turned on/off at compile time. We want to avoid adding tons of unnecessary data to particles. For instance, dark matter only simulations do not need dust physics. Therefore, if I implement new dust physics I need to add #ifdef statements to make sure particles only have dust attributes added if I have dust physics turned on.

Debugging

Turn off parallel processing by:
1. Use one process. With charmrun, add a +p 1 flag.
2. Use one tree piece. With ChaNGa, add -p 1 flag.
  For example: `charmrun +p 1 /path/to/ChaNGa -p 1 paramfile.param``
Add gdb debugging. In the Makefile add -g to the CXXFLAGS then run make. e.g., change the line CXXFLAGS += $(OPTS) ... to CXXFLAGS += -g $(OPTS) .... Note that re-running ./configure will re-make the Makefile.
Abort on floating point exceptions. In ParallelGravity.cpp uncomment #include <fenv.h> and uncomment line calling feenableexcept(). This will catch floating point errors by throwing an error when something calculated with a floating point error is Used. e.g.:
```
double a = 2.0;
double b;
b = a/0.0;
b = a * b;
```
will throw an error at b = a * b; But this is dangerous and may break other stuff, in particular vector instructions.
To print out quantities stored in a particle's extraData struct, you just need to cast it as the appropriate pointer type. For example, a gas particle uses extraSPHData for that struct, so to print all of the additional gas properties, you simply use the gdb command

p (extraSphData)(p->extraData)

Compatibility

Try to be C++03 compatible. Many high performance computing clusters use very old compilers and may not be C++11 compatible. Don't use lambda functions, be careful with templates, don't use auto pointers.

Add/Update wiki table of contents

To help keep this wiki up to date, new sections will require adding or updating a table of contents. This can be done easily with DocToc. Follow these steps:

npm install -g doctoc
git clone https://github.com/N-BodyShop/changa.wiki.git
cd changa.wiki
doctoc <filename>

This will install doctoc, clone this wiki, and create/update a table of contents in <filename>. If you are creating a new table of contents, this will create a line:

**Table of Contents**  *generated with [DocToc](https://github.com/thlorenz/doctoc)*

Remove the *generated with [DocToc](https://github.com/thlorenz/doctoc)*.

You can commit and push the updated file if you have permissions. Otherwise, you can copy/paste text to the GitHub editor and overwrite the text.

Quick Start

Flatiron Quickstart

Flatiron CFD‐summer school‐Moving‐mesh MANGA Quickstart

Obtaining ChaNGa

Building ChaNGa

Machine Specific Build Instructions

Running ChaNGa

More Details

Creating Initial Conditions

For Developers

Developer's Guide

Coding Guidelines

Documentation Guidelines

ChaNGa Benchmarks

Performance Analysis

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Developer's Guide

Tutorials

Add runtime parameters (as in a .param file)

Add a runtime parameters class

Add compile flags

Add SPH quantities

Write out an auxiliary array

Read in an auxiliary tipsy array

SPH quantity time derivatives.

Integrate SPH quantities

Use the smooth framework for neighbors

Loop over all particle data

Add a source file

General notes

Vectors

Implementing new modules

Debugging

Compatibility

Add/Update wiki table of contents

Quick Start

More Details

For Developers

Clone this wiki locally