Revert "had forgotten to add elements to pup"
[charm.git] / examples / ParFUM / 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
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   // Add a single ghost layer
131   //FEM_Add_ghost_layer(1, 1);  
132   //const static int tri2node[]={0,1,2};
133   //FEM_Add_ghost_elem(1, 3, tri2node);
134
135   CkPrintf("Finished with init (Reading took %.f s)\n",CmiWallTimer()-startTime);
136
137 }
138
139
140 // A driver() function 
141 // driver() is required in all FEM programs
142 extern "C" void
143 driver(void)
144 {
145   int nnodes,nelems,ignored;
146   int i, myId=FEM_My_partition();
147   myGlobals g;
148   FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
149   
150   
151   int mesh=FEM_Mesh_default_read(); // Tell framework we are reading data from the mesh
152   
153   // Get node data
154   nnodes=FEM_Mesh_get_length(mesh,FEM_NODE); // Get number of nodes
155   g.coord=new vector2d[nnodes];
156   // Get node positions
157   FEM_Mesh_data(mesh, FEM_NODE, FEM_DATA+0, (double*)g.coord, 0, nnodes, FEM_DOUBLE, 2);  
158
159
160   // Get element data
161   nelems=FEM_Mesh_get_length(mesh,FEM_ELEM+0); // Get number of elements
162   g.nnodes=nnodes; g.nelems=nelems;
163   g.conn=new connRec[nelems];
164   g.S11=new double[nelems];
165   g.S22=new double[nelems];
166   g.S12=new double[nelems];
167   // Get connectivity for elements
168   FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, nelems, FEM_INDEX_0, 3);  
169
170
171   //Initialize associated data
172   g.R_net=new vector2d[nnodes]; //Net force
173   g.d=new vector2d[nnodes];//Node displacement
174   g.v=new vector2d[nnodes];//Node velocity
175   g.a=new vector2d[nnodes];//Node acceleration
176   for (i=0;i<nnodes;i++)
177     g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
178
179 //Apply a small initial perturbation to positions
180   for (i=0;i<nnodes;i++) {
181           const double max=1.0e-10/15.0; //Tiny perturbation
182           g.d[i].x+=max*(i&15);
183           g.d[i].y+=max*((i+5)&15);
184   }
185
186   int fid=FEM_Create_simple_field(FEM_DOUBLE,2);
187
188   //Timeloop
189   if (myId==0)
190     CkPrintf("Entering timeloop\n");
191   int tSteps=5000;
192   double startTime, totalStart;
193   startTime=totalStart=CkWallTimer();
194   for (int t=0;t<tSteps;t++) {
195     if (1) { //Structural mechanics
196
197     //Compute forces on nodes exerted by elements
198         CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,nnodes,nelems,g.S11,g.S22,g.S12);
199
200     //Communicate net force on shared nodes
201         FEM_Update_field(fid,g.R_net);
202
203     //Advance node positions
204         advanceNodes(dt,nnodes,g.coord,g.R_net,g.a,g.v,g.d,0);
205     }
206
207     //Debugging/perf. output
208     double curTime=CkWallTimer();
209     double total=curTime-startTime;
210     startTime=curTime;
211     if (myId==0 && (t%64==0)) {
212             CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);
213             if (0) {
214               CkPrintf("    Triangle 0:\n");
215               for (int j=0;j<3;j++) {
216                     int n=g.conn[0][j];
217                     CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
218                              n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
219               }
220             }
221     }
222     /* perform migration-based load balancing */
223     if (t%1024==0)
224       FEM_Migrate();
225     
226     if (t%1024==0) { //Publish data to the net
227             NetFEM n=NetFEM_Begin(FEM_My_partition(),t,2,NetFEM_POINTAT);
228             
229             NetFEM_Nodes(n,nnodes,(double *)g.coord,"Position (m)");
230             NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
231             NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");
232             
233             NetFEM_Elements(n,nelems,3,(int *)g.conn,"Triangles");
234                 NetFEM_Scalar(n,g.S11,1,"X Stress (pure)");
235                 NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
236                 NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");
237             
238             NetFEM_End(n);
239     }
240   }
241
242   if (myId==0) {
243     double elapsed=CkWallTimer()-totalStart;
244     CkPrintf("Driver finished: average %.6f s/step\n",elapsed/tSteps);
245   }
246 }
247
248
249 // A PUP function to allow for migration and load balancing of mesh partitions.
250 // The PUP function is not needed if no migration or load balancing is desired.
251 void pup_myGlobals(pup_er p,myGlobals *g) 
252 {
253   FEM_Print("-------- called pup routine -------");
254   pup_int(p,&g->nnodes);
255   pup_int(p,&g->nelems);
256   int nnodes=g->nnodes, nelems=g->nelems;
257   if (pup_isUnpacking(p)) {
258     g->coord=new vector2d[nnodes];
259     g->conn=new connRec[nelems];
260     g->R_net=new vector2d[nnodes]; //Net force
261     g->d=new vector2d[nnodes];//Node displacement
262     g->v=new vector2d[nnodes];//Node velocity
263     g->a=new vector2d[nnodes];
264         g->S11=new double[nelems];
265         g->S22=new double[nelems];
266         g->S12=new double[nelems];
267   }
268   pup_doubles(p,(double *)g->coord,2*nnodes);
269   pup_ints(p,(int *)g->conn,3*nelems);
270   pup_doubles(p,(double *)g->R_net,2*nnodes);
271   pup_doubles(p,(double *)g->d,2*nnodes);
272   pup_doubles(p,(double *)g->v,2*nnodes);
273   pup_doubles(p,(double *)g->a,2*nnodes);
274   pup_doubles(p,(double *)g->S11,nelems);
275   pup_doubles(p,(double *)g->S22,nelems);
276   pup_doubles(p,(double *)g->S12,nelems);
277   if (pup_isDeleting(p)) {
278     delete[] g->coord;
279     delete[] g->conn;
280     delete[] g->R_net;
281     delete[] g->d;
282     delete[] g->v;
283     delete[] g->a;
284         delete[] g->S11;
285         delete[] g->S22;
286         delete[] g->S12;
287   }
288 }
289
290
291
292 //Update node position, velocity, acceleration based on net force.
293 void advanceNodes(const double dt,int nnodes,const vector2d *coord,
294                   vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,bool dampen)
295 {
296   const double nodeMass=1.0e-6; //Node mass, kilograms (HACK: hardcoded)
297   const double xm=1.0/nodeMass; //Inverse of node mass
298   const vector2d z(0,0);
299
300   const double shearForce=1.0e-6/(dt*dt*xm);
301
302   bool someNaNs=false;
303   int i;
304   for (i=0;i<nnodes;i++) {
305     vector2d R_n=R_net[i];
306 #if NANCHECK
307     if (((R_n.x-R_n.x)!=0)) {
308             CkPrintf("%d (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
309             someNaNs=true;
310     }
311 #endif
312     R_net[i]=z;
313 //Apply boundary conditions (HACK: hardcoded!)
314     if (1) {
315        if (coord[i].x<0.00001)
316                continue; //Left edge will NOT move
317        if (coord[i].y>0.02-0.00001)
318                R_n.x+=shearForce; //Top edge pushed hard to right
319     }
320 //Update displacement and velocity
321     vector2d aNew=R_n*xm;
322     v[i]+=(dt*0.5)*(aNew+a[i]);
323     d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
324     a[i]=aNew;   
325     //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
326   }
327   if (dampen)
328     for (i=0;i<nnodes;i++)
329           v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
330   if (someNaNs) {
331           CkPrintf("Nodes all NaN!\n");
332           CkAbort("Node forces NaN!");
333   }
334 }