A is not symmetric matrix
piyueh opened this issue · 14 comments
I thought A should be always symmetric after we use scaling matrices, but it turns out that, when we use non-uniform grids, it is not. Is there a bug? or just my understanding is wrong?
How to reproduce :
(I used PETSc 3.6 + branch AmgXSolver and PETSc 3.5 + branch develop, both cases produced the same results)
mpiexec -np 1 ${PETIBM_BINARY}/petibm2d -outputs -directory ${PETIBM_EXAMPLES}/2d/lidDrivenCavity/Re100NonUniform
Open the file ${PETIBM_EXAMPLES}/2d/lidDrivenCavity/Re100NonUniform/outputs/A.output, and then check the extries A(0, 31) and A(31, 0). They should have the same values, but actually they dont'.
Update: even in uniform-grid case, it is still not symmetric. Weird values can be observed in the rows representing velocity points on top and bottom boundaries. The locations of weird values correspond to the coefficients related to the next stencil points in the direction perpendicular to the boundaries.
Can you dump the matrix for a 3x3 uniform grid in a domain of size 3x3? I'll also look for bugs in the code.
This is the A matrix using 2D lid driven case with 3x3 grid in 3x3 domain:
Mat Object: 1 MPI processes
type: seqaij
row 0: (0, 50.03) (1, -0.005) (2, -0.00666667)
row 1: (0, -0.005) (1, 50.03) (3, -0.00666667)
row 2: (0, -0.005) (2, 50.02) (3, -0.005) (4, -0.005)
row 3: (1, -0.005) (2, -0.005) (3, 50.02) (5, -0.005)
row 4: (2, -0.00666667) (4, 50.03) (5, -0.005)
row 5: (3, -0.00666667) (4, -0.005) (5, 50.03)
row 6: (6, 50.03) (7, -0.00666667) (9, -0.005)
row 7: (6, -0.005) (7, 50.02) (8, -0.005) (10, -0.005)
row 8: (7, -0.00666667) (8, 50.03) (11, -0.005)
row 9: (6, -0.005) (9, 50.03) (10, -0.00666667)
row 10: (7, -0.005) (9, -0.005) (10, 50.02) (11, -0.005)
row 11: (8, -0.005) (10, -0.00666667) (11, 50.03)
I think line # 142-143, 164-165 and their corresponding parts in 3D version in generateA() should be modified in order to solve this problem. Especially the 0.5*dx[i] and 0.5*dy[i], which I guess should be the distances between the boundary points to ghost points (?). Also, I guess these values should be connected to the BC vector for velocity in line 85-86, 128-129 and corresponding part in 3D version in generateBC1(). Probably in the function updating ghost cells, some codes may also need modifications. Because all these functions are related to ghost points.
PetIBM and cuIBM both do not use ghost cells. They use the values of velocity directly from the boundary. So the stencils near the boundary are all slightly modified and do not resemble the stencils in the interior. If you want exact details about how this is done, you can refer to my thesis. That is why you see 0.5*dx instead of dx.
I am yet to figure out if this will in fact result in a non-symmetric matrix. I'll look at that in more detail later today.
For the Taylor-Green vortex problem, the domain is periodic is both x- and y-directions, and we end-up with a symmetric implicit matrix.
For the lid-driven cavity flow, the implicit matrix is not symmetric.
The implicit matrix A is not symmetric if the boundary conditions are not periodic and we do not use ghost-cells. The problem arises when computing the Laplacian coefficients.
Let's take the example of a uniform mesh with grid-spacing h, and we set a Dirichlet boundary condition on the bottom side of the domain.
We want to evaluate the second-order derivative of the u-velocity in the y-direction: d2udy2 at a bottom u-node (still in the fluid region, j=0). We get:
d2udy2 = a1*uSouth + b1*uPoint + c1*uNorth
where a1=2/dyMinus/(dyMinus+dyPlus) and c1=2/dyPlus/(dyMinus+dyPlus)
with dyMinus=0.5*h and dyPlus=h.
Now, let's compute the second-order derivative on the next u-node located on the north.
d2udy2 = a2*uSouth + b2*uPoint + c2*uNorth
where a2=2/dyMinus/(dyMinus+dyPlus) and c2=2/dyPlus/(dyMinus+dyPlus)
with dyMinus=h and dyPlus=h.
Because c1 != a2, the implicit matrix A is not symmetric.
So does this mean we make the matrix A nonsymmetric in purpose in our original design? If yes, I wonder the advantages of this. If no, then @mesnardo 's explanation explains why I think we may want to modify the part 0.5*dx[i] and 0.5*dy[i]. At least for uniform grid we can make A symmetric.
If we make A nonsymmetric in purpose, how can we assure that CG works all the time? We use CG in our all examples. I know there are some chances that CG works on nonsymmetric matrices and even non positive definite matrices.
The velocity value of the ghost-cell point can be extrapolated from the the velocity value at the physical boundary and the velocity value at the first point in the interior of the domain.
To make the implicit matrix symmetric, even for non-uniform grids, the distance between the ghost-cell point and the first interior point should be equal to the distance between the first interior point and the second interior point.
@piyueh We are not mathematically assured of it, but from a practical perspective I think CG will (and has) always worked on that system. It converges quickly and has never given us problems. I would recommend leaving it the way it is.
@gforsyth #53 is related to the Poisson system. That matrix is always going to be symmetric. But I believe the issue was that GAMG sometimes produces non-symmetric preconditioners. Which means it's probably best that that system is solved using BiCGStab by default.
If A is not symmetric or grid is non-uniform, and we use high-order approximation of inverse A as B (which is not implemented in our codes), will QTBNQ still be a symmetric matrix? Because I think is such case, B is not symmetric anymore (because the Laplacian operator may not be symmetric).
@piyueh You're right. But again, practically it may make no difference since most of the matrix is symmetric. In fact, the higher order approximations have been observed to converge with CG in the past. We have also noticed that there are no serious gains from using higher order approximations of QTBNQ. They take up too much memory to justify their use, especially when preconditioners are included. Moreover, if GAMG (which might be our PETSc preconditioner of choice) can produce asymmetric matrices anyway, then strict compatibility with CG criteria becomes even less important.
Taira & Colonius (2007) claim that there are ways to make the A matrix symmetric, and they use symmetric A matrix in their code. The real issue of this is that we claim to re-implement their method but are not able to create symmetric A matrix as what they do. This is why I felt confused at that time.
We are still not able to create a symmetric A matrix now, though this may not be important because the solvers still run well with unsymmetric A matrix.
I'm not sure if we should close this issue. We haven't understood why they can have a symmetric A matrix and how to do that. I feel this is still an open issue unless we all agree to ignore this.
But I want to remove the bug label first. It's not a bug in our code.