# How to avoid negative values of numerical solution of transport equation using FEM scheme?

The transport equation is actually an advection-diffussion-reaction equation, which has the form as

$$\frac{\partial C}{\partial t} + v_1 \frac{\partial C}{\partial x} + v_2 \frac{\partial C}{\partial x} = D \left(\frac{\partial^2C}{\partial x^2}+\frac{\partial^2C}{\partial y^2}\right)-R(x,y)\cdot C+S$$

where $C$ is the unknown substrate concentration, $v_1$ and $v_2$ are the fluid velocities in the $x$ and $y$ direction, respectively, $D$ is the diffusion coefficient, $R$ is the reaction term, and $S$ is the source term. I use the traditional FEM scheme with backward Euler for the time advancing. But it seems that there are some negative values appearing in the numerical solution I solved. Is there any method/scheme to avoid negativeness?

Typically you would use a slope limiter (or artificial diffusion and just cross your fingers) which detects where the solution has gone negative and modifies the solution to restore positivity (often by modifying the gradient of the solution in order to maintain conservation, at least in conservative schemes like Discontinuous Galerkin and Finite Volume).

There are also more general options - Ridzal, Bochev and Shashkov have a nice trick where they solve a bound-constrained optimization problem to minimize difference between the computed solution and a new positive solution. This turns out to decouple into a cheap iterative method which is computed independently for each solution coefficient. The iteration appears to converge very rapidly.

• Does finite volume method have any advantages in this context? – boyfarrell Mar 18 '15 at 15:07
• Some folks would prefer it, or maybe someone has an existing code for FV vs FEM. I just mentioned it since the method of modifying the gradient focuses on retaining local conservation, which FV and DG boast. – Jesse Chan Mar 18 '15 at 15:13
• Why are you so sure the "negativeness" is due to a high Peclet number? Doesn't that typically cause a numerical instability that essentially ruins the solution? I think your idea of artificial diffusion (just increase D) would be a good experiment to see if the solution becomes non-negative. – Bill Greene Mar 18 '15 at 17:15
• Yes. Increasing D will make the solution nonnegative. But our setting for the initial condition is a high peak of concentration in the center and we want maintain a high peak in the center. Higer diffusion will obviouly eliminate the difference between the center and where far away from the center. – winterfly Mar 18 '15 at 18:50
• That makes it sound less like issues with Peclet number (which as Bill mentioned is a stability problem, typically if you have boundary layers, small diffusion and large convection) and more like issues with Gibbs phenomena or representation. This is usually dealt with more with limiters or shock capturing. – Jesse Chan Mar 18 '15 at 21:59