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