I am solving an advection-diffusion equation using the FEM and am having trouble calculating my fluxes.

I start with the equation,

$$\frac{\partial n}{\partial t} = \frac{\partial j_{n}}{\partial x}\, ,$$

Where my flux (current density) is given by,

$$j_n = \mu E n + D \frac{\partial n}{\partial x}\, ,$$

where $\mu$ is my drift velocity and $D$ the diffusion coefficient.

Solving this equation using FEM-FCT to ensure stability I obtain my nodal values of n, which I then use to solve for my fluxes.

The issue I am having currently is that my steady state current density should be zero, where an in-built field causes the drift and diffusion terms to balance, however I am obtaining non-zero current densities.

In order to determine the cause of this, I plotted the two components of the current density in the following plot, and as is seen the two are not equal.

Drift and Diffusion Flux at steady state.

I am wondering if it has to do with how the drift flux uses the calculated $n$ directly while the diffusion flux requires a derivative which will introduce error.

Any help in generating more accurate fluxes would be appreciated.

A few notes: FEM-FCT limits me to linear basis functions, however I tried to create a higher order first derivative matrix just using finite differences however that did not seem to fix the problem.

  • $\begingroup$ Welcome to SciComp.SE, for future reference you can type your equations right in your question since it supports MathJax. $\endgroup$ – nicoguaro Mar 12 '19 at 13:54
  • 1
    $\begingroup$ Regarding your question, you mean that the two lines in your plot are different, right? At first glance, they look (about) the same to me. $\endgroup$ – nicoguaro Mar 12 '19 at 13:55
  • $\begingroup$ They look close but begin to deviate towards the boundary. The boundary itself isn't a problem as I know the flux on the boundary. The main issue is that I would like the fluxes to be as accurate as possible, which in this case should equate to equal drift and diffusion fluxes, I can simply increase the spatial resolution but this too only helps to a certain degree. $\endgroup$ – MGreen Mar 12 '19 at 15:24
  • $\begingroup$ In "regular" FEM derivatives are not continuous and you need to do some smoothing afterward if you want to visualize it as a smooth field. I don't know about FEM-FCT (Flux-corrected transport), though. $\endgroup$ – nicoguaro Mar 12 '19 at 15:48
  • $\begingroup$ Okay, would you be able to point me to some material that would help me implement this? $\endgroup$ – MGreen Mar 12 '19 at 16:05

If you are using a (classical) finite element formulation you have $C^0$ continuity, that is, your function is continuous but its derivatives are not. Thus, it is expected to have that behavior towards the end, since there the derivatives are larger. As you refine these two curves will get closer to each other.

As mentioned in my comment, you can do some smoothing of your derivatives. I have used the following approach (see [1]).:

  1. The variable is extrapolated from the Gauss point to nodes for each element.
  2. These values are averaged according to the number of elements that share that node.


  1. Zienkiewicz, Olgierd Cecil, and Jian Zhong Zhu. "The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique." International Journal for Numerical Methods in Engineering 33.7 (1992): 1331-1364.
| cite | improve this answer | |

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.