From 8bc27e4bfe0ad54538b02c6c6855f2c96cb7061d Mon Sep 17 00:00:00 2001 From: CallieHsu Date: Tue, 19 Sep 2023 10:41:01 -0400 Subject: [PATCH] modify the way to get eps_x value --- scripts/runcode.jl | 8 +++--- src/Chamber.jl | 2 -- src/plot_figs.jl | 61 ++++++++++++++++++++++++++++++++------------- src/runcode-func.jl | 7 +----- src/utils.jl | 1 - src/write_csv.jl | 14 +++++------ 6 files changed, 54 insertions(+), 39 deletions(-) diff --git a/scripts/runcode.jl b/scripts/runcode.jl index d4fdeb9..6821e1f 100644 --- a/scripts/runcode.jl +++ b/scripts/runcode.jl @@ -224,7 +224,6 @@ julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialCo # initialize vector to store quantities param_saved_var.storeTime = Vector{Float64}([0]) param_saved_var.storeTemp = Vector{Float64}([T_0]) - param_saved_var.storeEps_x = Vector{Float64}([eps_x0]) @info("IC_Finder parameters: $(param_IC_Finder)") @@ -266,14 +265,13 @@ julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialCo print_timer_log(io, to) close(io) df = DataFrame(sol) - write_csv(df, erupt_saved, path) + write_csv(df, erupt_saved, path, composition) if plotfig title = latexstringtitle(composition, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth) plot_combined_fig(df, path, title) - plot_dual_axis(df, param_saved_var.storeTime, param_saved_var.storeEps_x, path) + # plot_figs(df, path) + plot_dual_axis(df, path) end - # write_ϵx_csv(param_saved_var.storeTime, param_saved_var.storeEps_x, path) - # plot_ϵx(param_saved_var.storeTime, param_saved_var.storeEps_x, path) return df end return df diff --git a/src/Chamber.jl b/src/Chamber.jl index b0be27e..61c02d5 100644 --- a/src/Chamber.jl +++ b/src/Chamber.jl @@ -140,8 +140,6 @@ export @memoize, latexstringtitle, plot_combined_fig, plot_dual_axis, - plot_ϵx, - write_ϵx_csv, check_for_duplicates, chamber diff --git a/src/plot_figs.jl b/src/plot_figs.jl index d4ded76..212d771 100644 --- a/src/plot_figs.jl +++ b/src/plot_figs.jl @@ -5,6 +5,7 @@ function plot_sol(df, x_axis_col, y_axis_col) xaxis=x_axis_col, yaxis=y_axis_col, linewidth=2, + legend=false, marker=(:x, 3, Plots.stroke(2)), ) end @@ -21,6 +22,7 @@ function plot_sol_year_unit(df, xlabel, ylabel) "time" => "Time (yr)", "P+dP" => "Pressure (MPa)", "T" => "Temperature (K)", + "eps_g" => L"𝜺_g", "V" => "Volume (km³)", "total_mass" => "Total mass (×10¹² kg)", ) @@ -47,27 +49,50 @@ end function plot_combined_fig(df::DataFrame, path::String, title::LaTeXString)::Nothing ENV["GKSwstype"] = "100" # magic environmental variable for Plots + # plot_names = ["P+dP", "T", "eps_g", "V", "X_CO2", "total_mass"] + # plts = [plot_sol_year_unit(df, "time", l) for l in plot_names] + # combined_plot = plot(plts..., layout=(2, 3), size=(900, 500), dpi=200, margin = 3Plots.mm, plot_title=title, plot_titlefontsize=13) + # # add eps_x in eps_g plot + # color2 = RGB(255/255, 100/255, 0/255) + # twin_subplot3 = twinx(combined_plot[3]) + # plot!(twin_subplot3, df[:, 1]/31536000, df[!, 12];linewidth=2, legend=false, linecolor=color2,fg_color_guide=color2,fg_color_text=color2, yaxis=L"𝜺_X") + # savefig(combined_plot, joinpath(path, "combined_plots.png")) + + ### + combined_plot = plot(layout=(2, 3), size=(900, 500), dpi=400, margin=3Plots.mm, plot_title=title, plot_titlefontsize=13, legend=false) plot_names = ["P+dP", "T", "eps_g", "V", "X_CO2", "total_mass"] - plts = [plot_sol_year_unit(df, "time", l) for l in plot_names] - combined_plot = plot(plts..., layout=(2, 3), size=(900, 500), dpi=200, margin = 3Plots.mm, plot_title=title, plot_titlefontsize=13) + # Units of outputs: Time -> sec, Presure -> Pa, Volume -> m³ + divisors = Dict( + "time" => 31536000, + "P+dP" => 1e6, + "V" => 1e9, + "total_mass" => 1e12, + ) + labels = Dict( + "time" => "Time (yr)", + "P+dP" => "Pressure (MPa)", + "T" => "Temperature (K)", + "eps_g" => L"𝜺_g", + "V" => "Volume (km³)", + "total_mass" => "Total mass (×10¹² kg)", + ) + convert_unit(label) = df[:, label] / get(divisors, label, 1.0) + getaxis(label) = get(labels, label, label) + x_data = convert_unit("time") + for (i, label) in enumerate(plot_names) + subplot = combined_plot[i] + y_data = convert_unit(label) + plot!(subplot, x_data, y_data; xaxis=getaxis("time"), yaxis=getaxis(label), linewidth=2, legend=false) + end + # add eps_x in eps_g plot + color2 = RGB(255/255, 100/255, 0/255) + twin_subplot3 = twinx(combined_plot[3]) + plot!(twin_subplot3, x_data, df[!, "eps_x"]; linewidth=2, yaxis=L"𝜺_X", legend=false, linecolor=color2, fg_color_guide=color2, fg_color_text=color2) savefig(combined_plot, joinpath(path, "combined_plots.png")) return nothing end -function plot_ϵx(storeTime::Vector{Float64}, storeEps_x::Vector{Float64}, path::String)::Nothing - ENV["GKSwstype"] = "100" # magic environmental variable for Plots - savefig(plot( - storeTime, storeEps_x; - xaxis="time", - yaxis="eps_x", - label="Chamber.jl", - linewidth=2, - marker=(:x, 3, Plots.stroke(2))), - "$path/eps_x.png") - return nothing -end - -function plot_dual_axis(df::DataFrame, storeTime::Vector{Float64}, storeEps_x::Vector{Float64}, path::String)::Nothing +function plot_dual_axis(df::DataFrame, path::String)::Nothing color1 = RGB(20/255, 125/255, 204/255) color2 = RGB(255/255, 100/255, 0/255) ENV["GKSwstype"] = "100" # magic environmental variable for Plots @@ -115,8 +140,8 @@ function plot_dual_axis(df::DataFrame, storeTime::Vector{Float64}, storeEps_x::V ) p2 = plot!( twinx(), - storeTime/31536000, - storeEps_x; + df[!, 1]/31536000, + df[!, 12]; common_opt2..., yaxis=L"𝜺_X", ) diff --git a/src/runcode-func.jl b/src/runcode-func.jl index a706440..cd41dc8 100644 --- a/src/runcode-func.jl +++ b/src/runcode-func.jl @@ -31,7 +31,6 @@ function odeChamber( composition = param.composition storeTime = param_saved_var.storeTime storeTemp = param_saved_var.storeTemp - storeEps_x = param_saved_var.storeEps_x phase = param_saved_var.phase c_x, c_m = param.c_x, param.c_m L_e, L_m = param.L_e, param.L_m @@ -67,18 +66,15 @@ function odeChamber( =# if storeTime[end] == t storeTemp[end] = T - storeEps_x[end] = eps_x elseif t != 0 storeTime = [storeTime; t] storeTemp = [storeTemp; T] - storeEps_x = [storeEps_x; eps_x] end cross = findfirst(!=(0), diff(sign.(diff(storeTime)))) if cross !== nothing cross_time = storeTime[end] storeTemp = [storeTemp[storeTime .< cross_time]; storeTemp[end]] - storeEps_x = [storeEps_x[storeTime .< cross_time]; storeEps_x[end]] storeTime = [storeTime[storeTime .< cross_time]; cross_time] end @@ -204,7 +200,6 @@ function odeChamber( du[10] = Mdot_c_in - Mdot_c_out param_saved_var.storeTime = storeTime param_saved_var.storeTemp = storeTemp - param_saved_var.storeEps_x = storeEps_x return du end @@ -312,7 +307,6 @@ function affect!( ) composition = param.composition param_saved_var.storeTemp = param_saved_var.storeTemp[param_saved_var.storeTime .< int.t] - param_saved_var.storeEps_x = param_saved_var.storeEps_x[param_saved_var.storeTime .< int.t] param_saved_var.storeTime = param_saved_var.storeTime[param_saved_var.storeTime .< int.t] if param.dP_lit_dt_0 == 0 @@ -342,6 +336,7 @@ function affect!( record_erupt_end(int.t, erupt_saved, param) @info("*event idx: $idx\n time: $(int.t), Finished an eruption...") elseif idx == 6 || idx == 8 + println("idx = 8..") phase_here = param_saved_var.phase @info( "*event idx: $idx\n time: $(int.t), starting ic finder for conversion of phase, phase_here: $phase_here", diff --git a/src/utils.jl b/src/utils.jl index d82d1ad..f0a52bb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -137,7 +137,6 @@ end storeSumk_2::Vector{T} = [] storeSumk_old::Vector{T} = [] storeSumk_2_old::Vector{T} = [] - storeEps_x::Vector{T} = [] end @with_kw mutable struct EruptSaved{T} diff --git a/src/write_csv.jl b/src/write_csv.jl index 87c56a7..ae3197b 100644 --- a/src/write_csv.jl +++ b/src/write_csv.jl @@ -1,4 +1,4 @@ -function write_csv(df::DataFrame, erupt_saved::EruptSaved{Float64}, path::String)::Nothing +function write_csv(df::DataFrame, erupt_saved::EruptSaved{Float64}, path::String, composition::Union{Silicic,Mafic})::Nothing number_of_data = size(df, 1) println("number_of_data: $number_of_data") names = [ @@ -15,6 +15,12 @@ function write_csv(df::DataFrame, erupt_saved::EruptSaved{Float64}, path::String ] rename!(df, ["timestamp" => "time"]) rename!(df, ["value$i" => names[i] for i in 1:10]) + # add eps_x column + m_h2o = df[:, "total_mass_H2O"] ./ df[:, "total_mass"] + m_co2 = df[:, "total_mass_CO2"] ./ df[:, "total_mass"] + crystal_fraction_eps_x′(T, P, mH2O, mCO2) = crystal_fraction_eps_x(composition, T, P, mH2O, mCO2) + df[!, "eps_x"] = crystal_fraction_eps_x′.(df[:, "T"], df[:, "P+dP"], m_h2o, m_co2) + n = length(erupt_saved.time) df_erupt = DataFrame( time_of_eruption = erupt_saved.time, @@ -26,9 +32,3 @@ function write_csv(df::DataFrame, erupt_saved::EruptSaved{Float64}, path::String CSV.write("$path/eruptions.csv", df_erupt) return nothing end - -function write_ϵx_csv(storeTime::Vector{Float64}, storeEps_x::Vector{Float64}, path::String)::Nothing - df = DataFrame(Time = storeTime, eps_x = storeEps_x) - CSV.write("$path/eps_x.csv", df) - return nothing -end