Added this new manual for IFEM solvers.
[charm.git] / doc / ifem / manual.tex
1 \documentclass[10pt]{article}
2 \usepackage{../pplmanual}
3 \input{../pplmanual}
4
5 \makeindex
6
7 \title{\charmpp\\ Iterative Finite Element Matrix (IFEM) Library\\ Manual}
8 \version{1.2}
9 \credits{
10 Initial version of \charmpp{} Finite Element Framework was developed
11 by Orion Lawlor in the spring of 2003.
12 }
13
14 \begin{document}
15
16 \maketitle
17
18 \section{Introduction}
19
20 This manual presents the Iterative Finite Element Matrix (IFEM) library,
21 a library for easily solving matrix problems derived from 
22 finite-element formulations.  The library is designed to be matrix-free,
23 in that the only matrix operation required is matrix-vector product,
24 and hence the entire matrix need never be assembled
25
26 IFEM is built on the mesh and communication capabilities of the Charm++ 
27 FEM Framework, so for details on the basic runtime, problem setup, and 
28 partitioning see the FEM Framework manual.
29
30
31 \subsection{Terminology}
32
33 A FEM program manipulates elements and nodes. An \term{element} is a portion of
34 the problem domain, also known as a cell, and is typically some simple shape 
35 like a triangle, square, or hexagon in 2D; or tetrahedron or rectangular solid in 3D.  
36 A \term{node} is a point in the domain, and is often the vertex of several elements.  
37 Together, the elements and nodes form a \term{mesh}, which is the 
38 central data structure in the FEM framework.  See the FEM manual for details.
39
40 % --------------- Solvers -----------------
41 \section{Solvers}
42 \label{sec:solver}
43 A IFEM \term{solver} is a subroutine that controls the search for the solution.
44
45 Solvers often take extra parameters, which are listed in a type called 
46 in C \kw{ILSI\_Param}, which in Fortran is an array of \kw{ILSI\_PARAM} doubles.  
47 You initialize these solver parameters using the subroutine \kw{ILSI\_Param\_new}, 
48 which takes the parameters as its only argument.  The input and output
49 parameters in an \kw{ILSI\_Param} are listed in Table~\ref{table:solverIn}
50 and Table~\ref{table:solverOut}.
51
52 \begin{table}[hh]
53 \begin{center}
54 \begin{tabular}{|l|l|l|}\hline
55 C Field Name & Fortran Field Offset & Use\\\hline
56 maxResidual & 1 & If nonzero, termination criteria: magnitude of residual. \\
57 maxIterations & 2 & If nonzero, termination criteria: number of iterations. \\
58 solverIn[8] & 3-10 & Solver-specific input parameters. \\
59 \hline
60 \end{tabular}
61 \end{center}
62 \caption{\kw{ILSI\_Param} solver input parameters.}
63 \label{table:solverIn}
64 \end{table}
65
66 \begin{table}[hh]
67 \begin{center}
68 \begin{tabular}{|l|l|l|}\hline
69 C Field Name & Fortran Field Offset & Use\\\hline
70 residual & 11 & Magnitude of residual of final solution. \\
71 iterations & 12 & Number of iterations actually taken. \\
72 solverOut[8] & 13-20 & Solver-specific output parameters. \\
73 \hline
74 \end{tabular}
75 \end{center}
76 \caption{\kw{ILSI\_Param} solver output parameters.}
77 \label{table:solverOut}
78 \end{table}
79
80
81 \subsection{Conjugate Gradient Solver}
82
83 The only solver currently written using IFEM is the conjugate gradient solver.
84 This linear solver requires the matrix to be real, symmetric and positive definite.
85
86 Each iteration of the conjugate gradient solver requires one matrix-vector product
87 and two global dot products.  For well-conditioned problems, the solver typically 
88 converges in some small multiple of the diameter of the mesh--the number of elements
89 along the largest side of the mesh.
90
91 You access the conjugate gradient solver via the subroutine name \kw{ILSI\_CG\_Solver}.
92
93
94 % -------------- Shared-node ----------------
95 \section{Solving Shared-Node Systems}
96
97 Many problems encountered in FEM analysis place the entries of the
98 known and unknown vectors at the nodes--the vertices of the domain.
99 Elements provide linear relationships between the known and unknown node values, 
100 and the entire matrix expresses the combination of all these element relations.
101
102 For example, in a structural statics problem, we know the net force at
103 each node, $f$, and seek the displacements of each node, $u$.  Elements 
104 provide the relationship, often called a stiffness matrix $K$, between 
105 a nodes' displacments and its net forces:
106
107 \[
108         f=K u
109 \]
110
111 We normally label the known vector $b$ (in the example, the forces), the unknown
112 vector $x$ (in the example, the displacments), and the matrix $A$:
113
114 \[
115         b=A x
116 \]
117
118
119 IFEM provides two routines for solving problems of this type.  The first routine, \kw{IFEM\_Solve\_shared}, solves for the entire $x$ vector based on the known values of the $b$ vector.  The second, \kw{IFEM\_Solve\_shared\_bc}, allows certain entries in the $x$ vector to be given specific values before the problem is solved, creating values for the $b$ vector.
120
121
122 \newpage
123 \subsection{IFEM\_Solve\_shared}
124
125 \prototype{IFEM\_Solve\_shared}
126 \function{void IFEM\_Solve\_shared(ILSI\_Solver s,ILSI\_Param *p,
127         int fem\_mesh, int fem\_entity,int length,int width,
128         IFEM\_Matrix\_product\_c A, void *ptr, const double *b, double *x);}
129 \function{subroutine IFEM\_Solve\_shared(s,p,
130         fem\_mesh,fem\_entity,length,width,
131         A,ptr,b,x)}
132         \args{external solver subroutine :: s}
133         \args{double precision, intent(inout) :: p(ILSI\_PARAM)}
134         \args{integer, intent(in) :: fem\_mesh, fem\_entity, length,width}
135         \args{external matrix-vector product subroutine :: A}
136         \args{TYPE(varies), pointer :: ptr }
137         \args{double precision, intent(in) :: b(width,length)}
138         \args{double precision, intent(inout) :: x(width,length)}
139
140 This routine solves the linear system $A x = b$ for the unknown vector $x$.  \uw{s} and \uw{p} give the particular linear solver to use, and are described in more detail in Section~\ref{sec:solver}.  \uw{fem\_mesh} and \uw{fem\_entity} give the FEM framework mesh (often \kw{FEM\_Mesh\_default\_read()}) and entity (often \kw{FEM\_NODE}) with which the known and unknown vectors are listed.
141
142 \uw{width} gives the number of degrees of freedom (entries in the vector) per node. For example, if there is one degree of freedom per node, width is one.  \uw{length} should always equal the number of FEM nodes.
143
144 \uw{A} is a local matrix-vector product routine you must write.  Its interface is described in Section~\ref{sec:mvp}. \uw{ptr} is a pointer passed down to \uw{A}--it is not otherwise used by the framework.  
145
146 \uw{b} is the known vector.  \uw{x}, on input, is the initial guess for the unknown vector.  On output, \uw{x} is the final value for the unknown vector.  \uw{b} and \uw{x} should both have length * width entries.  In C, DOF $i$ of node $n$ should be indexed as $x[n*$\uw{width}$+i]$.  In Fortran, these arrays should be allocated like \uw{x(width,length)}.
147
148 When this routine returns, \uw{x} is the final value for the unknown vector, and the output values of the solver parameters \uw{p} will have been written.
149
150 \begin{alltt}
151 // C++ Example
152   int mesh=FEM_Mesh_default_read();
153   int nNodes=FEM_Mesh_get_length(mesh,FEM_NODE);
154   int width=3; //A 3D problem
155   ILSI_Param solverParam;
156   struct myProblemData myData;
157   
158   double *b=new double[nNodes*width];
159   double *x=new double[nNodes*width];
160   ... prepare solution target b and guess x ...
161   
162   ILSI_Param_new(&solverParam);
163   solverParam.maxResidual=1.0e-4;
164   solverParam.maxIterations=500; 
165   
166   IFEM_Solve_shared(IFEM_CG_Solver,&solverParam,
167          mesh,FEM_NODE, nNodes,width,
168          myMatrixVectorProduct, &myData, b,x);
169   
170 ! F90 Example
171   include 'ifem\_f.h'
172   INTEGER :: mesh, nNodes,width
173   DOUBLE PRECISION, ALLOCATABLE :: b(:,:), x(:,:)
174   DOUBLE PRECISION :: solverParam(ILSI_PARAM)
175   TYPE(myProblemData) :: myData
176   
177   mesh=FEM_Mesh_default_read()
178   nNodes=FEM_Mesh_get_length(mesh,FEM_NODE)
179   width=3   ! A 3D problem
180   
181   ALLOCATE(b(width,nNodes), x(width,nNodes))
182   ... prepare solution target b and guess x ..
183   
184   ILSI_Param_new(&solverParam);
185   solverParam(1)=1.0e-4;
186   solverParam(2)=500; 
187   
188   IFEM_Solve_shared(IFEM_CG_Solver,solverParam,
189          mesh,FEM_NODE, nNodes,width,
190          myMatrixVectorProduct, myData, b,x);
191
192 \end{alltt}
193
194 \newpage
195 \subsubsection{Matrix-vector product routine}
196 \label{sec:mvp}
197
198 IFEM requires you to write a matrix-vector product routine that will evaluate $A x$ for various vectors $x$.  You may use any subroutine name, but it must take these arguments:
199
200 \prototype{IFEM\_Matrix\_product}
201 \function{void IFEM\_Matrix\_product(void *ptr,int length,int width,
202         const double *src, double *dest);}
203 \function{subroutine IFEM\_Matrix\_product(ptr,length,width,src,dest)}
204   \args{TYPE(varies), pointer  :: ptr}
205   \args{integer, intent(in)  :: length,width}
206   \args{double precision, intent(in) :: src(width,length)}
207   \args{double precision, intent(out) :: dest(width,length)}
208
209 The framework calls this user-written routine when it requires a matrix-vector product.  This routine should compute $dest = A \, src$, interpreting $src$ and $dest$ as vectors.  \uw{length} gives the number of nodes and \uw{width} gives the number of degrees of freedom per node, as above.
210
211 In writing this routine, you are responsible for choosing a representation for the matrix $A$. For many problems, there is no need to represent $A$ explicitly--instead, you simply evaluate $dest$ by looping over local elements, taking into account the values of $src$.  This example shows how to write the matrix-vector product routine for simple 1D linear elastic springs, while solving for displacement given net forces.
212
213 After calling this routine, the framework will handle combining the overlapping portions of these vectors across processors to arrive at a consistent global matrix-vector product.
214
215 \begin{alltt}
216 // C++ Example
217 typedef struct \{
218   int nElements; //Number of local elements
219   int *conn; // Nodes adjacent to each element: 2*nElements entries
220   double k; //Uniform spring constant
221 \} myProblemData;
222
223 void myMatrixVectorProduct(void *ptr,int nNodes,int dofPerNode,
224           const double *src,double *dest) 
225 \{
226   myProblemData *d=(myProblemData *)ptr;
227   int n,e;
228   // Zero out output force vector:
229   for (n=0;n<nNodes;n++) dest[n]=0;
230   // Add in forces from local elements
231   for (e=0;e<d->nElements;e++) {
232     int n1=d->conn[2*e+0]; // Left node
233     int n2=d->conn[2*e+1]; // Right node
234     double f=d->k * (src[n2]-src[n1]); //Force
235     dest[n1]+=f;
236     dest[n2]-=f;
237   }
238 \}
239
240
241 ! F90 Example
242   TYPE(myProblemData) 
243     INTEGER :: nElements
244     INTEGER, ALLOCATABLE :: conn(2,:)
245     DOUBLE PRECISION :: k
246   END TYPE
247   
248 SUBROUTINE myMatrixVectorProduct(d,nNodes,dofPerNode,src,dest)
249   include 'ifem\_f.h'
250   TYPE(myProblemData), pointer :: d
251   INTEGER :: nNodes,dofPerNode
252   DOUBLE PRECISION :: src(dofPerNode,nNodes), dest(dofPerNode,nNodes)
253   INTEGER :: e,n1,n2
254   DOUBLE PRECISION :: f
255   
256   dest(:,:)=0.0
257   do e=1,d%nElements
258     n1=d%conn(1,e)
259     n2=d%conn(2,e)
260     f=d%k * (src(1,n2)-src(1,n1))
261     dest(1,n1)=dest(1,n1)+f
262     dest(1,n2)=dest(1,n2)+f
263   end do
264 END SUBROUTINE
265 \end{alltt}
266
267
268 \newpage
269 \subsection{IFEM\_Solve\_shared\_bc}
270
271 \prototype{IFEM\_Solve\_shared\_bc}
272 \function{void IFEM\_Solve\_shared\_bc(ILSI\_Solver s,ILSI\_Param *p,
273         int fem\_mesh, int fem\_entity,int length,int width,
274         int bcCount, const int *bcDOF, const double *bcValue,
275         IFEM\_Matrix\_product\_c A, void *ptr, const double *b, double *x);}
276 \function{subroutine IFEM\_Solve\_shared\_bc(s,p,
277         fem\_mesh,fem\_entity,length,width,
278         bcCount,bcDOF,bcValue,
279         A,ptr,b,x)}
280         \args{external solver subroutine :: s}
281         \args{double precision, intent(inout) :: p(ILSI\_PARAM)}
282         \args{integer, intent(in) :: fem\_mesh, fem\_entity, length,width}
283         \args{integer, intent(in) :: bcCount}
284         \args{integer, intent(in) :: bcDOF(bcCount)}
285         \args{double precision, intent(in) :: bcValue(bcCount)}
286         \args{external matrix-vector product subroutine :: A}
287         \args{TYPE(varies), pointer :: ptr }
288         \args{double precision, intent(in) :: b(width,length)}
289         \args{double precision, intent(inout) :: x(width,length)}
290
291 Like \kw{IFEM\_Solve\_shared}, this routine solves the linear system $A x = b$ for the unknown vector $x$.  This routine, however, adds support for boundary conditions associated with $x$. These so-called "essential" boundary conditions restrict the values of some unknowns. For example, in structural dynamics, a fixed displacment is such an essential boundary condition.  
292
293 The only form of boundary condition currently supported is to impose a fixed value on certain unknowns, listed by their degree of freedom--that is, their entry in the unknown vector.  In general, the $i$'th DOF of node $n$ has DOF number $n*width+i$ in C and $(n-1)*width+i$ in Fortran.  The framework guarantees that, on output, for all $bcCount$ boundary conditions, $x(bcDOF(f))=bcValue(f)$.  
294
295 For example, if $width$ is 3 in a 3d problem, we would set node $ny$'s y coordinate to 4.6 and node $nz$'s z coordinate to 7.3 like this:
296
297 \begin{alltt}
298 // C++ Example
299   int bcCount=2;
300   int bcDOF[bcCount];
301   double bcValue[bcCount];
302   // Fix node ny's y coordinate
303   bcDOF[0]=ny*width+1; // y is coordinate 1
304   bcValue[0]=4.6;
305   // Fix node nz's z coordinate
306   bcDOF[1]=nz*width+2; // z is coordinate 2
307   bcValue[1]=2.0;
308
309 ! F90 Example
310 // C++ Example
311   integer :: bcCount=2;
312   integer :: bcDOF(bcCount);
313   double precision :: bcValue(bcCount);
314   // Fix node ny's y coordinate
315   bcDOF(1)=(ny-1)*width+2; // y is coordinate 2
316   bcValue(1)=4.6;
317   // Fix node nz's z coordinate
318   bcDOF(2)=(nz-1)*width+3; // z is coordinate 3
319   bcValue(2)=2.0;
320 \end{alltt}
321
322
323
324 Mathematically, what is happening is we are splitting the partially unknown vector $x$ into a completely unknown portion $y$ and a known part $f$:
325 \[ A x = b \]
326 \[ A (y + f) = b \]
327 \[ A y = b - A f \]
328
329 We can then define a new right hand side vector $c=b-A f$ and solve the new linear system $A y=c$ normally.  Rather than renumbering, we do this by zeroing out the known portion of $x$ to make $y$.  The creation of the new linear system, and the substitution back to solve the original system are all done inside this subroutine.
330
331 One important missing feature is the ability to specify general linear constraints on the unknowns, rather than imposing specific values.
332
333
334 \section{Index}
335 \input{index}
336 \end{document}