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$.
radial diffusion equation
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
// compilation: g++ -Wall radial_diffusion.cpp -std=c++11 -fopenmp -O3 -o radial_diffusion.exe
#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?
