A new directory for ParFUM example programs, with its own simple2D example which...
[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   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   //Timeloop
184   if (myId==0)
185     CkPrintf("Entering timeloop\n");
186   int tSteps=5000;
187   double startTime, totalStart;
188   startTime=totalStart=CkWallTimer();
189   for (int t=0;t<tSteps;t++) {
190     if (1) { //Structural mechanics
191
192     //Compute forces on nodes exerted by elements
193         CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,nnodes,nelems,g.S11,g.S22,g.S12);
194
195     //Communicate net force on shared nodes
196         FEM_Update_field(fid,g.R_net);
197
198     //Advance node positions
199         advanceNodes(dt,nnodes,g.coord,g.R_net,g.a,g.v,g.d,0);
200     }
201
202     //Debugging/perf. output
203     double curTime=CkWallTimer();
204     double total=curTime-startTime;
205     startTime=curTime;
206     if (myId==0 && (t%64==0)) {
207             CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);
208             if (0) {
209               CkPrintf("    Triangle 0:\n");
210               for (int j=0;j<3;j++) {
211                     int n=g.conn[0][j];
212                     CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
213                              n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
214               }
215             }
216     }
217     /* perform migration-based load balancing */
218     if (t%1024==0)
219       FEM_Migrate();
220     
221     if (t%1024==0) { //Publish data to the net
222             NetFEM n=NetFEM_Begin(FEM_My_partition(),t,2,NetFEM_POINTAT);
223             
224             NetFEM_Nodes(n,nnodes,(double *)g.coord,"Position (m)");
225             NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
226             NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");
227             
228             NetFEM_Elements(n,nelems,3,(int *)g.conn,"Triangles");
229                 NetFEM_Scalar(n,g.S11,1,"X Stress (pure)");
230                 NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
231                 NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");
232             
233             NetFEM_End(n);
234     }
235   }
236
237   if (myId==0) {
238     double elapsed=CkWallTimer()-totalStart;
239     CkPrintf("Driver finished: average %.6f s/step\n",elapsed/tSteps);
240   }
241 }
242
243
244 // A PUP function to allow for migration and load balancing of mesh partitions.
245 // The PUP function is not needed if no migration or load balancing is desired.
246 void pup_myGlobals(pup_er p,myGlobals *g) 
247 {
248   FEM_Print("-------- called pup routine -------");
249   pup_int(p,&g->nnodes);
250   pup_int(p,&g->nelems);
251   int nnodes=g->nnodes, nelems=g->nelems;
252   if (pup_isUnpacking(p)) {
253     g->coord=new vector2d[nnodes];
254     g->conn=new connRec[nelems];
255     g->R_net=new vector2d[nnodes]; //Net force
256     g->d=new vector2d[nnodes];//Node displacement
257     g->v=new vector2d[nnodes];//Node velocity
258     g->a=new vector2d[nnodes];
259         g->S11=new double[nelems];
260         g->S22=new double[nelems];
261         g->S12=new double[nelems];
262   }
263   pup_doubles(p,(double *)g->coord,2*nnodes);
264   pup_ints(p,(int *)g->conn,3*nelems);
265   pup_doubles(p,(double *)g->R_net,2*nnodes);
266   pup_doubles(p,(double *)g->d,2*nnodes);
267   pup_doubles(p,(double *)g->v,2*nnodes);
268   pup_doubles(p,(double *)g->a,2*nnodes);
269   pup_doubles(p,(double *)g->S11,nelems);
270   pup_doubles(p,(double *)g->S22,nelems);
271   pup_doubles(p,(double *)g->S12,nelems);
272   if (pup_isDeleting(p)) {
273     delete[] g->coord;
274     delete[] g->conn;
275     delete[] g->R_net;
276     delete[] g->d;
277     delete[] g->v;
278     delete[] g->a;
279         delete[] g->S11;
280         delete[] g->S22;
281         delete[] g->S12;
282   }
283 }
284
285
286
287 //Update node position, velocity, acceleration based on net force.
288 void advanceNodes(const double dt,int nnodes,const vector2d *coord,
289                   vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,bool dampen)
290 {
291   const double nodeMass=1.0e-6; //Node mass, kilograms (HACK: hardcoded)
292   const double xm=1.0/nodeMass; //Inverse of node mass
293   const vector2d z(0,0);
294
295   const double shearForce=1.0e-6/(dt*dt*xm);
296
297   bool someNaNs=false;
298   int i;
299   for (i=0;i<nnodes;i++) {
300     vector2d R_n=R_net[i];
301 #if NANCHECK
302     if (((R_n.x-R_n.x)!=0)) {
303             CkPrintf("%d (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
304             someNaNs=true;
305     }
306 #endif
307     R_net[i]=z;
308 //Apply boundary conditions (HACK: hardcoded!)
309     if (1) {
310        if (coord[i].x<0.00001)
311                continue; //Left edge will NOT move
312        if (coord[i].y>0.02-0.00001)
313                R_n.x+=shearForce; //Top edge pushed hard to right
314     }
315 //Update displacement and velocity
316     vector2d aNew=R_n*xm;
317     v[i]+=(dt*0.5)*(aNew+a[i]);
318     d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
319     a[i]=aNew;   
320     //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
321   }
322   if (dampen)
323     for (i=0;i<nnodes;i++)
324           v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
325   if (someNaNs) {
326           CkPrintf("Nodes all NaN!\n");
327           CkAbort("Node forces NaN!");
328   }
329 }