diff --git a/Project.toml b/Project.toml index a487641..70b6616 100644 --- a/Project.toml +++ b/Project.toml @@ -1,10 +1,13 @@ name = "MAGEMin_C" uuid = "e5d170eb-415a-4524-987b-12f1bce1ddab" authors = ["Boris Kaus & Nicolas Riel "] -version = "1.4.7" +version = "1.4.8" [deps] CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82" +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" MAGEMin_jll = "763ebaa8-b0d2-5f6b-90ef-4fc23b5db1c4" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" @@ -12,6 +15,8 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] CEnum = "0.4" -MAGEMin_jll = "1.4.4 - 1.4.4" +DataFrames = "1.6.1" +CSV = "0.10.14" +MAGEMin_jll = "1.4.5 - 1.4.5" ProgressMeter = "1" julia = "1.6" diff --git a/gen/magemin_library.jl b/gen/magemin_library.jl index a6bc1b7..d42b796 100644 --- a/gen/magemin_library.jl +++ b/gen/magemin_library.jl @@ -582,6 +582,7 @@ mutable struct global_variables numPoint::Cint global_ite::Cint H2O_id::Cint + Al2O3_id::Cint K2O_id::Cint O_id::Cint TiO2_id::Cint @@ -1412,7 +1413,7 @@ mutable struct metapelite_datasets n_pp::Cint n_ss::Cint ox::NTuple{11, NTuple{20, Cchar}} - PP::NTuple{17, NTuple{20, Cchar}} + PP::NTuple{23, NTuple{20, Cchar}} SS::NTuple{16, NTuple{20, Cchar}} verifyPC::NTuple{16, Cint} n_SS_PC::NTuple{16, Cint} @@ -1438,7 +1439,7 @@ mutable struct metabasite_datasets n_pp::Cint n_ss::Cint ox::NTuple{10, NTuple{20, Cchar}} - PP::NTuple{18, NTuple{20, Cchar}} + PP::NTuple{24, NTuple{20, Cchar}} SS1::NTuple{14, NTuple{20, Cchar}} verifyPC1::NTuple{14, Cint} n_SS_PC1::NTuple{14, Cint} @@ -1468,7 +1469,7 @@ mutable struct igneous_datasets n_pp::Cint n_ss::Cint ox::NTuple{11, NTuple{20, Cchar}} - PP::NTuple{17, NTuple{20, Cchar}} + PP::NTuple{23, NTuple{20, Cchar}} SS::NTuple{15, NTuple{20, Cchar}} verifyPC::NTuple{15, Cint} n_SS_PC::NTuple{15, Cint} @@ -1494,7 +1495,7 @@ mutable struct ultramafic_datasets n_pp::Cint n_ss::Cint ox::NTuple{7, NTuple{20, Cchar}} - PP::NTuple{15, NTuple{20, Cchar}} + PP::NTuple{21, NTuple{20, Cchar}} SS::NTuple{12, NTuple{20, Cchar}} verifyPC::NTuple{12, Cint} n_SS_PC::NTuple{12, Cint} diff --git a/julia/MAGEMin_wrappers.jl b/julia/MAGEMin_wrappers.jl index f82dbe3..a85cbc4 100644 --- a/julia/MAGEMin_wrappers.jl +++ b/julia/MAGEMin_wrappers.jl @@ -3,6 +3,7 @@ import Base.show using Base.Threads: @threads using ProgressMeter +using DataFrames, Dates, CSV const VecOrMat = Union{Nothing, AbstractVector{Float64}, AbstractVector{<:AbstractVector{Float64}}} @@ -10,10 +11,111 @@ export retrieve_solution_phase_information, remove_phases, init_MAGEMin, finalize_MAGEMin, point_wise_minimization, convertBulk4MAGEMin, use_predefined_bulk_rock, define_bulk_rock, create_output, print_info, create_gmin_struct, pwm_init, pwm_run, single_point_minimization, multi_point_minimization, MAGEMin_Data, W_Data, + MAGEMin_data2dataframe, MAGEMin_dataTE2dataframe, Initialize_MAGEMin, Finalize_MAGEMin -export get_TE_database, compute_TE_partitioning, zirconium_saturation - +export adjust_chemical_system, TE_prediction, get_OL_KDs_database, adjust_bulk_4_zircon + + + +""" + structure that holds the result of the pointwise minimization +""" +struct gmin_struct{T,I} + MAGEMin_ver :: String + dataset :: String + G_system :: T # G of system + Gamma :: Vector{T} # Gamma + P_kbar :: T # Pressure in kbar + T_C :: T # Temperature in Celsius + X :: Vector{T} + M_sys :: T + + # 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} + + # Fractions: + # Solid, melt, fluid fractions + frac_M :: T + frac_S :: T + frac_F :: T + + frac_M_wt :: T + frac_S_wt :: T + frac_F_wt :: T + + # Solid, melt, fluid densities + alpha :: T + V :: T + cp :: T + s_cp :: Vector{T} + rho :: T + rho_M :: T + rho_S :: T + rho_F :: T + + # Oxygen fugacity + fO2 :: T + dQFM :: T + + # Activities + aH2O :: T + aSiO2 :: T + aTiO2 :: T + aAl2O3 :: T + aMgO :: T + aFeO :: T + + # Phase fractions and type: + n_PP :: Int64 # number of pure phases + n_SS :: Int64 # number of solid solutions + n_mSS :: Int64 # number of solid solutions + + ph_frac :: Vector{T} # phase fractions + ph_frac_wt :: Vector{T} # phase fractions + ph_frac_vol :: 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} + mSS_vec :: Vector{LibMAGEMin.mSS_data} + PP_vec :: Vector{LibMAGEMin.PP_data} + + oxides :: Vector{String} + + # Seismic velocity info + Vp :: T # P-wave velocity + Vs :: T # S-wave velocity + Vp_S :: T # P-wave velocity of solid aggregate + Vs_S :: T # S-wave velocity of solid aggregate + bulkMod :: T # Elastic bulk modulus + shearMod :: T # Elastic shear modulus + bulkModulus_M :: T # Elastic bulk modulus + bulkModulus_S :: T # Elastic bulk modulus + shearModulus_S :: T # Elastic shear modulus + + # thermodynamic properties + entropy :: T # entropy + enthalpy :: T # enthalpy + + # Numerics: + iter :: I # number of iterations required + bulk_res_norm :: T # bulk residual norm + time_ms :: T # computational time for this point + status :: I # status of calculations +end + + + """ Holds general information about solution phases """ @@ -64,7 +166,7 @@ end """ function retrieve_solution_phase_information(dtb) - db_inf = db_infos[db_infos("mp", "Metapelite (White et al., 2014)", ss_infos[ss_infos("liq", 8, 7, ["none", "q4L", "abL", "kspL", "anL", "slL", "fo2L", "fa2L", "h2oL"], ["none", "q", "fsp", "na", "an", "ol", "x", "h2o"]), ss_infos("fsp", 3, 2, ["none", "ab", "an", "san"], ["none", "ca", "k"]), ss_infos("bi", 7, 6, ["none", "phl", "annm", "obi", "east", "tbi", "fbi", "mmbi"], ["none", "x", "m", "y", "f", "t", "Q"]), ss_infos("g", 5, 4, ["none", "py", "alm", "spss", "gr", "kho"], ["none", "x", "z", "m", "f"]), ss_infos("ep", 3, 2, ["none", "cz", "ep", "fep"], ["none", "f", "Q"]), ss_infos("ma", 6, 5, ["none", "mut", "celt", "fcelt", "pat", "ma", "fmu"], ["none", "x", "y", "f", "n", "c"]), ss_infos("mu", 6, 5, ["none", "mut", "cel", "fcel", "pat", "ma", "fmu"], ["none", "x", "y", "f", "n", "c"]), ss_infos("opx", 7, 6, ["none", "en", "fs", "fm", "mgts", "fopx", "mnopx", "odi"], ["none", "x", "m", "y", "f", "c", "Q"]), ss_infos("sa", 5, 4, ["none", "spr4", "spr5", "fspm", "spro", "ospr"], ["none", "x", "y", "f", "Q"]), ss_infos("cd", 4, 3, ["none", "crd", "fcrd", "hcrd", "mncd"], ["none", "x", "m", "h"]), ss_infos("st", 5, 4, ["none", "mstm", "fst", "mnstm", "msto", "mstt"], ["none", "x", "m", "f", "t"]), ss_infos("chl", 8, 7, ["none", "clin", "afchl", "ames", "daph", "ochl1", "ochl4", "f3clin", "mmchl"], ["none", "x", "y", "f", "m", "QAl", "Q1", "Q4"]), ss_infos("ctd", 4, 3, ["none", "mctd", "fctd", "mnct", "ctdo"], ["none", "x", "m", "f"]), ss_infos("sp", 4, 3, ["none", "herc", "sp", "mt", "usp"], ["none", "x", "y", "z"]), ss_infos("ilmm", 5, 4, ["none", "oilm", "dilm", "dhem", "geik", "pnt"], ["none", "i", "g", "m", "Q"])], ["liq", "fsp", "bi", "g", "ep", "ma", "mu", "opx", "sa", "cd", "st", "chl", "ctd", "sp", "ilmm"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and", "ru", "sph", "O2", "H2O", "qfm", "qif", "nno", "hm", "cco"]), db_infos("mb", "Metabasite (Green et al., 2016)", ss_infos[ss_infos("sp", 4, 3, ["none", "herc", "sp", "mt", "usp"], ["none", "x", "y", "z"]), ss_infos("opx", 6, 5, ["none", "en", "fs", "fm", "mgts", "fopx", "odi"], ["none", "x", "y", "f", "c", "Q"]), ss_infos("fsp", 3, 2, ["none", "ab", "an", "san"], ["none", "ca", "k"]), ss_infos("liq", 9, 8, ["none", "q4L", "abL", "kspL", "wo1L", "sl1L", "fa2L", "fo2L", "watL", "anoL"], ["none", "q", "fsp", "na", "wo", "sil", "ol", "x", "yan"]), ss_infos("mu", 6, 5, ["none", "mu", "cel", "fcel", "pa", "mam", "fmu"], ["none", "x", "y", "f", "n", "c"]), ss_infos("ilm", 3, 2, ["none", "oilm", "dilm", "dhem"], ["none", "x", "Q"]), ss_infos("ol", 2, 1, ["none", "fo", "fa"], ["none", "x"]), ss_infos("hb", 11, 10, ["none", "tr", "tsm", "prgm", "glm", "cumm", "grnm", "a", "b", "mrb", "kprg", "tts"], ["none", "x", "y", "z", "a", "k", "c", "f", "t", "Q1", "Q2"]), ss_infos("ep", 3, 2, ["none", "cz", "ep", "fep"], ["none", "f", "Q"]), ss_infos("g", 4, 3, ["none", "py", "alm", "gr", "kho"], ["none", "x", "z", "f"]), ss_infos("chl", 7, 6, ["none", "clin", "afchl", "ames", "daph", "ochl1", "ochl4", "f3clin"], ["none", "x", "y", "f", "QAl", "Q1", "Q4"]), ss_infos("bi", 6, 5, ["none", "phl", "annm", "obi", "east", "tbi", "fbi"], ["none", "x", "y", "f", "t", "Q"]), ss_infos("dio", 7, 6, ["none", "jd", "di", "hed", "acmm", "om", "cfm", "jac"], ["none", "x", "j", "t", "c", "Qaf", "Qfm"]), ss_infos("abc", 2, 1, ["none", "abm", "anm"], ["none", "ca"])], ["sp", "opx", "fsp", "liq", "mu", "ilm", "ol", "hb", "ep", "g", "chl", "bi", "dio", "abc"], ["q", "crst", "trd", "coe", "law", "ky", "sill", "and", "ru", "sph", "sph", "ab", "H2O", "qfm", "qif", "nno", "hm", "cco"]), db_infos("ig", "Igneous (Holland et al., 2018)", ss_infos[ss_infos("spn", 8, 7, ["none", "nsp", "isp", "nhc", "ihc", "nmt", "imt", "pcr", "qndm"], ["none", "x", "y", "c", "t", "Q1", "Q2", "Q3"]), ss_infos("bi", 6, 5, ["none", "phl", "annm", "obi", "eas", "tbi", "fbi"], ["none", "x", "y", "f", "t", "Q"]), ss_infos("cd", 3, 2, ["none", "crd", "fcrd", "hcrd"], ["none", "x", "h"]), ss_infos("cpx", 10, 9, ["none", "di", "cfs", "cats", "crdi", "cess", "cbuf", "jd", "cen", "cfm", "kjd"], ["none", "x", "y", "o", "n", "Q", "f", "cr", "t", "k"]), ss_infos("ep", 3, 2, ["none", "cz", "ep", "fep"], ["none", "f", "Q"]), ss_infos("g", 6, 5, ["none", "py", "alm", "gr", "andr", "knom", "tig"], ["none", "x", "c", "f", "cr", "t"]), ss_infos("hb", 11, 10, ["none", "tr", "tsm", "prgm", "glm", "cumm", "grnm", "a", "b", "mrb", "kprg", "tts"], ["none", "x", "y", "z", "a", "k", "c", "f", "t", "Q1", "Q2"]), ss_infos("ilm", 3, 2, ["none", "oilm", "dilm", "dhem"], ["none", "x", "Q"]), ss_infos("liq", 12, 11, ["none", "q4L", "slL", "wo1L", "fo2L", "fa2L", "jdL", "hmL", "ekL", "tiL", "kjL", "ctL", "wat1L"], ["none", "wo", "sl", "fo", "fa", "jd", "hm", "ek", "ti", "kj", "yct", "h2o"]), ss_infos("ol", 4, 3, ["none", "mont", "fa", "fo", "cfm"], ["none", "x", "c", "Q"]), ss_infos("opx", 9, 8, ["none", "en", "fs", "fm", "odi", "mgts", "cren", "obuf", "mess", "ojd"], ["none", "x", "y", "c", "Q", "f", "t", "cr", "j"]), ss_infos("fsp", 3, 2, ["none", "ab", "an", "san"], ["none", "ca", "k"]), ss_infos("fl", 11, 10, ["none", "qfL", "slfL", "wofL", "fofL", "fafL", "jdfL", "hmfL", "ekfL", "tifL", "kjfL", "H2O"], ["none", "wo", "sl", "fo", "fa", "jd", "hm", "ek", "ti", "kj", "h2o"]), ss_infos("mu", 6, 5, ["none", "mu", "cel", "fcel", "pa", "mam", "fmu"], ["none", "x", "y", "f", "n", "c"]), ss_infos("fper", 2, 1, ["none", "per", "wu"], ["none", ""])], ["spn", "bi", "cd", "cpx", "ep", "g", "hb", "ilm", "liq", "ol", "opx", "fsp", "fl", "mu", "fper"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and", "ru", "sph", "O2", "qfm", "mw", "qif", "nno", "hm", "cco"]), db_infos("um", "Ultramafic (Evans & Frost., 2021)", ss_infos[ss_infos("fl", 2, 1, ["none", "H2", "H2O"], ["none", "x"]), ss_infos("ol", 2, 1, ["none", "fo", "fa"], ["none", "x"]), ss_infos("br", 2, 1, ["none", "br", "fbr"], ["none", "x"]), ss_infos("ch", 2, 1, ["none", "chum", "chuf"], ["none", "x"]), ss_infos("atg", 5, 4, ["none", "atgf", "fatg", "atgo", "aatg", "oatg"], ["none", "x", "y", "f", "t"]), ss_infos("g", 2, 1, ["none", "py", "alm"], ["none", "x"]), ss_infos("ta", 6, 5, ["none", "ta", "fta", "tao", "tats", "ota", "tap"], ["none", "x", "y", "f", "v", "Q"]), ss_infos("chl", 7, 6, ["none", "clin", "afchl", "ames", "daph", "ochl1", "ochl4", "f3clin"], ["none", "x", "y", "f", "m", "t", "QA1"]), ss_infos("spi", 3, 2, ["none", "herc", "sp", "mt"], ["none", "x", "y"]), ss_infos("opx", 5, 4, ["none", "en", "fs", "fm", "mgts", "fopx"], ["none", "x", "y", "f", "Q"]), ss_infos("po", 2, 1, ["none", "trov", "trot"], ["none", "y"]), ss_infos("anth", 5, 4, ["none", "anth", "gedf", "fant", "a", "b"], ["none", "x", "y", "z", "a"])], ["fl", "ol", "br", "ch", "atg", "g", "ta", "chl", "spi", "opx", "po", "anth"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and", "pyr", "O2", "qfm", "qif", "nno", "hm", "cco"])] + db_inf = db_infos[db_infos("mp", "Metapelite (White et al., 2014)", ss_infos[ss_infos("liq", 8, 7, ["none", "q4L", "abL", "kspL", "anL", "slL", "fo2L", "fa2L", "h2oL"], ["none", "q", "fsp", "na", "an", "ol", "x", "h2o"]), ss_infos("fsp", 3, 2, ["none", "ab", "an", "san"], ["none", "ca", "k"]), ss_infos("bi", 7, 6, ["none", "phl", "annm", "obi", "east", "tbi", "fbi", "mmbi"], ["none", "x", "m", "y", "f", "t", "Q"]), ss_infos("g", 5, 4, ["none", "py", "alm", "spss", "gr", "kho"], ["none", "x", "z", "m", "f"]), ss_infos("ep", 3, 2, ["none", "cz", "ep", "fep"], ["none", "f", "Q"]), ss_infos("ma", 6, 5, ["none", "mut", "celt", "fcelt", "pat", "ma", "fmu"], ["none", "x", "y", "f", "n", "c"]), ss_infos("mu", 6, 5, ["none", "mut", "cel", "fcel", "pat", "ma", "fmu"], ["none", "x", "y", "f", "n", "c"]), ss_infos("opx", 7, 6, ["none", "en", "fs", "fm", "mgts", "fopx", "mnopx", "odi"], ["none", "x", "m", "y", "f", "c", "Q"]), ss_infos("sa", 5, 4, ["none", "spr4", "spr5", "fspm", "spro", "ospr"], ["none", "x", "y", "f", "Q"]), ss_infos("cd", 4, 3, ["none", "crd", "fcrd", "hcrd", "mncd"], ["none", "x", "m", "h"]), ss_infos("st", 5, 4, ["none", "mstm", "fst", "mnstm", "msto", "mstt"], ["none", "x", "m", "f", "t"]), ss_infos("chl", 8, 7, ["none", "clin", "afchl", "ames", "daph", "ochl1", "ochl4", "f3clin", "mmchl"], ["none", "x", "y", "f", "m", "QAl", "Q1", "Q4"]), ss_infos("ctd", 4, 3, ["none", "mctd", "fctd", "mnct", "ctdo"], ["none", "x", "m", "f"]), ss_infos("sp", 4, 3, ["none", "herc", "sp", "mt", "usp"], ["none", "x", "y", "z"]), ss_infos("ilmm", 5, 4, ["none", "oilm", "dilm", "dhem", "geik", "pnt"], ["none", "i", "g", "m", "Q"])], ["liq", "fsp", "bi", "g", "ep", "ma", "mu", "opx", "sa", "cd", "st", "chl", "ctd", "sp", "ilmm"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and", "ru", "sph", "O2", "H2O", "qfm", "qif", "nno", "hm", "cco", "aH2O", "aO2", "aMgO", "aFeO", "aAl2O3", "aTiO2"]), db_infos("mb", "Metabasite (Green et al., 2016)", ss_infos[ss_infos("sp", 4, 3, ["none", "herc", "sp", "mt", "usp"], ["none", "x", "y", "z"]), ss_infos("opx", 6, 5, ["none", "en", "fs", "fm", "mgts", "fopx", "odi"], ["none", "x", "y", "f", "c", "Q"]), ss_infos("fsp", 3, 2, ["none", "ab", "an", "san"], ["none", "ca", "k"]), ss_infos("liq", 9, 8, ["none", "q4L", "abL", "kspL", "wo1L", "sl1L", "fa2L", "fo2L", "watL", "anoL"], ["none", "q", "fsp", "na", "wo", "sil", "ol", "x", "yan"]), ss_infos("mu", 6, 5, ["none", "mu", "cel", "fcel", "pa", "mam", "fmu"], ["none", "x", "y", "f", "n", "c"]), ss_infos("ilm", 3, 2, ["none", "oilm", "dilm", "dhem"], ["none", "x", "Q"]), ss_infos("ol", 2, 1, ["none", "fo", "fa"], ["none", "x"]), ss_infos("hb", 11, 10, ["none", "tr", "tsm", "prgm", "glm", "cumm", "grnm", "a", "b", "mrb", "kprg", "tts"], ["none", "x", "y", "z", "a", "k", "c", "f", "t", "Q1", "Q2"]), ss_infos("ep", 3, 2, ["none", "cz", "ep", "fep"], ["none", "f", "Q"]), ss_infos("g", 4, 3, ["none", "py", "alm", "gr", "kho"], ["none", "x", "z", "f"]), ss_infos("chl", 7, 6, ["none", "clin", "afchl", "ames", "daph", "ochl1", "ochl4", "f3clin"], ["none", "x", "y", "f", "QAl", "Q1", "Q4"]), ss_infos("bi", 6, 5, ["none", "phl", "annm", "obi", "east", "tbi", "fbi"], ["none", "x", "y", "f", "t", "Q"]), ss_infos("dio", 7, 6, ["none", "jd", "di", "hed", "acmm", "om", "cfm", "jac"], ["none", "x", "j", "t", "c", "Qaf", "Qfm"]), ss_infos("abc", 2, 1, ["none", "abm", "anm"], ["none", "ca"])], ["sp", "opx", "fsp", "liq", "mu", "ilm", "ol", "hb", "ep", "g", "chl", "bi", "dio", "abc"], ["q", "crst", "trd", "coe", "law", "ky", "sill", "and", "ru", "sph", "sph", "ab", "H2O", "qfm", "qif", "nno", "hm", "cco", "aH2O", "aO2", "aMgO", "aFeO", "aAl2O3", "aTiO2"]), db_infos("ig", "Igneous (Holland et al., 2018)", ss_infos[ss_infos("spn", 8, 7, ["none", "nsp", "isp", "nhc", "ihc", "nmt", "imt", "pcr", "qndm"], ["none", "x", "y", "c", "t", "Q1", "Q2", "Q3"]), ss_infos("bi", 6, 5, ["none", "phl", "annm", "obi", "eas", "tbi", "fbi"], ["none", "x", "y", "f", "t", "Q"]), ss_infos("cd", 3, 2, ["none", "crd", "fcrd", "hcrd"], ["none", "x", "h"]), ss_infos("cpx", 10, 9, ["none", "di", "cfs", "cats", "crdi", "cess", "cbuf", "jd", "cen", "cfm", "kjd"], ["none", "x", "y", "o", "n", "Q", "f", "cr", "t", "k"]), ss_infos("ep", 3, 2, ["none", "cz", "ep", "fep"], ["none", "f", "Q"]), ss_infos("g", 6, 5, ["none", "py", "alm", "gr", "andr", "knom", "tig"], ["none", "x", "c", "f", "cr", "t"]), ss_infos("hb", 11, 10, ["none", "tr", "tsm", "prgm", "glm", "cumm", "grnm", "a", "b", "mrb", "kprg", "tts"], ["none", "x", "y", "z", "a", "k", "c", "f", "t", "Q1", "Q2"]), ss_infos("ilm", 3, 2, ["none", "oilm", "dilm", "dhem"], ["none", "x", "Q"]), ss_infos("liq", 12, 11, ["none", "q4L", "slL", "wo1L", "fo2L", "fa2L", "jdL", "hmL", "ekL", "tiL", "kjL", "ctL", "wat1L"], ["none", "wo", "sl", "fo", "fa", "jd", "hm", "ek", "ti", "kj", "yct", "h2o"]), ss_infos("ol", 4, 3, ["none", "mont", "fa", "fo", "cfm"], ["none", "x", "c", "Q"]), ss_infos("opx", 9, 8, ["none", "en", "fs", "fm", "odi", "mgts", "cren", "obuf", "mess", "ojd"], ["none", "x", "y", "c", "Q", "f", "t", "cr", "j"]), ss_infos("fsp", 3, 2, ["none", "ab", "an", "san"], ["none", "ca", "k"]), ss_infos("fl", 11, 10, ["none", "qfL", "slfL", "wofL", "fofL", "fafL", "jdfL", "hmfL", "ekfL", "tifL", "kjfL", "H2O"], ["none", "wo", "sl", "fo", "fa", "jd", "hm", "ek", "ti", "kj", "h2o"]), ss_infos("mu", 6, 5, ["none", "mu", "cel", "fcel", "pa", "mam", "fmu"], ["none", "x", "y", "f", "n", "c"]), ss_infos("fper", 2, 1, ["none", "per", "wu"], ["none", ""])], ["spn", "bi", "cd", "cpx", "ep", "g", "hb", "ilm", "liq", "ol", "opx", "fsp", "fl", "mu", "fper"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and", "ru", "sph", "O2", "qfm", "mw", "qif", "nno", "hm", "cco", "aH2O", "aO2", "aMgO", "aFeO", "aAl2O3", "aTiO2"]), db_infos("um", "Ultramafic (Evans & Frost., 2021)", ss_infos[ss_infos("fl", 2, 1, ["none", "H2", "H2O"], ["none", "x"]), ss_infos("ol", 2, 1, ["none", "fo", "fa"], ["none", "x"]), ss_infos("br", 2, 1, ["none", "br", "fbr"], ["none", "x"]), ss_infos("ch", 2, 1, ["none", "chum", "chuf"], ["none", "x"]), ss_infos("atg", 5, 4, ["none", "atgf", "fatg", "atgo", "aatg", "oatg"], ["none", "x", "y", "f", "t"]), ss_infos("g", 2, 1, ["none", "py", "alm"], ["none", "x"]), ss_infos("ta", 6, 5, ["none", "ta", "fta", "tao", "tats", "ota", "tap"], ["none", "x", "y", "f", "v", "Q"]), ss_infos("chl", 7, 6, ["none", "clin", "afchl", "ames", "daph", "ochl1", "ochl4", "f3clin"], ["none", "x", "y", "f", "m", "t", "QA1"]), ss_infos("spi", 3, 2, ["none", "herc", "sp", "mt"], ["none", "x", "y"]), ss_infos("opx", 5, 4, ["none", "en", "fs", "fm", "mgts", "fopx"], ["none", "x", "y", "f", "Q"]), ss_infos("po", 2, 1, ["none", "trov", "trot"], ["none", "y"]), ss_infos("anth", 5, 4, ["none", "anth", "gedf", "fant", "a", "b"], ["none", "x", "y", "z", "a"])], ["fl", "ol", "br", "ch", "atg", "g", "ta", "chl", "spi", "opx", "po", "anth"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and", "pyr", "O2", "qfm", "qif", "nno", "hm", "cco", "aH2O", "aO2", "aMgO", "aFeO", "aAl2O3", "aTiO2"])] dbs = ["mp","mb","ig","um"] id = findall(dbs .== dtb)[1] @@ -235,9 +337,10 @@ function single_point_minimization( P :: T1, test :: Int64 = 0, # if using a build-in test case, X :: VecOrMat = nothing, B :: Union{Nothing, T1, Vector{T1}} = nothing, - scp = 0, - rm_list :: Union{Nothing, Vector{Int64}} = nothing, + scp :: Int64 = 0, + rm_list :: Union{Nothing, Vector{Int64}} = nothing, W :: Union{Nothing, W_Data} = nothing, + data_in :: Union{Nothing, gmin_struct{Float64, Int64}, Vector{gmin_struct{Float64, Int64}}} = nothing, Xoxides = Vector{String}, sys_in = "mol", progressbar = true # show a progress bar or not? @@ -249,6 +352,10 @@ function single_point_minimization( P :: T1, X = [X] end + if data_in isa gmin_struct{Float64, Int64} + data_in = [data_in] + end + Out_PT = multi_point_minimization( P, T, MAGEMin_db, @@ -257,12 +364,13 @@ function single_point_minimization( P :: T1, B = B, scp = scp, rm_list = rm_list, + data_in = data_in, W = W, Xoxides = Xoxides, sys_in = sys_in, progressbar = progressbar); - return Out_PT[1] + end @@ -335,8 +443,9 @@ function multi_point_minimization(P :: T2, test :: Int64 = 0, # if using a build-in test case, X :: VecOrMat = nothing, B :: Union{Nothing, T1, Vector{T1}} = nothing, - scp = 0, + scp :: Int64 = 0, rm_list :: Union{Nothing, Vector{Int64}} = nothing, + data_in :: Union{Nothing, Vector{gmin_struct{Float64, Int64}}} = nothing, W :: Union{Nothing, W_Data} = nothing, Xoxides = Vector{String}, sys_in = "mol", @@ -377,6 +486,7 @@ function multi_point_minimization(P :: T2, progr = Progress(length(P), desc="Computing $(length(P)) points...") # progress meter end @threads :static for i in eachindex(P) + # Get thread-local buffers. As of Julia v1.9, a dynamic scheduling of # the threads is the default setting. To avoid task migration and the # resulting concurrency issues, we restrict the loop to static scheduling. @@ -391,12 +501,21 @@ function multi_point_minimization(P :: T2, gv = define_bulk_rock(gv, X[i], Xoxides, sys_in, MAGEMin_db.db); end - # compute a new point using a ccall - if isnothing(B) - out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data; scp, rm_list) + if ~isnothing(data_in) + if isnothing(B) + out = point_wise_minimization_iguess(P[i], T[i], gv, z_b, DB, splx_data; scp, rm_list, data_in = data_in[i]) + else + out = point_wise_minimization_iguess(P[i], T[i], gv, z_b, DB, splx_data; buffer_n = B[i], W = W, scp, rm_list, data_in = data_in[i]) + end else - out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data; buffer_n = B[i], W = W, scp, rm_list) + if isnothing(B) + out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data; scp, rm_list) + else + out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data; buffer_n = B[i], W = W, scp, rm_list) + end end + + Out_PT[i] = deepcopy(out) if progressbar @@ -408,6 +527,7 @@ function multi_point_minimization(P :: T2, end return Out_PT + end """ @@ -420,34 +540,28 @@ function use_predefined_bulk_rock(gv, test=0, db="ig") if db == "ig" gv.test = test gv = LibMAGEMin.get_bulk_igneous(gv) - LibMAGEMin.norm_array(gv.bulk_rock, gv.len_ox) elseif db == "igd" gv.test = test gv = LibMAGEMin.get_bulk_igneous(gv) - LibMAGEMin.norm_array(gv.bulk_rock, gv.len_ox) elseif db == "ige" gv.test = test gv = LibMAGEMin.get_bulk_igneous(gv) - LibMAGEMin.norm_array(gv.bulk_rock, gv.len_ox) elseif db == "mp" gv.test = test gv = LibMAGEMin.get_bulk_metapelite(gv) - LibMAGEMin.norm_array(gv.bulk_rock, gv.len_ox) elseif db == "mb" gv.test = test gv = LibMAGEMin.get_bulk_metabasite(gv) - LibMAGEMin.norm_array(gv.bulk_rock, gv.len_ox) elseif db == "um" gv.test = test gv = LibMAGEMin.get_bulk_ultramafic(gv) - LibMAGEMin.norm_array(gv.bulk_rock, gv.len_ox) elseif db == "alk" gv.test = test gv = LibMAGEMin.get_bulk_igneous_alk(gv) - LibMAGEMin.norm_array(gv.bulk_rock, gv.len_ox) else print("Database not implemented...\n") end + LibMAGEMin.norm_array(gv.bulk_rock, gv.len_ox) return gv end @@ -644,7 +758,18 @@ julia> finalize_MAGEMin(gv,DB) ``` """ -function point_wise_minimization(P::Float64,T::Float64, gv, z_b, DB, splx_data; buffer_n = 0.0, scp = 0, rm_list = nothing, W = nothing) +function point_wise_minimization( P ::Float64, + T ::Float64, + gv, + z_b, + DB, + splx_data; + buffer_n = 0.0, + scp = 0, + rm_list = nothing, + data_in = nothing, + W = nothing ) + gv.buffer_n = buffer_n; input_data = LibMAGEMin.io_data(); # zero (not used actually) z_b.T = T + 273.15; # in K @@ -704,32 +829,271 @@ function point_wise_minimization(P::Float64,T::Float64, gv, z_b, DB, splx_data; # Transform results to a more convenient julia struct out = deepcopy(create_gmin_struct(DB, gv, time)); + # here we compute specific heat capacity using reactions if (scp == 1) mSS_vec = deepcopy(out.mSS_vec) dT = 2.0; - W = point_wise_minimization_with_guess(mSS_vec, P, T-dT, gv, z_b, DB, splx_data) - E = point_wise_minimization_with_guess(mSS_vec, P, T+dT, gv, z_b, DB, splx_data) - hcp = -(T+273.15)*(E + W - 2.0*out.G_system)/(dT*dT); + out_W = point_wise_minimization_with_guess(mSS_vec, P, T-dT, gv, z_b, DB, splx_data) + out_E = point_wise_minimization_with_guess(mSS_vec, P, T+dT, gv, z_b, DB, splx_data) + hcp = -(T+273.15)*(out_E.G_system + out_W.G_system - 2.0*out.G_system)/(dT*dT); s_cp = hcp/out.M_sys*1e6; out.s_cp .= s_cp + end + + return out + +end + + + + + +""" + in development +""" +function point_wise_minimization_iguess( P :: Number, + T :: Number, + gv, + z_b, + DB, + splx_data; + buffer_n :: Float64 = 0.0, + scp :: Int64 = 0, + rm_list :: Union{Nothing, Vector{Int64}} = nothing, + data_in :: Union{Nothing, gmin_struct{Float64, Int64}} = nothing, + W :: Union{Nothing, W_Data} = nothing ) + + mSS_vec = deepcopy(data_in.mSS_vec) + + gv.buffer_n = buffer_n; + input_data = LibMAGEMin.io_data(); # zero (not used actually) + z_b.T = T + 273.15; # in K + + if P < 0.001 + P = 0.001 + end + + z_b.P = P + gv.numPoint = 1; # the number of the current point */ + + # 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 ) + + LibMAGEMin.reset_simplex_A(pointer_from_objref(splx_data), z_b, gv) + LibMAGEMin.reset_simplex_B_em(pointer_from_objref(splx_data), gv) + + LibMAGEMin.reset_cp(gv,z_b, DB.cp) + LibMAGEMin.reset_SS(gv,z_b, DB.SS_ref_db) + LibMAGEMin.reset_sp(gv, DB.sp) - # print("E: $E W: $W G: $(out.G_system)\n") + gv = LibMAGEMin.ComputeG0_point(gv.EM_database, z_b, gv, DB.PP_ref_db,DB.SS_ref_db); + + if ~isnothing(rm_list) + SS_ref_db = unsafe_wrap(Vector{LibMAGEMin.SS_ref},DB.SS_ref_db,gv.len_ss); + + for i in eachindex(rm_list) + flags = zeros(Int32,5); + unsafe_copyto!(SS_ref_db[rm_list[i]].ss_flags,pointer(flags), 5) + end end - # LibMAGEMin.FreeDatabases(gv, DB, z_b); + # here we can over-ride default W's + if ~isnothing(W) + if gv.EM_database == W.database # check if the database fit + else + print(" Wrong database number, please make sure the custom Ws are linked to the right database\n") + end + end + + ############################################################################ + PP_ref_db = unsafe_wrap(Vector{LibMAGEMin.PP_ref},DB.PP_ref_db,gv.len_pp); + SS_ref_db = unsafe_wrap(Vector{LibMAGEMin.SS_ref},DB.SS_ref_db,gv.len_ss); + + np = z_b.nzEl_val + nzEl_array = unsafe_wrap(Vector{Cint},z_b.nzEl_array, gv.len_ox) .+ 1 + nzEl_array = nzEl_array[1:np] + + # Declare array to be copied in splx_data + A_jll = zeros(np,np) + g0_A_jll = zeros(np) + ph_id_A_jll = zeros(Int32,np,4) + + n_pc_ss = zeros(gv.len_ss) + + # fill the arrays to be copied in splx_data + for i = 1:np + if mSS_vec[i].ph_type == "pp" + ph_id = mSS_vec[i].ph_id+1 + g0_A_jll[i] = PP_ref_db[ph_id].gbase*PP_ref_db[ph_id].factor + A_jll[i,:] = mSS_vec[i].comp_Ppc[nzEl_array] + + ph_id_A_jll[i,1] = 1 + ph_id_A_jll[i,2] = ph_id-1 + ph_id_A_jll[i,3] = 0 + ph_id_A_jll[i,4] = 0 + elseif mSS_vec[i].ph_type == "ss" + ph_id = mSS_vec[i].ph_id+1 + ph = mSS_vec[i].ph_name + + unsafe_copyto!(SS_ref_db[ph_id].gb_lvl,SS_ref_db[ph_id].gbase, SS_ref_db[ph_id].n_em) + unsafe_copyto!(SS_ref_db[ph_id].iguess,pointer(mSS_vec[i].xeos_Ppc), SS_ref_db[ph_id].n_xeos) + + SS_ref_db[ph_id] = LibMAGEMin.PC_function(gv, SS_ref_db[ph_id], z_b, ph) + + g0_A_jll[i] = SS_ref_db[ph_id].df + A_jll[i,:] = mSS_vec[i].comp_Ppc[nzEl_array] + ph_id_A_jll[i,1] = 3 + ph_id_A_jll[i,2] = ph_id-1 + ph_id_A_jll[i,3] = 0 + ph_id_A_jll[i,4] = n_pc_ss[ph_id] + n_pc_ss[ph_id] += 1 + elseif mSS_vec[i].ph_type == "ss_em" + ph_id = mSS_vec[i].ph_id+1 + em_id = mSS_vec[i].em_id+1 + ape = unsafe_wrap(Vector{Cdouble},SS_ref_db[ph_id].ape, SS_ref_db[ph_id].n_em) + gbase = unsafe_wrap(Vector{Cdouble},SS_ref_db[ph_id].gbase, SS_ref_db[ph_id].n_em) + comp_ptr= unsafe_wrap(Vector{Ptr{Cdouble}},SS_ref_db[ph_id].Comp, SS_ref_db[ph_id].n_em) + Comp = unsafe_wrap(Vector{Cdouble},comp_ptr[em_id], gv.len_ox) + factor = z_b.fbc/ape[em_id] + ph = mSS_vec[i].ph_name + + g0_A_jll[i] = gbase[em_id]*factor; + A_jll[i,:] = Comp[nzEl_array]*factor + ph_id_A_jll[i,1] = 2 + ph_id_A_jll[i,2] = ph_id-1 + ph_id_A_jll[i,3] = 0 + ph_id_A_jll[i,4] = em_id-1 + end + end + + # copy to the appropriate places + ph_id_A = unsafe_wrap(Vector{Ptr{Int32}},splx_data.ph_id_A, np) + + for i=1:np + unsafe_copyto!(ph_id_A[i],pointer(ph_id_A_jll[i,:]),4) + end + + unsafe_copyto!(splx_data.A,pointer(vec(A_jll)),np*np) + unsafe_copyto!(splx_data.A1,pointer(vec(A_jll)),np*np) + unsafe_copyto!(splx_data.g0_A,pointer(g0_A_jll),np) + + # add pseudocompounds + n_mSS = length(mSS_vec) + for i = 1:n_mSS + + if mSS_vec[i].ph_type == "ss" + ph = mSS_vec[i].ph_name + ph_id = mSS_vec[i].ph_id+1 + n_xeos = SS_ref_db[ph_id].n_xeos + n_em = SS_ref_db[ph_id].n_em + + tot_pc = unsafe_wrap(Vector{Cint},SS_ref_db[ph_id].tot_pc, 1) + id_pc = unsafe_wrap(Vector{Cint},SS_ref_db[ph_id].id_pc, 1) + info = unsafe_wrap(Vector{Cint},SS_ref_db[ph_id].info, gv.max_n_mSS) + factor_pc = unsafe_wrap(Vector{Cdouble},SS_ref_db[ph_id].factor_pc, gv.max_n_mSS) + DF_pc = unsafe_wrap(Vector{Cdouble},SS_ref_db[ph_id].DF_pc, gv.max_n_mSS) + G_pc = unsafe_wrap(Vector{Cdouble},SS_ref_db[ph_id].G_pc, gv.max_n_mSS) + + m_pc = id_pc[1]+1; + ptr_comp_pc = unsafe_wrap(Vector{Ptr{Cdouble}},SS_ref_db[ph_id].comp_pc,gv.max_n_mSS) + ptr_p_pc = unsafe_wrap(Vector{Ptr{Cdouble}},SS_ref_db[ph_id].p_pc,gv.max_n_mSS) + ptr_xeos_pc = unsafe_wrap(Vector{Ptr{Cdouble}},SS_ref_db[ph_id].xeos_pc,gv.max_n_mSS) + + unsafe_copyto!(SS_ref_db[ph_id].gb_lvl,SS_ref_db[ph_id].gbase, SS_ref_db[ph_id].n_em) + xeos = mSS_vec[i].xeos_Ppc + + # retrieve bounds + bounds_ref = zeros( n_xeos,2) + ptr_bounds_ref = unsafe_wrap(Vector{Ptr{Cdouble}}, SS_ref_db[ph_id].bounds_ref, n_xeos) + + for k=1:n_xeos + bounds_ref[k,:] = unsafe_wrap(Vector{Cdouble}, ptr_bounds_ref[k], 2) + if xeos[k] < bounds_ref[k,1] + xeos[k] = bounds_ref[k,1] + elseif xeos[k] > bounds_ref[k,2] + xeos[k] = bounds_ref[k,2] + end + end + + # get solution phase information for given compositional variables + unsafe_copyto!(SS_ref_db[ph_id].iguess,pointer(xeos), n_xeos) + SS_ref_db[ph_id] = LibMAGEMin.PC_function(gv, SS_ref_db[ph_id], z_b, ph) + + # copy solution phase composition + ss_comp = unsafe_wrap(Vector{Cdouble}, SS_ref_db[ph_id].ss_comp, gv.len_ox) + comp_pc = unsafe_wrap(Vector{Cdouble}, ptr_comp_pc[m_pc], gv.len_ox) + comp_pc .= ss_comp .* SS_ref_db[ph_id].factor; + + # copy endmember fraction + p = unsafe_wrap(Vector{Cdouble}, SS_ref_db[ph_id].p, n_em) + p_pc = unsafe_wrap(Vector{Cdouble}, ptr_p_pc[m_pc], n_em) + p_pc .= p + + # copy compositional variables + xeos_pc = unsafe_wrap(Vector{Cdouble}, ptr_xeos_pc[m_pc], n_xeos) + xeos_pc .= xeos + + info[m_pc] = 1; + factor_pc[m_pc] = SS_ref_db[ph_id].factor; + DF_pc[m_pc] = SS_ref_db[ph_id].df; + G_pc[m_pc] = SS_ref_db[ph_id].df; + + tot_pc .+= 1; + id_pc .+= 1; + end + end + gv.leveling_mode = 1 + + out = deepcopy(pwm_run(gv, z_b, DB, splx_data)) + return out end + + + """ out = point_wise_minimization(P::Number,T::Number, data::MAGEMin_Data) Performs a point-wise optimization for a given pressure `P` and temperature `T` for the data specified in the MAGEMin database `MAGEMin_Data` (where also compoition is specified) """ -point_wise_minimization(P::Number,T::Number, gv, z_b, DB, splx_data; buffer_n::Float64 = 0.0, scp::Int64 = 0, rm_list::Union{Nothing, Vector{Int64}} = nothing, W::Union{Nothing, W_Data} = nothing) = point_wise_minimization(Float64(P),Float64(T), gv, z_b, DB, splx_data; buffer_n, scp, rm_list, W) - -point_wise_minimization(P::Number,T::Number, gv::LibMAGEMin.global_variables, z_b::LibMAGEMin.bulk_infos, DB::LibMAGEMin.Database, splx_data::LibMAGEMin.simplex_datas, sys_in::String; buffer_n::Float64 = 0.0, scp::Int64 = 0, rm_list::Union{Nothing, Vector{Int64}} = nothing, W::Union{Nothing, W_Data} = nothing) = point_wise_minimization(P,T, gv, z_b, DB, splx_data; buffer_n, scp, rm_list, W) - -point_wise_minimization(P::Number,T::Number, data::MAGEMin_Data; buffer_n::Float64 = 0.0, scp::Int64 = 0, rm_list::Union{Nothing, Vector{Int64}} = nothing, W::Union{Nothing, W_Data} = nothing) = point_wise_minimization(P,T, data.gv[1], data.z_b[1], data.DB[1], data.splx_data[1]; buffer_n, scp, rm_list, W) +point_wise_minimization(P :: Number, + T :: Number, + gv, + z_b, + DB, + splx_data; + buffer_n:: Float64 = 0.0, + scp :: Int64 = 0, + rm_list :: Union{Nothing, Vector{Int64}} = nothing, + data_in :: Union{Nothing, gmin_struct{Float64, Int64}, Vector{gmin_struct{Float64, Int64}}} = nothing, + W :: Union{Nothing, W_Data} = nothing) = + point_wise_minimization(Float64(P),Float64(T), gv, z_b, DB, splx_data; buffer_n, scp, rm_list, data_in, W) + +point_wise_minimization(P :: Number, + T :: Number, + gv :: LibMAGEMin.global_variables, + z_b :: LibMAGEMin.bulk_infos, + DB :: LibMAGEMin.Database, + splx_data:: LibMAGEMin.simplex_datas, + sys_in :: String; + buffer_n:: Float64 = 0.0, + scp :: Int64 = 0, + rm_list :: Union{Nothing, Vector{Int64}} = nothing, + data_in :: Union{Nothing, gmin_struct{Float64, Int64}, Vector{gmin_struct{Float64, Int64}}} = nothing, + W :: Union{Nothing, W_Data} = nothing) = + point_wise_minimization(Float64(P),Float64(T), gv, z_b, DB, splx_data; buffer_n, scp, rm_list, data_in, W) + +point_wise_minimization(P :: Number, + T :: Number, + data :: MAGEMin_Data; + buffer_n:: Float64 = 0.0, + scp :: Int64 = 0, + rm_list :: Union{Nothing, Vector{Int64}} = nothing, + data_in :: Union{Nothing, gmin_struct{Float64, Int64}, Vector{gmin_struct{Float64, Int64}}} = nothing, + W :: Union{Nothing, W_Data} = nothing) = + point_wise_minimization(Float64(P),Float64(T), data.gv[1], data.z_b[1], data.DB[1], data.splx_data[1]; buffer_n, scp, rm_list, data_in, W) """ @@ -788,105 +1152,6 @@ end # pwm_run(gv, z_b, DB, splx_data) = pwm_run( gv, z_b, DB, splx_data) - - -""" - structure that holds the result of the pointwise minisation - -""" -struct gmin_struct{T,I} - MAGEMin_ver :: String - dataset :: String - G_system :: T # G of system - Gamma :: Vector{T} # Gamma - P_kbar :: T # Pressure in kbar - T_C :: T # Temperature in Celsius - X :: Vector{T} - M_sys :: T - - # 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} - - # Fractions: - # Solid, melt, fluid fractions - frac_M :: T - frac_S :: T - frac_F :: T - - frac_M_wt :: T - frac_S_wt :: T - frac_F_wt :: T - - # Solid, melt, fluid densities - alpha :: T - V :: T - cp :: T - s_cp :: Vector{T} - rho :: T - rho_M :: T - rho_S :: T - rho_F :: T - - # Oxygen fugacity - fO2 :: T - dQFM :: T - - # Activities - aH2O :: T - aSiO2 :: T - aTiO2 :: T - aAl2O3 :: T - aMgO :: T - aFeO :: T - - # Phase fractions and type: - n_PP :: Int64 # number of pure phases - n_SS :: Int64 # number of solid solutions - n_mSS :: Int64 # number of solid solutions - - ph_frac :: Vector{T} # phase fractions - ph_frac_wt :: Vector{T} # phase fractions - ph_frac_vol :: 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} - mSS_vec :: Vector{LibMAGEMin.mSS_data} - PP_vec :: Vector{LibMAGEMin.PP_data} - - oxides :: Vector{String} - - # Seismic velocity info - Vp :: T # P-wave velocity - Vs :: T # S-wave velocity - Vp_S :: T # P-wave velocity of solid aggregate - Vs_S :: T # S-wave velocity of solid aggregate - bulkMod :: T # Elastic bulk modulus - shearMod :: T # Elastic shear modulus - bulkModulus_M :: T # Elastic bulk modulus - bulkModulus_S :: T # Elastic bulk modulus - shearModulus_S :: T # Elastic shear modulus - - # thermodynamic properties - entropy :: T # entropy - enthalpy :: T # enthalpy - - # Numerics: - iter :: I # number of iterations required - bulk_res_norm :: T # bulk residual norm - time_ms :: T # computational time for this point - status :: I # status of calculations -end - """ out = create_gmin_struct(gv, z_b) @@ -1027,8 +1292,6 @@ function show(io::IO, g::gmin_struct) end println(io, "Oxygen fugacity : $(g.fO2)") println(io, "Delta QFM : $(g.dQFM)") - - end """ @@ -1402,12 +1665,11 @@ function point_wise_minimization_with_guess(mSS_vec, P, T, gv, z_b, DB, splx_dat gv.leveling_mode = 1 out = deepcopy(pwm_run(gv, z_b, DB, splx_data)) - return out.G_system + return out end # The following section add post-processing routines - include("TE_partitioning.jl") include("Zircon_saturation.jl") - +include("export2CSV.jl") diff --git a/julia/TE_Zr_tests.jl b/julia/TE_Zr_tests.jl new file mode 100644 index 0000000..a6e86ce --- /dev/null +++ b/julia/TE_Zr_tests.jl @@ -0,0 +1,30 @@ +using MAGEMin_C + +# Initialize database - new way +dtb = "ig" +data = Initialize_MAGEMin(dtb, verbose=false); +test = 0 #KLB1 +data = use_predefined_bulk_rock(data, test); + +# Call optimization routine for given P & T & bulk_rock +P = 8.0 +T = 1400.0 +out = point_wise_minimization(P,T, data); + +elem_TE = ["SiO2","TiO2","Al2O3","O","FeO","MnO","MgO","CaO","K2O","Na2O","H2O","Li","Be","B","Sc","V","Cr","Ni","Cu","Zn","Rb","Sr","Y","Zr","Nb","Cs","Ba","La","Ce","Pr","Nd","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","Pb","Th","U"] +bulk_TE = [50.60777553,0.953497243,13.70435413,0.19,11.28130762,0.202560796,8.496024312,9.502380068,0.700881685,2.07434927,4,29.14258603,0.434482763,29.69200003,38.23663423,257.4346716,529.333066,208.2057375,88.87615683,91.7592182,16.60777308,163.4533209,20.74016207,66.90677472,3.808354064,1.529226981,122.8449739,6.938172601,16.04827796,2.253943183,10.18276823,3.3471043,0.915941652,3.28230146,1.417695298,3.851230952,0.914966282,2.20425,0.343734976,2.136202593,0.323405135,1.841502082,0.330971265,5.452969044,1.074692058,0.290233271] + +KDs_dtb = get_OL_KDs_database(); + +C0 = adjust_chemical_system( KDs_dtb,bulk_TE,elem_TE); + +out_TE = TE_prediction(C0,KDs_dtb, "CB",out,dtb); + + +if ~isnothing(out_TE.bulk_cor_wt) #then we can recompute the equilibrium after removing the SiO2 entering zircon + sys_in = "wt" + out = single_point_minimization(P, T, data, X=out_TE.bulk_cor_wt, Xoxides=out.oxides, sys_in=sys_in) + out_TE = TE_prediction(C0,KDs_dtb, "CB",out,dtb) +end + +# add parallel test \ No newline at end of file diff --git a/julia/TE_partitioning.jl b/julia/TE_partitioning.jl index 6dce6d7..e458f7a 100644 --- a/julia/TE_partitioning.jl +++ b/julia/TE_partitioning.jl @@ -2,41 +2,10 @@ # Note that only metapelite, metabasite and igneous database can be used for trace element prediction # NR 11/04/2023 -# using MAGEMin_C - -""" - Holds the partitioning coefficient database -""" -mutable struct TE_database - conditions :: String - infos :: String - n_element :: Int64 - n_phase :: Int64 - element_name :: Vector{String} - phase_name :: Vector{String} - Kds :: Matrix{Float64} -end - -""" - This function converts mineral names from the TE database into names that can be compared with MAGEMin mineral classification -""" -function mineral_name_convertor( phase_name :: Vector{String} ) - - ph = copy(phase_name) - ph = replace.(ph,"All" => "all"); ph = replace.(ph,"Ap" => "ap"); ph = replace.(ph,"Ttn" => "ttn"); ph = replace.(ph,"Amp" => "hb"); - ph = replace.(ph,"Bt" => "bi"); ph = replace.(ph,"Crd" => "cd"); ph = replace.(ph,"Kfs" => "ksp"); ph = replace.(ph,"Pl" => "pl"); - ph = replace.(ph,"Qtz" => "q"); ph = replace.(ph,"Rt" => "ru"); ph = replace.(ph,"Grt" => "g"); ph = replace.(ph,"Ep" => "ep"); - ph = replace.(ph,"Ol" => "ol"); ph = replace.(ph,"Opx" => "opx"); ph = replace.(ph,"Mica" => "mu"); ph = replace.(ph,"Mt" => "mt"); - ph = replace.(ph,"Zrn" => "zrn"); ph = replace.(ph,"Cpx" => "cpx"); ph = replace.(ph,"Spl" => "sp"); - - return ph -end - - """ Classify the mineral output from MAGEMin to be able to be compared with partitioning coefficient database """ -function mineral_classification( out :: gmin_struct{Float64, Int64}, +function mineral_classification( out :: MAGEMin_C.gmin_struct{Float64, Int64}, dtb :: String ) ph = Array{String}(undef, out.n_SS + out.n_PP) @@ -98,119 +67,205 @@ function mineral_classification( out :: gmin_struct{Float64, Int6 end + """ - This routine stores the TE partitioning coefficients + Holds the partitioning coefficient database """ -function get_TE_database(tedb :: String) - - if tedb == "TE_OL_felsic" - conditions = "SiO2 of the melt > 63 wt" - infos = "Laurent, O. (2012). Les changements géodynamiques à la transition Archéen-Protérozoïque : étude des granitoïdes de la marge Nord du craton du Kaapvaal (Afrique du Sud). PhD, Université Blaise Pascal, Clermont-Ferrand." - element_name = ["Rb", "Ba", "Th", "U", "Nb", "Ta", "La", "Ce", "Pb", "Pr", "Sr", "Nd", "Zr", "Hf", "Sm", "Eu", "Gd", "Tb", "Dy", "Y", "Ho", "Er", "Tm", "Yb", "Lu", "V", "Sc"] - phase_name_OL = ["All", "Amp", "Ap", "Bt", "Crd", "Cpx", "FeTiOx", "Grt", "Kfs", "Mt", "Ol", "Opx", "Pl", "Qtz", "Rt", "Spl", "Ttn", "Zrn", "Ep", "and", "sill", "Mica"] - n_element = length(element_name); - n_phase = length(phase_name_OL); - Kds = [0.0632455532033676 3.46410161513776 424.264068711929 20 0.447213595499958 2.73861278752583 1549.19333848297 1224.74487139159 0.316227766016838 1005.98210719674 1 793.725393319378 0.173205080756888 12.2474487139159 447.213595499958 97.9795897113271 254.950975679639 169.558249578132 118.321595661992 95.3939201416947 84.8528137423857 63.2455532033676 43.301270189222 26.8328157299975 17.7482393492988 1 57.0087712549569; 0.0547722557505166 0.189736659610103 0.212132034355964 0.821583836257749 2 0.4 0.559016994374947 1 0.632455532033676 1.59687194226713 0.387298334620742 2.64575131106459 0.632455532033676 0.707106781186547 3.46410161513776 3.46410161513776 4.74341649025257 5.37354631505117 5.42217668469038 5.74456264653803 5.24404424085076 4.96990945591567 4.5 3.74165738677394 2.95803989154981 6.70820393249937 14.142135623731; 0.00316227766016838 0.273861278752583 1 0.948683298050514 0.0316227766016838 0.0316227766016838 13.228756555323 21.2132034355964 0.158113883008419 30.5777697028413 7.07106781186548 40.3112887414928 0.387298334620742 0.387298334620742 47.4341649025257 21.2132034355964 52.9150262212919 52.9150262212919 45.8257569495584 38.7298334620742 37.0809924354783 30 22.9128784747792 15.8113883008419 12.5499003980111 0.316227766016838 0.316227766016838; 4 6.70820393249937 0.316227766016838 0.316227766016838 6.32455532033676 1.89736659610103 1.26491106406735 0.948683298050514 0.316227766016838 0.790569415042095 0.223606797749979 0.790569415042095 0.632455532033676 0.632455532033676 0.632455532033676 0.316227766016838 0.51234753829798 0.521536192416212 0.519615242270663 0.774596669241483 0.53851648071345 0.554977477020464 0.557225268630201 0.558569601750758 0.559016994374947 3.16227766016838 9.21954445729289; 0.106066017177982 0.0223606797749979 0.193649167310371 0.632455532033676 0.0316227766016838 0.0316227766016838 0.0741619848709566 0.0848528137423857 0.0316227766016838 0.0989949493661167 0.187082869338697 0.12 0.0612372435695794 0.0774596669241483 0.173205080756888 0.0316227766016838 0.403112887414927 0.670820393249937 1.14564392373896 0.987420882906575 1.58113883008419 2.29128784747792 2.80624304008046 3.16227766016838 3.51781181986757 0.447213595499958 2.23606797749979; 0.0316227766016838 0.122474487139159 0.14142135623731 0.154919333848297 0.109544511501033 0.154919333848297 0.632455532033676 1.09544511501033 0.223606797749979 1.48323969741913 0.335410196624968 1.93649167310371 0.387298334620742 0.547722557505166 2.32379000772445 2.77488738510232 3.06186217847897 3.46410161513775 3.57770876399966 3.24037034920393 3.35410196624969 3.1859064644148 3.07408522978788 2.89827534923789 2.77488738510232 4.74341649025257 26.4575131106459; 0.0223606797749979 0.0223606797749979 0.387298334620742 0.387298334620742 28.2842712474619 44.7213595499958 3.35410196624968 2.82842712474619 0.447213595499958 2.57875939164553 0.193649167310371 2.50998007960223 1.58113883008419 1.58113883008419 2.36643191323985 0.632455532033676 2.36643191323985 1.89736659610103 1.58113883008419 0.447213595499958 1.58113883008419 1.58113883008419 1.42302494707577 1.26491106406735 1.26491106406735 31.6227766016838 4.47213595499958; 0.00790569415042094 0.0173205080756888 0.0774596669241483 0.158113883008419 0.0316227766016838 0.0707106781186547 0.1 0.237170824512628 0.0316227766016838 0.460977222864644 0.0223606797749979 0.866025403784439 0.447213595499958 0.707106781186547 1.93649167310371 1.58113883008419 6.12372435695795 11.180339887499 22.3606797749979 25.4950975679639 32.4037034920393 41.2310562561766 44.7213595499958 44.7213595499958 37.4165738677394 4.47213595499958 22.3606797749979; 0.948683298050514 9.48683298050514 0.0144913767461894 0.0158113883008419 0.0447213595499958 0.00707106781186547 0.14142135623731 0.0866025403784439 0.948683298050514 0.0547722557505166 3.87298334620742 0.0316227766016838 0.0193649167310371 0.0158113883008419 0.0187082869338697 3.87298334620742 0.0187082869338697 0.0167332005306815 0.02 0.0316227766016838 0.02 0.0223606797749979 0.0234520787991171 0.0234520787991171 0.0244948974278318 0.223606797749979 0.0316227766016838; 0.0223606797749979 0.0223606797749979 0.221359436211787 0.447213595499958 1.67332005306815 2.12132034355964 1.02469507659596 1 0.547722557505166 1.06066017177982 0.0264575131106459 1.09544511501033 0.670820393249937 0.707106781186547 1.09544511501033 0.948683298050514 1.18321595661992 1.18321595661992 1.10679718105893 0.948683298050514 1.02469507659596 0.866025403784439 0.866025403784439 0.707106781186547 0.547722557505166 36.0555127546399 8.66025403784439; 0.0223606797749979 0.03 0.02 0.0154919333848297 0.031224989991992 0.0335410196624968 0.0670820393249937 0.0774596669241483 0.158113883008419 0.0916515138991168 0.13228756555323 0.106066017177982 0.025298221281347 0.0284604989415154 0.121963109176505 0.14142135623731 0.13228756555323 0.154919333848297 0.175499287747842 0.126491106406735 0.193649167310371 0.232379000772445 0.277488738510232 0.316227766016838 0.387298334620742 0.244948974278318 0.447213595499958; 0.0158113883008419 0.05 0.158113883008419 0.158113883008419 0.25298221281347 0.316227766016838 0.158113883008419 0.239791576165636 0.0894427190999916 0.289827534923789 0.0447213595499958 0.353553390593274 0.0547722557505166 0.0948683298050513 0.424264068711928 0.316227766016838 0.58309518948453 0.724568837309472 0.848528137423857 0.774596669241483 0.924662100445346 1 1.16189500386223 1.2747548783982 1.54919333848297 3.16227766016838 22.3606797749979; 0.126491106406735 0.632455532033676 0.0547722557505166 0.1 0.1 0.0670820393249937 0.387298334620742 0.284604989415154 0.670820393249937 0.221359436211786 4.47213595499958 0.158113883008419 0.0894427190999916 0.0670820393249937 0.126491106406735 2.73861278752583 0.116189500386222 0.102469507659596 0.0948683298050513 0.1 0.0774596669241483 0.0692820323027551 0.0574456264653803 0.0524404424085076 0.0447213595499958 0.158113883008419 0.0316227766016838; 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838; 0.0118321595661992 0.01 0.335410196624968 0.335410196624968 61.2372435695795 27.3861278752583 0.0353553390593274 0.0458257569495584 0.0316227766016838 0.0494974746830583 0.0707106781186547 0.0565685424949238 3.46410161513776 5.29150262212918 0.0565685424949238 0.00111803398874989 0.0354964786985977 0.0278208554864871 0.0256904651573302 0.0707106781186547 0.0241867732448956 0.0217944947177034 0.0185202591774521 0.0173205080756888 0.0158113883008419 47.4341649025257 0.316227766016838; 0.00547722557505166 0.00223606797749979 0.00790569415042094 0.013228756555323 0.0707106781186547 0.0707106781186547 0.00244948974278318 0.00387298334620742 0.000316227766016838 0.006 0.00346410161513776 0.00812403840463596 0.0790569415042095 0.1 0.0116189500386222 0.0106066017177982 0.0144913767461894 0.0149666295470958 0.015 0.0158113883008419 0.0141421356237309 0.0132664991614216 0.0130384048104053 0.012369316876853 0.0116189500386222 8.36660026534076 0.316227766016838; 0.4 1.36930639376292 0.223606797749979 0.2 3.46410161513776 8.66025403784439 6 7.41619848709567 0.223606797749979 8.83176086632785 2.23606797749979 10.2469507659596 1.34164078649987 1.73205080756888 10.9544511501033 8.06225774829855 10.2469507659596 9.2951600308978 7.88035532193822 6 6.557438524302 5.17010638188423 3.73496987939662 2.56904651573303 1.67332005306815 5.47722557505166 2.23606797749979; 0.632455532033676 0.632455532033676 18.02775637732 31.6227766016838 22.3606797749979 22.3606797749979 1.26491106406735 3.46410161513776 0.1 2 4.47213595499958 2.82842712474619 948.683298050515 948.683298050515 7.74596669241484 2.23606797749979 20.4939015319192 28.4604989415154 44.1588043316392 67.0820393249938 73.4846922834954 110 139.64240043769 173.205080756888 223.606797749979 0.1 63.2455532033676; 0.0045 0.408 156 1.29 0.226 0.226 2.05 2.44 0.5 2.85475042692001 2 3.34 0.1 10 4.22 3.78 4.67 4.58421203698084 4.5 4.3 4.12431812546026 3.78 3.34496636754392 2.96 2.61933874283863 0.1 0.0001; 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5; 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5; 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5]; - phase_name = mineral_name_convertor(phase_name_OL); - TE_dtb = TE_database(conditions,infos,n_element,n_phase,element_name,phase_name,Kds) - - elseif tedb == "TE_OL_intermediate" - conditions = "SiO2 of the melt < 63 wt && SiO2 of the melt > 52 wt" - infos = "Laurent, O. (2012). Les changements géodynamiques à la transition Archéen-Protérozoïque : étude des granitoïdes de la marge Nord du craton du Kaapvaal (Afrique du Sud). PhD, Université Blaise Pascal, Clermont-Ferrand." - element_name = ["Rb", "Ba", "Th", "U", "Nb", "Ta", "La", "Ce", "Pb", "Pr", "Sr", "Nd", "Zr", "Hf", "Sm", "Eu", "Gd", "Tb", "Dy", "Y", "Ho", "Er", "Tm", "Yb", "Lu", "V", "Sc"] - phase_name_OL = ["All"; "Amp"; "Ap"; "Bt"; "Crd"; "Cpx"; "FeTiOx"; "Grt"; "Kfs"; "Mt"; "Ol"; "Opx"; "Pl"; "Qtz"; "Rt"; "Spl"; "Ttn"; "Zrn"] - n_element = length(element_name); - n_phase = length(phase_name_OL); - Kds = [0.0632455532033676 3.46410161513776 324.037034920393 13.228756555323 0.836660026534076 0.836660026534076 1414.2135623731 1174.73401244707 1 948.683298050515 1 612.372435695795 0.173205080756888 11.180339887499 300 34.6410161513775 116.189500386222 58.309518948453 31.6227766016838 31.6227766016838 18.9736659610103 13.4164078649987 10.2469507659596 7.74596669241484 6.70820393249937 1 57.0087712549569; 0.122474487139159 0.223606797749979 0.14142135623731 0.1 0.559016994374947 0.273861278752583 0.335410196624968 0.412310562561766 0.223606797749979 0.574456264653803 0.273861278752583 0.866025403784439 0.4 0.547722557505166 1.36930639376292 1 1.80277563773199 2.21359436211787 2.29128784747792 2.12132034355964 2.08566536146142 1.8165902124585 1.51657508881031 1.40712472794703 1.30384048104053 5.47722557505166 10; 0.00316227766016838 0.0790569415042095 0.790569415042095 0.387298334620742 0.00353553390593274 0.00707106781186547 3.87298334620742 7.07106781186548 0.316227766016838 11.0679718105893 1.93649167310371 15 0.193649167310371 0.0707106781186547 17.3205080756888 11.180339887499 18.3711730708738 17.8885438199983 15.6524758424985 12.2474487139159 13.4164078649987 11.180339887499 8.94427190999916 6.70820393249937 4.47213595499958 0.223606797749979 0.316227766016838; 3.16227766016838 5.47722557505166 0.0474341649025257 0.0707106781186547 0.474341649025257 0.316227766016838 0.0591607978309961 0.06 0.0316227766016838 0.0648074069840786 0.234520787991171 0.0670820393249937 0.229128784747792 0.316227766016838 0.0707106781186547 0.13228756555323 0.0741619848709566 0.0774596669241483 0.0793725393319377 0.111803398874989 0.0774596669241483 0.0747663025700749 0.0734846922834953 0.0726636084983398 0.0707106781186547 3.16227766016838 9.21954445729289; 0.106066017177982 0.0223606797749979 0.193649167310371 0.632455532033676 0.0316227766016838 0.0316227766016838 0.0741619848709566 0.0848528137423857 0.0316227766016838 0.0989949493661167 0.187082869338697 0.12 0.0612372435695794 0.0774596669241483 0.173205080756888 0.0316227766016838 0.403112887414927 0.670820393249937 1.14564392373896 0.987420882906575 1.58113883008419 2.29128784747792 2.80624304008046 3.16227766016838 3.51781181986757 0.447213595499958 2.23606797749979; 0.0316227766016838 0.0547722557505166 0.0707106781186547 0.0447213595499958 0.111803398874989 0.0707106781186547 0.13228756555323 0.212132034355964 0.316227766016838 0.290688837074973 0.287228132326901 0.401870625948202 0.234520787991171 0.31224989991992 0.5 0.692820323027551 0.725603197346869 0.848528137423857 0.953939201416946 1.09544511501033 1.06348483769163 1.14105214604767 1.18215904175369 1.1861703081767 1.18321595661992 1.93649167310371 7.74596669241484; 0.00316227766016838 0.00447213595499958 0.0158113883008419 0.0237170824512628 4.47213595499958 6.70820393249937 0.0223606797749979 0.0306186217847897 0.0158113883008419 0.0367423461417476 0.00632455532033676 0.0412310562561766 1.22474487139159 1.5 0.0547722557505166 0.0447213595499958 0.0790569415042095 0.09 0.0989949493661167 0.0707106781186547 0.101980390271856 0.111803398874989 0.117260393995586 0.117473401244707 0.109544511501033 11.180339887499 1.73205080756888; 0.00836660026534076 0.0154919333848297 0.0474341649025257 0.0935414346693485 0.0264575131106459 0.05 0.0264575131106459 0.0519615242270663 0.0316227766016838 0.102469507659596 0.0173205080756888 0.205548047910945 0.547722557505166 0.387298334620742 0.433012701892219 0.5 1.58113883008419 3.16227766016838 5.45435605731786 6.12372435695795 8.2915619758885 9.89949493661167 10.0623058987491 9.35414346693486 7.74596669241484 2.64575131106459 3.16227766016838; 0.948683298050514 5.47722557505166 0.0173205080756888 0.0273861278752583 0.0158113883008419 0.00316227766016838 0.111803398874989 0.0670820393249937 0.223606797749979 0.0447213595499958 2.23606797749979 0.03 0.0223606797749979 0.0223606797749979 0.0212132034355964 3.35410196624969 0.0167332005306815 0.015 0.0144913767461894 0.0254950975679639 0.0162018517460197 0.0189736659610103 0.02 0.0232379000772445 0.0273861278752583 0.111803398874989 0.02; 0.0387298334620742 0.158113883008419 0.0935414346693485 0.13228756555323 0.136930639376291 0.122474487139159 0.14142135623731 0.189736659610103 0.316227766016838 0.234520787991171 0.0223606797749979 0.267394839142419 0.220794021658196 0.158113883008419 0.3 0.254950975679639 0.324037034920393 0.31224989991992 0.301662062579967 0.335410196624968 0.268328157299975 0.245967477524977 0.234520787991171 0.223606797749979 0.212132034355964 22.3606797749979 3.16227766016838; 0.013228756555323 0.0122474487139159 0.0154919333848297 0.0223606797749979 0.0193649167310371 0.0273861278752583 0.00273861278752583 0.00489897948556636 0.0187082869338697 0.00790569415042094 0.01 0.0116189500386222 0.0387298334620742 0.0223606797749979 0.0152970585407783 0.0187616630392937 0.0228035085019828 0.0270554985169374 0.0324037034920393 0.0346410161513775 0.0369864840178138 0.04 0.0439545219516718 0.0474341649025257 0.0529150262212918 0.16431676725155 0.433012701892219; 0.0244948974278318 0.0316227766016838 0.1 0.0316227766016838 0.05 0.0670820393249937 0.00894427190999916 0.0204939015319192 0.111803398874989 0.0363318042491699 0.0193649167310371 0.0574456264653803 0.0346410161513775 0.0547722557505166 0.0812403840463596 0.114017542509914 0.157480157480236 0.207846096908265 0.256124969497314 0.282842712474619 0.311126983722081 0.351425667816112 0.4 0.458257569495584 0.5 1.5 5.47722557505166; 0.0447213595499958 0.387298334620742 0.0316227766016838 0.0547722557505166 0.0612372435695794 0.0447213595499958 0.2 0.154919333848297 0.5 0.122474487139159 2.44948974278318 0.106066017177982 0.0223606797749979 0.0324037034920393 0.0774596669241483 1.1180339887499 0.0447213595499958 0.0367423461417476 0.03 0.0292403830344269 0.0256904651573302 0.0212132034355964 0.0182345825288105 0.0169705627484771 0.0154919333848297 0.158113883008419 0.0316227766016838; 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838; 0.01 0.00866025403784438 0.162018517460196 0.2 33.5410196624969 22.3606797749979 0.0111803398874989 0.0141421356237309 0.0273861278752583 0.0156524758424985 0.0316227766016838 0.0172481883106603 2.23606797749979 3.74165738677394 0.0184932420089069 0.000632455532033676 0.0194935886896179 0.0204939015319192 0.022248595461287 0.0273861278752583 0.024 0.0264575131106459 0.0291204395571221 0.0314642654451045 0.033466401061363 3.16227766016838 0.273861278752583; 0.00547722557505166 0.00223606797749979 0.00790569415042094 0.013228756555323 0.0707106781186547 0.0707106781186547 0.00244948974278318 0.00387298334620742 0.000316227766016838 0.006 0.00346410161513776 0.00812403840463596 0.0790569415042095 0.1 0.0116189500386222 0.0106066017177982 0.0144913767461894 0.0149666295470958 0.015 0.0158113883008419 0.0141421356237309 0.0132664991614216 0.0130384048104053 0.012369316876853 0.0116189500386222 8.36660026534076 0.316227766016838; 0.4 1.36930639376292 0.223606797749979 0.2 3.46410161513776 8.66025403784439 6 7.41619848709567 0.223606797749979 8.83176086632785 2.23606797749979 10.2469507659596 1.34164078649987 1.73205080756888 10.9544511501033 8.06225774829855 10.2469507659596 9.2951600308978 7.88035532193822 6 6.557438524302 5.17010638188423 3.73496987939662 2.56904651573303 1.67332005306815 5.47722557505166 2.23606797749979; 0.632455532033676 0.632455532033676 18.02775637732 31.6227766016838 22.3606797749979 22.3606797749979 1.26491106406735 3.46410161513776 0.1 2 4.47213595499958 2.82842712474619 948.683298050515 948.683298050515 7.74596669241484 2.23606797749979 20.4939015319192 28.4604989415154 44.1588043316392 67.0820393249938 73.4846922834954 110 139.64240043769 173.205080756888 223.606797749979 0.1 63.2455532033676] - phase_name = mineral_name_convertor(phase_name_OL); - TE_dtb = TE_database(conditions,infos,n_element,n_phase,element_name,phase_name,Kds) - - elseif tedb == "TE_OL_mafic" - conditions = "SiO2 of the melt < 52 wt" - infos = "Laurent, O. (2012). Les changements géodynamiques à la transition Archéen-Protérozoïque : étude des granitoïdes de la marge Nord du craton du Kaapvaal (Afrique du Sud). PhD, Université Blaise Pascal, Clermont-Ferrand." - element_name = ["Rb", "Ba", "Th", "U", "Nb", "Ta", "La", "Ce", "Pb", "Pr", "Sr", "Nd", "Zr", "Hf", "Sm", "Eu", "Gd", "Tb", "Dy", "Y", "Ho", "Er", "Tm", "Yb", "Lu", "V", "Sc"] - phase_name_OL = ["All"; "Amp"; "Ap"; "Bt"; "Crd"; "Cpx"; "FeTiOx"; "Grt"; "Kfs"; "Mt"; "Ol"; "Opx"; "Pl"; "Qtz"; "Rt"; "Spl"; "Ttn"; "Zrn"; "Ep"] - n_element = length(element_name); - n_phase = length(phase_name_OL); - Kds = [0.0632455532033676 3.46410161513776 324.037034920393 13.228756555323 0.836660026534076 0.836660026534076 1414.2135623731 1174.73401244707 1 948.683298050515 1 612.372435695795 0.173205080756888 11.180339887499 300 34.6410161513775 116.189500386222 58.309518948453 31.6227766016838 31.6227766016838 18.9736659610103 13.4164078649987 10.2469507659596 7.74596669241484 6.70820393249937 1 57.0087712549569; 0.223606797749979 0.31224989991992 0.0141421356237309 0.0193649167310371 0.158113883008419 0.220794021658196 0.0724568837309472 0.144913767461894 0.0632455532033676 0.234520787991171 0.3 0.346410161513775 0.335410196624968 0.391152144312159 0.447213595499958 0.524404424085076 0.589915248150105 0.634822809924155 0.66932802122726 0.612372435695794 0.681175454637056 0.648074069840786 0.603738353924943 0.53619026473818 0.458257569495584 4.89897948556636 3.87298334620742; 0.00316227766016838 0.0790569415042095 0.790569415042095 0.387298334620742 0.00353553390593274 0.00707106781186547 3.87298334620742 7.07106781186548 0.316227766016838 11.0679718105893 1.93649167310371 15 0.193649167310371 0.0707106781186547 17.3205080756888 11.180339887499 18.3711730708738 17.8885438199983 15.6524758424985 12.2474487139159 13.4164078649987 11.180339887499 8.94427190999916 6.70820393249937 4.47213595499958 0.223606797749979 0.316227766016838; 2 4.47213595499958 0.0474341649025257 0.0707106781186547 0.273861278752583 0.187082869338697 0.0144913767461894 0.0189736659610103 0.0316227766016838 0.0244948974278318 0.234520787991171 0.031224989991992 0.0547722557505166 0.0547722557505166 0.0387298334620742 0.0418330013267038 0.0447213595499958 0.0519615242270663 0.0591607978309961 0.066332495807108 0.0692820323027551 0.0747663025700749 0.0793725393319377 0.0819756061276768 0.0866025403784439 3.16227766016838 9.21954445729289; 0.106066017177982 0.0223606797749979 0.193649167310371 0.632455532033676 0.0316227766016838 0.0316227766016838 0.0741619848709566 0.0848528137423857 0.0316227766016838 0.0989949493661167 0.187082869338697 0.12 0.0612372435695794 0.0774596669241483 0.173205080756888 0.0316227766016838 0.403112887414927 0.670820393249937 1.14564392373896 0.987420882906575 1.58113883008419 2.29128784747792 2.80624304008046 3.16227766016838 3.51781181986757 0.447213595499958 2.23606797749979; 0.00948683298050514 0.005 0.00935414346693485 0.00707106781186547 0.00894427190999916 0.0180277563773199 0.0474341649025257 0.0866025403784439 0.0223606797749979 0.14456832294801 0.122474487139159 0.216333076527839 0.0707106781186547 0.187082869338697 0.301662062579967 0.387298334620742 0.460977222864644 0.522494019104525 0.551361950083609 0.6 0.585662018573853 0.608276253029822 0.608276253029822 0.585662018573853 0.559910707166777 1.58113883008419 4.74341649025257; 0.00316227766016838 0.00447213595499958 0.0158113883008419 0.0237170824512628 4.47213595499958 6.70820393249937 0.0223606797749979 0.0306186217847897 0.0158113883008419 0.0367423461417476 0.00632455532033676 0.0412310562561766 1.22474487139159 1.5 0.0547722557505166 0.0447213595499958 0.0790569415042095 0.09 0.0989949493661167 0.0707106781186547 0.101980390271856 0.111803398874989 0.117260393995586 0.117473401244707 0.109544511501033 11.180339887499 1.73205080756888; 0.000866025403784438 0.000387298334620741 0.002 0.00612372435695794 0.013228756555323 0.0180277563773199 0.00894427190999916 0.0219089023002067 0.0122474487139159 0.0565685424949238 0.013228756555323 0.122474487139159 0.3 0.252487623459052 0.264575131106459 0.424264068711928 0.821583836257749 1.3228756555323 1.93649167310371 2.32379000772445 2.54950975679639 3.1224989991992 3.51781181986757 4.02336923485777 4.47213595499958 2.44948974278318 3.74165738677394; 0.948683298050514 5.47722557505166 0.0173205080756888 0.0273861278752583 0.0158113883008419 0.00316227766016838 0.111803398874989 0.0670820393249937 0.223606797749979 0.0447213595499958 2.23606797749979 0.03 0.0223606797749979 0.0223606797749979 0.0212132034355964 3.35410196624969 0.0167332005306815 0.015 0.0144913767461894 0.0254950975679639 0.0162018517460197 0.0189736659610103 0.02 0.0232379000772445 0.0273861278752583 0.111803398874989 0.02; 0.0158113883008419 0.0158113883008419 0.0387298334620742 0.0387298334620742 0.1 0.1 0.0126491106406735 0.0273861278752583 0.0591607978309961 0.0474341649025257 0.0122474487139159 0.062449979983984 0.150831031289984 0.2 0.0724568837309471 0.0360555127546399 0.0774596669241483 0.0702851335632223 0.0666333249958307 0.0591607978309961 0.0620483682299543 0.0565685424949238 0.0547722557505166 0.05 0.0447213595499958 6.32455532033676 1.4142135623731; 0.00223606797749979 0.00223606797749979 0.000707106781186547 0.000707106781186547 0.00273861278752583 0.005 0.000223606797749979 0.000447213595499958 0.000193649167310371 0.000935414346693485 0.000109544511501033 0.00158113883008419 0.00316227766016838 0.00547722557505166 0.00250998007960223 0.00316227766016838 0.00424264068711928 0.0062928530890209 0.00908295106229247 0.00774596669241483 0.013228756555323 0.0193649167310371 0.0241867732448956 0.0273861278752583 0.0282842712474619 0.0774596669241483 0.287228132326901; 0.00547722557505166 0.00173205080756888 0.000316227766016838 0.0005 0.00346410161513776 0.00632455532033676 0.000836660026534075 0.00264575131106459 0.02 0.00547722557505166 0.0141421356237309 0.0109544511501033 0.0291547594742265 0.0547722557505166 0.0180277563773199 0.0223606797749979 0.0324037034920393 0.0412310562561766 0.05 0.0670820393249937 0.0619677335393187 0.0734846922834953 0.0916515138991168 0.107238052947636 0.120415945787923 1 1.22474487139159; 0.0547722557505166 0.335410196624968 0.1 0.0632455532033676 0.0447213595499958 0.0447213595499958 0.111803398874989 0.0894427190999916 0.591607978309962 0.0714142842854285 1.6583123951777 0.0591607978309961 0.01 0.0154919333848297 0.0489897948556635 0.707106781186547 0.0387298334620742 0.0338230690505755 0.03 0.0357071421427143 0.0267394839142419 0.0240831891575846 0.0222261107708929 0.0203469899493758 0.0178885438199983 0.0670820393249937 0.0316227766016838; 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838; 0.01 0.00866025403784438 0.162018517460196 0.2 33.5410196624969 22.3606797749979 0.0111803398874989 0.0141421356237309 0.0273861278752583 0.0156524758424985 0.0316227766016838 0.0172481883106603 2.23606797749979 3.74165738677394 0.0184932420089069 0.000632455532033676 0.0194935886896179 0.0204939015319192 0.022248595461287 0.0273861278752583 0.024 0.0264575131106459 0.0291204395571221 0.0314642654451045 0.033466401061363 3.16227766016838 0.273861278752583; 0.00547722557505166 0.00223606797749979 0.00790569415042094 0.013228756555323 0.0707106781186547 0.0707106781186547 0.00244948974278318 0.00387298334620742 0.000316227766016838 0.006 0.00346410161513776 0.00812403840463596 0.0790569415042095 0.1 0.0116189500386222 0.0106066017177982 0.0144913767461894 0.0149666295470958 0.015 0.0158113883008419 0.0141421356237309 0.0132664991614216 0.0130384048104053 0.012369316876853 0.0116189500386222 8.36660026534076 0.316227766016838; 0.4 1.36930639376292 0.223606797749979 0.2 3.46410161513776 8.66025403784439 6 7.41619848709567 0.223606797749979 8.83176086632785 2.23606797749979 10.2469507659596 1.34164078649987 1.73205080756888 10.9544511501033 8.06225774829855 10.2469507659596 9.2951600308978 7.88035532193822 6 6.557438524302 5.17010638188423 3.73496987939662 2.56904651573303 1.67332005306815 5.47722557505166 2.23606797749979; 0.632455532033676 0.632455532033676 18.02775637732 31.6227766016838 22.3606797749979 22.3606797749979 1.26491106406735 0.1 0.316227766016838 2 4.47213595499958 2.82842712474619 948.683298050515 948.683298050515 7.74596669241484 2.23606797749979 20.4939015319192 28.4604989415154 44.1588043316392 67.0820393249938 73.4846922834954 110 139.64240043769 173.205080756888 223.606797749979 0.1 63.2455532033676; 0.0045 0.408 156 1.29 0.226 0.226 2.05 2.44 0.5 2.89 2 3.78 0.1 10 4.22 3.78 4.67 4.59 4.5 4.3 4.14 3.78 3.37 2.96 2.55 0.1 0.0001] - phase_name = mineral_name_convertor(phase_name_OL); - TE_dtb = TE_database(conditions,infos,n_element,n_phase,element_name,phase_name,Kds) +struct KDs_database + infos :: String + element_name :: Vector{String} + conditions :: Tuple{String, Vector{Vector{Float64}}} + phase_name :: Tuple{Vector{String}, Vector{String}, Vector{String}} + KDs :: Tuple{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} +end + + +function get_OL_KDs_database() + infos = "Laurent, O. (2012). Les changements géodynamiques à la transition Archéen-Protérozoïque : étude des granitoïdes de la marge Nord du craton du Kaapvaal (Afrique du Sud). PhD, Université Blaise Pascal, Clermont-Ferrand." + + element_name = ["Rb", "Ba", "Th", "U", "Nb", "Ta", "La", "Ce", "Pb", "Pr", "Sr", "Nd", "Zr", "Hf", "Sm", "Eu", "Gd", "Tb", "Dy", "Y", "Ho", "Er", "Tm", "Yb", "Lu", "V", "Sc"] + conditions = ("SiO2",[[0.0,52.0],[52.0,63.0],[63.0,100.0]]) + + ph_1 = ["all", "hb", "ap", "bi", "cd", "cpx", "FeTiOx", "g", "ksp", "mt", "ol", "opx", "pl", "q", "ru", "sp", "ttn", "zrn", "ep", "and", "sill", "mu"] + ph_2 = ["all"; "hb"; "ap"; "bi"; "cd"; "cpx"; "FeTiOx"; "g"; "ksp"; "mt"; "ol"; "opx"; "pl"; "q"; "ru"; "sp"; "ttn"; "zrn"] + ph_3 = ["all"; "hb"; "ap"; "bi"; "cd"; "cpx"; "FeTiOx"; "g"; "ksp"; "mt"; "ol"; "opx"; "pl"; "q"; "ru"; "sp"; "ttn"; "zrn"; "ep"] + + KDs_1 = [0.0632455532033676 3.46410161513776 424.264068711929 20 0.447213595499958 2.73861278752583 1549.19333848297 1224.74487139159 0.316227766016838 1005.98210719674 1 793.725393319378 0.173205080756888 12.2474487139159 447.213595499958 97.9795897113271 254.950975679639 169.558249578132 118.321595661992 95.3939201416947 84.8528137423857 63.2455532033676 43.301270189222 26.8328157299975 17.7482393492988 1 57.0087712549569; 0.0547722557505166 0.189736659610103 0.212132034355964 0.821583836257749 2 0.4 0.559016994374947 1 0.632455532033676 1.59687194226713 0.387298334620742 2.64575131106459 0.632455532033676 0.707106781186547 3.46410161513776 3.46410161513776 4.74341649025257 5.37354631505117 5.42217668469038 5.74456264653803 5.24404424085076 4.96990945591567 4.5 3.74165738677394 2.95803989154981 6.70820393249937 14.142135623731; 0.00316227766016838 0.273861278752583 1 0.948683298050514 0.0316227766016838 0.0316227766016838 13.228756555323 21.2132034355964 0.158113883008419 30.5777697028413 7.07106781186548 40.3112887414928 0.387298334620742 0.387298334620742 47.4341649025257 21.2132034355964 52.9150262212919 52.9150262212919 45.8257569495584 38.7298334620742 37.0809924354783 30 22.9128784747792 15.8113883008419 12.5499003980111 0.316227766016838 0.316227766016838; 4 6.70820393249937 0.316227766016838 0.316227766016838 6.32455532033676 1.89736659610103 1.26491106406735 0.948683298050514 0.316227766016838 0.790569415042095 0.223606797749979 0.790569415042095 0.632455532033676 0.632455532033676 0.632455532033676 0.316227766016838 0.51234753829798 0.521536192416212 0.519615242270663 0.774596669241483 0.53851648071345 0.554977477020464 0.557225268630201 0.558569601750758 0.559016994374947 3.16227766016838 9.21954445729289; 0.106066017177982 0.0223606797749979 0.193649167310371 0.632455532033676 0.0316227766016838 0.0316227766016838 0.0741619848709566 0.0848528137423857 0.0316227766016838 0.0989949493661167 0.187082869338697 0.12 0.0612372435695794 0.0774596669241483 0.173205080756888 0.0316227766016838 0.403112887414927 0.670820393249937 1.14564392373896 0.987420882906575 1.58113883008419 2.29128784747792 2.80624304008046 3.16227766016838 3.51781181986757 0.447213595499958 2.23606797749979; 0.0316227766016838 0.122474487139159 0.14142135623731 0.154919333848297 0.109544511501033 0.154919333848297 0.632455532033676 1.09544511501033 0.223606797749979 1.48323969741913 0.335410196624968 1.93649167310371 0.387298334620742 0.547722557505166 2.32379000772445 2.77488738510232 3.06186217847897 3.46410161513775 3.57770876399966 3.24037034920393 3.35410196624969 3.1859064644148 3.07408522978788 2.89827534923789 2.77488738510232 4.74341649025257 26.4575131106459; 0.0223606797749979 0.0223606797749979 0.387298334620742 0.387298334620742 28.2842712474619 44.7213595499958 3.35410196624968 2.82842712474619 0.447213595499958 2.57875939164553 0.193649167310371 2.50998007960223 1.58113883008419 1.58113883008419 2.36643191323985 0.632455532033676 2.36643191323985 1.89736659610103 1.58113883008419 0.447213595499958 1.58113883008419 1.58113883008419 1.42302494707577 1.26491106406735 1.26491106406735 31.6227766016838 4.47213595499958; 0.00790569415042094 0.0173205080756888 0.0774596669241483 0.158113883008419 0.0316227766016838 0.0707106781186547 0.1 0.237170824512628 0.0316227766016838 0.460977222864644 0.0223606797749979 0.866025403784439 0.447213595499958 0.707106781186547 1.93649167310371 1.58113883008419 6.12372435695795 11.180339887499 22.3606797749979 25.4950975679639 32.4037034920393 41.2310562561766 44.7213595499958 44.7213595499958 37.4165738677394 4.47213595499958 22.3606797749979; 0.948683298050514 9.48683298050514 0.0144913767461894 0.0158113883008419 0.0447213595499958 0.00707106781186547 0.14142135623731 0.0866025403784439 0.948683298050514 0.0547722557505166 3.87298334620742 0.0316227766016838 0.0193649167310371 0.0158113883008419 0.0187082869338697 3.87298334620742 0.0187082869338697 0.0167332005306815 0.02 0.0316227766016838 0.02 0.0223606797749979 0.0234520787991171 0.0234520787991171 0.0244948974278318 0.223606797749979 0.0316227766016838; 0.0223606797749979 0.0223606797749979 0.221359436211787 0.447213595499958 1.67332005306815 2.12132034355964 1.02469507659596 1 0.547722557505166 1.06066017177982 0.0264575131106459 1.09544511501033 0.670820393249937 0.707106781186547 1.09544511501033 0.948683298050514 1.18321595661992 1.18321595661992 1.10679718105893 0.948683298050514 1.02469507659596 0.866025403784439 0.866025403784439 0.707106781186547 0.547722557505166 36.0555127546399 8.66025403784439; 0.0223606797749979 0.03 0.02 0.0154919333848297 0.031224989991992 0.0335410196624968 0.0670820393249937 0.0774596669241483 0.158113883008419 0.0916515138991168 0.13228756555323 0.106066017177982 0.025298221281347 0.0284604989415154 0.121963109176505 0.14142135623731 0.13228756555323 0.154919333848297 0.175499287747842 0.126491106406735 0.193649167310371 0.232379000772445 0.277488738510232 0.316227766016838 0.387298334620742 0.244948974278318 0.447213595499958; 0.0158113883008419 0.05 0.158113883008419 0.158113883008419 0.25298221281347 0.316227766016838 0.158113883008419 0.239791576165636 0.0894427190999916 0.289827534923789 0.0447213595499958 0.353553390593274 0.0547722557505166 0.0948683298050513 0.424264068711928 0.316227766016838 0.58309518948453 0.724568837309472 0.848528137423857 0.774596669241483 0.924662100445346 1 1.16189500386223 1.2747548783982 1.54919333848297 3.16227766016838 22.3606797749979; 0.126491106406735 0.632455532033676 0.0547722557505166 0.1 0.1 0.0670820393249937 0.387298334620742 0.284604989415154 0.670820393249937 0.221359436211786 4.47213595499958 0.158113883008419 0.0894427190999916 0.0670820393249937 0.126491106406735 2.73861278752583 0.116189500386222 0.102469507659596 0.0948683298050513 0.1 0.0774596669241483 0.0692820323027551 0.0574456264653803 0.0524404424085076 0.0447213595499958 0.158113883008419 0.0316227766016838; 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838; 0.0118321595661992 0.01 0.335410196624968 0.335410196624968 61.2372435695795 27.3861278752583 0.0353553390593274 0.0458257569495584 0.0316227766016838 0.0494974746830583 0.0707106781186547 0.0565685424949238 3.46410161513776 5.29150262212918 0.0565685424949238 0.00111803398874989 0.0354964786985977 0.0278208554864871 0.0256904651573302 0.0707106781186547 0.0241867732448956 0.0217944947177034 0.0185202591774521 0.0173205080756888 0.0158113883008419 47.4341649025257 0.316227766016838; 0.00547722557505166 0.00223606797749979 0.00790569415042094 0.013228756555323 0.0707106781186547 0.0707106781186547 0.00244948974278318 0.00387298334620742 0.000316227766016838 0.006 0.00346410161513776 0.00812403840463596 0.0790569415042095 0.1 0.0116189500386222 0.0106066017177982 0.0144913767461894 0.0149666295470958 0.015 0.0158113883008419 0.0141421356237309 0.0132664991614216 0.0130384048104053 0.012369316876853 0.0116189500386222 8.36660026534076 0.316227766016838; 0.4 1.36930639376292 0.223606797749979 0.2 3.46410161513776 8.66025403784439 6 7.41619848709567 0.223606797749979 8.83176086632785 2.23606797749979 10.2469507659596 1.34164078649987 1.73205080756888 10.9544511501033 8.06225774829855 10.2469507659596 9.2951600308978 7.88035532193822 6 6.557438524302 5.17010638188423 3.73496987939662 2.56904651573303 1.67332005306815 5.47722557505166 2.23606797749979; 0.632455532033676 0.632455532033676 18.02775637732 31.6227766016838 22.3606797749979 22.3606797749979 1.26491106406735 3.46410161513776 0.1 2 4.47213595499958 2.82842712474619 948.683298050515 948.683298050515 7.74596669241484 2.23606797749979 20.4939015319192 28.4604989415154 44.1588043316392 67.0820393249938 73.4846922834954 110 139.64240043769 173.205080756888 223.606797749979 0.1 63.2455532033676; 0.0045 0.408 156 1.29 0.226 0.226 2.05 2.44 0.5 2.85475042692001 2 3.34 0.1 10 4.22 3.78 4.67 4.58421203698084 4.5 4.3 4.12431812546026 3.78 3.34496636754392 2.96 2.61933874283863 0.1 0.0001; 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5; 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5; 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5 1.0e-5]; + KDs_2 = [0.0632455532033676 3.46410161513776 324.037034920393 13.228756555323 0.836660026534076 0.836660026534076 1414.2135623731 1174.73401244707 1 948.683298050515 1 612.372435695795 0.173205080756888 11.180339887499 300 34.6410161513775 116.189500386222 58.309518948453 31.6227766016838 31.6227766016838 18.9736659610103 13.4164078649987 10.2469507659596 7.74596669241484 6.70820393249937 1 57.0087712549569; 0.122474487139159 0.223606797749979 0.14142135623731 0.1 0.559016994374947 0.273861278752583 0.335410196624968 0.412310562561766 0.223606797749979 0.574456264653803 0.273861278752583 0.866025403784439 0.4 0.547722557505166 1.36930639376292 1 1.80277563773199 2.21359436211787 2.29128784747792 2.12132034355964 2.08566536146142 1.8165902124585 1.51657508881031 1.40712472794703 1.30384048104053 5.47722557505166 10; 0.00316227766016838 0.0790569415042095 0.790569415042095 0.387298334620742 0.00353553390593274 0.00707106781186547 3.87298334620742 7.07106781186548 0.316227766016838 11.0679718105893 1.93649167310371 15 0.193649167310371 0.0707106781186547 17.3205080756888 11.180339887499 18.3711730708738 17.8885438199983 15.6524758424985 12.2474487139159 13.4164078649987 11.180339887499 8.94427190999916 6.70820393249937 4.47213595499958 0.223606797749979 0.316227766016838; 3.16227766016838 5.47722557505166 0.0474341649025257 0.0707106781186547 0.474341649025257 0.316227766016838 0.0591607978309961 0.06 0.0316227766016838 0.0648074069840786 0.234520787991171 0.0670820393249937 0.229128784747792 0.316227766016838 0.0707106781186547 0.13228756555323 0.0741619848709566 0.0774596669241483 0.0793725393319377 0.111803398874989 0.0774596669241483 0.0747663025700749 0.0734846922834953 0.0726636084983398 0.0707106781186547 3.16227766016838 9.21954445729289; 0.106066017177982 0.0223606797749979 0.193649167310371 0.632455532033676 0.0316227766016838 0.0316227766016838 0.0741619848709566 0.0848528137423857 0.0316227766016838 0.0989949493661167 0.187082869338697 0.12 0.0612372435695794 0.0774596669241483 0.173205080756888 0.0316227766016838 0.403112887414927 0.670820393249937 1.14564392373896 0.987420882906575 1.58113883008419 2.29128784747792 2.80624304008046 3.16227766016838 3.51781181986757 0.447213595499958 2.23606797749979; 0.0316227766016838 0.0547722557505166 0.0707106781186547 0.0447213595499958 0.111803398874989 0.0707106781186547 0.13228756555323 0.212132034355964 0.316227766016838 0.290688837074973 0.287228132326901 0.401870625948202 0.234520787991171 0.31224989991992 0.5 0.692820323027551 0.725603197346869 0.848528137423857 0.953939201416946 1.09544511501033 1.06348483769163 1.14105214604767 1.18215904175369 1.1861703081767 1.18321595661992 1.93649167310371 7.74596669241484; 0.00316227766016838 0.00447213595499958 0.0158113883008419 0.0237170824512628 4.47213595499958 6.70820393249937 0.0223606797749979 0.0306186217847897 0.0158113883008419 0.0367423461417476 0.00632455532033676 0.0412310562561766 1.22474487139159 1.5 0.0547722557505166 0.0447213595499958 0.0790569415042095 0.09 0.0989949493661167 0.0707106781186547 0.101980390271856 0.111803398874989 0.117260393995586 0.117473401244707 0.109544511501033 11.180339887499 1.73205080756888; 0.00836660026534076 0.0154919333848297 0.0474341649025257 0.0935414346693485 0.0264575131106459 0.05 0.0264575131106459 0.0519615242270663 0.0316227766016838 0.102469507659596 0.0173205080756888 0.205548047910945 0.547722557505166 0.387298334620742 0.433012701892219 0.5 1.58113883008419 3.16227766016838 5.45435605731786 6.12372435695795 8.2915619758885 9.89949493661167 10.0623058987491 9.35414346693486 7.74596669241484 2.64575131106459 3.16227766016838; 0.948683298050514 5.47722557505166 0.0173205080756888 0.0273861278752583 0.0158113883008419 0.00316227766016838 0.111803398874989 0.0670820393249937 0.223606797749979 0.0447213595499958 2.23606797749979 0.03 0.0223606797749979 0.0223606797749979 0.0212132034355964 3.35410196624969 0.0167332005306815 0.015 0.0144913767461894 0.0254950975679639 0.0162018517460197 0.0189736659610103 0.02 0.0232379000772445 0.0273861278752583 0.111803398874989 0.02; 0.0387298334620742 0.158113883008419 0.0935414346693485 0.13228756555323 0.136930639376291 0.122474487139159 0.14142135623731 0.189736659610103 0.316227766016838 0.234520787991171 0.0223606797749979 0.267394839142419 0.220794021658196 0.158113883008419 0.3 0.254950975679639 0.324037034920393 0.31224989991992 0.301662062579967 0.335410196624968 0.268328157299975 0.245967477524977 0.234520787991171 0.223606797749979 0.212132034355964 22.3606797749979 3.16227766016838; 0.013228756555323 0.0122474487139159 0.0154919333848297 0.0223606797749979 0.0193649167310371 0.0273861278752583 0.00273861278752583 0.00489897948556636 0.0187082869338697 0.00790569415042094 0.01 0.0116189500386222 0.0387298334620742 0.0223606797749979 0.0152970585407783 0.0187616630392937 0.0228035085019828 0.0270554985169374 0.0324037034920393 0.0346410161513775 0.0369864840178138 0.04 0.0439545219516718 0.0474341649025257 0.0529150262212918 0.16431676725155 0.433012701892219; 0.0244948974278318 0.0316227766016838 0.1 0.0316227766016838 0.05 0.0670820393249937 0.00894427190999916 0.0204939015319192 0.111803398874989 0.0363318042491699 0.0193649167310371 0.0574456264653803 0.0346410161513775 0.0547722557505166 0.0812403840463596 0.114017542509914 0.157480157480236 0.207846096908265 0.256124969497314 0.282842712474619 0.311126983722081 0.351425667816112 0.4 0.458257569495584 0.5 1.5 5.47722557505166; 0.0447213595499958 0.387298334620742 0.0316227766016838 0.0547722557505166 0.0612372435695794 0.0447213595499958 0.2 0.154919333848297 0.5 0.122474487139159 2.44948974278318 0.106066017177982 0.0223606797749979 0.0324037034920393 0.0774596669241483 1.1180339887499 0.0447213595499958 0.0367423461417476 0.03 0.0292403830344269 0.0256904651573302 0.0212132034355964 0.0182345825288105 0.0169705627484771 0.0154919333848297 0.158113883008419 0.0316227766016838; 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838; 0.01 0.00866025403784438 0.162018517460196 0.2 33.5410196624969 22.3606797749979 0.0111803398874989 0.0141421356237309 0.0273861278752583 0.0156524758424985 0.0316227766016838 0.0172481883106603 2.23606797749979 3.74165738677394 0.0184932420089069 0.000632455532033676 0.0194935886896179 0.0204939015319192 0.022248595461287 0.0273861278752583 0.024 0.0264575131106459 0.0291204395571221 0.0314642654451045 0.033466401061363 3.16227766016838 0.273861278752583; 0.00547722557505166 0.00223606797749979 0.00790569415042094 0.013228756555323 0.0707106781186547 0.0707106781186547 0.00244948974278318 0.00387298334620742 0.000316227766016838 0.006 0.00346410161513776 0.00812403840463596 0.0790569415042095 0.1 0.0116189500386222 0.0106066017177982 0.0144913767461894 0.0149666295470958 0.015 0.0158113883008419 0.0141421356237309 0.0132664991614216 0.0130384048104053 0.012369316876853 0.0116189500386222 8.36660026534076 0.316227766016838; 0.4 1.36930639376292 0.223606797749979 0.2 3.46410161513776 8.66025403784439 6 7.41619848709567 0.223606797749979 8.83176086632785 2.23606797749979 10.2469507659596 1.34164078649987 1.73205080756888 10.9544511501033 8.06225774829855 10.2469507659596 9.2951600308978 7.88035532193822 6 6.557438524302 5.17010638188423 3.73496987939662 2.56904651573303 1.67332005306815 5.47722557505166 2.23606797749979; 0.632455532033676 0.632455532033676 18.02775637732 31.6227766016838 22.3606797749979 22.3606797749979 1.26491106406735 3.46410161513776 0.1 2 4.47213595499958 2.82842712474619 948.683298050515 948.683298050515 7.74596669241484 2.23606797749979 20.4939015319192 28.4604989415154 44.1588043316392 67.0820393249938 73.4846922834954 110 139.64240043769 173.205080756888 223.606797749979 0.1 63.2455532033676] + KDs_3 = [0.0632455532033676 3.46410161513776 324.037034920393 13.228756555323 0.836660026534076 0.836660026534076 1414.2135623731 1174.73401244707 1 948.683298050515 1 612.372435695795 0.173205080756888 11.180339887499 300 34.6410161513775 116.189500386222 58.309518948453 31.6227766016838 31.6227766016838 18.9736659610103 13.4164078649987 10.2469507659596 7.74596669241484 6.70820393249937 1 57.0087712549569; 0.223606797749979 0.31224989991992 0.0141421356237309 0.0193649167310371 0.158113883008419 0.220794021658196 0.0724568837309472 0.144913767461894 0.0632455532033676 0.234520787991171 0.3 0.346410161513775 0.335410196624968 0.391152144312159 0.447213595499958 0.524404424085076 0.589915248150105 0.634822809924155 0.66932802122726 0.612372435695794 0.681175454637056 0.648074069840786 0.603738353924943 0.53619026473818 0.458257569495584 4.89897948556636 3.87298334620742; 0.00316227766016838 0.0790569415042095 0.790569415042095 0.387298334620742 0.00353553390593274 0.00707106781186547 3.87298334620742 7.07106781186548 0.316227766016838 11.0679718105893 1.93649167310371 15 0.193649167310371 0.0707106781186547 17.3205080756888 11.180339887499 18.3711730708738 17.8885438199983 15.6524758424985 12.2474487139159 13.4164078649987 11.180339887499 8.94427190999916 6.70820393249937 4.47213595499958 0.223606797749979 0.316227766016838; 2 4.47213595499958 0.0474341649025257 0.0707106781186547 0.273861278752583 0.187082869338697 0.0144913767461894 0.0189736659610103 0.0316227766016838 0.0244948974278318 0.234520787991171 0.031224989991992 0.0547722557505166 0.0547722557505166 0.0387298334620742 0.0418330013267038 0.0447213595499958 0.0519615242270663 0.0591607978309961 0.066332495807108 0.0692820323027551 0.0747663025700749 0.0793725393319377 0.0819756061276768 0.0866025403784439 3.16227766016838 9.21954445729289; 0.106066017177982 0.0223606797749979 0.193649167310371 0.632455532033676 0.0316227766016838 0.0316227766016838 0.0741619848709566 0.0848528137423857 0.0316227766016838 0.0989949493661167 0.187082869338697 0.12 0.0612372435695794 0.0774596669241483 0.173205080756888 0.0316227766016838 0.403112887414927 0.670820393249937 1.14564392373896 0.987420882906575 1.58113883008419 2.29128784747792 2.80624304008046 3.16227766016838 3.51781181986757 0.447213595499958 2.23606797749979; 0.00948683298050514 0.005 0.00935414346693485 0.00707106781186547 0.00894427190999916 0.0180277563773199 0.0474341649025257 0.0866025403784439 0.0223606797749979 0.14456832294801 0.122474487139159 0.216333076527839 0.0707106781186547 0.187082869338697 0.301662062579967 0.387298334620742 0.460977222864644 0.522494019104525 0.551361950083609 0.6 0.585662018573853 0.608276253029822 0.608276253029822 0.585662018573853 0.559910707166777 1.58113883008419 4.74341649025257; 0.00316227766016838 0.00447213595499958 0.0158113883008419 0.0237170824512628 4.47213595499958 6.70820393249937 0.0223606797749979 0.0306186217847897 0.0158113883008419 0.0367423461417476 0.00632455532033676 0.0412310562561766 1.22474487139159 1.5 0.0547722557505166 0.0447213595499958 0.0790569415042095 0.09 0.0989949493661167 0.0707106781186547 0.101980390271856 0.111803398874989 0.117260393995586 0.117473401244707 0.109544511501033 11.180339887499 1.73205080756888; 0.000866025403784438 0.000387298334620741 0.002 0.00612372435695794 0.013228756555323 0.0180277563773199 0.00894427190999916 0.0219089023002067 0.0122474487139159 0.0565685424949238 0.013228756555323 0.122474487139159 0.3 0.252487623459052 0.264575131106459 0.424264068711928 0.821583836257749 1.3228756555323 1.93649167310371 2.32379000772445 2.54950975679639 3.1224989991992 3.51781181986757 4.02336923485777 4.47213595499958 2.44948974278318 3.74165738677394; 0.948683298050514 5.47722557505166 0.0173205080756888 0.0273861278752583 0.0158113883008419 0.00316227766016838 0.111803398874989 0.0670820393249937 0.223606797749979 0.0447213595499958 2.23606797749979 0.03 0.0223606797749979 0.0223606797749979 0.0212132034355964 3.35410196624969 0.0167332005306815 0.015 0.0144913767461894 0.0254950975679639 0.0162018517460197 0.0189736659610103 0.02 0.0232379000772445 0.0273861278752583 0.111803398874989 0.02; 0.0158113883008419 0.0158113883008419 0.0387298334620742 0.0387298334620742 0.1 0.1 0.0126491106406735 0.0273861278752583 0.0591607978309961 0.0474341649025257 0.0122474487139159 0.062449979983984 0.150831031289984 0.2 0.0724568837309471 0.0360555127546399 0.0774596669241483 0.0702851335632223 0.0666333249958307 0.0591607978309961 0.0620483682299543 0.0565685424949238 0.0547722557505166 0.05 0.0447213595499958 6.32455532033676 1.4142135623731; 0.00223606797749979 0.00223606797749979 0.000707106781186547 0.000707106781186547 0.00273861278752583 0.005 0.000223606797749979 0.000447213595499958 0.000193649167310371 0.000935414346693485 0.000109544511501033 0.00158113883008419 0.00316227766016838 0.00547722557505166 0.00250998007960223 0.00316227766016838 0.00424264068711928 0.0062928530890209 0.00908295106229247 0.00774596669241483 0.013228756555323 0.0193649167310371 0.0241867732448956 0.0273861278752583 0.0282842712474619 0.0774596669241483 0.287228132326901; 0.00547722557505166 0.00173205080756888 0.000316227766016838 0.0005 0.00346410161513776 0.00632455532033676 0.000836660026534075 0.00264575131106459 0.02 0.00547722557505166 0.0141421356237309 0.0109544511501033 0.0291547594742265 0.0547722557505166 0.0180277563773199 0.0223606797749979 0.0324037034920393 0.0412310562561766 0.05 0.0670820393249937 0.0619677335393187 0.0734846922834953 0.0916515138991168 0.107238052947636 0.120415945787923 1 1.22474487139159; 0.0547722557505166 0.335410196624968 0.1 0.0632455532033676 0.0447213595499958 0.0447213595499958 0.111803398874989 0.0894427190999916 0.591607978309962 0.0714142842854285 1.6583123951777 0.0591607978309961 0.01 0.0154919333848297 0.0489897948556635 0.707106781186547 0.0387298334620742 0.0338230690505755 0.03 0.0357071421427143 0.0267394839142419 0.0240831891575846 0.0222261107708929 0.0203469899493758 0.0178885438199983 0.0670820393249937 0.0316227766016838; 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838 0.00316227766016838; 0.01 0.00866025403784438 0.162018517460196 0.2 33.5410196624969 22.3606797749979 0.0111803398874989 0.0141421356237309 0.0273861278752583 0.0156524758424985 0.0316227766016838 0.0172481883106603 2.23606797749979 3.74165738677394 0.0184932420089069 0.000632455532033676 0.0194935886896179 0.0204939015319192 0.022248595461287 0.0273861278752583 0.024 0.0264575131106459 0.0291204395571221 0.0314642654451045 0.033466401061363 3.16227766016838 0.273861278752583; 0.00547722557505166 0.00223606797749979 0.00790569415042094 0.013228756555323 0.0707106781186547 0.0707106781186547 0.00244948974278318 0.00387298334620742 0.000316227766016838 0.006 0.00346410161513776 0.00812403840463596 0.0790569415042095 0.1 0.0116189500386222 0.0106066017177982 0.0144913767461894 0.0149666295470958 0.015 0.0158113883008419 0.0141421356237309 0.0132664991614216 0.0130384048104053 0.012369316876853 0.0116189500386222 8.36660026534076 0.316227766016838; 0.4 1.36930639376292 0.223606797749979 0.2 3.46410161513776 8.66025403784439 6 7.41619848709567 0.223606797749979 8.83176086632785 2.23606797749979 10.2469507659596 1.34164078649987 1.73205080756888 10.9544511501033 8.06225774829855 10.2469507659596 9.2951600308978 7.88035532193822 6 6.557438524302 5.17010638188423 3.73496987939662 2.56904651573303 1.67332005306815 5.47722557505166 2.23606797749979; 0.632455532033676 0.632455532033676 18.02775637732 31.6227766016838 22.3606797749979 22.3606797749979 1.26491106406735 0.1 0.316227766016838 2 4.47213595499958 2.82842712474619 948.683298050515 948.683298050515 7.74596669241484 2.23606797749979 20.4939015319192 28.4604989415154 44.1588043316392 67.0820393249938 73.4846922834954 110 139.64240043769 173.205080756888 223.606797749979 0.1 63.2455532033676; 0.0045 0.408 156 1.29 0.226 0.226 2.05 2.44 0.5 2.89 2 3.78 0.1 10 4.22 3.78 4.67 4.59 4.5 4.3 4.14 3.78 3.37 2.96 2.55 0.1 0.0001] + + # phase_name = (mineral_name_convertor(ph_1),mineral_name_convertor(ph_2),mineral_name_convertor(ph_3)) + phase_name = (ph_1,ph_2,ph_3) + KDs = (KDs_1,KDs_2,KDs_3) + + KDs_dtb = KDs_database(infos, element_name,conditions,phase_name,KDs) + + return KDs_dtb +end + + +function adjust_chemical_system( KDs_dtb :: KDs_database, + bulk_TE :: Vector{Float64}, + elem_TE :: Vector{String}) + + C0_TE_idx = [findfirst(isequal(x), elem_TE) for x in KDs_dtb.element_name] + C0_TE = bulk_TE[C0_TE_idx] + + return C0_TE +end + + +struct out_tepm + elements :: Union{Nothing, Vector{String}} + C0 :: Union{Nothing, Vector{Float64}} + Cliq :: Union{Nothing, Vector{Float64}} + Csol :: Union{Nothing, Vector{Float64}} + Cmin :: Union{Nothing, Matrix{Float64}} + Czrc :: Union{Nothing, Vector{Float64}} + ph_TE :: Union{Nothing, Vector{String}} + ph_wt_norm :: Union{Nothing, Vector{Float64}} + liq_wt_norm :: Union{Nothing, Float64} + Cliq_Zr :: Union{Nothing, Float64} + Sat_zr_liq :: Union{Nothing, Float64} + zrc_wt :: Union{Nothing, Float64} + bulk_cor_wt :: Union{Nothing, Vector{Float64}} +end + + +function compute_partitioning( cond :: Int64, + C0 :: Vector{Float64}, + ph :: Vector{String}, + ph_wt :: Vector{Float64}, + liq_wt :: Float64, + KDs_dtb :: KDs_database) + + + TE_ph = intersect(ph,KDs_dtb.phase_name[cond]); + + # get indexes of the phase with respect to MAGEMin output and TE_database + MM_ph_idx = [findfirst(isequal(x), ph) for x in TE_ph]; + TE_ph_idx = [findfirst(isequal(x), KDs_dtb.phase_name[cond]) for x in TE_ph]; + + # normalize phase fractions + sum_ph_frac = sum(ph_wt[MM_ph_idx]); + liq_wt_norm = liq_wt/(sum_ph_frac+liq_wt); + ph_wt_norm = ph_wt[MM_ph_idx]./sum_ph_frac; + ph_TE = ph[MM_ph_idx]; + + # compute bulk distributiion coefficient + D = KDs_dtb.KDs[cond][TE_ph_idx,:]'*ph_wt_norm; + + Cliq = C0 ./ (D .+ liq_wt_norm.*(1.0 .- D)); + Csol = (C0 .- Cliq .* liq_wt_norm) ./ (1.0 .- liq_wt_norm) + Cmin = similar(KDs_dtb.KDs[cond][TE_ph_idx,:]); + + for i = 1:length(ph_wt_norm) + Cmin[i,:] = KDs_dtb.KDs[cond][TE_ph_idx[i],:] .* Cliq; end - return TE_dtb + return Cliq, Cmin, Csol, ph_TE, ph_wt_norm, liq_wt_norm end -""" - Routine to compute the TE partitioning -""" -function compute_TE_partitioning( C0 :: Vector{Float64}, - out :: gmin_struct{Float64, Int64}, - dtb :: String; - TE_db :: String = "OL" ) - - if out.frac_M > 0.0 && out.frac_S > 0.0 - SiO2_M = out.bulk_M_wt[1] - liq_wt = out.frac_M - if TE_db == "OL" - if SiO2_M > 63.0 - TE_dtb = get_TE_database("TE_OL_felsic") - elseif SiO2_M > 52.0 && SiO2_M <= 63.0 - TE_dtb = get_TE_database("TE_OL_intermediate") - elseif SiO2_M <= 52.0 - TE_dtb = get_TE_database("TE_OL_mafic") + +function TE_prediction( C0 :: Vector{Float64}, + KDs_dtb :: KDs_database, + ZrSat_model:: String, + out :: MAGEMin_C.gmin_struct{Float64, Int64}, + dtb :: String ) + + ox_id = findfirst(out.oxides .== KDs_dtb.conditions[1])[1] + ox_non_H2O = findall(out.oxides .!= "H2O") + + ox_M = out.bulk_M_wt[ox_id]./sum(out.bulk_M_wt[ox_non_H2O]) + liq_wt = out.frac_M_wt + sol_wt = out.frac_S_wt + elements = KDs_dtb.element_name + if liq_wt > 0.0 && liq_wt < 1.0 && sol_wt > 0.0 + n_cond = length(KDs_dtb.conditions[2]) + cond = -1 + for i=1:n_cond + if ox_M > KDs_dtb.conditions[2][i][1] && ox_M < KDs_dtb.conditions[2][i][2] + cond = i + break end - else - print("Trace element database $TE_db is not implemented!\n") end ph, ph_wt = mineral_classification(out, dtb); - TE_ph = intersect(ph,TE_dtb.phase_name); + Cliq, Cmin, Csol, ph_TE, ph_wt_norm, liq_wt_norm, = compute_partitioning( cond, + C0, + ph, + ph_wt, + liq_wt, + KDs_dtb) + + id_Zr = findfirst(KDs_dtb.element_name .== "Zr")[1] + Cliq_Zr = Cliq[id_Zr] + + Sat_zr_liq = zirconium_saturation( out; + model = ZrSat_model) + + if Cliq_Zr > Sat_zr_liq + zrc_wt, SiO2_zrc_wt, O_wt = adjust_bulk_4_zircon(Cliq_Zr, Sat_zr_liq) + + push!(ph,"zrn") + push!(ph_wt,zrc_wt/100.0) + + ph_wt = ph_wt ./(sum(ph_wt)) - # get indexes of the phase with respect to MAGEMin output and TE_database - MM_ph_idx = [findfirst(isequal(x), ph) for x in TE_ph]; - TE_ph_idx = [findfirst(isequal(x), TE_dtb.phase_name) for x in TE_ph]; + Cliq, Cmin, Csol, ph_TE, ph_wt_norm, liq_wt_norm, = compute_partitioning( cond, + C0, + ph, + ph_wt, + liq_wt, + KDs_dtb) - # normalize phase fractions - sum_ph_frac = sum(ph_wt[MM_ph_idx]); - liq_wt_norm = liq_wt/(sum_ph_frac+liq_wt); - ph_wt_norm = ph_wt[MM_ph_idx]./sum_ph_frac; - ph_TE = ph[MM_ph_idx]; + Cliq_Zr = Cliq[id_Zr] + Sat_zr_liq = zirconium_saturation( out; + model = ZrSat_model) - # compute bulk distributiion coefficient - D = TE_dtb.Kds[TE_ph_idx,:]'*ph_wt_norm; + SiO2_id = findall(out.oxides .== "SiO2")[1] - Cliq = C0 ./ (D .+ liq_wt_norm.*(1.0 .- D)); - Cmin = similar(TE_dtb.Kds[TE_ph_idx,:]); + bulk_cor_wt = copy(out.bulk_wt) + bulk_cor_wt[SiO2_id] = out.bulk_wt[SiO2_id] - SiO2_zrc_wt + bulk_cor_wt ./= sum(bulk_cor_wt) - for i = 1:length(ph_wt_norm) - Cmin[i,:] = TE_dtb.Kds[TE_ph_idx[i],:] .* Cliq; + else + zrc_wt, bulk_cor_wt = nothing, nothing end - end + elseif liq_wt == 0.0 + Csol = C0 + Cliq, Cmin, ph_TE, ph_wt_norm, liq_wt_norm, Cliq_Zr, Sat_zr_liq, zrc_wt, bulk_cor_wt = nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing + elseif liq_wt == 1.0 || (sol_wt == 0.0 && liq_wt > 0.0) #latter means there is fluid + melt + Cliq = C0 + Csol = nothing + id_Zr = findfirst(KDs_dtb.element_name .== "Zr")[1] - return Cliq, Cmin, ph_TE, ph_wt_norm, liq_wt_norm -end + Cliq_Zr = Cliq[id_Zr] + Cmin, ph_TE, ph_wt_norm, zrc_wt = nothing, nothing, nothing, nothing + Sat_zr_liq = zirconium_saturation( out; + model = ZrSat_model) + if Cliq_Zr > Sat_zr_liq + zrc_wt, SiO2_zrc_wt, O_wt = adjust_bulk_4_zircon(Cliq_Zr, Sat_zr_liq) -TE_db = get_TE_database("TE_OL_felsic") -components = ["SiO2","TiO2","Al2O3","O","FeO","MnO","MgO","CaO","K2O","Na2O","H2O","Li","Be","B","Sc","V","Cr","Ni","Cu","Zn","Rb","Sr","Y","Zr","Nb","Cs","Ba","La","Ce","Pr","Nd","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","Pb","Th","U"] -Ton = [69.20111492,0.385480479,15.36802994,0.08,2.835169639,0.052720082,1.173738821,3.375444597,1.767078881,4.503466421,2,5.524762752,39.87762745,49.26129576,15.04133712,15.15293785,63.74940359,64.5373756,484.5753262,7.459369222,150.2402264,5.430265866,2.474580356,499.0093971,29.06202556,53.78324684,5.797620808,20.54095673,3.317738791,0.961548784,2.418869313,0.315058021,1.593968618,0.296925128,0.799927928,0.128677572,0.713354049,0.133287066,3.95645511,0.722386615,15.17205921,5.69152854,1.064824497] -Bas = [50.60777553,0.953497243,13.70435413,0.19,11.28130762,0.202560796,8.496024312,9.502380068,0.700881685,2.07434927,4,29.14258603,0.434482763,29.69200003,38.23663423,257.4346716,529.333066,208.2057375,88.87615683,91.7592182,16.60777308,163.4533209,20.74016207,66.90677472,3.808354064,1.529226981,122.8449739,6.938172601,16.04827796,2.253943183,10.18276823,3.3471043,0.915941652,3.28230146,1.417695298,3.851230952,0.914966282,2.20425,0.343734976,2.136202593,0.323405135,1.841502082,0.330971265,5.452969044,1.074692058,0.290233271] + SiO2_id = findall(out.oxides .== "SiO2")[1] -C0_TE_idx = [findfirst(isequal(x), components) for x in TE_db.element_name] -C0_TE = Bas[C0_TE_idx] + bulk_cor_wt = copy(out.bulk_wt) + bulk_cor_wt[SiO2_id] = out.bulk_wt[SiO2_id] - SiO2_zrc_wt + bulk_cor_wt ./= sum(bulk_cor_wt) + else + zrc_wt, bulk_cor_wt = nothing, nothing + end -# using MAGEMin_C -# dtb = "mb" -# data = Initialize_MAGEMin(dtb, verbose=true); -# test = 0 -# data = use_predefined_bulk_rock(data, test); -# P = 10.0 -# T = 1000.0 -# out = point_wise_minimization(P,T, data); + liq_wt_norm = 1.0 + else + print("unrecognized case!\n") + Cliq, Csol, Cmin, ph_TE, ph_wt_norm, liq_wt_norm, Cliq_Zr, zrc_wt, bulk_cor_wt, Sat_zr_liq = nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing + end -# Cliq, Cmin, ph_TE, ph_wt_norm, liq_wt_norm = compute_TE_partitioning(C0_TE,out,dtb) + out_TE = out_tepm(elements, C0, Cliq, Csol, Cmin, nothing, ph_TE, ph_wt_norm, liq_wt_norm, Cliq_Zr, Sat_zr_liq, zrc_wt, bulk_cor_wt) -# C_zr_liq = zirconium_saturation( out; -# model="CB") \ No newline at end of file + return out_TE +end diff --git a/julia/Zircon_saturation.jl b/julia/Zircon_saturation.jl index 1e02bd7..eb71ba6 100644 --- a/julia/Zircon_saturation.jl +++ b/julia/Zircon_saturation.jl @@ -1,11 +1,7 @@ # Routine to compute zircon saturation and adjust bulk-rock composition when zircon crystallizes # NR 12/04/2023 -# include("TE_routine.jl") - -# using MAGEMin_C - -function zirconium_saturation( out :: gmin_struct{Float64, Int64}; +function zirconium_saturation( out :: MAGEMin_C.gmin_struct{Float64, Int64}; model :: String = "WH" ) if out.frac_M > 0.0 @@ -63,16 +59,18 @@ function zirconium_saturation( out :: gmin_struct{Float64, Int64}; return C_zr_liq end +function adjust_bulk_4_zircon( zr_liq :: Float64, + sat_liq :: Float64 ) + SiO2_wt = 0.0 + O2_wt = 0.0 + zircon_wt = 0.0 + zircon_excess = (zr_liq - sat_liq)/1e4 -# dtb = "ig" -# data = Initialize_MAGEMin(dtb, verbose=true); -# test = 0 -# data = use_predefined_bulk_rock(data, test); -# P = 2.0 -# T = 1800.0 -# out = point_wise_minimization(P,T, data); + zircon_wt = zircon_excess*0.497644 + SiO2_wt = (zircon_wt *0.327765) + O2_wt = (zircon_wt *0.174570) -# C_zr_liq = zirconium_saturation( out; -# model="CB") + return zircon_wt ,SiO2_wt ,O2_wt +end diff --git a/julia/export2CSV.jl b/julia/export2CSV.jl new file mode 100644 index 0000000..6a17fda --- /dev/null +++ b/julia/export2CSV.jl @@ -0,0 +1,352 @@ + +""" +MAGEMin_dataTE2dataframe( out:: Union{Vector{out_tepm}, out_tepm},dtb,fileout) + + Transform MAGEMin trace-element output into a dataframe for quick(ish) save + +""" +function MAGEMin_dataTE2dataframe( out :: Union{Vector{gmin_struct{Float64, Int64}}, gmin_struct{Float64, Int64}}, + out_te :: Union{Vector{out_tepm}, out_tepm}, + dtb, + fileout) + + # here we fill the dataframe with the all minimized point entries + if typeof(out) == MAGEMin_C.gmin_struct{Float64, Int64} + out = [out] + end + if typeof(out_te) == MAGEMin_C.out_tepm + out_te = [out_te] + end + np = length(out) + + db_in = retrieve_solution_phase_information(dtb) + datetoday = string(Dates.today()) + rightnow = string(Dates.Time(Dates.now())) + + metadata = "# MAGEMin_TE " * " $(out[1].MAGEMin_ver);" * datetoday * ", " * rightnow * "; using database " * db_in.db_info * "\n" + + # Here we create the dataframe's header: + MAGEMin_db = DataFrame( Symbol("point[#]") => Int64[], + Symbol("X[0.0-1.0]") => Float64[], + Symbol("P[kbar]") => Float64[], + Symbol("T[°C]") => Float64[], + Symbol("phase") => String[], + Symbol("mode[wt%]") => Float64[], + Symbol("Zr_sat[μg/g]") => Float64[], + ) + for i in out[1].oxides + col = i*"_cor[wt%]" + MAGEMin_db[!, col] = Float64[] + end + + for i in out_te[1].elements + col = i*"_[μg/g]" + MAGEMin_db[!, col] = Float64[] + end + + print("\noutput path: $(pwd())\n") + @showprogress "Saving data to csv..." for k=1:np + + part_1 = Dict( "point[#]" => k, + "X[0.0-1.0]" => out[k].X[1], + "P[kbar]" => out[k].P_kbar, + "T[°C]" => out[k].T_C, + "phase" => "system", + "mode[wt%]" => 100.0, + "Zr_sat[μg/g]" => "-") + + if ~isnothing(out_te[k].bulk_cor_wt) + part_2 = Dict( (out[1].oxides[j]*"_cor[wt%]" => out_te[k].bulk_cor_wt[j]*100.0) + for j in eachindex(out[1].oxides)) + else + part_2 = Dict( (out[1].oxides[j]*"_cor[wt%]" => "-") + for j in eachindex(out[1].oxides)) + end + + part_3 = Dict( ( out_te[1].elements[j]*"_[μg/g]" => out_te[k].C0[j]) + for j in eachindex(out_te[1].elements)) + + row = merge(part_1,part_2,part_3) + + push!(MAGEMin_db, row, cols=:union) + + # liquid + if ~isnothing(out_te[k].liq_wt_norm) + part_1 = Dict( "point[#]" => k, + "X[0.0-1.0]" => out[k].X[1], + "P[kbar]" => out[k].P_kbar, + "T[°C]" => out[k].T_C, + "phase" => "liq", + "mode[wt%]" => out_te[k].liq_wt_norm, + "Zr_sat[μg/g]" => out_te[k].Sat_zr_liq) + + part_2 = Dict( (out[1].oxides[j]*"_cor[wt%]" => "-") + for j in eachindex(out[1].oxides)) + + part_3 = Dict( ( out_te[1].elements[j]*"_[μg/g]" => out_te[k].Cliq[j]) + for j in eachindex(out_te[1].elements)) + + row = merge(part_1,part_2,part_3) + + push!(MAGEMin_db, row, cols=:union) + end + + # solid + if ~isnothing(out_te[k].Csol) + if ~isnothing(out_te[k].liq_wt_norm) + sol_mode = (1.0-out_te[k].liq_wt_norm)*100.0 + else + sol_mode = 100.0 + end + part_1 = Dict( "point[#]" => k, + "X[0.0-1.0]" => out[k].X[1], + "P[kbar]" => out[k].P_kbar, + "T[°C]" => out[k].T_C, + "phase" => "sol", + "mode[wt%]" => sol_mode, + "Zr_sat[μg/g]" => "-") + + part_2 = Dict( (out[1].oxides[j]*"_cor[wt%]" => "-") + for j in eachindex(out[1].oxides)) + + part_3 = Dict( ( out_te[1].elements[j]*"_[μg/g]" => out_te[k].Csol[j]) + for j in eachindex(out_te[1].elements)) + + row = merge(part_1,part_2,part_3) + + push!(MAGEMin_db, row, cols=:union) + end + + if ~isnothing(out_te[k].ph_TE) + nph = length(out_te[k].ph_TE) + for i=1:nph + part_1 = Dict( "point[#]" => k, + "X[0.0-1.0]" => out[k].X[1], + "P[kbar]" => out[k].P_kbar, + "T[°C]" => out[k].T_C, + "phase" => out_te[k].ph_TE, + "mode[wt%]" => out_te[k].ph_wt_norm, + "Zr_sat[μg/g]" => "-") + + part_2 = Dict( (out[1].oxides[j]*"_cor[wt%]" => "-") + for j in eachindex(out[1].oxides)) + + part_3 = Dict( ( out_te[1].elements[j]*"_[μg/g]" => out_te[k].Cmin[i,j]) + for j in eachindex(out_te[1].elements)) + + row = merge(part_1,part_2,part_3) + + push!(MAGEMin_db, row, cols=:union) + + end + end + + end + + + meta = fileout*"_metadata.txt" + filename = fileout*".csv" + CSV.write(filename, MAGEMin_db) + + open(meta, "w") do file + write(file, metadata) + end + + return nothing +end + +""" +MAGEMin_data2dataframe( out:: Union{Vector{MAGEMin_C.gmin_struct{Float64, Int64}}, MAGEMin_C.gmin_struct{Float64, Int64}}) + + Transform MAGEMin output into a dataframe for quick(ish) save + +""" +function MAGEMin_data2dataframe( out:: Union{Vector{gmin_struct{Float64, Int64}}, gmin_struct{Float64, Int64}},dtb,fileout) + + # here we fill the dataframe with the all minimized point entries + if typeof(out) == MAGEMin_C.gmin_struct{Float64, Int64} + out = [out] + end + np = length(out) + + db_in = retrieve_solution_phase_information(dtb) + datetoday = string(Dates.today()) + rightnow = string(Dates.Time(Dates.now())) + + metadata = "# MAGEMin " * " $(out[1].MAGEMin_ver);" * datetoday * ", " * rightnow * "; using database " * db_in.db_info * "\n" + + # Here we create the dataframe's header: + MAGEMin_db = DataFrame( Symbol("point[#]") => Int64[], + Symbol("X[0.0-1.0]") => Float64[], + Symbol("P[kbar]") => Float64[], + Symbol("T[°C]") => Float64[], + Symbol("phase") => String[], + Symbol("mode[mol%]") => Float64[], + Symbol("mode[wt%]") => Float64[], + Symbol("log10(fO2)") => Float64[], + Symbol("log10(dQFM)") => Float64[], + Symbol("aH2O") => Float64[], + Symbol("aSiO2") => Float64[], + Symbol("aTiO2") => Float64[], + Symbol("aAl2O3") => Float64[], + Symbol("aMgO") => Float64[], + Symbol("aFeO") => Float64[], + Symbol("density[kg/m3]") => Float64[], + Symbol("volume[cm3/mol]") => Float64[], + Symbol("heatCapacity[kJ/K]")=> Float64[], + Symbol("alpha[1/K]") => Float64[], + Symbol("Entropy[J/K]") => Float64[], + Symbol("Enthalpy[J]") => Float64[], + Symbol("Vp[km/s]") => Float64[], + Symbol("Vs[km/s]") => Float64[], + Symbol("BulkMod[GPa]") => Float64[], + Symbol("ShearMod[GPa]") => Float64[], + ) + + for i in out[1].oxides + col = i*"[mol%]" + MAGEMin_db[!, col] = Float64[] + end + + for i in out[1].oxides + col = i*"[wt%]" + MAGEMin_db[!, col] = Float64[] + end + + + print("\noutput path: $(pwd())\n") + @showprogress "Saving data to csv..." for k=1:np + np = length(out[k].ph) + nss = out[k].n_SS + npp = np-nss + + part_1 = Dict( "point[#]" => k, + "X[0.0-1.0]" => out[k].X[1], + "P[kbar]" => out[k].P_kbar, + "T[°C]" => out[k].T_C, + "phase" => "system", + "mode[mol%]" => 100.0, + "mode[wt%]" => 100.0, + "log10(fO2)" => out[k].fO2[1], + "log10(dQFM)" => out[k].dQFM[1], + "aH2O" => out[k].aH2O, + "aSiO2" => out[k].aSiO2, + "aTiO2" => out[k].aTiO2, + "aAl2O3" => out[k].aAl2O3, + "aMgO" => out[k].aMgO, + "aFeO" => out[k].aFeO, + "density[kg/m3]" => out[k].rho, + "volume[cm3/mol]" => out[k].V, + "heatCapacity[kJ/K]"=> out[k].cp, + "alpha[1/K]" => out[k].alpha, + "Entropy[J/K]" => out[k].entropy, + "Enthalpy[J]" => out[k].enthalpy, + "Vp[km/s]" => out[k].Vp, + "Vs[km/s]" => out[k].Vs, + "BulkMod[GPa]" => out[k].bulkMod, + "ShearMod[GPa]" => out[k].shearMod ) + + part_2 = Dict( (out[1].oxides[j]*"[mol%]" => out[k].bulk[j]*100.0) + for j in eachindex(out[1].oxides)) + + part_3 = Dict( (out[1].oxides[j]*"[wt%]" => out[k].bulk_wt[j]*100.0) + for j in eachindex(out[1].oxides)) + + row = merge(part_1,part_2,part_3) + + push!(MAGEMin_db, row, cols=:union) + + for i=1:nss + part_1 = Dict( "point[#]" => k, + "X[0.0-1.0]" => out[k].X[1], + "P[kbar]" => out[k].P_kbar, + "T[°C]" => out[k].T_C, + "phase" => out[k].ph[i], + "mode[mol%]" => out[k].ph_frac[i].*100.0, + "mode[wt%]" => out[k].ph_frac_wt[i].*100.0, + "log10(fO2)" => "-", + "log10(dQFM)" => "-", + "aH2O" => "-", + "aSiO2" => "-", + "aTiO2" => "-", + "aAl2O3" => "-", + "aMgO" => "-", + "aFeO" => "-", + "density[kg/m3]" => out[k].SS_vec[i].rho, + "volume[cm3/mol]" => out[k].SS_vec[i].V, + "heatCapacity[kJ/K]"=> out[k].SS_vec[i].cp, + "alpha[1/K]" => out[k].SS_vec[i].alpha, + "Entropy[J/K]" => out[k].SS_vec[i].entropy, + "Enthalpy[J]" => out[k].SS_vec[i].enthalpy, + "Vp[km/s]" => out[k].SS_vec[i].Vp, + "Vs[km/s]" => out[k].SS_vec[i].Vs, + "BulkMod[GPa]" => out[k].SS_vec[i].bulkMod, + "ShearMod[GPa]" => out[k].SS_vec[i].shearMod ) + + part_2 = Dict( (out[1].oxides[j]*"[mol%]" => out[k].SS_vec[i].Comp[j]*100.0) + for j in eachindex(out[1].oxides)) + + part_3 = Dict( (out[1].oxides[j]*"[wt%]" => out[k].SS_vec[i].Comp_wt[j]*100.0) + for j in eachindex(out[1].oxides)) + + row = merge(part_1,part_2,part_3) + + push!(MAGEMin_db, row, cols=:union) + + end + + if npp > 0 + for i=1:npp + pos = i + nss + + part_1 = Dict( "point[#]" => k, + "X[0.0-1.0]" => out[k].X[1], + "P[kbar]" => out[k].P_kbar, + "T[°C]" => out[k].T_C, + "phase" => out[k].ph[pos], + "mode[mol%]" => out[k].ph_frac[pos].*100.0, + "mode[wt%]" => out[k].ph_frac_wt[pos].*100.0, + "log10(fO2)" => "-", + "log10(dQFM)" => "-", + "aH2O" => "-", + "aSiO2" => "-", + "aTiO2" => "-", + "aAl2O3" => "-", + "aMgO" => "-", + "aFeO" => "-", + "density[kg/m3]" => out[k].PP_vec[i].rho, + "volume[cm3/mol]" => out[k].PP_vec[i].V, + "heatCapacity[kJ/K]"=> out[k].PP_vec[i].cp, + "alpha[1/K]" => out[k].PP_vec[i].alpha, + "Entropy[J/K]" => out[k].PP_vec[i].entropy, + "Enthalpy[J]" => out[k].PP_vec[i].enthalpy, + "Vp[km/s]" => out[k].PP_vec[i].Vp, + "Vs[km/s]" => out[k].PP_vec[i].Vs, + "BulkMod[GPa]" => out[k].PP_vec[i].bulkMod, + "ShearMod[GPa]" => out[k].PP_vec[i].shearMod ) + + part_2 = Dict( (out[1].oxides[j]*"[mol%]" => out[k].PP_vec[i].Comp[j]*100.0) + for j in eachindex(out[1].oxides)) + + part_3 = Dict( (out[1].oxides[j]*"[wt%]" => out[k].PP_vec[i].Comp_wt[j]*100.0) + for j in eachindex(out[1].oxides)) + + row = merge(part_1,part_2,part_3) + + push!(MAGEMin_db, row, cols=:union) + + end + end + + end + + meta = fileout*"_metadata.txt" + filename = fileout*".csv" + CSV.write(filename, MAGEMin_db) + + open(meta, "w") do file + write(file, metadata) + end + + return nothing +end \ No newline at end of file diff --git a/src/MAGEMin.h b/src/MAGEMin.h index df1659e..774cbf3 100644 --- a/src/MAGEMin.h +++ b/src/MAGEMin.h @@ -109,6 +109,7 @@ typedef struct global_variables { int numPoint; /** the number of the current point */ int global_ite; /** global iteration increment */ int H2O_id; + int Al2O3_id; int K2O_id; int O_id; int TiO2_id; diff --git a/src/dump_function.c b/src/dump_function.c index 0ea14b8..20e7a7d 100644 --- a/src/dump_function.c +++ b/src/dump_function.c @@ -155,11 +155,18 @@ void fill_output_struct( global_variable gv, sum_Molar_mass_bulk = 0.0; sp[0].M_sys = 0.0; for (i = 0; i < nox; i++){ - sp[0].M_sys += z_b.bulk_rock[i]*z_b.masspo[i]; + sp[0].M_sys += z_b.bulk_rock[i]*z_b.masspo[i]; sp[0].bulk[i] = z_b.bulk_rock[i]; sp[0].gamma[i] = gv.gam_tot[i]; sp[0].bulk_wt[i] = z_b.bulk_rock[i]*z_b.masspo[i]; sum_Molar_mass_bulk += sp[0].bulk_wt[i]; + + sp[0].bulk_S[i] = 0.0; + sp[0].bulk_S_wt[i] = 0.0; + sp[0].bulk_M[i] = 0.0; + sp[0].bulk_M_wt[i] = 0.0; + sp[0].bulk_F[i] = 0.0; + sp[0].bulk_F_wt[i] = 0.0; } for (i = 0; i < nox; i++){ sp[0].bulk_wt[i] /= sum_Molar_mass_bulk; @@ -424,9 +431,11 @@ void fill_output_struct( global_variable gv, /* The following section normalizes the entries for S (solid), M (melt) and F (fluid) which are entries useful for geodynamic coupling */ // normalize rho_S and bulk_S - sp[0].rho_S /= sp[0].frac_S; - for (j = 0; j < gv.len_ox; j++){ - sp[0].bulk_S[j] /= sp[0].frac_S; + if (sp[0].frac_S > 0.0){ + sp[0].rho_S /= sp[0].frac_S; + for (j = 0; j < gv.len_ox; j++){ + sp[0].bulk_S[j] /= sp[0].frac_S; + } } sum = 0.0; @@ -487,11 +496,11 @@ void fill_output_struct( global_variable gv, sum = sp[0].frac_F_wt + sp[0].frac_M_wt + sp[0].frac_S_wt; + sp[0].frac_F_wt /= sum; sp[0].frac_M_wt /= sum; sp[0].frac_S_wt /= sum; - /* compute cp as J/K/kg for given bulk-rock composition */ double MolarMass_system = 0.0; for (int i = 0; i < gv.len_ox; i++){ diff --git a/src/gss_function.c b/src/gss_function.c index 2d4593f..da0c9b1 100644 --- a/src/gss_function.c +++ b/src/gss_function.c @@ -7653,36 +7653,34 @@ SS_ref G_SS_ig_EM_function( global_variable gv, } if (strcmp( name, "bi") == 0 ){ - // if no H2O, deactivate - if (z_b.bulk_rock[gv.H2O_id] == 0. || z_b.bulk_rock[gv.K2O_id] == 0.){ + if (z_b.bulk_rock[gv.H2O_id] == 0. || z_b.bulk_rock[gv.K2O_id] == 0. || z_b.bulk_rock[gv.Al2O3_id] == 0.){ SS_ref_db.ss_flags[0] = 0; } SS_ref_db = G_SS_ig_bi_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "cd") == 0){ - // if no H2O, deactivate - if (z_b.bulk_rock[gv.H2O_id] == 0.){ + if (z_b.bulk_rock[gv.H2O_id] == 0. || z_b.bulk_rock[gv.Al2O3_id] == 0. ){ SS_ref_db.ss_flags[0] = 0; } SS_ref_db = G_SS_ig_cd_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "cpx") == 0){ SS_ref_db = G_SS_ig_cpx_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "ep") == 0){ - // if no h2O, deactivate - if (z_b.bulk_rock[gv.H2O_id] == 0.){ + if (z_b.bulk_rock[gv.H2O_id] == 0. || z_b.bulk_rock[gv.Al2O3_id] == 0.){ SS_ref_db.ss_flags[0] = 0; } SS_ref_db = G_SS_ig_ep_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "fl") == 0){ - // if no H2O, deactivate if (z_b.bulk_rock[gv.H2O_id] == 0.){ SS_ref_db.ss_flags[0] = 0; } SS_ref_db = G_SS_ig_fl_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "g") == 0){ + if (z_b.bulk_rock[gv.Al2O3_id] == 0.){ + SS_ref_db.ss_flags[0] = 0; + } SS_ref_db = G_SS_ig_g_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "hb") == 0){ - // if no H2O, deactivate - if (z_b.bulk_rock[gv.H2O_id] == 0.){ + if (z_b.bulk_rock[gv.H2O_id] == 0. || z_b.bulk_rock[gv.Al2O3_id] == 0. ){ SS_ref_db.ss_flags[0] = 0; } SS_ref_db = G_SS_ig_hb_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } @@ -7692,13 +7690,11 @@ SS_ref G_SS_ig_EM_function( global_variable gv, } SS_ref_db = G_SS_ig_ilm_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "liq") == 0){ - /* turn of liquid when T < 600°C) */ if ( T < gv.min_melt_T){ SS_ref_db.ss_flags[0] = 0; } SS_ref_db = G_SS_ig_liq_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "mu") == 0){ - // if no H2O, deactivate if (z_b.bulk_rock[gv.H2O_id] == 0. || z_b.bulk_rock[gv.K2O_id] == 0.){ SS_ref_db.ss_flags[0] = 0; } @@ -7711,6 +7707,9 @@ SS_ref G_SS_ig_EM_function( global_variable gv, else if (strcmp( name, "fper") == 0 ){ SS_ref_db = G_SS_ig_fper_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps);} else if (strcmp( name, "fsp") == 0){ + if (z_b.bulk_rock[gv.Al2O3_id] == 0.){ + SS_ref_db.ss_flags[0] = 0; + } SS_ref_db = G_SS_ig_fsp_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } else if (strcmp( name, "spn") == 0){ SS_ref_db = G_SS_ig_spn_function(SS_ref_db, EM_database, gv.len_ox, z_b, eps); } diff --git a/src/initialize.h b/src/initialize.h index 145e5f0..797e296 100644 --- a/src/initialize.h +++ b/src/initialize.h @@ -180,7 +180,7 @@ global_variable global_variable_alloc( bulk_info *z_b ){ } strcpy(gv.outpath,"./output/"); /** define the outpath to save logs and final results file */ - strcpy(gv.version,"1.4.4 [28/04/2024]"); /** MAGEMin version */ + strcpy(gv.version,"1.4.5 [17/05/2024]"); /** MAGEMin version */ /* generate parameters */ strcpy(gv.buffer,"none"); @@ -288,7 +288,7 @@ typedef struct metapelite_datasets { int n_pp; int n_ss; char ox[11][20]; - char PP[17][20]; + char PP[23][20]; char SS[16][20]; int verifyPC[16]; @@ -314,11 +314,11 @@ typedef struct metapelite_datasets { metapelite_dataset metapelite_db = { 256, /* number of endmembers */ 11, /* number of oxides */ - 17, /* number of pure phases */ + 23, /* number of pure phases */ 15, /* number of solution phases */ - {"SiO2" ,"Al2O3","CaO" ,"MgO" ,"FeO" ,"K2O" ,"Na2O" ,"TiO2" ,"O" ,"MnO" ,"H2O" }, + {"SiO2" ,"Al2O3","CaO" ,"MgO" ,"FeO" ,"K2O" ,"Na2O" ,"TiO2" ,"O" ,"MnO" ,"H2O" }, {"q" ,"crst" ,"trd" ,"coe" ,"stv" ,"ky" ,"sill" ,"and" ,"ru" ,"sph" ,"O2" ,"H2O" , - "qfm" ,"qif" ,"nno" ,"hm" ,"cco" }, + "qfm" ,"qif" ,"nno" ,"hm" ,"cco" ,"aH2O" , "aO2" ,"aMgO" ,"aFeO" ,"aAl2O3" ,"aTiO2" }, {"liq" ,"fsp" ,"bi" ,"g" ,"ep" ,"ma" ,"mu" ,"opx" ,"sa" ,"cd" ,"st" ,"chl" ,"ctd" ,"sp" ,"ilmm" ,"aq17" }, {1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 }, // allow solvus? @@ -349,7 +349,7 @@ typedef struct metabasite_datasets { int n_pp; int n_ss; char ox[10][20]; - char PP[18][20]; + char PP[24][20]; char SS1[14][20]; int verifyPC1[14]; @@ -381,11 +381,11 @@ typedef struct metabasite_datasets { metabasite_dataset metabasite_db = { 256, /* number of endmembers */ 10, /* number of oxides */ - 18, /* number of pure phases */ + 24, /* number of pure phases */ 14, /* number of solution phases */ {"SiO2" ,"Al2O3","CaO" ,"MgO" ,"FeO" ,"K2O" ,"Na2O" ,"TiO2" ,"O" ,"H2O" }, {"q" ,"crst" ,"trd" ,"coe" ,"law" ,"ky" ,"sill" ,"and" ,"ru" ,"sph" ,"sph" ,"ab" ,"H2O" , - "qfm" ,"qif" ,"nno" ,"hm" ,"cco" }, + "qfm" ,"qif" ,"nno" ,"hm" ,"cco" ,"aH2O" , "aO2" ,"aMgO" ,"aFeO" ,"aAl2O3" ,"aTiO2" }, {"sp" ,"opx" ,"fsp" ,"liq" ,"mu" ,"ilm" ,"ol" ,"hb" ,"ep" ,"g" ,"chl" ,"bi" ,"dio" ,"abc" }, {1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 }, // allow solvus? @@ -421,7 +421,7 @@ typedef struct igneous_datasets { int n_pp; int n_ss; char ox[11][20]; - char PP[17][20]; + char PP[23][20]; char SS[15][20]; int verifyPC[15]; @@ -447,11 +447,11 @@ typedef struct igneous_datasets { igneous_dataset igneous_db = { 291, /* number of endmembers */ 11, /* number of oxides */ - 17, /* number of pure phases */ + 23, /* number of pure phases */ 15, /* number of solution phases */ {"SiO2" ,"Al2O3","CaO" ,"MgO" ,"FeO" ,"K2O" ,"Na2O" ,"TiO2" ,"O" ,"Cr2O3","H2O" }, {"q" ,"crst" ,"trd" ,"coe" ,"stv" ,"ky" ,"sill" ,"and" ,"ru" ,"sph" ,"O2" , - "qfm" ,"mw" ,"qif" ,"nno" ,"hm" ,"cco" }, + "qfm" ,"mw" ,"qif" ,"nno" ,"hm" ,"cco" ,"aH2O" , "aO2" ,"aMgO" ,"aFeO" ,"aAl2O3" ,"aTiO2" }, {"spn" ,"bi" ,"cd" ,"cpx" ,"ep" ,"g" ,"hb" ,"ilm" ,"liq" ,"ol" ,"opx" ,"fsp" ,"fl" ,"mu" ,"fper" }, {1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 }, // allow solvus? @@ -482,7 +482,7 @@ typedef struct ultramafic_datasets { int n_pp; int n_ss; char ox[7][20]; - char PP[15][20]; + char PP[21][20]; char SS[12][20]; int verifyPC[12]; @@ -509,12 +509,12 @@ typedef struct ultramafic_datasets { ultramafic_dataset ultramafic_db = { 256, /* number of endmembers */ 7, /* number of oxides */ - 15, /* number of pure phases */ + 21, /* number of pure phases */ 12, /* number of solution phases */ {"SiO2" ,"Al2O3","MgO" ,"FeO" ,"O" ,"H2O" ,"S" }, {"q" ,"crst" ,"trd" ,"coe" ,"stv" ,"ky" ,"sill" ,"and" ,"pyr" ,"O2" , - "qfm" ,"qif" ,"nno" ,"hm" ,"cco" }, - {"fl", "ol" ,"br" ,"ch" ,"atg" ,"g" ,"ta" ,"chl" ,"spi" ,"opx" ,"po" ,"anth" }, + "qfm" ,"qif" ,"nno" ,"hm" ,"cco" ,"aH2O" , "aO2" ,"aMgO" ,"aFeO" ,"aAl2O3" ,"aTiO2" }, + {"fl" ,"ol" ,"br" ,"ch" ,"atg" ,"g" ,"ta" ,"chl" ,"spi" ,"opx" ,"po" ,"anth" }, {1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 }, // allow solvus? {11 ,10 ,10 ,10 ,489 ,10 ,985 ,2691 ,100 ,196 ,10 ,274 }, // No. of pseudocompound @@ -838,6 +838,9 @@ global_variable global_variable_init( global_variable gv, if (strcmp( gv.ox[i], "H2O") == 0){ gv.H2O_id = i; } + else if (strcmp( gv.ox[i], "Al2O3") == 0){ + gv.Al2O3_id = i; + } else if (strcmp( gv.ox[i], "K2O") == 0){ gv.K2O_id = i; } diff --git a/src/pp_min_function.c b/src/pp_min_function.c index 19b2a63..d00f6ae 100644 --- a/src/pp_min_function.c +++ b/src/pp_min_function.c @@ -131,7 +131,7 @@ global_variable init_em_db( int EM_database, strcpy(PP_ref_db[i].Name, gv.PP_list[i]); for (int j = 0; j < gv.len_ox; j++){ - PP_ref_db[i].Comp[j] = 2.*NiO.Comp[j] -2.*Ni.Comp[j]; + PP_ref_db[i].Comp[j] = 2.0*NiO.Comp[j] -2.0*Ni.Comp[j]; } /* Calculate the number of atoms in the bulk-rock composition */ @@ -148,11 +148,263 @@ global_variable init_em_db( int EM_database, /* Calculate normalizing factor */ double factor = fbc/ape; - PP_ref_db[i].gbase = 2.*NiO.gbase - 2.*Ni.gbase + z_b.T*0.019145*gv.buffer_n; + PP_ref_db[i].gbase = 2.0*NiO.gbase - 2.0*Ni.gbase + z_b.T*0.019145*gv.buffer_n; PP_ref_db[i].factor = factor; - PP_ref_db[i].phase_shearModulus = 2.*NiO.phase_shearModulus - 2.*Ni.phase_shearModulus; + PP_ref_db[i].phase_shearModulus = 2.0*NiO.phase_shearModulus - 2.0*Ni.phase_shearModulus; gv.pp_flags[i][4] = 1; } + else if (strcmp( gv.PP_list[i], "aH2O") == 0){ + + PP_ref H2O = G_EM_function( EM_database, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "H2O", + state ); + + strcpy(PP_ref_db[i].Name, gv.PP_list[i]); + for (int j = 0; j < gv.len_ox; j++){ + PP_ref_db[i].Comp[j] = H2O.Comp[j]; + } + + /* Calculate the number of atoms in the bulk-rock composition */ + double fbc = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + fbc += z_b.bulk_rock[j]*z_b.apo[j]; + } + + /* Calculate the number of atom in the solution */ + double ape = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + ape += PP_ref_db[i].Comp[j]*z_b.apo[j]; + } + /* Calculate normalizing factor */ + double factor = fbc/ape; + if (gv.buffer_n < 1e-8){ + gv.buffer_n = 1e-8; + } + else if (gv.buffer_n >= 1.0){ + gv.buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + H2O.gbase; + PP_ref_db[i].factor = factor; + PP_ref_db[i].phase_shearModulus = H2O.phase_shearModulus; + gv.pp_flags[i][4] = 1; + } + else if (strcmp( gv.PP_list[i], "aO2") == 0){ + + PP_ref O2 = G_EM_function( EM_database, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "O2", + state ); + + strcpy(PP_ref_db[i].Name, gv.PP_list[i]); + for (int j = 0; j < gv.len_ox; j++){ + PP_ref_db[i].Comp[j] = O2.Comp[j]; + } + + /* Calculate the number of atoms in the bulk-rock composition */ + double fbc = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + fbc += z_b.bulk_rock[j]*z_b.apo[j]; + } + + /* Calculate the number of atom in the solution */ + double ape = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + ape += PP_ref_db[i].Comp[j]*z_b.apo[j]; + } + /* Calculate normalizing factor */ + double factor = fbc/ape; + if (gv.buffer_n <= 1e-60){ + gv.buffer_n = 1e-60; + } + else if (gv.buffer_n >= 1.0){ + gv.buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + O2.gbase; + PP_ref_db[i].factor = factor; + PP_ref_db[i].phase_shearModulus = O2.phase_shearModulus; + gv.pp_flags[i][4] = 1; + } + else if (strcmp( gv.PP_list[i], "aMgO") == 0){ + + PP_ref MgO = G_EM_function( EM_database, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "per", + state ); + + strcpy(PP_ref_db[i].Name, gv.PP_list[i]); + for (int j = 0; j < gv.len_ox; j++){ + PP_ref_db[i].Comp[j] = MgO.Comp[j]; + } + + /* Calculate the number of atoms in the bulk-rock composition */ + double fbc = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + fbc += z_b.bulk_rock[j]*z_b.apo[j]; + } + + /* Calculate the number of atom in the solution */ + double ape = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + ape += PP_ref_db[i].Comp[j]*z_b.apo[j]; + } + /* Calculate normalizing factor */ + double factor = fbc/ape; + if (gv.buffer_n <= 1e-8){ + gv.buffer_n = 1e-8; + } + else if (gv.buffer_n >= 1.0){ + gv.buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + MgO.gbase; + PP_ref_db[i].factor = factor; + PP_ref_db[i].phase_shearModulus = MgO.phase_shearModulus; + gv.pp_flags[i][4] = 1; + } + else if (strcmp( gv.PP_list[i], "aFeO") == 0){ + + PP_ref FeO = G_EM_function( EM_database, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "fper", + state ); + + strcpy(PP_ref_db[i].Name, gv.PP_list[i]); + for (int j = 0; j < gv.len_ox; j++){ + PP_ref_db[i].Comp[j] = FeO.Comp[j]; + } + + /* Calculate the number of atoms in the bulk-rock composition */ + double fbc = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + fbc += z_b.bulk_rock[j]*z_b.apo[j]; + } + + /* Calculate the number of atom in the solution */ + double ape = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + ape += PP_ref_db[i].Comp[j]*z_b.apo[j]; + } + /* Calculate normalizing factor */ + double factor = fbc/ape; + if (gv.buffer_n <= 1e-8){ + gv.buffer_n = 1e-8; + } + else if (gv.buffer_n >= 1.0){ + gv.buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + FeO.gbase; + PP_ref_db[i].factor = factor; + PP_ref_db[i].phase_shearModulus = FeO.phase_shearModulus; + gv.pp_flags[i][4] = 1; + } + else if (strcmp( gv.PP_list[i], "aAl2O3") == 0){ + + PP_ref Al2O3 = G_EM_function( EM_database, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "cor", + state ); + + strcpy(PP_ref_db[i].Name, gv.PP_list[i]); + for (int j = 0; j < gv.len_ox; j++){ + PP_ref_db[i].Comp[j] = Al2O3.Comp[j]; + } + + /* Calculate the number of atoms in the bulk-rock composition */ + double fbc = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + fbc += z_b.bulk_rock[j]*z_b.apo[j]; + } + + /* Calculate the number of atom in the solution */ + double ape = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + ape += PP_ref_db[i].Comp[j]*z_b.apo[j]; + } + /* Calculate normalizing factor */ + double factor = fbc/ape; + if (gv.buffer_n <= 1e-8){ + gv.buffer_n = 1e-8; + } + else if (gv.buffer_n >= 1.0){ + gv.buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + Al2O3.gbase; + PP_ref_db[i].factor = factor; + PP_ref_db[i].phase_shearModulus = Al2O3.phase_shearModulus; + gv.pp_flags[i][4] = 1; + } + else if (strcmp( gv.PP_list[i], "aTiO2") == 0){ + + PP_ref TiO2 = G_EM_function( EM_database, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "ru", + state ); + + strcpy(PP_ref_db[i].Name, gv.PP_list[i]); + for (int j = 0; j < gv.len_ox; j++){ + PP_ref_db[i].Comp[j] = TiO2.Comp[j]; + } + + /* Calculate the number of atoms in the bulk-rock composition */ + double fbc = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + fbc += z_b.bulk_rock[j]*z_b.apo[j]; + } + + /* Calculate the number of atom in the solution */ + double ape = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + ape += PP_ref_db[i].Comp[j]*z_b.apo[j]; + } + /* Calculate normalizing factor */ + double factor = fbc/ape; + if (gv.buffer_n <= 1e-8){ + gv.buffer_n = 1e-8; + } + else if (gv.buffer_n >= 1.0){ + gv.buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + TiO2.gbase; + PP_ref_db[i].factor = factor; + PP_ref_db[i].phase_shearModulus = TiO2.phase_shearModulus; + gv.pp_flags[i][4] = 1; + } else if (strcmp( gv.PP_list[i], "mw") == 0){ PP_ref mt = G_EM_function( EM_database, diff --git a/test/tests.jl b/test/tests.jl index 90bdc5f..4a58291 100644 --- a/test/tests.jl +++ b/test/tests.jl @@ -29,6 +29,45 @@ print_info(out) Finalize_MAGEMin(data) +@testset "test activity buffers" begin + # Initialize database - new way + data = Initialize_MAGEMin("mp", verbose=true, buffer="aH2O"); + test = 0 + data = use_predefined_bulk_rock(data, test); + + # Call optimization routine for given P & T & bulk_rock + P = 8.0 + T = 400.0 + out = point_wise_minimization(P,T, data, buffer_n=0.6); + @test sort(out.ph) == ["aH2O", "chl", "ep", "fsp", "mu", "mu", "q", "ru", "sp"] + Finalize_MAGEMin(data) + + # Initialize database - new way + data = Initialize_MAGEMin("mp", verbose=true, buffer="aTiO2"); + test = 0 + data = use_predefined_bulk_rock(data, test); + + # Call optimization routine for given P & T & bulk_rock + P = 8.0 + T = 400.0 + out = point_wise_minimization(P,T, data, buffer_n=0.6); + @test sort(out.ph) == ["H2O", "aTiO2", "chl", "ep", "fsp", "ilmm", "mu", "mu", "q"] + Finalize_MAGEMin(data) + + # Initialize database - new way + data = Initialize_MAGEMin("ig", verbose=true, buffer="aTiO2"); + test = 0 + data = use_predefined_bulk_rock(data, test); + + # Call optimization routine for given P & T & bulk_rock + P = 8.0 + T = 1200.0 + out = point_wise_minimization(P,T, data, buffer_n=0.1); + @test sort(out.ph) == ["aTiO2", "cpx", "fsp", "liq", "ol", "opx"] + + Finalize_MAGEMin(data) +end + @testset "test normalization" begin @@ -69,6 +108,15 @@ Finalize_MAGEMin(data) @test sum(out.bulk_S) ≈ 1.0 @test sum(out.bulk_S_wt) ≈ 1.0 + + P = 8.0 + T = 1900.0 + out = point_wise_minimization(P,T, data); + @test out.frac_M + out.frac_S + out.frac_F ≈ 1.0 + @test out.frac_M_wt + out.frac_S_wt + out.frac_F_wt ≈ 1.0 + @test sum(out.bulk_M) ≈ 1.0 + @test sum(out.bulk_M_wt) ≈ 1.0 + Finalize_MAGEMin(data) end