diff --git a/julia/MAGEMin_wrappers.jl b/julia/MAGEMin_wrappers.jl index c06e68fe..4da8a600 100755 --- a/julia/MAGEMin_wrappers.jl +++ b/julia/MAGEMin_wrappers.jl @@ -1,15 +1,12 @@ # The full functionality of MAGEMin is wrapped in ../gen/magemin_library.jl # Yet, the routines here make it more convenient to use this from julia import Base.show -using Base.Threads: @threads +using Base.Threads: @threads using ProgressMeter export init_MAGEMin, finalize_MAGEMin, point_wise_minimization, convertBulk4MAGEMin, use_predefined_bulk_rock, define_bulk_rock, create_output, - print_info, create_gmin_struct, - single_point_minimization, - multi_point_minimization, MAGEMin_Data, - Initialize_MAGEMin, Finalize_MAGEMin - + print_info, create_gmin_struct, single_point_minimization, multi_point_minimization, MAGEMin_Data, Initialize_MAGEMin, Finalize_MAGEMin + """ Holds the MAGEMin databases & required structures for every thread @@ -25,7 +22,7 @@ end """ Dat = Initialize_MAGEMin(db = "ig"; verbose::Union{Bool, Int64} = true) -Initializes MAGEMin on one or more threads, for the database `db`. You can surpress all output with `verbose=false`. `verbose=true` will give a brief summary of the result, whereas `verbose=1` will give more details about the computations. +Initializes MAGEMin on one or more threads, for the database `db`. You can suppress all output with `verbose=false`. `verbose=true` will give a brief summary of the result, whereas `verbose=1` will give more details about the computations. """ function Initialize_MAGEMin(db = "ig"; verbose::Union{Int64,Bool} = 0) gv, z_b, DB, splx_data = init_MAGEMin(db); @@ -82,10 +79,10 @@ function init_MAGEMin(db="ig") z_b = LibMAGEMin.bulk_infos() gv = LibMAGEMin.global_variables() - splx_data = LibMAGEMin.simplex_data(); + splx_data = LibMAGEMin.simplex_data(); DB = LibMAGEMin.Database() gv = LibMAGEMin.global_variable_alloc( pointer_from_objref(z_b)) - + if db == "ig" gv.EM_database = 2 elseif db == "mp" @@ -114,7 +111,7 @@ end # wrapper for single point minimization function single_point_minimization( P::Float64, T::Float64, - MAGEMin_db::MAGEMin_Data; + MAGEMin_db::MAGEMin_Data; test::Int64 = 0, # if using a build-in test case X::Union{Nothing, Vector{_T}, Vector{Vector{_T}}} = nothing, Xoxides = Vector{String}, @@ -170,7 +167,7 @@ julia> P = fill(10.0,n) julia> T = fill(1100.0,n) julia> Xoxides = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "Cr2O3"; "H2O"]; julia> X = [48.43; 15.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; -julia> sys_in = "wt" +julia> sys_in = "wt" julia> out = multi_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in) julia> Finalize_MAGEMin(data) ``` @@ -185,7 +182,7 @@ julia> Xoxides = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; julia> X1 = [48.43; 15.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; julia> X2 = [49.43; 14.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; julia> X = [X1,X2] -julia> sys_in = "wt" +julia> sys_in = "wt" julia> out = multi_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in) julia> Finalize_MAGEMin(data) ``` @@ -201,12 +198,12 @@ which will automatically use all threads on your machine. Alternatively, use `ju If you are interested to see what you can do on your machine type: ``` julia> versioninfo() -``` +``` """ function multi_point_minimization( P::Vector{Float64}, T::Vector{Float64}, - MAGEMin_db::MAGEMin_Data; + MAGEMin_db::MAGEMin_Data; test = 0, # if using a build-in test case X::Union{Nothing, Vector{_T}, Vector{Vector{_T}}}=nothing, Xoxides = Vector{String}, @@ -214,10 +211,10 @@ function multi_point_minimization( P::Vector{Float64}, progressbar = true # show a progress bar or not? ) where _T <: Float64 - # Set the compositional info + # Set the compositional info CompositionType::Int64 = 0; if isnothing(X) - # Use one of the build-in tests + # Use one of the build-in tests # Create thread-local data for i in 1:Threads.nthreads() MAGEMin_db.gv[i] = use_predefined_bulk_rock(MAGEMin_db.gv[i], test, MAGEMin_db.db) @@ -278,7 +275,7 @@ function multi_point_minimization( P::Vector{Float64}, # compute a new point using a ccall out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data) - Out_PT[i] = deepcopy(out) + Out_PT[i] = deepcopy(out) if progressbar next!(progr) @@ -349,32 +346,43 @@ end """ convertBulk4MAGEMin( bulk_in, bulk_in_ox, sys_in) - + receives bulk-rock composition in [mol,wt] fraction and associated oxide list and sends back bulk-rock composition converted for MAGEMin use - + """ function convertBulk4MAGEMin(bulk_in::Vector{Float64},bulk_in_ox::Vector{String},sys_in::String,db::String); - ref_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "O"; "Cr2O3"; "MnO"; "H2O"]; - ref_MolarMass = [60.08; 101.96; 56.08; 40.30; 71.85; 159.69; 94.2; 61.98; 79.88; 16.0; 151.99; 70.937; 18.015]; #Molar mass of oxides + ref_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "O"; "Cr2O3"; "MnO"; "H2O"; "S"]; + ref_MolarMass = [60.08; 101.96; 56.08; 40.30; 71.85; 159.69; 94.2; 61.98; 79.88; 16.0; 151.99; 70.937; 18.015; 32.06]; #Molar mass of oxides if db == "ig" MAGEMin_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "K2O"; "Na2O"; "TiO2"; "O"; "Cr2O3"; "H2O"]; + elseif db == "igd" + MAGEMin_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "K2O"; "Na2O"; "TiO2"; "O"; "Cr2O3"; "H2O"]; + elseif db == "alk" + MAGEMin_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "K2O"; "Na2O"; "TiO2"; "O"; "Cr2O3"; "H2O"]; + elseif db == "mb" + MAGEMin_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "K2O"; "Na2O"; "TiO2"; "O"; "H2O"]; + elseif db == "um" + MAGEMin_ox = ["SiO2"; "Al2O3"; "MgO" ;"FeO"; "O"; "H2O"; "S"]; elseif db == "mp" MAGEMin_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "K2O"; "Na2O"; "TiO2"; "O"; "MnO"; "H2O"]; else print("Database not implemented...\n") end - MAGEMin_bulk = zeros(11); - bulk = zeros(11); + MAGEMin_bulk = zeros(length(MAGEMin_ox)); + bulk = zeros(length(MAGEMin_ox)); # convert to mol, if system unit = wt if sys_in == "wt" for i=1:length(bulk_in_ox) id = findall(ref_ox .== bulk_in_ox[i]); bulk[i] = bulk_in[i]/ref_MolarMass[id[1]]; end + else + bulk .= bulk_in; end + bulk = normalize(bulk); for i=1:length(MAGEMin_ox) @@ -404,6 +412,7 @@ function convertBulk4MAGEMin(bulk_in::Vector{Float64},bulk_in_ox::Vector{String} MAGEMin_bulk .= normalize(MAGEMin_bulk)*100.0; + # print("MAGEMin_bulk $MAGEMin_bulk\n") return MAGEMin_bulk; end @@ -413,9 +422,9 @@ end """ point_wise_minimization(P::Float64,T::Float64, gv, z_b, DB, splx_data, sys_in::String="mol") - + Computes the stable assemblage at `P` [kbar], `T` [C] and for a given bulk rock composition - + # Example 1 @@ -431,24 +440,24 @@ julia> T = 800.0; julia> gv.verbose = -1; # switch off any verbose julia> out = point_wise_minimization(P,T, gv, z_b, DB, splx_data, sys_in) Pressure : 8.0 [kbar] -Temperature : 800.0 [Celcius] - Stable phase | Fraction (mol 1 atom basis) - opx 0.24229 - ol 0.58808 - cpx 0.14165 - spn 0.02798 - Stable phase | Fraction (wt fraction) - opx 0.23908 - ol 0.58673 - cpx 0.14583 - spn 0.02846 +Temperature : 800.0 [Celsius] + Stable phase | Fraction (mol 1 atom basis) + opx 0.24229 + ol 0.58808 + cpx 0.14165 + spn 0.02798 + Stable phase | Fraction (wt fraction) + opx 0.23908 + ol 0.58673 + cpx 0.14583 + spn 0.02846 Gibbs free energy : -797.749183 (26 iterations; 94.95 ms) Oxygen fugacity : 9.645393319147175e-12 ``` # Example 2 -And here a case in which you specify your own bulk rock composition. -We convert that in the correct format, using the `convertBulk4MAGEMin` function. +And here a case in which you specify your own bulk rock composition. +We convert that in the correct format, using the `convertBulk4MAGEMin` function. ```julia julia> db = "ig" # database: ig, igneous (Holland et al., 2018); mp, metapelite (White et al 2014b) julia> gv, z_b, DB, splx_data = init_MAGEMin(db); @@ -460,17 +469,17 @@ julia> P,T = 10.0, 1100.0; julia> gv.verbose = -1; # switch off any verbose julia> out = point_wise_minimization(P,T, gv, z_b, DB, splx_data, sys_in) Pressure : 10.0 [kbar] -Temperature : 1100.0 [Celcius] - Stable phase | Fraction (mol 1 atom basis) - pl4T 0.01114 - liq 0.74789 - cpx 0.21862 - opx 0.02154 - Stable phase | Fraction (wt fraction) - pl4T 0.01168 - liq 0.72576 - cpx 0.23872 - opx 0.02277 +Temperature : 1100.0 [Celsius] + Stable phase | Fraction (mol 1 atom basis) + pl4T 0.01114 + liq 0.74789 + cpx 0.21862 + opx 0.02154 + Stable phase | Fraction (wt fraction) + pl4T 0.01168 + liq 0.72576 + cpx 0.23872 + opx 0.02277 Gibbs free energy : -907.27887 (47 iterations; 187.11 ms) Oxygen fugacity : 0.02411835177808492 julia> finalize_MAGEMin(gv,DB) @@ -478,7 +487,7 @@ julia> finalize_MAGEMin(gv,DB) """ function point_wise_minimization(P::Float64,T::Float64, gv, z_b, DB, splx_data) - + input_data = LibMAGEMin.io_data(); # zero (not used actually) z_b.T = T + 273.15 # in K @@ -488,7 +497,7 @@ function point_wise_minimization(P::Float64,T::Float64, gv, z_b, DB, splx_data) # Perform the point-wise minimization after resetting variables gv = LibMAGEMin.reset_gv(gv,z_b, DB.PP_ref_db, DB.SS_ref_db) - z_b = LibMAGEMin.reset_z_b_bulk( gv, z_b ) + z_b = LibMAGEMin.reset_z_b_bulk( gv, z_b ) LibMAGEMin.reset_simplex_A(pointer_from_objref(splx_data), z_b, gv) LibMAGEMin.reset_simplex_B_em(pointer_from_objref(splx_data), gv) @@ -505,12 +514,12 @@ function point_wise_minimization(P::Float64,T::Float64, gv, z_b, DB, splx_data) # Fill structure LibMAGEMin.fill_output_struct( gv, z_b, DB.PP_ref_db,DB.SS_ref_db, DB.cp, DB.sp ); - # Print output to screen - LibMAGEMin.PrintOutput(gv, 0, 1, DB, time, z_b); + # Print output to screen + LibMAGEMin.PrintOutput(gv, 0, 1, DB, time, z_b); # Transform results to a more convenient julia struct out = create_gmin_struct(DB, gv, time); - + # LibMAGEMin.FreeDatabases(gv, DB); return out @@ -526,44 +535,44 @@ point_wise_minimization(P::Number,T::Number, gv::MAGEMin_C.LibMAGEMin.global_var Performs a point-wise optimization for a given pressure `P` and temperature `T` foir the data specified in the MAGEMin database `MAGEMin_Data` (where also compoition is specified) """ -point_wise_minimization(P::Number,T::Number, data::MAGEMin_Data) = point_wise_minimization(P,T, data.gv[1], data.z_b[1], data.DB[1], data.splx_data[1]) +point_wise_minimization(P::Number,T::Number, data::MAGEMin_Data) = point_wise_minimization(P,T, data.gv[1], data.z_b[1], data.DB[1], data.splx_data[1]) """ - structure that holds the result of the pointwise minisation + structure that holds the result of the pointwise minisation """ struct gmin_struct{T,I} G_system::T # G of system Gamma::Vector{T} # Gamma P_kbar::T # Pressure in kbar - T_C::T # Temperature in Celcius + T_C::T # Temperature in Celsius # bulk rock composition: - bulk::Vector{T} - bulk_M::Vector{T} - bulk_S::Vector{T} - bulk_F::Vector{T} - - bulk_wt::Vector{T} - bulk_M_wt::Vector{T} - bulk_S_wt::Vector{T} - bulk_F_wt::Vector{T} - + bulk::Vector{T} + bulk_M::Vector{T} + bulk_S::Vector{T} + bulk_F::Vector{T} + + bulk_wt::Vector{T} + bulk_M_wt::Vector{T} + bulk_S_wt::Vector{T} + bulk_F_wt::Vector{T} + # Fractions: # Solid, melt, fluid fractions - frac_M::T + frac_M::T frac_S::T frac_F::T - frac_M_wt::T + frac_M_wt::T frac_S_wt::T frac_F_wt::T - + # Solid, melt, fluid densities rho::T - rho_M::T + rho_M::T rho_S::T rho_F::T @@ -573,16 +582,16 @@ struct gmin_struct{T,I} # Phase fractions and type: n_PP::Int64 # number of pure phases n_SS::Int64 # number of solid solutions - + ph_frac::Vector{T} # phase fractions ph_frac_wt::Vector{T} # phase fractions ph_type::Vector{I} # type of phase (SS or PP) ph_id::Vector{I} # id of phase ph::Vector{String} # Name of phase - + SS_vec::Vector{LibMAGEMin.SS_data} PP_vec::Vector{LibMAGEMin.PP_data} - + oxides::Vector{String} # Seismic velocity info @@ -619,8 +628,8 @@ function create_gmin_struct(DB, gv, time) G_system = stb.G; Gamma = unsafe_wrap(Vector{Cdouble},stb.gamma,gv.len_ox) P_kbar = stb.P - T_C = stb.T-273.15 - + T_C = stb.T-273.15 + # Bulk rock info (total, melt, solid, fluid) bulk = unsafe_wrap(Vector{Cdouble},stb.bulk, gv.len_ox) bulk_M = unsafe_wrap(Vector{Cdouble},stb.bulk_M, gv.len_ox) @@ -632,20 +641,20 @@ function create_gmin_struct(DB, gv, time) bulk_M_wt = unsafe_wrap(Vector{Cdouble},stb.bulk_M_wt, gv.len_ox) bulk_S_wt = unsafe_wrap(Vector{Cdouble},stb.bulk_S_wt, gv.len_ox) bulk_F_wt = unsafe_wrap(Vector{Cdouble},stb.bulk_F_wt, gv.len_ox) - + # Solid, melt, fluid fractions - frac_M = stb.frac_M + frac_M = stb.frac_M frac_S = stb.frac_S frac_F = stb.frac_F # Solid, melt, fluid fractions - frac_M_wt = stb.frac_M_wt + frac_M_wt = stb.frac_M_wt frac_S_wt = stb.frac_S_wt frac_F_wt = stb.frac_F_wt # Solid, melt, fluid densities rho = stb.rho - rho_M = stb.rho_M + rho_M = stb.rho_M rho_S = stb.rho_S rho_F = stb.rho_F @@ -660,67 +669,67 @@ function create_gmin_struct(DB, gv, time) n_ph = stb.n_ph # total # of stable phases n_PP = stb.n_PP # number of pure phases n_SS = stb.n_SS # number of solid solutions - + ph_frac = unsafe_wrap(Vector{Cdouble},stb.ph_frac, n_ph) ph_frac_wt = unsafe_wrap(Vector{Cdouble},stb.ph_frac_wt, n_ph) ph_type = unsafe_wrap(Vector{Cint}, stb.ph_type, n_ph) ph_id = unsafe_wrap(Vector{Cint}, stb.ph_id , n_ph) ph = unsafe_string.(unsafe_wrap(Vector{Ptr{Int8}}, stb.ph, n_ph)) # stable phases - + # extract info about compositional variables of the solution models: SS_vec = convert.(LibMAGEMin.SS_data, unsafe_wrap(Vector{LibMAGEMin.stb_SS_phase},stb.SS,n_SS)) - + # Info about the endmembers: - PP_vec = convert.(LibMAGEMin.PP_data, unsafe_wrap(Vector{LibMAGEMin.stb_PP_phase},stb.PP,n_PP)) + PP_vec = convert.(LibMAGEMin.PP_data, unsafe_wrap(Vector{LibMAGEMin.stb_PP_phase},stb.PP,n_PP)) # Names of oxides: oxides = unsafe_string.(unsafe_wrap(Vector{Ptr{Int8}}, stb.oxides, gv.len_ox)) - + # Numerics bulk_res_norm = gv.BR_norm - iter = gv.global_ite + iter = gv.global_ite time_ms = time*1000.0 - # Store all in output struct - out = gmin_struct{Float64,Int64}( G_system, Gamma, P_kbar, T_C, + # Store all in output struct + out = gmin_struct{Float64,Int64}( G_system, Gamma, P_kbar, T_C, bulk, bulk_M, bulk_S, bulk_F, - bulk_wt, bulk_M_wt, bulk_S_wt, bulk_F_wt, - frac_M, frac_S, frac_F, - frac_M_wt, frac_S_wt, frac_F_wt, - rho, rho_M, rho_S, rho_F, - fO2, + bulk_wt, bulk_M_wt, bulk_S_wt, bulk_F_wt, + frac_M, frac_S, frac_F, + frac_M_wt, frac_S_wt, frac_F_wt, + rho, rho_M, rho_S, rho_F, + fO2, n_PP, n_SS, ph_frac, ph_frac_wt, ph_type, ph_id, ph, - SS_vec, PP_vec, - oxides, + SS_vec, PP_vec, + oxides, stb.Vp, stb.Vs, stb.Vp_S, stb.Vs_S, stb.bulkMod, stb.shearMod, stb.bulkModulus_M, stb.bulkModulus_S, stb.shearModulus_S, entropy, enthalpy, iter, bulk_res_norm, time_ms, stb.status) - + return out end -# Print brief info about pointwise calculation result -function show(io::IO, g::gmin_struct) - println(io, "Pressure : $(g.P_kbar) [kbar]") - println(io, "Temperature : $(round(g.T_C,digits=4)) [Celcius]") - - println(io, " Stable phase | Fraction (mol 1 atom basis) ") +# Print brief info about pointwise calculation result +function show(io::IO, g::gmin_struct) + println(io, "Pressure : $(g.P_kbar) [kbar]") + println(io, "Temperature : $(round(g.T_C,digits=4)) [Celsius]") + + println(io, " Stable phase | Fraction (mol 1 atom basis) ") for i=1:length(g.ph) - println(io, " $(lpad(g.ph[i],14," ")) $( round(g.ph_frac[i], digits=5)) ") + println(io, " $(lpad(g.ph[i],14," ")) $( round(g.ph_frac[i], digits=5)) ") end - println(io, " Stable phase | Fraction (wt fraction) ") + println(io, " Stable phase | Fraction (wt fraction) ") for i=1:length(g.ph) - println(io, " $(lpad(g.ph[i],14," ")) $( round(g.ph_frac_wt[i], digits=5)) ") + println(io, " $(lpad(g.ph[i],14," ")) $( round(g.ph_frac_wt[i], digits=5)) ") end - println(io, "Gibbs free energy : $(round(g.G_system,digits=6)) ($(g.iter) iterations; $(round(g.time_ms,digits=2)) ms)") + println(io, "Gibbs free energy : $(round(g.G_system,digits=6)) ($(g.iter) iterations; $(round(g.time_ms,digits=2)) ms)") if g.status>0 - println(io, "WARNING: calculation did not converge ----------------------------") + println(io, "WARNING: calculation did not converge ----------------------------") end - println(io, "Oxygen fugacity : $(g.fO2)") + println(io, "Oxygen fugacity : $(g.fO2)") + - end """ @@ -729,19 +738,19 @@ end Prints a more extensive overview of the simulation results """ function print_info(g::gmin_struct) - + println("Stable phases @ {$(round(g.P_kbar,digits=4)), $(round(g.T_C,digits=4))} kbar/°C :") for i=1:length(g.ph) - print("$(lpad(g.ph[i],6," ")) ") + print("$(lpad(g.ph[i],6," ")) ") end - print(" \n \n") - + print(" \n \n") + # == println("Compositional variables (solution phase):") for i=1:g.n_SS - print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(g.ph[i],15," ")) ") for j=1:length(g.SS_vec[i].compVariables) - print("$(lpad(round(g.SS_vec[i].compVariables[j],digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].compVariables[j],digits=5),8," ")) ") end print("\n") end @@ -754,12 +763,12 @@ function print_info(g::gmin_struct) print(" ") for j=1:length(g.SS_vec[i].emNames) - print("$(lpad(g.SS_vec[i].emNames[j],8," ")) ") + print("$(lpad(g.SS_vec[i].emNames[j],8," ")) ") end print("\n") - print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(g.ph[i],15," ")) ") for j=1:length(g.SS_vec[i].emFrac) - print("$(lpad(round(g.SS_vec[i].emFrac[j],digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].emFrac[j],digits=5),8," ")) ") end print("\n") end @@ -770,12 +779,12 @@ function print_info(g::gmin_struct) print(" ") for j=1:length(g.SS_vec[i].emNames) - print("$(lpad(g.SS_vec[i].emNames[j],8," ")) ") + print("$(lpad(g.SS_vec[i].emNames[j],8," ")) ") end print("\n") - print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(g.ph[i],15," ")) ") for j=1:length(g.SS_vec[i].emFrac_wt) - print("$(lpad(round(g.SS_vec[i].emFrac_wt[j],digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].emFrac_wt[j],digits=5),8," ")) ") end print("\n") end @@ -783,30 +792,30 @@ function print_info(g::gmin_struct) # == # == println("Oxide compositions [mol% 1 atom basis] (normalized):") - print(" ") + print(" ") for i=1:length(g.oxides) - print("$(lpad(g.oxides[i],8," ")) ") + print("$(lpad(g.oxides[i],8," ")) ") end print("\n") - - print("$(lpad("SYS",15)) ") + + print("$(lpad("SYS",15)) ") for i=1:length(g.oxides) - print("$(lpad(round(g.bulk[i],digits=5),8," ")) ") + print("$(lpad(round(g.bulk[i],digits=5),8," ")) ") end print("\n") for i=1:g.n_SS - print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(g.ph[i],15," ")) ") for j=1:length(g.oxides) - print("$(lpad(round(g.SS_vec[i].Comp[j],digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].Comp[j],digits=5),8," ")) ") end print("\n") end for i=1:g.n_PP - print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(g.ph[i],15," ")) ") for j=1:length(g.oxides) - print("$(lpad(round(g.PP_vec[i].Comp[j],digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].Comp[j],digits=5),8," ")) ") end print("\n") end @@ -818,31 +827,31 @@ function print_info(g::gmin_struct) # == # == println("Oxide compositions [wt%] (normalized):") - print(" ") + print(" ") for i=1:length(g.oxides) - print("$(lpad(g.oxides[i],8," ")) ") + print("$(lpad(g.oxides[i],8," ")) ") end print("\n") - - print("$(lpad("SYS",15)) ") + + print("$(lpad("SYS",15)) ") for i=1:length(g.oxides) - print("$(lpad(round(g.bulk_wt[i],digits=5),8," ")) ") + print("$(lpad(round(g.bulk_wt[i],digits=5),8," ")) ") end print("\n") for i=1:g.n_SS - print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(g.ph[i],15," ")) ") for j=1:length(g.oxides) - print("$(lpad(round(g.SS_vec[i].Comp_wt[j],digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].Comp_wt[j],digits=5),8," ")) ") end print("\n") end print("\n") for i=1:g.n_PP - print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(g.ph[i],15," ")) ") for j=1:length(g.oxides) - print("$(lpad(round(g.PP_vec[i].Comp_wt[j],digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].Comp_wt[j],digits=5),8," ")) ") end print("\n") end @@ -850,85 +859,85 @@ function print_info(g::gmin_struct) # == # == - + println("Stable mineral assemblage:") println(" phase mode[mol1at] mode[wt] f G V Cp rho[kg/m3] Thermal_Exp Entropy[J/K] Enthalpy[J] BulkMod[GPa] ShearMod[GPa] Vp[km/s] Vs[km/s]") for i=1:g.n_SS - print("$(lpad(g.ph[i],15," ")) ") - print("$(lpad(round(g.ph_frac[i],digits=5),13," ")) ") - print("$(lpad(round(g.ph_frac_wt[i],digits=5),8," ")) ") - print("$(lpad(round(g.SS_vec[i].f,digits=5),8," ")) ") - print("$(lpad(round(g.SS_vec[i].G,digits=5),8," ")) ") - print("$(lpad(round(g.SS_vec[i].V,digits=5),8," ")) ") - print("$(lpad(round(g.SS_vec[i].cp,digits=5),8," ")) ") - print("$(lpad(round(g.SS_vec[i].rho,digits=5),11," ")) ") - print("$(lpad(round(g.SS_vec[i].alpha,digits=5),8," ")) ") - print("$(lpad(round(g.SS_vec[i].entropy,digits=5),8," ")) ") - print("$(lpad(round(g.SS_vec[i].enthalpy,digits=5),8," ")) ") - print("$(lpad(round(g.SS_vec[i].bulkMod,digits=5),12," ")) ") - print("$(lpad(round(g.SS_vec[i].shearMod,digits=5),13," ")) ") - print("$(lpad(round(g.SS_vec[i].Vp,digits=5),10," ")) ") - print("$(lpad(round(g.SS_vec[i].Vs,digits=5),10," ")) ") - + print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(round(g.ph_frac[i],digits=5),13," ")) ") + print("$(lpad(round(g.ph_frac_wt[i],digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].f,digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].G,digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].V,digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].cp,digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].rho,digits=5),11," ")) ") + print("$(lpad(round(g.SS_vec[i].alpha,digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].entropy,digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].enthalpy,digits=5),8," ")) ") + print("$(lpad(round(g.SS_vec[i].bulkMod,digits=5),12," ")) ") + print("$(lpad(round(g.SS_vec[i].shearMod,digits=5),13," ")) ") + print("$(lpad(round(g.SS_vec[i].Vp,digits=5),10," ")) ") + print("$(lpad(round(g.SS_vec[i].Vs,digits=5),10," ")) ") + print("\n") end for i=1:g.n_PP - print("$(lpad(g.ph[i],15," ")) ") - print("$(lpad(round(g.ph_frac[i],digits=5),13," ")) ") - print("$(lpad(round(g.ph_frac_wt[i],digits=5),8," ")) ") - print("$(lpad(round(g.PP_vec[i].f,digits=5),8," ")) ") - print("$(lpad(round(g.PP_vec[i].G,digits=5),8," ")) ") - print("$(lpad(round(g.PP_vec[i].V,digits=5),8," ")) ") - print("$(lpad(round(g.PP_vec[i].cp,digits=5),8," ")) ") - print("$(lpad(round(g.PP_vec[i].rho,digits=5),11," ")) ") - print("$(lpad(round(g.PP_vec[i].alpha,digits=5),8," ")) ") - print("$(lpad(round(g.PP_vec[i].entropy,digits=5),8," ")) ") - print("$(lpad(round(g.PP_vec[i].enthalpy,digits=5),8," ")) ") - print("$(lpad(round(g.PP_vec[i].bulkMod,digits=5),12," ")) ") - print("$(lpad(round(g.PP_vec[i].shearMod,digits=5),13," ")) ") - print("$(lpad(round(g.PP_vec[i].Vp,digits=5),10," ")) ") - print("$(lpad(round(g.PP_vec[i].Vs,digits=5),10," ")) ") - + print("$(lpad(g.ph[i],15," ")) ") + print("$(lpad(round(g.ph_frac[i],digits=5),13," ")) ") + print("$(lpad(round(g.ph_frac_wt[i],digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].f,digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].G,digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].V,digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].cp,digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].rho,digits=5),11," ")) ") + print("$(lpad(round(g.PP_vec[i].alpha,digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].entropy,digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].enthalpy,digits=5),8," ")) ") + print("$(lpad(round(g.PP_vec[i].bulkMod,digits=5),12," ")) ") + print("$(lpad(round(g.PP_vec[i].shearMod,digits=5),13," ")) ") + print("$(lpad(round(g.PP_vec[i].Vp,digits=5),10," ")) ") + print("$(lpad(round(g.PP_vec[i].Vs,digits=5),10," ")) ") + print("\n") end - - print("$(lpad("SYS",15," ")) ") - print("$(lpad(round(sum(g.ph_frac),digits=5),13," ")) ") - print("$(lpad(round(sum(g.ph_frac_wt),digits=5),8," ")) ") - print("$(lpad(round(g.G_system,digits=5),20," ")) ") - print("$(lpad(round(g.rho,digits=5),29," ")) ") - print("$(lpad(round(g.entropy,digits=5),21," ")) ") - print("$(lpad(round(g.enthalpy,digits=5),13," ")) ") - print("$(lpad(round(g.bulkMod,digits=5),13," ")) ") - print("$(lpad(round(g.shearMod,digits=5),13," ")) ") - print("$(lpad(round(g.Vp,digits=5),10," ")) ") - print("$(lpad(round(g.Vs,digits=5),10," ")) ") + + print("$(lpad("SYS",15," ")) ") + print("$(lpad(round(sum(g.ph_frac),digits=5),13," ")) ") + print("$(lpad(round(sum(g.ph_frac_wt),digits=5),8," ")) ") + print("$(lpad(round(g.G_system,digits=5),20," ")) ") + print("$(lpad(round(g.rho,digits=5),29," ")) ") + print("$(lpad(round(g.entropy,digits=5),21," ")) ") + print("$(lpad(round(g.enthalpy,digits=5),13," ")) ") + print("$(lpad(round(g.bulkMod,digits=5),13," ")) ") + print("$(lpad(round(g.shearMod,digits=5),13," ")) ") + print("$(lpad(round(g.Vp,digits=5),10," ")) ") + print("$(lpad(round(g.Vs,digits=5),10," ")) ") print("\n") print("\n") # == - println("Gamma (chemical potential of oxides):") + println("Gamma (chemical potential of oxides):") for i=1:length(g.oxides) - println(" $(lpad(g.oxides[i],6," ")) $(rpad(round(g.Gamma[i],digits=5),15," ")) ") + println(" $(lpad(g.oxides[i],6," ")) $(rpad(round(g.Gamma[i],digits=5),15," ")) ") end print("\n") - println("delta Gibbs energy (G-hyperplane distance):") + println("delta Gibbs energy (G-hyperplane distance):") for i=1:g.n_SS - println(" $(lpad(g.ph[i],6," ")) $(rpad(g.SS_vec[i].deltaG,15," ")) ") + println(" $(lpad(g.ph[i],6," ")) $(rpad(g.SS_vec[i].deltaG,15," ")) ") end for i=1:g.n_PP - println(" $(lpad(g.ph[i],6," ")) $(rpad(g.PP_vec[i].deltaG,15," ")) ") + println(" $(lpad(g.ph[i],6," ")) $(rpad(g.PP_vec[i].deltaG,15," ")) ") end - + print("\n") - println("mass residual : $(g.bulk_res_norm)") - println("# iterations : $(g.iter)") - println("status : $(g.status)") - + println("mass residual : $(g.bulk_res_norm)") + println("# iterations : $(g.iter)") + println("status : $(g.status)") + print("\n") diff --git a/test/gen_tests.jl b/test/gen_tests.jl index e74b7ee6..d5915b41 100755 --- a/test/gen_tests.jl +++ b/test/gen_tests.jl @@ -1,21 +1,21 @@ -# This script helps to generate a lsit of points for testing MAGEMin using reference built-in bulk-rock compositions +# This script helps to generate a list of points for testing MAGEMin using reference built-in bulk-rock compositions -cur_dir = pwd(); +cur_dir = pwd(); if cur_dir[end-3:end]=="test" cd("../") # change to main directory if we are in /test end using MAGEMin_C -# Initialize database +# Initialize database db = "ig" # database: ig, igneous (Holland et al., 2018); mp, metapelite (White et al 2014b) gv, z_b, DB, splx_data = init_MAGEMin(db); sys_in = "mol" #default is mol, if wt is provided conversion will be done internally (MAGEMin works on mol basis) -mutable struct outP{ _T } +mutable struct outP{ _T } P :: _T - T :: _T + T :: _T test :: Int64 G :: _T diff --git a/test/runtests.jl b/test/runtests.jl index ca580763..7ff2b7fa 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,217 +1,34 @@ # this tests the julia interface to MAGEMin using Test -cur_dir = pwd(); +# MAGEMin_C loads a locally compiled version of the library if it detects +# a file with the appropriate name in the current working directory. +# To make use of this also for tests, we change the working directory to +# the parent directory of the package iff necessary. +cur_dir = pwd() -if cur_dir[end-3:end]=="test" +if endswith(cur_dir, "test") cd("../") # change to main directory if we are in /test end -using MAGEMin_C # load MAGEMin (needs to be loaded from main directory to pick up correct library in case it is locally compiled) -# Initialize database - new way -data = Initialize_MAGEMin("ig", verbose=true); -test = 0 #KLB1 -data = use_predefined_bulk_rock(data, test); -# Call optimization routine for given P & T & bulk_rock -P = 8.0 -T = 800.0 -out = point_wise_minimization(P,T, data); - -@show out - -@test out.G_system ≈ -797.7491824683576 -@test out.ph == ["opx", "ol", "cpx", "spn"] -@test all(abs.(out.ph_frac - [ 0.24226960158631541, 0.5880694152724345, 0.1416697366114075, 0.027991246529842587]) .< 1e-4) - -# print more detailed info about this point: -print_info(out) -Finalize_MAGEMin(data) - - - -# previous way we defined this (left here for backwards compatibility) -db = "ig" -gv, z_b, DB, splx_data = init_MAGEMin(db); -sys_in = "mol" #default is mol, if wt is provided conversion will be done internally (MAGEMin works on mol basis) -test = 0 #KLB1 -gv = use_predefined_bulk_rock(gv, test, db); -gv.verbose=-1 -P = 8.0 -T = 800.0 -out = point_wise_minimization(P,T, gv, z_b, DB, splx_data, sys_in); -@test out.G_system ≈ -797.7491824683576 -@test out.ph == ["opx", "ol", "cpx", "spn"] -@test all(abs.(out.ph_frac - [ 0.24226960158631541, 0.5880694152724345, 0.1416697366114075, 0.027991246529842587]) .< 1e-4) -finalize_MAGEMin(gv,DB) - - -@testset "pointwise tests " begin - n = 100; - P = fill(8.0,n) - T = fill(800.0,n) - db = "ig" - data = Initialize_MAGEMin(db, verbose=false); - out = multi_point_minimization(P, T, data, test=0); - @test out[end].G_system ≈ -797.7491824683576 - @test out[end].ph == ["opx", "ol", "cpx", "spn"] - @test all(abs.(out[end].ph_frac - [ 0.24226960158631541, 0.5880694152724345, 0.1416697366114075, 0.027991246529842587]) .< 1e-4) - - Finalize_MAGEMin(data) -end - -@testset "specify bulk rock" begin - - - data = Initialize_MAGEMin("ig", verbose=false); - - # One bulk rock for all points - P,T = 10.0, 1100.0 - Xoxides = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "Cr2O3"; "H2O"]; - X = [48.43; 15.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; - sys_in = "wt" - out = single_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in) - - @test abs(out.G_system + 917.008526)/abs(917.008526) < 2e-4 - - - # different bulk rock per point - P = [10.0, 10.0] - T = [1100.0, 1100.0] - Xoxides = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "Cr2O3"; "H2O"]; - X1 = [48.43; 15.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; - X2 = [49.43; 14.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; - X = [X1,X2] - sys_in = "wt" - out = multi_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in) - - @test out[1].G_system ≈ -917.0085258258146 rtol=2e-4 - @test out[2].G_system ≈ -912.7754865001111 rtol=2e-4 - - Finalize_MAGEMin(data) -end - -@testset "convert bulk rock" begin - - bulk_in_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "Cr2O3"; "H2O"]; - bulk_in = [48.43; 15.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; - bulk_rock = convertBulk4MAGEMin(bulk_in,bulk_in_ox,"wt","ig" ); - - @test bulk_rock ≈ [46.12136551106488,8.52404156868422,11.80437413592006,14.382090296468109,6.470772496142926,0.3583593339444149,1.7262657202608955,0.487066733471891,0.5876026409473947,0.009999000099989998,9.528062562995213] - - bulk_in_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "MnO"; "H2O"]; - bulk_in = [69.64; 13.76; 1.77; 1.73; 4.32; 0.4; 2.61; 2.41; 0.80; 0.07; 0.0]; - bulk_rock = convertBulk4MAGEMin(bulk_in,bulk_in_ox,"wt","mp" ); - - @test bulk_rock ≈ [76.57038397179574,8.914984523583415,2.0849576977131403,2.835783318610597,4.30275071755529,1.8302970975627948,2.568605789798099,0.6615823604771729,0.16546809116073818,0.06518643174302832,0.0] +@testset "serial" begin + include(joinpath(@__DIR__, "tests.jl")) end -@testset "test Seismic velocities & modulus" begin - # Call optimization routine for given P & T & bulk_rock - data = Initialize_MAGEMin("ig", verbose=false); - test = 0; - data = use_predefined_bulk_rock(data, test) - P = 8.0 - T = 1200.0 - out = point_wise_minimization(P,T, data) +@testset "threaded" begin + # Do a dummy `@test true`: + # If the process errors out the testset would error out as well + @test true - tol = 5e-2; - @test abs(out.bulkMod - 95.35222421341481 < tol) - @test abs(out.shearMod - 29.907907390690557 < tol) - @test abs(out.Vs - 3.056253320843246 < tol) - @test abs(out.Vp - 6.498781717400121 < tol) - @test abs(out.Vs_S - 4.30872049030154 < tol) - @test abs(out.Vp_S - 7.392153167537697 < tol) - # @test abs(out.bulkModulus_M - 27.260603902167567 < tol) - @test abs(out.bulkModulus_S - 95.74343528580735 < tol) - @test abs(out.shearModulus_S - 59.4665150508297 < tol) - - Finalize_MAGEMin(data) -end - -# Stores data of tests -mutable struct outP{ _T } - P :: _T - T :: _T - test :: Int64 - - G :: _T - ph :: Vector{String} - ph_frac :: Vector{Float64} -end - -print_error_msg(i,out) = println("ERROR for point $i with test=$(out.test); P=$(out.P); T=$(out.T); stable phases=$(out.ph), fractions=$(out.ph_frac)") - - -# Automatic testing of all points -function TestPoints(list, data::MAGEMin_Data) - - # Compute all points - P = [ l.P for l in list] - T = [ l.T for l in list] - test = [ l.test for l in list] - out_vec = multi_point_minimization(P, T, data, test = test[1]); - - # Check if the points this fit - for (i,out) in enumerate(out_vec) - VerifyPoint(out, list[i], i) + # Only run additional tests if we are running with a single thread right now + if Threads.nthreads() == 1 + # We explicitly disable code coverage tracking with multiple threads since + # this is expensive, see https://github.com/JuliaLang/julia/issues/36142 + run(`$(Base.julia_cmd()) --threads=2 --check-bounds=yes --code-coverage=none $(abspath(joinpath(@__DIR__, "tests.jl")))`) end - return nothing end -# This checks whether a single point agrees with precomputed values & prints a message if not -function VerifyPoint(out, list, i) - - # We need to sort the phases (sometimes they are ordered differently) - ind_sol = sortperm(list.ph) - ind_out = sortperm(out.ph) - - result1 = @test out.G_system ≈ list.G rtol=1e-3 - result2 = @test out.ph[ind_out] == list.ph[ind_sol] - result3 = @test out.ph_frac[ind_out] ≈ list.ph_frac[ind_sol] atol=5e-2 # ok, this is really large (needs fixing for test6!) - - # print more info about the point if one of the tests above fails - if isa(result1,Test.Fail) || isa(result2,Test.Fail) || isa(result3,Test.Fail) - print_error_msg(i,list) - end - - return nothing -end - -# load reference for built-in tests -println("Testing points from the reference diagrams:") -@testset verbose = true "Total tests" begin - println(" Starting KLB-1 peridotite tests") - db = "ig" # database: ig, igneous (Holland et al., 2018); mp, metapelite (White et al 2014b) - data = Initialize_MAGEMin(db, verbose=false); - - gv.verbose=-1; - @testset "KLB-1 peridotite tests" begin - include("test_diagram_test0.jl") - TestPoints(list, data) - end - Finalize_MAGEMin(data) - - println(" Starting RE-46 icelandic basalt tests") - db = "ig" # database: ig, igneous (Holland et al., 2018); mp, metapelite (White et al 2014b) - data = Initialize_MAGEMin(db, verbose=false); - gv.verbose=-1; - @testset "RE-46 icelandic basalt tests" begin - include("test_diagram_test1.jl") - TestPoints(list, data) - end - - - println(" Starting Wet MORB tests") - db = "ig" # database: ig, igneous (Holland et al., 2018); mp, metapelite (White et al 2014b) - data = Initialize_MAGEMin(db, verbose=false); - @testset "Wet MORB tests" begin - include("test_diagram_test6.jl") - TestPoints(list, data) - end - Finalize_MAGEMin(data) -end - - +# Change back to the working directory we used when starting to run the tests cd(cur_dir) diff --git a/test/tests.jl b/test/tests.jl new file mode 100644 index 00000000..40b06e2d --- /dev/null +++ b/test/tests.jl @@ -0,0 +1,222 @@ +using Test + +# Load MAGEMin (needs to be loaded from main directory to pick up correct +# library in case it is locally compiled). This is handled by the logic in +# runtests.jl +using MAGEMin_C + +# Initialize database - new way +data = Initialize_MAGEMin("ig", verbose=true); +test = 0 #KLB1 +data = use_predefined_bulk_rock(data, test); + +# Call optimization routine for given P & T & bulk_rock +P = 8.0 +T = 800.0 +out = point_wise_minimization(P,T, data); + +@show out + +@test out.G_system ≈ -797.7491824683576 +@test out.ph == ["opx", "ol", "cpx", "spn"] +@test all(abs.(out.ph_frac - [ 0.24226960158631541, 0.5880694152724345, 0.1416697366114075, 0.027991246529842587]) .< 1e-4) + +# print more detailed info about this point: +print_info(out) +Finalize_MAGEMin(data) + + + +# previous way we defined this (left here for backwards compatibility) +db = "ig" +gv, z_b, DB, splx_data = init_MAGEMin(db); +sys_in = "mol" #default is mol, if wt is provided conversion will be done internally (MAGEMin works on mol basis) +test = 0 #KLB1 +gv = use_predefined_bulk_rock(gv, test, db); +gv.verbose=-1 +P = 8.0 +T = 800.0 +out = point_wise_minimization(P,T, gv, z_b, DB, splx_data, sys_in); +@test out.G_system ≈ -797.7491824683576 +@test out.ph == ["opx", "ol", "cpx", "spn"] +@test all(abs.(out.ph_frac - [ 0.24226960158631541, 0.5880694152724345, 0.1416697366114075, 0.027991246529842587]) .< 1e-4) +finalize_MAGEMin(gv,DB) + + +@testset "pointwise tests " begin + n = 100; + P = fill(8.0,n) + T = fill(800.0,n) + db = "ig" + data = Initialize_MAGEMin(db, verbose=false); + out = multi_point_minimization(P, T, data, test=0); + @test out[end].G_system ≈ -797.7491824683576 + @test out[end].ph == ["opx", "ol", "cpx", "spn"] + @test all(abs.(out[end].ph_frac - [ 0.24226960158631541, 0.5880694152724345, 0.1416697366114075, 0.027991246529842587]) .< 1e-4) + + Finalize_MAGEMin(data) +end + + +@testset "ultramafic database" begin + P,T = 10.0, 500.0 + data = Initialize_MAGEMin("um", verbose=false); + out = single_point_minimization(P, T, data) + @test abs(out.G_system + 556.243321)/abs(556.243321) < 2e-4 + @test all(abs.(out.ph_frac - [0.5753192548298169, 0.024582756340200455, 0.024184726745689656, 0.0554904375970385, 0.0016172391886508995, 0.3184197873654363, 0.0003857953985807551]) .< 1e-4) + @test out.ph == ["atg", "ol", "spi", "chl", "po", "fluid", "pyr"] + +end + + +@testset "specify bulk rock" begin + + data = Initialize_MAGEMin("ig", verbose=false); + + # One bulk rock for all points + P,T = 10.0, 1100.0 + Xoxides = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "Cr2O3"; "H2O"]; + X = [48.43; 15.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; + sys_in = "wt" + out = single_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in) + + @test abs(out.G_system + 917.008526)/abs(917.008526) < 2e-4 + + + # different bulk rock per point + P = [10.0, 10.0] + T = [1100.0, 1100.0] + Xoxides = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "Cr2O3"; "H2O"]; + X1 = [48.43; 15.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; + X2 = [49.43; 14.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; + X = [X1,X2] + sys_in = "wt" + out = multi_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in) + + @test out[1].G_system ≈ -917.0085258258146 rtol=2e-4 + @test out[2].G_system ≈ -912.7754865001111 rtol=2e-4 + + Finalize_MAGEMin(data) +end + +@testset "convert bulk rock" begin + + bulk_in_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "Cr2O3"; "H2O"]; + bulk_in = [48.43; 15.19; 11.57; 10.13; 6.65; 1.64; 0.59; 1.87; 0.68; 0.0; 3.0]; + bulk_rock = convertBulk4MAGEMin(bulk_in,bulk_in_ox,"wt","ig" ); + + @test bulk_rock ≈ [46.12136551106488,8.52404156868422,11.80437413592006,14.382090296468109,6.470772496142926,0.3583593339444149,1.7262657202608955,0.487066733471891,0.5876026409473947,0.009999000099989998,9.528062562995213] + + bulk_in_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "MnO"; "H2O"]; + bulk_in = [69.64; 13.76; 1.77; 1.73; 4.32; 0.4; 2.61; 2.41; 0.80; 0.07; 0.0]; + bulk_rock = convertBulk4MAGEMin(bulk_in,bulk_in_ox,"wt","mp" ); + + @test bulk_rock ≈ [76.57038397179574,8.914984523583415,2.0849576977131403,2.835783318610597,4.30275071755529,1.8302970975627948,2.568605789798099,0.6615823604771729,0.16546809116073818,0.06518643174302832,0.0] +end + + +@testset "Seismic velocities & modulus" begin + # Call optimization routine for given P & T & bulk_rock + data = Initialize_MAGEMin("ig", verbose=false); + test = 0; + data = use_predefined_bulk_rock(data, test) + P = 8.0 + T = 1200.0 + out = point_wise_minimization(P,T, data) + + tol = 5e-2; + @test abs(out.bulkMod - 95.35222421341481 < tol) + @test abs(out.shearMod - 29.907907390690557 < tol) + @test abs(out.Vs - 3.056253320843246 < tol) + @test abs(out.Vp - 6.498781717400121 < tol) + @test abs(out.Vs_S - 4.30872049030154 < tol) + @test abs(out.Vp_S - 7.392153167537697 < tol) + # @test abs(out.bulkModulus_M - 27.260603902167567 < tol) + @test abs(out.bulkModulus_S - 95.74343528580735 < tol) + @test abs(out.shearModulus_S - 59.4665150508297 < tol) + + Finalize_MAGEMin(data) +end + +# Stores data of tests +mutable struct outP{ _T } + P :: _T + T :: _T + test :: Int64 + + G :: _T + ph :: Vector{String} + ph_frac :: Vector{Float64} +end + +print_error_msg(i,out) = println("ERROR for point $i with test=$(out.test); P=$(out.P); T=$(out.T); stable phases=$(out.ph), fractions=$(out.ph_frac)") + + +# Automatic testing of all points +function TestPoints(list, data::MAGEMin_Data) + + # Compute all points + P = [ l.P for l in list] + T = [ l.T for l in list] + test = [ l.test for l in list] + out_vec = multi_point_minimization(P, T, data, test = test[1]); + + # Check if the points this fit + for (i,out) in enumerate(out_vec) + VerifyPoint(out, list[i], i) + end + return nothing +end + +# This checks whether a single point agrees with precomputed values & prints a message if not +function VerifyPoint(out, list, i) + + # We need to sort the phases (sometimes they are ordered differently) + ind_sol = sortperm(list.ph) + ind_out = sortperm(out.ph) + + result1 = @test out.G_system ≈ list.G rtol=1e-3 + result2 = @test out.ph[ind_out] == list.ph[ind_sol] + result3 = @test out.ph_frac[ind_out] ≈ list.ph_frac[ind_sol] atol=5e-2 # ok, this is really large (needs fixing for test6!) + + # print more info about the point if one of the tests above fails + if isa(result1,Test.Fail) || isa(result2,Test.Fail) || isa(result3,Test.Fail) + print_error_msg(i,list) + end + + return nothing +end + +# load reference for built-in tests +println("Testing points from the reference diagrams:") +@testset verbose = true "Total tests" begin + println(" Starting KLB-1 peridotite tests") + db = "ig" # database: ig, igneous (Holland et al., 2018); mp, metapelite (White et al 2014b) + data = Initialize_MAGEMin(db, verbose=false); + + gv.verbose=-1; + @testset "KLB-1 peridotite tests" begin + include("test_diagram_test0.jl") + TestPoints(list, data) + end + Finalize_MAGEMin(data) + + println(" Starting RE-46 icelandic basalt tests") + db = "ig" # database: ig, igneous (Holland et al., 2018); mp, metapelite (White et al 2014b) + data = Initialize_MAGEMin(db, verbose=false); + gv.verbose=-1; + @testset "RE-46 icelandic basalt tests" begin + include("test_diagram_test1.jl") + TestPoints(list, data) + end + + + println(" Starting Wet MORB tests") + db = "ig" # database: ig, igneous (Holland et al., 2018); mp, metapelite (White et al 2014b) + data = Initialize_MAGEMin(db, verbose=false); + @testset "Wet MORB tests" begin + include("test_diagram_test6.jl") + TestPoints(list, data) + end + Finalize_MAGEMin(data) +end diff --git a/utilities/indent b/utilities/indent deleted file mode 100755 index 54cf1af5..00000000 --- a/utilities/indent +++ /dev/null @@ -1,69 +0,0 @@ -#!/bin/bash - -# collect all header and source files and process them in batches of 50 files -# with up to 10 in parallel - -# remove execute permission on source files: -find src \( -name '*.c' -o -name '*.h' \) -print | xargs -n 50 -P 10 chmod -x -find julia \( -name '*.jl' \) -print | xargs -n 50 -P 10 chmod -x -find matlab \( -name '*.m' \) -print | xargs -n 50 -P 10 chmod -x - -# convert tabs to four spaces: -tab_to_space() -{ - f=$1 - # awkward tab replacement because of OSX sed, do not change unless you test it on OSX - TAB=$'\t' - sed -e "s/$TAB/ /g" $f >$f.tmp - diff -q $f $f.tmp >/dev/null || mv $f.tmp $f - rm -f $f.tmp -} -export -f tab_to_space -find src \( -name '*.c' -o -name '*.h' \) -print | xargs -n 1 -P 10 -I {} bash -c 'tab_to_space "$@"' _ {} -find julia \( -name '*.jl' \) -print | xargs -n 1 -P 10 -I {} bash -c 'tab_to_space "$@"' _ {} -find matlab \( -name '*.m' \) -print | xargs -n 1 -P 10 -I {} bash -c 'tab_to_space "$@"' _ {} - -# Remove trailing whitespace from files -remove_trailing_whitespace() -{ - f=$1 - # awkward tab replacement because of OSX sed, do not change unless you test it on OSX - TAB=$'\t' - sed -e "s/[ $TAB]*$//" $f >$f.tmp - diff -q $f $f.tmp >/dev/null || mv $f.tmp $f - rm -f $f.tmp -} -export -f remove_trailing_whitespace -find src \( -name '*.c' -o -name '*.h' \) -print | xargs -n 1 -P 10 -I {} bash -c 'remove_trailing_whitespace "$@"' _ {} -find julia \( -name '*.jl' \) -print | xargs -n 1 -P 10 -I {} bash -c 'remove_trailing_whitespace "$@"' _ {} -find matlab \( -name '*.m' \) -print | xargs -n 1 -P 10 -I {} bash -c 'remove_trailing_whitespace "$@"' _ {} - -# Ensure only a single newline at end of files -ensure_single_trailing_newline() -{ - f=$1 - - # Remove newlines at end of file - # Check that the current line only contains newlines - # If it doesn't match, print it - # If it does match and we're not at the end of the file, - # append the next line to the current line and repeat the check - # If it does match and we're at the end of the file, - # remove the line. - sed -e :a -e '/^\n*$/{$d;N;};/\n$/ba' $f >$f.tmpi - - # Then add a newline to the end of the file - # '$' denotes the end of file - # 'a\' appends the following text (which in this case is nothing) - # on a new line - sed -e '$a\' $f.tmpi >$f.tmp - - diff -q $f $f.tmp >/dev/null || mv $f.tmp $f - rm -f $f.tmp $f.tmpi -} -export -f ensure_single_trailing_newline -find src \( -name '*.c' -o -name '*.h' \) -print | xargs -n 1 -P 10 -I {} bash -c 'ensure_single_trailing_newline "$@"' _ {} -find julia \( -name '*.jl' \) -print | xargs -n 1 -P 10 -I {} bash -c 'ensure_single_trailing_newline "$@"' _ {} -find matlab \( -name '*.m' \) -print | xargs -n 1 -P 10 -I {} bash -c 'ensure_single_trailing_newline "$@"' _ {} - -exit 0