1c84cc247875c84992af42ef6b096a5f14f5f9c8
[charm.git] / doc / fem / manual.tex
1 \documentclass[10pt]{article}
2 \usepackage{../pplmanual}
3 \input{../pplmanual}
4
5 \makeindex
6
7 \title{\charmpp\\ Finite Element Framework\\ Manual}
8 \version{1.0}
9 \credits{
10 Initial version of \charmpp{} Finite Element Framework was developed
11 by Milind Bhandarkar with inputs from Timothy Hinrichs and Kathikeyan
12 Mahesh. The current version is almost completely rewritten by
13 Orion Lawlor.
14 }
15
16 \begin{document}
17
18 \maketitle
19
20 \section{Motivation}
21
22 The Finite Element Method (FEM) approach is used in many engineering
23 applications with irregular domains, from elastic deformation problems to
24 crack propagation to fluid flow.  \charmpp{} is a free message-passing parallel
25 runtime system for machines from clusters of workstations to tightly-coupled
26 SMPs.  The \charmpp{} FEM framework allows you to write a parallel FEM program,
27 in C or Fortran 90, that closely resembles a serial version but includes
28 a few framework calls.
29 Using the FEM framework also allows you to take advantage of all the
30 features of \charmpp, including run-time load balancing,  performance
31 monitoring and visualization, and checkpoint/restart, with no additional
32 effort.
33
34
35 \section{Introduction/Terminology}
36
37 A FEM program manipulates elements and nodes. An element is a portion of
38 the problem domain, typically in the shape of a triangle, square, or hexagon
39 in 2D; or tetrahedron or rectangular solid in 3D.  A node is a point in the
40 domain.  An element knows which nodes surround it via the connectivity
41 table, which lists adjacent nodes.
42
43 \begin{figure}[h]
44 \begin{center}
45 \includegraphics[width=4in]{simple_mesh}
46 \end{center}
47 \caption{3-element, 5 node mesh.}
48 \label{fig:simplemesh}
49 \end{figure}
50
51 \begin{table}[h]
52 \begin{center}
53 \begin{tabular}{||l||r|r|r||}\hline
54 Element & \multicolumn{3}{c||}{Adjacent Nodes} \\\hline
55 e1 & n1 & n3 & n4 \\
56 e2 & n1 & n2 & n4 \\
57 e3 & n2 & n4 & n5 \\
58 \hline
59 \end{tabular}
60 \end{center}
61 \caption{Connectivity table for mesh in figure~\ref{fig:simplemesh}.}
62 \label{table:simplemesh}
63 \end{table}
64
65 A typical FEM program performs some element-by-element calculations which
66 update adjacent node values; then some node-by-node calculations.  For
67 example, a material dynamics program has the structure:
68
69 \begin{alltt}
70      time loop
71           element loop-- Element deformation applies forces to
72           surrounding nodes
73           node loop-- Forces and boundary conditions change node
74           positions
75      end time loop
76 \end{alltt}
77
78 We can parallelize such FEM programs by partitioning the serial mesh
79 elements into several chunks  (at least one chunk per processor; perhaps
80 even more).  During partitioning, we give nodes and elements new,
81 chunk-local numbers.  Below, we partition the mesh above into two chunks, A
82 and B.
83
84 \begin{figure}[h]
85 \begin{center}
86 \includegraphics[width=4in]{partitioned_mesh}
87 \end{center}
88 \caption{Partitioned mesh.}
89 \label{fig:partitionedmesh}
90 \end{figure}
91
92 \begin{table}[h]
93 \begin{center}
94 \begin{tabular}{||l||r|r|r||}\hline
95 Element & \multicolumn{3}{c||}{Adjacent Nodes} \\\hline
96 e1 & n1 & n3 & n4 \\
97 e2 & n1 & n2 & n4 \\
98 \hline
99 \end{tabular}
100 \end{center}
101 \caption{Connectivity table for chunk A in figure~\ref{fig:partitionedmesh}.}
102 \label{table:chunkA}
103 \end{table}
104
105 \begin{table}[h]
106 \begin{center}
107 \begin{tabular}{||l||r|r|r||}\hline
108 Element & \multicolumn{3}{c||}{Adjacent Nodes}\\\hline
109 e1 & n1 & n2 & n3 \\
110 \hline
111 \end{tabular}
112 \end{center}
113 \caption{Connectivity table for chunk B in figure~\ref{fig:partitionedmesh}.}
114 \label{table:chunkB}
115 \end{table}
116
117 Note that chunk A's node n2 and B's node n1 were actually the same node in
118 the original mesh-- partitioning split this single node into two shared
119 copies (one on each chunk).  However, since adding forces is associative, we
120 can handle shared nodes by computing the forces normally (ignoring the
121 existence of the other chunk), then adding both chunks' net force for the
122 shared node together.  This ``node update'' will give us the same resulting
123 force on each shared node as we would get without partitioning, thus the
124 same positions, thus the same final result.  Hence, each chunk's time loop
125 has the structure:
126
127 \begin{alltt}
128      chunk time loop
129           element loop-- Element deformation applies forces to
130           surrounding nodes
131           <update forces on shared nodes>
132           node loop-- Forces and boundary conditions change node
133           positions
134      end time loop
135 \end{alltt}
136
137 This is exactly the form of the time loop for a \charmpp{} FEM framework
138 program.  The framework will accept a serial mesh, partition it, distribute
139 the chunks to each processor, allow the user to run their time loop, and
140 handle the node-updates.
141
142
143 \section{Structure of a FEM Framework Program}
144
145 A FEM framework program consists of three subroutines: \kw{init}, \kw{driver},
146 and \kw{finalize}.  \kw{init} and \kw{finalize} are called by the FEM framework
147 only on the first processor -- these routines typically do specialized I/O,
148 startup and shutdown tasks.  \kw{driver} is called for every chunk on every
149 processor, and does the main work of the program.
150
151 \begin{alltt}
152      subroutine init
153           read the serial mesh and configuration data
154      end subroutine
155
156      subroutine driver
157           get local mesh chunk
158           time loop
159                FEM computations
160                update shared node fields
161                more FEM computations
162           end time loop
163      end subroutine
164
165      subroutine mesh_updated
166           write intermediate results; modify serial mesh
167      end subroutine
168      subroutine finalize
169            write results
170      end subroutine
171 \end{alltt}
172
173 \section{Compilation and Execution}
174
175 A FEM framework program is a \charmpp\ program, so you must begin by
176 downloading the latest source version of \charmpp\ from
177 {\tt http://charm.cs.uiuc.edu/}.  Build the source with 
178 {\tt ./SUPER\_INSTALL FEM version} or {\tt cd} into the build directory, 
179 {\tt version/tmp}, and type {\tt make FEM}.
180 To compile a FEM program, pass the {\tt -language fem} (for C) or 
181 {\tt -language femf} (for Fortran) option to {\tt charmc}.
182
183 In a charm installation, see charm/version/pgms/charm++/fem/
184 for several example and test programs.
185
186
187 \section{FEM Framework API Reference}
188
189 \subsection{Utility}
190
191 \function{int FEM\_Num\_Partitions();}
192 \function{function integer :: FEM\_Num\_Partitions()}
193
194      Return the number of mesh chunks in the current computation.  Can
195      only be called from the driver routine.
196
197 \function{int FEM\_My\_Partition();}
198 \function{function integer :: FEM\_My\_Partition()}
199
200      Return the number of the current chunk, from 0 to
201      \kw{num\_partitions}-1.  Can only be called from the driver routine.
202
203 \function{double FEM\_Timer();}
204 \function{function double precision :: FEM\_Timer()}
205
206      Return the current wall clock time, in seconds.  Resolution is
207      machine-dependent, but is at worst 10ms.
208
209 \function{void FEM\_Print\_Partition();}
210 \function{subroutine FEM\_Print\_Partition()}
211
212      Print a debugging representation of the current chunk's mesh.
213      Prints the entire connectivity array, and data associated with
214      each local node and element.
215
216 \function{void FEM\_Print(const char *str);}
217 \function{subroutine FEM\_Print(str)}
218 \args{  character*, intent(in) :: str}
219
220      Print the given string.  Works on all machines; unlike \kw{printf} or
221      \kw{print *}, which may not work on all parallel machines.
222
223 \subsection{Mesh}
224
225 These routines describe and retreive the finite element mesh for this
226 computation.  A mesh, from the framework's perspective, is a list of
227 elements, nodes, uninterpreted data associated with each, and the
228 connectivity table.  Elements and nodes have both a global number (number in
229 the serial mesh) as well as a chunk-local number (number in the partitioned
230 mesh chunk).  The FEM framework currently uses the free Metis package for
231 partitioning.
232
233 A simple program would set the serial mesh in init, get the partitioned
234 chunk in driver, and work on that chunk.  A more complex program would set
235 the initial mesh in init; get, work on, update and repartition the mesh
236 several times in driver; and perform post-processing on the reassembled,
237 modified mesh in finalize.
238
239 From \kw{init()}, the \kw{FEM\_Set\_} routines describe the serial mesh, which
240 will be partitioned into chunks.
241
242 From \kw{driver()}, the \kw{FEM\_Get\_} routines ask for the current mesh
243 chunk; the \kw{FEM\_Set\_} routines describe a new partitioned mesh chunk.  The
244 new chunk need not have the same number of elements or nodes as the old chunk;
245 but any new added nodes are assumed private (not shared).
246 \kw{FEM\_Update\_Mesh} will reassemble a serial version of the mesh from the
247 new pieces and optionally repartition the serial mesh.
248
249 From \kw{mesh\_updated()}, the \kw{FEM\_Get} and \kw{FEM\_Set} routines
250 manipulate the serial mesh.   \kw{Mesh\_updated} is only executed if you call
251 \kw{FEM\_Update\_Mesh} during driver-- otherwise, \kw{mesh\_updated} can be an
252 empty procedure.  The parameter \kw{callMeshUpdated} is passed down to
253 \kw{mesh\_update}.
254
255 From \kw{finalize(),} the \kw{FEM\_Get\_} routines ask for the reassembled
256 serial mesh-- for this, you must have previously called \kw{FEM\_Update\_Mesh}
257 during driver.
258
259 \function{void FEM\_Set\_Mesh(int nElem, int nNodes, int nodePerEl,const int* conn);}
260
261      This is a convenience routine equivalent to:
262 \begin{alltt}
263           FEM\_Set\_Node(nNodes,0);
264           FEM\_Set\_Elem(0,nElem,0,nodePerEl);
265           FEM\_Set\_Elem\_Conn(0,conn);
266 \end{alltt}
267
268 \function{subroutine FEM\_Set\_Mesh(nElem,nNodes,nodePerEl,conn)}
269     \args{integer, intent(in) :: nElem, nNodes, nodePerEl}
270     \args{integer, intent(in), dimention(nElem,nodePerEl) :: conn;}
271
272      This is a convenience routine equivalent to:
273 \begin{alltt}
274           CALL FEM\_Set\_Node(nNodes,0)
275           CALL FEM\_Set\_Elem(1,nElem,0,nodePerEl)
276           CALL FEM\_Set\_Elem\_Conn\_c(1,conn)
277 \end{alltt}
278
279 \function{void FEM\_Set\_Elem(int elType,int  nEl,int  doublePerEl,int  nodePerEl);}
280 \function{void FEM\_Get\_Elem(int elType,int *nEl,int *doublePerEl,int *nodePerEl);}
281 \function{subroutine FEM\_Set\_Elem(elType,nEl,doublePerEl,nodePerEl)}
282   \args{integer, intent(in)  :: elType,nEl,doublePerEl,nodePerEl}
283 \function{subroutine FEM\_Get\_Elem(elType,nEl,doublePerEl,nodePerEl)}
284   \args{integer, intent(in)  :: elType}
285   \args{integer, intent(out) :: nEl,doublePerEl,nodePerEl}
286
287      Describe/retreive the number and type of elements.  \kw{ElType} is the
288 number of element types registered so far (the first element type is 1, then 2,
289 etc.).  \kw{nEl} is the number of elements being registered.  \kw{doublesPerEl}
290 and \kw{nodePerEl} are the number of doubles of user data, and nodes
291 (respectively) associated with each element.
292
293      \kw{doublePerEl} or \kw{nodePerEl} may be zero, indicating that no user
294 data or connectivity data (respectively) is associated with the element.
295
296 \function{void FEM\_Set\_Elem\_Conn(int elType,const int *conn);}
297 \function{void FEM\_Get\_Elem\_Conn(int elType,int *conn);}
298 \function{subroutine FEM\_Set\_Elem\_Conn\_r(elType,conn)}
299   \args{integer, intent(in)  :: elType}
300   \args{integer, intent(in),  dimention(nodePerEl,nEl) :: conn}
301 \function{subroutine FEM\_Get\_Elem\_Conn\_r(elType,conn)}
302   \args{integer, intent(in)  :: elType}
303   \args{integer, intent(out), dimention(nodePerEl,nEl) :: conn}
304 \function{subroutine FEM\_Set\_Elem\_Conn\_c(elType,conn)}
305   \args{integer, intent(in)  :: elType}
306   \args{integer, intent(in),  dimention(nEl,nodePerEl) :: conn}
307 \function{subroutine FEM\_Get\_Elem\_Conn\_c(elType,conn)}
308   \args{integer, intent(in)  :: elType}
309   \args{integer, intent(out), dimention(nEl,nodePerEl) :: conn}
310
311      Describe/retreive the element connectivity array for this element
312      type.  The connectivity array is indexed by the element number,
313      and gives the indices of the nodes surrounding the element.  It is
314      hence \kw{nodePerEl*nEl} integers long.
315
316      The C version array indices are zero-based, and must be stored in
317      row-major order (a given element's surrounding nodes are stored
318      contiguously in the conn array).  The Fortran version indices are
319      one-based, and are available in row-major (named \_r) and
320      column-major (named \_c) versions.  We recommend row-major storage
321      because it results in better cache utilization (because the nodes
322      around an element are stored contiguously).
323
324 \function{void FEM\_Set\_Node(int  nNode,int  doublePerNode);}
325 \function{void FEM\_Get\_Node(int *nNode,int *doublePerNode);}
326 \function{subroutine FEM\_Set\_Node(nNode,doublePerNode)}
327   \args{integer, intent(in)  :: nNode,doublePerNode}
328 \function{subroutine FEM\_Get\_Node(nNode,doublePerNode)}
329   \args{integer, intent(out) :: nNode,doublePerNode}
330
331      Describe/retreive the number of nodes and doubles of user data
332      associated with each node.  There is only one type of node, so no
333      \kw{nodeType} identifier is needed.
334
335      \kw{doublePerNode} may be zero, indicating that no user data is
336      associated with each node.
337
338 \function{void FEM\_Set\_Node\_Data(const double *data);}
339 \function{void FEM\_Get\_Node\_Data(double *data);}
340 \function{void FEM\_Set\_Elem\_Data(int elType,const double *data);}
341 \function{void FEM\_Get\_Elem\_Data(int elType,double *data);}
342 \function{subroutine FEM\_Set\_Node\_Data\_r(data)}
343   \args{REAL*8, intent(in),  dimention(doublePerNode,nNode)  :: data}
344 \function{subroutine FEM\_Get\_Node\_Data\_r(data)}
345   \args{REAL*8, intent(out), dimention(doublePerNode,nNode)  :: data}
346 \function{subroutine FEM\_Set\_Elem\_Data\_r(data)}
347   \args{REAL*8, intent(in),  dimention(doublePerElem,nElem)  :: data}
348 \function{subroutine FEM\_Get\_Elem\_Data\_r(data)}
349   \args{REAL*8, intent(out), dimention(doublePerElem,nElem)  :: data}
350 \function{subroutine FEM\_Set\_Node\_Data\_c(data)}
351   \args{REAL*8, intent(in),  dimention(nNode,doublePerNode)  :: data}
352 \function{subroutine FEM\_Get\_Node\_Data\_c(data)}
353   \args{REAL*8, intent(out), dimention(nNode,doublePerNode)  :: data}
354 \function{subroutine FEM\_Set\_Elem\_Data\_c(data)}
355   \args{REAL*8, intent(in),  dimention(nElem,doublePerElem)  :: data}
356 \function{subroutine FEM\_Get\_Elem\_Data\_c(data)}
357   \args{REAL*8, intent(out), dimention(nElem,doublePerElem)  :: data}
358
359      Describe/retrieve the optional, uninterpreted user data associated with
360 each node and element.  This user data is partitioned and reassembled along
361 with the connectivity matrix, and may include initial conditions, boundary
362 values, or any other data needed or produced by the program.   The Fortran
363 arrays can be row- or column- major (see \kw{FEM\_Set\_Elem\_Conn} for
364 details).  The row-major form is preferred.
365
366 \function{void FEM\_Update\_Mesh(int callMeshUpdated,int doRepartition);}
367 \function{subroutine FEM\_Update\_Mesh(callMeshUpdated,doRepartition)}
368     \args{integer, intent(in) :: callMeshUpdated,doRepartition}
369
370      Reassemble the mesh chunks from each partition into a single serial mesh.
371 Can only be called from driver; and must be called by the driver routine for
372 every chunk.  The call is blocking only if \kw{doRepartition} is true;
373 otherwise the call will immediately return.  \kw{FEM\_Get} calls from
374 \kw{driver()} will only return the new mesh after a \kw{FEM\_Update\_Mesh} call
375 where \kw{doRepartition} is true; otherwise \kw{FEM\_Get} returns the old mesh.
376
377      If \kw{callMeshUpdated} is not zero, \kw{mesh\_updated(callMeshUpdated)}
378 will be called on the first processor after the mesh is reassembled.
379
380      If \kw{doRepartition} is 1, the reassembled serial mesh will be
381 repartitioned back into chunks and redistributed.  Otherwise,
382 \kw{FEM\_Update\_Mesh} will return immediately.
383
384      If both \kw{doRepartition} and \kw{callMeshUpdated} are nonzero, the
385 serial mesh will be reassembled on the first processor, \kw{mesh\_updated} will
386 be called, and the resulting mesh (which \kw{mesh\_updated} may change) is
387 repartitioned into chunks and redistributed.
388
389      If both \kw{doRepartition} and \kw{callMeshUpdated} are 0, the serial mesh
390 will be reassembled on the first processor for use in the \kw{finalize()}
391 routine.  Note that to use the serial mesh in \kw{finalize(),} you must call
392 this routine from \kw{driver()}.
393
394      \kw{FEM\_Update\_Mesh} reassembles the serial mesh with an attempt to
395      preserve the element and node global numbering.  If the new mesh
396      has the same number and type of elements and nodes, the global
397      numbers (and hence serial mesh) will be unchanged.  If new
398      elements or nodes are added at each chunk, they will be assigned
399      new unique global numbers.  If elements or nodes are removed,
400      their global numbers are not re-used-- you can detect the
401      resulting holes in the serial mesh since the user data associated
402      with the deleted elements will be all zero.
403
404      It may be easier to perform major mesh modifications from
405      \kw{mesh\_updated}, since the entire serial mesh is available there.
406
407 \subsection{Node Fields}
408
409 The FEM framework handles the updating of the values of shared nodes-- that
410 is, it combines shared nodes' values across all processors.  The basic
411 mechanism to do this update is the field-- numeric data items associated
412 with each node. We make no assumptions about the meaning of the node data,
413 and allow various data types and non-communicated data associated with each
414 node.  To do this, the framework must be able to find the data items
415 associated with each node in memory.
416
417 Each field represents a (set of) node data items stored in a contiguous array,
418 indexed by node number.  You create a field once, with \kw{FEM\_Create\_Field},
419 then pass the resulting field ID to \kw{FEM\_Update\_Field} (which does the
420 shared node communication) and/or \kw{FEM\_Reduce\_Field} (which applies a
421 reduction over node values).
422
423 \function{int FEM\_Create\_Field(int base\_type,int vec\_len,int offset,int dist);}
424 \function{function integer :: FEM\_Create\_Field(base\_type, vec\_len, offset, dist)}
425   \args{integer, intent(in)  :: base\_type, vec\_len, offset, dist}
426
427      Creates and returns a FEM field ID, which can be passed to
428 \kw{FEM\_Update\_Field} and \kw{FEM\_Reduce\_Field.}  Can only be called from
429 \kw{driver().}  A field is a range of values associated with each local node--
430 the FEM framework uses the information you pass to find the values associated
431 with shared nodes (for \kw{FEM\_Update\_Field}) and primary nodes (for
432 \kw{FEM\_Reduce\_Field}).
433
434      \kw{base\_type} describes the kind of data item associated with each
435      node, one of:
436
437      \begin{itemize}
438         \item \kw{FEM\_BYTE}-- unsigned char, INTEGER*1, or CHARACTER*1
439         \item \kw{FEM\_INT}-- int or INTEGER*4
440         \item \kw{FEM\_REAL}-- float or REAL*4
441         \item \kw{FEM\_DOUBLE}-- double, DOUBLE PRECISION, or REAL*8
442      \end{itemize}
443
444      \kw{vec\_len} describes the number of data items associated with each
445      node, an integer at least 1.
446
447      \kw{offset} is the byte offset from the start of the nodes array to the
448      data items, a non-negative integer.
449
450      \kw{dist} is the byte offset from the first node's data item to the
451      second, a positive integer.
452
453 \begin{figure}[h]
454 \begin{center}
455 \includegraphics[width=4in]{create_field}
456 \end{center}
457 \caption{Creating a Node Field.}
458 \label{fig:createfield}
459 \end{figure}
460
461      For example, if you store a 3D force for each node \kw{n}, in an array
462 indexed by 3*\kw{n}, \kw{vec\_len} is 3, \kw{offset} is 0, and \kw{dist} is
463 3*8=24.  You can register the node force for update with:
464
465 \begin{alltt}
466           /* C */
467           double *nodeForce;
468           ... allocate nodeForce as 3*n_nodes...
469           int fid=FEM_Create_Field(FEM_DOUBLE,3,0,24);
470  
471           ! - Fortran90
472           REAL*8 ALLOCATABLE, DIMENTION(:) :: nodeForce
473           INTEGER :: fid
474           ... allocate nodeForce as 3*n_nodes...
475           fid=FEM_Create_Field(FEM_DOUBLE,3,0,24)
476 \end{alltt}
477
478      If the 3D force is a member \kw{fXYZ} of a structure (in C) or named type
479 (in Fortran 90) \kw{node\_type}, in an array called \kw{nodes}, you can
480 register this node force for update with:
481
482 \begin{alltt}
483           /* C */
484           node_type *nodes;
485           ...allocate nodes array as n_nodes...
486           int fid=FEM_Create_Field(FEM_DOUBLE,3,
487               (int)((char *)\&nodes[0].fXYZ-(char *)nodes),
488               (int)((char *)\&nodes[1]-(char *)\&nodes[0]) );
489  
490           ! - Fortran90
491           TYPE(node_type), ALLOCATABLE, DIMENTION(:) :: nodes
492           INTEGER :: fid
493           ...allocate nodes array as n_nodes...
494           fid=FEM_Create_Field(FEM_DOUBLE,3,
495               offsetof(nodes(1), nodes(1)%fXYZ),
496               offsetof(nodes(1), nodes(2)) )
497 \end{alltt}
498
499      This example uses the Fortran-only helper routine \kw{offsetof}, which
500      returns the offset in bytes of memory between its two given
501      variables.  The C version uses pointer arithmetic to achieve the
502      same result.
503
504 \function{void FEM\_Update\_Field(int fid,void *nodes);}
505 \function{subroutine FEM\_Update\_Field(fid,nodes)}
506   \args{integer, intent(in)  :: fid}
507   \args{varies, intent(inout) :: nodes}
508
509      Combine a field of all shared nodes with the other chunks.  Sums
510      the value of the given field across all chunks that share each
511      node.  For the example above, once each chunk has computed the net
512      force on each local node, this routine will sum the net force
513      across all shared nodes.
514
515      \kw{FEM\_Update\_Field} can only be called from driver, and to be useful,
516      must be called from every chunk's driver routine.
517
518      After this routine returns, the given field of each shared node
519      will be the same across all processors that share the node.
520
521 \function{void FEM\_Read\_Field(int fid,void *nodes,char *fName);}
522 \function{subroutine FEM\_Read\_Field(fid,nodes,fName)}
523   \args{integer, intent(in)  :: fid}
524   \args{varies, intent(out) :: nodes}
525   \args{character*, intent(in) :: fName}
526
527      Read a field out of the given serial input file.  The serial input
528      file is line-oriented ASCII-- each line begins with the global
529      node number (which must match the line order in the file),
530      followed by the data to be read into the node field.  The
531      remainder of each line is unread.  If called from Fortran, the
532      first line must be numbered 1; if called from C, the first line
533      must be numbered zero.  All fields are separated by white space
534      (any number of tabs or spaces).
535
536      For example, if we have called \kw{Create\_Field} to describe 3 doubles,
537      the input file could begin with
538
539 \begin{alltt}
540           1    0.2    0.7    -0.3      First node
541           2    0.4    1.12   -17.26    another node
542           ...
543 \end{alltt}
544
545      \kw{FEM\_Read\_Field} must be called from driver at any time, independent
546      of other chunks.
547
548 \function{void FEM\_Reduce\_Field(int fid,const void *nodes,void *out,int op);}
549 \function{subroutine FEM\_Reduce\_Field(fid,nodes,outVal,op)}
550   \args{integer, intent(in)  :: fid,op}
551   \args{varies, intent(in) :: nodes}
552   \args{varies, intent(out) :: outVal}
553
554      Combine a field from each node, according to op, across all chunks.
555 Shared nodes are not double-counted-- only once copy will contribute to the
556 reduction.  After \kw{Reduce\_Field} returns, all chunks will have identical
557 values in \kw{outVal,} which must be \kw{vec\_len} copies of \kw{base\_type}.
558
559      May only be called from driver, and to complete, must be called
560      from every chunk's driver routine.
561
562      \kw{op} must be one of:
563
564 \begin{itemize}
565         \item \kw{FEM\_SUM}-- each element of \kw{outVal} will be the sum 
566 of the corresponding fields of all nodes
567         \item \kw{FEM\_MIN}-- each element of \kw{outVal} will be the 
568 smallest value among the corresponding field of all nodes
569         \item \kw{FEM\_MAX}-- each element of \kw{outVal} will be the largest 
570 value among the corresponding field of all nodes
571 \end{itemize}
572
573 \function{void FEM\_Reduce(int fid,const void *inVal,void *outVal,int op);}
574 \function{subroutine FEM\_Reduce(fid,inVal,outVal,op)}
575   \args{integer, intent(in)  :: fid,op}
576   \args{varies, intent(in) :: inVal}
577   \args{varies, intent(out) :: outVal}
578
579      Combine a field from each chunk, acoording to \kw{op}, across all chunks.
580 \kw{Fid} is only used for the \kw{base\_type} and \kw{vec\_len}-- offset and
581 \kw{dist} are not used.  After this call returns, all chunks will have
582 identical values in \kw{outVal}.  Op has the same values and meaning as
583 \kw{FEM\_Reduce\_Field}.
584
585      May only be called from driver, and to complete, must be called
586      from every chunk's driver routine.
587
588 \subsection{Migration}
589
590 The \charmpp\ runtime framework includes an automated, run-time load balancer,
591 which will automatically monitor the performance of your parallel program.
592 If needed, the load balancer can ``migrate'' mesh chunks from heavily-loaded
593 processors to more lightly-loaded processors, improving the load balance and
594 speeding up the program.  For this to be useful, pass the \kw{+vpN} argument
595 with a larger number of chunks \kw{N} than processors-- a few thousand
596 finite elements per chunk acheives a good trade-off between overhead and
597 ability to load balance.  Because this is somewhat involved, you
598 may refrain from calling \kw{FEM\_Migrate} and migration will never take place.
599
600 The runtime system can automatically move your thread stack to the new
601 processor, but you must write a PUP function to move any global or
602 heap-allocated data to the new processor (global data is declared at file scope
603 or \kw{static} in C and \kw{COMMON} in Fortran77; heap allocated data comes
604 from C \kw{malloc}, C++ \kw{new}, or Fortran90 \kw{ALLOCATE}).  A PUP
605 (Pack/UnPack) function performs both packing (converting heap data into a
606 message) and unpacking (converting a message back into heap data).  All your
607 global and heap data must be collected into a single block (\kw{struct} in C;
608 user-defined \kw{TYPE} in Fortran) so the PUP function can access it all.
609
610 Your PUP function will be passed a pointer to your heap data block and a
611 special handle called a ``pupper'', which contains the network message to be
612 sent.  Your PUP function returns a pointer to your heap data block.  In a PUP
613 function, you pass all your heap data to routines named \kw{pup\_type}, where
614 type is either a basic type (such as int, char, float, or double) or an array
615 type (as before, but with a ``s'' suffix).  Depending on the direction of
616 packing, the pupper will either read from or write to the values you pass--
617 normally, you shouldn't even know which.  The only time you need to know the
618 direction is when you are leaving a processor or just arriving.
619 Correspondingly, the pupper passed to you may be deleting (indicating that you
620 are leaving the processor, and should delete your heap storage after packing),
621 unpacking (indicating you've just arrived on a processor, and should allocate
622 your heap storage before unpacking), or neither (indicating the system is
623 merely sizing a buffer, or checkpointing your values).
624
625 PUP functions are much easier to write than explain-- a simple C heap block
626 and the corresponding PUP function is:
627
628 \begin{alltt}
629      typedef struct {
630        int n1;/*Length of first array below*/
631        int n2;/*Length of second array below*/
632        double *arr1; /*Some doubles, allocated on the heap*/
633        int *arr2; /*Some ints, allocated on the heap*/
634      } my_block;
635  
636      my_block *pup_my_block(pup_er p,my_block *m)
637      {
638        if (pup_isUnpacking(p)) m=malloc(sizeof(my_block));
639        pup_int(p,\&m->n1);
640        pup_int(p,\&m->n2);
641        if (pup_isUnpacking(p)) {
642          m->arr1=malloc(m->n1*sizeof(double));
643          m->arr2=malloc(m->n2*sizeof(int));
644        }
645        pup_doubles(p,m->arr1,m->n1);
646        pup_ints(p,m->arr2,m->n2);
647        if (pup_isDeleting(p)) {
648          free(m->arr1);
649          free(m->arr2);
650          free(m);
651        }
652        return m;
653      }
654 \end{alltt}
655
656 This single PUP function can be used to copy the \kw{my\_block} data into a
657 message buffer and free the old heap storage (deleting pupper); allocate
658 storage on the new processor and copy the message data back (unpacking pupper);
659 or save the heap data for debugging or checkpointing.
660
661 A Fortran block TYPE and corresponding PUP routine is as follows:
662
663 \begin{alltt}
664      MODULE my_block_mod
665        TYPE my_block
666          INTEGER :: n1,n2x,n2y
667          REAL*8, POINTER, DIMENSION(:) :: arr1
668          INTEGER, POINTER, DIMENSION(:,:) :: arr2
669        END TYPE
670      END MODULE
671  
672      SUBROUTINE pup_my_block(p,m)
673        IMPLICIT NONE
674        USE my_block_mod
675        USE pupmod
676        INTEGER :: p
677        TYPE(my_block) :: m
678        call pup_int(p,m%n1)
679        call pup_int(p,m%n2x)
680        call pup_int(p,m%n2y)
681        IF (pup_isUnpacking(p)) THEN
682          ALLOCATE(m%arr1(m%n1))
683          ALLOCATE(m%arr2(m%n2x,m%n2y))
684        END IF
685        call pup_doubles(p,m%arr1,m%n1)
686        call pup_ints(p,m%arr2,m%n2x*m%n2y)
687        IF (pup_isDeleting(p)) THEN
688          DEALLOCATE(m%arr1)
689          DEALLOCATE(m%arr2)
690        END IF
691      END SUBROUTINE
692 \end{alltt}
693
694 \function{int FEM\_Register(void *block, FEM\_PupFn pup\_ud)}
695 \function{function integer :: FEM\_Register(block,pup\_ud)}
696     \args{TYPE(varies), POINTER :: block}
697     \args{SUBROUTINE :: pup\_ud}
698
699      Associates the given data block and PUP function.  Returns a block
700      ID, which can be passed to \kw{FEM\_Get\_Userdata} later.  Can only be
701      called from driver.  For the declarations above, you call
702      \kw{FEM\_Register} as:
703
704 \begin{alltt}
705           /*C/C++ driver() function*/
706           my_block *m=malloc(sizeof(my_block));
707           int myId=FEM_Register(m,(FEM_PupFn)pup_my_block);
708  
709           !- Fortran driver subroutine
710           use my_block_mod
711           interface
712             subroutine pup_my_block(p,m)
713               use my_block_mod
714               INTEGER :: p
715               TYPE(my_block) :: m
716             end subroutine
717           end interface
718           TYPE(my_block) :: m
719           INTEGER :: myId
720           myId=FEM_Register(m,pup_my_block)
721 \end{alltt}
722
723      Note that Fortran blocks must be allocated on the stack in driver;
724      while C/C++ blocks may be allocated on the heap.
725
726 \function{void FEM\_Migrate()}
727 \function{subroutine FEM\_Migrate()}
728
729      Informs the load balancing system that you are ready to be
730      migrated, if needed.  If the system decides to migrate you, the
731      PUP function passed to \kw{FEM\_Register} will be called with a sizing
732      pupper, then a packing, deleting pupper.  Your stack (and pupped
733      data) will then be sent to the destination machine, where your PUP
734      function will be called with an unpacking pupper.  \kw{FEM\_Migrate}
735      will then return, whereupon you should call \kw{FEM\_Get\_Userdata} to
736      get your unpacked data block.  Can only be called from driver.
737
738 \function{void *FEM\_Get\_Userdata(int n)}
739
740      Return your unpacked userdata after migration-- that is, the
741      return value of the unpacking call to your PUP function.  Takes
742      the userdata ID returned by \kw{FEM\_Register}.  Can be called from
743      driver at any time.
744
745      Since Fortran blocks are always allocated on the stack, the system
746      migrates them to the same location on the new processor, so no
747      \kw{Get\_Userdata} call is needed from Fortran.
748 \input{index}
749 \end{document}