search

SLACKHA / pyJac

Creates C and CUDA analytical Jacobians for chemical kinetics ODE systems
http://slackha.github.io/pyJac/
MIT License
52 stars 23 forks source link

about 98% of the elements are not zero? #24

Closed ghost closed 6 years ago

ghost commented 6 years ago

about 98% of the elements are not zero?

import pyjac as pj
import cantera as ct
import numpy as np
import pyjacob 
from numpy.linalg import eig
from numpy import argmax,count_nonzero

#create gas from original mechanism file gri30.cti
gas = ct.Solution('gri30.cti')
#reorder the gas to match pyJac
n2_ind = gas.species_index('N2')
specs = gas.species()[:]
specs2 = specs[:n2_ind] + specs[n2_ind + 1:] + [specs[n2_ind]]

names = []
print("index{}:{}".format(0,'T'))
names.append('T')
for i,s in enumerate(specs2[:-1]):
    print("index{}:{}".format(i+1,s.name))
    names.append(s.name)
gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
        species=specs2,
        reactions=gas.reactions())

# generate the c code to ./out folder
pj.create_jacobian(lang='c',gas=gas)

T = 1500
P = ct.one_atm
gas.TPY = T, P, "CH4:1.0, O2:2, N2:7.52"
r = ct.IdealGasConstPressureReactor(gas, name='R1')
sim = ct.ReactorNet([r])
t=0.0
print('time, temp, press, internal energy')

t=0.0
print('time, temp, press, internal energy')

for t in np.arange(0,2e-3,5e-6):
    sim.advance(t)
    print('%10.3e %10.3f %10.3f %14.6e' %(sim.time, r.T, r.thermo.P, r.thermo.u))
    #setup the state vector
    y = np.zeros(gas.n_species)
    y[0] = r.T
    y[1:] = gas.Y[:-1]

    #create a dydt vector
    dydt = np.zeros_like(y)
    pyjacob.py_dydt(0, P, y, dydt)

    #create a jacobian vector
    jac = np.zeros(gas.n_species * gas.n_species)

    #evaluate the Jacobian
    pyjacob.py_eval_jacobian(0, P, y, jac)

    # reshape
    jac = jac.reshape((gas.n_species,gas.n_species)).T

    # get eigen
    ev,eV = eig(jac)
    print(max(ev.real),
          min(ev.real),
          ev[argmax(ev.imag)],
          norm(jac),
          count_nonzero(jac)*1.0/np.prod(jac.shape))

output looks like:

time, temp, press, internal energy
 0.000e+00   1500.000 101325.000   7.871386e+05
4930.71181537 -439421094.945 (-0.529376424808+0.124547878742j) 3.03427476721e+12 0.223211107156
 5.000e-06   1499.995 101325.000   7.871402e+05
4827.60926329 -439421090.26 (-0.914795644697+0.232602551787j) 3.03428268642e+12 0.981132075472
 1.000e-05   1499.989 101325.000   7.871418e+05
4727.37493927 -439421124.204 (-1.28101057506+0.325071212446j) 3.03429057428e+12 0.981132075472
 1.500e-05   1499.984 101325.000   7.871433e+05
4627.38763267 -439421204.107 (-1.66572622569+0.404929343758j) 3.03429850855e+12 0.981132075472
 2.000e-05   1499.978 101325.000   7.871448e+05
4527.45372387 -439421334.374 (-2.06901445963+0.46448002886j) 3.03430641503e+12 0.981132075472
 2.500e-05   1499.973 101325.000   7.871463e+05
4427.57335683 -439421517.527 (-2.49128133606+0.490524785492j) 3.03431419575e+12 0.981132075472
 3.000e-05   1499.968 101325.000   7.871477e+05
4327.82812694 -439421755.309 (-2.93308784638+0.45862546959j) 3.03432174209e+12 0.981132075472
 3.500e-05   1499.962 101325.000   7.871491e+05
4228.33622426 -439422049.121 (-3.39510035357+0.296650050187j) 3.03432894018e+12 0.981132075472
 4.000e-05   1499.957 101325.000   7.871504e+05
4129.23455195 -439422400.186 (-0.51433469421+0.0857915491796j) 3.03433567328e+12 0.981132075472
 4.500e-05   1499.952 101325.000   7.871516e+05
4030.67147569 -439422809.612 (-0.563803708368+0.0943895508633j) 3.03434182281e+12 0.981132075472
 5.000e-05   1499.948 101325.000   7.871528e+05
3932.80361075 -439423278.409 (-0.61372794618+0.102058569703j) 3.03434726902e+12 0.981132075472
 5.500e-05   1499.943 101325.000   7.871538e+05
3835.79397938 -439423807.493 (-0.663938985189+0.108979897334j) 3.03435189142e+12 0.981132075472
 6.000e-05   1499.939 101325.000   7.871548e+05
3739.81046582 -439424397.682 (-0.714256793073+0.115260242347j) 3.03435556916e+12 0.981132075472
 6.500e-05   1499.935 101325.000   7.871557e+05
3645.02415979 -439425049.694 (-0.76450892722+0.120966937934j) 3.03435818138e+12 0.981132075472
 7.000e-05   1499.932 101325.000   7.871564e+05
3551.60745544 -439425764.148 (-0.814545373793+0.126140875325j) 3.03435960762e+12 0.981132075472
 7.500e-05   1499.928 101325.000   7.871571e+05
3459.73189239 -439426541.562 (-0.864246217858+0.130797984995j) 3.03435972819e+12 0.981132075472
 8.000e-05   1499.926 101325.000   7.871576e+05
3369.56579507 -439427382.359 (-0.913522015029+0.134926788484j) 3.03435842453e+12 0.981132075472
 8.500e-05   1499.923 101325.000   7.871580e+05
3281.27179126 -439428286.868 (-15.218755323+0.926382878656j) 3.03435557957e+12 0.981132075472
 9.000e-05   1499.921 101325.000   7.871582e+05
3195.00431014 -439429255.333 (-15.9759663086+1.84158794689j) 3.03435107813e+12 0.981132075472
 9.500e-05   1499.919 101325.000   7.871583e+05
3110.90716522 -439430287.92 (-16.7525128817+2.36871575832j) 3.03434480722e+12 0.981132075472
 1.000e-04   1499.918 101325.000   7.871583e+05
3029.11132546 -439431384.728 (-17.5472350549+2.72877112317j) 3.03433665636e+12 0.981132075472
 1.050e-04   1499.917 101325.000   7.871580e+05
2949.73296806 -439432545.798 (-18.3585188522+2.97072343796j) 3.03432651791e+12 0.981132075472
 1.100e-04   1499.917 101325.000   7.871577e+05
2872.87189069 -439433771.123 (-19.1844613156+3.11044242447j) 3.03431428728e+12 0.981132075472
 1.150e-04   1499.917 101325.000   7.871571e+05
2798.61034056 -439435060.66 (-20.0229175483+3.15009172636j) 3.03429986317e+12 0.981132075472
 1.200e-04   1499.918 101325.000   7.871564e+05
2727.0122924 -439436414.342 (-20.8715059653+3.08226882463j) 3.03428314779e+12 0.981132075472
 1.250e-04   1499.920 101325.000   7.871555e+05
2658.12318804 -439437832.083 (-21.7275837764+2.88820619856j) 3.03426404699e+12 0.981132075472
 1.300e-04   1499.922 101325.000   7.871544e+05
2591.97011928 -439439313.795 (-22.5881773319+2.52754975503j) 3.03424247039e+12 0.981132075472
 1.350e-04   1499.924 101325.000   7.871531e+05
2528.56242112 -439440859.39 (-23.4497856317+1.89467421874j) 3.03421833146e+12 0.981132075472
 1.400e-04   1499.927 101325.000   7.871516e+05
2467.89262137 -439442468.791 (-1.46178141469+0.0948506986798j) 3.03419154758e+12 0.981132075472
 1.450e-04   1499.931 101325.000   7.871500e+05
2409.9376855 -439444141.936 (-1.50406695282+0.0690416476984j) 3.0341620401e+12 0.981132075472
 1.500e-04   1499.935 101325.000   7.871481e+05
2354.66048342 -439445878.784 (-22.8016930027+0.249607323251j) 3.03412973426e+12 0.981132075472
...
kyleniemeyer commented 6 years ago

Hi @chengdi123000, yes, for the Jacobian based on a constant-pressure system, the Jacobian is dense—almost all elements are nonzero. This happens because all species mass fractions contribute to density, which shows up in the time derivative of species mass fraction, and thus all species are "connected" to each other in the Jacobian. You can read a bit more about this in Section 3.3 of our article describing pyJac:

Kyle E. Niemeyer, Nicholas J. Curtis, and Chih-Jen Sung. 2017. “pyJac: analytical Jacobian generator for chemical kinetics.” Computer Physics Communications, 215:188–203. https://doi.org/10.1016/j.cpc.2017.02.004 | Available at https://arxiv.org/abs/1605.03262

We are currently working on v2 of pyJac which will enable a sparse formulation—this is almost ready to share! For updates, please follow this repository, or join the SLACKHA user group here: https://groups.io/g/slackha-users/join (the mailing list is extremely low volume).

michael-a-hansen commented 6 years ago

You could try counting nonzeros by a threshold. Something like np.sum(np.abs(J) > 1.e-12) ought to do the trick (I haven't tested that).

kyleniemeyer commented 6 years ago

I should also point out that we currently formulate the exact Jacobian; it may be possible to increase the sparsity of the current constant-pressure formulation by introducing some approximations, but we haven't done this yet.

kahilah commented 6 years ago

I once played around and tested to increase sparsity of the Jacobian by easily applying different threshold levels after which the element is forced to zero. I used gri3 and jacobian was created with the Cp assumption in pyjac. I noticed that the ODE solver was able to solve the problem but it took quite a bit more internal steps, making the overall cpu-time higher than with a dense and accurate jacobian. For the sparse tests I applied sparse linear algebra (KLU package). The greater number of integration steps + medium sparsity levels were not enough to make any speedup.

So lessons learned: zeroing jacobian elements by thresholding may not be a good idea for a priori dense jacobian. At least I'm doubting it's generality. Story could be different with larger mechanisms...

Great to hear that pyjac V2 is soon ready.

Cheers, HK

ghost commented 6 years ago

@kyleniemeyer Thank you very much for your reply!

@michael-a-hansen Thresholding is not a good idea. It depends on what the Jacobian is used for. For quasi-Newton methods, inexact Jacobian is useful with the cost of more steps. For eigen value and vectors, sometime it becomes sensitive to the small numbers for chemical kinetics. Considering energy in the system, small perturbation of radical may influence heavily.

michael-a-hansen commented 6 years ago

@chengdi123000 agreed. If the goal is to estimate sparsity, however, it might be useful to use a threshold to eliminate numbers on the order of machine epsilon. That's all I was suggesting.

skyreflectedinmirrors commented 6 years ago

@chengdi123000 I'm going to close this issue as it appears to be resolved

ghost commented 6 years ago

@arghdos OK