Transport Phenomena problem coupled by Chemical Reactions: Negative values in Finite-element functions lead to numerical divergence

Hello FEniCS community,

I am solving a mathematical system of Transport Phenomena Coupled by Chemical Reactions as follows:

The PDEs system and thus, the variational formulation is composed by 7 Finite-Element Functions (5 species - mass transport equation, Temperature - energy transport equation, and the last accounting for pressure drop - Ergun equation). Besides, the coupled chemical system considers three reversible reactions inside a catalytic fixed-bed reactor.

  1. Problem:

Solver fails at some point during the dynamic simulation with the following message error:

  Process 0: Newton iteration 0: r (abs) = 2.798e+12 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Process 0: Newton iteration 1: r (abs) = 9.642e+18 (tol = 1.000e-10) r (rel) = 3.446e+06 (tol = 1.000e-09)
  Process 0: Newton iteration 2: r (abs) = 2.936e+31 (tol = 1.000e-10) r (rel) = 1.049e+19 (tol = 1.000e-09)
  Process 0: Newton iteration 3: r (abs) = 6.031e+36 (tol = 1.000e-10) r (rel) = 2.156e+24 (tol = 1.000e-09)

*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PCSETUP_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 1
*** DOLFIN version: 2018.1.0

The FEniCS - Python code is too extensive to post it over here. However, the command line that computes the numerical solution is:

    solve(F == 0, u, bcs, solver_parameters={"newton_solver": {
            "relative_tolerance": 1e-5}, "newton_solver": {
                    "maximum_iterations": 250}})

After diagnosing the behaviour of the system, It was evidenced that at some point, the functions which account for the reaction rates started to fluctuate and over the course, those oscillations became stronger resulting in negative concentrations for some species. Which, of course, leads to numerical divergence.

  1. Questioning:
  • Do mass- and energy-transport assure a “physical auto-control” (non-negative compositions)?
  • If not: How can I restrict my system so that, non-negative numerical solutions may be computed for all Finite-Element Functions?
  • Any suggestion about the Non-linear solver? Should I implement another or specify different arguments to the one above?

Best Regards,


it is a well known phenomenon that standard Finite Element discretizations and Advection-Diffusion equations don’t get along well, where problems including artificial oscillations and numerical damping may arise.
There are many techniques to somewhat obviate these issues. A nice overview is this work. There are also a couple of papers on this topic by G. Wells, one of the lead developers of Dolfin, e.g. this one.

Hi Klunkean,

Thank you for your fast response. I will be checking the documentation you shared and I will let you know.

Best regards,


You need to make it sure that the global mass balance is fulfilled as well. Try to get rid of alpha=5 and solve mass balance instead. As sum of partial masses is 1 you can alway get alpha=5 in postprocessing. All difficulties about FEM in Navier–Stokes equation needs to be handled with care in your case. A mathematician calls Navier–Stokes if it is an incompressible flow with a Newton fluid. This fact is definitely not your case, you are solving a compressible flow. The biggest problem is the ideal gas equation if you use for pressure, there are also Grüneisen type of energy/pressure relations increasing the robustness.

For your specific question, FEniCS has a nice

conditional(le(partialmass,0.),0.,1.) * partialmass

function, if partialmass is less or equal than 0 then the outcome is 0 else 1 multiplied by partialmass.

All the best, Emek