Correction on var order in ThermoFunction and improve serialization method in jfb script

Author: jfbarthelemy

data/aq17.jls data/aq17-thermofun.jlsdata/aq17.jls renamed to data/aq17-thermofun.jls

@@ -48,7 +48,8 @@ try CemSpecies(Species("Ca(OH)")) catch; "ERROR: Ca(OH) cannot be decomposed in
 CemSpecies(Species("CaCO3"; name="Calcite", aggregate_state=AS_CRYSTAL, class=SC_COMPONENT)) # ok here
 
 # Thermofun input
-function extract_database(json_file, jls_file)
+function extract_database(json_file)
+    jls_file = splitext(json_file)[1] * ".jls"
     df_substances, df_reactions = nothing, nothing
     try
         df_substances, df_reactions = deserialize(jls_file)
@@ -62,9 +63,9 @@ function extract_database(json_file, jls_file)
     dict_reactions = Dict(zip(df_reactions.symbol, df_reactions.reaction))
     return df_substances, df_reactions, dict_species, dict_reactions
 end
-df_substances, df_reactions, dict_species, dict_reactions = extract_database("data/cemdata18-merged.json", "data/cemdata18.jls")
-df_substances_psi, df_reactions_psi, dict_species_psi, dict_reactions_psi = extract_database("data/psinagra-12-07-thermofun.json", "data/psinagra.jls")
-df_substances_aq17, df_reactions_aq17, dict_species_aq17, dict_reactions_aq17 = extract_database("data/aq17-thermofun.json", "data/aq17.jls")
+df_substances, df_reactions, dict_species, dict_reactions = extract_database("data/cemdata18-merged.json")
+df_substances_psi, df_reactions_psi, dict_species_psi, dict_reactions_psi = extract_database("data/psinagra-12-07-thermofun.json")
+df_substances_aq17, df_reactions_aq17, dict_species_aq17, dict_reactions_aq17 = extract_database("data/aq17-thermofun.json")
 
 # Extraction of primaries from .dat
 df_primaries = extract_primary_species("data/CEMDATA18-31-03-2022-phaseVol.dat")
@@ -307,14 +308,18 @@ CSM = CanonicalStoichMatrix([H₂O, H⁺, OH⁻, CO₂, HCO₃⁻, CO₃²⁻]);
 SM = StoichMatrix([H₂O, H⁺, OH⁻, CO₂, HCO₃⁻, CO₃²⁻]); pprint(SM)
 
 # Example from Miron et al. 2023, https://doi.org/10.21105/joss.04624
-plot(xlabel="Temperature [°C]", ylabel="log₁₀K⁰", xticks=0:50:250, yticks=-20:2:10)
+p1 = plot(xlabel="Temperature [°C]", ylabel="Cp⁰ [J K⁻¹]", xticks=0:50:250, yticks=-2000:200:2000)
+for sp ∈ getindex.(Ref(dict_species_aq17), split("Na+ Ca+2 SiO2@ CO3-2 OH-"))
+    plot!(p1, θ->sp.Cp⁰(273.15+θ), 0:0.1:250, label=unicode(sp))
+end
+p2 = plot(xlabel="Temperature [°C]", ylabel="log₁₀K⁰", xticks=0:50:250, yticks=-20:2:10)
 for lsp ∈ [
             "Calcite Ca+2 CO3-2",
             "H2O@ H+ OH-",
             "NaCl@ Na+ Cl-",
             "Al+3 H2O@ AlOH+2 H+",
           ]
     rr = Reaction(getindex.(Ref(dict_species_aq17), split(lsp)))
-    plot!(θ->rr.logK⁰(273.15+θ), 0:0.1:250, label=rr.equation)
+    plot!(p2, θ->rr.logK⁰(273.15+θ), 0:0.1:250, label=rr.equation)
 end
-plot!()
+plot(p1, p2, layout = (1, 2))

data/cemdata18.jls data/cemdata18-merged.jlsdata/cemdata18.jls renamed to data/cemdata18-merged.jls

@@ -101,7 +101,7 @@ const thermo_function_library = Dict(
 )
 
 """
-    ThermoFunction(expr::Union{Symbol,Expr}, params=Pair[], vars=[:T, :P, :t]; ref=[])
+    ThermoFunction(expr::Union{Symbol,Expr}, params=Pair[]; vars=[:T, :P, :t, :x, :y, :z], ref=[])
 
 Construct a ThermoFunction from an expression, parameters, variables, and reference conditions.
 
@@ -127,11 +127,7 @@ julia> params = [:α => 210.0u"J/K/mol", :β => 0.0u"J/mol/K^2", :γ => -3.07e6u
  :β => 0.0 m² kg s⁻² K⁻² mol⁻¹
  :γ => -3.07e6 m² kg s⁻² K⁻¹ mol⁻¹
 
-julia> vars = [:T]
-1-element Vector{Symbol}:
- :T
-
-julia> tf = ThermoFunction(expr, params, vars; ref=[:T => 298.15u"K"])
+julia> tf = ThermoFunction(expr, params; vars=[:T], ref=[:T => 298.15u"K"])
 210.0 - 3.07e6log(T) ♢ unit=[m² kg s⁻² K⁻¹ mol⁻¹] ♢ ref=[T=298.15 K]
 
 julia> tf() # default evaluation at reference variable here Tref=298.15K
@@ -143,7 +139,7 @@ julia> tf(300.0u"K")
 julia> tf(300.0)
 -1.7510402197194535e7
 
-julia> tf_nounits = ThermoFunction(expr, [:α => 210.0, :β => 0.0, :γ => -3.07e6], vars; ref=[:T => 298.15])
+julia> tf_nounits = ThermoFunction(expr, [:α => 210.0, :β => 0.0, :γ => -3.07e6]; vars=[:T], ref=[:T => 298.15])
 210.0 - 3.07e6log(T) ♢ unit=[] ♢ ref=[T=298.15]
 
 julia> tf_nounits(300.0u"K")
@@ -158,8 +154,8 @@ julia> tf_nounits(300.)
 """
 function ThermoFunction(
     expr::Union{Symbol,Expr},
-    params=Pair[],
-    vars=[:T, :P, :t, :x, :y, :z];
+    params=Pair[];
+    vars=[:T, :P, :t, :x, :y, :z],
     ref=[],
 )
     expr = get(thermo_function_library, expr, expr)
@@ -172,6 +168,8 @@ function ThermoFunction(
                     OrderedDict(:T => 298.15, :P => 100000., :t => 0.)
     givenref = OrderedDict(ref)
     vecvars = filter(x -> Symbol(x) ∈ vars, varofexpr)
+    pos = Dict(v => i for (i, v) in enumerate(vars))  # ordre de référence
+    sort!(vecvars, by = x -> pos[Symbol(x)])
     dictvars = OrderedDict{Symbol,Num}(zip(Symbol.(vecvars), vecvars))
     dictref = OrderedDict(v=>get(givenref, v, get(default_ref, v, 0)) for v in union(keys(dictvars), keys(givenref)))
     dictvarref = OrderedDict(dictvars[v]=>get(givenref, v, get(default_ref, v, 0)) for v in keys(dictvars))