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