Discretizing Elliptic Equations

Discretizing elliptic equations generally involve the discretization of the Laplacian operator. While discretizing the Laplacian operator on a fixed grid is somewhat trivial, the discretization on an nested grid type mesh is not. Here we look at some of the difficulties and how they are handled in Astrobear.

Outline

• Fixed Grid 1D Discretization
• Non-Fixed 1D Stencils
  • Arbitrary stencils
• 2D Fixed Grid Stencils
• 2D Non-Fixed Grid Stencils
  • Using higher level cells
  • Using lower level cells
  • One possible approach
  • One other approach
• 3D AMR

Fixed Grid 1D Discretization

In 1D the Laplacian is simply d²/dx² and when acting on a variable y_i centered at a set of points x_i separated by DX can be written as D²/Dx² (y) = (y_i+1-2y_i+y_i-1)/(DX²) where we have used a capital D for the discrete operator and DX for the grid spacing. Of course any discretization will introduce errors and we can look to see what errors are introduced by performing a taylor expansion of y about x__i. y__i+1=y_i+dy/dx_i DX + ½ d²y/dx²_i DX² + 1/6 d³/dx³_i DX³ + 1/24 d⁴/dy⁴_i DX⁴ and so on. If we plug these taylor expansions into the discretized laplacian expression we find that D²/Dx²=d²/dy²_i+ 1/12 d⁴/dy⁴_i DX² + … so the discretization is effectively second order (even though it is the fourth order term in the taylor expansion since we are interested in the second order term)

It is worth taking a moment to see how we could have calculated this second order expression since we will use it later on. Most simply we can write our discretization as a linear combination of surrounding values (y_i-1 ,y_i ,y_i+1). In principle we could use (y_i-2 and y_i+2) if we wanted an expression that was fourth order but second order is sufficient. If we write D²/Dx² (y_i) as A_i-1 y_i-1 + A_i y_i + A_i+1 y_i+1 , and replace the y's with there taylor expansion (and set DX=1 for simplicity) we will end up with something like (A_i-1+A_i+A_i+1) y_i + (A_i-1-A_i+1) dy/dx_i + ½(A_i-1+A_i+1) d²y/dx²_i + 1/6(A_i+1-A_i-1) d³y/dx³_i + 1/24(A_i+1+A_i-1) dy⁴/dx⁴_i + …

Since we would like this discretization to be 0 y_i + 0 dy/dx_i + 1 d²y/dx²_i + 0 d³y/dx³_i + 0 d⁴y/dx⁴_i + … we can write this as a system of equations for A_i's (where we have set DX to 1)

   1    1    1                              0
  -1    0    1          A_{i-1}             0
  1/2   0   1/2         A_{i}       =       1
 -1/6   0   1/6         A_{i+1}             0
  1/24  0   1/24                            0

Of course this system is over determined so we will have to drop the last two rows (not worrying about third or fourth order terms in the expansion)

   1    1    1          A_{i-1}             0
  -1    0    1          A_{i}       =       0
  1/2   0   1/2         A_{i+1}             1

and and we are left with a system of equations with one solution A__i-1 = 1, A_i=-2, A_i+1=1. Plugging this back in we see that the coefficient for the d³/dx³ term is zero so this method with only three points is accidentally second order. If our points are not evenly spaced we will only achieve first order accuracy with a three point stencil.

Non-Fixed 1D Stencils

Before considering arbitrary spacing we should look at one other option in fixed grid. Instead of using (y_i-1 , y_i , & y_i+1) what if we use (y_i , y_i+1 , & y_i+2)? Our matrix now becomes

   1    1    1          A_{i}               0
   0    1    2          A_{i+1}     =       0
   0   1/2   2          A_{i+2}             1

and our solution vector is A__i = 1, A_i+1=-2, A_i+2=1 (Note if d²y/dx² is a constant this is accurate) - but what about the third and fourth order terms? Before the coefficients were 0 & 1/12. Now they are 1 & 7/12. So we have introduced a large first order error by not centering our stencil. What if we stretch our stencil but keep it centered? Use (y_i-2 , y_i , & y_i+2)? Our third order term again vanishes, but our fourth order term now has a coefficient of 4/3 instead of 1/12. Finally what if we use (y_i-2 , y_i , & y_i+1)? Our matrix equation becomes:

    1    1    1          A_{i}               0
   -2    0    1          A_{i+1}     =       0
    2    0   1/2         A_{i+2}             1

And our solution vector is A__i-2 = 1/3, A_i=-1, A_i+1=2/3. Since we are no longer centered we have a third and fourth order coefficients of -1/3 & ¼. If we don't have data at y_i-1 what is the best option? Here are the stencils and associated third and fourth order terms

 [  0  1  2 ]  ->  [   1   7/12 ]
 [ -2  0  2 ]  ->  [   0   4/3  ]
 [ -2  0  1 ]  ->  [ -1/3  1/4  ]

If we assume y is a fourth order polynomial with a local expansion of a + bx + cx² + dx³ + ex⁴ then d³y/dx³ is 6dx+24ex² and d⁴y/dx⁴ is 24ex. The error for each option (divided by DX³) is therefore

 [  0  1  2 ]  ->  [   1   7/12 ] -> 6dx + 24ex^2^ + 14ex^2^ = 6dx+38ex^2^
 [ -2  0  2 ]  ->  [   0   4/3  ] -> 32ex^2^
 [ -2  0  1 ]  ->  [ -1/3  1/4  ] -> -2dx - 8ex^2^ + 6ex^2^ = -2dx-ex^2^

Clearly there are situations (some contrived) when each of these approximations would be better but in general it is likely that using closer data will result in less error even if it might reduce the order of the approximation.

Arbitrary stencils

By now we have more or less described how to deal with arbitrary three point stencils, but here is a summary…

• Given a set of three points whose relative position to the cell i is given by x_i , calculate the coefficients A_i by solving the matrix equation

        1           1           1          A_1       0
       x_1         x_2         x_3    x    A_2   =   0
   1/2 x_1^2   1/2 x_2^2   1/2 x_3^2       A_3       1
  

2D Fixed Grid Stencils

In 2D the Laplacian operator becomes d²/dx² + d²/dy². In a cartesian fixed grid, the points (f_i-1,j f_i,j f_i+1,j) can be used to calculate the coefficients for the D²f/Dx² and the points (f_i,j-1 f_i,j f_i,j+1) for D²f/Dy². And our coefficients in each direction are {1,-2,1} and since y_i,j is used in both directions the complete list of coefficients is given by {A_i,j = -4, A_i-1,j = 1, A_i+1,j = 1, A_i,j-1 = 1, A_i,j+1 = 1}. Of course we could proceed as we did in the 1D case and create the 2D taylor expansion of each of the points. Recall that the 2D taylor expansion of f(x,y) is f_i + df/dx_i x + df/dy_i y + ½ d²f/dx²_i x² + ½ d²f/dx²_i y² + d²f/dxdy_i xy + 1/6 d³f/dx³ x³ + ½ d³f/dx²dy x²y + ½ d³f/dxdy² xy² + 1/6d³f/dy³ y³ + O(x⁴). Our matrix equation (out to second order terms) then becomes

   1    1    1    1    1          A_{i,i}               0    (f)
   0   -1    1    0    0          A_{i-1,j}             0    (df/dx)
   0    0    0   -1    1          A_{i+1,j}     =       0    (df/dy)
   0   1/2  1/2   0    0          A_{i,j-1}             1    (d^2^f/dx^2^)
   0    0    0   1/2  1/2         A_{i,j+1}             1    (d^2^f/dy^2^)
   0    0    0    0    0                                0    (d^2^f/dxdy)

Note the matrix is no longer square and we would normally need a sixth point to make a 6x6 matrix - however since the bottom row is all zeroes, any solution to the first five equations (the first five rows) will also satisfy the last equation in the sixth row, and our reduced matrix equation becomes:

   1    1    1    1    1          A_{i,i}               0    (f)
   0   -1    1    0    0          A_{i-1,j}             0    (df/dx)
   0    0    0   -1    1          A_{i+1,j}     =       0    (df/dy)
   0   1/2  1/2   0    0          A_{i,j-1}             1    (d^2^f/dx^2^)
   0    0    0   1/2  1/2         A_{i,j+1}             1    (d^2^f/dy^2^)

which has the solution [ -4 1 1 1 1 ] as we have already seen. Note since none of the points used were diagonally displaced the product of x*y was zero for each point. This prevented the introduction of any cross terms in the taylor expansions and removed the need for a sixth data point. The lack of a any cross terms also cancels any higher order errors involving cross terms leaving only the pure 1D error terms we have already seen ( 1/12 d^2f/dx2 DX^{2 + 1/12 d}2f/dy^{2 DY}2) giving a second order stencil with only five data points! We will see next that a generalized 2D stencil needs 6 points to be first order and 10 points to achieve second order accuracy.

2D Non-Fixed Grid Stencils

In nested grid type AMR, cells will often surrounding data points on many different amr levels - determining which data points to use in the discretization of the laplacian operator is therefore not a trivial task. We can however, proceed in a similar manner to the 1D case by expanding the function in a taylor series and collecting the various terms… If we expand the potential in a taylor series at the surrounding points

[ \phi(x_i,y_i)=\phi +\frac{\partial \phi}{\partial x} x_i +\frac{\partial \phi}{\partial y} y_i + \frac{1}{2} \frac{\partial^2 \phi}{\partial x^2} x_i^2 + \frac{1}{2} \frac{\partial^2 \phi}{\partial y^2} y_i^2 + \frac{\partial^2 \phi}{\partial x \partial y} x_i y_i + O(\Delta x^3) \]

And solve for the linear combination of surrounding cells that leaves the right combination of the derivatives

[ \left [ \begin{array}{cccccc} 1 & 1 & 1 & 1 & 1 & 1 \\ x_1 & x_2 & x_3 & x_4 & x_5 & x_6 \\ y_1 & y_2 & y_3 & y_4 & y_5 & y_6 \\ \frac{x_1^2}{2} & \frac{x_2^2}{2} & \frac{x_3^2}{2} & \frac{x_4^2}{2} & \frac{x_5^2}{2} & \frac{x_6^2}{2} \\ \frac{y_1^2}{2} & \frac{y_2^2}{2} & \frac{y_3^2}{2} & \frac{y_4^2}{2} & \frac{y_5^2}{2} & \frac{y_6^2}{2} \\ x_1y_1 & x_2y_2 & x_3y_3 & x_4y_4 & x_5y_5 & x_6y_6 \\ \end{array} \right ] \times \left [ \begin{array}{c} A_1 \\ A_2 \\ A_3 \\ A_4 \\ A_5 \\ A_6 \\ \end{array} \right ] = \left [ \begin{array}{c} 0 \\ 0 \\ 0 \\ 1 \\ 1 \\ 0 \\ \end{array} \right ] \mbox{then} \sum{A_i \phi_i} = \nabla^2\phi_i+O(\Delta_x^3) \]

Using higher level cells

When a neighbor cell is refined there is the possibility of using some of the neighbors children instead of the neighbor itself in the stencil. Some of these neighboring children cells will be physically closer, but we have seen how nice the 5pt fixed grid stencil is at efficiently reducing errors. While the closer data would allow for a more accurate expression, it is likely not worth the added cost of potentially breaking the left-right symmetry and introducing diagonal stencil components as we will need at least 10 stencil points to achieve better accuracy instead of just 5.

Using lower level cells

While we never have to use higher level cells in the stencil (since the higher level cells will have a neighboring parent cell) we will have to use lower level cells at the edges of grid interfaces. These parent cell centers will be off-diagonal and will at most have an offset of (+-3/2, +-3/2). Of course with improper nesting grandparent cells could be at (+-5/2, +-5/2) and great grandparent cells at (+-9/2, +-9/2). Since cells next to parents etc… are likely to not have large gradients in the flow variables (ie density), the gravitational potential (phi) will likely have small third and higher derivatives. However, the second derivatives could still be large - including the cross terms - so it is important to construct a stencil that is at least 1st order accurate Since we need 6pts to achieve this, and there are 5 obvious choices (the cells to the left, right, top, bottom, and the cell itself) we need to select a sixth (and perhaps 7th etc…) cell to use in our stencil.

One possible approach

• If the left cell is a parent at (-3/2,-½) then use the parent at (-3/2,3/2) to balance the diagonal term. If that parent is absent use the grandparent etc…
• If the left cell is a grandparent at (-5/2,-3/2) then use the grandparent at (-5/2,5/2) to balance. If that grandparent is absent use the greatgrandparent etc…
• and so on…

One other approach

Of course adding more cells to the stencil will only improve the accuracy and since cells are fined in 2x2 regions, there are always two adjacent neighbors on the same level and one other neighbor on the diagonal that could also be used to balance. This would be algorithmically simpler to implement though it is not clear whether this would lead to smaller errors. This approach would be summarized as follows:

• If any of the adjacent cells are ancestors, link to them but also link to the diagonal neighbor within the same parent cell.

Of course coefficients would have to be calculated for the range of possibilities, but since only two adjacent edges can be ancestors, there would be a small number of scenarios to calculate coefficients for (MaxLevel²)

3D AMR

In 3D the taylor expansion has 1+3+6 = 10 terms second order and below and so in general a 10 pt stencil will be needed. (Of those 10 terms, 3 are cross terms so a fixed grid stencil would only need the expected 7 points.) An extension to the second approach above would be to include the 3 diagonal neighbors that share an edge to balance any parents/grandparents. (Each amr parent cell would have 8 child cells so any child cell will have 7 surrounding neighbors - three that share a face, three that share an edge, and 1 that shares a corner. ) Since a cell has 6 adjacent neighbors (on any level), only three additional cells are needed to cancel any cross terms and the three that share an edge seem to be the best bet. And there are only MaxLevel³ sets of coefficients to calculate.