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