# Spurious oscillations in diffusion-reaction problems with finite volume

I have successfully solved the multi-species diffusion-reaction equation $$\frac{\partial c_i}{\partial t} = \nabla \cdot (d_i(x)\nabla c_i) + s_i(x,t), \quad \quad (1)$$ with dicontinuous source term $$s(x,t) = \left\{ \begin{array}{rcl} s_1(x,t) & \text{for} & 0<x\le x' , \\ s_2(x,t) &\text{for} & x'<x\le 1. \end{array}\right.$$ In fact, the variable diffusion coefficient also takes slightly different forms in these two region. The equation is discretized with the conservative central scheme $$c_j'(t) = \frac{1}{h^2}\biggl( d(x_{j+1/2})(c_{j+1}(t)-c_j(t))- d(x_{j-1/2})(c_{j}(t)-c_{j-1}(t))\biggr) + s(x_j,t).$$ I used SUNDIALS package for time integration.

Now I'm trying to solve the problem in a concentrated solution $$\frac{\partial c_i}{\partial t} = \nabla \cdot \sum_k D_{ik}(x)\nabla c_k + s(x,t),$$ where $D_{ik}$ may be negative. In dilute solution approximation (solvent concentration $c_0 \rightarrow \infty$), we have $D_{ik} \rightarrow d_{i}\delta_{ik}$. However, the solution exhibit spurious oscillations.

The blue line is sloved with $c_0=10^{20}$ or equation (1); purple and green lines with $c_0=10^{10}$. The difference between purple and green line is the average scheme used for the diffusion coefficient at the cell interface, the former is with simple average and the latter with harmonic average. As shown in the figures below, there are oscillations near the boundaries and the point $x'$. (I should mention that the vary first osscilations began at the left boundary before they built up at the right boundary and the point $x'$).

It is odd that even with the oscillations, the scheme with harmonic average never fails to converge at each timestep and follows closely the result of analytical dilution scheme (equation (1)); where as the simple average scheme fails at the one tenth of the time domain. At much later time the solutions look like

Are there discretizatons within the finite volume scheme that avoid such oscillations?

• is your mesh aligned with the discontinuity? the harmonic mean flux is meant to deal with discontinuities in the diffusion coefficient $d$, see this reference. How do you deal with the boundaries? – GoHokies Mar 9 '17 at 12:18
• I'm using the cell centered grid, so I aligned the discontinuity at a control-volume face, say point $j+1/2$. The coefficient is then the harmonic mean of the value at $j$, $j+1$. Since I have pure Neumann boundary conditions, they fit in to the cell centered grid naturally. – Henry Mar 9 '17 at 13:54
• Pure Neumann BCs may be problematic here, you may need to check that you're satisfying the compatibility condition (so the solution does not blow up). I suggest you consider normalizing your variables, since these large numbers ($10^{10}, 10^{20}$) are more susceptible to significant precision errors. – Charles Mar 9 '17 at 18:27
• Initial conditions are consistent. And I have tried to rewrite the code with gmp, but it did not help. – Henry Mar 10 '17 at 3:25
• @GoHokies are there special techniques to align with the discontinuity? I though negative diffusion coefficients may as well cause some problems as discussed here. – Henry Mar 10 '17 at 6:18