1d diffusion equation

Integrating the diffusion equation,

$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2},$$

with a constant diffusion coefficient D using forward Euler for time and the finite difference approximation for space,

$$u_i^{t+1} = u_i^t + D\frac{\Delta t}{\Delta x^2} ( u_{i+1}^{t} + u_{i-1}^{t} - 2u_i^{t} ),$$

leads to the conservation of $\bar{u}=\sum_i \Delta x\, u_i$ over time (see animation 1 and figure 1), because reflective Neumann boundary conditions,$\partial_x u=0$, are employed at the borders (forward differences: $u_i = u_{i \pm 1}$):

$$u_i^{t+1} = u_i^t + D\frac{\Delta t}{\Delta x^2} \cdot ( u_{i\pm1}^{t} - u_{i}^{t} )$$

Space and time are discretized with $\Delta x=0.01$ and $\Delta t=10^{-7}$. The diffusion coefficient is $D=10$.

In 2d polar coordinates $(r,\phi)$, the Laplacian is given by:

$$\nabla^2 u = \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2} \frac{\partial u^2}{\partial \phi^2}.$$

In case of an axisymmetric distribution $u(r,\phi)=u(r)$, the Laplacian reduces to:

$$\nabla^2 u = \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r}$$

Forward Euler and finite difference approximation with central differences for the advective term leads to (derivation):

$$u_i^{t+1} = u_i^t + D\frac{\Delta t}{\Delta r^2} \big( u_{i+1}^{t} + u_{i-1}^{t} - 2u_i^{t} \big) + D\frac{\Delta t}{2 \Delta r} \big( u_{i+1}^{t} - u_{i-1}^{t} \big) \\ = u_i^t + D\frac{\Delta t}{\Delta r^2} \big( (1+0.5/i)u_{i+1}^{t} + (1-0.5/i)u_{i-1}^{t} -2u_i^{t} \big).$$

The same approximation is obtained with the Finite Volume Method (derivation). The boundary condition at the origin exploits the rotational symmetry and thus is $\partial_r u = 0$, leading to (central differences: $u_{i-1}=u_{i+1}$)

$$u_i^{t+1} = u_i^t + D\frac{\Delta t}{\Delta r^2} \cdot 4( u_{i+1}^{t} - u_i^{t} )$$

At the other boundary, Neumann conditions $\partial_r u = 0$ are employed and realized as (central differences: $u_{i-1}=u_{i+1}$):

$$u_i^{t+1} = u_i^t + D\frac{\Delta t}{\Delta r^2} \cdot 2( u_{i-1}^{t} - u_{i}^{t} )$$

However, there is a problem: Conservation of u is violated (see animation 2 and figure 2).

Over time the total amount of u grows as calculated via

$$2\pi \int_\Omega u(r) r dr = 2\pi \sum_{i=0}^{n} u_i\; (i \,\Delta r) \; \cdot \; \Delta r,$$

where $(i \Delta r)$ is the discretized volume element from 2d polar coordinates. The dashed line gives the analytical result and the regular line the numerical result. Numerical parameters $\Delta r=\Delta x$, $\Delta t$ and $D$ are as before.

Question

Why is conservation of u violated when going from the 1d system to the 1d radial (reduced from 2d polar coordinates)? What can I do to regain the conservation?

Finite volume methods lead to the same discretization scheme as shown above. Robin boundaries (as used for conservation in advection equations) are not applicable, since the flux $j$ at the boundaries is given by $j=r\partial_r u$, which leads to the already employed Neumann boundary conditions $\partial_r u=0$ at $r=0$ and $r=r_{end}$.

Edit:

While testing a few things, I came up with an MWE in c++ for others interested to try. Compilation instructions are given at the top. It reproduces the conservation violation problem :(

// test program to check conservation in radial diffusion

#include <fstream>                              // std::ofstream
#include <string>                               // std::string
#include <cmath>                                // exp

int main(){

// initialize variables
std::string name = "/tmp/u_sum.dat";
int nx = 1000;
size_t nsteps=100000;
double d = 10;
double dr = 0.01;
double dt = 1e-7;
double diffcoeff = d*dt/(dr*dr);
double *u = new double[nx];
double *unew = new double[nx];
double *usum = new double[nsteps]();
double usum0;

// initial condition
float mu_x = 0*nx;  // mean x
float xsigma = 0.01*nx; // variance x
for(int x=0; x<nx; x++) u[x] = exp(-0.25*((x-mu_x)*(x-mu_x)/(xsigma*xsigma)));

// time evolution
for(size_t step=0; step<nsteps; step++){

#pragma omp parallel for
for(int x=0; x<nx; x++){

int left=x-1;
int right=x+1;

// calculation, central difference with lHospital at r->0 with du/dr -> 0
if(x==0){
unew[x] = u[x] + diffcoeff*( 4*(u[right] - u[x]) );
}else if(x==nx-1){
unew[x] = u[x] + diffcoeff*( 2*(u[left] - u[x]) );
}else{
unew[x] = u[x] + diffcoeff*( u[right]*(1+0.5/x) + u[left]*(1-0.5/x) - 2*u[x] );
}
}

// sum up to check conservation
usum0 = 0;
for(int x=0; x<nx; x++) usum0 += 2*M_PI*dr*dr*x*unew[x];
if(!(step%10000)) printf("%12.12f\n",usum0);

for(int x=0; x<nx; x++) usum[step] += 2*M_PI*dr*dr*x*unew[x];

// update u,unew
std::swap(u,unew);
}

// save stuff
FILE *fp;
if((fp=fopen(name.c_str(), "w"))==NULL){ printf("Cannot open file.\n"); }
for(size_t step=0; step<nsteps; step++) fprintf(fp,"%12.12f\n",usum[step]);
fclose(fp);

// free memory
delete[] u;
delete[] unew;
delete[] usum;

}

Update

In the paper "High-order schemes for cylindrical/spherical geometries with cylindrical/spherical symmetry" by Wang et al. (2013), the authors propose the following scheme:

$$\frac{d u_i}{dt} = \frac{1}{\Delta V} (r_{i+\frac{1}{2}}j_{i+\frac{1}{2}} - r_{i-\frac{1}{2}}j_{i-\frac{1}{2}}), \\ \Delta V = \frac{1}{2}(r^2_{i+\frac{1}{2}} - r^2_{i-\frac{1}{2}}), \\ j_i = \nabla u_i = \frac{1}{2\Delta r}(u_{i+1} - u_{i-1}).$$ The half index, $i+\frac{1}{2}$, implies an arithmetic mean. Is this scheme applicable here? Is the gradient discretized correctly?

• When you say "conservation of $u$" in the radial case, do you think of $\int_\Omega u(r,t) \; dr$, or $\int_\Omega u(r,t) \; 2\pi r \; dr$? Jan 24 '18 at 20:26
• For the radial case, I think of the second integral with the volume element. For the 1d case on a line, your first integral is appropriate. Jan 24 '18 at 20:50
• Are you integrating in a finite volume, or finite difference sense? It might just be notation, but your quadrature formula appears off by either the half size first cell, or by using $i\Delta r$ rather than $(i+1/2)\Delta r$ . Jan 25 '18 at 11:52
• I integrate in a finite difference sense. I left out the equation I am starting from for brevity but now included it again for clarity. The $1/2$ originates from the discretizing the first derivative. Jan 26 '18 at 15:14

You may want a method that works right down to the coordinate singularity at $r=0$. I will do the spherical case, but the cylindrical case is similar. We want to avoid ever dividing by $r$. To think of a Finite Volume method, consider the control volumes to be concentric shells. The inner and outer radii of shell i are $(r_i-\Delta r/2, r_i+\Delta r/2)$, containing an amount of heat $$4\pi\left((r_i+\Delta r/2)^3-(r_i-\Delta r/2)^3\right)u_i=4\pi\Delta r(3r_i^2+\Delta r^2/4) \;u_i$$

Now take the first "shell" to have $r_0=\Delta r/2$ so that it is actually a little sphere of radius $\Delta r$, whose volume is given correctly by $4\pi\Delta r(3\Delta r^2/4+\Delta r^2/4)=4\pi\Delta r^3$. The control volumes are separated by interfaces that have surface area $$S_{i+1/2}=4\pi/3(r_i+\Delta r/2)^2$$ carrying heat fluxes $D(u_{i+1}-u_i)/\Delta r$.

You should find that this all works out. The secret of conservation in non-Cartesian coordinates is to to define carefully the non-overlapping control volumes and then to do the geometry without any approximation.

I think that this, and some related topics, are not well treated in the literature,so I give my version. You always have conservation in some approximate sense, if you solve the equations consistently. But as a guard against accumulation of secular errors, we usually require conservation in some exact sense. That is, if you specify the quadrature rule that you will use to check for conservation, it is guaranteed to hold exactly (It might not hold for a different quadrature) This property makes no reference to taking any infinitesimal limit. It is a property of the discrete method alone. It is a mistake to derive it by going to the limit (writing a pde) and then moving away from the limit. The limit is not involved, so dont go there.

$$(Volume)x(New \;average\; conserved \;value)=(Volume)x(Old\; average \;conserved \;variable) = \Sigma(Rate\; of \;flux \;through \;faces)$$ This works for arbitrary control volumes in any number of dimensions.

One way to see the importance of exact geometry is to consider a uniform flow, that should not change. The flux is a constant vector, and the RHS is $$\Sigma(Rate\;of\;flow\;through\;surfaces)=\Sigma(Flux \;vector)x(Normal\;to\; Surface) = (Flux Vector)x\Sigma(Normal\;to\;Surface)$$ If the geometry is calculated exactly then $\Sigma(Normal\; to\; Surface)$ vanishes for any closed body and there will be no changes. But if the body does not close there will be changes. This is a rather drastic loss of conservation.

You can have control volumes that move and distort, but you need to think of this in four dimensions.

• The approach of avoiding division by $r$ sounds quite interesting. Could you elaborate a little more on how this translates into the variable update equation - or add some link for further details/tutorial/book? Jan 26 '18 at 15:12
• @Oscillon. I have expanded my answer a bit. I can enlarge again if you wish Jan 27 '18 at 17:49
• What is not clear to me is, how is the volume accounted for in the update equation? Can you write that down explicitly? Jan 29 '18 at 18:30
• @PhilipRoe Thanks for the finite volume idea - I asked a similar question here before I saw this. Could you add the update formulas in the finite volume case? Also what does the CFL condition look like? Nov 8 '21 at 2:33

Your diffusive equation leads always to the conservation of energy in your spatial domain if Neumann BC are imposed. This is natural because there is no heat flux through walls (analogy to heat equation). Therefore integrating our equation whatever coordinate system is used leads to this conservation. For example in cylindrical coordinates: $$\frac{d}{d t} \int_{\Omega}{ur\,dr}= r\frac{\partial u}{\partial r}\hat{n}\bigg|_{\partial \Omega}=0$$

Finite element/volume schemes in 1D leads to finite differences schemes. Make sure it coincides with: $$\frac{\partial u_i}{\partial t}=\frac{D}{\Delta r^2 }[u_{i-1}-2u_{i}+u_{i+1}+\frac{\Delta r}{2r_i}(u_{i+1}-u_{i-1})] \qquad i=2,...,N-1$$ And for the first and last nodes: $$\frac{\partial u_i}{\partial t}=\frac{2D}{\Delta r^2}(u_{i+1}-u_i) \qquad i=1$$ $$\frac{\partial u_i}{\partial t}=-\frac{2D}{\Delta r^2}(u_{i}-u_{i-1})\qquad i=N$$

• I completely agree with your point that the diffusion / heat equation shows conservation when employed with Neumann boundaries. I am using the same update scheme you show (updated my question to reflect that). For the inner boundary at $r=0$, I think the correct relation is $4(u_{i+1}-u_i)$ due to l'Hospitals rule (see the link for derivation in the question). However, I tried it out and also it does not lead to conservation. Jan 26 '18 at 15:09