NumMethodsPDEs - Revision history

GethinWilliams: /* Practice: Guassian Elimination */

2011-11-08T17:30:11Z

Practice: Guassian Elimination

GethinWilliams: /* Practice: Guassian Elimination */

2011-11-08T17:28:22Z

Practice: Guassian Elimination

GethinWilliams: /* Practice: Guassian Elimination */

2011-11-08T17:26:38Z

Practice: Guassian Elimination

GethinWilliams: /* Practice: Solving Eliptical Problems using LAPACK */

2011-10-31T13:02:59Z

Practice: Solving Eliptical Problems using LAPACK

GethinWilliams: /* Practice */

2011-10-31T12:59:45Z

Practice

GethinWilliams: /* Practice: Solving Parabolic Problems using LAPACK */

2011-10-31T12:57:54Z

Practice: Solving Parabolic Problems using LAPACK

GethinWilliams: /* Practice: Guassian Elimination */

2011-10-31T12:56:30Z

Practice: Guassian Elimination

GethinWilliams at 12:54, 31 October 2011

2011-10-31T12:54:39Z

GethinWilliams: /* Practice: Solving Parabolic Problems using LAPACK */

2011-10-31T12:51:21Z

Practice: Solving Parabolic Problems using LAPACK

GethinWilliams at 12:47, 31 October 2011

2011-10-31T12:47:03Z

@@ Line 85: / Line 85: @@
 # top row x2, subtract from 3rd row
 # dived the 2nd row by 4
-# multiply the 2nd row by 2, add the the 3rd row
+# multiply the 2nd row by 2, add to the 3rd row
 etc:

@@ Line 84: / Line 84: @@
 # top row x3, subtract from 2nd row
 # top row x2, subtract from 3rd row
 etc:

@@ Line 81: / Line 81: @@
 </math>
-The following computation shows Gauss-Jordan elimination applied a further simplified (assumed the x, y and z) matrix version ('''Ax'''='''b''') of the above system:
+The following computation shows Gauss-Jordan elimination applied a further simplified (assumed the x, y and z) matrix version ('''Ax'''='''b''') of the above system, where the operations include:
 :<math>\left[\begin{array}{rrr|r}
 & 3 & -2 & 5 \\
@@ Line 122: / Line 126: @@
 & 0 & 1 & 2
 \end{array}\right].</math>
 OK, so far so good?  Let's try our first code example:

@@ Line 289: / Line 289: @@
 ==Practice: Solving Eliptical Problems using LAPACK==
-Next, take a look at '''simple-laplace.f90'''.  This program solves the 2D laplace equation for a 4x4 domain, as described on [[NumMethodsPDEs]].  This program only differs from the previous one in way in which we initialise the matrix '''A''' and vector '''b'''.  To run the program type:
+Next, take a look at '''simple-laplace.f90'''.  This program solves the 2D laplace equation for a 4x4 domain, as described above.  This program only differs from the previous one in way in which we initialise the matrix '''A''' and vector '''b'''.  To run the program type:
 <pre>
-./simple-laplace-expert.exe
+./simple-laplace.exe
 </pre>
@@ Line 304: / Line 304: @@
 <pre>
-./mxm-laplace.exe
+./mxm-laplace-expert.exe
 </pre>

@@ Line 16: / Line 16: @@
 ==Practice==
-In this tutorial we'll take a look at a number of example programs which solve equations using linear algebra packages.  These packages include LAPACK, PETSc, PLASMA and MAGMA.  The examples will be based upon the solution of a PDE using the Finite Difference Method (FDM) (although we'll also touch upon the theory for an alternative approach--the Finite Element Method (FEM).
+In this tutorial we'll take a look at a number of example programs which solve equations using linear algebra packages.  These packages include LAPACK, PETSc, PLASMA and MAGMA.  The examples will be based upon the solution of a PDE using the Finite Difference Method (FDM), although we'll also touch upon the theory for an alternative approach--the Finite Element Method (FEM).
 To download the code for the examples, type:

@@ Line 797: / Line 797: @@
 In '''simple-crank-nicolson.f90''' we see that the initial call to '''dgesv''' replaces the entries in the matrix A with an encoding of the LU decompostion (the unit diagonal is assumed).  For subsequent timesteps we can thus call '''dgetrs''' instead:
-<pre>
+<source lang="fortran">
 do ii=2,4
    ...
@@ Line 803: / Line 803: @@
    ...
 end do
-</pre>
+</source>
 where we have remembered to take account of the boundary conditions at each time-step too:
-<pre>
+<source lang="fortran">
 do ii=2,4
    ...
@@ Line 828: / Line 828: @@
    ...
 end do
-</pre>
+</source>
 ----

@@ Line 136: / Line 136: @@
 To see how to solve the above matrix example using LAPACK take a look inside '''dgesv-example.f90'''.  First the declarations:
-<pre>
+<source lang="fortran">
    integer, parameter    :: N  = 3
    integer, parameter    :: LDA  = N  ! leading dimension of A
@@ Line 146: / Line 146: @@
    real(kind=8), dimension(LDA,N)    :: A ! LDAxN matrix
    real(kind=8), dimension(LDB,NRHS) :: B ! LDBxNRHS matrix
-</pre>
+</source>
 Where, we see that the ''leading dimension'' of a matrix is used to specify the number of rows of storage.  Note also, that we could simultaneously solve for multiple vectors B, by specifying NRHS > 1.
 Next, we set up the matrix A:
-<pre>
+<source lang="fortran">
    ! ( 1  3 -2)
    ! ( 3  5  6)
@@ Line 164: / Line 164: @@
    A(3,2) =  4.0d+0
    A(3,3) =  3.0d+0
-</pre>
+</source>
 And the vector B:
-<pre>
+<source lang="fortran">
    ! ( 5)
    ! ( 7)
@@ Line 174: / Line 174: @@
    B(2,1) =  7.0d+0
    B(3,1) =  8.0d+0
-</pre>
+</source>
 Once we have all that in place, all we need to do is call the solver routine:
-<pre>
+<source lang="fortran">
    call dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
-</pre>
+</source>
 Note that the solution vector ''replaces'' the values in B.  (The contents of A will be changed in place too, but more of that later.)
-<pre>
+<source lang="fortran">
    ! Print the result vector X.
    ! It should be:
@@ Line 193: / Line 193: @@
       write(*,*) B(ii,1)
    end do
-</pre>
+</source>
 Run the program to satisfy yourself that it works:

@@ Line 914: / Line 914: @@
 For further examples see INSTALL_DIR/petsc-3.1-p1/src/ksp/ksp/examples/tutorials.
-=Theory: Other Numerical Methods=
+=Other Numerical Methods=
 We have looked at the ''finite difference'' method (FDM) above, but other (more mathsy) approaches exist, such as the ''finite element'' and ''finite volume'' methods.  Each approach has its pros and cons.  For example, the FDM is limited to regular grids spread over rectangular domains, but is easy to program.  In contrast the FEM is applicable to intricately shaped domains and the resolution of the individual elements can be varied, but is requires more involved programming and a deeper level of understand of the relevant mathematics.  It should also be noted, however, that under certain circumstances, the FEM and FVM boil down to be ''exactly the same'' as the FDM, even though we've approached the task of solving a PDE from a different direction.
@@ Line 920: / Line 920: @@
 One thing which is ''always'' common to the approaches is that ''the problem boils down to a system of linear equations to solve, and that a library of linear algebra routines will be the workhorse when finding solutions to PDEs on a computer''.  Repeat after me, "a fast linear algebra library is key!"
-==The Finite Element Method (FEM)==
+==Theory: The Finite Element Method (FEM)==
 The general approach here is to create a potential function from our PDE and to find a solution to it via minimising the potential function.  Key steps include:
@@ Line 999: / Line 999: @@
 [[Image:Elmer.jpg|500px|centre|thumbnail|A FEM mesh used by Elmer to solve a steady state heat problem.]]
 [[Image:Elmer-solution.jpg|500px|centre|thumbnail|The solution to the steady state heat problem as computed by Elmer.]]
 =Further Reading=
 * [https://source.ggy.bris.ac.uk/wiki/A_Good_Read#Numerical_Methods] Peter J. Olver has made available some excellent notes on the topic of solving PDEs, including a description of the FEM.

← Older revision		Revision as of 12:51, 31 October 2011
Line 835:		Line 835:

	* Experiment with values of <math>r=\frac{k}{h^2}</math> for the 3 schemes, so as to highlight the effects of the CFL condition.		* Experiment with values of <math>r=\frac{k}{h^2}</math> for the 3 schemes, so as to highlight the effects of the CFL condition.
		+
		+	==Practice: Banded Matrices==
		+
		+	In all my examples, thus far, I've called routines from the '''dge'''xxx series offered by LAPACK. The '''d''' stands for '''d'''ouble and the '''ge''' for '''ge'''nerel purpose. For banded matrices, LAPACK offers a much more efficient format and routines to match. Take a look at http://www.netlib.org/lapack/double/ for a list of driver routines using matrices of double precision numbers, or http://www.netlib.org/lapack for more documentation in general.
		+
		+	==Practice: PETSc & Iterative solutions and sparse matrices==
		+
		+	In order to compute the solution of a linear system via the direct method of an LU-decomposition, we know that it will cost <math>O(n^3)</math> operations. Can we compute a solution more cheaply? The class of iterative methods strive to compute a cheaper solution by first guessing a value for the solution vector <math>x</math> and then iteratively refining it. Thus we will get a series of (hopefully!) improving approximations:
		+
		+	:<math>
		+	x^0, x^1, x^2, x^3, \cdots, x^j \rightarrow x
		+	</math>
		+
		+	By way of an illustration of such a method, consider splitting matrix A into 3 others:
		+
		+	:<math>
		+	A = L + D + U
		+	</math>
		+
		+	Note the summation, and so the L and the U matrices are '''not the same as those considered in the LU-decomposition''' used by the direct approach. With this split, '''Ax'''='''b''', becomes:
		+
		+	:<math>
		+	Lx + Dx + Ux = b
		+	</math>
		+
		+	The game now is to replace some occurrences of <math>x</math> by <math>x^k</math> and others by <math>x^{k+1}</math>, and thereby obtaining a relationship which tells us how to get an approximation at one iteration step, from the approximation at the previous step. If the resulting iterative scheme converges, then we're in business!
		+
		+	One classic form of such an iteration is the Jacobi method (which is useful when A is diagonally dominant):
		+
		+	:<math>
		+	Lx^j + Dx^{j+1} + Ux^j = b
		+	</math>
		+
		+	:<math>
		+	x^{j+1} = D^{-1}(b - (L + U)x^j)
		+	</math>
		+
		+	To take a look at a illustrative implementation of the Jacobi method using LAPACK:
		+
		+	<pre>
		+	cd ../example3
		+	make
		+	</pre>
		+
		+	and take a look at '''jacobi-example.f90'''. To run the executable, type:
		+
		+	<pre>
		+	./jacobi-example.exe
		+	</pre>
		+
		+	Happily, however, we do not need to write our own iterative schemes, as the good folk at Argonne have written the '''PETSc''' library.
		+
		+	http://www.mcs.anl.gov/petsc/petsc-as/
		+
		+	Take a look inside '''ex1f.F''' to see how to solve our laplacian example using PETSc. To run the program you will obviously need PETSc installed on the computer you are using. If you have, then edit Makefile and set PETSC_DIR to the appropriate value, and then type:
		+
		+	<pre>
		+	make ex1f
		+	qsub mpi_submit
		+	</pre>
		+
		+	Look in the file called '''OUT''' for the results of running PETSc.
		+
		+	We can again plot a graph of timings given varying N:
		+
		+	[[Image:gmres-vs-LU.jpg\|450px\|centre\|thumbnail\|Comparing the time taken by an iterative method to that taken by LAPACK's LU decomposition for a sparse matrix.]]
		+
		+	We can see that an iterative method can scale approximately <math>O(n^2)</math>. For large matrices, this can mean that a solution is found several orders of magnitude faster. Now that is a speed-up worth having! It's liking buying a supercomputer than runs several orders of magnitude faster than your old one..
		+
		+	Better yet, you can try different schemes for iterative solution, just by varying some command-line arguments. For example:
		+
		+	<pre>
		+	./ex1f -ksp_type bcgs -pc_type jacobi
		+	</pre>
		+
		+	will invoke the biconjugate gradient stabalised scheme together with a jacobi preconditioner matrix. (This in fact runs faster than the GMRES method for the above case--808ms for N=9694. Note, however, that it cannot scale better than <math>O(N^2)</math>.) For those interested, the PETSc manual is the best next port of call for more details on the various methods made available by the package.
		+
		+	For further examples see INSTALL_DIR/petsc-3.1-p1/src/ksp/ksp/examples/tutorials.

	=Theory: Other Numerical Methods=		=Theory: Other Numerical Methods=

@@ Line 641: / Line 641: @@
 (*The two difference schemes for either side of the equation were chosen so that the approximation error characteristics for either side are approximately balanced and so the errors on one side or the other do not begin to dominate the calculations.)
-==Solution Stability and the CFL Condition==
+==Theory: Solution Stability and the CFL Condition==
 The Forward Euler scheme has the benefit of simplicity, but is unstable for large time steps.  The CFL (Courant–Friedrichs–Lewy) condition states that the following must hold:
@@ Line 649: / Line 649: @@
 for stability.
-==Parabolic Equations: Backward (Implicit) Euler==
+==Theory: Parabolic Equations & Backward (Implicit) Euler==
 Another way to formulate our finite difference equations for the very same problem is to use backward-differences instead of forward-differences.  The backward Euler scheme requires more computation per time-step than the forward scheme, but has the distinct advantage that it is '''unconditionally stable''' with regard to the relative size between the time and spatial step sizes; i.e. we are not limited to <math>r \le 0.5</math> when we use the backward Euler scheme.
@@ Line 719: / Line 719: @@
 </math>
-==Parabolic Equations: Crank-Nicolson==
+==Theory: Parabolic Equations & Crank-Nicolson==
 Now, forward and backward Euler schemes are by far from the only ways in which to approach this 1D heat problem.  Yet another scheme is the to create finite difference equations using the Crank-Nicolson stencil:
@@ Line 786: / Line 786: @@
 \begin{bmatrix} ru_{11} + (2-2r)u_{12} + ru_{13} \\ \vdots \end{bmatrix}
 </math>
 ----
-'''Exercise'''
+'''Exercises'''
-* The examples on [[LinAlgebraPacks#Solving_Parabolic_Problems_using_LAPACK | solving parabolic problems]] contains an exercise in which you see what happens if the CFL condition is not satisfied.  Then read on through until the end of that page.
+* Experiment with values of <math>r=\frac{k}{h^2}</math> for the 3 schemes, so as to highlight the effects of the CFL condition.
-=Other Numerical Methods=
+=Theory: Other Numerical Methods=
 We have looked at the ''finite difference'' method (FDM) above, but other (more mathsy) approaches exist, such as the ''finite element'' and ''finite volume'' methods.  Each approach has its pros and cons.  For example, the FDM is limited to regular grids spread over rectangular domains, but is easy to program.  In contrast the FEM is applicable to intricately shaped domains and the resolution of the individual elements can be varied, but is requires more involved programming and a deeper level of understand of the relevant mathematics.  It should also be noted, however, that under certain circumstances, the FEM and FVM boil down to be ''exactly the same'' as the FDM, even though we've approached the task of solving a PDE from a different direction.