Extending the solution components causes error "Jacobian is singular"

[Ewald Kleefstra]

Dear everyone

I have made changes to the Cantera source code to extend the solution components in the solution array with a new parameter "soot" for counterflow diffusion flames.

When running the simulation via Python, I get the following error:

Attempt Newton solution of steady-state problem...
***********************************************************************
CanteraError thrown by MultiNewton::step:
Jacobian is singular for domain flame, component soot at point 0
(Matrix row 5)
***********************************************************************

Traceback (most recent call last):
  File "diffusion_b.py", line 67, in <module>
    f.solve(loglevel, auto=True)
  File "C:\Users\Ewald\Anaconda3\envs\diffusion\lib\site-packages\cantera\onedim.py", line 914, in solve
    super().solve(loglevel, refine_grid, auto)
  File "interfaces\cython\cantera\onedim.pyx", line 992, in cantera._cantera.Sim1D.solve
cantera._cantera.CanteraError: Could not find a solution for the 1D problem

What could cause this problem? I am using Cantera v2.4.0 and Anaconda 3 on Windows 10. My code changes can be found here.

My Python file "diffusion_b.py" follows the onedim/diffusion_flame.py example file with my added toggle to enable soot formation calculation.

I am rather new to computational equation solving and C++, so any help would be appreciated.

Kind regards

Ewald Kleefstra

[Ewald Kleefstra]

I should add that the error is thrown regardless of the soot toggle f.soot_enabled in my "diffusion_b.py" example being set to true or not.

[Ray Speth]

Ewald,

What is the equation you are intending to solve at the left and right boundaries? The error is because the equation that you have right now,

rsd[index(c_offset_S,0)] = m_soot_left;

does not depend on any value in the state vector x. If you want to impose a fixed boundary value, then your residual equation should be for the difference between the desired and actual values at the boundary. This would look something more like

rsd[index(c_offset_S, 0)] = m_soot_left - N_p(x, 0);

And of course, you would need to do something similar for the right boundary.

Regards,
Ray

[Ewald Kleefstra]

Dear Ray

Thank you for your reply. The Jacobian Singularity error is now not occuring anymore. While Cantera tries to solve the newly implemented equations, I now receive the following damping error:

Re-evaluating Jacobian, since no damping coefficient
could be found with this Jacobian.

0    1.00000     0.00000    19.00104    10.59517    10.45733  103  1/20
At limits.

Re-evaluating Jacobian, since no damping coefficient
could be found with this Jacobian.

At limits.
...failure.
    0    1.49e-13       4.743

What have I done wrong? I have attached a pdf with the equations I tried to implement into Cantera. The updated code can be found again on my Github.

Kind regards

Ewald

[Ewald Kleefstra]

Dear Ray,

It seems that my m_offset_S solution vector index remained 0 because my soot toggle is disabled by default. Therefore m_offset_S is not refreshed in the if clause at line 58 of the StFlow.cpp file during StFlow class initialization. The m_offset_S index being 0 caused overwriting of the original solution vector index 0, which probably caused the damping error. I will try to fix this and give an update.

Kind regards,
Ewald

[Ewald Kleefstra]

Good evening Ray

In my latest commit m_offset_S is now correctly set and Newton solution attempts are being made. However, during each solution attempt I get the following two errors until the code at step "Solving on 24 point grid with energy equation disabled" abruptly comes to a halt.

Damped Newton iteration:
-----------------------------------------------------------------
m     F_damp     F_bound     log10(ss)   log10(s0)   log10(s1)  N_jac   Age
-----------------------------------------------------------------

ERROR: solution out of bounds.
domain 1:               eField(1) =  1.000e+00 ( 0.000e+00,  0.000e+00)

ERROR: solution out of bounds.
domain 1:                 soot(2) =  3.000e+02 ( 0.000e+00,  0.000e+00)


***********************************************************************
CanteraError thrown by Phase::setTemperature:
temperature must be positive. T = nan

***********************************************************************

Traceback (most recent call last):

  File "DIFFUSION_FLAME.py", line 67, in <module>

    f.solve(loglevel, auto=True)
  File "C:\Users\Ewald\Anaconda3\envs\diffusion\lib\site-packages\cantera\onedim.py", line 914, in solve
    super().solve(loglevel, refine_grid, auto)

  File "interfaces\cython\cantera\onedim.pyx", line 935, in cantera._cantera.Sim1D.solve
  File "C:\Users\Ewald\Anaconda3\envs\diffusion\lib\site-packages\cantera\onedim.py", line 855, in set_initial_guess
    self.gas.equilibrate('HP')

I do not understand why the soot value is out of bounds even though it should be in between the set boundaries, whereas eField should not even be affected by my code changes. And where could a temperature value nan be passed?
I assume not all affected arrays/vectors have the proper length. I try to resize and set boundaries within the initSoot() function in StFlow.cpp at line 103.

Kind regards

Ewald

[Ewald Kleefstra]

The setBounds error was caused because the solution array was not resized yet. The T = nan error is caused by a wrong equation implementation in the StFlow::updateParticleDiffFluxes() method causing a division by zero at StFlow.cpp line 640, where N_p(x,0) is zero. I am currently looking into it.

Best regards

Ewald

[Ewald Kleefstra]

Dear everyone

I have fixed the division by zero by defaulting the solution vector x[index(m_offset_S, j)] for each j to 1 at StFlow::_getInitialSol().

Now I get my first post's "Jacobian is singular" error again. Unlike my first post however, the residual vector's entries should all be dependent on the state vector x.

As soon as I comment out line 642 of StFlow.cpp the error message disappears and Cantera tries to solve the equations.

Below is a snippet of the console error output.

************ Solving on 6 point grid with energy equation enabled ************

..............................................................................

Attempt Newton solution of steady-state problem...

***********************************************************************
CanteraError thrown by MultiNewton::step:

Jacobian is singular for domain flame, component soot at point 3
(Matrix row 467)
***********************************************************************


******** Initial solve failed; Retrying with energy equation disabled ********

..............................................................................

Attempt Newton solution of steady-state problem...

Maximum Jacobian age reached (20)

***********************************************************************
CanteraError thrown by MultiNewton::step:

Jacobian is singular for domain flame, component soot at point 3
(Matrix row 467)
***********************************************************************


*********** Solving on 12 point grid with energy equation enabled ************

..............................................................................

Attempt Newton solution of steady-state problem...

***********************************************************************
CanteraError thrown by MultiNewton::step:

Jacobian is singular for domain flame, component soot at point 1
(Matrix row 233)
***********************************************************************

I am close, but what is missing?

Kind regards

Ewald

[Ewald Kleefstra]

Dear everyone

To specify the problem of my previous post, it is probably best to have a look at lines 476 to 485 of StFlow.cpp, where the residual of my new equations is being calculated.

The there mentioned variable m_particle_flux[j] is set in the method updateParticleDiffFluxes().

When line 644 in method updateParticleDiffFluxes is not commented out, the "domain singular" error in my previous post is being thrown.

If it is commented out however, Cantera tries to endlessly solve while throwing the "damping" error below.

Re-evaluating Jacobian, since no damping coefficient
could be found with this Jacobian.


Any help and hints are appreciated.

Kind regards

Ewald