Cactus on Unstructured Meshes with CaPHG

Jian Tao (
2039 Digital Media Center
Center for Computation and Technology
Louisiana State University,
Baton Rouge, LA 70803

(last updated: 07/04/2014)

                 o-o     o--o  o  o  o-o 
                /        |   | |  | o 
               O      oo O--o  O--O |  -o 
                \    | | |     |  | o   |
                 o-o o-o-o     o  o  o-o

             Cactus Parallel Hierarchical Grid

                 (c) Copyright The Authors 
                 GNU Licensed. No Warranty

Background of CaPHG

The goal of CaPHG is to integrate Cactus with PHG, an FEM library to help the application development with both FEM and DG methods on unstructure meshes. The motivation is to take advantage of the programmability of the Cactus computational framework and the performance and scalability of the PHG library to provide an integrated problem solving environment to solve large scale scientific problems on high-end computing facilities.

The primary task that I have been working on is to properly map Cactus GF to PHG DOF (degree of freedom) while extending the structured cGH (Cactus Grid Hierarchy) to unstructured PHG (Parallel Hierarchical Grid). The trick to minimize the changes in Cactus flesh is that I limited the variable type of GF in Cactus to DOF in PHG. Namely, I kept the name GF as the group type, but the variable type now can only be DOF instead of previously used CCTK_REAL, CCTK_INT, etc. By default, we will use CCTK_REAL for all numerical values of the unknowns.

DOF is defined in PHG to handle variables and relevant operations on a triangular grid. The DOF_TYPE in the tag refers to different types of DOF used. Different basis functions will be automatically associated with the type of DOF defined. The built-in integrator, differential operators will take that into consideration automatically. Unlike structured codes where a user specifies the number of grid points to decide the problem size, CaPHG gets that information from a mesh file.

Note: Adding a namespace, such as **DOMAIN** in interface.ccl might be a better solution, which can potentially support the co-existence of both unstructured and structured meshes and codes.

Grid Abstractions in Cactus

The Cactus flesh together with the Cactus computational toolkit contain a collection of data structures and functions that to support operations on a dynamically distributed grid hierarchy and its grid functions (fields living on the grid hierarchy). These structures and functions can be categorized into the following three grid abstractions, which usually appear in high level programming frameworks for parallel block-structured

A simple maping between Cactus and PHG would be

PHG Data Structures

PHG DOF Structure

The DOF structure holds the information about the DOFs on a given grid. It has a built-in caching block to store some relavant data, such as the gradient of the DOF, etc. Functions defined in PHG could directly carry out algebraic, differential, and integral operations on DOFs.

typedef struct DOF_ {
    int         magic;          /* magic number */
    void        *cache_func;    /* cached function values */
    char        *name;          /* name of the dof */
    const char  *srcfile;       /* src file name in which DOF is created */
    GRID        *g;             /* the mesh */

    /* Note: exactly one of 'type' and 'hp' members is not NULL,
     * and 'hp' != NULL means a variable order DOF */
    DOF_TYPE    *type;          /* type of the DOF */
    HP_TYPE     *hp;            /* variable order DOF (h-p adaptivity) */

    /* data pointers */
    FLOAT       *data;          /* length = dim * (
                                        nvert * np_vert + nedge * np_edge +
                                        nface * np_face + nelem * np_elem ) */
    FLOAT       *data_vert;     /* pointer to vertex data */
    FLOAT       *data_edge;     /* pointer to edge data */
    FLOAT       *data_face;     /* pointer to face data */
    FLOAT       *data_elem;     /* pointer to volume data */

    /* (x,y,z) function attached to DOF */
    DOF_USER_FUNC userfunc;
    DOF_USER_FUNC_LAMBDA userfunc_lambda;

    /* some constants */
    size_t      count_vert;     /* data count per vertex */
    size_t      count_edge;     /* data count per edge */
    size_t      count_face;     /* data count per face */
    size_t      count_elem;     /* data count per volume */
    int         srcline;        /* line no in src file */

    BOOLEAN     invariant;      /* TRUE ==> same values on all elems,
                                 * default FALSE */
    SHORT       dim;            /* # of variables */
    BTYPE       DB_mask;        /* Dirichlet boundary mask */
    BTYPE       *DB_masks;      /* Dirichlet boundary masks for each dof dim */
} DOF;

PHG DOF Type Structure

DOF_TYPE defines the type of a DOF. Each DOF_TYPE object describes the properties that certain type of DOF shall have. For instance, it defines for a DOF what type of basis functions and their corresponding order, what interpolator to use to switch between the coarse and fine grids, etc.

struct DOF_TYPE_ {
    /*----------------------- pointers in DofReserved -----------------------*/
    /* caches for 3D quadrature */
    void    *cache_basfunc, *cache_basgrad; /* cached bases, Dbas/Dlambda */
    void    *cache_gradient, *cache_curl;   /* cached Dbas/Dxyz */

    struct MASS_LUMPING_ *mass_lumping;     /* mass lumping struct */
    void (*Initialize)(struct DOF_TYPE_ *); /* Initialization function */

    const char  *name;      /* name (or description) */
    FLOAT   *points;    /* barycentric coordinates of the points,
                   (np_vert + np_edge*2 + np_face*3 + np_elem*4)
                   entries */
    BYTE    *orders;    /* polynomial order of individual basis
                   functions (may be NULL), should contain
                   np_vert + np_edge + np_face + np_elem
                   entries */

    DOF_TYPE    *grad_type; /* DOF type of the gradients */
    DOF_TYPE    *base_type; /* for DGn: the corresponding Pn/HBn type */
    HP_INFO *hp_info;   /* NULL for non-hierarchic bases */

    /* pointers to DOF functions */
    DOF_INTERP_FUNC InterpC2F; /* coarse-fine interpolation function */
    DOF_INTERP_FUNC InterpF2C; /* fine-coarse interpolation function */
    DOF_INIT_FUNC   InitFunc;  /* function for assigning DOF values */
    DOF_BASIS_FUNC  BasFuncs;  /* function for evaluating basis funcs */
    DOF_BASIS_GRAD  BasGrads;  /* function for evaluating gradients */
    DOF_MAP_FUNC    DofMap;   /* function for ordering face data */

    FE_SPACE    fe_space;   /* the FE space the functions belong to */

    BOOLEAN invariant;  /* TRUE ==> same bas functions for all elems */
    BOOLEAN free_after_use; /* TRUE ==> DOF_TYPE is dynamically allocated
                        and should be freed after last
                        DOF using it is freed */

    SHORT   id;     /* id (assigned when referenced) */

    SHORT   nbas;       /* number of basis functions per element:
                    NVert * np_vert + NEdge * np_edge +
                    NFace * np_face + np_elem */
    BYTE    order;      /* polynomials order */
    CHAR    continuity; /* continuity across element boundaries:
                   < 0 means discontinuous, 0 means C^0,
                   1 means C^1, etc. */
    SHORT   dim;        /* space dimension (# of values) */

    SHORT   np_vert;    /* # of DOFs per vertex (must be 0 or 1) */
    SHORT   np_edge;    /* # of DOFs per edge */
    SHORT   np_face;    /* # of DOFs per face */
    SHORT   np_elem;    /* # of DOFs per element */

PHG GRID Structure

This structure holds all the information PHG needs to partition and distribute unstructured meshes. It also helps to balance the load and handle the coarsening and refinement of the grids adaptively via an error estimator from the user.

typedef struct GRID_ {
    FLOAT   lif;        /**< Load imbalance factor */
    FLOAT   bbox[2][Dim];   /**< BoundingBox of the mesh */
    FLOAT   volume;     /**< Volume of the mesh */
    char    *filename;  /**< name of the mesh file */
    COORD   *verts;     /**< List of vertices */
    ELEMENT *roots;     /**< Root elements */
    ELEMENT **elems;    /**< Array of size g->nlem, pointers to leaf
                     elements (for mesh traversal), NULL if
                     an element index is unreferenced.
                     Note for any i, i>=0 && i<nelems,
                     the following assertion holds:
                     assert(elems[i] == NULL ||
                        elems[i]->index == i) */
    GRID_ELEM_TYPE  elem_type;  /**< Grid element types */  
    struct MAPPING_ *mapping;/**< Mapping from reference coord to
                      real world coord */
    struct DOF_ *geom;      /**< DOF containing the geometric data,
                     maintained by geom.c */
    struct DOF_ **dof;      /**< NULL terminated list of DOFs
                     defined on the mesh */
    struct HP_TYPE_ **hp;   /**< NULL terminated list of HP_TYPEs
                     defined on the mesh */

    EXPR    **bdry_funcs;   /**< NULL terminated list of projection
                   functions for curved boundaries */
    BTYPE   *types_vert;    /**< Types of vertices (bit flags) */
    BTYPE   *types_edge;    /**< Types of edges (bit flags) */
    BTYPE   *types_face;    /**< Types of faces (bit flags) */
    BTYPE   *types_elem;    /**< Types of elements (bit flags) */

    /* inter-grid links */
    struct GRID_ *alien;    /**< Pointer to another grid */
    INT     *alien_map; /**< Map of elements between the two grids */

    /* struct for handling periodic BC */
    struct PERIOD_ *period;

    /* Arrays containing ranks of, and local indices in the owner process
     * for vertices. They are non nil only when needed and are updated by
     * phgUpdateBoundaryTypes in utils.c and used by build_L2Vmap() in map.c.
     * Note: an entry has value -1 iff it's uniquely owned. */
    INT     *owner_index_vert;
    SHORT   *owner_rank_vert;

    struct {
    /* elements with neighbours on other processes are stored here */
    int     *counts;    /* Counts of remote faces */
    int     *displs;    /* Start positions in list */
    RNEIGHBOUR  *list;      /* List of remote faces */
    INT     count;      /* Total number of remote faces */
    INT     allocated;  /* Allocated size of list */
    }       neighbours; /**< List of remote faces */
    /* Note: don't access L2Gmaps directly, use the GlobalXxxx macros instead */
    INT     *L2Gmap_vert;   /**< Local to global map of vertices */
    INT     *L2Gmap_edge;   /**< Local to global map of edges */
    INT     *L2Gmap_face;   /**< Local to global map of faces */
    INT     *L2Gmap_elem;   /**< Local to global map of element indices */

    /* Arrays containing ranks of, and local indices in the owner process
     * for edges, faces and elements. They are non nil only when
     * needed and are updated by phgUpdateBoundaryTypes in utils.c
     * and used by build_L2Vmap() in map.c.
     * Note: an entry has value -1 iff it's uniquely owned. */
    INT     *owner_index_edge;
    SHORT   *owner_rank_edge;

    /* Note: owner_xxxx_face and owner_xxxx_elem are currently not used */
    INT     *owner_index_face;
    SHORT   *owner_rank_face;

    INT     *owner_index_elem;
    SHORT   *owner_rank_elem;

    MPI_Comm    comm;       /**< The mesh's communicator */
#else   /* USE_MPI */
    int     comm;
#endif  /* USE_MPI */
    INT     nleaf;      /**< Number of leaf elements in the submesh */
    INT     nvert;      /**< Number of vertices in the submesh */
    INT     nedge;      /**< Number of edges in the submesh */
#if Dim == 3
    INT     nface;      /**< Number of faces in the submesh */
    INT     nelem;      /**< Number of element indices in the submesh */

    /* counts of owned vertices/edges/faces/elements, they are updated by
     * the function phgUpdateBoundaryTypes */
    INT     nvert_owned;
    INT     nedge_owned;
    INT     nface_owned;
    INT     nelem_owned;

    /* global counters. 2D, V:E:T \approx 1:3:2. 3D, V:E:F:T \approx 1:6:10:5 */
    INT     nvert_global;   /**< Number of vertices in the global mesh */
    INT     nedge_global;   /**< Number of edges in the global mesh */
#if Dim == 3
    INT     nface_global;   /**< Number of faces in the global mesh */
    INT     nelem_global;   /**< Number of tetrahedra in the global mesh */

    INT     nroot;      /**< Number of root elements */
    INT     ntree;      /**< Total number of elements in the tree */

    int     serial_no;  /**< Number of grid changes */
    int     last_partitioner; /**< Last partitioner used */

    int     rank;       /**< Process rank (submesh no) */
    int     nprocs;     /**< Number of processes for the mesh */

    int     bc_alloc;   /**< Allocated size of bc_index and bc_rmap */
    int     bc_n;       /**< Actual size of bc_index and bc_rmap */
    int     *bc_list;   /**< sorted list of BDRY types */
    int     *bc_rmap;   /**< PHG:bc_list[i] <==> Medit:bc_rmap[i] */

    BYTE    flags;      /**< {VERT,EDGE,FACE,ELEM,GEOM}_FLAG bits */

As one may guess, some changes in Cactus flesh are necessary, but the changes have been kept to minimal and most changes are constrained in the unstructured mesh driver CaPHG, which handles the memory allocation and creates a new set of time bins that could be better mapped to solvers built on top of unstructured meshes. CaPHG registers an extension and overloads a couple of functions to mimic what a structured mesh driver will do in Cactus. Following the basic design ideas of Cactus 4.0, CaPHG makes it relatively easy for Cactus users to switch to unstructured meshes without worrying about the handling of unstructured meshes and refinement in a distributed computing environment at all.

Implementation Details

In Chapter C2 (Infrastructure Thorns) of the Cactus Users Guide (, you may find instructions on how to write a Cactus driver. It is actually surprisingly easy to get started. However, when you begin to think about supporting real world applications, as well as issues on performance and scalability, writing a driver will quickly become complicated. Fortunately, one of Cactus design considerations to detach Grid Function from Grid Hiearchy make it possible to switch to an external library, PHG in this case, to handle a completely different type of Grid Hiearchy without dramatically modifying the flesh.

CaPHG depends on PHG for parallelization, load balancing, mesh partitioning, mesh refinement, internal error checking, and performance measurement. To facilitate the error estimation, CaPHG defines a DOF to hold the error that an application developer needs to calculate during the adaptive iteration. The L1 norm of this DOF will be calculated in CaPHG as an indicator for adaptive mesh refinement.

Here I list some of the key functions implemented in CaPHG. Functions like CCTK_MyProc and CCTK_nProcs could be easily implemented and I will skip those functions here.

There are two main structures defined in CaPHG: CaPHGGroup and CaPHGExt.


typedef struct CaPHGGroup_
    int vartypesize; /* number of bytes of variable type */
    int numtimelevels; /* number of timelevels in group */
    int numvars; /* number of variables in group */
    int first_var; /* index of first variable in group */
    int storage; /* must be an int, because returned by EnableGroupStorage() */
    int tagstable; /* tag table for storing extra information of the DOF groups*/
    CCTK_DOF ***data; /* *data[var][tl] */

/* the following information shall be passed in via code generation at the CST stage */
    DOF_TYPE *dof_type; /* DOF_TYPE of the DOF */
    SOLVER *solver; /* solver for the group */
    DOF_USER_FUNC user_func; /* user defined function to set the DOF */
  } CaPHGGroup;


  typedef struct CaPHGExt_
    GRID *grid;
    int numgroups;
    CaPHGGroup *groups;
  } CaPHGExt;

Here are some of the main functions implemented in CaPHG to fulfill its role as a driver for Cactus.

TODO: The default Evolve function is currently used. I am working on a customized version of it to take advantage of the extra information I get from the tag table to integrate an explicit integrator similar to MoL in the driver CaPHG.

startup.c : StartUp

int CaPHG_StartUp(void)
  /* Register the GH extension */
  GHEx_Handle = CCTK_RegisterGHExtension(DRIVERNAME);
  CCTK_RegisterGHExtensionSetupGH(GHEx_Handle, CaPHG_SetupGH);

  /* Overload necessary functions */
  CCTK_OverloadExit((int (*)(cGH*, int)) CaPHG_Exit);
  CCTK_OverloadAbort((int (*)(cGH*, int)) CaPHG_Abort);

caphg.c : SetupGH

void *CaPHG_SetupGH(tFleshConfig *config, int conv_level, cGH *GH)
  g = phgNewGrid(-1);

  phgSetPeriodicity(g, periodicity);

  if (!phgImport(g, input_mesh_file, TRUE))
    /* PHG shall abort the program if phgImport fails, the error message below shall never be called.*/
    CCTK_ERROR("Failed to import the mesh file!");
  GHExt->grid = g;
  GHExt->numgroups = CCTK_NumGroups();

      GHExt->groups = malloc(GHExt->numgroups * sizeof(CaPHGGroup)));

  /* iterate through all the groups */
  for (gi = 0; gi < GHExt->numgroups; gi++)
    CCTK_GroupData(gi, &groupdata);
    GHExt->groups[gi].vartypesize = CCTK_VarTypeSize(groupdata.vartype);
    GHExt->groups[gi].numtimelevels = groupdata.numtimelevels;
    GHExt->groups[gi].numvars = groupdata.numvars;
    GHExt->groups[gi].first_var = first_var;
    GHExt->groups[gi].storage = 0; /* off by default */
    GHExt->groups[gi].tagstable = CCTK_GroupTagsTableI(gi);
    GHExt->groups[gi].dof_type = groupdata.dof_type;
    GHExt->groups[gi].solver = NULL;
    GHExt->groups[gi].user_func = groupdata.user_func;

    first_var += groupdata.numvars;
        GHExt->groups[gi].data = malloc(groupdata.numvars * sizeof(void **)));
    for (vi = 0; vi < groupdata.numvars; vi++)
          GHExt->groups[gi].data[vi] = malloc(
              groupdata.numtimelevels * sizeof(void *)));
      for (tl = 0; tl < groupdata.numtimelevels; tl++)
        GHExt->groups[gi].data[vi][tl] = NULL;
  return GHExt;

caphg.c : TraverseGH

int CaPHG_ScheduleTraverseGH(cGH *GH, const char *where)
  /* these two numbers won't affect CaPHG. */
  GH->cctk_convlevel = 1;
  GH->cctk_convfac = 2;

  /* timefac wil be used in CaPHG */
  GH->cctk_timefac = 1;
  GH->grid = GHExt->grid;
  /*TODO: the default shutdown takes care of convergence levels which
   * is not necessary for adaptive unstructured grids. Here we skip
   * the CCTK_SHUTDOWN time bin. Need to get back to this later on if
   * necessary. */

  if (strcmp("CCTK_SHUTDOWN", where))
    /* variable groups and relations */
    for (gi = 0; gi < GHExt->numgroups; gi++)
      for (vi = 0; vi < GHExt->groups[gi].numvars; vi++)
        for (tl = 0; tl < GHExt->groups[gi].numtimelevels; tl++)
          GH->data[GHExt->groups[gi].first_var + vi][tl] =

caphg.c : EnableStorage

int CaPHG_EnableGroupStorageI(const cGH *GH, int gi)
  /* enable the storage only when it is not enabled yet. */
  if (retval != 1)
    /* in CaPHG, Grid Functions will be PHG DOFs. */
    for (vi = 0; vi < g.numvars; vi++)
      gvi = g.first_var + vi;
      vtype = CCTK_VarTypeI(gvi);
      if (gtype == CCTK_GF)
        if (vtype == CCTK_VARIABLE_DOF)
          for (tl = 0; tl < g.numtimelevels; tl++)
            gf = phgDofNew(GHExt->grid, g.dof_type, 1, CCTK_VarName(gvi),

  [vi][tl] = gf;
          } // for (tl = 0; tl < g.numtimelevels; tl++)
        /* TODO: this shall be checked at the parsing stage. */
              "Failed to create a DOF! Please check your tag table in interface.ccl.");
      else if (gtype == CCTK_SCALAR)
        if (vtype == CCTK_VARIABLE_DOF)
          /* TODO: this shall be checked at the parsing stage. */
          CCTK_ERROR("CCTK_DOF only works with grid functions in CaPHG!");
          for (vi = 0; vi < g.numvars; vi++)
            for (tl = 0; tl < g.numtimelevels; tl++)
              gvi = g.first_var + vi;
        [vi][tl] = malloc(
                      g.numtimelevels * CCTK_VarTypeSize(CCTK_VarTypeI(gvi))));
        CCTK_ERROR("CaPHG now only support CCTK_GF and CCTK_SCALAR");
    } = 1;
  } //if (retval != 1)
  return retval;

caphg.c : QueryGroupStorage

int CaPHG_QueryGroupStorage(const cGH *GH, int group, const char *groupname)
  CaPHGExt *GHExt;
  int gi;

  if (groupname)
    gi = CCTK_GroupIndex(groupname);
    gi = -1;

  if (gi < 0)
    /* Invalid groupname, try group index */
    if ((gi = group) < 0 || group >= CCTK_NumGroups())
      gi = -10;

  if (gi >= 0)
    /* gi is valid, grab the GHExt */
    GHExt = CCTK_GHExtension(GH, DRIVERNAME);

  return gi < 0 ? gi : GHExt->groups->storage;

caphg.c : GroupDynamicData

int CaPHG_GroupDynamicData(const cGH *GH, int group, cGroupDynamicData *data)
  CaPHGExt *GHExt;

  /* Grab GH extension */

  /* Dimension has no meaning for unstructured meshes. */
  data->dim = 1;

  /* not relevant to CaPHG */
  data->gsh = NULL;
  data->lsh = NULL;
  data->nghostzones = NULL;
  data->lbnd = NULL;
  data->ubnd = NULL;
  data->bbox = NULL;
  data->activetimelevels = GHExt->groups[group].storage
      * GHExt->groups[group].numtimelevels;
  return EXIT_SUCCESS;

A Simple Example with CaPHG

The standalone versoin of this example using PHG directly can be found in the appendix at the end of this article.

Based on a standalone PHG code, one example using CaPHG was implemented. The goal is to reduce the size of the standalone code as much as possible by moving the common components into the driver.

Here comes the interface.ccl and schedule.ccl file of the testing thorn that solves the Helmholtz equation: Δu+au=f.


implements: phg_examples
inherits: caphg

CCTK_DOF main type=GF timelevels=1 tags='dof_type="DOF_DEFAULT" solver="SOLVER_DEFAULT" dof_user_func="DofInterpolation"'
} "main unknowns of a PDE"

CCTK_DOF analytic_solution type=GF timelevels=1 tags='dof_type="DOF_ANALYTIC" dof_user_func="myu"'
} "analytical solution"

CCTK_DOF rhs_function type=GF timelevels=1 tags='dof_type="DOF_DEFAULT" dof_user_func="myf"'
} "right hand side function"


STORAGE: main analytic_solution rhs_function
schedule simplest in CUMT_ITERATION
} "Simplest testing case for CaPHG"

schedule estimate_error IN CUMT_ESTIMATE_ERROR
} "estimate error"

Now, only three functions (will be two after CaPHG_Evolve is ready) need the input from an application developer to build a linear system, estimate error, and schedule a controlling function (to be removed) to call the linear system constructor and solver.

linear system builder

void build_linear_system(SOLVER *solver, DOF *u_h, DOF *f_h)
/* application specific part */

/* this routine just needs to calculate the error, a DOF whose L1 Norm
 * will be used to compare against the iteration_tolerance set in the
 * parameter file */

error estimator

void estimate_error( CCTK_ARGUMENTS)
/* application specific part */

controlling function

/* Main routine: construct and call a linear solver for each adaptive iteration */

void simplest( CCTK_ARGUMENTS)
  SOLVER *solver;
  solver = phgSolverCreate(SOLVER_DEFAULT, u_h, NULL);
  build_linear_system(solver, u_h, f_h);
  phgSolverSolve(solver, TRUE, u_h, NULL);
  phgPrintf("Final mesh written to \"%s\".\n",
  phgExportVTK(cctkGH->grid, "simplest.vtk", u_h, error, NULL));

Some special (x,y,z,value)-type numerical functions shall be implemented by the user. The function names could be directly used in the interface.ccl file to initialize the DOF. For instance, myu and myf are user-defined functions. Some changes in one of the CST scripts, is necessary to make this magic happen.

void myu(FLOAT x, FLOAT y, FLOAT z, FLOAT *value)
    *value = Cos(2. * M_PI * x) * Cos(2. * M_PI * y) * Cos(2. * M_PI * z);

void myf(FLOAT x, FLOAT y, FLOAT z, FLOAT *value)
  myu(x, y, z, value);
    *value = 12. * M_PI * M_PI * *value + a * *value;

Parameter file

ActiveThorns = "phg caphg phgexamples"

# flesh parameters
Cactus::cctk_run_title = "CaPHG Example (Solving the Helmholtz equation -\Delta u + a u = f)"
Cactus::cctk_show_schedule = "yes"
Cactus::cctk_show_banners = "yes"

# time_dependent is default to no. so time dependent problems need to set it to yes. 
# cctk_itlast shall be set accordingly. CaPHG resets cctk_itlast to 1 if time_dependent = "no"

CaPHG::time_dependent = "no"
CaPHG::input_mesh_file = "cube4.dat"
CaPHG::periodicity = 0
CaPHG::verbose= "no"
CaPHG::mem_max = 2000
CaPHG::iteration_tolerance = 0.4

Mesh file

DIM: 3
number of vertices: 125
number of elements: 384

vertex coordinates:
 0 0 0
 1 0 0
 0.5 0.75 1
 0.25 0.5 1

element vertices:
 0 27 35 71
 15 27 35 71
 26 34 69 124
 12 34 69 124

element boundaries:
 0 1 0 0
 0 1 0 0
 0 0 0 0
 0 0 0 0

element type:

element neighbours:
 1 -1 64 192
 0 -1 65 3
 383 374 379 342
 382 375 378 381

Standard Output

##### MPI_L2setup: setting up two-level topology.
##### MPI_L2setup: 1 node or 1 proc/node, do nothing. #####
Parallel Hierarchical Grid (version 0.9.1).
number of process(es): 1
SPMD mode.
phgMaxThreads = 1

  1   0101       ************************  
  01  1010 10      The Cactus Code V4.2.3    
 1010 1101 011     
  1001 100101    ************************  
     100011     (c) Copyright The Authors  
      0100      GNU Licensed. No Warranty  

Cactus version:    4.2.3
Compile date:      Jul 03 2014 (16:37:25)
Run date:          Jul 04 2014 (02:17:51-0500)
Run host: (pid=13059)
Working directory: /home/jtao/Development/HOPAGH/Cactus/exe
Executable:        /home/jtao/Development/HOPAGH/Cactus/exe/cactus_phg
Parameter file:    /home/jtao/Development/HOPAGH/Cactus/exe/phg.par

Activating thorn Cactus...Success -> active implementation Cactus
Activation requested for 
--->phg caphg phgexamples <---
Activating thorn caphg...Success -> active implementation caphg
Activating thorn phg...Success -> active implementation phg
Activating thorn phgexamples...Success -> active implementation phg_examples
Scheduled CaPHG_IterationStart at CUMT_PREITERATION
Scheduled CaPHG_BalanceGrid at CUMT_PREITERATION
Scheduled CaPHG_RefineMesh at CUMT_REFINE_MESH
Scheduled simplest at CUMT_ITERATION
Scheduled estimate_error at CUMT_ESTIMATE_ERROR
  if (recover initial data)
    Recover parameters

  Startup routines
      GROUP CUMT_STARTUP: Cactus Unstructured Mesh Toolkit Start Up
        CaPHG::DRIVER_BANNER: Register CaPH Banner

  Startup routines which need an existing grid hierarchy
      GROUP CUMT_WRAGH: Cactus Unstructured Mesh Toolkit Driver Initialization
  Parameter checking routines
      GROUP CUMT_PARAMCHECK: Parameter Check
        CaPHG::DRIVER_PARAMCHECK: Check CaPHG Parameters

    if (NOT (recover initial data AND recovery_mode is 'strict'))
      Set up grid hierarchy
        GROUP CUMT_BASEGRID: Set Up Coordinate Systems
        GROUP CUMT_INITIAL: Initialize Variables
        GROUP CUMT_POSTINITIAL: Modify Initial Data or Calculate Data That Depend on the Initial Data
      Initialise finer grids recursively
      Restrict from finer grids
        GROUP CUMT_POSTSTEP: After the Main Evolution Step
    if (recover initial data)
        GROUP CUMT_BASEGRID: Set Up Coordinate Systems
    if (checkpoint initial data)
    if (analysis)
        GROUP CUMT_ANALYSIS: Carry Out Analysis
  Output grid variables

  do loop over timesteps
    Change grid hierarchy
    Rotate timelevels
    iteration = iteration+1
    t = t+dt
      GROUP CUMT_PRESTEP: Before the Main Evolution Step
      GROUP CUMT_EVOL: Main Evolution Loop
        while (CaPHG::iterate)
          GROUP CUMTi_ITERATION: The Internal Adaptive Iteration Loop
            GROUP CUMT_PREITERATION: Before the Adaptive Iteration Loop
              CaPHG::CaPHG_IterationStart: Set iteration indicator
              CaPHG::CaPHG_BalanceGrid: Repartition the Mesh to Balance the Load
            GROUP CUMT_ITERATION: Inside the Adaptive Iteration Loop
              PHGExamples::simplest: Simplest testing case for CaPHG
            GROUP CUMT_POSTITERATION: Post the Adaptive Iteration Loop
              GROUP CUMT_ESTIMATE_ERROR: Estimate the Error for Adaptive Iteration
                PHGExamples::estimate_error: estimate error
              GROUP CUMT_REFINE_MESH: Refine the Mesh Based on Estimated Error
                CaPHG::CaPHG_RefineMesh: Set iteration indicator
        end while
    Evolve finer grids recursively
    Restrict from finer grids
      GROUP CUMT_POSTSTEP: After the Main Evolution Step
  if (checkpoint)
  if (analysis)
      GROUP CUMT_ANALYSIS: Carry Out Analysis
  Output grid variables

  Termination routines
      GROUP CUMT_TERMINATE: Cactus Unstructured Mesh Toolkit Terminates
        CaPHG::DRIVER_TERMINATE: Destroy Grid Hierarchy and Remove Grid Variables

  Shutdown routines
      GROUP CUMT_SHUTDOWN: Cactus Unstructured Mesh Toolkit Shut Down

  Routines run after changing the grid hierarchy:
CCTK_RegisterGHExtension: Extension CaPHG gets handle 0
                 o-o     o--o  o  o  o-o 
                /        |   | |  | o 
               O      oo O--o  O--O |  -o 
                \    | | |     |  | o   |
                 o-o o-o-o     o  o  o-o

             Cactus Parallel Hierarchical Grid

                 (c) Copyright The Authors 
                 GNU Licensed. No Warranty

INFO (CaPHG): Successfully imported the mesh file cube4.dat !

CCTKi_SetupGHExtensions: Set up extension for handle 0
Final mesh written to "simplest.vtk".
indicator = 1.20e+01, mem = 8.25MB
Final mesh written to "simplest.vtk".
indicator = 9.16e+00, mem = 9.02MB
Final mesh written to "simplest.vtk".
indicator = 4.96e+00, mem = 9.64MB
Final mesh written to "simplest.vtk".
indicator = 4.48e+00, mem = 10.41MB
Final mesh written to "simplest.vtk".
indicator = 3.55e+00, mem = 11.45MB
Final mesh written to "simplest.vtk".
indicator = 2.85e+00, mem = 14.28MB
Final mesh written to "simplest.vtk".
indicator = 2.21e+00, mem = 17.93MB
Final mesh written to "simplest.vtk".
indicator = 1.56e+00, mem = 23.90MB
Final mesh written to "simplest.vtk".
indicator = 1.06e+00, mem = 24.41MB
Final mesh written to "simplest.vtk".
indicator = 8.29e-01, mem = 43.05MB
Final mesh written to "simplest.vtk".
indicator = 6.49e-01, mem = 68.61MB
Final mesh written to "simplest.vtk".
indicator = 4.82e-01, mem = 106.51MB
Final mesh written to "simplest.vtk".
indicator = 3.22e-01, mem = 159.57MB

The output of could be dumpt to vtk or other formats. I use Visit ( to visualize the output and mesh. The work to add an I/O extenson is planned though not urgent at this point since the default standalone functions works well so far.

Here comes the output using Visit:

enter image description here


simplest.c : Standalone Version

/* Parallel Hierarchical Grid -- an adaptive finite element library.
 * Copyright (C) 2005-2010 State Key Laboratory of Scientific and
 * Engineering Computing, Chinese Academy of Sciences. */

/* This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * Lesser General Public License for more details.
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 * MA  02110-1301  USA */

/* $Id: simplest.c,v 1.115 2012/11/20 08:35:55 zlb Exp $
 * This sample code solves the Helmholtz equation:
 *   -\Delta u + a u = f
 * with Dirichlet or periodic boundary conditions */

#include "phg.h"

# undef ELEMENT

#include <string.h>
#include <math.h>

static FLOAT a = 1.0;

static void
func_u(FLOAT x, FLOAT y, FLOAT z, FLOAT *value)
    *value = Cos(2. * M_PI * x) * Cos(2. * M_PI * y) * Cos(2. * M_PI * z);

static void
func_f(FLOAT x, FLOAT y, FLOAT z, FLOAT *value)
    func_u(x, y, z, value);
    *value = 12. * M_PI * M_PI * *value + a * *value;

static void
build_linear_system(SOLVER *solver, DOF *u_h, DOF *f_h)
    int i, j;
    GRID *g = u_h->g;
    ELEMENT *e;

    assert(u_h->dim == 1);
    ForAllElements(g, e) {
    int N = DofGetNBas(u_h, e); /* number of basises in the element */
    FLOAT A[N][N], rhs[N], buffer[N];
    INT I[N];

    /* compute \int \grad\phi_j \cdot \grad\phi_i making use of symmetry */
    for (i = 0; i < N; i++) {
        I[i] = phgSolverMapE2L(solver, 0, e, i);
        for (j = 0; j <= i; j++) {
        A[j][i] = A[i][j] =
            /* stiffness */
            phgQuadGradBasDotGradBas(e, u_h, j, u_h, i, QUAD_DEFAULT) +
            /* mass */
            a * phgQuadBasDotBas(e, u_h, j, u_h, i, QUAD_DEFAULT);

    /* loop on basis functions */
    for (i = 0; i < N; i++) {
        if (phgDofDirichletBC(u_h, e, i, func_u, buffer, rhs+i,
                    DOF_PROJ_NONE)) {
        phgSolverAddMatrixEntries(solver, 1, I + i, N, I, buffer); 
        else {  /* interior node */
        /* right hand side = \int f * phi_i */
        phgQuadDofTimesBas(e, f_h, u_h, i, QUAD_DEFAULT, rhs + i);
        phgSolverAddMatrixEntries(solver, 1, I + i, N, I, A[i]); 
    phgSolverAddRHSEntries(solver, N, I, rhs);

static void
estimate_error(DOF *u_h, DOF *f_h, DOF *grad_u, DOF *error)
/* compute H1 error indicator */
    GRID *g = u_h->g;
    ELEMENT *e;
    DOF *jump, *residual, *tmp;

    jump = phgQuadFaceJump(grad_u, DOF_PROJ_DOT, NULL, QUAD_DEFAULT);
    tmp = phgDofDivergence(grad_u, NULL, NULL, NULL);
    residual = phgDofGetSameOrderDG(u_h, -1, NULL);
    phgDofCopy(f_h, &residual, NULL, NULL); 
    phgDofAXPY(-a, u_h, &residual);
    phgDofAXPY(1.0, tmp, &residual);
    ForAllElements(g, e) {
    int i;
    FLOAT eta, h;
    FLOAT diam = phgGeomGetDiameter(g, e);
    eta = 0.0;
    /* for each face F compute [grad_u \cdot n] */
    for (i = 0; i < NFace; i++) {
        if (e->bound_type[i] & (DIRICHLET | NEUMANN))
        continue;   /* boundary face */
        h = phgGeomGetFaceDiameter(g, e, i);
        eta += *DofFaceData(jump, e->faces[i]) * h;
    eta = eta * .5 + diam * diam * phgQuadDofDotDof(e, residual, residual,
    *DofElementData(error, e->index) = eta;


main(int argc, char *argv[])
    INT periodicity = 0 /* X_MASK | Y_MASK | Z_MASK */;
    INT mem_max = 2000;
    FLOAT tol = 0.1;
    char *fn = "../test/cube4.dat";
    GRID *g;
    ELEMENT *e;
    DOF *u_h, *f_h, *grad_u, *error, *u;
    SOLVER *solver;
    FLOAT L2error, indicator;
    size_t mem_peak;
    phgOptionsRegisterFloat("a", "Coefficient", &a);
    phgOptionsRegisterFloat("tol", "Tolerance", &tol);
    phgOptionsRegisterInt("mem_max", "Maximum memory (MB)", &mem_max);
    phgOptionsRegisterInt("periodicity", "Set periodicity", &periodicity);

    phgInit(&argc, &argv);

    g = phgNewGrid(-1);
    phgSetPeriodicity(g, periodicity);

    if (!phgImport(g, fn, FALSE))
    phgError(1, "can't read file \"%s\".\n", fn);

    /* The discrete solution */
    if (FALSE) {
    /* u_h is h-p type */
    HP_TYPE *hp = phgHPNew(g, HP_HB);
    ForAllElements(g, e)
    e->hp_order = DOF_DEFAULT->order + GlobalElement(g, e->index) % 3;
    phgHPSetup(hp, FALSE);
    u_h = phgHPDofNew(g, hp, 1, "u_h", DofInterpolation);
    else {
    /* u_h is non h-p type */
    u_h = phgDofNew(g, DOF_DEFAULT, 1, "u_h", DofInterpolation);
    phgDofSetDataByValue(u_h, 0.0);

    /* RHS function */
    f_h = phgDofNew(g, DOF_DEFAULT, 1, "f_h",  func_f);

    /* DOF for storing a posteriori error estimates */
    error = phgDofNew(g, DOF_P0, 1, "error indicator", DofNoAction);

    /* The analytic solution */
    u = phgDofNew(g, DOF_ANALYTIC, 1, "u", func_u);

    while (TRUE) {
    double t0 = phgGetTime(NULL);
    if (phgBalanceGrid(g, 1.2, 1, NULL, 0.))
        phgPrintf("Repartition mesh, load imbalance: %lg\n",
    phgPrintf("%"dFMT" DOF, %"dFMT" elements, %"dFMT
          " submeshes, load imbalance: %lg\n",
            DofGetDataCountGlobal(u_h), g->nleaf_global, g->nprocs,
    solver = phgSolverCreate(SOLVER_DEFAULT, u_h, NULL);
    phgPrintf("  DOF: %"dFMT", unknowns: %"dFMT
          ", Dirichlet bdry: %"dFMT"\n",
        DofGetDataCountGlobal(u_h), solver->rhs->map->nglobal,
    build_linear_system(solver, u_h, f_h);
    phgSolverSolve(solver, TRUE, u_h, NULL);
    phgPrintf("  nits = %d, ", solver->nits);
    grad_u = phgDofGradient(u_h, NULL, NULL, NULL);
    estimate_error(u_h, f_h, grad_u, error);
    indicator = Sqrt(phgDofNormL1Vec(error));
    grad_u = phgDofCopy(u_h, NULL, NULL, NULL);
    phgDofAXPY(-1.0, u, &grad_u);
    L2error = phgDofNormL2(grad_u);
    phgMemoryUsage(g, &mem_peak);
    phgPrintf("L2 error = %0.3le, indicator = %0.3le, mem = %0.2lfMB\n",
            (double)L2error, (double)indicator,
            (double)mem_peak / (1024.0 * 1024.0));
    phgPrintf("  Wall time: %0.3le\n", phgGetTime(NULL) - t0);
    if (indicator < tol || mem_peak >= 1024 * (size_t)mem_max * 1024)
    phgMarkRefine(MARK_DEFAULT, error, Pow(0.8,2), NULL, 0., 1,
            Pow(tol, 2) / g->nleaf_global);

#if 0
    phgPrintf("Final mesh written to \"%s\".\n",
    phgExportDX(g, "simplest.dx", u_h, error, NULL));
#elif 1
    phgPrintf("Final mesh written to \"%s\".\n",
    phgExportVTK(g, "simplest.vtk", u_h, error, NULL));




    return 0;

Written with StackEdit.