I used pretty much all my expertise with finite differences to solve an advection-diffusion equation with space-dependent coefficient with a grid in the $x$,$z$ domain with regular spacings. Something like this:

enter image description here

This is nice because it is easy to initialize with some experimental measurements that I have. However, near the sloping boundary, I have to carefully choose the grid spacings or sometimes my code behaves funny. This also limits the vertical resolution that I can have very close to the boundary.

To solve this problem I was thinking to change reference system by rotating it by an angle $\theta$. Something like this:

enter image description here

so that I could easily increase the resolution near the boundary. What I do not know is if this imposes some limitations on the physics. Is it sufficient to multiply my equations by a rotation matrix in order to be consistent with the new reference system?

Alternatively, is it possible (and fairly easy to do) to add triangular cells near the slope and still use finite difference approximations?

  • $\begingroup$ is "finite-element" tag there by mistake? $\endgroup$ – Anton Menshov Mar 22 '18 at 0:07
  • $\begingroup$ @AntonMenshov more likely due to my ignorance. I thought that having rectangular and triangular cells in the same domain would imply passing to finite element techniques. Is it so? $\endgroup$ – shamalaia Mar 22 '18 at 0:36
  • $\begingroup$ I am not a "finite-difference evangelist", so I would be biased towards saying if you have a complex boundary, you should look towards the finite-element method. However, your boundary is not that complex and you might have reasons to wish staying limited to the FD-realm. $\endgroup$ – Anton Menshov Mar 22 '18 at 1:09
  • $\begingroup$ A couple of notes: (i)At the moment you appear to only show one boundary in the second figure. If your other boundaries are grid aligned, you do realise you're solving on a different domain, right? (ii) in drawing boxes/cells rather than points, you're implying you might actually be dealing in finite volumes, rather than finite differences, is this true? $\endgroup$ – origimbo Mar 22 '18 at 12:02
  • $\begingroup$ How are you choosing the finite difference equations at the sloped boundary? $\endgroup$ – H. Rittich Mar 22 '18 at 12:57

In case you can easily implement the roated grid, this is possibly the easist fix, so let me answer your question about the invariance of the convection diffusion equation under rotation of the coordinate system.


Your physical problem does not depend on the coordinate system that you choose. Only the definition of the Laplace operator and the gradient require that you choose a cartesian coordinate system. Thus, rotating the coordinate system will change nothing, as long as you rotate the velocities apropriately.


Consider the convection diffusion equation $$ \underbrace{\sum_i \frac{\partial^2 u}{\partial x_i^2}}_{\text{diffusion}} + \underbrace{\sum_i \frac{\partial (v_i u)}{ \partial x_i}}_{\text{convection}} = f \,. $$ A rotation is a linear transformation. Hence, assume the new coordinates $x_i'$ are given by $$ x_i = \sum_j a_{ij} x_j', $$ which is the multiplication by a matrix $A$. Furthermore, assume that we can obtain the new coordinates given the old ones by $$x'_i = \sum_j b_{ij} x_j .$$ The numbers $b_{ij}$ are the entries of the matrix $B := A^{-1}$.

We denote the function $u$ in the new coordinates also by $u$: $$ u(x') = u(x(x'))\,. $$

Computation of the Diffusion Term

To compute the diffusion term, we need to use the chain rule. We start with the first partial derivative. We have that $$ \frac{\partial u}{\partial x_i} = \sum_k \frac{\partial u}{\partial x_k'} \frac{\partial x_k'}{\partial x_i} = \sum_k \frac{\partial u}{\partial x_k'} b_{ki} \,. $$ We use this equation to compute the second partial derivative. We obtain that \begin{align*} \frac{\partial}{\partial x_i} \left(\frac{\partial u}{\partial x_i} \right) &= \frac{\partial}{\partial x_i} \left( \sum_k \frac{\partial u}{\partial x_k'} b_{ki} \right) \\ &= \sum_k \frac{\partial}{\partial x_i} \left(\frac{\partial u}{\partial x_k'} b_{ki} \right) \\ &= \sum_k \sum_\ell \frac{\partial}{\partial x_\ell'} \left( \frac{\partial u}{\partial x_k'} b_{ki} \right) \frac{\partial x_\ell'}{\partial x_i} \\ &= \sum_k \sum_\ell \frac{\partial^2 u}{\partial x_k' \partial x_\ell'} b_{ki} b_{li} \,. \end{align*} The Laplace operator is the sum over all second partial derivatives. Hence, we compute \begin{align*} \sum_i \sum_k \sum_\ell \frac{\partial^2 u}{\partial x_k' \partial x_\ell'} b_{ki} b_{\ell i} &= \sum_k \sum_\ell \frac{\partial^2 u}{\partial x_k' \partial x_\ell'} \sum_i b_{ki} b_{\ell i} \\ &= \sum_k \sum_\ell \frac{\partial^2 u}{\partial x_k' \partial x_\ell'} c_{k\ell} \,, \end{align*} where $$ c_{kl} := \sum_i b_{ki} b_{\ell i} \,. $$

Computation of the Convection Term

To compute the convection term, we have to use the chain rule, again. Using this rule, \begin{align*} \frac{\partial (v_i u)}{\partial x_i} &= \frac{\partial (v_i u)}{\partial x_i} \\ &= \sum_k \frac{\partial (v_i u)}{\partial x_k'} \frac{\partial x_k'}{\partial x_i} \\ &= \sum_k \frac{\partial (v_i u)}{\partial x_k'} b_{ki} \\ \,. \end{align*} We obtain the gradient by summing over all partial derivatives. Thus, $$ \sum_i \frac{\partial (v_i u)}{\partial x_i} = \sum_i \sum_k \frac{\partial (v_i u)}{\partial x_k'} b_{ki} = \sum_k \frac{\partial (v_k' u)}{\partial x_k'} \,, $$ where $$ v_k' = \sum_i b_{ki} v_i \,. $$ The numbers $v_k'$ are just the velocities in the new coordinate system.

Rewriting the PDE

We can now give the PDE in the new coordinate system. The PDE is $$ \sum_k \sum_\ell \frac{\partial^2 u}{\partial x_k' \partial x_\ell'} c_{k\ell} + \sum_k \frac{\partial (v_k' u)}{\partial x_k'} b_{ki} = f \,. $$ However, since we have assumed that $A$ is a rotation matrix, we can simplify this equation.


A rotation is a special linear transformation; the matrix $A$ is orthogonal, i.e., the inverse $B$ of $A$ is the matrix $A$ transposed ($b_{ij} = a_{ji}$). Thus, $$ \delta_{ij} = \sum_{\ell} b_{i\ell} a_{\ell j} = \sum_{\ell} b_{i \ell} b_{j \ell} = c_{ij} $$

Hence, the diffusion term simplifies to $$ \sum_k \sum_\ell \frac{\partial^2 u}{\partial x_k' \partial x_\ell'} c_{k\ell} = \sum_k \frac{\partial^2 u}{\partial x_k' \partial x_k'} \,. $$ And the PDE becomes $$ \sum_k \frac{\partial^2 u}{\partial {x_k'}^2} + \sum_k \frac{\partial (v_k' u)}{\partial x_k'} = f \,, $$ which is just the convection diffusion equation with rotated velocities.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.