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