Cactus on Unstructured Meshes with CaPHG
Jian Tao (jtao@cct.lsu.edu)
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
applications:
- Grid Hierarchy (GH) represents the distributed adaptive GH. The abstraction enables application developers to create, operate and destroy hierarchical grid structures. The regridding and partitioning of a grid structure are done automatically whenever necessary. In Cactus, grid operations are handled by a driver thorn which is a special module in Cactus.
- Grid Function (GF) represents a distributed data structure that represents the variables in an application. Storage, synchronization, arithmetic, and reduction operations are implemented for the GF by standard thorns. The application developers are responsible for providing proper routines for initialization, boundary updates, etc.
- Grid Geometry (GG) represents the coordinates, bounding boxes, and bounding box lists of the computational domain. Operations on the GG, such as union, intersection, refine, and coarsen are usually implemented in a driver thorn as well.
A simple maping between Cactus and PHG would be
- Cactus GH + GG
→ PHG GRID - Cactus GF
→ PHG DOF
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 */
#if ALLOW_CURVED_BOUNDARY
EXPR **bdry_funcs; /**< NULL terminated list of projection
functions for curved boundaries */
#endif
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;
#if USE_MPI
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 */
#endif
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 */
#endif
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 */
} GRID;
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 (http://cactuscode.org/documentation/UsersGuide.pdf), 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.
CaPHGGroup
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;
CaPHGExt
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);
CCTK_RegisterGHExtensionScheduleTraverseGH(GHEx_Handle,
CaPHG_ScheduleTraverseGH);
/* Overload necessary functions */
CCTK_OverloadExit((int (*)(cGH*, int)) CaPHG_Exit);
CCTK_OverloadAbort((int (*)(cGH*, int)) CaPHG_Abort);
CCTK_OverloadMyProc(CaPHG_MyProc);
CCTK_OverloadnProcs(CaPHG_nProcs);
CCTK_OverloadBarrier(CaPHG_Barrier);
CCTK_OverloadEnableGroupStorage(CaPHG_EnableGroupStorage);
CCTK_OverloadQueryGroupStorageB(CaPHG_QueryGroupStorage);
CCTK_OverloadGroupDynamicData(CaPHG_GroupDynamicData);
...
}
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();
MALLOC_SAFE_CALL(
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;
MALLOC_SAFE_CALL(
GHExt->groups[gi].data = malloc(groupdata.numvars * sizeof(void **)));
for (vi = 0; vi < groupdata.numvars; vi++)
{
MALLOC_SAFE_CALL(
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] =
GHExt->groups[gi].data[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),
g.user_func);
g.data[vi][tl] = gf;
} // for (tl = 0; tl < g.numtimelevels; tl++)
}
/* TODO: this shall be checked at the parsing stage. */
else
{
CCTK_ERROR(
"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!");
}
else
{
for (vi = 0; vi < g.numvars; vi++)
{
for (tl = 0; tl < g.numtimelevels; tl++)
{
gvi = g.first_var + vi;
MALLOC_SAFE_CALL(
g.data[vi][tl] = malloc(
g.numtimelevels * CCTK_VarTypeSize(CCTK_VarTypeI(gvi))));
}
}
}
}
else
{
CCTK_ERROR("CaPHG now only support CCTK_GF and CCTK_SCALAR");
}
}
g.storage = 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);
}
else
{
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 */
GHExt = CCTK_GHExtension(GH, DRIVERNAME);
/* 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:
interface.ccl
implements: phg_examples
inherits: caphg
public:
CCTK_DOF main type=GF timelevels=1 tags='dof_type="DOF_DEFAULT" solver="SOLVER_DEFAULT" dof_user_func="DofInterpolation"'
{
u_h
} "main unknowns of a PDE"
CCTK_DOF analytic_solution type=GF timelevels=1 tags='dof_type="DOF_ANALYTIC" dof_user_func="myu"'
{
u
} "analytical solution"
CCTK_DOF rhs_function type=GF timelevels=1 tags='dof_type="DOF_DEFAULT" dof_user_func="myf"'
{
f_h
} "right hand side function"
schedule.ccl
STORAGE: main analytic_solution rhs_function
schedule simplest in CUMT_ITERATION
{
LANG: C
} "Simplest testing case for CaPHG"
schedule estimate_error IN CUMT_ESTIMATE_ERROR
{
LANG: C
} "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)
{
DECLARE_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);
phgSolverDestroy(&solver);
phgPrintf("Final mesh written to \"%s\".\n",
phgExportVTK(cctkGH->grid, "simplest.vtk", u_h, error, NULL));
return;
}
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, GridFuncStuff.pl 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
DIM_OF_WORLD: 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:
0
0
...
0
0
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
--------------------------------------------------------------------------------
10
1 0101 ************************
01 1010 10 The Cactus Code V4.2.3
1010 1101 011 www.cactuscode.org
1001 100101 ************************
00010101
100011 (c) Copyright The Authors
0100 GNU Licensed. No Warranty
0101
--------------------------------------------------------------------------------
Cactus version: 4.2.3
Compile date: Jul 03 2014 (16:37:25)
Run date: Jul 04 2014 (02:17:51-0500)
Run host: jtao1-2.lsu.edu (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 CUMT_STARTUP at CCTK_STARTUP
Scheduled DRIVER_BANNER at CUMT_STARTUP
Scheduled CUMT_WRAGH at CCTK_WRAGH
Scheduled CUMT_PARAMCHECK at CCTK_PARAMCHECK
Scheduled DRIVER_PARAMCHECK at CUMT_PARAMCHECK
Scheduled CUMT_BASEGRID at CCTK_BASEGRID
Scheduled CUMT_INITIAL at CCTK_INITIAL
Scheduled CUMT_POSTINITIAL at CCTK_POSTINITIAL
Scheduled CUMT_PRESTEP at CCTK_PRESTEP
Scheduled CUMT_EVOL at CCTK_EVOL
Scheduled CUMTi_ITERATION at CUMT_EVOL
Scheduled CUMT_PREITERATION at CUMTi_ITERATION
Scheduled CUMT_ITERATION at CUMTi_ITERATION
Scheduled CUMT_POSTITERATION at CUMTi_ITERATION
Scheduled CaPHG_IterationStart at CUMT_PREITERATION
Scheduled CaPHG_BalanceGrid at CUMT_PREITERATION
Scheduled CUMT_ESTIMATE_ERROR at CUMT_POSTITERATION
Scheduled CUMT_REFINE_MESH at CUMT_POSTITERATION
Scheduled CaPHG_RefineMesh at CUMT_REFINE_MESH
Scheduled CUMT_POSTSTEP at CCTK_POSTSTEP
Scheduled CUMT_ANALYSIS at CCTK_ANALYSIS
Scheduled CUMT_TERMINATE at CCTK_TERMINATE
Scheduled DRIVER_TERMINATE at CUMT_TERMINATE
Scheduled CUMT_SHUTDOWN at CCTK_SHUTDOWN
Scheduled simplest at CUMT_ITERATION
Scheduled estimate_error at CUMT_ESTIMATE_ERROR
--------------------------------------------------------------------------------
if (recover initial data)
Recover parameters
endif
Startup routines
[CCTK_STARTUP]
GROUP CUMT_STARTUP: Cactus Unstructured Mesh Toolkit Start Up
CaPHG::DRIVER_BANNER: Register CaPH Banner
Startup routines which need an existing grid hierarchy
[CCTK_WRAGH]
GROUP CUMT_WRAGH: Cactus Unstructured Mesh Toolkit Driver Initialization
Parameter checking routines
[CCTK_PARAMCHECK]
GROUP CUMT_PARAMCHECK: Parameter Check
CaPHG::DRIVER_PARAMCHECK: Check CaPHG Parameters
Initialisation
if (NOT (recover initial data AND recovery_mode is 'strict'))
[CCTK_PREREGRIDINITIAL]
Set up grid hierarchy
[CCTK_POSTREGRIDINITIAL]
[CCTK_BASEGRID]
GROUP CUMT_BASEGRID: Set Up Coordinate Systems
[CCTK_INITIAL]
GROUP CUMT_INITIAL: Initialize Variables
[CCTK_POSTINITIAL]
GROUP CUMT_POSTINITIAL: Modify Initial Data or Calculate Data That Depend on the Initial Data
Initialise finer grids recursively
Restrict from finer grids
[CCTK_POSTRESTRICTINITIAL]
[CCTK_POSTPOSTINITIAL]
[CCTK_POSTSTEP]
GROUP CUMT_POSTSTEP: After the Main Evolution Step
endif
if (recover initial data)
[CCTK_BASEGRID]
GROUP CUMT_BASEGRID: Set Up Coordinate Systems
[CCTK_RECOVER_VARIABLES]
[CCTK_POST_RECOVER_VARIABLES]
endif
if (checkpoint initial data)
[CCTK_CPINITIAL]
endif
if (analysis)
[CCTK_ANALYSIS]
GROUP CUMT_ANALYSIS: Carry Out Analysis
endif
Output grid variables
do loop over timesteps
[CCTK_PREREGRID]
Change grid hierarchy
[CCTK_POSTREGRID]
Rotate timelevels
iteration = iteration+1
t = t+dt
[CCTK_PRESTEP]
GROUP CUMT_PRESTEP: Before the Main Evolution Step
[CCTK_EVOL]
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
[CCTK_POSTRESTRICT]
[CCTK_POSTSTEP]
GROUP CUMT_POSTSTEP: After the Main Evolution Step
if (checkpoint)
[CCTK_CHECKPOINT]
endif
if (analysis)
[CCTK_ANALYSIS]
GROUP CUMT_ANALYSIS: Carry Out Analysis
endif
Output grid variables
enddo
Termination routines
[CCTK_TERMINATE]
GROUP CUMT_TERMINATE: Cactus Unstructured Mesh Toolkit Terminates
CaPHG::DRIVER_TERMINATE: Destroy Grid Hierarchy and Remove Grid Variables
Shutdown routines
[CCTK_SHUTDOWN]
GROUP CUMT_SHUTDOWN: Cactus Unstructured Mesh Toolkit Shut Down
Routines run after changing the grid hierarchy:
[CCTK_POSTREGRID]
--------------------------------------------------------------------------------
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
--------------------------------------------------------------------------------
Done.
I/O
The output of could be dumpt to vtk or other formats. I use Visit (https://wci.llnl.gov/codes/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:
Appendix
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
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* 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"
#if (PHG_VERSION_MAJOR <= 0 && PHG_VERSION_MINOR < 9)
# undef ELEMENT
typedef SIMPLEX ELEMENT;
#endif
#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);
phgDofFree(&tmp);
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,
QUAD_DEFAULT);
*DofElementData(error, e->index) = eta;
}
phgDofFree(&jump);
phgDofFree(&residual);
return;
}
int
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;
phgSetVerbosity(1);
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);
phgHPFree(&hp);
}
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",
(double)g->lif);
phgPrintf("%"dFMT" DOF, %"dFMT" elements, %"dFMT
" submeshes, load imbalance: %lg\n",
DofGetDataCountGlobal(u_h), g->nleaf_global, g->nprocs,
(double)g->lif);
solver = phgSolverCreate(SOLVER_DEFAULT, u_h, NULL);
phgPrintf(" DOF: %"dFMT", unknowns: %"dFMT
", Dirichlet bdry: %"dFMT"\n",
DofGetDataCountGlobal(u_h), solver->rhs->map->nglobal,
solver->rhs->map->bdry_nglobal);
build_linear_system(solver, u_h, f_h);
phgSolverSolve(solver, TRUE, u_h, NULL);
phgPrintf(" nits = %d, ", solver->nits);
phgSolverDestroy(&solver);
grad_u = phgDofGradient(u_h, NULL, NULL, NULL);
estimate_error(u_h, f_h, grad_u, error);
indicator = Sqrt(phgDofNormL1Vec(error));
phgDofFree(&grad_u);
grad_u = phgDofCopy(u_h, NULL, NULL, NULL);
phgDofAXPY(-1.0, u, &grad_u);
L2error = phgDofNormL2(grad_u);
phgDofFree(&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)
break;
phgMarkRefine(MARK_DEFAULT, error, Pow(0.8,2), NULL, 0., 1,
Pow(tol, 2) / g->nleaf_global);
phgRefineMarkedElements(g);
}
#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));
#endif
phgDofFree(&u);
phgDofFree(&u_h);
phgDofFree(&f_h);
phgDofFree(&error);
phgFreeGrid(&g);
phgFinalize();
return 0;
}
Written with StackEdit.