# Symmetrization of Laplacian Matrix Operator (finite volumes)

The aim is to construct symmetric Laplacian matrix $$L$$ for 2D square domain.

I have the following discretization of second derivatives on non-uniform grid (will skip the steps of derivation, any textbook works). velocity vectors $$u$$ and $$v$$ are defined on cell faces, $$x_i$$ corresponds to $$u_{i+mk}$$, where m is amount of unknown $$u$$'s and k is an integer.

$$\left(\frac{\partial^2 u}{\partial x^2}\right)_i \approx \frac{u_{i+1}\left(x_i-x_{i-1}\right)-u_i\left(x_{i+1}-x_{i-1}\right)+u_{i-1}\left(x_{i+1}-x_i\right)}{\left(\frac{x_{i+1}-x_{i-1}}{2}\right)\left(x_i-x_{i-1}\right)\left(x_{i+1}-x_i\right)}+O(\Delta x).$$ This scheme becomes second order when $$\Delta x_i\approx\Delta x_{i-1}$$. The above discretization leads to matrix $$\hat{L}$$.

Now, consider the symmetrizer matrix $$\hat{M} =\left[\begin{array}{cc} \frac{1}{2}\left(x_{i+1}+x_{i-1}\right) & 0 \\ 0 & \frac{1}{2}\left(y_{j+1}+ y_{j-1}\right) \end{array}\right],$$ which is a diagonal.

The product $$\hat{M}\hat{L}=L$$ makes L symmetric for Dirichlet BC and Neumann (1st order) BC.

Now, if I want to use the second order scheme for Neumann BC, discretization becomes 3 point. Therefore, it starts affecting the off diagonal elements in $$\hat{L}$$ (first oder boundary condition scheme changes only the diagonal element, which does not kill the symmetricity). Hence, $$L$$ is not symmetric anymore.

Is there a way to modify the diagonal symmetrizer $$\hat{M}$$ without inducing a lot of computational time? $$\hat{M}$$ is multiplied by the rhs vector at every time step. I tried the simple $$\hat{M}=(\frac{\hat{L}+\hat{L}'}{2})\hat{L}^{-1}$$ but it is too dense (see pic below).

Edit: I figured it out. Will post the solution when finish typesetting it in LaTeX.

I finally figured it out.

We need to move to the new variables.

(Idk why the typeset result formulas look different than the preview.)

Let us have $$m$$ intervals in $$x$$ direction and $$n$$ intervals in $$y$$ direction. The number of unknowns is equal to the number of face velocities in the interior of the domain excluding the boundaries. Rewrite $$\vec{u}$$ as a vector of unknowns

$$\vec{u} = [\underbrace{u_1;u_2;\dotsc;u_{m-1};u_1;u_2;\dotsc;u_{m-1};\dotsc;u_1;u_2;\dotsc;u_{m-1}}_{n\text{ times for each }1...m-1\text{ block}}, \underbrace{v_1;v_2;\dotsc;v_{m};v_1;v_2;\dotsc;v_{m};\dotsc;v_1;v_2;\dotsc;v_{m}}_{n-1 \text{ of each }1,\dotsc,m\text{ block}}].$$

To maintain consistency - indexation goes left to right from bottom to top.

Define matrices

$$\Delta y_j=diag([\underbrace{\Delta y_1; \Delta y_1; \dotsc;\Delta y_1}_{m-1\text{ times}}; \underbrace{\Delta y_2; \Delta y_2; \dotsc;\Delta y_2}_{m-1\text{ times}};\dotsc;\underbrace{\Delta y_n; \Delta y_n; \dotsc;\Delta y_n}_{m-1\text{ times}}]),$$

$$\Delta x_i=diag([\underbrace{\Delta x_1; \Delta x_2;\dotsc\Delta x_m;\Delta x_1; \Delta x_2;\dotsc\Delta x_m; \dotsc;\Delta x_1; \Delta x_2;\dotsc;\Delta x_m}_{n\text{ times}}])$$

$$$$R =\left[\begin{array}{cc} \Delta y_j & 0 \\ 0 & \Delta x_i \end{array}\right],$$$$

Let the new vector variable of unknowns (something like mass flux) be defined as

$$\vec{q} = R\vec{u},$$

where $$R$$ is a diagonal matrix of size $$(m-1)n\times m(n-1)$$ defined above.

Next, we need to define the difference matrices to cancel out the $$\frac{x_{i+1}−x_{i-1}}{2}$$ central term in diffusion discretization.

Let

$${\Delta x_i+\Delta x_{i-1}} = diag[\underbrace{\Delta x_2+\Delta x_1; \Delta x_3+\Delta x_2; \dotsc \Delta x_m-\Delta x_{m-1};...;\Delta x_2+\Delta x_1; \Delta x_3+\Delta x_2; \dotsc \Delta x_m-\Delta x_{m-1}}_{n\text{ times for each block of } \Delta x_2+\Delta x_1\text{ to }\Delta x_m-\Delta x_{m-1}}]$$

$${\Delta y_j+\Delta y_{j-1}} = diag[\underbrace{\Delta y_2+\Delta y_1; \dotsc \Delta y_2+\Delta y_1;}_{m\text{ times}} \underbrace{\Delta y_3+\Delta y_2; \dotsc \Delta y_3+\Delta y_2;}_{m\text{ times}} \underbrace{\Delta y_n+\Delta y_{n-1}\dotsc \Delta y_n+\Delta y_{n-1}}_{m\text{ times}}],$$

$$$$\hat{M} =\left[\begin{array}{cc} \frac{1}{2}\left(\Delta x_i+\Delta x_{i-1}\right) & 0 \\ 0 & \frac{1}{2}\left(\Delta y_j+\Delta y_{j-1}\right) \end{array}\right].$$$$

Now, since $$\vec{q}=R\vec{u}$$ then $$\vec{u}=R^{-1}\vec{q}$$, Going back to Laplacian matrix

$$L\vec{u}= LR^{-1}\vec{q}.$$

Premultiply by $$\hat M$$, the expression will become

$$\hat MLR^{-1}\vec{q},$$

where $$\hat M L R^{-1}$$ is symmetric by construction.