Intro to Hypre and Useful Linear Solvers

1. Intro to Hypre and Useful Linear Solvers Shule Li Computational Astrophysics Group University of Rochester

2. Part I. Intro to HYPRE (high performance preconditioners) What is hypre? Hypre is a software library of high performance preconditioners for solutions of large and sparse linear systems on massively parallel computers. What are the features of hypre? Scalability. Provides a suite of common iterative methods. Intuitive grid-centric interfaces. (Structured grid, Semi-Structured grid.). Wide choice of options. Language flexibility What can hypre solve particularly in computational astrophysics? Problem that requires additional step of solving sparse matrix can be dealt with in hypre. How to run with hypre in Astrobear? 1. Turn on the elliptic solver flags in domain.data. 2. Turn on the corresponding linear problem flag in physics.data. When working on a problem with diffusion turned on, assign values to kappa and ndif. kappa = 0.0 means no diffusion, ndif = 0.0 means linear diffusion. 3. When compiling, turn on the hypre flag in Make.inc. 4. Compile and run!

3. Work on fixed grid in Hypre How do I build new solver in hypre? The major hypre steps are contained in the subroutine Hypre in the hypreBEAR.f90 file. The solver is built within this subroutine. Hypre calls subroutine StructHypre or SStructHypre, depending on which interface is chosen. For fixed grid, work in StructInterface. How do I pass parameters and variables into hypre? A struct “hypreinfo” is used to define parameters when calling hypre. It has the following members: Interface: choose which interface is used in hypre: struct or semi-struct. Solver: choose what kind of solver to use. Currently PCG is used. qvarIn: which variable in q is passed in. qvarOut: which variable in q is passed out. InterpOrder: order of interpolation between AMR levels Tolerance: tolerance of iteration steps. PrintLevel: Print on screen informations. maxIters: maximum number of iterations hVerbosity: how verbose hypre should be. These parameters can be assigned in bearez.f90 file. The rest of the physical constants should be defined in globaldeclarations.f90 or physics.data.

4. What are the main steps in StructHypre? StructHypre works under fixed grid. The main steps are: Set up the grid. Set up the stencil. Set up the matrix. Set up the right vector. Call the solver, get the result. When building a new solver in hypre, just follow these steps in StructHypre subroutine. matrix right(variable) vector A S = B solution vector How to set up the grid? The following functions are called: (usually not necessary to change.)‏ 1. StructGridCreate 2,3. StructGridSetExtents Apply BC 4.StructGridAssemble For periodic boundary, call StructGridSetPeriodic.

5. What is a stencil? A stencil defines the adjacent points around a center. In 3D, the labels used in hypre are: iEntry 0:center 1:left 2:right 3:bottom 4:top 5:back 6:front offsets3D (3, 7 array)‏ 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1 For instance, when solving equation We assign stencil values: u(0) = -4, u(1) = 1, u(2) = 1, u(3) = 1, u(4) = 1. How do I assign stencil entries? In the code, it is realized by the following lines. Usually we do not need to change anything here when solving a new problem. The actual stencil values are assigned when setting up the matrix. npts=2*ndim+1 do ientry=0:npts-1 // call StructStencilSetElement enddo

6. How do I set the matrix values? The matrix is set up by making the following calls: StructMatrixCreate no need to change StructMatrixInitialize no need to change assign matrix coefficients here. StructMatrixSetBoxValues StructMatrixAssemble no need to change As the graph shows, for each new solver, a new StructMatrixSetBoxValues should be written.

7. How do I write StructMatrixSetBoxValues? Just initialize an 1D array called matrixValues, assign the coefficients and then call the subroutine StructMatrixSetBoxValues. Why I assign coefficients to an array not a matrix? What is the rule of assigning its values? The 1D array holds all the information for grid points and their stencils. It works as (if 2D): matrixValues(0) 　⇒　central point A matrixValues(1) 　⇒　stencil 1 of A matrixValues(2) 　⇒　stencil 2 of A matrixValues(3) 　⇒　stencil 3 of A matrixValues(4) 　⇒　stencil 4 of A matrixValues(5) 　⇒　stencil 5 of A matrixValues(6) 　⇒　stencil 6 of A matrixValues(7) 　⇒　central point B matrixValues(8) 　⇒　stencil 1 of B matrixValues(9) 　⇒　stencil 2 of B matrixValues(10)　⇒　stencil 3 of B matrixValues(11)　⇒　stencil 4 of B matrixValues(12)　⇒　stencil 5 of B matrixValues(12)　⇒　stencil 6 of B . . . . So if the stencil has the npoints points on it where npoints = 2*ndim+1 The dimension of this array should be npoints*ncells, where ncells = mx*my*mz. In the example on page 5, the values should be: m0 = -4, m5 = -4, m11 = -4, m17 = -4 … all the rest values are 1.

8. What if the matrix central values and stencil values are not the same, say, dependent on the location? You do the same thing as before, the following lines can be used: Allocate(matrixValues(:))‏ m = 0 do k = 1, mz do j = 1, my do i = 1, mx matrixValues(m) = f(i,j,k) // central value do n = m+1, m+npoints-1 nn = n-m // label of stencil select case(nn)‏ case(1) // stencil 1 matrixValues(n) = f(i,j,k)‏ . . . Function (i,j,k) can be defined arbitrarily. The above case assigns coefficients in a order of x-y-z. It is better to start m index from 0 (tried diffusion with m starting from 1, also works) . How do I set up the vectors? The subroutines on page 8 should be called. Usually you only modify the subroutine StructVectorSetBoxValues. It works similar to the StructMatrixSetBoxValues, but here you need to setup two vectors, the right (variable) vector and the solution vector.

9. How do I assign right vector values in StructVectorSetBoxValues? Allocate an array called vectorValues, assign values in the same manner as of the matrix is assigned. For instance, in 2D solver, if matrixValues(0:4) is assigned to point (1,1) in the grid, then vectorValues(0) should also correspond to point (1,1). The dimension of the vectorValues should equal to ncells, where ncells = mx*my*mz. Then set it to variable vector: call C_StructVectorSetBoxValues(SvariableVector, …)‏ StructVectorCreate StructVectorIntialize How do I assign solution vector values in StructVectorSetBoxValues? Allocate an array called vectorValues, assign values in the same manner as of the matrix is assigned. In 2D solver, if matrixValues(0:4) is assigned to point (1,1) in the grid, then vectorValues(0) should also correspond to point (1,1). The dimension of the vectorValues should equal to ncells, where ncells = mx*my*mz. Here, since the solution vector is unknown, we can simply assign it to be the solution vector in the previous time step. For instance, in solving the density diffusion, we can set: vectorValues(m) = Info%q(i,j,k,qvarOut). Then set it to solution vector: call C_StructVectorSetBoxValues(SsolutionVector, …)‏ StructVectorSetBoxValues assign vector values here. StructVectorAssemble How do I call the solver? Use the following steps to call the solver and return the solution vector. StructPCGCreate tol=HypreInfo%tolerance // set parameters StructPCGSetup StructPCGSolve

10. How do I get the solution vector? After calling the solver, call StructVectorGetBoxValues. Initialize vectorValues(m) with m from 0 to ncells-1, assign the solution vector to vectorValues(m), and then pass the values of vectorValues(m) out. Example: allocate(vectorValues(0:ncells-1))‏ call C_StructVectorGetBoxValues(SsolutionVector, …)‏ m = 0 do k = 1, mz do j = 1, my do i = 1, mx Info%q(i,j,k,1,qvarOut) = vectorValues(m)‏ . . . What should I do after getting the solution vector? Once the solution is passed out, simply destroy the grid, stencil, matrix, variable and solution vector. (#・∀・)ノ finished!

11. Semi-Struct Interface Semi-Struct Interface is if the grids are not entirely structured, e.g. block-structured grids, AMR applications. to be continued …

12. Part II. Useful Linear Solvers - Diffusion Problem An Implicit Scheme The 1D diffusion equation reads: where D is the thermal conductivity: To solve this equation, we first define: or equivalently, Substitute D into the diffusion equation we have: It is just:

13. Now write the differentials as finite differences on both sides, we have: which is equivalent to: where This simple explicit method has a stability problem. Let’s look at the simplest case: linear diffusion. In this case Thus we have the differencing scheme to be: Consider the method suggested by Von Neumann, we assume the variables are single mode planar waves on a periodic grid: Then the differencing scheme becomes:

14. Regrouping terms we get: The stability requires: Thus for all ka value we require: we thus require h > 2: In other words dt should be smaller than the thermal conduction time for each grid. For applications with finer grid spacing, this requirement is devastatingly slow and inefficient. We are thus in need of a scheme with looser stability requirement. This kind of method is usually an implicit method. Now look back at the differencing scheme for linear diffusion on page 13, we can make it implicit by substituting the variables at step n on the right hand side by variables at step “n+1/2”:

15. Now we can check the stability. After some math we have: Obviously this is always stable if h > 0. This scheme is called Crank Nicholson scheme. For nonlinear diffusion, we use z instead of T on the right hand side. Here To avoid having to solve a set of nonlinear equations we do the following expansion: Substitute the above expression for z into the differencing scheme, and regroup terms we finally have:

16. This is a tri-diagonal matrix, which is easy to solve by backward iteration. We can assign values: Then the scheme becomes: Boundary Condition in Crank Nicholson scheme. When the boundaries are transmissive, the equations are not closed. Thus when i=1, we should include the ghost region values to close the equations: We thus have to extend the computation domain and compute ghost region values explicitly during the calculation, which is not desired. But obviously, this can be avoided by combining the above two equations:

17. Similarly, at the other end we have: In hypre, we can simply set the matrix with values a, b, c and the right vector with value d, and then pass them into the linear solver. 2D Crank Nicholson Scheme The above 1D scheme can be extended to higher dimensions by using sweep method easily. We first do calculation on x direction with conductivity Dx for time step dt pretending no diffusion on y, and then do calculation on y direction with conductivity Dy for time step dt pretending no diffusion on x. What you should do is defining a public variable, for example, hflag, and then pass it around in the “SetBoxValues “ subroutines. For instance, in Hypre call, define: integer hflag common /sweep/hflag select case(HypreInfo%Interface)‏ case(StructInterface)‏ hflag = 0 call StructHypre(HypreInfo)‏ hflag = 1 call StructHypre(HypreInfo)‏ . In the diffusion solver, the matrix and vector values should change according to which direction you are solving. But obviously, we do not need a sweep method for multi-dimension problem since the matrix is solved by hypre. We can simply decompose the 2D equation and solve it at once. For example, if the eqution is:

18. We can directly write it into difference form with Crank Nicholson method: Grouping terms we finally get: In terms of stencil, we have: except for the boundaries.

19. At the boundaries, we use the same method as before. For instance, if i=1, we set: the rest elements remain the same. To Verify the code, we can use an arbitrary function, fft it, take each frequency components and do a corresponding exponential decay. But remember: it's NOT a Gaussian! Gaussian is the solution for linear diffusion only when the Dirichlet boundaries are placed at -inf to +inf! More Complicated Boundary Conditions need fft method. coming soon.