From d49c39f77c809e46c1f0d819768380bb1a4523f0 Mon Sep 17 00:00:00 2001 From: RSuryaNarayan <52670594+RSuryaNarayan@users.noreply.github.com> Date: Sun, 30 May 2021 02:52:25 +0530 Subject: [PATCH] Add files via upload --- CEMA.ipynb | 846 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 846 insertions(+) create mode 100644 CEMA.ipynb diff --git a/CEMA.ipynb b/CEMA.ipynb new file mode 100644 index 0000000..3ca6a21 --- /dev/null +++ b/CEMA.ipynb @@ -0,0 +1,846 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "21cc6738", + "metadata": {}, + "source": [ + "# A 1-D example of Chemical Explosive Mode Analysis (CEMA) for $CH_4$ ignition" + ] + }, + { + "cell_type": "markdown", + "id": "c984d88d", + "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" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "9989f948", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "PyObject " + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using Arrhenius\n", + "using ForwardDiff\n", + "using DifferentialEquations\n", + "using LinearAlgebra\n", + "using Plots\n", + "using PyCall\n", + "ct = pyimport(\"cantera\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "10536c49", + "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]))" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gas_ct = ct.Solution(\"gri30.xml\")\n", + "gas = CreateSolution(\"../mechanism/gri30.yaml\")" + ] + }, + { + "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", + "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": 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" + ] + }, + "execution_count": 50, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "plt2 = plot(sol.t, sol[end,:],lw=2, label =\"Temperature\");\n", + "plot(plt2)" + ] + }, + { + "cell_type": "code", + "execution_count": 72, + "id": "067127ea", + "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" + ] + }, + { + "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" + ] + }, + "execution_count": 72, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#compute the jacobian \n", + "function dudt(u)\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\n", + "@time J_w = ForwardDiff.jacobian(dudt,sol[end])" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "id": "33e12125", + "metadata": {}, + "outputs": [ + { + "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" + ] + }, + "execution_count": 73, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "eig_vals = eigvals(J_w)" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "id": "488bb2c0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2-element Vector{Int64}:\n", + " 53\n", + " 54" + ] + }, + "execution_count": 78, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#now extract all the positive eigen values: \n", + "findall(x->(x > 0), eig_vals) " + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "id": "19874e03", + "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", + " 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" + ] + }, + "execution_count": 79, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "eigvecs(J_w)" + ] + }, + { + "cell_type": "code", + "execution_count": 82, + "id": "927e9dc8", + "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" + ] + }, + "execution_count": 82, + "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))" + ] + }, + { + "cell_type": "code", + "execution_count": 86, + "id": "1a49e810", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "54.000000000000014" + ] + }, + "execution_count": 86, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "EP = Diagonal(Q*QL)\n", + "sum(EP)" + ] + }, + { + "cell_type": "code", + "execution_count": 87, + "id": "6db2f983", + "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" + ] + }, + "execution_count": 87, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "EI = Diagonal(Q*QL)/sum(EP)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.6.0", + "language": "julia", + "name": "julia-1.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.6.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}