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