The beginning of a test of the TOPS layer.
[charm.git] / examples / ParFUM / TOPS_simple2D / pgm.C
1 /*
2  Charm++ Finite-Element Framework Program:
3    Performs simple 2D structural simulation on Triangle-style inputs.
4    
5  init: 
6     read node input file
7     pass nodes into ParFUM framework
8     read element input file
9     pass elements into ParFUM framework
10  
11  driver:
12     extract mesh chunk from framework
13     calculate masses of nodes
14     timeloop
15       compute forces within my chunk
16       communicate
17       apply boundary conditions and integrate
18       pass data to NetFEM
19  
20  Among the hideous things about this program are:
21    -Hardcoded material properties, timestep, and boundary conditions
22    -Limited to 2D
23  
24
25  Converted from f90 Structural Materials program by 
26         Orion Sky Lawlor, 2001
27
28  Updated to new ParFUM interface by
29     Isaac Dooley, 2005
30
31  */
32
33 #include "pgm.h"
34 #include "ParFUM_TOPS.h"
35
36 extern "C" void
37 init(void)
38 {
39   CkPrintf("init started\n");
40   double startTime=CmiWallTimer();
41   const char *eleName="xxx.1.ele";
42   const char *nodeName="xxx.1.node";
43   int nPts=0; //Number of nodes
44   vector2d *pts=0; //Node coordinates
45
46   CkPrintf("Reading node coordinates from %s\n",nodeName);
47   //Open and read the node coordinate file
48   {
49     char line[1024];
50     FILE *f=fopen(nodeName,"r");
51     if (f==NULL) die("Can't open node file!");
52     fgets(line,1024,f);
53     if (1!=sscanf(line,"%d",&nPts)) die("Can't read number of points!");
54     pts=new vector2d[nPts];
55     for (int i=0;i<nPts;i++) {
56       int ptNo;
57       if (NULL==fgets(line,1024,f)) die("Can't read node input line!");
58       if (3!=sscanf(line,"%d%lf%lf",&ptNo,&pts[i].x,&pts[i].y)) 
59         die("Can't parse node input line!");
60     }
61     fclose(f);
62   }
63  
64   int nEle=0;
65   connRec *ele=NULL;
66   CkPrintf("Reading elements from %s\n",eleName);
67   //Open and read the element connectivity file
68   {
69     char line[1024];
70     FILE *f=fopen(eleName,"r");
71     if (f==NULL) die("Can't open element file!");
72     fgets(line,1024,f);
73     if (1!=sscanf(line,"%d",&nEle)) die("Can't read number of elements!");
74     ele=new connRec[nEle];
75     for (int i=0;i<nEle;i++) {
76       int elNo;
77       if (NULL==fgets(line,1024,f)) die("Can't read element input line!");
78       if (4!=sscanf(line,"%d%d%d%d",&elNo,&ele[i][0],&ele[i][1],&ele[i][2])) 
79         die("Can't parse element input line!");  
80       ele[i][0]--; //Fortran to C indexing
81       ele[i][1]--; //Fortran to C indexing
82       ele[i][2]--; //Fortran to C indexing
83     }
84     fclose(f);
85   }
86  
87
88
89   int mesh=FEM_Mesh_default_write(); // Tell framework we are writing to the mesh
90
91   CkPrintf("Passing node coords to framework\n");
92
93   /*   Old versions used FEM_Set_node() and FEM_Set_node_data()
94    *   New versions use the more flexible FEM_Set_Data()
95    */
96
97   FEM_Mesh_data(mesh,        // Add nodes to the current mesh
98                 FEM_NODE,        // We are registering nodes
99                 FEM_DATA+0,      // Register the point locations which are normally 
100                                  // the first data elements for an FEM_NODE
101                 (double *)pts,   // The array of point locations
102                 0,               // 0 based indexing
103                 nPts,            // The number of points
104                 FEM_DOUBLE,      // Coordinates are doubles
105                 2);              // Points have dimension 2 (x,y)
106  
107
108   CkPrintf("Passing elements to framework\n");
109
110   /*   Old versions used FEM_Set_elem() and FEM_Set_elem_conn() 
111    *   New versions use the more flexible FEM_Set_Data()
112    */
113
114   FEM_Mesh_data(mesh,      // Add nodes to the current mesh
115                 FEM_ELEM+0,      // We are registering elements with type 0
116                                  // The next type of element could be registered with FEM_ELEM+1
117                 FEM_CONN,        // Register the connectivity table for this
118                                  // data elements for this type of FEM entity
119                 (int *)ele,      // The array of point locations
120                 0,               // 0 based indexing
121                 nEle,            // The number of elements
122                 FEM_INDEX_0,     // We use zero based node numbering
123                 3);              // Elements have degree 3, since triangles are defined 
124                                  // by three nodes
125  
126    
127   delete[] ele;
128   delete[] pts;
129
130   CkPrintf("Finished with init (Reading took %.f s)\n",CmiWallTimer()-startTime);
131
132 }
133
134
135 // A driver() function 
136 // driver() is required in all FEM programs
137 extern "C" void
138 driver(void)
139 {
140   int nnodes,nelems,ignored;
141   int i, myId=FEM_My_partition();
142   myGlobals g;
143   FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
144   
145   
146   int mesh=FEM_Mesh_default_read(); // Tell framework we are reading data from the mesh
147   
148   // Get node data
149   nnodes=FEM_Mesh_get_length(mesh,FEM_NODE); // Get number of nodes
150   g.coord=new vector2d[nnodes];
151   // Get node positions
152   FEM_Mesh_data(mesh, FEM_NODE, FEM_DATA+0, (double*)g.coord, 0, nnodes, FEM_DOUBLE, 2);  
153
154
155   // Get element data
156   nelems=FEM_Mesh_get_length(mesh,FEM_ELEM+0); // Get number of elements
157   g.nnodes=nnodes; g.nelems=nelems;
158   g.conn=new connRec[nelems];
159   g.S11=new double[nelems];
160   g.S22=new double[nelems];
161   g.S12=new double[nelems];
162   // Get connectivity for elements
163   FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, nelems, FEM_INDEX_0, 3);  
164
165
166   //Initialize associated data
167   g.R_net=new vector2d[nnodes]; //Net force
168   g.d=new vector2d[nnodes];//Node displacement
169   g.v=new vector2d[nnodes];//Node velocity
170   g.a=new vector2d[nnodes];//Node acceleration
171   for (i=0;i<nnodes;i++)
172     g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
173
174 //Apply a small initial perturbation to positions
175   for (i=0;i<nnodes;i++) {
176           const double max=1.0e-10/15.0; //Tiny perturbation
177           g.d[i].x+=max*(i&15);
178           g.d[i].y+=max*((i+5)&15);
179   }
180
181   int fid=FEM_Create_simple_field(FEM_DOUBLE,2);
182
183
184
185     
186     /* Do some TOPS style stuff here */
187     TopModel *m = topModel_Create();
188     assert(m!=NULL);
189
190
191     return;
192
193
194
195
196
197
198
199
200
201   //Timeloop
202   if (myId==0)
203     CkPrintf("Entering timeloop\n");
204   int tSteps=5000;
205   double startTime, totalStart;
206   startTime=totalStart=CkWallTimer();
207   for (int t=0;t<tSteps;t++) {
208     if (1) { //Structural mechanics
209
210     //Compute forces on nodes exerted by elements
211         CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,nnodes,nelems,g.S11,g.S22,g.S12);
212
213     //Communicate net force on shared nodes
214         FEM_Update_field(fid,g.R_net);
215
216     //Advance node positions
217         advanceNodes(dt,nnodes,g.coord,g.R_net,g.a,g.v,g.d,0);
218     }
219
220     //Debugging/perf. output
221     double curTime=CkWallTimer();
222     double total=curTime-startTime;
223     startTime=curTime;
224     if (myId==0 && (t%64==0)) {
225             CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);
226             if (0) {
227               CkPrintf("    Triangle 0:\n");
228               for (int j=0;j<3;j++) {
229                     int n=g.conn[0][j];
230                     CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
231                              n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
232               }
233             }
234     }
235     /* perform migration-based load balancing */
236     if (t%1024==0)
237       FEM_Migrate();
238     
239     if (t%1024==0) { //Publish data to the net
240             NetFEM n=NetFEM_Begin(FEM_My_partition(),t,2,NetFEM_POINTAT);
241             
242             NetFEM_Nodes(n,nnodes,(double *)g.coord,"Position (m)");
243             NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
244             NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");
245             
246             NetFEM_Elements(n,nelems,3,(int *)g.conn,"Triangles");
247                 NetFEM_Scalar(n,g.S11,1,"X Stress (pure)");
248                 NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
249                 NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");
250             
251             NetFEM_End(n);
252     }
253   }
254
255   if (myId==0) {
256     double elapsed=CkWallTimer()-totalStart;
257     CkPrintf("Driver finished: average %.6f s/step\n",elapsed/tSteps);
258   }
259 }
260
261
262 // A PUP function to allow for migration and load balancing of mesh partitions.
263 // The PUP function is not needed if no migration or load balancing is desired.
264 void pup_myGlobals(pup_er p,myGlobals *g) 
265 {
266   FEM_Print("-------- called pup routine -------");
267   pup_int(p,&g->nnodes);
268   pup_int(p,&g->nelems);
269   int nnodes=g->nnodes, nelems=g->nelems;
270   if (pup_isUnpacking(p)) {
271     g->coord=new vector2d[nnodes];
272     g->conn=new connRec[nelems];
273     g->R_net=new vector2d[nnodes]; //Net force
274     g->d=new vector2d[nnodes];//Node displacement
275     g->v=new vector2d[nnodes];//Node velocity
276     g->a=new vector2d[nnodes];
277         g->S11=new double[nelems];
278         g->S22=new double[nelems];
279         g->S12=new double[nelems];
280   }
281   pup_doubles(p,(double *)g->coord,2*nnodes);
282   pup_ints(p,(int *)g->conn,3*nelems);
283   pup_doubles(p,(double *)g->R_net,2*nnodes);
284   pup_doubles(p,(double *)g->d,2*nnodes);
285   pup_doubles(p,(double *)g->v,2*nnodes);
286   pup_doubles(p,(double *)g->a,2*nnodes);
287   pup_doubles(p,(double *)g->S11,nelems);
288   pup_doubles(p,(double *)g->S22,nelems);
289   pup_doubles(p,(double *)g->S12,nelems);
290   if (pup_isDeleting(p)) {
291     delete[] g->coord;
292     delete[] g->conn;
293     delete[] g->R_net;
294     delete[] g->d;
295     delete[] g->v;
296     delete[] g->a;
297         delete[] g->S11;
298         delete[] g->S22;
299         delete[] g->S12;
300   }
301 }
302
303
304
305 //Update node position, velocity, acceleration based on net force.
306 void advanceNodes(const double dt,int nnodes,const vector2d *coord,
307                   vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,bool dampen)
308 {
309   const double nodeMass=1.0e-6; //Node mass, kilograms (HACK: hardcoded)
310   const double xm=1.0/nodeMass; //Inverse of node mass
311   const vector2d z(0,0);
312
313   const double shearForce=1.0e-6/(dt*dt*xm);
314
315   bool someNaNs=false;
316   int i;
317   for (i=0;i<nnodes;i++) {
318     vector2d R_n=R_net[i];
319 #if NANCHECK
320     if (((R_n.x-R_n.x)!=0)) {
321             CkPrintf("%d (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
322             someNaNs=true;
323     }
324 #endif
325     R_net[i]=z;
326 //Apply boundary conditions (HACK: hardcoded!)
327     if (1) {
328        if (coord[i].x<0.00001)
329                continue; //Left edge will NOT move
330        if (coord[i].y>0.02-0.00001)
331                R_n.x+=shearForce; //Top edge pushed hard to right
332     }
333 //Update displacement and velocity
334     vector2d aNew=R_n*xm;
335     v[i]+=(dt*0.5)*(aNew+a[i]);
336     d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
337     a[i]=aNew;   
338     //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
339   }
340   if (dampen)
341     for (i=0;i<nnodes;i++)
342           v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
343   if (someNaNs) {
344           CkPrintf("Nodes all NaN!\n");
345           CkAbort("Node forces NaN!");
346   }
347 }