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
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 0center 1left 2right 3bottom 4top
5back 6front
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. npts2ndim1 do
ientry0npts-1 // call StructStencilSetElement en
ddo
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 2ndim1 The dimension of this array
should be npointsncells, where ncells
mxmymz. 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 m1, mnpoints-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(04) 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 mxmymz. 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(04) 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 mxmymz. 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) Infoq(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 tolHypreInfotolerance
// 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(vectorValue
s(0ncells-1))? call C_StructVectorGetBoxValues(Ss
olutionVector, )? m 0 do k 1, mz do j 1,
my do i 1, mx Infoq(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,

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. Lets 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 gt 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 n1/2
15
Now we can check the stability. After some math
we have Obviously this is always stable if
h gt 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 S
ubstitute 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 i1, 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(HypreInfoInterface)? cas
e(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 i1, 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.