# How far is a non-symmetric discretization of an elliptic operator from the continuous operator itself?

I am investigating the accuracy and stability properties of a non-symmetric discretization of a Poisson problem. The problem originates from a ghost fluid discretization of the projection step of a fractional step Navier-Stokes solver. In this type of approach, grid points away from the geometry surface are discretized using standard 5-point stencils. For points adjacent to the geometry, a boundary condition preserving extrapolation is used to fill in a set of ghost points right outside the geometry. These points are then used to fill in the 5 point stencil of the points immediately interior to the geometry surface. The act of extrapolating to the 'ghost points' introduces a non-symmetric character to the matrix, and I'm pretty much trying to figure out how much doing this ruins the matrix properties. One particular example of such an approach is shown in the paper by Mittal et al.

Mittal R, Dong H, Bozkurttas M, Najjar FM, Vargas A, von Loebbecke A., A VERSATILE SHARP INTERFACE IMMERSED BOUNDARY METHOD FOR INCOMPRESSIBLE FLOWS WITH COMPLEX BOUNDARIES. J Comput Phys. 2008;227(10):4825-4852.

• Could you clarify what you mean "non-symmetric discretization"? Do you have specific discretization schemes in mind?
– Paul
Jul 24 '12 at 17:39
• Can you elaborate an example, including the choice of spaces and the motivation to use a "non-symmetric discretization"? It is absolutely not clear what you are talking about. Jul 24 '12 at 18:15
• @JohnMousel Welcome to scicomp! I recommend that you edit your question directly rather than putting essential information in a comment. Also, you're more likely to get help if you write down an example of the discretizations you're considering, or at least provide a link to details. Jul 24 '12 at 19:25
• Edited. Hopefully it's more clear now. Jul 24 '12 at 20:02
• Yes, it's better. You say " I'm pretty much trying to figure out how much doing this ruins the matrix properties." Which matrix properties? Just symmetry? Jul 25 '12 at 5:20

There is a field of research that can be identified by the keywords "mimetic discretizations" and "compatible spatial discretizations" that starts from your same assumption: algebraic equations arising from the discretization of partial differential equations should mimic the properties of the continuum operators.

A good starting point could be Principles of Mimetic Discretizations of Differential Operators by Bochev and Hyman, which appeared in this book.

It is an interesting framework, but I would not focus to much on "matrix properties" per se: finally what matters is obtaining a "good solution" (robust, accurate, not to much expensive, ...) and not having a nice matrix at hand.

I would suggest you to start asking more narrow questions: e.g. how can I deal efficiently with the unsymmetric terms? Can I trade accuracy with speed by using a symmetric solver?

Jed Brown correctly points out that sometimes lack of symmetry is only apparent.

Let me explain this point with a very trivial example. Suppose that $\mathcal{L}u=f$ + b.c (continuous, with $\mathcal{L}$ self adjoint) becomes $Ax=b$ (discrete, with $A=A^T$) where the b.c. are expressed as $u = Bv$ (see this answer for a concrete example, sorry for the shameless self citation). Now by substitution we have $ABv = b$ where $AB$ is not symmetric, but symmetry can be restored by premultiplication by $B^T$ to obtain $B^TABv = B^Tf$. If matrix $B$ is square (and note that usually it is an identity matrix plus some low rank correction), problem $ABv=f$ is perfectly equivalent to $B^TABv = B^Tf$. If premultiplication by $B^T$ is not feasible, sometimes a symmetric problem can be obtained by augmenting the system, i.e. by introducing extra unknown variables.

[disclaimer: what follows is only a conjecture since I had no time to work out your problem in detail]

In the case at hand b.c. are introduced by adding ghost points outside the problem domain: unknown fields at the ghost points are obtained by extrapolation from the domain points, taking into account boundary geometry and b.c. (The algorithm actually goes the other way: from the ghost point an image point in the problem domain is derived, and its value is obtained by interpolation of the surrounding problem points. For this discussion however the approaches are equivalent.) The source of "unsymmetry" is clear, and the situation is slightly different from my trivial example above, so that there is no easy recipe for fixing the symmetry.

My suggestion is to investigate in these two directions:

1. Is this a sensible way of introducing the boundary conditions? If the answer is yes, lack of symmetry is acceptable at the formulation level.

2. Is lack of symmetry acceptable at the computational level? I.e. are you able to solve the unsymmetric problem? If no, try to modify the system for restoring symmetry in a way that does not introduces further approximations (like the trivial example above) or trade accuracy with speed by introducing some extra approx.

• I think this might be a lead into what I was looking for. I'm using GAMG in PETSc as a preconditioner for this system which symmetrizes the input graph. I agree that the performance is what matters, I just didn't want to put something out there without fully understanding it! Aug 2 '12 at 14:20
• @JohnMousel It's common for nonsymmetric discretizations to still have a symmetric graph. For a self-adjoint operator, lack of discrete symmetry usually comes from boundary conditions (often "fixable" without deeply changing the method) and nonuniform grid spacing or variable coefficients (usually requiring deeper changes, with tradeoffs, to make symmetric). Aug 2 '12 at 15:30
• @JedBrown Jed, can you explain what you mean by 'fixable'? Aug 2 '12 at 19:13
• @JedBrown Jed, thanks for your comments. I've taken the liberty of trying to incorporate your point in my reply. Aug 3 '12 at 8:16
• @StefanoM Thanks for the lucid edit. I never thought of it from the point of view of factoring it into a product of matrices. Currently, to make the system solvable, I've been explicitly enforcing a solvability constraint on the rhs by projecting it out of the null space of Atranspose. If I symmetrize like you suggest, then this won't be necessary anymore. Aug 3 '12 at 13:51