Adding connectivity information to fix a bug in accumulating to R_int
[charm.git] / src / libs / ck-libs / ParFUM-Tops / ParFUM_TOPS.cc
1 /**
2
3    @file
4    @brief Implementation of ParFUM-TOPS layer, except for Iterators
5
6    @author Isaac Dooley
7
8    @todo add code to generate ghost layers
9    @todo Support multiple models
10
11    @note FEM_DATA+0 holds the elemAttr or nodeAtt data
12    @note FEM_DATA+1 holds the id
13    @note FEM_DATA+2 holds nodal coordinates
14    @note FEM_CONN holds the element connectivity
15
16 */
17
18 #include "ParFUM_TOPS.h"
19 #include "ParFUM.decl.h"
20 #include "ParFUM_internals.h"
21 #ifdef CUDA
22     #include <cuda.h>
23     #include <cuda_runtime.h>
24 #endif
25
26 void setTableReferences(TopModel* model){
27   model->ElemConn_T = &((FEM_IndexAttribute*)model->mesh->elem[0].lookup(FEM_CONN,""))->get();
28   model->node_id_T = &((FEM_DataAttribute*)model->mesh->node.lookup(FEM_DATA+1,""))->getInt();
29   model->elem_id_T = &((FEM_DataAttribute*)model->mesh->elem[0].lookup(FEM_DATA+1,""))->getInt();
30   model->ElemData_T = &((FEM_DataAttribute*)model->mesh->elem[0].lookup(FEM_DATA+0,""))->getChar();
31   model->NodeData_T = &((FEM_DataAttribute*)model->mesh->node.lookup(FEM_DATA+0,""))->getChar();
32 #ifdef FP_TYPE_FLOAT
33   model->coord_T = &((FEM_DataAttribute*)model->mesh->node.lookup(FEM_DATA+2,""))->getFloat();
34 #else
35   model->coord_T = &((FEM_DataAttribute*)model->mesh->node.lookup(FEM_DATA+2,""))->getDouble();
36   model->n2eConn_T = &((FEM_DataAttribute*)model->mesh->elem[0].lookup(FEM_DATA+2,""))->getInt();
37 #endif
38 }
39
40
41 TopModel* topModel_Create_Init(int elem_attr_sz, int node_attr_sz, int model_attr_sz){
42
43   CkAssert(elem_attr_sz > 0);
44   CkAssert(node_attr_sz > 0);
45   TopModel *model = new TopModel;
46   memset(model, 0, sizeof(TopModel));
47   model->elem_attr_size = elem_attr_sz;
48   model->node_attr_size = node_attr_sz;
49   model->model_attr_size = model_attr_sz;
50
51   // This only uses a single mesh, so better not create multiple ones of these
52   int which_mesh=FEM_Mesh_default_write();
53   model->mesh = FEM_Mesh_lookup(which_mesh,"TopModel::TopModel()");
54
55
56 /** @note   Here we allocate the arrays with a single
57             initial node and element, which are set as
58             invalid. If no initial elements were provided.
59             the AllocTable2d's would not ever get allocated,
60             and insertNode or insertElement would fail.
61  */
62   char temp_array[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
63
64   // Allocate node coords
65 #ifdef FP_TYPE_FLOAT
66   FEM_Mesh_data(which_mesh,FEM_NODE,FEM_DATA+2,temp_array, 0, 1, FEM_FLOAT, 3);
67 #else
68   FEM_Mesh_data(which_mesh,FEM_NODE,FEM_DATA+2,temp_array, 0, 1, FEM_DOUBLE, 3);
69 #endif
70
71   // Allocate element connectivity
72   FEM_Mesh_data(which_mesh,FEM_ELEM+0,FEM_CONN,temp_array, 0, 1, FEM_INDEX_0, 4);
73
74   // Allocate element attributes
75   FEM_Mesh_data(which_mesh,FEM_ELEM+0,FEM_DATA+0,temp_array, 0, 1, FEM_BYTE, model->elem_attr_size);
76   // Allocate element Id array
77   FEM_Mesh_data(which_mesh,FEM_ELEM+0,FEM_DATA+1,temp_array, 0, 1, FEM_INT, 1);
78   // Allocate n2e connectivity
79   FEM_Mesh_data(which_mesh,FEM_ELEM+0,FEM_DATA+2,temp_array, 0, 1, FEM_INT, 4);
80
81
82   // Allocate node attributes
83   FEM_Mesh_data(which_mesh,FEM_NODE+0,FEM_DATA+0,temp_array, 0, 1, FEM_BYTE, model->node_attr_size);
84   // Allocate node Id array
85   FEM_Mesh_data(which_mesh,FEM_NODE+0,FEM_DATA+1,temp_array, 0, 1, FEM_INT, 1);
86
87   FEM_Mesh_allocate_valid_attr(which_mesh, FEM_NODE);
88   FEM_Mesh_allocate_valid_attr(which_mesh, FEM_ELEM+0);
89
90   FEM_set_entity_invalid(which_mesh, FEM_NODE, 0);
91   FEM_set_entity_invalid(which_mesh, FEM_ELEM+0, 0);
92
93   setTableReferences(model);
94
95   return model;
96 }
97
98 TopModel* topModel_Create_Driver(int elem_attr_sz, int node_attr_sz, int model_attr_sz, void *mAtt){
99     // This only uses a single mesh, so don't create multiple TopModels of these
100     CkAssert(elem_attr_sz > 0);
101     CkAssert(node_attr_sz > 0);
102     int which_mesh=FEM_Mesh_default_read();
103     TopModel *model = new TopModel;
104     memset(model, 0, sizeof(TopModel));
105
106     model->elem_attr_size = elem_attr_sz;
107     model->node_attr_size = node_attr_sz;
108     model->model_attr_size = model_attr_sz;
109
110     model->mesh = FEM_Mesh_lookup(which_mesh,"TopModel::TopModel()");
111
112     model->mAtt = mAtt;
113
114     model->num_local_elem = model->mesh->elem[0].size();
115     model->num_local_node = model->mesh->node.size();
116
117     setTableReferences(model);
118
119 #if CUDA
120     /** copy number/sizes of nodes and elements to device structure */
121     model->device_model.elem_attr_size = elem_attr_sz;
122     model->device_model.node_attr_size = node_attr_sz;
123     model->device_model.model_attr_size = model_attr_sz;
124     model->device_model.num_local_node = model->num_local_node;
125     model->device_model.num_local_elem = model->num_local_elem;
126
127
128     /** Create n2e connectivity array and copy to device global memory */
129     {
130         FEM_Mesh_create_node_elem_adjacency(which_mesh);
131         FEM_Mesh* mesh = FEM_Mesh_lookup(which_mesh, "topModel_Create_Driver");
132         FEM_DataAttribute * at = (FEM_DataAttribute*) 
133             model->mesh->elem[0].lookup(FEM_DATA+2,"topModel_Create_Driver");
134         int* n2eTable  = at->getInt().getData();
135
136         FEM_IndexAttribute * iat = (FEM_IndexAttribute*) 
137             model->mesh->elem[0].lookup(FEM_CONN,"topModel_Create_Driver");
138         int* connTable  = iat->get().getData();
139
140         int* adjElements;
141         int size;
142         for (int i=0; i<model->num_local_node; ++i) {
143             mesh->n2e_getAll(i, adjElements, size);
144             for (int j=0; j<size; ++j) {
145                 for (int k=0; k<5; ++k) {
146                     if (connTable[4*adjElements[j]+k] == i) {
147                         n2eTable[4*adjElements[j]+k] = j;
148                         break;
149                     }
150                     assert(k != 4);
151                 }
152             }
153             delete[] adjElements;
154         }
155
156         //for (int i=0; i<model->num_local_elem*4; ++i) {
157         //    printf("%d ", connTable[i]);
158         //    if ((i+1)%4 == 0) printf("\n");
159         //}
160         //printf("\n\n");
161         //for (int i=0; i<model->num_local_elem*4; ++i) {
162         //    printf("%d ", n2eTable[i]);
163         //    if ((i+1)%4 == 0) printf("\n");
164         //}
165
166         size = model->num_local_elem * 4 *sizeof(int);
167         cudaMalloc((void**)&(model->device_model.n2eConnDevice), size);
168         cudaMemcpy(model->device_model.n2eConnDevice,n2eTable,size,
169                 cudaMemcpyHostToDevice);
170
171     }
172
173     /** Copy element Attribute array to device global memory */
174     {
175         FEM_DataAttribute * at = (FEM_DataAttribute*) model->mesh->elem[0].lookup(FEM_DATA+0,"topModel_Create_Driver");
176         AllocTable2d<unsigned char> &dataTable  = at->getChar();
177         unsigned char *ElemData = dataTable.getData();
178         int size = dataTable.size()*dataTable.width();
179         cudaMalloc((void**)&(model->device_model.ElemDataDevice), size);
180         cudaMemcpy(model->device_model.ElemDataDevice,ElemData,size,
181                 cudaMemcpyHostToDevice);
182     }
183
184     /** Copy node Attribute array to device global memory */
185     {
186         FEM_DataAttribute * at = (FEM_DataAttribute*) model->mesh->node.lookup(FEM_DATA+0,"topModel_Create_Driver");
187         AllocTable2d<unsigned char> &dataTable  = at->getChar();
188         unsigned char *NodeData = dataTable.getData();
189         int size = dataTable.size()*dataTable.width();
190         cudaMalloc((void**)&(model->device_model.NodeDataDevice), size);
191         cudaMemcpy(model->device_model.NodeDataDevice,NodeData,size,
192                 cudaMemcpyHostToDevice);
193     }
194
195     /** Copy elem connectivity array to device global memory */
196     {
197         FEM_IndexAttribute * at = (FEM_IndexAttribute*) model->mesh->elem[0].lookup(FEM_CONN,"topModel_Create_Driver");
198         AllocTable2d<int> &dataTable  = at->get();
199         int *data = dataTable.getData();
200         int size = dataTable.size()*dataTable.width()*sizeof(int);
201         cudaMalloc((void**)&(model->device_model.ElemConnDevice), size);
202         cudaMemcpy(model->device_model.ElemConnDevice,data,size,
203                 cudaMemcpyHostToDevice);
204     }
205
206
207
208     /** Copy model Attribute to device global memory */
209     {
210       printf("Copying model attribute of size %d\n", model->model_attr_size);
211         cudaMalloc((void**)&(model->device_model.mAttDevice),
212                 model->model_attr_size);
213         cudaMemcpy(model->device_model.mAttDevice,mAtt,model->model_attr_size,
214                 cudaMemcpyHostToDevice);
215     }
216
217 #endif
218
219     return model;
220 }
221
222 /** Copy node attribute array from CUDA device back to the ParFUM attribute */
223 void top_retrieve_node_data(TopModel* m){ 
224 #if CUDA
225   cudaMemcpy(m->NodeData_T->getData(),
226             m->device_model.NodeDataDevice,
227             m->num_local_node * m->node_attr_size,
228             cudaMemcpyDeviceToHost);
229 #endif
230 }
231
232
233 /** Copy element attribute array from CUDA device back to the ParFUM attribute */
234 void top_retrieve_elem_data(TopModel* m){
235 #if CUDA
236   cudaMemcpy(m->ElemData_T->getData(),
237             m->device_model.ElemDataDevice,
238             m->num_local_elem * m->elem_attr_size,
239             cudaMemcpyDeviceToHost);
240 #endif
241 }
242
243
244
245 /** Cleanup a model */
246 void topModel_Destroy(TopModel* m){
247 #if CUDA
248     cudaFree(m->device_model.mAttDevice);
249     cudaFree(m->device_model.NodeDataDevice);
250     cudaFree(m->device_model.ElemDataDevice);
251 #endif
252     delete m;
253 }
254
255
256 TopNode topModel_InsertNode(TopModel* m, double x, double y, double z){
257   int newNode = FEM_add_node_local(m->mesh,false,false,false);
258   (*m->coord_T)(newNode,0)=x;
259   (*m->coord_T)(newNode,1)=y;
260   (*m->coord_T)(newNode,2)=z;
261   return newNode;
262 }
263
264 TopNode topModel_InsertNode(TopModel* m, float x, float y, float z){
265   int newNode = FEM_add_node_local(m->mesh,false,false,false);
266   (*m->coord_T)(newNode,0)=x;
267   (*m->coord_T)(newNode,1)=y;
268   (*m->coord_T)(newNode,2)=z;
269   return newNode;
270 }
271
272
273 /** Set id of a node
274 @todo Make this work with ghosts
275 */
276 void topNode_SetId(TopModel* m, TopNode n, TopID id){
277   CkAssert(n>=0);
278   (*m->node_id_T)(n,0)=id;
279 }
280
281 /** Insert an element */
282 TopElement topModel_InsertElem(TopModel*m, TopElemType type, TopNode* nodes){
283   CkAssert(type ==  TOP_ELEMENT_TET4);
284   int conn[4];
285   conn[0] = nodes[0];
286   conn[1] = nodes[1];
287   conn[2] = nodes[2];
288   conn[3] = nodes[3];
289   int newEl = FEM_add_element_local(m->mesh, conn, 4, 0, 0, 0);
290   return newEl;
291 }
292
293 /** Set id of an element
294 @todo Make this work with ghosts
295 */
296 void topElement_SetId(TopModel* m, TopElement e, TopID id){
297   CkAssert(e>=0);
298   (*m->elem_id_T)(e,0)=id;
299 }
300
301
302
303 /**
304         @brief Set attribute of a node
305
306         The attribute passed in must be a contiguous data structure with size equal to the value node_attr_sz passed into topModel_Create_Driver() and topModel_Create_Init()
307
308         The supplied attribute will be copied into the ParFUM attribute array "FEM_DATA+0". Then ParFUM will own this data. The function topNode_GetAttrib() will return a pointer to the copy owned by ParFUM. If a single material parameter attribute is used for multiple nodes, each node will get a separate copy of the array. Any subsequent modifications to the data will only be reflected at a single node.
309
310         The user is responsible for deallocating parameter d passed into this function.
311
312 */
313 void topNode_SetAttrib(TopModel* m, TopNode n, void* d){
314   CkAssert(n>=0);
315   unsigned char *data = m->NodeData_T->getData();
316   memcpy(data + n*m->node_attr_size, d, m->node_attr_size);
317 }
318
319 /** @brief Set attribute of an element
320 See topNode_SetAttrib() for description
321 */
322 void topElement_SetAttrib(TopModel* m, TopElement e, void* d){
323   CkAssert(e>=0);
324   unsigned char *data = m->ElemData_T->getData();
325   memcpy(data + e*m->elem_attr_size, d, m->elem_attr_size);
326 }
327
328
329 /** @brief Get elem attribute
330 See topNode_SetAttrib() for description
331 */
332 void* topElement_GetAttrib(TopModel* m, TopElement e){
333   if(! m->mesh->elem[0].is_valid_any_idx(e))
334         return NULL;
335   unsigned char *data = m->ElemData_T->getData();
336   return (data + e*m->elem_attr_size);
337 }
338
339 /** @brief Get nodal attribute
340 See topNode_SetAttrib() for description
341 */
342 void* topNode_GetAttrib(TopModel* m, TopNode n){
343   if(! m->mesh->node.is_valid_any_idx(n))
344         return NULL;
345   unsigned char *data = m->NodeData_T->getData();
346   return (data + n*m->node_attr_size);
347 }
348
349
350
351 /**
352         Get node via id
353         @todo Re-implement this function with some hash to make it fast.
354         @note In the explicit FEA example, this is just used during initialization, so speed is not too important.
355         @todo Does not work with ghosts yet.
356 */
357 TopNode topModel_GetNodeAtId(TopModel* m, TopID id){
358   // lookup node via global ID
359   for(int i=0;i<m->node_id_T->size();++i){
360         if((*m->node_id_T)(i,0)==id){
361           return i;
362         }
363   }
364   return -1;
365 }
366
367
368 /**
369         Get elem via id
370  */
371 TopElement topModel_GetElemAtId(TopModel*m,TopID id){
372     // lookup node via global ID
373     for(int i=0;i<m->elem_id_T->size();++i){
374         if( m->mesh->elem[0].is_valid(i) && (*m->elem_id_T)(i,0)==id){
375                   return i;
376         }
377     }
378     return -1;
379 }
380
381
382 TopNode topElement_GetNode(TopModel* m,TopElement e,int idx){
383   CkAssert(e>=0);
384   const AllocTable2d<int> &conn = m->mesh->elem[0].getConn();
385   CkAssert(idx>=0 && idx<conn.width());
386
387   int node = conn(e,idx);
388
389   return conn(e,idx);
390 }
391
392 int topNode_GetId(TopModel* m, TopNode n){
393   CkAssert(n>=0);
394   return (*m->node_id_T)(n,0);
395 }
396
397
398 /** @todo handle ghost nodes as appropriate */
399 int topModel_GetNNodes(TopModel *model){
400   return model->mesh->node.count_valid();
401 }
402
403 /** @todo Fix to return the width of the conn array */
404 int topElement_GetNNodes(TopModel* model, TopElement elem){
405   return 4;
406 }
407
408 /** @todo make sure we are in a getting mesh */
409 void topNode_GetPosition(TopModel*model, TopNode node,double*x,double*y,double*z){
410   CkAssert(node>=0);
411   *x = (*model->coord_T)(node,0);
412   *y = (*model->coord_T)(node,1);
413   *z = (*model->coord_T)(node,2);
414 }
415
416 /** @todo make sure we are in a getting mesh */
417 void topNode_GetPosition(TopModel*model, TopNode node,float*x,float*y,float*z){
418   CkAssert(node>=0);
419   *x = (*model->coord_T)(node,0);
420   *y = (*model->coord_T)(node,1);
421   *z = (*model->coord_T)(node,2);
422 }
423
424 void topModel_Sync(TopModel*m){
425   MPI_Barrier(MPI_COMM_WORLD);
426
427
428   //  CkPrintf("%d: %d local, %d ghost elements\n", FEM_My_partition(), m->mesh->elem[0].size(),m->mesh->elem[0].ghost->size() );
429   //  CkPrintf("%d: %d local, %d ghost valid elements\n", FEM_My_partition(), m->mesh->elem[0].count_valid(),m->mesh->elem[0].ghost->count_valid() );
430
431 }
432
433 void topModel_TestIterators(TopModel*m){
434   CkAssert(m->mesh->elem[0].ghost!=NULL);
435   CkAssert(m->mesh->node.ghost!=NULL);
436
437   int expected_elem_count = m->mesh->elem[0].count_valid() + m->mesh->elem[0].ghost->count_valid();
438   int iterated_elem_count = 0;
439
440   int expected_node_count = m->mesh->node.count_valid() + m->mesh->node.ghost->count_valid();
441   int iterated_node_count = 0;
442
443   int myId = FEM_My_partition();
444
445
446   TopNodeItr* itr = topModel_CreateNodeItr(m);
447   for(topNodeItr_Begin(itr);topNodeItr_IsValid(itr);topNodeItr_Next(itr)){
448         iterated_node_count++;
449         TopNode node = topNodeItr_GetCurr(itr);
450         void* na = topNode_GetAttrib(m,node);
451         CkAssert(na != NULL);
452   }
453
454   TopElemItr* e_itr = topModel_CreateElemItr(m);
455   for(topElemItr_Begin(e_itr);topElemItr_IsValid(e_itr);topElemItr_Next(e_itr)){
456         iterated_elem_count++;
457         TopElement elem = topElemItr_GetCurr(e_itr);
458         void* ea = topElement_GetAttrib(m,elem);
459         CkAssert(ea != NULL);
460   }
461
462   CkAssert(iterated_node_count == expected_node_count);
463   CkAssert(iterated_elem_count==expected_elem_count);
464
465 }
466
467
468 #include "ParFUM_TOPS.def.h"