GGCM library
From OpenGGCM
Contents |
Intro
The GGCM library is a collection of C routines for reading and working with the 'MHD' file type, as output by Professor Joachim Raeder's OpenGGCM model.
This page used to be on Tom's artemis site, but hopefully everything that is there is now moved to here.
Downloads
Media:ggcm-2.2.0.tar.gz -- consts some new things, code cleanups and maintainability fixes. Final release by Tom Fogal.
Media:ggcm-2.1.1.tar.gz -- fixes a problem with the VTK reader which may affect some users.
Media:ggcm-2.1.0.tar.gz -- can now read timecode information from 3df files. Bug Fix related to GSE conversion. Improved documentation.
Media:ggcm-2.0.1.tar.gz -- cleanups, symbol renaming, some documentation.
Compilation
The build system is straightforward and noninteresting; simply run
./configure make sudo make install
on just about any UNIX clone. The library is developed and used on mostly Linux systems.
The library was the basis for developing the Visit Reader. Hopefully you'll be able to get by with just premade utilities which wrap the library, but if you cannot some documentation is included in the form of UNIX manual pages. Of note is the ggcm_read_mhd manual page, which contains an example on how to read in and manipulate GGCM data.
Also included with the library is a VTK reader for the data format. Most people will simply want to copy the vtkGGCMReader.[cxx|h] files and include them in their own project. However it is possible to build the reader as part of VTK, which has the advantage that it gets wrapped into Python / TCL / Java / whatever VTK is wrapped with. See the file 'INSTALL', included with the distribution.
Documentation
Online Doxygen
Importing doxygen into a wiki is probably not feasible, so this documentation is posted on Tom's artemis site.
There is documentation for the C library which includes documentation for the VTK reader.
C library Usage
The following example program reads in a 3df file and prints out the 'BY' value at an arbitrary location:
/** sample program to read in data from a `3df' file. It takes two arguments,
* the grid file followed by the data file, reads in all of the data, and
* prints out an element of it. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ggcm.h>
struct jdata {
char *data_fn; /** data filename */
char *grid_fn; /** grid filename */
float *grid[3];
float ***b[3];
float ***v[3];
float ***xj[3];
float ***rr;
float ***pp;
float ***resis;
unsigned int dim[3];
};
static void read_grid(struct jdata *, const char *fn);
static void read_data(struct jdata *, const char *fn);
static void do_your_algorithm(const struct jdata *);
static void usage(const char *s);
static void free_data(struct jdata *);
int main(int argc, char *argv[])
{
struct jdata data;
if(argc != 3) {
usage(argv[0]);
exit(1);
}
read_grid(&data, argv[1]);
read_data(&data, argv[2]);
do_your_algorithm(&data);
free_data(&data);
return 0;
}
static void usage(const char *s)
{
fprintf(stderr, "usage: %s grid_file data_file\n", s);
}
/** reads in grid data from the specified filename */
static void read_grid(struct jdata *jd, const char *fn)
{
MHDdata *grid;
const char *fields[] = {"gridx", "gridy", "gridz", NULL};
float *grd[3]; /* non-GSE data */
/* keep a record of the filename. */
jd->grid_fn = malloc(strlen(fn)+1);
strcpy(jd->grid_fn, fn);
grid = ggcm_read_mhd(fn, fields);
grd[0] = ggcm_grid(grid, &jd->dim[0], "gridx");
grd[1] = ggcm_grid(grid, &jd->dim[1], "gridy");
grd[2] = ggcm_grid(grid, &jd->dim[2], "gridz");
ggcm_free_data(grid);
/* at this point, `grd' is the grid data as written to the grid file -- aka
* in jimmy's coordinate system, NOT GSE. Also, jd->dim is initialized to
* be the dimensions of the grid. */
jd->grid[0] = malloc(sizeof(float) * jd->dim[0]);
jd->grid[1] = malloc(sizeof(float) * jd->dim[1]);
jd->grid[2] = malloc(sizeof(float) * jd->dim[2]);
ggcm_grid_gse(grd[0], grd[1], grd[2],
jd->grid[0], jd->grid[1], jd->grid[2],
jd->dim[0], jd->dim[1], jd->dim[2]);
/* now, jd->grid contains the GSE grid data. We won't need the old data
* anymore. */
free(grd[0]);
free(grd[1]);
free(grd[2]);
}
/** reads field data from the specified filename */
static void read_data(struct jdata *jd, const char *fn)
{
MHDdata *data;
const char *fields[] = {"bx","by","bz",
"vx","vy","vz",
"xjx","xjy","xjz",
"rr","pp","resis", NULL};
/* non-GSE temporary buffers: */
float ***b[3], ***v[3], ***xj[3];
float ***rr, ***pp, ***resis;
/* keep a copy of the filename */
jd->data_fn = malloc(strlen(fn) + 1);
strcpy(jd->data_fn, fn);
/* read in the data */
data = ggcm_read_mhd(fn, fields);
/* now pull that data into a temporary array; temporary because its in a
* Jimmy-coordinate system, not GSE. */
b[0] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "bx");
b[1] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "by");
b[2] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "bz");
v[0] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "vx");
v[1] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "vy");
v[2] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "vz");
xj[0] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "xjx");
xj[1] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "xjy");
xj[2] = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "xjz");
rr = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "rr");
pp = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "pp");
resis = ggcm_field_matrix(data, jd->dim[0], jd->dim[1], jd->dim[2], "resis");
/* allocate space for our GSE coordinate data */
jd->b[0] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->b[1] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->b[2] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->v[0] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->v[1] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->v[2] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->xj[0] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->xj[1] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->xj[2] = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->rr = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->pp = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
jd->resis = m_alloc(jd->dim[0], jd->dim[1], jd->dim[2]);
/* convert to GSE */
ggcm_mhd_gse(b[0], jd->b[0], jd->dim[0], jd->dim[1], jd->dim[2], "bx");
ggcm_mhd_gse(b[1], jd->b[1], jd->dim[0], jd->dim[1], jd->dim[2], "by");
ggcm_mhd_gse(b[2], jd->b[2], jd->dim[0], jd->dim[1], jd->dim[2], "bz");
ggcm_mhd_gse(v[0], jd->v[0], jd->dim[0], jd->dim[1], jd->dim[2], "vx");
ggcm_mhd_gse(v[1], jd->v[1], jd->dim[0], jd->dim[1], jd->dim[2], "vy");
ggcm_mhd_gse(v[2], jd->v[2], jd->dim[0], jd->dim[1], jd->dim[2], "vz");
ggcm_mhd_gse(xj[0], jd->xj[0], jd->dim[0], jd->dim[1], jd->dim[2], "xjx");
ggcm_mhd_gse(xj[1], jd->xj[1], jd->dim[0], jd->dim[1], jd->dim[2], "xjy");
ggcm_mhd_gse(xj[2], jd->xj[2], jd->dim[0], jd->dim[1], jd->dim[2], "xjz");
ggcm_mhd_gse(rr, jd->rr, jd->dim[0], jd->dim[1], jd->dim[2], "rr");
ggcm_mhd_gse(pp, jd->pp, jd->dim[0], jd->dim[1], jd->dim[2], "pp");
ggcm_mhd_gse(resis, jd->resis, jd->dim[0], jd->dim[1], jd->dim[2], "resis");
/* we won't need our old, non-GSE data anymore */
m_free(b[0]); m_free(b[1]); m_free(b[2]);
m_free(v[0]); m_free(v[1]); m_free(v[2]);
m_free(xj[0]); m_free(xj[1]); m_free(xj[2]);
m_free(rr); m_free(pp); m_free(resis);
}
/** replace this with your code. */
static void do_your_algorithm(const struct jdata *jd)
{
printf("The BY value at (%g,%g,%g) is [%g]\n",
jd->grid[0][12], jd->grid[1][6], jd->grid[2][4], jd->b[1][12][6][4]);
}
/** frees all the memory allocated previously. */
static void free_data(struct jdata *jd)
{
free(jd->data_fn);
free(jd->grid_fn);
/* grid data was allocated with system's malloc */
free(jd->grid[0]); free(jd->grid[1]); free(jd->grid[2]);
/* data was allocated via m_alloc; free with m_free */
m_free(jd->b[0]); m_free(jd->b[1]); m_free(jd->b[2]);
m_free(jd->v[0]); m_free(jd->v[1]); m_free(jd->v[2]);
m_free(jd->xj[0]); m_free(jd->xj[1]); m_free(jd->xj[2]);
/* scalar data: */
m_free(jd->rr); m_free(jd->pp); m_free(jd->resis);
}
The above code reads all the available data in a file, which may or may not be required in your application. Also note in read_data that we call a block of ggcm_field_matrix (which allocates memory), followed by a block m_allocs, followed by a block of ggcm_mhd_gse. In a real application, you should be concerned with 1) that any of these routines can fail due to lack of memory, returning NULL, and 2) the memory usage of the application. It would make more sense to do:
b[0] = ggcm_field_matrix(...); jd->b[0] = m_alloc(...); ggcm_mhd_gse(b[0], jd->b[0], ..., "bx"); m_free(b[0]) b[1] = ggcm_field_matrix(...); jd->b[1] = m_alloc(...); ggcm_mhd_gse(b[1], jd->b[1], ..., "by"); m_free(b[1]);
Note how this fragment reorders the operations so that we can m_free memory as soon as possible -- immediately after converting to GSE coordinates -- and before other allocations occur. Our non-GSE b[0] (BX) data is thus deallocated before we allocate any memory for the BY field. With large data files, this can be very important -- the system might have enough RAM to fit a whole data file in memory, but if you wait to read and GSE-ify the whole data file before freeing the non-GSE data arrays, then the system will need enough RAM to fit two of the data files in memory.
VTK Reader
The following is the example program included with the GGCM library distribution. You can find it in /vtk_mhd/objtest.cc.
#include "vtkGGCMReader.h"
#include "vtkContourFilter.h"
#include "vtkActor.h"
#include "vtkRenderWindow.h"
#include "vtkRenderer.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkPolyDataMapper.h"
#include "vtkRectilinearGridGeometryFilter.h"
#include "vtkProperty.h"
#include "vtkCamera.h"
int main()
{
char *grid_list[] = {"gridx", "gridy", "gridz", NULL};
vtkGGCMReader *mhd = vtkGGCMReader::New();
vtkPolyDataMapper *map = vtkPolyDataMapper::New();
vtkActor *actor = vtkActor::New();
vtkRenderer *ren = vtkRenderer::New();
vtkRenderWindow *ren_wind = vtkRenderWindow::New();
vtkRenderWindowInteractor *ir = vtkRenderWindowInteractor::New();
vtkRectilinearGridGeometryFilter *plane = vtkRectilinearGridGeometryFilter::New();
mhd->SetGridFile("/home/tfogal/data/jcd0006.grid");
mhd->SetDataFile("/home/tfogal/data/jcd0006.3df.004800");
/* Don't do this type of thing in your programs. I just do it this way
* so I can test both methods of setting fields */
mhd->SetGridList(grid_list);
mhd->AddField("bx");
mhd->AddField("by");
mhd->AddField("bz");
mhd->Read();
plane->SetInput(dynamic_cast<vtkRectilinearGrid*>(mhd->GetOutput()));
plane->SetExtent(0,391, 0,111, 0,111);
map->SetInput(plane->GetOutput());
actor->SetMapper(map);
actor->GetProperty()->SetRepresentationToWireframe();
actor->GetProperty()->SetColor(0,0,0);
ren_wind->AddRenderer(ren);
ir->SetRenderWindow(ren_wind);
ren->AddActor(actor);
ren->SetBackground(1,1,1);
ren_wind->Render();
ir->Start();
return 0;
}
I will not explain VTK here, but this describes the basic usage of our VTK reader. Simply allocate a New() reader, set the grid and data files, set which fields to read, and call Read(). The reader will then behave like a normal VTK reader and Update() appropriately.
