Merge pull request #5 from ReactionMechanismGenerator/pyrms

Pyrms and improved flux diagrams

Author: mjohnson541

Manifest.toml

@@ -91,6 +91,12 @@ git-tree-sha1 = "70885e5e038cba1c4c17a84ad6c40756e10a4fb5"
 uuid = "19ecbf4d-ef7c-5e4b-b54a-0a0ff23c5aed"
 version = "0.5.0"
 
+[[Colors]]
+deps = ["ColorTypes","FixedPointNumbers","InteractiveUtils","Printf","Reexport"]
+name = "Colors"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.9.4"
+
 [[ColorTypes]]
 deps = ["FixedPointNumbers", "Random", "Test"]
 git-tree-sha1 = "f73b0e10f2a5756de7019818a41654686da06b09"

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

REQUIRE

@@ -13,3 +13,4 @@ SpecialFunctions
 Images
 PyPlot
 Dierckx
+Colors

src/Domain.jl

@@ -32,9 +32,8 @@ export AbstractVariableKDomain
     jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
     t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
 end
-function ConstantTPDomain(;phase::E2,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{String,E},constantspecies::Array{String,1}=Array{String,1}(),
-    sparse::Bool=false,sensitivity::Bool=false) where {E<:Real,E2<:AbstractPhase,Q<:AbstractInterface,W<:Real}
-
+function ConstantTPDomain(;phase::E2,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,X2},constantspecies::Array{X3,1}=Array{String,1}(),
+    sparse::Bool=false,sensitivity::Bool=false) where {E<:Real,E2<:AbstractPhase,Q<:AbstractInterface,W<:Real,X,X2,X3}
     #set conditions and initialconditions
     T = 0.0
     P = 0.0
@@ -106,8 +105,8 @@ export ConstantTPDomain
     jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
     t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
 end
-function ConstantVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{String,E},constantspecies::Array{String,1}=Array{String,1}(),
-    sparse::Bool=false,sensitivity::Bool=false) where {E<:Real,Z<:IdealGas,Q<:AbstractInterface}
+function ConstantVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,E},constantspecies::Array{X2,1}=Array{String,1}(),
+    sparse::Bool=false,sensitivity::Bool=false) where {E,X,X2,Z<:IdealGas,Q<:AbstractInterface}
 
     #set conditions and initialconditions
     T = 0.0
@@ -174,8 +173,8 @@ export ConstantVDomain
     jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
     t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
 end
-function ParametrizedTPDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{String,Any},constantspecies::Array{String,1}=Array{String,1}(),
-    sparse::Bool=false,sensitivity::Bool=false) where {Z<:IdealGas,Q<:AbstractInterface}
+function ParametrizedTPDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,Any},constantspecies::Array{X2,1}=Array{String,1}(),
+    sparse::Bool=false,sensitivity::Bool=false) where {X,X2,Z<:IdealGas,Q<:AbstractInterface}
 
     #set conditions and initialconditions
     T = 0.0
@@ -257,8 +256,8 @@ export ParametrizedTPDomain
     jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
     t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
 end
-function ParametrizedVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{String,Any},constantspecies::Array{String,1}=Array{String,1}(),
-    sparse::Bool=false,sensitivity::Bool=false) where {E<:Real,Z<:IdealGas,Q<:AbstractInterface}
+function ParametrizedVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,Any},constantspecies::Array{X2,1}=Array{String,1}(),
+    sparse::Bool=false,sensitivity::Bool=false) where {X,X2,E<:Real,Z<:IdealGas,Q<:AbstractInterface}
 
     #set conditions and initialconditions
     T = 0.0
@@ -340,8 +339,8 @@ export ParametrizedVDomain
     jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
     t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
 end
-function ConstantTVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{String,E},constantspecies::Array{String,1}=Array{String,1}(),
-    sparse=false,sensitivity=false) where {E<:Real, Z<:AbstractPhase,Q<:AbstractInterface,W<:Real}
+function ConstantTVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,E},constantspecies::Array{X2,1}=Array{String,1}(),
+    sparse=false,sensitivity=false) where {E,X,X2, Z<:AbstractPhase,Q<:AbstractInterface,W<:Real}
     #set conditions and initialconditions
     T = 0.0
     V = 0.0
@@ -415,8 +414,8 @@ export ConstantTVDomain
     jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
     t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
 end
-function ParametrizedTConstantVDomain(;phase::IdealDiluteSolution,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{String,Any},constantspecies::Array{String,1}=Array{String,1}(),
-    sparse::Bool=false,sensitivity::Bool=false) where {Q<:AbstractInterface}
+function ParametrizedTConstantVDomain(;phase::IdealDiluteSolution,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,X3},constantspecies::Array{X2,1}=Array{String,1}(),
+    sparse::Bool=false,sensitivity::Bool=false) where {X,X2,X3,Q<:AbstractInterface}
     #set conditions and initialconditions
     T = 0.0
     P = 0.0

src/ReactionMechanismSimulator.jl

@@ -19,15 +19,15 @@ module ReactionMechanismSimulator
         copy!(Chem,pyimport_conda("rdkit.Chem","rdkit","rmg"))
         copy!(molecule,pyimport_conda("rmgpy.molecule","rmgpy","rmg"))
         copy!(pydot,pyimport_conda("pydot","pydot","rmg"))
-        copy!(chemkin,pyimport_conda("rmgpy.chemkin","rmgpy","rmg"))
-        copy!(species,pyimport_conda("rmgpy.species","rmgpy","rmg"))
-        copy!(reaction,pyimport_conda("rmgpy.reaction","rmgpy","rmg"))
-        copy!(nasa,pyimport_conda("rmgpy.thermo.nasa","rmgpy","rmg"))
-        copy!(wilhoit,pyimport_conda("rmgpy.thermo.wilhoit","rmgpy","rmg"))
-        copy!(arrhenius,pyimport_conda("rmgpy.kinetics.arrhenius","rmgpy","rmg"))
-        copy!(falloff,pyimport_conda("rmgpy.kinetics.falloff","rmgpy","rmg"))
-        copy!(chebyshev,pyimport_conda("rmgpy.kinetics.chebyshev","rmgpy","rmg"))
-        copy!(solvation,pyimport_conda("rmgpy.data.solvation","rmgpy","rmg"))
+        copy!(chemkin,pyimport_conda("rmgpy.chemkin","rmg","rmg"))
+        copy!(species,pyimport_conda("rmgpy.species","rmg","rmg"))
+        copy!(reaction,pyimport_conda("rmgpy.reaction","rmg","rmg"))
+        copy!(nasa,pyimport_conda("rmgpy.thermo.nasa","rmg","rmg"))
+        copy!(wilhoit,pyimport_conda("rmgpy.thermo.wilhoit","rmg","rmg"))
+        copy!(arrhenius,pyimport_conda("rmgpy.kinetics.arrhenius","rmg","rmg"))
+        copy!(falloff,pyimport_conda("rmgpy.kinetics.falloff","rmg","rmg"))
+        copy!(chebyshev,pyimport_conda("rmgpy.kinetics.chebyshev","rmg","rmg"))
+        copy!(solvation,pyimport_conda("rmgpy.data.solvation","rmg","rmg"))
         copy!(yaml,pyimport_conda("yaml","yaml"))
         copy!(os,pyimport_conda("os","os"))
     end

src/fluxdiagrams.jl

@@ -1,6 +1,7 @@
 using PyCall
 using SparseArrays
 using Images
+using Colors
 import Base: length
 
 struct FluxDiagram{T<:Real}
@@ -63,16 +64,18 @@ export drawspecies
 
 """
 generates and returns the image of a single flux diagram at the given time point
+all PyPlot colorscheme names are valid inputs for colorscheme
 """
 function getfluxdiagram(bsol,t;centralspecieslist=Array{String,1}(),superimpose=false,
     maximumnodecount=50, maximumedgecount=50, concentrationtol=1e-6, speciesratetolerance=1e-6,
-    maximumnodepenwidth=10.0,maximumedgepenwidth=10.0,radius=1,centralreactioncount=-1,outputdirectory="fluxdiagrams")
+    maximumnodepenwidth=10.0,maximumedgepenwidth=10.0,radius=1,centralreactioncount=-1,outputdirectory="fluxdiagrams",
+    colorscheme="viridis")
 
     fd = makefluxdiagrams(bsol,[t]; centralspecieslist=centralspecieslist,superimpose=superimpose,
         maximumnodecount=maximumnodecount, maximumedgecount=maximumedgecount, concentrationtol=concentrationtol,
         speciesratetolerance=speciesratetolerance,maximumnodepenwidth=maximumnodepenwidth,
         maximumedgepenwidth=maximumedgepenwidth,radius=radius,centralreactioncount=centralreactioncount,
-        outputdirectory=outputdirectory)
+        outputdirectory=outputdirectory,colorscheme=colorscheme)
 
     return getdiagram(fd,1)
 end
@@ -81,10 +84,12 @@ export getfluxdiagram
 """
 generates a series of flux diagrams at the time points indicated
 each flux diagram will have the same nodes and edges as determined by the options
+all PyPlot colorscheme names are valid inputs for colorscheme
 """
 function makefluxdiagrams(bsol,ts;centralspecieslist=Array{String,1}(),superimpose=false,
     maximumnodecount=50, maximumedgecount=50, concentrationtol=1e-6, speciesratetolerance=1e-6,
-    maximumnodepenwidth=10.0,maximumedgepenwidth=10.0,radius=1,centralreactioncount=-1,outputdirectory="fluxdiagrams")
+    maximumnodepenwidth=10.0,maximumedgepenwidth=10.0,radius=1,centralreactioncount=-1,outputdirectory="fluxdiagrams",
+    colorscheme="viridis")
 
     specieslist = bsol.domain.phase.species
     speciesnamelist = getfield.(specieslist,:name)
@@ -262,6 +267,15 @@ function makefluxdiagrams(bsol,ts;centralspecieslist=Array{String,1}(),superimpo
         end
 
         slope = -maximumedgepenwidth / log10(speciesratetolerance)
+        minspeciesrate = Inf
+        for index in 1:length(edges)
+            reactantindex,productindex = edges[index]
+            sprate = abs(speciesrates[reactantindex,productindex,t])
+            if sprate > speciesratetolerance && minspeciesrate > sprate
+                minspeciesrate = sprate
+            end
+        end
+
         for index in 1:length(edges)
             reactantindex,productindex = edges[index]
             if reactantindex in nodes && productindex in nodes
@@ -298,6 +312,7 @@ function makefluxdiagrams(bsol,ts;centralspecieslist=Array{String,1}(),superimpo
                 end
 
                 edge[:set_penwidth](penwidth)
+                edge[:set_color](getcolor(speciesrates[reactantindex,productindex,t],maxspeciesrate,minspeciesrate,colorscheme))
 
             end
         end
@@ -312,6 +327,7 @@ function makefluxdiagrams(bsol,ts;centralspecieslist=Array{String,1}(),superimpo
         graph[:set_label](label)
         graph[:write_dot](joinpath(outputdirectory,"flux_diagram_$t.dot"))
         graph[:write_png](joinpath(outputdirectory,"flux_diagram_$t.png"))
+        graph[:write_svg](joinpath(outputdirectory,"flux_diagram_$t.svg"))
     end
     return FluxDiagram(ts,outputdirectory)
 end
@@ -369,3 +385,14 @@ function addadjacentnodes!(targetnodeindex,nodes,edges,phase,maxreactionrates,ma
         end
     end
 end
+
+function getcolor(speciesrate,maxspeciesrate,minspeciesrate,colorscheme="viridis")
+    """
+    gives the color corresponding to the scaled log species rate
+    for a given PyPlot color scheme
+    """
+    scale = log(maxspeciesrate)-log(minspeciesrate)
+    value = (log(abs(speciesrate))-log(minspeciesrate))/scale
+    out = PyPlot.get_cmap(colorscheme)(value)[1:3]
+    return "#"*hex(RGB(out...))
+end