Now that you have got a feel for running Cactus you must be able to write your own thorns. Below you'll encounter a very simple but instructive example of a thorn. The purpose is to solve the wave equation in three dimensions
This work has been already done when you run WaveDemo or WaveToy, the difference is that here we have made some simplifications to make the code shorter and simpler.
At this point you have read at least once the section Application Thorn Writing of the Cactus user's manual. Here you'll find just an example of the principal features of Cactus.
We use finite differences to integrate the wave equation. Initial conditions are Gaussian shape and zero velocity. Boundary conditions is taken as fixed. Thinking in one dimension, this is the well known problem of the string fixed at its ends with some initial shape and velocity. You can take a look at the code here. Below is a step-by-step explanation of the code concerning cactus issues.
#include "cctk.h" #include "cctk_Arguments.h"
Cactus header files that must be included in every source code file.
#define PI 4*atan(1)
We define our own Pi.
void wave3D_generate_ID( CCTK_ARGUMENTS )
This function sets the initial shape and velocity of the wave. Its arguments are defined in the CCTK_ARGUMENTS macro, which is expanded when the code is compiled. Every function that is scheduled to run (in schedule.ccl) needs to take these arguments.
DECLARE_CCTK_ARGUMENTS
Defines the arguments passed by cactus infrastructure.
ierr = CCTK_CoordRegisterSystem( dim, "cart3d" );
Assigns a coordinate system with a chosen name and dimension.
iend = cctk_lsh[0] - 1; jend = cctk_lsh[1] - 1; kend = cctk_lsh[2] - 1;
Depending on the number of processors we use, the work is distributed among all of them. cctk_lsh gives the local grid size on this processor. For instance, if you have a 100x100x100 grid and 8 processors the local size on each processor would be 50x50x50. We substract 1 because arrays in C start from zero.
/* initial conditions are zero velocity and gaussian shape */
for(k=kstart-1; k< kend+1; k++)
{
for(j=jstart-1; j< jend+1; j++)
{
for(i=istart-1; i< iend+1; i++)
{
index = CCTK_GFINDEX3D( cctkGH, i, j, k );
/* ID with CartGrid3D variables */
phi_p[ index ] = exp( -(x[index]*x[index] +
y[index]*y[index] + z[index]*z[index]) );
phi[ index ] = phi_p[ index ];
}
}
}
These loops assign the initial data. The macro CCTK_GFINDEX3D makes the convertion between a 3D array indices and 1D array index. This is due to the fact that arrays are created as a single 1D array, i.e. a continous memory block.
Here appears for the first time our grid function phi. A grid function is like an array. They are distributed in all the processors. We know that phi depends on both space and time, here is where phi_p comes to play. The suffix _p represents the value of phi at the previuos time level. That is phi = phi(x, y, z, t); phi_p = phi(x, y, z, t - /Delta t): phi_p_p = phi(x, y, z, t - 2 /Delta t); and so on. These are referred as timelevels which are rotated after every iteration.
The arrays x[], y[] and z[] are defined in thorn CartGrid3d. As you will see you can use the variables of other thorns if you specify it in the interface.ccl (see Interface.ccl below).
CCTK_REAL dx, dt, dx2, dt2, factor; CCTK_REAL dx2i;
CCTK_REAL is one of the Cactus variables types.
dx = CCTK_DELTA_SPACE( 0 ); dt = CCTK_DELTA_TIME;
We extract the information about the space and time intervals. For simplicity we have done that the spacing is equal in all directions, i.e. dx = dy = dz. That's the reason why we use only CCTK_DELTA_SPACE( 0 ).
factor = 2*( 1 - 3*dt2*dx2i );
for(k=kstart; k< kend; k++)
{
for(j=jstart; j< jend; j++)
{
for(i=istart; i< iend; i++)
{
index = CCTK_GFINDEX3D( cctkGH, i, j, k );
phi[ index ] = factor*phi_p[ index ] - phi_p_p[ index ] +
dt2*dx2i * ( phi_p[ CCTK_GFINDEX3D( cctkGH, i+1, j , k ) ] +
phi_p[ CCTK_GFINDEX3D( cctkGH, i-1, j , k ) ] +
phi_p[ CCTK_GFINDEX3D( cctkGH, i , j+1, k ) ] +
phi_p[ CCTK_GFINDEX3D( cctkGH, i , j-1, k ) ] +
phi_p[ CCTK_GFINDEX3D( cctkGH, i , j , k+1) ] +
phi_p[ CCTK_GFINDEX3D( cctkGH, i , j , k-1) ] );
}
}
}
This part make the evolution of the wave. As you can see phi depends on the values of phi_p and phi_p_p. After each iteration these timelevels are rotated by the driver (PUGH), i.e. the values of phi are assigned to phi_p, those of phi_p to phi_p_p and so on, depending on the number of timelevels you need.
void wave3D_Boundary( CCTK_ARGUMENTS )
{
DECLARE_CCTK_ARGUMENTS
/* planes normal to x axis */
i = 0;
for(j=jstart-1; j< jend+1; j++)
{
for(k=kstart-1; k< kend+1; k++)
{
index = CCTK_GFINDEX3D( cctkGH, i, j, k );
phi[ index ] = phi_p[ index ];
}
}
.
.
.
Finally this functions sets the proper boundary conditions. In this case we have chosen a fixed one. Here we just apply the previous concepts.
The last step is to make that your source file be compiled. To do this just edit the file make.code.defn that is in the src directory of your thorn. Complete the line that starts with 'SRC =' with the name of your source file. This file must look as follows
# Main make.code.defn file for thorn wave3D # Source files in this directory SRCS = wave3D.c # Subdirectories containing source files SUBDIRS =
To make a new thorn go to the Cactus directory and type
gmake newthorn
The program will ask questions like the thorn name, the arrangement, author name, etc. and will create a group of new files that are present in any thorn. Supose that the name of your thorn is wave3D and you put it in an arrangement that you call CactusMythorns. Now go to arrangements/CactusMythorns there will be a directory called wave3D. Inside this directory there are a REAME file where you give a brief explanation of what you thorn is. There are also the interface.ccl, param.ccl and schedule.ccl. CCL stands for Cactus Configuration Language. CCL files are case independent, and may contain comments introduced by the '#' character, which indicates that the rest of the line is a comment. If the last non-blank character of a line in a CCL file is a backslash '\' the following line is treated as a continuation of the current line.
The code for interface.ccl is as follows
implements: wave3D
inherits: Grid
private:
CCTK_REAL scalarevolve type=GF dim=3 timelevels=3
{
phi
} "The wave function"
The first statement specifies the name of the thorn you are implementing, i.e. wave3D. The second line means that your thorn inherits all the public variables of thorn grid. You do this in order to use the public variables of any other thorn.
Now we declare our only 'Cactus variable' in our thorn. It is a private variable, which means that it can not be inherited by another thorn. You specify that this variable to be CCTK_REAL, that belongs to a group called 'scalarevolve', that is a grid function (type=GF), it has 3 dimensions (dim=3) and 3 timelevels (timelevels=3). You give the name of your variable inside curly brackets, the name is 'phi' in this case. Finally is strongly recommended that you put a small comment of what this variable is. It is done after you close the curly brackets.
There are many more options you can use to declare you variables.
STORAGE: scalarevolve[3]
schedule wave3D_generate_ID at INITIAL
{
LANG: C
} "Boundary and initial conditions"
schedule wave3D_Evolution at EVOL
{
LANG: C
} "Evolution of 3D wave equation"
schedule wave3D_Boundary at EVOL after wave3D_Evolution
{
LANG: C
SYNC: scalarevolve
} "Set up fixed boundary conditions"
The first line is needed to allocate memory to the group 'scalarevolve'. The number in brackets is the number of timelevels your variable has. (Remember that 'phi' belongs to 'scalarevolve' group).
What follows is very interesting. Here you tell when Cactus will execute the funtions that are contained in your source file. Cactus has several time bins. You must specify what function must be in what time bin. In this case we are scheduling 'wave3D_generate_ID' at time bin 'INITIAL'. This take care of the initial and boundary conditions. The main loop is performed by 'EVOL' time bin, here we schedule 'wave3D_Evolution'; which is the funtion that integrates the wave equation. After doing this we schedule 'wave3D_Boundary' again in 'EVOL' but now we add 'after wave3D_Evolution'. In this way we make that these two functions be executed one after each other.
You must specify the language (LANG: C) you used to write your code. In this case everything is C code. The SYNC keyword is used to specify the groups of variables that need to be synchonized, i.e. their ghostzones need to be exchanged. This is very important because if you forget to do it, you'll notice that your results depend on the number of processors you used to run the program!
Last but not least, it's a good idea to write some comment about what each funtion does. Comments are at the end of each 'schedule' statement. They are printed at runtime along with all other Cactus information.
Now you are almost ready to compile you code. What you need is a thornlist. It is a file in which you specify the thorns to be compiled. For your thorn to work you will need a thorn list like this
# Thorn list for wave3D CactusMythorns/wave3D CactusPUGH/PUGH CactusPUGH/PUGHReduce CactusPUGH/PUGHSlab CactusBase/Time CactusBase/IOASCII CactusBase/IOBasic CactusBase/IOUtil CactusBase/CartGrid3D
The format is < arrangement > / < thorn > . Copy and paste this to a file with name 'wave3D.th'. Now go to the Cactus directory and proceed in the same way you did with Cactus WaveDemo.
gmake wave3D-config THORNLIST=wave3D.th options=MyConf.conf
Now compile
gmake wave3D
If everything went all right by now your program is ready to run. At this time you'll need a parameter file. A parameter file (or parfile) contains information needed by thorn in order to operate. A possible parfile could be
ActiveThorns = "wave3D time pugh pughreduce pughslab ioutil ioascii iobasic cartgrid3d" # grid cactus::cctk_itlast = 50 #number of iteration to run driver::global_nx = 10 #number of points in the x direction driver::global_ny = 10 #number of points in the y direction driver::global_nz = 10 #number of points in the z direction grid::type = "byspacing" grid::domain = "full" grid::dxyz = 0.1 time::dtfac = 0.1 # output IO::out_dir = "wave3D" IOBasic::outInfo_every = 1 IOBasic::outInfo_vars = "wave3D::phi" #maximum and minimum of phi IOASCII::out_format = ".13e" IOASCII::out1D_every = 1 IOASCII::out1D_vars = "wave3D::phi" cactus::cctk_timer_output="full" pugh::timer_output="yes"
The first line specifies the thorns that need to be activated. Note that they are the same that those of the thornlist. A program may contain much more thorns compiled but if you don't use them you don't have to activate them. If you try to activate a thorn that is not in your thornlist then you'll get an error message. If this happens add the thorn you want to activate to your thorn list. Make the configuration and compile again.
In the parameter file you specify things like the size of the grid, spacing on each direction, how many iteration you want to do, how often to output the results, etc. For a complete explanation of these parameters see the documentation of each particular thorn, you'll find it in the 'doc' directory of each thorn.
To run your program in a single processor type
mpirun -np 1 ./exe/cactus_wave3D wave3D.par
have fun!