Added this new manual for IFEM solvers.
authorOrion Lawlor <olawlor@acm.org>
Mon, 10 Mar 2003 19:49:40 +0000 (19:49 +0000)
committerOrion Lawlor <olawlor@acm.org>
Mon, 10 Mar 2003 19:49:40 +0000 (19:49 +0000)
doc/ifem/Makefile [new file with mode: 0644]
doc/ifem/manual.tex [new file with mode: 0644]

diff --git a/doc/ifem/Makefile b/doc/ifem/Makefile
new file mode 100644 (file)
index 0000000..bc5944d
--- /dev/null
@@ -0,0 +1,8 @@
+# Stub makefile for LaTeX PPL manual
+FILE=manual
+TEX=$(FILE).tex index.tex
+DEST=ifem
+LATEX2HTML=$(L2H) -split 4
+PROJECT_LINK='<a href="http://charm.cs.uiuc.edu/research/fem">FEM Homepage</a><br>'
+
+include ../Makefile.common
diff --git a/doc/ifem/manual.tex b/doc/ifem/manual.tex
new file mode 100644 (file)
index 0000000..07be4da
--- /dev/null
@@ -0,0 +1,336 @@
+\documentclass[10pt]{article}
+\usepackage{../pplmanual}
+\input{../pplmanual}
+
+\makeindex
+
+\title{\charmpp\\ Iterative Finite Element Matrix (IFEM) Library\\ Manual}
+\version{1.2}
+\credits{
+Initial version of \charmpp{} Finite Element Framework was developed
+by Orion Lawlor in the spring of 2003.
+}
+
+\begin{document}
+
+\maketitle
+
+\section{Introduction}
+
+This manual presents the Iterative Finite Element Matrix (IFEM) library,
+a library for easily solving matrix problems derived from 
+finite-element formulations.  The library is designed to be matrix-free,
+in that the only matrix operation required is matrix-vector product,
+and hence the entire matrix need never be assembled
+
+IFEM is built on the mesh and communication capabilities of the Charm++ 
+FEM Framework, so for details on the basic runtime, problem setup, and 
+partitioning see the FEM Framework manual.
+
+
+\subsection{Terminology}
+
+A FEM program manipulates elements and nodes. An \term{element} is a portion of
+the problem domain, also known as a cell, and is typically some simple shape 
+like a triangle, square, or hexagon in 2D; or tetrahedron or rectangular solid in 3D.  
+A \term{node} is a point in the domain, and is often the vertex of several elements.  
+Together, the elements and nodes form a \term{mesh}, which is the 
+central data structure in the FEM framework.  See the FEM manual for details.
+
+% --------------- Solvers -----------------
+\section{Solvers}
+\label{sec:solver}
+A IFEM \term{solver} is a subroutine that controls the search for the solution.
+
+Solvers often take extra parameters, which are listed in a type called 
+in C \kw{ILSI\_Param}, which in Fortran is an array of \kw{ILSI\_PARAM} doubles.  
+You initialize these solver parameters using the subroutine \kw{ILSI\_Param\_new}, 
+which takes the parameters as its only argument.  The input and output
+parameters in an \kw{ILSI\_Param} are listed in Table~\ref{table:solverIn}
+and Table~\ref{table:solverOut}.
+
+\begin{table}[hh]
+\begin{center}
+\begin{tabular}{|l|l|l|}\hline
+C Field Name & Fortran Field Offset & Use\\\hline
+maxResidual & 1 & If nonzero, termination criteria: magnitude of residual. \\
+maxIterations & 2 & If nonzero, termination criteria: number of iterations. \\
+solverIn[8] & 3-10 & Solver-specific input parameters. \\
+\hline
+\end{tabular}
+\end{center}
+\caption{\kw{ILSI\_Param} solver input parameters.}
+\label{table:solverIn}
+\end{table}
+
+\begin{table}[hh]
+\begin{center}
+\begin{tabular}{|l|l|l|}\hline
+C Field Name & Fortran Field Offset & Use\\\hline
+residual & 11 & Magnitude of residual of final solution. \\
+iterations & 12 & Number of iterations actually taken. \\
+solverOut[8] & 13-20 & Solver-specific output parameters. \\
+\hline
+\end{tabular}
+\end{center}
+\caption{\kw{ILSI\_Param} solver output parameters.}
+\label{table:solverOut}
+\end{table}
+
+
+\subsection{Conjugate Gradient Solver}
+
+The only solver currently written using IFEM is the conjugate gradient solver.
+This linear solver requires the matrix to be real, symmetric and positive definite.
+
+Each iteration of the conjugate gradient solver requires one matrix-vector product
+and two global dot products.  For well-conditioned problems, the solver typically 
+converges in some small multiple of the diameter of the mesh--the number of elements
+along the largest side of the mesh.
+
+You access the conjugate gradient solver via the subroutine name \kw{ILSI\_CG\_Solver}.
+
+
+% -------------- Shared-node ----------------
+\section{Solving Shared-Node Systems}
+
+Many problems encountered in FEM analysis place the entries of the
+known and unknown vectors at the nodes--the vertices of the domain.
+Elements provide linear relationships between the known and unknown node values, 
+and the entire matrix expresses the combination of all these element relations.
+
+For example, in a structural statics problem, we know the net force at
+each node, $f$, and seek the displacements of each node, $u$.  Elements 
+provide the relationship, often called a stiffness matrix $K$, between 
+a nodes' displacments and its net forces:
+
+\[
+       f=K u
+\]
+
+We normally label the known vector $b$ (in the example, the forces), the unknown
+vector $x$ (in the example, the displacments), and the matrix $A$:
+
+\[
+       b=A x
+\]
+
+
+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.
+
+
+\newpage
+\subsection{IFEM\_Solve\_shared}
+
+\prototype{IFEM\_Solve\_shared}
+\function{void IFEM\_Solve\_shared(ILSI\_Solver s,ILSI\_Param *p,
+       int fem\_mesh, int fem\_entity,int length,int width,
+       IFEM\_Matrix\_product\_c A, void *ptr, const double *b, double *x);}
+\function{subroutine IFEM\_Solve\_shared(s,p,
+       fem\_mesh,fem\_entity,length,width,
+       A,ptr,b,x)}
+       \args{external solver subroutine :: s}
+       \args{double precision, intent(inout) :: p(ILSI\_PARAM)}
+       \args{integer, intent(in) :: fem\_mesh, fem\_entity, length,width}
+       \args{external matrix-vector product subroutine :: A}
+       \args{TYPE(varies), pointer :: ptr }
+        \args{double precision, intent(in) :: b(width,length)}
+        \args{double precision, intent(inout) :: x(width,length)}
+
+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.
+
+\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.
+
+\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.  
+
+\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)}.
+
+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.
+
+\begin{alltt}
+// C++ Example
+  int mesh=FEM_Mesh_default_read();
+  int nNodes=FEM_Mesh_get_length(mesh,FEM_NODE);
+  int width=3; //A 3D problem
+  ILSI_Param solverParam;
+  struct myProblemData myData;
+  
+  double *b=new double[nNodes*width];
+  double *x=new double[nNodes*width];
+  ... prepare solution target b and guess x ...
+  
+  ILSI_Param_new(&solverParam);
+  solverParam.maxResidual=1.0e-4;
+  solverParam.maxIterations=500; 
+  
+  IFEM_Solve_shared(IFEM_CG_Solver,&solverParam,
+         mesh,FEM_NODE, nNodes,width,
+         myMatrixVectorProduct, &myData, b,x);
+  
+! F90 Example
+  include 'ifem\_f.h'
+  INTEGER :: mesh, nNodes,width
+  DOUBLE PRECISION, ALLOCATABLE :: b(:,:), x(:,:)
+  DOUBLE PRECISION :: solverParam(ILSI_PARAM)
+  TYPE(myProblemData) :: myData
+  
+  mesh=FEM_Mesh_default_read()
+  nNodes=FEM_Mesh_get_length(mesh,FEM_NODE)
+  width=3   ! A 3D problem
+  
+  ALLOCATE(b(width,nNodes), x(width,nNodes))
+  ... prepare solution target b and guess x ..
+  
+  ILSI_Param_new(&solverParam);
+  solverParam(1)=1.0e-4;
+  solverParam(2)=500; 
+  
+  IFEM_Solve_shared(IFEM_CG_Solver,solverParam,
+         mesh,FEM_NODE, nNodes,width,
+         myMatrixVectorProduct, myData, b,x);
+
+\end{alltt}
+
+\newpage
+\subsubsection{Matrix-vector product routine}
+\label{sec:mvp}
+
+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:
+
+\prototype{IFEM\_Matrix\_product}
+\function{void IFEM\_Matrix\_product(void *ptr,int length,int width,
+       const double *src, double *dest);}
+\function{subroutine IFEM\_Matrix\_product(ptr,length,width,src,dest)}
+  \args{TYPE(varies), pointer  :: ptr}
+  \args{integer, intent(in)  :: length,width}
+  \args{double precision, intent(in) :: src(width,length)}
+  \args{double precision, intent(out) :: dest(width,length)}
+
+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.
+
+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.
+
+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.
+
+\begin{alltt}
+// C++ Example
+typedef struct \{
+  int nElements; //Number of local elements
+  int *conn; // Nodes adjacent to each element: 2*nElements entries
+  double k; //Uniform spring constant
+\} myProblemData;
+
+void myMatrixVectorProduct(void *ptr,int nNodes,int dofPerNode,
+          const double *src,double *dest) 
+\{
+  myProblemData *d=(myProblemData *)ptr;
+  int n,e;
+  // Zero out output force vector:
+  for (n=0;n<nNodes;n++) dest[n]=0;
+  // Add in forces from local elements
+  for (e=0;e<d->nElements;e++) {
+    int n1=d->conn[2*e+0]; // Left node
+    int n2=d->conn[2*e+1]; // Right node
+    double f=d->k * (src[n2]-src[n1]); //Force
+    dest[n1]+=f;
+    dest[n2]-=f;
+  }
+\}
+
+
+! F90 Example
+  TYPE(myProblemData) 
+    INTEGER :: nElements
+    INTEGER, ALLOCATABLE :: conn(2,:)
+    DOUBLE PRECISION :: k
+  END TYPE
+  
+SUBROUTINE myMatrixVectorProduct(d,nNodes,dofPerNode,src,dest)
+  include 'ifem\_f.h'
+  TYPE(myProblemData), pointer :: d
+  INTEGER :: nNodes,dofPerNode
+  DOUBLE PRECISION :: src(dofPerNode,nNodes), dest(dofPerNode,nNodes)
+  INTEGER :: e,n1,n2
+  DOUBLE PRECISION :: f
+  
+  dest(:,:)=0.0
+  do e=1,d%nElements
+    n1=d%conn(1,e)
+    n2=d%conn(2,e)
+    f=d%k * (src(1,n2)-src(1,n1))
+    dest(1,n1)=dest(1,n1)+f
+    dest(1,n2)=dest(1,n2)+f
+  end do
+END SUBROUTINE
+\end{alltt}
+
+
+\newpage
+\subsection{IFEM\_Solve\_shared\_bc}
+
+\prototype{IFEM\_Solve\_shared\_bc}
+\function{void IFEM\_Solve\_shared\_bc(ILSI\_Solver s,ILSI\_Param *p,
+       int fem\_mesh, int fem\_entity,int length,int width,
+       int bcCount, const int *bcDOF, const double *bcValue,
+       IFEM\_Matrix\_product\_c A, void *ptr, const double *b, double *x);}
+\function{subroutine IFEM\_Solve\_shared\_bc(s,p,
+       fem\_mesh,fem\_entity,length,width,
+       bcCount,bcDOF,bcValue,
+       A,ptr,b,x)}
+       \args{external solver subroutine :: s}
+       \args{double precision, intent(inout) :: p(ILSI\_PARAM)}
+       \args{integer, intent(in) :: fem\_mesh, fem\_entity, length,width}
+       \args{integer, intent(in) :: bcCount}
+       \args{integer, intent(in) :: bcDOF(bcCount)}
+       \args{double precision, intent(in) :: bcValue(bcCount)}
+       \args{external matrix-vector product subroutine :: A}
+       \args{TYPE(varies), pointer :: ptr }
+        \args{double precision, intent(in) :: b(width,length)}
+        \args{double precision, intent(inout) :: x(width,length)}
+
+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.  
+
+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)$.  
+
+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:
+
+\begin{alltt}
+// C++ Example
+  int bcCount=2;
+  int bcDOF[bcCount];
+  double bcValue[bcCount];
+  // Fix node ny's y coordinate
+  bcDOF[0]=ny*width+1; // y is coordinate 1
+  bcValue[0]=4.6;
+  // Fix node nz's z coordinate
+  bcDOF[1]=nz*width+2; // z is coordinate 2
+  bcValue[1]=2.0;
+
+! F90 Example
+// C++ Example
+  integer :: bcCount=2;
+  integer :: bcDOF(bcCount);
+  double precision :: bcValue(bcCount);
+  // Fix node ny's y coordinate
+  bcDOF(1)=(ny-1)*width+2; // y is coordinate 2
+  bcValue(1)=4.6;
+  // Fix node nz's z coordinate
+  bcDOF(2)=(nz-1)*width+3; // z is coordinate 3
+  bcValue(2)=2.0;
+\end{alltt}
+
+
+
+Mathematically, what is happening is we are splitting the partially unknown vector $x$ into a completely unknown portion $y$ and a known part $f$:
+\[ A x = b \]
+\[ A (y + f) = b \]
+\[ A y = b - A f \]
+
+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.
+
+One important missing feature is the ability to specify general linear constraints on the unknowns, rather than imposing specific values.
+
+
+\section{Index}
+\input{index}
+\end{document}