itaps-parallel updated codes
Mark Miller
miller86 at llnl.gov
Mon Dec 17 09:51:39 CST 2007
Hi All,
I did not get as far on this over weekend as I liked. I did take a minute to make the code
I wrote a little better documented and changed terminology in code to patch parallel
termonilogy a bit better (e.g. use 'part' instead of 'block'). I have a strategy to make
this code simulate reading the parallel decomposed mesh such that each processor
gets a different piece by reading all parts onto all procs and then removing the parts
on each proc that are not owned by that proc. However, the rmv method does not
seem to be working.
Mark
--
Mark C. Miller, Lawrence Livermore National Laboratory
email: mailto:miller86 at llnl.gov
(M/T/W) (925)-423-5901 (!!LLNL BUSINESS ONLY!!)
(Th/F) (530)-753-8511 (!!LLNL BUSINESS ONLY!!)
-------------- next part --------------
#include <mpi.h>
#include <iMesh.h>
#include <iBase.h>
#include <stdarg.h>
#include <stdio.h>
// end of list
#define EoL (void*)-1
// no list
#define NoL (0,EoL)
// supports passing variable length list of void pointers that if non-zero
// need to be free'd.
static void
ITAPSErrorCleanupHelper(int dummy, ...)
{
va_list ap;
va_start(ap, dummy);
void *ptr = va_arg(ap, void*);
while (ptr != EoL)
{
if (ptr)
free(ptr);
ptr = va_arg(ap, void*);
}
va_end(ap);
}
// ****************************************************************************
// Macro: CheckITAPSError2
//
// Purpose: Very useful macro for checking error condition of ITAPS
// implementation after a call into the implementation. This macro does work
// to get an error description, issues a VisIt warning, keeps track of how
// many times a given warning has been issued and suppresses above 5 aborts
// (via 'goto') to the end of the function from which it is called and
// ensures that all pointers allocated prior to the abort get freed.
//
//
// ****************************************************************************
#define CheckITAPSError2(IMI, ERR, ARGS, THELINE, THEFILE) \
if (ERR != 0) \
{ \
int i; \
char msg[1024]; \
char desc[256]; \
for (i = 0; i < sizeof(desc); i++) desc[i] = '\0'; \
int dummyError = ERR; \
iMesh_getDescription(IMI, desc, &dummyError, sizeof(desc)); \
snprintf(msg, sizeof(msg), "%d: Encountered ITAPS error (%d) at line %d in file \"%s\"\n" \
"The description is...\n" \
" \"%s\"\n", commRank, ERR, THELINE, THEFILE, desc); \
fprintf(stderr, "%s\n", msg); \
ITAPSErrorCleanupHelper ARGS; \
goto funcEnd; \
}
#define CheckITAPSError(IMI, ARGS) CheckITAPSError2(IMI, iMerr, ARGS, __LINE__, __FILE__)
int main(int argc, char **argv)
{
int i, j, k;
char *imeshFileName = "./myTest.h5m";
int commSize = 1;
int commRank = 0;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &commSize);
MPI_Comm_rank(MPI_COMM_WORLD, &commRank);
/* Create a new, empty iMesh instance to laod stuff into */
iMesh_Instance iMi;
int iMerr;
iMesh_newMesh(0, &iMi, &iMerr, 0);
iBase_EntitySetHandle iMroot;
iMesh_getRootSet(iMi, &iMroot, &iMerr);
CheckITAPSError(iMi, NoL);
/* Load the entire contents of the h5m file. It would be nice if we
could load only the EntSet (not entities themselves, just ent sets)
hierarchy here, onto every proc, and then follo that up with another
load to load ALL the stuff rooted at a particular ent set in the
hierarchy. Alternatively, if we could specify tag and/or tag values
to indicate which portions of the ent set hierarchy to load onto this
processor in this call, that would be very helpful. Still, another
thing might be useful is to 'leave' everything on disk here and return
enough state internal to iMesh instance that suceeding calls on it
can do I/O to the file, to get information necessary to satisfy the
query. */
/* bogus: loads whole mesh onto every processor. Its just an example. */
iMesh_load(iMi,
iMroot, /*in const iBase_EntitySetHandle entity_set_handle,*/
imeshFileName, /*in const char *name,*/
0, /*in const char *options,*/
&iMerr, /*out int *err,*/
strlen(imeshFileName)+2,/*in int name_len,*/
0 /*in int options_len*/);
/* Get information on the tags we'll use for parts and bcs */
iBase_TagHandle partTag, bcTag;
iMesh_getTagHandle(iMi, "part", &partTag, &iMerr, 7);
CheckITAPSError(iMi, NoL);
iMesh_getTagHandle(iMi, "bogus_bc", &bcTag, &iMerr, 10);
CheckITAPSError(iMi, NoL);
/* Get the top-level list of ent sets from the hierarchy */
/* Why is this returning many more sets than I 'rooted' in root in the
silo_to_itaps tool that I used to create it? It seems to return a
number equal to the number of times I called createEntSet plus 2
(the latter 2 being I think universe and empty sets). I was expecting
it to return a number of sets equal to the number I actually rooted
at the root set via addEntSet calls. */
iBase_EntityHandle *top_sets = 0;
int top_sets_allocated = 0;
int top_sets_occupied = 0;
iMesh_getEntSets(iMi, iMroot, 1, &top_sets, &top_sets_allocated, &top_sets_occupied, &iMerr);
CheckITAPSError(iMi, NoL);
if (commRank == 0)
printf("top set count = %d\n", top_sets_occupied);
/* Ok, paw through (Carl's jargon) this collection of entities, looking
for ones that are the 'part' parts (e.g. have partTag). I use
failure to find tag as indication that it is not present. Alternatively,
I could query for all tags defined on the ent set, top_sets[k]. */
int nBlocks = 0;
iBase_EntityHandle part_sets[256];
for (k = 0; k < top_sets_occupied; k++)
{
int partId;
iMesh_getIntData(iMi, top_sets[k], partTag, &partId, &iMerr);
if (iMerr == 0)
part_sets[nBlocks++] = top_sets[k];
else if (iMerr != iBase_TAG_NOT_FOUND)
CheckITAPSError(iMi, NoL);
}
if (commRank == 0)
printf("part count = %d\n", nBlocks);
/* Ok, make a second pass and look for ent sets with boundary condition
(bc) tag. Note, given how I originally constructed this, I would have
expected to have to do this by drilling down underneath the 'part'
ent sets I created instead of working directly from top-level sets. */
int nBCsets = 0;
iBase_EntityHandle bc_sets[256];
for (k = 0; k < top_sets_occupied; k++)
{
int bcId;
iMesh_getIntData(iMi, top_sets[k], bcTag, &bcId, &iMerr);
if (iMerr == 0)
bc_sets[nBCsets++] = top_sets[k];
else if (iMerr != iBase_TAG_NOT_FOUND)
CheckITAPSError(iMi, NoL);
}
/* Divide top-level parts over processors */
/* If commSize does not divide nBlocks evenly, add remainder parts
such that proc 0...#remainder get an extra part */
int baseBlocksPerProc = nBlocks / commSize;
int remainderBlocks = nBlocks % commSize;
int myStart, myEnd;
if (commRank < remainderBlocks)
{
myStart = commRank * (baseBlocksPerProc+1);
myEnd = myStart + baseBlocksPerProc+1;
}
else
{
myStart = commRank * baseBlocksPerProc + remainderBlocks;
myEnd = myStart + baseBlocksPerProc;
}
/* ok, to simulate what it would really be like working with
this parallel mesh, we can remove parts from this instance
that we don't own on this processor */
for (i = 0; i < nBlocks; i++)
{
if (i >= myStart && i < myEnd)
continue;
/* This call fails to remove set due to entity not found error */
iMesh_rmvEntSet(iMi, part_sets[i], &iMroot, &iMerr);
CheckITAPSError(iMi, NoL);
}
/* 1) find all parts adjacent to a given part
2) find all parts adjacent to the bc
3) find all ents adjacent to bc and collect them to
a processor
4) find all face entities shared between parts
5) find all node entities shared by 2, 3, 4, 8 parts
6) find all parts that have things on the 'external boundary'
*/
if (nBCsets > 0)
{
for (k = 0; k < nBCsets; k++)
{
int bcId;
iMesh_getIntData(iMi, bc_sets[k], bcTag, &bcId, &iMerr);
if (myStart <= bcId && bcId < myEnd)
{
/* get hex ents adjacent to the bc nodes */
iBase_EntityHandle *ents_adj_to_bc_nodes = 0;
int adj_allocated = 0;
int adj_occupied = 0;
int *offsets = 0;
int offsets_allocated = 0;
int offsets_occupied = 0;
int *in_es = 0;
int in_es_alloc = 0;
int in_es_size = 0;
iMesh_getAdjEntities(iMi,
bc_sets[k], /*in const iBase_EntityHandle entity_set_handle,*/
iBase_VERTEX, /*in const int entity_type_requestor, */
iMesh_POINT, /*in const int entity_topology_requestor,*/
iBase_REGION, /*in const int entity_type_requested,*/
&ents_adj_to_bc_nodes, /*inout iBase_EntityHandle** adj_entity_handles,*/
&adj_allocated, /*inout int* adj_entity_handles_allocated,*/
&adj_occupied, /*out int* adj_entity_handles_size,*/
&offsets, /*inout int** offset,*/
&offsets_allocated, /*inout int* offset_allocated,*/
&offsets_occupied, /*out int* offset_size,*/
&in_es, /*inout int** in_entity_set,*/
&in_es_alloc, /*inout int* in_entity_set_allocated,*/
&in_es_size, /*out int* in_entity_set_size,*/
&iMerr);
CheckITAPSError(iMi, (0,ents_adj_to_bc_nodes,offsets,in_es,EoL));
printf("on processor %d, not %d ents adjacent to bc nodes for id %d\n",
commRank, adj_occupied, bcId);
}
}
}
for (k = myStart; k < myEnd; k++)
{
int hex=-100;
int point=-200;
/* this returns zero for ALL_TOPOLOGIES */
iMesh_getNumOfTopo(iMi, part_sets[k], iMesh_HEXAHEDRON, &hex, &iMerr);
iMesh_getNumOfTopo(iMi, part_sets[k], iMesh_POINT, &point, &iMerr);
printf("processor %d reading part %d, has %d,%d points,hexes\n", commRank, k, point,hex);
CheckITAPSError(iMi, NoL);
}
funcEnd:
MPI_Finalize();
return 0;
}
-------------- next part --------------
#include <silo.h>
#include <iMesh.h>
#include <iBase.h>
#include <stdarg.h>
// end of list
#define EoL (void*)-1
// no list
#define NoL (0,EoL)
// supports passing variable length list of void pointers that if non-zero
// need to be free'd.
static void
ITAPSErrorCleanupHelper(int dummy, ...)
{
va_list ap;
va_start(ap, dummy);
void *ptr = va_arg(ap, void*);
while (ptr != EoL)
{
if (ptr)
free(ptr);
ptr = va_arg(ap, void*);
}
va_end(ap);
}
// ****************************************************************************
// Macro: CheckITAPSError2
//
// Purpose: Very useful macro for checking error condition of ITAPS
// implementation after a call into the implementation. This macro does work
// to get an error description, issues a warning, and aborts
// (via 'goto') to the end of the function from which it is called and
// ensures that all pointers listed in ARGS get freed.
//
// ****************************************************************************
#define CheckITAPSError2(IMI, ERR, ARGS, THELINE, THEFILE) \
if (ERR != 0) \
{ \
int i; \
char msg[1024]; \
char desc[256]; \
for (i = 0; i < sizeof(desc); i++) desc[i] = '\0'; \
int dummyError = ERR; \
iMesh_getDescription(IMI, desc, &dummyError, sizeof(desc)); \
snprintf(msg, sizeof(msg), "Encountered ITAPS error (%d) at line %d in file \"%s\"\n" \
"The description is...\n" \
" \"%s\"\n", ERR, THELINE, THEFILE, desc); \
fprintf(stderr, "%s\n", msg); \
ITAPSErrorCleanupHelper ARGS; \
goto funcEnd; \
}
#define CheckITAPSError(IMI, ARGS) CheckITAPSError2(IMI, iMerr, ARGS, __LINE__, __FILE__)
int main(int argc, char **argv)
{
int i, j, k;
char *siloFileName = "/usr/gapps/visit/data/multi_ucd3d.silo";
char *imeshFileName = "./myTest.h5m";
/* tell Silo to load up data in single precision only */
DBForceSingle(1);
/* open the silo file for input data */
DBfile *sfile = DBOpen(siloFileName, DB_UNKNOWN, DB_READ);
/* create the iMesh instance and get its root set */
/* Set-theoretically speacking, is the root set the same notion as 'universe' */
iMesh_Instance iMi;
int iMerr;
iMesh_newMesh(0, &iMi, &iMerr, 0);
iBase_EntitySetHandle iMroot;
iMesh_getRootSet(iMi, &iMroot, &iMerr);
CheckITAPSError(iMi, NoL);
/* create tags well use to specify parts of the mesh
and boundary conditions (bc) */
iBase_TagHandle partTag, bcTag;
iMesh_createTag(iMi, "part", 1, iBase_INTEGER, &partTag, &iMerr, 7);
iMesh_createTag(iMi, "bogus_bc", 1, iBase_INTEGER, &bcTag, &iMerr, 10);
/* get the Silo multimesh object: a list of mesh pieces */
/* Note: For now, this code assumes all of Silo's mesh pieces are
composed entirely of Hex elements. This is not true of Silo in
general but is true of the particular Silo file we are reading
here. */
DBmultimesh *mm = DBGetMultimesh(sfile, "mesh1");
/* iterate over pieces of mesh from Silo */
for (k = 0; k < mm->nblocks; k++)
{
char tmpName[256];
/* Create EntSet for the part of mesh */
/* ent/es creation always seems to occur in a global space and then
objects are assigned to sets in that space as needed */
iBase_EntityHandle partes;
iMesh_createEntSet(iMi, 0, &partes, &iMerr);
CheckITAPSError(iMi, NoL);
/* Tag the part ent set with partTag value of part number, k */
iMesh_setEntSetIntData(iMi, partes, partTag, k, &iMerr);
CheckITAPSError(iMi, NoL);
/* Add the entset handle for this part to root set */
iMesh_addEntSet(iMi, partes, &iMroot, &iMerr);
/* Build name of the Silo object to read */
sprintf(tmpName, "block%d/mesh1", k);
printf ("processing %s\n", tmpName);
/* Read this block of mesh from Silo */
DBucdmesh *um = DBGetUcdmesh(sfile, tmpName);
/* Convert three sep. coord arrays Silo returns to a single array
of interleaved x,y,z values for iMesh. Also, keep track of which
nodes, if any, are on the z=0 bc. */
double *vtx_coords = (double *) malloc(3*um->nnodes*sizeof(double));
int *bc_node_ids = (int *) malloc(um->nnodes*sizeof(int));
int bc_nodes = 0;
j = 0;
for (i = 0; i < um->nnodes; i++)
{
vtx_coords[j++] = ((float*) um->coords[0])[i];
vtx_coords[j++] = ((float*) um->coords[1])[i];
vtx_coords[j++] = ((float*) um->coords[2])[i];
if ((um->coords[2])[i] <= 0.0)
bc_node_ids[bc_nodes++] = i;
}
/* Create iMesh vertex entities and assign coord values for them.
By setting vtx_ents to zero, iMesh allocates for us. I forget,
do we require all implementations to use malloc for the memory
allocation? If not, we need some way of knowing so that we can
later free with either malloc or delete []. I guess since its
a C interface, we can safely assume malloc is used? */
/* these seem to be created in some kind of a 'global' space.
How do I create vertex in an entity set */
iBase_EntityHandle *vtx_ents = 0;
int vtx_allocated = 0;
int vtx_occupied = 0;
iMesh_createVtxArr(iMi, um->nnodes, /*in const int num_verts, */
iBase_INTERLEAVED, /*in const int storage_order, */
vtx_coords, /*in const double* new_coords, */
3*um->nnodes, /*in const int new_coords_size, */
&vtx_ents, /*inout iBase_EntityHandle** new_vertex_handles,*/
&vtx_allocated, /*inout int* new_vertex_handles_allocated, */
&vtx_occupied, /*inout int* new_vertex_handles_size, */
&iMerr);
CheckITAPSError(iMi, (0,vtx_ents,EoL));
/* Ok, at this point, iMesh knows about the verts in terms of
iBase_EntityHandles. Now, we need to tell iMesh about the
hex elements of Silo's mesh block too. Above, we assume that
the vtx_ents array returned is 1:1 with the vtx_coords array
passed in. Is this a requirement of all iMesh implementations?
If not, we'd have an additional step to figure out which
iBase_EntityHandle in vtx_ents goes with which 'node' in
vtx_coords passed in.
Each of the hex elements in Silo's mesh is specified by 8
integer references into vtx_coords. So, I simply create a
analgous array of iMesh entity handles. I also assume that
Silo's node order over a hex is equivalent to iMesh. However,
I don't recall seeing anywhere where this node order is either
specified or an iMesh client is supposed to specify it. */
/* create refs to vtx ents for hexes */
iBase_EntityHandle *vtx_ent_refs = (iBase_EntityHandle *) malloc(um->zones->lnodelist*sizeof(iBase_EntityHandle));
for (i = 0; i < um->zones->lnodelist; i++)
vtx_ent_refs[i] = vtx_ents[um->zones->nodelist[i]];
/* Create all the hexes for this part in iMesh */
iBase_EntityHandle *hex_ents = 0;
int hex_allocated = 0;
int hex_occupied = 0;
int *hex_status = 0;
int hex_status_alloc = 0;
int hex_status_occ = 0;
iMesh_createEntArr(iMi,
iMesh_HEXAHEDRON, /*in const int new_entity_topology,*/
vtx_ent_refs, /*in const iBase_EntityHandle* lower_order_entity_handles,*/
um->zones->lnodelist, /*in const int lower_order_entity_handles_size,*/
&hex_ents, /*out iBase_EntityHandle** new_entity_handles,*/
&hex_allocated, /*out int* new_entity_handles_allocated,*/
&hex_occupied, /*out int* new_entity_handles_size,*/
&hex_status, /*inout int** status,*/
&hex_status_alloc, /*inout int* status_allocated,*/
&hex_status_occ, /*out int* status_size,*/
&iMerr);
CheckITAPSError(iMi, (0,vtx_ent_refs,hex_ents,hex_status,EoL));
/* Now, add the hexes created above to the part ent set */
iMesh_addEntArrToSet(iMi, hex_ents, hex_occupied, &partes, &iMerr);
CheckITAPSError(iMi, NoL);
/* Create entity set for any bc nodes we might have */
if (bc_nodes > 0)
{
iBase_EntityHandle bc_nodes_es;
iMesh_createEntSet(iMi, 0, &bc_nodes_es, &iMerr);
CheckITAPSError(iMi, NoL);
/* tag the bc_nodes_es ent set with bcTag value of part number */
iMesh_setEntSetIntData(iMi, bc_nodes_es, bcTag, k, &iMerr);
CheckITAPSError(iMi, NoL);
/* add an es handle for these bc nodes to the part ent set */
iMesh_addEntSet(iMi, bc_nodes_es, &partes, &iMerr);
CheckITAPSError(iMi, NoL);
/* Add the node entities to this bc_nodes_es */
iBase_EntityHandle *bc_node_ents = (iBase_EntityHandle *) malloc(bc_nodes*sizeof(iBase_EntityHandle));
for (i = 0; i < bc_nodes; i++)
bc_node_ents[i] = vtx_ents[bc_node_ids[i]];
iMesh_addEntArrToSet(iMi, bc_node_ents, bc_nodes, &bc_nodes_es, &iMerr);
CheckITAPSError(iMi, (0,bc_node_ents,EoL));
free(bc_node_ids);
}
free(vtx_ent_refs);
free(hex_ents);
free(hex_status);
free(vtx_coords);
DBFreeUcdmesh(um);
}
iMesh_save(iMi,
iMroot, /*in const iBase_EntitySetHandle entity_set_handle, */
imeshFileName, /*in const char *name,*/
0, /*in const char *options,*/
&iMerr, /*out int *err,*/
strlen(imeshFileName)+2, /*in const int name_len, */
0 /*in int options_len*/);
CheckITAPSError(iMi, NoL);
funcEnd:
return 0;
}
More information about the itaps-parallel
mailing list