I want to verify and compare different Discretizations of the anisotropic diffusion equation in 2D / 3D.

enter image description here

In order to both test the timestepping and the spatial discretisations I had a look at using the heat kernels as an analytical solution to diffusion equations. The one dimensional heat kernel (wiki-link) looks like this: $$ u(x,t) = \frac{1}{(4\pi t)^{d/2}} e^{-x^2/4t}\ $$ and solves following diffusion equation of dimension $d$ analytically: $$ \frac{\partial u}{\partial t} = -\Delta u $$

The idea is that for an analytical solution I can verify my convergence both spatially and temporal. It is important in my case, that the diffusion is anisotropic, and might not be aligned with the underlying grids.

My anisotropic diffusion equation reads:

$$ \frac{\partial u}{\partial t} = \nabla (D \nabla u) $$

In order to expand the analytical solution to my case, where the diffusion is not scalar and unity, but a full (symmetric positive definite) tensor D, I had a look at multivariate gaußian distributions (wiki-link).

My question is the following:

  1. Is there somewhere in the literature a heat kernel for anisotropic (but homogeneous) diffusion which I can use to verify my numerics?

By analogy I would come up with something like:

$$ u(\vec{x},\vec{\mu},t) = \frac{1}{(4\pi~ |\mathbf{D}| t)^{d/2}} \exp \large(\frac{-(\vec{x}-\vec{\mu})^{T} \mathbf{D}^{-1}(\vec{x}-\vec{\mu})}{4~|\mathbf{D}|~t}\large). $$

  1. Is the above equation a proper analytical solution of the anisotropic diffusion equation?

I am not sure if the exponent $d/2$ is correct for the $|D|$. I am aware that the question might touch both physics, and numerics, but I thought this might be useful to other people trying to verify their codes.

Any hints or help is appreciated!

  • 1
    $\begingroup$ I'd start by looking for a coordinate transformation (using the eigenvalues and eigenvectors of $D$) which transforms your anisotropic diffusion into isotropic diffusion. Then you can transform the classical heat kernel solution back into the original anisotropic space. But I haven't worked out the details. $\endgroup$ – user3883 Jul 20 '18 at 10:46
  • $\begingroup$ This may be a stupid comment, but have you thought about just manufacturing an analytical solution using the method of manufactured solutions? This is generally what I do for anisotropic equations... $\endgroup$ – BlaB Aug 22 '18 at 14:07

Since Laplacian is an elliptic operator you are looking for the Cholesky decomposition of the assumed constant diffusion matrix $D$: $$D = LL^T$$

Therefore the parabolic equation may be written as: $$\partial_t u=-\nabla^T(D\nabla u)=-(L^T\nabla)^T(L^T\nabla u)=-\tilde{\nabla}^T(\tilde{\nabla} u)=-\tilde{\nabla}^2u$$

Naming $\tilde{\nabla}=L^T\nabla$ and therefore $\tilde{x}=L^{-1}x$. This last PDE is totally decoupled (using Fourier transform and separation of variables).

The solution for the unnormalised heat kernel will be (using $\tilde{x}=L^{-1}x$): $$\Phi=\prod_{n=1}^d\frac{\exp{(-\tilde{x}_n^2/(4t))}}{\sqrt{4\pi t}}=\frac{\exp{(-\tilde{x}^T\tilde{x}/(4t))}}{(4\pi t)^{d/2}}=\frac{\exp{(-x^TD^{-1}x/(4t))}}{(4\pi t)^{d/2}}$$

| cite | improve this answer | |
  • $\begingroup$ There is one step missing here in getting kernel in $x$ space. The solution involves a convolution so you have to transform the $d$-dimensional measure $d \tilde{x}$ to $dx$. In $x$-space the kernel becomes $\exp(-x^\top D^{-1}x/(4t))/(4\pi t |D|)^{d/2}$. $\endgroup$ – cfdlab Aug 22 '18 at 11:00

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.