diff --git a/CEMA.ipynb b/CEMA.ipynb index 3ca6a21..c869da1 100644 --- a/CEMA.ipynb +++ b/CEMA.ipynb @@ -2,275 +2,130 @@ "cells": [ { "cell_type": "markdown", - "id": "21cc6738", + "id": "f080cd4e", "metadata": {}, "source": [ - "# A 1-D example of Chemical Explosive Mode Analysis (CEMA) for $CH_4$ ignition" + "# 1-D Chemical Explosive Mode Analysis (CEMA) validation (Lu et.al.) for stoichiometric $H_2-air$ combustion ($\\phi$ =1) " ] }, { "cell_type": "markdown", - "id": "c984d88d", + "id": "5a9db857", "metadata": {}, "source": [ - "Let us implement an sample CEMA for methane combustion case.To do so we first generate the 1-D methane combustion data using Arrhenius.jl" + "Chemical Explosive Mode Analysis (CEMA) has proven to be a useful computational diagnostic tool for analysis combustion data. It can help identify regimes of auto-ignition and flame propagation using an eigen analysis on the chemical Jacobian. Arrhenius.jl can be exploited to perform eigen-analysis thanks to the built in tools for calculating jacobians and eigenvalues/eigenvectors of a matrix. " ] }, { "cell_type": "code", - "execution_count": 2, - "id": "9989f948", + "execution_count": 29, + "id": "0b874589", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "PyObject " + "171×11 Matrix{Float64}:\n", + " 0.01 0.295858 0.147929 … 1.97804e-13 0.556213 1001.0\n", + " 0.02 0.295858 0.147929 1.0614e-12 0.556213 1001.0\n", + " 0.03 0.295858 0.147929 3.31954e-12 0.556213 1001.0\n", + " 0.04 0.295858 0.147929 8.50956e-12 0.556213 1001.0\n", + " 0.05 0.295858 0.147929 1.98756e-11 0.556213 1001.0\n", + " 0.06 0.295858 0.147929 … 4.42671e-11 0.556213 1001.0\n", + " 0.07 0.295858 0.147929 9.61617e-11 0.556213 1001.0\n", + " 0.08 0.295858 0.147929 2.06229e-10 0.556213 1001.0\n", + " 0.09 0.295858 0.147929 4.39757e-10 0.556213 1001.0\n", + " 0.1 0.295857 0.147928 9.37095e-10 0.556213 1001.0\n", + " 0.11 0.295856 0.147928 … 2.00602e-9 0.556213 1001.01\n", + " 0.12 0.295854 0.147926 4.34759e-9 0.556214 1001.01\n", + " 0.13 0.29585 0.147923 9.66824e-9 0.556215 1001.03\n", + " ⋮ ⋱ ⋮\n", + " 0.894688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 0.904688 0.0356073 0.0126109 … 3.45191e-7 0.631019 2691.91\n", + " 0.914688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 0.924688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 0.934688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 0.944688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 0.954688 0.0356073 0.0126109 … 3.45191e-7 0.631019 2691.91\n", + " 0.964688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 0.974688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 0.984688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 0.994688 0.0356073 0.0126109 3.45191e-7 0.631019 2691.91\n", + " 1.00469 0.0356073 0.0126109 … 3.45191e-7 0.631019 2691.91" ] }, - "execution_count": 2, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Arrhenius\n", - "using ForwardDiff\n", - "using DifferentialEquations\n", "using LinearAlgebra\n", + "using DifferentialEquations\n", + "using ForwardDiff\n", + "using DiffEqSensitivity\n", + "using Sundials\n", "using Plots\n", - "using PyCall\n", - "ct = pyimport(\"cantera\")" + "using DelimitedFiles\n", + "using Profile\n", + "cantera_data = readdlm(\"ct_data.dat\")" ] }, { "cell_type": "code", - "execution_count": 5, - "id": "10536c49", + "execution_count": 31, + "id": "c24a80eb", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "Arrhenius.Solution(53, 325, [2.016, 1.008, 15.999, 31.998, 17.007, 18.015, 33.006, 34.014, 12.011, 13.018999999999998 … 43.025, 43.025, 43.025, 42.016999999999996, 28.014, 39.95, 43.089, 44.097, 43.045, 44.053], [\"H2\", \"H\", \"O\", \"O2\", \"OH\", \"H2O\", \"HO2\", \"H2O2\", \"C\", \"CH\" … \"HCNO\", \"HOCN\", \"HNCO\", \"NCO\", \"N2\", \"AR\", \"C3H7\", \"C3H8\", \"CH2CHO\", \"CH3CHO\"], [\"O\", \"H\", \"C\", \"N\", \"Ar\"], [0.0 0.0 … 1.0 1.0; 2.0 1.0 … 3.0 4.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], Arrhenius.IdealGasThermo([2.34433112 0.00798052075 … -917.935173 0.683010238; 2.5 7.05332819e-13 … 25473.6599 -0.446682853; … ; 3.409062 0.010738574 … 1521.4766 9.55829; 4.7294595 -0.0031932858 … -21572.878 4.1030159], [3.3372792 -4.94024731e-5 … -950.158922 -3.20502331; 2.50000001 -2.30842973e-11 … 25473.6599 -0.446682914; … ; 5.97567 0.008130591 … 490.3218 -5.045251; 5.4041108 0.011723059 … -22593.122 -3.4807917], [200.0 1000.0 3500.0; 200.0 1000.0 3500.0; … ; 300.0 1000.0 5000.0; 200.0 1000.0 6000.0], false), Arrhenius.Transport(7, [-2.366640508787296e-10 -1.2299135731442892e-9 … 8.693592178691943e-8 8.794793602234737e-8; 2.1192301406670526e-8 6.668825847991541e-8 … -3.271068954511353e-6 -3.309147211220777e-6; … ; 5.7448300967228756e-5 0.0001298839327866297 … -0.0049138820317118 -0.00497108414638971; -6.404100632763197e-5 -0.00013572755355049403 … 0.0051027272253661325 0.005162127672107992], [-0.0016279891577270268 7.206309409434899e-5 … -0.00011924922331411204 -0.00015757657556546702; 0.06178549433963653 -0.0021239618978405594 … 0.004243684109832842 0.005587634198577984; … ; 96.65522170663003 -2.0666287999557595 … 5.836219561810777 7.641767878764346; -102.92280467937866 1.8536466578657944 … -5.97664208790569 -7.812749492011645], [3.348485446466861e-5 5.169493728638671e-5 … 3.2144250898910778e-6 3.195770414642371e-6; -0.0012612725705465348 -0.0019481129061467706 … -0.00012272386962752324 -0.00012201165084214482; … ; -2.0626595269712937 -3.1895799986447386 … -0.20978320877300424 -0.20856574763356261; 2.2014391913047877 3.404772178136124 … 0.22590225587600743 0.2245912490539997]), Arrhenius.Reaction(\n", - "⢛⣛⢺⣻⠚⡛⣍⢙⡋⢉⡟⣿⡓⢓⢶⡍⣙⣼⣧⣋⣓⢛⠛⠛⠛⢛⢛⣁⣙⢌⠙⡫⣚⣋⣊⡋⣺⣟⢓⡐\n", - "⠈⠉⠱⠥⠁⠀⠈⠛⠹⠦⠅⠙⠉⡶⠍⠻⠷⠆⠮⠭⠎⠗⠀⣀⢀⡀⢀⣀⣀⣀⡁⠈⣀⣀⡀⠙⠂⠩⠁⠭\n", - "⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠀⠀⠀⢈⣉⠙⣒⣑⢑⢛⡲⡲⡶⣞⣂⠶⡛⠀⠀⠀⣀⠀\n", - "⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠈⠉⠁⠉⠉, \n", - "⠹⣉⠉⢙⠛⡻⣍⠉⠉⠉⡳⢖⡒⠒⢖⣷⣼⣜⣂⣆⡒⠓⠛⠑⠚⢛⠛⠓⠙⢥⢤⡙⠙⢉⣋⣷⠛⡛⣳⢓\n", - "⠀⠈⠓⠦⠁⠀⠈⠉⠓⠦⡄⠀⠙⠲⠌⠀⠅⠄⠠⠘⠝⢦⡀⣀⣀⠀⠀⢀⠀⡀⠀⠀⠀⢠⣠⠐⠂⠄⠤⠀\n", - "⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠀⠀⠀⠈⠚⠁⠉⠉⠒⣊⡓⣊⠉⠥⠤⠌⢋⠀⠀⠀⢀⣀\n", - "⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠀, \n", - "⠹⣉⠉⢙⠛⡻⣍⠉⠉⠉⡳⢖⡒⠒⢖⣷⣼⣜⣂⣆⡒⠓⠛⠑⠚⢛⠛⠓⠙⢥⢤⡙⠙⢉⣋⣷⠛⡛⣳⢓\n", - "⠀⠈⠓⠦⠁⠀⠈⠉⠓⠦⡄⠀⠙⠲⠌⠀⠅⠄⠠⠘⠝⢦⡀⣀⣀⠀⠀⢀⠀⡀⠀⠀⠀⢠⣠⠐⠂⠄⠤⠀\n", - "⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠀⠀⠀⠈⠚⠁⠉⠉⠒⣊⡓⣊⠉⠥⠤⠌⢋⠀⠀⠀⢀⣀\n", - "⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠀, Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 0, 1], [1.2e11 -1.0 0.0; 5.0e11 -1.0 0.0; … ; 2.41e10 0.0 0.0; 1.927e10 -0.32 0.0], [6.02e8 0.0 3000.0; 1.04e20 -2.76 1600.0; … ; 3.0000000000000003e57 -14.6 18170.0; 4.42e55 -13.545 11357.0], [0.562 5836.0 8552.0 91.0; 0.783 2941.0 6964.0 74.0; … ; 0.1894 8748.0 7891.0 277.0; 0.315 3285.0 6667.0 369.0], [1, 2, 33, 39, 43, 167, 187, 205, 212, 227, 230, 269], [12, 50, 52, 54, 56, 57, 59, 63, 70, 71 … 158, 174, 185, 237, 241, 289, 304, 312, 318, 320], [-1, 1, 2, 3, 4, 5, 6, 7, 8, 9 … 19, 20, -1, -1, 21, 22, 23, 24, 25, 26], \n", - "⡇⢸⠀⠀⣿⡇⣿⣿⢸⡇⡇⢸⠀⠀⠀⠀⡇⡇⡇⡇⡇⡇⢸⠀⠀⣿⠀⢸⡇⣿⠀⠀⢸⠀⠀⡇⠀⡇⣿⡇\n", - "⡇⢸⠀⠀⣿⡇⣿⣿⢸⡇⡇⢸⠀⠀⠀⠀⡇⡇⡇⡇⡇⡇⢸⠀⠀⣿⠀⢸⡇⣿⠀⠀⢸⠀⠀⡇⠀⡇⣿⡇\n", - "⡇⢸⠀⠀⣿⡇⣿⣿⢸⡇⡇⢸⠀⠀⠀⠀⡇⡇⡇⡇⡇⡇⢸⠀⠀⣿⠀⢸⡇⣿⠀⠀⢸⠀⠀⡇⠀⡇⣿⡇\n", - "⠁⠈⠀⠀⠉⠁⠉⠉⠈⠁⠁⠈⠀⠀⠀⠀⠁⠁⠁⠁⠁⠁⠈⠀⠀⠉⠀⠈⠁⠉⠀⠀⠈⠀⠀⠁⠀⠁⠉⠁, [[3], [2, 3], [1, 3], [3, 7], [3, 8], [3, 10], [3, 11], [3, 12], [3, 12], [3, 13] … [8, 50], [13, 51], [13, 25], [3, 50], [2, 50], [2, 50], [5, 50], [7, 50], [7, 50], [13, 50]], [[4], [5], [2, 5], [4, 5], [5, 7], [2, 15], [2, 17], [1, 15], [2, 17], [2, 18] … [7, 51], [14, 50], [50], [18, 26], [51], [13, 26], [19, 26], [4, 51], [5, 18, 26], [26]], 325, \n", - "⢻⣛⢻⣻⠛⠻⣍⢙⡋⢉⡿⣿⡓⢓⢶⣿⣽⣼⣧⣏⣑⢛⠛⠛⠛⢛⢛⣓⣙⢭⢽⡻⣛⣋⣋⣿⣻⣟⣳⣓\n", - "⠈⠉⠳⠧⠁⠀⠈⠛⠻⠦⡅⠙⠙⡶⠍⠻⠷⠆⠮⠽⠟⢷⡀⣀⣀⡀⢀⣀⣀⣀⡁⠈⣀⣠⣠⠙⠂⠭⠥⠭\n", - "⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢈⣛⠙⣛⣙⢓⣛⡳⣺⡿⣿⣦⠾⣛⠀⠀⠀⣀⣀\n", - "⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠈⠉⠉⠉⠉, [-1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0]))" + "171×9 Matrix{Float64}:\n", + " 0.295858 0.147929 1.67685e-11 … 2.10484e-10 1.97804e-13 0.556213\n", + " 0.295858 0.147929 5.47466e-11 6.50633e-10 1.0614e-12 0.556213\n", + " 0.295858 0.147929 1.34852e-10 1.57473e-9 3.31954e-12 0.556213\n", + " 0.295858 0.147929 3.03816e-10 3.51957e-9 8.50956e-12 0.556213\n", + " 0.295858 0.147929 6.60211e-10 7.61734e-9 1.98756e-11 0.556213\n", + " 0.295858 0.147929 1.41198e-9 … 1.6256e-8 4.42671e-11 0.556213\n", + " 0.295858 0.147929 2.99783e-9 3.44714e-8 9.61617e-11 0.556213\n", + " 0.295858 0.147929 6.3437e-9 7.28822e-8 2.06229e-10 0.556213\n", + " 0.295858 0.147929 1.34051e-8 1.53871e-7 4.39757e-10 0.556213\n", + " 0.295857 0.147928 2.83176e-8 3.24585e-7 9.37095e-10 0.556213\n", + " 0.295856 0.147928 5.9854e-8 … 6.84178e-7 2.00602e-9 0.556213\n", + " 0.295854 0.147926 1.26739e-7 1.44052e-6 4.34759e-9 0.556214\n", + " 0.29585 0.147923 2.69458e-7 3.02627e-6 9.66824e-9 0.556215\n", + " ⋮ ⋱ \n", + " 0.0356073 0.0126109 0.00388351 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.00388351 … 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.00388351 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.00388351 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.00388351 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.0038835 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.0038835 … 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.0038835 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.0038835 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.0038835 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.0038835 4.74015e-6 3.45191e-7 0.631019\n", + " 0.0356073 0.0126109 0.0038835 … 4.74015e-6 3.45191e-7 0.631019" ] }, - "execution_count": 5, + "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "gas_ct = ct.Solution(\"gri30.xml\")\n", - "gas = CreateSolution(\"../mechanism/gri30.yaml\")" + "ts= cantera_data[:, 1]\n", + "u = cantera_data[:,2:end]\n", + "T = cantera_data[:,end]\n", + "Y = cantera_data[:,2:end-1]" ] }, { "cell_type": "code", - "execution_count": 6, - "id": "b59c86ef", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "53\n", - "325" - ] - } - ], - "source": [ - "ns = gas.n_species\n", - "print(ns);\n", - "nr = gas.n_reactions\n", - "print(\"\\n\",nr);" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "bb5a6252", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "53-element Vector{String}:\n", - " \"H2\"\n", - " \"H\"\n", - " \"O\"\n", - " \"O2\"\n", - " \"OH\"\n", - " \"H2O\"\n", - " \"HO2\"\n", - " \"H2O2\"\n", - " \"C\"\n", - " \"CH\"\n", - " \"CH2\"\n", - " \"CH2(S)\"\n", - " \"CH3\"\n", - " ⋮\n", - " \"H2CN\"\n", - " \"HCNN\"\n", - " \"HCNO\"\n", - " \"HOCN\"\n", - " \"HNCO\"\n", - " \"NCO\"\n", - " \"N2\"\n", - " \"AR\"\n", - " \"C3H7\"\n", - " \"C3H8\"\n", - " \"CH2CHO\"\n", - " \"CH3CHO\"" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "species_arr = gas.species_names" - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "id": "2205b286", - "metadata": {}, - "outputs": [], - "source": [ - "#get equivalence ratio\n", - "gas_ct.TP = 900, ct.one_atm\n", - "phi = 1\n", - "gas_ct.set_equivalence_ratio(phi, \"CH4\",\"O2:1, N2:3.76\")" - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "id": "91ef853c", - "metadata": {}, - "outputs": [ - { - "ename": "LoadError", - "evalue": "UndefVarError: find not defined", - "output_type": "error", - "traceback": [ - "UndefVarError: find not defined", - "", - "Stacktrace:", - " [1] top-level scope", - " @ In[21]:2", - " [2] eval", - " @ ./boot.jl:360 [inlined]", - " [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)", - " @ Base ./loading.jl:1094" - ] - } - ], - "source": [ - "specs = gas_ct.X\n", - "find(x->(x>0),specs)" - ] - }, - { - "cell_type": "code", - "execution_count": 46, - "id": "c8b39202", - "metadata": {}, - "outputs": [], - "source": [ - "Y0 = zeros(ns)\n", - "Y0[species_index(gas, \"CH4\")] = 0.05\n", - "Y0[species_index(gas, \"O2\")] = 0.22\n", - "Y0[species_index(gas,\"N2\")] = 0.73\n", - "T0 = 1300 #K\n", - "P = one_atm\n", - "u0 = vcat(Y0, T0);" - ] - }, - { - "cell_type": "code", - "execution_count": 47, - "id": "9458bb24", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "dudt! (generic function with 1 method)" - ] - }, - "execution_count": 47, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "@inbounds function dudt!(du, u, p, t)\n", - " T = u[end]\n", - " Y = @view(u[1:ns])\n", - " mean_MW = 1. / dot(Y, 1 ./ gas.MW)\n", - " ρ_mass = P / R / T * mean_MW\n", - " X = Y2X(gas, Y, mean_MW)\n", - " C = Y2C(gas, Y, ρ_mass)\n", - " cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)\n", - " h_mole = get_H(gas, T, Y, X)\n", - " S0 = get_S(gas, T, P, X)\n", - " wdot = wdot_func(gas.reaction, T, C, S0, h_mole)\n", - " Ydot = wdot / ρ_mass .* gas.MW\n", - " Tdot = -dot(h_mole, wdot) / ρ_mass / cp_mass\n", - " du .= vcat(Ydot, Tdot)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 48, - "id": "a10b1e0d", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[0.0007372047769254919, 0.00015783813441258581, 0.0024684764564995667, 0.030240894023154722, 0.009514030316890703, 0.09925295177222965, 5.999200322512603e-6, 2.745199996984539e-7, 1.1655762192849003e-14, 1.1195793281653313e-15, 9.169500260625093e-16, 7.417468071139761e-17, 1.1526296348853825e-15, 1.6559665608689453e-16, 0.02639335196179212, 0.09569054288713183, 1.3926994094453778e-8, 1.1347285657166752e-10, 1.3074647243442918e-15, 2.6681749783590557e-17, 3.8453520809023416e-17, 5.6420499565209635e-21, 1.4785698813585826e-19, 3.2237756053507546e-24, 7.281778600979616e-25, 1.6873468039224525e-29, 3.6881881095716126e-31, 1.2623809037968665e-17, 4.2550975151592306e-18, 1.2943563840775742e-20, 3.8472657287692335e-7, 4.0801441084301696e-8, 7.14521858643824e-9, 4.924599456330184e-9, 1.0282015957786918e-8, 0.010380621594693025, 4.3757920969141466e-6, 8.035663360454821e-7, 3.6552718100736904e-7, 7.382181252022507e-12, 2.950729500815938e-10, 3.6453197064588955e-16, 4.785953872629456e-19, 8.234820493660357e-15, 1.9286533429595644e-11, 2.8854025065189156e-9, 5.447650914165151e-10, 0.7251518037980842, 0.0, 1.2457688091893588e-42, 3.1329280487968883e-43, 4.898003839005981e-23, 7.896716066200371e-25, 2625.321240165389]" - ] - } - ], - "source": [ - "tspan = [0.0, 0.07];\n", - "prob = ODEProblem(dudt!, u0, tspan);\n", - "sol = solve(prob, TRBDF2(), reltol=1e-6, abstol=1e-9);\n", - "print(sol[end])" - ] - }, - { - "cell_type": "code", - "execution_count": 49, - "id": "6dcc99bf", + "execution_count": 32, + "id": "0ae9d944", "metadata": {}, "outputs": [ { @@ -279,299 +134,290 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", - " \n", + " \n", + " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" - ] - }, - "execution_count": 49, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "#plot stuff\n", - "plt1 = plot(sol.t, sol[species_index(gas, \"CH4\"), :], lw=2, label=\"CH4 mass fraction\");\n", - "plot(plt1)" - ] - }, - { - "cell_type": "code", - "execution_count": 50, - "id": "743f60a0", - "metadata": {}, - "outputs": [ - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", "\n", - " \n", - " \n", + " \n", + " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" + "\n" ] }, - "execution_count": 50, + "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "plt2 = plot(sol.t, sol[end,:],lw=2, label =\"Temperature\");\n", - "plot(plt2)" + "plt1 = plot(ts,T,lw=2,label=\"Temperature\",legend=:bottomright,linecolor = \"red\")\n", + "plt2 = plot(ts,Y[:,1],lw=2,label=\"H2\",legend =:best,linecolor=\"orange\")\n", + "plot(plt1,plt2)" ] }, { "cell_type": "code", - "execution_count": 72, - "id": "067127ea", + "execution_count": 36, + "id": "d815e223", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - " 0.230798 seconds (405.00 k allocations: 25.309 MiB, 16.89% gc time, 99.58% compilation time)\n" + "index 1 => H2\n", + "index 2 => O2\n", + "index 3 => O\n", + "index 4 => OH\n", + "index 5 => H2O\n", + "index 6 => H\n", + "index 7 => HO2\n", + "index 8 => H2O2\n", + "index 9 => N2\n" ] - }, + } + ], + "source": [ + "gas = CreateSolution(\"../mechanism/LiDryer.yaml\")\n", + "P = one_atm\n", + "T = 1000\n", + "ns = gas.n_species\n", + "nr = gas.n_reactions\n", + "spec_names = gas.species_names\n", + "for i in 1:ns\n", + " print(\"index \", i,\" => \",spec_names[i],\"\\n\")\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "3a111882", + "metadata": {}, + "outputs": [ { "data": { "text/plain": [ - "54×54 Matrix{Float64}:\n", - " -1.52695e6 7.13195e6 … 21485.2 0.824522\n", - " 7.63502e5 -4.49871e6 -10742.5 -0.584406\n", - " -3.05046e6 2.3922e7 -28991.5 -5.04239\n", - " -5439.7 -1.94202e7 -1416.48 -3.77348\n", - " -6.39097e6 4.54834e7 -2.15116e5 24.7637\n", - " 1.0205e7 -4.76365e7 … 2.60533e5 -15.4171\n", - " 5391.86 19308.6 1455.83 0.150413\n", - " -36.2059 -1744.06 -6.56466 0.210893\n", - " -4.62871e-6 0.000280186 -1.68825e-18 1.31532e-9\n", - " 8.09969e-7 5.60372e-5 5.85536e-12 1.36516e-10\n", - " 4.41289e-6 2.61133e-6 … 2.77507e-12 2.2412e-11\n", - " -2.26432e-7 5.40157e-6 -1.70388e-10 9.57e-12\n", - " 2.85661e-7 1.08957e-5 3.11854e5 7.702e-12\n", - " ⋮ ⋱ \n", - " -2.08883e-10 -2.43301e-9 -9.57743e-12 -2.04416e-14\n", - " -2.88522e-17 -2.72726e-6 -1.32036e-18 4.8977e-12\n", - " -1.95895e-13 -0.154448 -8.96475e-15 8.09177e-8\n", - " 0.222192 -3.58519 … -1.56636e-5 4.71834e-7\n", - " -0.219206 -2.72015 -8.62846e-5 2.93628e-6\n", - " 55.0283 -31106.4 3.12279 -0.00718459\n", - " 0.0 0.0 0.0 0.0\n", - " -7.04006e-35 -9.53325e-33 2.72048e-36 6.30171e-39\n", - " 1.58833e-34 -6.73867e-33 … 1.72678e-36 -3.24201e-39\n", - " -7.55654e-15 1.09565e-11 2.68386e5 -4.46754e-18\n", - " 1.08791e-15 -2.34898e-15 -1.18843e6 7.88067e-22\n", - " 2.09748e10 -1.38905e11 1.24381e9 -48623.4" + "dudt (generic function with 1 method)" ] }, - "execution_count": 72, + "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "#compute the jacobian \n", "function dudt(u)\n", " T = u[end]\n", " Y = @view(u[1:ns])\n", @@ -586,245 +432,601 @@ " Ydot = wdot / ρ_mass .* gas.MW\n", " Tdot = -dot(h_mole, wdot) / ρ_mass / cp_mass\n", " du = vcat(Ydot, Tdot)\n", - "end\n", - "@time J_w = ForwardDiff.jacobian(dudt,sol[end])" + "end" ] }, { "cell_type": "code", - "execution_count": 73, - "id": "33e12125", + "execution_count": 41, + "id": "e65c60f8", "metadata": {}, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.000169 seconds (102 allocations: 45.594 KiB)\n" + ] + }, { "data": { "text/plain": [ - "54-element Vector{Float64}:\n", - " -5.1402940277116776e8\n", - " -5.13152406186228e8\n", - " -1.047488154189249e8\n", - " -4.016960521435924e7\n", - " -3.712423735420916e7\n", - " -3.656225824242137e7\n", - " -3.465295356192211e7\n", - " -3.458235019079888e7\n", - " -3.222971058255094e7\n", - " -3.218474998638545e7\n", - " -2.1424526246316772e7\n", - " -1.6560164920316156e7\n", - " -1.5350392042396849e7\n", - " ⋮\n", - " -376970.8910503711\n", - " -263703.93432325067\n", - " -159718.90651112387\n", - " -123941.40921549303\n", - " -6630.883336015616\n", - " -304.2359055369181\n", - " -3.951849578020171e-8\n", - " -5.965913016854841e-9\n", - " -1.1194473550679741e-10\n", - " 0.0\n", - " 2.555719020746241e-8\n", - " 1.1448487553316728e-7" + "10×10 Matrix{Float64}:\n", + " -8.29518e-6 -1.4236e-5 -2.90133e5 … 3.57898e-6 -1.15115e-7\n", + " 0.000601601 -0.00173505 0.0300868 6.81274e-5 -1.28866e-6\n", + " -0.000151941 0.000309684 -2.3025e6 -1.54301e-6 1.32556e-7\n", + " -0.000161496 0.000329189 2.44756e6 -1.64839e-6 3.54018e-7\n", + " 4.19916e-5 -1.58617e-5 -3.38605e-5 -1.81174e-5 1.93986e-7\n", + " 2.372e-5 -3.90892e-5 1.45066e5 … 8.88356e-8 5.18766e-8\n", + " -0.000345645 0.00116539 -0.030487 -5.04583e-5 6.69769e-7\n", + " 6.50107e-8 -2.4557e-8 -1.1923e-6 -2.80493e-8 1.57299e-9\n", + " 0.0 0.0 0.0 0.0 0.0\n", + " -0.718704 0.446667 -1.65482e8 -0.0725448 -0.00226572" ] }, - "execution_count": 73, + "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "eig_vals = eigvals(J_w)" + "#compute Jacobian for one\n", + "@time J_w = ForwardDiff.jacobian(dudt,u[1,:])" ] }, { "cell_type": "code", - "execution_count": 78, - "id": "488bb2c0", + "execution_count": 46, + "id": "463b5616", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "2-element Vector{Int64}:\n", - " 53\n", - " 54" + "10-element Vector{Float64}:\n", + " -1.3649840528402898e7\n", + " -2.339248624761963e6\n", + " -1281.7349093373377\n", + " -1014.242369146887\n", + " -0.0030068068408325276\n", + " -2.3356973247393255e-6\n", + " -4.567161198433734e-10\n", + " 0.0\n", + " 4.915698386826457e-9\n", + " 10526.224724477977" ] }, - "execution_count": 78, + "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "#now extract all the positive eigen values: \n", - "findall(x->(x > 0), eig_vals) " + " lambda = eigvals(J_w)" ] }, { "cell_type": "code", - "execution_count": 79, - "id": "19874e03", + "execution_count": 48, + "id": "52451f89", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "54×54 Matrix{Float64}:\n", - " 2.70714e-8 -7.1509e-7 -2.51146e-6 … 3.68468e-6 2.64177e-6\n", - " 1.72087e-9 7.05283e-5 2.57949e-6 1.02975e-6 9.2112e-7\n", - " -2.65337e-7 -7.3539e-6 3.41717e-6 9.43655e-6 1.23394e-5\n", - " 1.77707e-5 1.80682e-6 -3.41017e-5 -4.11591e-5 2.88098e-5\n", - " 0.000175692 9.22861e-6 8.56997e-5 2.33229e-5 2.75997e-5\n", - " -0.000274375 -1.75474e-5 -8.70614e-5 … -1.07149e-5 -3.82003e-5\n", - " -0.000537388 -5.1194e-5 -4.49305e-7 7.61265e-9 1.7135e-8\n", - " 0.000618398 6.21162e-5 -2.38254e-8 3.1115e-10 5.47463e-10\n", - " 4.26755e-18 -4.03337e-17 6.37716e-7 2.2017e-16 2.0083e-16\n", - " -8.61843e-18 -8.09984e-18 -1.97865e-5 2.15168e-17 1.88548e-17\n", - " -6.28377e-19 -1.0581e-19 0.000285173 … 1.68573e-17 1.3881e-17\n", - " -2.15674e-19 -1.08432e-18 -0.000372372 1.40081e-18 1.16429e-18\n", - " -4.37923e-19 -2.01582e-18 7.11077e-5 1.87427e-17 1.41786e-17\n", - " ⋮ ⋱ \n", - " -4.35165e-22 5.23838e-22 -1.60644e-10 7.95851e-21 6.79772e-21\n", - " 5.54854e-21 3.90167e-19 2.17009e-7 9.22833e-17 8.15808e-17\n", - " -7.22635e-16 2.17072e-14 -7.47314e-10 1.49925e-13 1.25346e-13\n", - " 1.67162e-14 -1.21096e-11 5.10093e-7 … 1.74148e-11 1.37376e-11\n", - " 2.14649e-14 -7.77908e-11 4.40618e-7 4.34532e-12 4.04681e-12\n", - " -6.00393e-7 0.00195975 -9.3441e-8 -1.7715e-6 -1.00397e-5\n", + "(10526.224724477977, 10)" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(lambda_exp,index) = findmax(lambda)" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "id": "db4f2202", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "10×10 Matrix{Float64}:\n", + " 0.000170858 -0.000361498 5.88285e-5 … 0.000574718 -9.19066e-5\n", + " -8.17434e-6 0.000102412 0.000964555 0.00383013 -0.00116079\n", + " 2.00963e-6 -0.00133222 1.09587e-7 4.07793e-14 1.07975e-6\n", + " 0.00144242 -0.000284022 5.3767e-8 -1.15539e-14 3.89019e-7\n", + " -0.00152641 0.00175381 -0.000605946 -0.00533396 0.000533834\n", + " -8.56866e-5 0.000183975 9.47585e-7 … -3.3314e-12 9.38443e-6\n", + " 4.98554e-6 -6.24748e-5 -0.00114298 1.46373e-10 0.000668433\n", + " -2.33532e-10 1.70778e-8 0.000724427 3.37723e-10 3.95748e-5\n", " 0.0 0.0 0.0 0.0 0.0\n", - " -1.60856e-46 1.27799e-45 -1.33364e-34 0.0 9.10949e-22\n", - " 1.67265e-47 9.5001e-46 -2.28985e-36 … 2.48673e-21 2.04963e-21\n", - " -1.70227e-25 -1.63193e-24 5.22568e-14 1.61638e-21 5.76934e-21\n", - " 2.00756e-29 3.84882e-28 -9.66608e-18 -8.08188e-22 4.40292e-21\n", - " -1.0 0.999996 1.0 1.0 1.0" + " -0.999998 0.999997 -0.999998 -0.999978 0.999999" + ] + }, + "execution_count": 49, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "right_vecs = eigvecs(J_w)" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "id": "bb367a73", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "10-element Vector{Float64}:\n", + " -9.190658370980618e-5\n", + " -0.0011607883812806672\n", + " 1.0797534510935694e-6\n", + " 3.8901928168201305e-7\n", + " 0.0005338340983246179\n", + " 9.384430492797351e-6\n", + " 0.0006684328938853956\n", + " 3.9574769554508826e-5\n", + " 0.0\n", + " 0.9999989553427485" + ] + }, + "execution_count": 50, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a_e = right_vecs[:,index]" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "id": "3eb45354", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "10×10 Matrix{Float64}:\n", + " 6.5555e-9 -1.3354e-8 … 6.68957e-11 -1.64749e-11\n", + " -4.85195e-8 9.85554e-8 -4.87198e-10 3.91395e-11\n", + " 6.88239e-5 -0.000670348 -7.99384e-5 -1.37255e-5\n", + " -6.69216e-5 0.000556992 5.81798e-5 1.0395e-5\n", + " 155.349 273.753 2.27638 1.09586\n", + " 15276.0 1003.7 … 1198.23 0.0956694\n", + " -1077.86 26.3587 -49.5947 -7.73812e-7\n", + " 0.0 0.0 1.94056 0.0\n", + " -14042.9 -756.306 -1144.79 0.000192246\n", + " -4.85132e-5 0.000187584 -2.57503e-6 7.97522e-7" ] }, - "execution_count": 79, + "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "eigvecs(J_w)" + "left_vecs = inv(eigvecs(J_w))" ] }, { "cell_type": "code", - "execution_count": 82, - "id": "927e9dc8", + "execution_count": 55, + "id": "d37fde01", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "54×54 Diagonal{Float64, Vector{Float64}}:\n", - " -5.14029e8 ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ -5.13152e8 ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ -1.04749e8 ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋮ ⋱ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ 2.55572e-8 ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ 1.14485e-7" + "10-element Vector{Float64}:\n", + " -1.6474922337819664e-11\n", + " 3.913951714812347e-11\n", + " -1.372545613520515e-5\n", + " 1.0394973473992764e-5\n", + " 1.0958552265203763\n", + " 0.09566936427761297\n", + " -7.73811907217592e-7\n", + " 0.0\n", + " 0.000192246436911514\n", + " 7.975222917622204e-7" ] }, - "execution_count": 82, + "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "F = eigen(J_w)\n", - "Q = eigvecs(F) # right eigenvectors \n", - "QL = inv(eigvecs(F)) # left eigenvectors \n", - "Λ = Diagonal(eigvals(F))" + "b_e = left_vecs[:,index]" ] }, { "cell_type": "code", - "execution_count": 86, - "id": "1a49e810", + "execution_count": 61, + "id": "f5b239f7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "54.000000000000014" + "10-element Vector{Float64}:\n", + " 1.5141538289533787e-15\n", + " -4.543269675447716e-14\n", + " -1.4820108629821165e-11\n", + " 4.043845113956244e-12\n", + " 0.000585004886743825\n", + " 8.978024993533688e-7\n", + " -5.172413324644322e-10\n", + " 0.0\n", + " 0.0\n", + " 7.975214586247751e-7" ] }, - "execution_count": 86, + "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "EP = Diagonal(Q*QL)\n", - "sum(EP)" + "EP = a_e.*b_e" ] }, { "cell_type": "code", - "execution_count": 87, - "id": "6db2f983", + "execution_count": 64, + "id": "cb627116", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "54×54 Diagonal{Float64, Vector{Float64}}:\n", - " 0.0185185 ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ 0.0185185 ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ 0.0185185 ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋮ ⋱ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ 0.0185185 ⋅ ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ 0.0185185 ⋅ \n", - " ⋅ ⋅ ⋅ ⋅ ⋅ 0.0185185" + "10-element Vector{Float64}:\n", + " 2.580798786423959e-12\n", + " 7.743773875932433e-11\n", + " 2.5260127230897995e-8\n", + " 6.892529915404038e-9\n", + " 0.9971113059259952\n", + " 0.0015302590506152027\n", + " 8.816117474901688e-7\n", + " 0.0\n", + " 0.0\n", + " 0.0013593350775915512" + ] + }, + "execution_count": 64, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "EI = abs.(EP)/sum(EP)" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "id": "f9b3d7a6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(171, 10)" + ] + }, + "execution_count": 69, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nrow,ncol = size(u)" + ] + }, + { + "cell_type": "code", + "execution_count": 88, + "id": "d25a295f", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "execution_count": 88, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "EP = []\n", + "EI = []\n", + "lambda_exp_arr = zeros(nrow)\n", + "for i = 1:nrow\n", + " u_node = u[i,:]\n", + " J_w = ForwardDiff.jacobian(dudt,u_node)\n", + " lambda = eigvals(J_w)\n", + " (lambda_exp,index) = findmax(real(lambda))\n", + " if(all(<=(0),real(lambda)))\n", + " (lambda_exp,index) = findmin(real(lambda))\n", + " end\n", + " lambda_exp_arr[i] = lambda_exp\n", + " right_vecs = eigvecs(J_w)\n", + " a_e = right_vecs[:,index]\n", + " left_vecs = inv(eigvecs(J_w))\n", + " b_e = left_vecs[:,index]\n", + " EP = [EP,a_e.*b_e]\n", + " EI = [EI,abs.(a_e.*b_e)/sum(a_e.*b_e)]\n", + "end\n", + "plot(ts,lambda_exp_arr)" + ] + }, + { + "cell_type": "code", + "execution_count": 94, + "id": "fb441897", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" ] }, - "execution_count": 87, + "execution_count": 94, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "EI = Diagonal(Q*QL)/sum(EP)" + "index = findall(x->(x>0),lambda_exp_arr)\n", + "t = zeros(length(index));\n", + "l = zeros(length(index));\n", + "for i =1:length(index)\n", + " t[i] = ts[index[i]]\n", + " l[i] = lambda_exp_arr[index[i]]\n", + "end\n", + "plot(t,l)" ] } ], diff --git a/Isobaric reactor H2 autoignition.ipynb b/Isobaric reactor H2 autoignition.ipynb new file mode 100644 index 0000000..e5de25a --- /dev/null +++ b/Isobaric reactor H2 autoignition.ipynb @@ -0,0 +1,1460 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 9, + "id": "1ef5f90d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Initializing reactor network.\n", + "Reactor 0: 11 variables.\n", + " 0 sensitivity params.\n", + "Number of equations: 11\n", + "Maximum time step: 0\n", + "t [s] T [K] P [Pa] u [J/kg] \n", + " 1.000e-05 1001.000 101325.000 627729.037756\n", + " 2.000e-05 1001.000 101325.000 627729.037334\n", + " 3.000e-05 1001.000 101325.000 627729.035564\n", + " 4.000e-05 1001.000 101325.000 627729.030951\n", + " 5.000e-05 1001.000 101325.000 627729.020340\n", + " 6.000e-05 1001.000 101325.000 627728.997077\n", + " 7.000e-05 1001.000 101325.000 627728.947117\n", + " 8.000e-05 1001.001 101325.000 627728.840804\n", + " 9.000e-05 1001.001 101325.000 627728.615444\n", + " 1.000e-04 1001.003 101325.000 627728.138141\n", + " 1.100e-04 1001.005 101325.000 627727.125684\n", + " 1.200e-04 1001.012 101325.000 627724.967729\n", + " 1.300e-04 1001.025 101325.000 627720.318786\n", + " 1.400e-04 1001.054 101325.000 627710.078787\n", + " 1.500e-04 1001.121 101325.000 627686.511834\n", + " 1.600e-04 1001.286 101325.000 627627.715416\n", + " 1.700e-04 1001.746 101325.000 627460.902190\n", + " 1.800e-04 1003.234 101325.000 626912.222255\n", + " 1.900e-04 1008.368 101325.000 625000.947293\n", + " 2.000e-04 1025.800 101325.000 618508.642613\n", + " Limiting global state vector component 1 (dt = 1e-05): 65.3607 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 21.7315 > 20\n", + " 2.025e-04 1034.853 101325.000 615141.344949\n", + " Limiting global state vector component 1 (dt = 1e-05): 95.8519 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 30.8315 > 20\n", + " 2.050e-04 1047.551 101325.000 610423.578960\n", + " Limiting global state vector component 1 (dt = 1e-05): 146.424 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 44.8852 > 20\n", + " 2.075e-04 1065.728 101325.000 603680.965945\n", + " Limiting global state vector component 1 (dt = 1e-05): 238.702 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 67.9702 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 26.815 > 20\n", + " 2.088e-04 1077.777 101325.000 599218.274728\n", + " Limiting global state vector component 1 (dt = 1e-05): 316.097 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 85.4499 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 33.0973 > 20\n", + " 2.100e-04 1092.548 101325.000 593755.232848\n", + " Limiting global state vector component 1 (dt = 1e-05): 428.544 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 109.372 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 41.4314 > 20\n", + " 2.113e-04 1110.882 101325.000 586985.970807\n", + " Limiting global state vector component 1 (dt = 1e-05): 602.411 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 143.283 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 52.7847 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 23.1108 > 20\n", + " 2.119e-04 1121.745 101325.000 582981.273789\n", + " Limiting global state vector component 1 (dt = 1e-05): 720.839 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 165.521 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 60.0406 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 26.1272 > 20\n", + " 2.125e-04 1133.994 101325.000 578471.280591\n", + " Limiting global state vector component 1 (dt = 1e-05): 862.393 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 192.287 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 68.6806 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 29.6917 > 20\n", + " 2.131e-04 1147.873 101325.000 573367.851198\n", + " Limiting global state vector component 1 (dt = 1e-05): 1037.86 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 224.852 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 79.0433 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 33.934 > 20\n", + " 2.138e-04 1163.686 101325.000 567562.740425\n", + " Limiting global state vector component 1 (dt = 1e-05): 1237.73 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 263.907 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 91.5093 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 39.0143 > 20\n", + " 2.144e-04 1181.808 101325.000 560922.744261\n", + " Limiting global state vector component 1 (dt = 1e-05): 1435.55 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 309.518 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 106.489 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 45.1242 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 20.8928 > 20\n", + " 2.147e-04 1191.875 101325.000 557239.922243\n", + " Limiting global state vector component 1 (dt = 1e-05): 1512.58 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 334.246 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 115.033 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 48.6313 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 22.4838 > 20\n", + " 2.150e-04 1202.700 101325.000 553284.493883\n", + " Limiting global state vector component 1 (dt = 1e-05): 1553.69 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 359.482 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 124.31 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 52.4749 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 24.2299 > 20\n", + " 2.153e-04 1214.358 101325.000 549030.624234\n", + " Limiting global state vector component 1 (dt = 1e-05): 1533.61 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 384.155 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 134.306 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 56.6777 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 26.1446 > 20\n", + " 2.156e-04 1226.930 101325.000 544450.264233\n", + " Limiting global state vector component 1 (dt = 1e-05): 1407.84 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 406.337 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 144.95 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 61.2556 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 28.2402 > 20\n", + " 2.159e-04 1240.503 101325.000 539513.355591\n", + " Limiting global state vector component 1 (dt = 1e-05): 1100.3 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 422.624 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 156.062 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 66.21 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 30.5255 > 20\n", + " 2.162e-04 1255.170 101325.000 534188.302332\n", + " Limiting global state vector component 1 (dt = 1e-05): 612.248 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 431.745 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 167.471 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 71.5274 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 33.004 > 20\n", + " 2.166e-04 1271.028 101325.000 528442.840312\n", + " Limiting global state vector component 1 (dt = 1e-05): 153.495 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 428.521 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 178.762 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 77.1596 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 35.6699 > 20\n", + " 2.169e-04 1288.173 101325.000 522245.474757\n", + " Limiting global state vector component 1 (dt = 1e-05): 1252.64 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 408.421 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 189.406 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 83.0185 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 38.5031 > 20\n", + " 2.172e-04 1306.697 101325.000 515567.685630\n", + " Limiting global state vector component 1 (dt = 1e-05): 2689.02 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 367.979 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 198.753 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 88.9627 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 41.4635 > 20\n", + " 2.175e-04 1326.676 101325.000 508387.099957\n", + " Limiting global state vector component 1 (dt = 1e-05): 4354.84 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 307.153 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 206.11 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 94.7896 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 44.4849 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 21.4842 > 20\n", + " 2.177e-04 1337.228 101325.000 504603.975172\n", + " Limiting global state vector component 1 (dt = 1e-05): 5187.36 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 270.964 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 208.861 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 97.5795 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 45.9892 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 22.244 > 20\n", + " 2.178e-04 1348.160 101325.000 500691.757009\n", + " Limiting global state vector component 1 (dt = 1e-05): 5977.39 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 232.003 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 210.896 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 100.237 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 47.4693 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 23.0001 > 20\n", + " 2.180e-04 1359.472 101325.000 496651.445412\n", + " Limiting global state vector component 1 (dt = 1e-05): 6489.14 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 198.017 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 212.391 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 102.729 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 48.9083 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 23.7448 > 20\n", + " 2.181e-04 1371.160 101325.000 492485.388175\n", + " Limiting global state vector component 1 (dt = 1e-05): 6742.94 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 168.616 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 213.255 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 105.013 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 50.2868 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 24.4693 > 20\n", + " 2.183e-04 1383.217 101325.000 488197.462722\n", + " Limiting global state vector component 1 (dt = 1e-05): 6433.01 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 154.284 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 213.79 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 107.059 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 51.5845 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 25.1639 > 20\n", + " 2.184e-04 1395.629 101325.000 483793.239421\n", + " Limiting global state vector component 1 (dt = 1e-05): 6187.7 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 136.257 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 213.389 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 108.809 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 52.7795 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 25.8184 > 20\n", + " 2.186e-04 1408.381 101325.000 479280.111351\n", + " Limiting global state vector component 1 (dt = 1e-05): 4740.58 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 155.535 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 213.365 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 110.269 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 53.8517 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 26.4217 > 20\n", + " 2.187e-04 1421.448 101325.000 474667.373565\n", + " Limiting global state vector component 1 (dt = 1e-05): 2640.95 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 197.13 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 213.355 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 111.404 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 54.7807 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 26.9632 > 20\n", + " 2.189e-04 1434.802 101325.000 469966.234156\n", + " Limiting global state vector component 1 (dt = 1e-05): 423.973 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 245.283 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 212.924 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 112.175 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 55.5472 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 27.4324 > 20\n", + " 2.191e-04 1448.411 101325.000 465189.741212\n", + " Limiting global state vector component 1 (dt = 1e-05): 2172.22 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 310.18 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 212.547 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 112.591 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 56.1367 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 27.8198 > 20\n", + " 2.192e-04 1462.235 101325.000 460352.613599\n", + " Limiting global state vector component 1 (dt = 1e-05): 4974.03 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 387.712 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 212.241 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 112.652 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 56.5381 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 28.1175 > 20\n", + " 2.194e-04 1476.231 101325.000 455470.970534\n", + " Limiting global state vector component 1 (dt = 1e-05): 7734.93 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 470.923 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 211.93 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 112.365 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 56.745 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 28.3198 > 20\n", + " 2.195e-04 1490.352 101325.000 450561.964316\n", + " Limiting global state vector component 1 (dt = 1e-05): 10501 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 561.647 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 211.813 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 111.759 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 56.7571 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 28.4233 > 20\n", + " 2.197e-04 1504.551 101325.000 445643.331536\n", + " Limiting global state vector component 1 (dt = 1e-05): 12144.9 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 622.836 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 210.744 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 110.82 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 56.5773 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 28.4274 > 20\n", + " 2.198e-04 1518.776 101325.000 440732.889006\n", + " Limiting global state vector component 1 (dt = 1e-05): 12958.5 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 662.604 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 209.044 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 109.591 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 56.216 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 28.3345 > 20\n", + " 2.200e-04 1532.978 101325.000 435848.009488\n", + " Limiting global state vector component 1 (dt = 1e-05): 12847 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 675.936 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 206.545 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 108.104 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 55.6873 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 28.1496 > 20\n", + " 2.202e-04 1547.110 101325.000 431005.117142\n", + " Limiting global state vector component 1 (dt = 1e-05): 11840.5 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 661.535 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 203.167 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 106.39 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 55.0093 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 27.8801 > 20\n", + " 2.203e-04 1561.128 101325.000 426219.242033\n", + " Limiting global state vector component 1 (dt = 1e-05): 9300.93 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 596.393 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 198.069 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 104.456 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 54.2007 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 27.5352 > 20\n", + " 2.205e-04 1574.990 101325.000 421503.666782\n", + " Limiting global state vector component 1 (dt = 1e-05): 6975.04 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 535.197 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 193.019 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 102.398 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 53.2858 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 27.1257 > 20\n", + " 2.206e-04 1588.663 101325.000 416869.687530\n", + " Limiting global state vector component 1 (dt = 1e-05): 4508.27 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 465.189 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 187.537 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 100.23 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 52.2854 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 26.6626 > 20\n", + " 2.208e-04 1602.116 101325.000 412326.497706\n", + " Limiting global state vector component 1 (dt = 1e-05): 2106.69 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 392.459 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 181.798 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 97.987 > 20\n", + " Limiting global state vector component 1 (dt = 6.25e-07): 51.2205 > 20\n", + " Limiting global state vector component 1 (dt = 3.125e-07): 26.1571 > 20\n", + " 2.209e-04 1615.325 101325.000 407881.189334\n", + " 2.309e-04 2069.410 101325.000 262435.182754\n", + " Limiting global state vector component 1 (dt = 1e-05): 173.067 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 105.153 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 57.3676 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 30.0311 > 20\n", + " 2.316e-04 2084.790 101325.000 257725.554812\n", + " Limiting global state vector component 1 (dt = 1e-05): 167.909 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 100.752 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 54.7896 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 28.6359 > 20\n", + " 2.322e-04 2099.442 101325.000 253251.179041\n", + " Limiting global state vector component 1 (dt = 1e-05): 162.346 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 96.6273 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 52.4068 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 27.351 > 20\n", + " 2.328e-04 2113.426 101325.000 248992.381685\n", + " Limiting global state vector component 1 (dt = 1e-05): 157.419 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 92.8155 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 50.2026 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 26.1646 > 20\n", + " 2.334e-04 2126.794 101325.000 244931.770655\n", + " Limiting global state vector component 1 (dt = 1e-05): 152.3 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 89.2341 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 48.1547 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 25.0657 > 20\n", + " 2.341e-04 2139.591 101325.000 241053.902501\n", + " Limiting global state vector component 1 (dt = 1e-05): 147.208 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 85.8727 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 46.248 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 24.0454 > 20\n", + " 2.347e-04 2151.860 101325.000 237345.007727\n", + " Limiting global state vector component 1 (dt = 1e-05): 142.698 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 82.748 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 44.4714 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 23.0958 > 20\n", + " 2.353e-04 2163.637 101325.000 233792.761803\n", + " Limiting global state vector component 1 (dt = 1e-05): 138.399 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 79.8164 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 42.8108 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 22.21 > 20\n", + " 2.359e-04 2174.956 101325.000 230386.093937\n", + " Limiting global state vector component 1 (dt = 1e-05): 134.558 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 77.0789 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 41.2569 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 21.3821 > 20\n", + " 2.366e-04 2185.847 101325.000 227115.026032\n", + " Limiting global state vector component 1 (dt = 1e-05): 130.577 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 74.48 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 39.7973 > 20\n", + " Limiting global state vector component 1 (dt = 1.25e-06): 20.6064 > 20\n", + " 2.372e-04 2196.338 101325.000 223970.536384\n", + " Limiting global state vector component 1 (dt = 1e-05): 126.781 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 72.0292 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 38.4251 > 20\n", + " 2.384e-04 2216.217 101325.000 218029.308276\n", + " Limiting global state vector component 1 (dt = 1e-05): 119.72 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 67.5265 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 35.9145 > 20\n", + " 2.397e-04 2234.766 101325.000 212505.357245\n", + " Limiting global state vector component 1 (dt = 1e-05): 113.249 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 63.4856 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 33.6743 > 20\n", + " 2.409e-04 2252.133 101325.000 207351.036617\n", + " Limiting global state vector component 1 (dt = 1e-05): 107.229 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 59.8351 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 31.6636 > 20\n", + " 2.422e-04 2268.442 101325.000 202526.183662\n", + " Limiting global state vector component 1 (dt = 1e-05): 101.846 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 56.5369 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 29.8505 > 20\n", + " 2.434e-04 2283.798 101325.000 197996.659176\n", + " Limiting global state vector component 1 (dt = 1e-05): 96.8224 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 53.5307 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 28.207 > 20\n", + " 2.447e-04 2298.294 101325.000 193733.223315\n", + " Limiting global state vector component 1 (dt = 1e-05): 92.1894 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 50.7842 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 26.711 > 20\n", + " 2.459e-04 2312.006 101325.000 189710.658968\n", + " Limiting global state vector component 1 (dt = 1e-05): 87.9052 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 48.2655 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 25.3438 > 20\n", + " 2.472e-04 2325.005 101325.000 185907.080618\n", + " Limiting global state vector component 1 (dt = 1e-05): 83.9357 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 45.9479 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 24.0897 > 20\n", + " 2.484e-04 2337.351 101325.000 182303.383804\n", + " Limiting global state vector component 1 (dt = 1e-05): 80.2452 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 43.8085 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 22.9354 > 20\n", + " 2.497e-04 2349.096 101325.000 178882.802449\n", + " Limiting global state vector component 1 (dt = 1e-05): 76.8372 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 41.8296 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 21.8698 > 20\n", + " 2.509e-04 2360.287 101325.000 175630.549780\n", + " Limiting global state vector component 1 (dt = 1e-05): 73.6123 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 39.9894 > 20\n", + " Limiting global state vector component 1 (dt = 2.5e-06): 20.8828 > 20\n", + " 2.522e-04 2370.966 101325.000 172533.524769\n", + " Limiting global state vector component 1 (dt = 1e-05): 70.5775 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 38.2751 > 20\n", + " 2.547e-04 2390.932 101325.000 166759.773115\n", + " Limiting global state vector component 1 (dt = 1e-05): 65.1517 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 35.1857 > 20\n", + " 2.572e-04 2409.249 101325.000 161482.235511\n", + " Limiting global state vector component 1 (dt = 1e-05): 60.323 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 32.4701 > 20\n", + " 2.597e-04 2426.122 101325.000 156636.744373\n", + " Limiting global state vector component 1 (dt = 1e-05): 56.3284 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 30.0858 > 20\n", + " 2.622e-04 2441.722 101325.000 152170.710009\n", + " Limiting global state vector component 1 (dt = 1e-05): 52.1351 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 27.92 > 20\n", + " 2.647e-04 2456.189 101325.000 148040.530821\n", + " Limiting global state vector component 1 (dt = 1e-05): 48.6427 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 25.9951 > 20\n", + " 2.672e-04 2469.643 101325.000 144209.680683\n", + " Limiting global state vector component 1 (dt = 1e-05): 45.4775 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 24.2579 > 20\n", + " 2.697e-04 2482.186 101325.000 140647.270032\n", + " Limiting global state vector component 1 (dt = 1e-05): 42.5929 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 22.6817 > 20\n", + " 2.722e-04 2493.902 101325.000 137326.949959\n", + " Limiting global state vector component 1 (dt = 1e-05): 39.9423 > 20\n", + " Limiting global state vector component 1 (dt = 5e-06): 21.2442 > 20\n", + " 2.747e-04 2504.868 101325.000 134226.065620\n", + " Limiting global state vector component 1 (dt = 1e-05): 37.5209 > 20\n", + " 2.797e-04 2524.798 101325.000 128606.627363\n", + " Limiting global state vector component 1 (dt = 1e-05): 33.2206 > 20\n", + " 2.847e-04 2542.405 101325.000 123659.686257\n", + " Limiting global state vector component 1 (dt = 1e-05): 29.5128 > 20\n", + " 2.897e-04 2558.026 101325.000 119284.625474\n", + " Limiting global state vector component 1 (dt = 1e-05): 26.3031 > 20\n", + " 2.947e-04 2571.928 101325.000 115401.393286\n", + " Limiting global state vector component 1 (dt = 1e-05): 23.4915 > 20\n", + " 2.997e-04 2584.334 101325.000 111945.053717\n", + " Limiting global state vector component 1 (dt = 1e-05): 21.0152 > 20\n", + " 3.047e-04 2595.423 101325.000 108862.026574\n", + " 3.147e-04 2614.248 101325.000 103643.204895\n", + " 3.247e-04 2629.387 101325.000 99459.507291\n", + " 3.347e-04 2641.584 101325.000 96097.518634\n", + " 3.447e-04 2651.416 101325.000 93392.794713\n", + " 3.547e-04 2659.342 101325.000 91216.054335\n", + " 3.647e-04 2665.730 101325.000 89464.365984\n", + " 3.747e-04 2670.873 101325.000 88055.164187\n", + " 3.847e-04 2675.013 101325.000 86921.952903\n", + " 3.947e-04 2678.343 101325.000 86011.078636\n", + " 4.047e-04 2681.020 101325.000 85279.227992\n", + " 4.147e-04 2683.171 101325.000 84691.443161\n", + " 4.247e-04 2684.899 101325.000 84219.525164\n", + " 4.347e-04 2686.286 101325.000 83840.744791\n", + " 4.447e-04 2687.399 101325.000 83536.795780\n", + " 4.547e-04 2688.293 101325.000 83292.945071\n", + " 4.647e-04 2689.010 101325.000 83097.343340\n", + " 4.747e-04 2689.584 101325.000 82940.465970\n", + " 4.847e-04 2690.045 101325.000 82814.661018\n", + " 4.947e-04 2690.415 101325.000 82713.783497\n", + " 5.047e-04 2690.712 101325.000 82632.900367\n", + " 5.147e-04 2690.949 101325.000 82568.052649\n", + " 5.247e-04 2691.140 101325.000 82516.063852\n", + " 5.347e-04 2691.293 101325.000 82474.385818\n", + " 5.447e-04 2691.415 101325.000 82440.974799\n", + " 5.547e-04 2691.513 101325.000 82414.191597\n", + " 5.647e-04 2691.592 101325.000 82392.721898\n", + " 5.747e-04 2691.655 101325.000 82375.511861\n", + " 5.847e-04 2691.706 101325.000 82361.716531\n", + " 5.947e-04 2691.746 101325.000 82350.658511\n", + " 6.047e-04 2691.779 101325.000 82341.794730\n", + " 6.147e-04 2691.805 101325.000 82334.689836\n", + " 6.247e-04 2691.826 101325.000 82328.994836\n", + " 6.347e-04 2691.842 101325.000 82324.429972\n", + " 6.447e-04 2691.856 101325.000 82320.770988\n", + " 6.547e-04 2691.867 101325.000 82317.838138\n", + " 6.647e-04 2691.875 101325.000 82315.487358\n", + " 6.747e-04 2691.882 101325.000 82313.603093\n", + " 6.847e-04 2691.888 101325.000 82312.092733\n", + " 6.947e-04 2691.892 101325.000 82310.882108\n", + " 7.047e-04 2691.896 101325.000 82309.911737\n", + " 7.147e-04 2691.899 101325.000 82309.133937\n", + " 7.247e-04 2691.901 101325.000 82308.510491\n", + " 7.347e-04 2691.903 101325.000 82308.010771\n", + " 7.447e-04 2691.904 101325.000 82307.610223\n", + " 7.547e-04 2691.905 101325.000 82307.289163\n", + " 7.647e-04 2691.906 101325.000 82307.031820\n", + " 7.747e-04 2691.907 101325.000 82306.825554\n", + " 7.847e-04 2691.908 101325.000 82306.660241\n", + " 7.947e-04 2691.908 101325.000 82306.527748\n", + " 8.047e-04 2691.908 101325.000 82306.421544\n", + " 8.147e-04 2691.909 101325.000 82306.336403\n", + " 8.247e-04 2691.909 101325.000 82306.268147\n", + " 8.347e-04 2691.909 101325.000 82306.213434\n", + " 8.447e-04 2691.909 101325.000 82306.169583\n", + " 8.547e-04 2691.910 101325.000 82306.134434\n", + " 8.647e-04 2691.910 101325.000 82306.106258\n", + " 8.747e-04 2691.910 101325.000 82306.083670\n", + " 8.847e-04 2691.910 101325.000 82306.065564\n", + " 8.947e-04 2691.910 101325.000 82306.051056\n", + " 9.047e-04 2691.910 101325.000 82306.039434\n", + " 9.147e-04 2691.910 101325.000 82306.030137\n", + " 9.247e-04 2691.910 101325.000 82306.022681\n", + " 9.347e-04 2691.910 101325.000 82306.016718\n", + " 9.447e-04 2691.910 101325.000 82306.011919\n", + " 9.547e-04 2691.910 101325.000 82306.008069\n", + " 9.647e-04 2691.910 101325.000 82306.004970\n", + " 9.747e-04 2691.910 101325.000 82306.002482\n", + " 9.847e-04 2691.910 101325.000 82306.000488\n", + " 9.947e-04 2691.910 101325.000 82305.998892\n", + " 1.005e-03 2691.910 101325.000 82305.997631\n" + ] + } + ], + "source": [ + "import cantera as ct\n", + "import sys\n", + "import os\n", + "import csv\n", + "\n", + "gas = ct.Solution('LiDryer.yaml')\n", + "gas.TP = 1001.0, ct.one_atm\n", + "gas.set_equivalence_ratio(1,\"H2\",\"O2:1, N2:3.76\")\n", + "r = ct.IdealGasConstPressureReactor(gas)\n", + "\n", + "sim = ct.ReactorNet([r])\n", + "sim.verbose = True\n", + "\n", + "# limit advance when temperature difference is exceeded\n", + "delta_T_max = 20.\n", + "r.set_advance_limit('temperature', delta_T_max)\n", + "\n", + "dt_max = 1.e-5\n", + "t_end = 100 * dt_max\n", + "states = ct.SolutionArray(gas, extra=['t'])\n", + "time = []\n", + "\n", + "print('{:10s} {:10s} {:10s} {:14s}'.format(\n", + " 't [s]', 'T [K]', 'P [Pa]', 'u [J/kg]'))\n", + "while sim.time < t_end:\n", + " sim.advance(sim.time + dt_max)\n", + " states.append(r.thermo.state, t=sim.time*1e3)\n", + " time.append(sim.time*1e3)\n", + " print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(\n", + " sim.time, r.T, r.thermo.P, r.thermo.u))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "fd5843ca", + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support. ' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
');\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('