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.

**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.

**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*,

*Santiago*