From 7aea2048a3fabf53061e75686e2e44072d8c7171 Mon Sep 17 00:00:00 2001 From: Nicolas Riel <50235786+NicolasRiel@users.noreply.github.com> Date: Thu, 7 Nov 2024 20:37:30 +0100 Subject: [PATCH] MAGEMin_C v1.6.1 (#59) - Added 1atom basis to output structure - Added buffer and buffer_n to output structure - Corrected reduced compositional system for metabasite database - Re-added consideration of pure-end member to the minimization, improves performance and accuracy --- Project.toml | 8 +- gen/magemin_library.jl | 38 +- julia/MAGEMin_wrappers.jl | 14 +- src/MAGEMin.c | 59 +- src/MAGEMin.h | 8 +- src/SB_database/SB_endmembers.h | 6 +- src/SB_database/SB_gem_function.c | 22 +- src/SB_database/SB_init_database.c | 56 +- src/SB_database/SB_init_database.h | 2 +- src/SB_database/SB_solution_phases.h | 21 + src/SB_database/sb_gss_function.c | 1066 ++++++++++++++++++++++++ src/SB_database/sb_gss_function.h | 22 + src/SB_database/sb_gss_init_function.c | 12 +- src/SB_database/sb_gss_init_function.h | 6 +- src/TC_database/TC_init_database.c | 14 +- src/TC_database/objective_functions.c | 3 +- src/TC_database/tc_gss_function.c | 16 +- src/all_endmembers.h | 1 + src/all_solution_phases.h | 2 +- src/dump_function.c | 14 +- src/gem_function.c | 4 +- src/initialize.c | 13 +- src/pp_min_function.c | 277 +++++- src/pp_min_function.h | 6 + src/toolkit.c | 45 + src/toolkit.h | 2 +- test/test_diagram_test0_mb.jl | 98 +-- 27 files changed, 1652 insertions(+), 183 deletions(-) create mode 100644 src/SB_database/SB_solution_phases.h create mode 100644 src/SB_database/sb_gss_function.c create mode 100644 src/SB_database/sb_gss_function.h diff --git a/Project.toml b/Project.toml index e737939..c810dfd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MAGEMin_C" uuid = "e5d170eb-415a-4524-987b-12f1bce1ddab" -authors = ["Boris Kaus & Nicolas Riel "] -version = "1.6.0" +authors = ["Nicolas Riel & Boris Kaus "] +version = "1.6.1" [deps] CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82" @@ -15,8 +15,8 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] CEnum = "0.4" -DataFrames = "1.6.1" CSV = "0.10.14" -MAGEMin_jll = "1.5.6 - 1.5.6" +DataFrames = "1.6.1" +MAGEMin_jll = "1.5.7 - 1.5.7" ProgressMeter = "1" julia = "1.6" diff --git a/gen/magemin_library.jl b/gen/magemin_library.jl index 364524d..1826116 100644 --- a/gen/magemin_library.jl +++ b/gen/magemin_library.jl @@ -131,9 +131,9 @@ function rho_wat_calc(wat, Pbar, TK, opt) end mutable struct EM_db_sb_ - Name::NTuple{20, Cchar} - FullName::NTuple{50, Cchar} - Equation::NTuple{50, Cchar} + Name::NTuple{50, Cchar} + FullName::NTuple{80, Cchar} + Equation::NTuple{90, Cchar} Comp::NTuple{6, Cdouble} input_1::NTuple{10, Cdouble} input_2::NTuple{3, Cdouble} @@ -146,6 +146,10 @@ function Access_SB_EM_DB(id, EM_dataset) ccall((:Access_SB_EM_DB, libMAGEMin), EM_db_sb, (Cint, Cint), id, EM_dataset) end +function SB_G_EM_function(EM_database, len_ox, id, bulk_rock, apo, P, T, name, state) + ccall((:SB_G_EM_function, libMAGEMin), PP_ref, (Cint, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, Ptr{Cchar}, Ptr{Cchar}), EM_database, len_ox, id, bulk_rock, apo, P, T, name, state) +end + # typedef double ( * nlopt_func ) ( unsigned n , const double * x , double * gradient , /* NULL if not needed */ void * func_data ) const nlopt_func = Ptr{Cvoid} @@ -1114,6 +1118,7 @@ const stb_PP_phase = stb_PP_phases struct stb_systems MAGEMin_ver::Ptr{Cchar} dataset::Ptr{Cchar} + database::Ptr{Cchar} bulk_res_norm::Cdouble n_iterations::Cint status::Cint @@ -1125,6 +1130,8 @@ struct stb_systems X::Cdouble bulk::Ptr{Cdouble} bulk_wt::Ptr{Cdouble} + buffer::Ptr{Cchar} + buffer_n::Cdouble gamma::Ptr{Cdouble} G::Cdouble M_sys::Cdouble @@ -1175,6 +1182,7 @@ struct stb_systems ph::Ptr{Ptr{Cchar}} ph_frac::Ptr{Cdouble} ph_frac_wt::Ptr{Cdouble} + ph_frac_1at::Ptr{Cdouble} ph_frac_vol::Ptr{Cdouble} ph_type::Ptr{Cint} ph_id::Ptr{Cint} @@ -1506,7 +1514,7 @@ mutable struct stx11_datasets n_pp::Cint n_ss::Cint ox::NTuple{6, NTuple{20, Cchar}} - PP::NTuple{10, NTuple{20, Cchar}} + PP::NTuple{9, NTuple{20, Cchar}} SS::NTuple{14, NTuple{20, Cchar}} verifyPC::NTuple{14, Cint} n_SS_PC::NTuple{14, Cint} @@ -2196,6 +2204,14 @@ function SS_mtl_pc_init_function(SS_pc_xeos, iss, name) ccall((:SS_mtl_pc_init_function, libMAGEMin), Cvoid, (Ptr{PC_ref}, Cint, Ptr{Cchar}), SS_pc_xeos, iss, name) end +function SB_SS_init(SS_init, gv) + ccall((:SB_SS_init, libMAGEMin), Cvoid, (Ptr{SS_init_type}, global_variable), SS_init, gv) +end + +function G_SS_sb11_EM_function(gv, SS_ref_db, EM_dataset, z_b, name) + ccall((:G_SS_sb11_EM_function, libMAGEMin), SS_ref, (global_variable, SS_ref, Cint, bulk_info, Ptr{Cchar}), gv, SS_ref_db, EM_dataset, z_b, name) +end + function PGE(z_b, gv, PC_read, SS_objective, NLopt_opt, splx_data, PP_ref_db, SS_ref_db, cp) ccall((:PGE, libMAGEMin), global_variable, (bulk_info, global_variable, Ptr{PC_type}, Ptr{obj_type}, Ptr{NLopt_type}, Ptr{simplex_data}, Ptr{PP_ref}, Ptr{SS_ref}, Ptr{csd_phase_set}), z_b, gv, PC_read, SS_objective, NLopt_opt, splx_data, PP_ref_db, SS_ref_db, cp) end @@ -2287,7 +2303,7 @@ function read_in_data(gv, input_data, n_points) ccall((:read_in_data, libMAGEMin), Cvoid, (global_variable, Ptr{io_data}, Cint), gv, input_data, n_points) end -mutable struct ketopt_t +struct ketopt_t ind::Cint opt::Cint arg::Ptr{Cchar} @@ -2295,14 +2311,12 @@ mutable struct ketopt_t i::Cint pos::Cint n_args::Cint - ketopt_t() = new() end -mutable struct ko_longopt_t +struct ko_longopt_t name::Ptr{Cchar} has_arg::Cint val::Cint - ko_longopt_t() = new() end function ketopt_permute(argv, j, n) @@ -2351,6 +2365,10 @@ function init_em_db(EM_database, z_b, gv, PP_ref_db) ccall((:init_em_db, libMAGEMin), global_variable, (Cint, bulk_info, global_variable, Ptr{PP_ref}), EM_database, z_b, gv, PP_ref_db) end +function init_em_db_sb(EM_database, z_b, gv, PP_ref_db) + ccall((:init_em_db_sb, libMAGEMin), global_variable, (Cint, bulk_info, global_variable, Ptr{PP_ref}), EM_database, z_b, gv, PP_ref_db) +end + function update_dG(splx_data) ccall((:update_dG, libMAGEMin), Cvoid, (Ptr{simplex_data},), splx_data) end @@ -2451,6 +2469,10 @@ function convert_system_comp(gv, sys_in, z_b) ccall((:convert_system_comp, libMAGEMin), Cvoid, (global_variable, Ptr{Cchar}, bulk_info), gv, sys_in, z_b) end +function get_tests_bulks(gv) + ccall((:get_tests_bulks, libMAGEMin), global_variable, (global_variable,), gv) +end + function get_act_sf_id(result, A, n) ccall((:get_act_sf_id, libMAGEMin), Cvoid, (Ptr{Cint}, Ptr{Cdouble}, Cint), result, A, n) end diff --git a/julia/MAGEMin_wrappers.jl b/julia/MAGEMin_wrappers.jl index 47c06f9..e4762cb 100644 --- a/julia/MAGEMin_wrappers.jl +++ b/julia/MAGEMin_wrappers.jl @@ -40,6 +40,9 @@ end struct gmin_struct{T,I} MAGEMin_ver :: String dataset :: String + database :: String + buffer :: String + buffer_n :: T G_system :: T # G of system Gamma :: Vector{T} # Gamma P_kbar :: T # Pressure in kbar @@ -97,6 +100,7 @@ struct gmin_struct{T,I} ph_frac :: Vector{T} # phase fractions ph_frac_wt :: Vector{T} # phase fractions + ph_frac_1at :: 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 @@ -183,7 +187,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", "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", "h2oL", "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"]), ss_infos("spn", 3, 2, ["none", "herc", "sp", "usp"], ["none", "x", "y"])], ["sp", "opx", "fsp", "liq", "mu", "ilm", "ol", "hb", "ep", "g", "chl", "bi", "dio", "abc", "spn"], ["q", "crst", "trd", "coe", "law", "ky", "sill", "and", "ru", "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", "h2o1L"], ["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", "x"])], ["spn", "bi", "cd", "cpx", "ep", "g", "hb", "ilm", "liq", "ol", "opx", "fsp", "fl", "mu", "fper"], ["ne", "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("igad", "Igneous alkaline dry (Weller et al., 2024)", ss_infos[ss_infos("spn", 8, 7, ["none", "nsp", "isp", "nhc", "ihc", "nmt", "imt", "pcr", "usp"], ["none", "x", "y", "c", "t", "Q1", "Q2", "Q3"]), 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("g", 6, 5, ["none", "py", "alm", "gr", "andr", "knr", "tig"], ["none", "x", "c", "f", "cr", "t"]), ss_infos("ilm", 5, 4, ["none", "oilm", "dilm", "hm", "ogk", "dgk"], ["none", "i", "m", "Q", "Qt"]), ss_infos("liq", 14, 13, ["none", "q3L", "sl1L", "wo1L", "fo2L", "fa2L", "nmL", "hmL", "ekL", "tiL", "kmL", "anL", "ab1L", "enL", "kfL"], ["none", "wo", "sl", "fo", "fa", "ns", "hm", "ek", "ti", "ks", "yan", "yab", "yen", "ykf"]), ss_infos("ol", 4, 3, ["none", "mnt", "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("lct", 2, 1, ["none", "nlc", "klc"], ["none", "n"]), ss_infos("mel", 5, 4, ["none", "geh", "ak", "fak", "nml", "fge"], ["none", "x", "n", "y", "f"]), ss_infos("ness", 6, 5, ["none", "neN", "neS", "neK", "neO", "neC", "neF"], ["none", "s", "k", "Q", "f", "c"]), ss_infos("kals", 2, 1, ["none", "nks", "kls"], ["none", "k"])], ["spn", "cpx", "g", "ilm", "liq", "ol", "opx", "fsp", "lct", "mel", "ness", "kals"], ["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"]), db_infos("ume", "Ultramafic extended (Evans & Frost., 2021) with pl, hb and aug from Green et al., 2016", 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"]), ss_infos("pl4tr", 2, 1, ["none", "ab", "an"], ["none", "ca"]), ss_infos("hb", 9, 8, ["none", "tr", "tsm", "prgm", "glm", "cumm", "grnm", "a", "b", "mrb"], ["none", "x", "y", "z", "a", "c", "f", "Q1", "Q2"]), ss_infos("aug", 8, 7, ["none", "di", "cenh", "cfs", "jdm", "acmm", "ocats", "dcats", "fmc"], ["none", "x", "y", "f", "z", "j", "Qfm", "Qa1"])], ["fl", "ol", "br", "ch", "atg", "g", "ta", "chl", "spi", "opx", "po", "anth", "pl4tr", "hb", "aug"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and", "pyr", "O2", "qfm", "qif", "nno", "hm", "cco", "aH2O", "aO2", "aMgO", "aFeO", "aAl2O3", "aTiO2"]), db_infos("mtl", "Mantle (Holland et al., 2013)", ss_infos[ss_infos("g", 6, 5, ["none", "py", "alm", "gr", "maj", "gfm", "nagt"], ["none", "x", "c", "y", "Q", "n"]), ss_infos("fp", 2, 1, ["none", "per", "fper"], ["none", "x"]), ss_infos("mpv", 5, 4, ["none", "mpv", "fpvm", "cpvm", "apv", "npvm"], ["none", "x", "y", "c", "n"]), ss_infos("cpv", 5, 4, ["none", "mpv", "fpvm", "cpvm", "apv", "npvm"], ["none", "x", "y", "c", "n"]), ss_infos("crn", 3, 2, ["none", "cor", "mcor", "fcor"], ["none", "x", "y"]), ss_infos("cf", 6, 5, ["none", "macf", "cacf", "mscf", "fscf", "oscf", "nacfm"], ["none", "y", "x", "Q", "c", "n"]), ss_infos("nal", 7, 6, ["none", "nanal", "canal", "manal", "msnal", "fsnal", "o1nal", "o2nal"], ["none", "y", "x", "Q1", "Q2", "c", "n"]), ss_infos("aki", 3, 2, ["none", "aak", "mak", "fak"], ["none", "x", "y"]), ss_infos("ol", 2, 1, ["none", "fo", "fa"], ["none", "x"]), ss_infos("wad", 2, 1, ["none", "mwd", "fwd"], ["none", "x"]), ss_infos("ring", 2, 1, ["none", "mrw", "frw"], ["none", "x"]), ss_infos("cpx", 6, 5, ["none", "di", "cfs", "cats", "jd", "cen", "cfm"], ["none", "x", "y", "o", "n", "Q"]), ss_infos("opx", 5, 4, ["none", "en", "fs", "fm", "odi", "mgts"], ["none", "x", "y", "c", "Q"]), ss_infos("hpx", 5, 4, ["none", "en", "fs", "fm", "odi", "hmts"], ["none", "x", "y", "c", "Q"])], ["g", "fp", "mpv", "cpv", "crn", "cf", "nal", "aki", "ol", "wad", "ring", "cpx", "opx", "hpx"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and"])]; + 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", "h2oL", "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("ilmm", 4, 3, ["none", "oilm", "dilm", "dhem", "geik"], ["none", "c", "t", "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"]), ss_infos("spn", 3, 2, ["none", "herc", "sp", "usp"], ["none", "x", "y"])], ["sp", "opx", "fsp", "liq", "mu", "ilmm", "ol", "hb", "ep", "g", "chl", "bi", "dio", "abc", "spn"], ["q", "crst", "trd", "coe", "law", "ky", "sill", "and", "ru", "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", "h2o1L"], ["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", "x"])], ["spn", "bi", "cd", "cpx", "ep", "g", "hb", "ilm", "liq", "ol", "opx", "fsp", "fl", "mu", "fper"], ["ne", "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("igad", "Igneous alkaline dry (Weller et al., 2024)", ss_infos[ss_infos("spn", 8, 7, ["none", "nsp", "isp", "nhc", "ihc", "nmt", "imt", "pcr", "usp"], ["none", "x", "y", "c", "t", "Q1", "Q2", "Q3"]), 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("g", 6, 5, ["none", "py", "alm", "gr", "andr", "knr", "tig"], ["none", "x", "c", "f", "cr", "t"]), ss_infos("ilm", 5, 4, ["none", "oilm", "dilm", "hm", "ogk", "dgk"], ["none", "i", "m", "Q", "Qt"]), ss_infos("liq", 14, 13, ["none", "q3L", "sl1L", "wo1L", "fo2L", "fa2L", "nmL", "hmL", "ekL", "tiL", "kmL", "anL", "ab1L", "enL", "kfL"], ["none", "wo", "sl", "fo", "fa", "ns", "hm", "ek", "ti", "ks", "yan", "yab", "yen", "ykf"]), ss_infos("ol", 4, 3, ["none", "mnt", "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("lct", 2, 1, ["none", "nlc", "klc"], ["none", "n"]), ss_infos("mel", 5, 4, ["none", "geh", "ak", "fak", "nml", "fge"], ["none", "x", "n", "y", "f"]), ss_infos("ness", 6, 5, ["none", "neN", "neS", "neK", "neO", "neC", "neF"], ["none", "s", "k", "Q", "f", "c"]), ss_infos("kals", 2, 1, ["none", "nks", "kls"], ["none", "k"])], ["spn", "cpx", "g", "ilm", "liq", "ol", "opx", "fsp", "lct", "mel", "ness", "kals"], ["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"]), db_infos("ume", "Ultramafic extended (Evans & Frost., 2021) with pl, hb and aug from Green et al., 2016", 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"]), ss_infos("pl4tr", 2, 1, ["none", "ab", "an"], ["none", "ca"]), ss_infos("hb", 9, 8, ["none", "tr", "tsm", "prgm", "glm", "cumm", "grnm", "a", "b", "mrb"], ["none", "x", "y", "z", "a", "c", "f", "Q1", "Q2"]), ss_infos("aug", 8, 7, ["none", "di", "cenh", "cfs", "jdm", "acmm", "ocats", "dcats", "fmc"], ["none", "x", "y", "f", "z", "j", "Qfm", "Qa1"])], ["fl", "ol", "br", "ch", "atg", "g", "ta", "chl", "spi", "opx", "po", "anth", "pl4tr", "hb", "aug"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and", "pyr", "O2", "qfm", "qif", "nno", "hm", "cco", "aH2O", "aO2", "aMgO", "aFeO", "aAl2O3", "aTiO2"]), db_infos("mtl", "Mantle (Holland et al., 2013)", ss_infos[ss_infos("g", 6, 5, ["none", "py", "alm", "gr", "maj", "gfm", "nagt"], ["none", "x", "c", "y", "Q", "n"]), ss_infos("fp", 2, 1, ["none", "per", "fper"], ["none", "x"]), ss_infos("mpv", 5, 4, ["none", "mpv", "fpvm", "cpvm", "apv", "npvm"], ["none", "x", "y", "c", "n"]), ss_infos("cpv", 5, 4, ["none", "mpv", "fpvm", "cpvm", "apv", "npvm"], ["none", "x", "y", "c", "n"]), ss_infos("crn", 3, 2, ["none", "cor", "mcor", "fcor"], ["none", "x", "y"]), ss_infos("cf", 6, 5, ["none", "macf", "cacf", "mscf", "fscf", "oscf", "nacfm"], ["none", "y", "x", "Q", "c", "n"]), ss_infos("nal", 7, 6, ["none", "nanal", "canal", "manal", "msnal", "fsnal", "o1nal", "o2nal"], ["none", "y", "x", "Q1", "Q2", "c", "n"]), ss_infos("aki", 3, 2, ["none", "aak", "mak", "fak"], ["none", "x", "y"]), ss_infos("ol", 2, 1, ["none", "fo", "fa"], ["none", "x"]), ss_infos("wad", 2, 1, ["none", "mwd", "fwd"], ["none", "x"]), ss_infos("ring", 2, 1, ["none", "mrw", "frw"], ["none", "x"]), ss_infos("cpx", 6, 5, ["none", "di", "cfs", "cats", "jd", "cen", "cfm"], ["none", "x", "y", "o", "n", "Q"]), ss_infos("opx", 5, 4, ["none", "en", "fs", "fm", "odi", "mgts"], ["none", "x", "y", "c", "Q"]), ss_infos("hpx", 5, 4, ["none", "en", "fs", "fm", "odi", "hmts"], ["none", "x", "y", "c", "Q"])], ["g", "fp", "mpv", "cpv", "crn", "cf", "nal", "aki", "ol", "wad", "ring", "cpx", "opx", "hpx"], ["q", "crst", "trd", "coe", "stv", "ky", "sill", "and"])]; dbs = ["mp","mb","ig","igad","um","ume","mtl"] id = findall(dbs .== dtb)[1] @@ -1212,6 +1216,9 @@ function create_gmin_struct(DB, gv, time) MAGEMin_ver = unsafe_string(stb.MAGEMin_ver) dataset = unsafe_string(stb.dataset) + database = unsafe_string(stb.database) + buffer = unsafe_string(stb.buffer) + buffer_n = stb.buffer_n G_system = stb.G Gamma = unsafe_wrap(Vector{Cdouble},stb.gamma,gv.len_ox) P_kbar = stb.P @@ -1275,6 +1282,7 @@ function create_gmin_struct(DB, gv, time) ph_frac = unsafe_wrap(Vector{Cdouble},stb.ph_frac, n_ph) ph_frac_wt = unsafe_wrap(Vector{Cdouble},stb.ph_frac_wt, n_ph) + ph_frac_1at = unsafe_wrap(Vector{Cdouble},stb.ph_frac_1at, n_ph) ph_frac_vol = unsafe_wrap(Vector{Cdouble},stb.ph_frac_vol, n_ph) ph_type = unsafe_wrap(Vector{Cint}, stb.ph_type, n_ph) ph_id = unsafe_wrap(Vector{Cint}, stb.ph_id , n_ph) @@ -1299,7 +1307,7 @@ function create_gmin_struct(DB, gv, time) time_ms = time*1000.0 # Store all in output struct - out = gmin_struct{Float64,Int64}( MAGEMin_ver, dataset, G_system, Gamma, P_kbar, T_C, X, M_sys, + out = gmin_struct{Float64,Int64}( MAGEMin_ver, dataset, database, buffer, buffer_n, G_system, Gamma, P_kbar, T_C, X, M_sys, bulk, bulk_M, bulk_S, bulk_F, bulk_wt, bulk_M_wt, bulk_S_wt, bulk_F_wt, frac_M, frac_S, frac_F, @@ -1308,7 +1316,7 @@ function create_gmin_struct(DB, gv, time) rho, rho_M, rho_S, rho_F, fO2, dQFM, aH2O, aSiO2, aTiO2, aAl2O3, aMgO, aFeO, n_PP, n_SS, n_mSS, - ph_frac, ph_frac_wt, ph_frac_vol, ph_type, ph_id, ph, + ph_frac, ph_frac_wt, ph_frac_1at, ph_frac_vol, ph_type, ph_id, ph, SS_vec, mSS_vec, PP_vec, oxides, elements, stb.Vp, stb.Vs, stb.Vp_S, stb.Vs_S, stb.bulkMod, stb.shearMod, stb.bulkModulus_M, stb.bulkModulus_S, stb.shearModulus_S, diff --git a/src/MAGEMin.c b/src/MAGEMin.c index 87c8755..4633ddb 100644 --- a/src/MAGEMin.c +++ b/src/MAGEMin.c @@ -117,20 +117,17 @@ int runMAGEMin( int argc, argc, argv ); - /* initialize global structure to store shared variables (e.g. Gamma, SS and PP list, ...) */ gv = global_variable_init( gv, &z_b ); - /* Allocate both pure and solid-solution databases */ DB = InitializeDatabases( gv, gv.EM_database ); - /* initialize simplex (levelling stage using pseudocompounds) @@ -157,31 +154,7 @@ int runMAGEMin( int argc, /* get bulk rock composition parsed from args */ - - if (gv.EM_database == 0){ - gv = get_bulk_metapelite( gv ); - } - else if (gv.EM_database == 1){ - gv = get_bulk_metabasite( gv ); - } - else if (gv.EM_database == 2){ - gv = get_bulk_igneous( gv ); - } - else if (gv.EM_database == 3){ - gv = get_bulk_igneous_igad( gv ); - } - else if (gv.EM_database == 4){ - gv = get_bulk_ultramafic( gv ); - } - else if (gv.EM_database == 5){ - gv = get_bulk_ultramafic_ext( gv ); - } - else if (gv.EM_database == 6){ - gv = get_bulk_mantle( gv ); - } - else{ - printf(" Wrong database...\n"); - } + gv = get_tests_bulks( gv ); /****************************************************************************************/ /** LAUNCH MINIMIZATION ROUTINE **/ @@ -352,10 +325,18 @@ int runMAGEMin( int argc, SS_ref *SS_ref_db ){ /* initialize endmember database for given P-T point */ - gv = init_em_db( EM_database, - z_b, /** bulk rock informations */ - gv, /** global variables (e.g. Gamma) */ - PP_ref_db ); + if ( strcmp(gv.research_group, "tc") == 0 ){ + gv = init_em_db( EM_database, + z_b, /** bulk rock informations */ + gv, /** global variables (e.g. Gamma) */ + PP_ref_db ); + } + else if ( strcmp(gv.research_group, "sb") == 0 ){ + gv = init_em_db_sb( EM_database, + z_b, /** bulk rock informations */ + gv, /** global variables (e.g. Gamma) */ + PP_ref_db ); + } /* Calculate solution phase data at given P-T conditions (G0 based on G0 of endmembers) */ gv = init_ss_db( EM_database, @@ -928,8 +909,15 @@ Databases InitializeDatabases( global_variable gv, SS_init_type SS_init[gv.len_ss]; - TC_SS_init( SS_init, - gv ); + + if ( strcmp(gv.research_group, "tc") == 0 ){ + TC_SS_init( SS_init, + gv ); + } + else if (strcmp(gv.research_group, "sb") == 0 ){ + SB_SS_init( SS_init, + gv ); + } DB.SS_ref_db = malloc ((gv.len_ss) * sizeof(SS_ref)); for (int iss = 0; iss < gv.len_ss; iss++){ @@ -1078,12 +1066,15 @@ void FreeDatabases( global_variable gv, free(DB.sp[0].bulk_F_wt); free(DB.sp[0].ph_frac); free(DB.sp[0].ph_frac_wt); + free(DB.sp[0].ph_frac_1at); free(DB.sp[0].ph_frac_vol); free(DB.sp[0].ph_id); free(DB.sp[0].ph_type); free(DB.sp[0].MAGEMin_ver); + free(DB.sp[0].buffer); free(DB.sp[0].dataset); + free(DB.sp[0].database); /* ==================== CP ============================== */ diff --git a/src/MAGEMin.h b/src/MAGEMin.h index c726ae6..259d879 100644 --- a/src/MAGEMin.h +++ b/src/MAGEMin.h @@ -702,7 +702,8 @@ typedef struct stb_systems { char *MAGEMin_ver; char *dataset; - + char *database; + double bulk_res_norm; int n_iterations; int status; @@ -716,7 +717,9 @@ typedef struct stb_systems { double X; double *bulk; double *bulk_wt; - + char *buffer; + double buffer_n; + double *gamma; double G; double M_sys; @@ -766,6 +769,7 @@ typedef struct stb_systems { char **ph; /* phases names */ double *ph_frac; /* phase fractions */ double *ph_frac_wt; /* phase fractions in wt fraction */ + double *ph_frac_1at; /* phase fractions in wt fraction */ double *ph_frac_vol; /* phase fractions in wt fraction */ int *ph_type; /* 0 -> Solution phases; 1 -> Pure phases */ int *ph_id; /* position in the related stb_SS_phase or stb_PP_phase structure arrays */ diff --git a/src/SB_database/SB_endmembers.h b/src/SB_database/SB_endmembers.h index 3a7d91c..3a45e07 100644 --- a/src/SB_database/SB_endmembers.h +++ b/src/SB_database/SB_endmembers.h @@ -13,9 +13,9 @@ /** store endmember database **/ typedef struct EM_db_sb_ { - char Name[20]; /** pure species name */ - char FullName[50]; /** pure species name */ - char Equation[50]; /** pure species name */ + char Name[50]; /** pure species name */ + char FullName[80]; /** pure species name */ + char Equation[90]; /** pure species name */ double Comp[6]; /** pure species composition [0-10] + number of atom [11] */ double input_1[10]; /** second line of the thermodynamics datable */ double input_2[3]; /** second line of the thermodynamics datable */ diff --git a/src/SB_database/SB_gem_function.c b/src/SB_database/SB_gem_function.c index e40849c..ada06be 100644 --- a/src/SB_database/SB_gem_function.c +++ b/src/SB_database/SB_gem_function.c @@ -22,7 +22,7 @@ #include "../toolkit.h" #include "SB_gem_function.h" -#define eps 1e-12 +#define eps 1e-10 double plg(double t) { @@ -49,6 +49,7 @@ double plg(double t) { i += 1; } + return plg; } @@ -73,7 +74,7 @@ PP_ref SB_G_EM_function( int EM_dataset, for (i = 0; i < len_ox; i ++){ composition[i] = EM_return.Comp[id[i]]; } - + double kbar2bar = 1e3; double RTlnf = 0.0; double T0 = 298.15; @@ -93,7 +94,7 @@ PP_ref SB_G_EM_function( int EM_dataset, double F0 = EM_return.input_1[0]; double n = EM_return.input_1[1]; - double V0 = EM_return.input_1[2]; + double V0 = -EM_return.input_1[2]; double K0 = EM_return.input_1[3]; double Kp = EM_return.input_1[4]; double z00 = EM_return.input_1[5]; @@ -140,7 +141,7 @@ PP_ref SB_G_EM_function( int EM_dataset, bad = 1; itic = 0; - while (itic < 100){ + while (itic < max_ite){ itic += 1; @@ -206,24 +207,24 @@ PP_ref SB_G_EM_function( int EM_dataset, V -= dv; - if (itic > max_ite || fabs(f1) > 1e40){ + if (itic >= max_ite || fabs(f1) > 1e40){ if (fabs(f1 / P) < 0.0){ ibad = 5; printf("ERROR abs(f1/p)\n"); } - else if( fabs(dv / (1.0 + V)) < eps ){ + } + else if( fabs(dv / (1.0 + V)) < eps ){ bad = 0; break; - } } + } if (bad == 1){ printf("ERROR bad\n"); } - /* get helmoltz energy:*/ f = 0.5 * pow((V0 / V),r23) - 0.5; z = 1.0 + (aii + aiikk2 * f) * f; @@ -263,12 +264,13 @@ PP_ref SB_G_EM_function( int EM_dataset, for (i = 0; i < len_ox; i++){ PP_ref_db.Comp[i] = composition[i]; } - PP_ref_db.gbase = gbase; + PP_ref_db.gbase = gbase/kbar2bar; PP_ref_db.factor = factor; PP_ref_db.phase_shearModulus = (EM_return.input_2[0]*kbar2bar + (P - P0)*(EM_return.input_2[1])*kbar2bar + (T - T0)*(EM_return.input_2[2]))/kbar2bar; - // printf(" %4s %+10f\n",name,PP_ref_db.gbase); + // printf("gbase %4s %+10f\n",name,PP_ref_db.gbase); + // printf("phase_shearModulus %4s %+10f\n",name,PP_ref_db.phase_shearModulus); // for (i = 0; i < len_ox; i++){ // printf("%+10f",PP_ref_db.Comp[i]*PP_ref_db.factor); // } diff --git a/src/SB_database/SB_init_database.c b/src/SB_database/SB_init_database.c index d11974f..e5e56b4 100644 --- a/src/SB_database/SB_init_database.c +++ b/src/SB_database/SB_init_database.c @@ -22,23 +22,23 @@ oxide_data oxide_info_sb = { 15, /* number of endmembers */ - {"SiO2" ,"Al2O3","CaO" ,"MgO" ,"FeO" ,"K2O" ,"Na2O" ,"TiO2" ,"O" ,"MnO" ,"Cr2O3","H2O" ,"CO2" ,"S" ,"Cl" }, - {"Si" ,"Al" ,"Ca" ,"Mg" ,"Fe" ,"K" ,"Na" ,"Ti" ,"O" ,"Mn" ,"Cr" ,"H" ,"C" ,"S" ,"Cl" }, - {60.08 ,101.96 ,56.08 ,40.30 ,71.85 ,94.2 ,61.98 ,79.88 ,16.0 ,70.94 ,151.99 ,18.015 ,44.01 , 32.06 ,35.453 }, - {3.0 ,5.0 ,2.0 ,2.0 ,2.0 ,3.0 ,3.0 ,3.0 ,1.0 ,2.0 ,5.0 ,3.0 ,3.0 , 1.0 ,1.0 }, - {66.7736,108.653,42.9947,40.3262,38.7162,69.1514,61.1729,70.3246,30.5827,40.1891,106.9795,69.5449,62.8768,9.5557,33.2556 }, - {2,3,1,1,1,1,1,2,1,1,3,1,2,0,0}, - {1,2,1,1,1,2,2,1,1,1,2,2,1,1,1} + {"SiO2" ,"CaO" ,"Al2O3","FeO" ,"MgO" ,"Na2O" ,"K2O" ,"TiO2" ,"O" ,"MnO" ,"Cr2O3","H2O" ,"CO2" ,"S" ,"Cl"}, + {"Si" ,"Ca" ,"Al" ,"Fe" ,"Mg" ,"Na" ,"K" ,"Ti" ,"O" ,"Mn" ,"Cr" ,"H" ,"C" ,"S" ,"Cl"}, + {60.08 ,56.08 ,101.96 ,71.85 ,40.30 ,61.98 ,94.2 ,79.88 ,16.0 ,70.94 ,151.99 ,18.015 ,44.01 , 32.06 ,35.453}, + {3.0 ,2.0 ,5.0 ,2.0 ,2.0 ,3.0 ,3.0 ,3.0 ,1.0 ,2.0 ,5.0 ,3.0 ,3.0 , 1.0 ,1.0}, + {66.7736,42.9947,108.653,38.7162,40.3262,61.1729,69.1514,70.3246,30.5827,40.1891,106.9795,69.5449,62.8768,9.5557,33.2556}, + {2,1,3,1,1,1,1,2,1,1,3,1,2,0,0}, + {1,1,2,1,1,2,2,1,1,1,2,2,1,1,1} // for the standard molar entropy the values are already normalized by the reference temperature = 298.15K (25°C) and expressed in kJ }; // SiO2:0 CaO:1 Al2O3:2 FeO:3 MgO:4 Na2O:5 */ stx11_dataset stx11_db = { - 2011, /* Endmember default dataset number */ + 2011, /* Endmember default dataset number */ 6, /* number of oxides */ - 10, /* number of pure phases */ + 9, /* number of pure phases */ 14, /* number of solution phases */ {"SiO2" ,"CaO" ,"Al2O3","FeO" ,"MgO" ,"Na2O" }, - {"neph" ,"ky" ,"seif" ,"st" ,"coe" ,"qtz" ,"capv" ,"aMgO" ,"aFeO" ,"aAl2O3" }, + {"neph" ,"ky" ,"st" ,"coe" ,"qtz" ,"capv" ,"aMgO" ,"aFeO" ,"aAl2O3" }, {"plg" ,"sp" ,"ol" ,"wa" ,"ri" ,"opx" ,"cpx" ,"hpcpx","ak" ,"gtmj" ,"pv" ,"ppv" ,"mw" ,"cf" }, {1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 }, // allow solvus? @@ -66,6 +66,7 @@ global_variable global_variable_SB_init( global_variable gv, if (gv.EM_database == 0){ stx11_dataset db = stx11_db; + gv.EM_dataset = db.ds_version; gv.len_pp = db.n_pp; gv.len_ss = db.n_ss; gv.len_ox = db.n_ox; @@ -107,7 +108,6 @@ global_variable global_variable_SB_init( global_variable gv, gv.SS_PC_stp[i] = db.SS_PC_stp[i]; } } - /** ALLOCATE MEMORY OF OTHER GLOBAL VARIABLES */ @@ -206,10 +206,19 @@ global_variable global_variable_SB_init( global_variable gv, z_b->cpo = malloc (gv.len_ox * sizeof (double) ); z_b->ElEntropy = malloc (gv.len_ox * sizeof (double) ); z_b->id = malloc (gv.len_ox * sizeof (int) ); - + z_b->elName = malloc (gv.len_ox * sizeof (char*) ); + for (i = 0; i < (gv.len_ox); i++){ + z_b->elName[i] = malloc(20 * sizeof(char)); + } /** retrieve the right set of oxide and their informations */ + gv.Al2O3_id = -1; + gv.CaO_id = -1; + gv.Na2O_id = -1; + gv.FeO_id = -1; + gv.MgO_id = -1; + oxide_data ox_in = oxide_info_sb; for (i = 0; i < gv.len_ox; i++){ for (j = 0; j < ox_in.n_ox; j++){ @@ -217,12 +226,24 @@ global_variable global_variable_SB_init( global_variable gv, if (strcmp( gv.ox[i], "Al2O3") == 0){ gv.Al2O3_id = i; } - else if (strcmp( gv.ox[i], "TiO2") == 0){ - gv.TiO2_id = i; + // else if (strcmp( gv.ox[i], "TiO2") == 0){ + // gv.TiO2_id = i; + // } + // else if (strcmp( gv.ox[i], "O") == 0){ + // gv.O_id = i; + // } + else if (strcmp( gv.ox[i], "CaO") == 0){ + gv.CaO_id = i; + } + else if (strcmp( gv.ox[i], "Na2O") == 0){ + gv.Na2O_id = i; } - else if (strcmp( gv.ox[i], "O") == 0){ - gv.O_id = i; - } + else if (strcmp( gv.ox[i], "MgO") == 0){ + gv.MgO_id = i; + } + else if (strcmp( gv.ox[i], "FeO") == 0){ + gv.FeO_id = i; + } z_b->apo[i] = ox_in.atPerOx[j]; z_b->masspo[i] = ox_in.oxMass[j]; z_b->opo[i] = ox_in.OPerOx[j]; @@ -267,6 +288,7 @@ global_variable get_bulk_stx11( global_variable gv) { printf(" - No predefined bulk provided -> user custom bulk (if none provided, will run default KLB1)\n"); } } + if (gv.test == 0){ //KLB1 /* SiO2 Al2O3 CaO MgO FeO K2O Na2O TiO2 O Cr2O3 H2O */ /* Bulk rock composition of Peridotite from Holland et al., 2018, given by E. Green */ diff --git a/src/SB_database/SB_init_database.h b/src/SB_database/SB_init_database.h index a7fea06..b6905b9 100644 --- a/src/SB_database/SB_init_database.h +++ b/src/SB_database/SB_init_database.h @@ -35,7 +35,7 @@ int n_pp; int n_ss; char ox[6][20]; - char PP[10][20]; + char PP[9][20]; char SS[14][20]; int verifyPC[14]; diff --git a/src/SB_database/SB_solution_phases.h b/src/SB_database/SB_solution_phases.h new file mode 100644 index 0000000..e5cf88d --- /dev/null +++ b/src/SB_database/SB_solution_phases.h @@ -0,0 +1,21 @@ +/*@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ** + ** Project : MAGEMin + ** License : GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 + ** Developers : Nicolas Riel, Boris Kaus + ** Contributors : Dominguez, H., Assunção J., Green E., Berlie N., and Rummel L. + ** Organization : Institute of Geosciences, Johannes-Gutenberg University, Mainz + ** Contact : nriel[at]uni-mainz.de, kaus[at]uni-mainz.de + ** + ** ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @*/ +#ifndef __SB_sol_H_ +#define __SB_sol_H_ + + #include "sb_gss_init_function.h" + #include "sb_gss_function.h" + // #include "objective_functions.h" + // #include "NLopt_opt_function.h" + + /* include pseudocompounds */ + // #include "SS_xeos_PC_mtl.h" +#endif \ No newline at end of file diff --git a/src/SB_database/sb_gss_function.c b/src/SB_database/sb_gss_function.c new file mode 100644 index 0000000..92863e6 --- /dev/null +++ b/src/SB_database/sb_gss_function.c @@ -0,0 +1,1066 @@ +/*@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ** + ** Project : MAGEMin + ** License : GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 + ** Developers : Nicolas Riel, Boris Kaus + ** Contributors : Dominguez, H., Assunção J., Green E., Berlie N., and Rummel L. + ** Organization : Institute of Geosciences, Johannes-Gutenberg University, Mainz + ** Contact : nriel[at]uni-mainz.de, kaus[at]uni-mainz.de + ** + ** ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @*/ +/** +Function to calculate the reference chemical potential of solid-solutions + +Stixrude thermodynamic database for the mantle minerals +*/ + +#include +#include +#include +#include +#include +#include + +#include "../MAGEMin.h" +#include "../initialize.h" +#include "../all_solution_phases.h" +#include "../simplex_levelling.h" +#include "../toolkit.h" + +/** + Solution phase data for sb11_plg +*/ +SS_ref G_SS_sb11_plg_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"an","ab",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 26000.0; + + + em_data an = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "an", + "equilibrium"); + + em_data ab = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "ab", + "equilibrium"); + + SS_ref_db.gbase[0] = an.gb; + SS_ref_db.gbase[1] = ab.gb; + + SS_ref_db.ElShearMod[0] = an.ElShearMod; + SS_ref_db.ElShearMod[1] = ab.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = an.C[i]; + SS_ref_db.Comp[1][i] = ab.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_sp +*/ +SS_ref G_SS_sb11_sp_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"sp","hc",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 5000.0; + + + em_data sp = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "sp", + "equilibrium"); + + em_data hc = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "hc", + "equilibrium"); + + SS_ref_db.gbase[0] = sp.gb; + SS_ref_db.gbase[1] = hc.gb; + + SS_ref_db.ElShearMod[0] = sp.ElShearMod; + SS_ref_db.ElShearMod[1] = hc.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = sp.C[i]; + SS_ref_db.Comp[1][i] = hc.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_ol +*/ +SS_ref G_SS_sb11_ol_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"fa","fo",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 7600.0; + + + em_data fa = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "fa", + "equilibrium"); + + em_data fo = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "fo", + "equilibrium"); + + SS_ref_db.gbase[0] = fa.gb; + SS_ref_db.gbase[1] = fo.gb; + + SS_ref_db.ElShearMod[0] = fa.ElShearMod; + SS_ref_db.ElShearMod[1] = fo.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = fa.C[i]; + SS_ref_db.Comp[1][i] = fo.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_wa +*/ +SS_ref G_SS_sb11_wa_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"fewa","mgwa",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 16500.0; + + + em_data fewa = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "fewa", + "equilibrium"); + + em_data mgwa = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "mgwa", + "equilibrium"); + + SS_ref_db.gbase[0] = fewa.gb; + SS_ref_db.gbase[1] = mgwa.gb; + + SS_ref_db.ElShearMod[0] = fewa.ElShearMod; + SS_ref_db.ElShearMod[1] = mgwa.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = fewa.C[i]; + SS_ref_db.Comp[1][i] = mgwa.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_ri +*/ +SS_ref G_SS_sb11_ri_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"mgri","feri",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 9100.0; + + + em_data mgri = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "mgri", + "equilibrium"); + + em_data feri = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "feri", + "equilibrium"); + + SS_ref_db.gbase[0] = mgri.gb; + SS_ref_db.gbase[1] = feri.gb; + + SS_ref_db.ElShearMod[0] = mgri.ElShearMod; + SS_ref_db.ElShearMod[1] = feri.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = mgri.C[i]; + SS_ref_db.Comp[1][i] = feri.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_opx +*/ +SS_ref G_SS_sb11_opx_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"mgts","fs","en","odi",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 0.0; + SS_ref_db.W[1] = 0.0; + SS_ref_db.W[2] = 32100.0; + SS_ref_db.W[3] = 0.0; + SS_ref_db.W[4] = 0.0; + SS_ref_db.W[5] = 48000.0; + + + em_data mgts = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "mgts", + "equilibrium"); + + em_data fs = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "fs", + "equilibrium"); + + em_data en = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "en", + "equilibrium"); + + em_data odi = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "odi", + "equilibrium"); + + SS_ref_db.gbase[0] = mgts.gb; + SS_ref_db.gbase[1] = fs.gb; + SS_ref_db.gbase[2] = en.gb; + SS_ref_db.gbase[3] = odi.gb; + + SS_ref_db.ElShearMod[0] = mgts.ElShearMod; + SS_ref_db.ElShearMod[1] = fs.ElShearMod; + SS_ref_db.ElShearMod[2] = en.ElShearMod; + SS_ref_db.ElShearMod[3] = odi.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = mgts.C[i]; + SS_ref_db.Comp[1][i] = fs.C[i]; + SS_ref_db.Comp[2][i] = en.C[i]; + SS_ref_db.Comp[3][i] = odi.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + SS_ref_db.bounds_ref[2][0] = 0.0+eps; SS_ref_db.bounds_ref[2][1] = 1.0-eps; + SS_ref_db.bounds_ref[3][0] = 0.0+eps; SS_ref_db.bounds_ref[3][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_cpx +*/ +SS_ref G_SS_sb11_cpx_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"he","jd","cen","cats","di",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 24700.0; + SS_ref_db.W[1] = 0.0; + SS_ref_db.W[2] = 26000.0; + SS_ref_db.W[3] = 24700.0; + SS_ref_db.W[4] = 24300.0; + SS_ref_db.W[5] = 0.0; + SS_ref_db.W[6] = 60600.0; + SS_ref_db.W[7] = 10000.0; + SS_ref_db.W[8] = 0.0; + SS_ref_db.W[9] = 0.0; + + SS_ref_db.v[0] = 1.0; + SS_ref_db.v[1] = 1.0; + SS_ref_db.v[2] = 1.0; + SS_ref_db.v[3] = 3.5; + SS_ref_db.v[4] = 1.0; + + em_data he = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "he", + "equilibrium"); + + em_data jd = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "jd", + "equilibrium"); + + em_data cen = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "cen", + "equilibrium"); + + em_data cats = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "cats", + "equilibrium"); + + em_data di = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "di", + "equilibrium"); + + SS_ref_db.gbase[0] = he.gb; + SS_ref_db.gbase[1] = jd.gb; + SS_ref_db.gbase[2] = cen.gb; + SS_ref_db.gbase[3] = cats.gb; + SS_ref_db.gbase[4] = di.gb; + + SS_ref_db.ElShearMod[0] = he.ElShearMod; + SS_ref_db.ElShearMod[1] = jd.ElShearMod; + SS_ref_db.ElShearMod[2] = cen.ElShearMod; + SS_ref_db.ElShearMod[3] = cats.ElShearMod; + SS_ref_db.ElShearMod[4] = di.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = he.C[i]; + SS_ref_db.Comp[1][i] = jd.C[i]; + SS_ref_db.Comp[2][i] = cen.C[i]; + SS_ref_db.Comp[3][i] = cats.C[i]; + SS_ref_db.Comp[4][i] = di.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + SS_ref_db.bounds_ref[2][0] = 0.0+eps; SS_ref_db.bounds_ref[2][1] = 1.0-eps; + SS_ref_db.bounds_ref[3][0] = 0.0+eps; SS_ref_db.bounds_ref[3][1] = 1.0-eps; + SS_ref_db.bounds_ref[4][0] = 0.0+eps; SS_ref_db.bounds_ref[4][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_hpcpx +*/ +SS_ref G_SS_sb11_hpcpx_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"hpcen","hpcfs",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 0.0; + + + em_data hpcen = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "hpcen", + "equilibrium"); + + em_data hpcfs = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "hpcfs", + "equilibrium"); + + SS_ref_db.gbase[0] = hpcen.gb; + SS_ref_db.gbase[1] = hpcfs.gb; + + SS_ref_db.ElShearMod[0] = hpcen.ElShearMod; + SS_ref_db.ElShearMod[1] = hpcfs.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = hpcen.C[i]; + SS_ref_db.Comp[1][i] = hpcfs.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_ak +*/ +SS_ref G_SS_sb11_ak_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"co","mgak","feak",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 66000.0; + SS_ref_db.W[1] = 0.0; + SS_ref_db.W[2] = 0.0; + + + em_data co = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "co", + "equilibrium"); + + em_data mgak = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "mgak", + "equilibrium"); + + em_data feak = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "feak", + "equilibrium"); + + SS_ref_db.gbase[0] = co.gb; + SS_ref_db.gbase[1] = mgak.gb; + SS_ref_db.gbase[2] = feak.gb; + + SS_ref_db.ElShearMod[0] = co.ElShearMod; + SS_ref_db.ElShearMod[1] = mgak.ElShearMod; + SS_ref_db.ElShearMod[2] = feak.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = co.C[i]; + SS_ref_db.Comp[1][i] = mgak.C[i]; + SS_ref_db.Comp[2][i] = feak.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + SS_ref_db.bounds_ref[2][0] = 0.0+eps; SS_ref_db.bounds_ref[2][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_gtmj +*/ +SS_ref G_SS_sb11_gtmj_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"alm","jdmj","mgmj","py","gr",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 30000.0; + SS_ref_db.W[1] = 0.0; + SS_ref_db.W[2] = 0.0; + SS_ref_db.W[3] = 0.0; + SS_ref_db.W[4] = 0.0; + SS_ref_db.W[5] = 0.0; + SS_ref_db.W[6] = 21300.0; + SS_ref_db.W[7] = 0.0; + SS_ref_db.W[8] = 0.0; + SS_ref_db.W[9] = 58000.0; + + + em_data alm = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "alm", + "equilibrium"); + + em_data jdmj = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "jdmj", + "equilibrium"); + + em_data mgmj = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "mgmj", + "equilibrium"); + + em_data py = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "py", + "equilibrium"); + + em_data gr = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "gr", + "equilibrium"); + + SS_ref_db.gbase[0] = alm.gb; + SS_ref_db.gbase[1] = jdmj.gb; + SS_ref_db.gbase[2] = mgmj.gb; + SS_ref_db.gbase[3] = py.gb; + SS_ref_db.gbase[4] = gr.gb; + + SS_ref_db.ElShearMod[0] = alm.ElShearMod; + SS_ref_db.ElShearMod[1] = jdmj.ElShearMod; + SS_ref_db.ElShearMod[2] = mgmj.ElShearMod; + SS_ref_db.ElShearMod[3] = py.ElShearMod; + SS_ref_db.ElShearMod[4] = gr.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = alm.C[i]; + SS_ref_db.Comp[1][i] = jdmj.C[i]; + SS_ref_db.Comp[2][i] = mgmj.C[i]; + SS_ref_db.Comp[3][i] = py.C[i]; + SS_ref_db.Comp[4][i] = gr.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + SS_ref_db.bounds_ref[2][0] = 0.0+eps; SS_ref_db.bounds_ref[2][1] = 1.0-eps; + SS_ref_db.bounds_ref[3][0] = 0.0+eps; SS_ref_db.bounds_ref[3][1] = 1.0-eps; + SS_ref_db.bounds_ref[4][0] = 0.0+eps; SS_ref_db.bounds_ref[4][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_pv +*/ +SS_ref G_SS_sb11_pv_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"alpv","fepv","mgpv",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 0.0; + SS_ref_db.W[1] = 116000.0; + SS_ref_db.W[2] = 0.0; + + SS_ref_db.v[0] = 1.0; + SS_ref_db.v[1] = 1.0; + SS_ref_db.v[2] = 0.39; + + em_data alpv = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "alpv", + "equilibrium"); + + em_data fepv = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "fepv", + "equilibrium"); + + em_data mgpv = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "mgpv", + "equilibrium"); + + SS_ref_db.gbase[0] = alpv.gb; + SS_ref_db.gbase[1] = fepv.gb; + SS_ref_db.gbase[2] = mgpv.gb; + + SS_ref_db.ElShearMod[0] = alpv.ElShearMod; + SS_ref_db.ElShearMod[1] = fepv.ElShearMod; + SS_ref_db.ElShearMod[2] = mgpv.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = alpv.C[i]; + SS_ref_db.Comp[1][i] = fepv.C[i]; + SS_ref_db.Comp[2][i] = mgpv.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + SS_ref_db.bounds_ref[2][0] = 0.0+eps; SS_ref_db.bounds_ref[2][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_ppv +*/ +SS_ref G_SS_sb11_ppv_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"feppv","mgppv","alppv",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 0.0; + SS_ref_db.W[1] = 0.0; + SS_ref_db.W[2] = 60000.0; + + + em_data feppv = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "feppv", + "equilibrium"); + + em_data mgppv = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "mgppv", + "equilibrium"); + + em_data alppv = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "alppv", + "equilibrium"); + + SS_ref_db.gbase[0] = feppv.gb; + SS_ref_db.gbase[1] = mgppv.gb; + SS_ref_db.gbase[2] = alppv.gb; + + SS_ref_db.ElShearMod[0] = feppv.ElShearMod; + SS_ref_db.ElShearMod[1] = mgppv.ElShearMod; + SS_ref_db.ElShearMod[2] = alppv.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = feppv.C[i]; + SS_ref_db.Comp[1][i] = mgppv.C[i]; + SS_ref_db.Comp[2][i] = alppv.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + SS_ref_db.bounds_ref[2][0] = 0.0+eps; SS_ref_db.bounds_ref[2][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_mw +*/ +SS_ref G_SS_sb11_mw_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"wu","pe",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 13000.0; + + + em_data wu = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "wu", + "equilibrium"); + + em_data pe = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "pe", + "equilibrium"); + + SS_ref_db.gbase[0] = wu.gb; + SS_ref_db.gbase[1] = pe.gb; + + SS_ref_db.ElShearMod[0] = wu.ElShearMod; + SS_ref_db.ElShearMod[1] = pe.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = wu.C[i]; + SS_ref_db.Comp[1][i] = pe.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + + + return SS_ref_db; +} + +/** + Solution phase data for sb11_cf +*/ +SS_ref G_SS_sb11_cf_function(SS_ref SS_ref_db, char* research_group, int EM_dataset, int len_ox, bulk_info z_b, double eps){ + + int i, j; + int n_em = SS_ref_db.n_em; + + char *EM_tmp[] = {"mgcf","nacf","fecf",}; + for (int i = 0; i < SS_ref_db.n_em; i++){ + strcpy(SS_ref_db.EM_list[i],EM_tmp[i]); + }; + + SS_ref_db.W[0] = 0.0; + SS_ref_db.W[1] = 0.0; + SS_ref_db.W[2] = 0.0; + + + em_data mgcf = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "mgcf", + "equilibrium"); + + em_data nacf = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "nacf", + "equilibrium"); + + em_data fecf = get_em_data( research_group, EM_dataset, + len_ox, + z_b, + SS_ref_db.P, + SS_ref_db.T, + "fecf", + "equilibrium"); + + SS_ref_db.gbase[0] = mgcf.gb; + SS_ref_db.gbase[1] = nacf.gb; + SS_ref_db.gbase[2] = fecf.gb; + + SS_ref_db.ElShearMod[0] = mgcf.ElShearMod; + SS_ref_db.ElShearMod[1] = nacf.ElShearMod; + SS_ref_db.ElShearMod[2] = fecf.ElShearMod; + + for (i = 0; i < len_ox; i++){ + SS_ref_db.Comp[0][i] = mgcf.C[i]; + SS_ref_db.Comp[1][i] = nacf.C[i]; + SS_ref_db.Comp[2][i] = fecf.C[i]; + } + + for (i = 0; i < n_em; i++){ + SS_ref_db.z_em[i] = 1.0; + }; + + SS_ref_db.bounds_ref[0][0] = 0.0+eps; SS_ref_db.bounds_ref[0][1] = 1.0-eps; + SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; + SS_ref_db.bounds_ref[2][0] = 0.0+eps; SS_ref_db.bounds_ref[2][1] = 1.0-eps; + + + return SS_ref_db; +} + + SS_ref G_SS_sb11_EM_function( global_variable gv, + SS_ref SS_ref_db, + int EM_dataset, + bulk_info z_b, + char *name){ + + double eps = gv.bnd_val; + double P = SS_ref_db.P; + double T = SS_ref_db.T; + + SS_ref_db.ss_flags[0] = 1; + + /* Associate the right solid-solution data */ + for (int FD = 0; FD < gv.n_Diff; FD++){ + + if (FD == 8 || FD == 9){ + SS_ref_db.P = 1.+ gv.gb_P_eps*gv.pdev[0][FD]; + SS_ref_db.T = T + gv.gb_T_eps*gv.pdev[1][FD]; + } + else{ + SS_ref_db.P = P + gv.gb_P_eps*gv.pdev[0][FD]; + SS_ref_db.T = T + gv.gb_T_eps*gv.pdev[1][FD]; + } + + if (strcmp( name, "plg") == 0 ){ + SS_ref_db = G_SS_sb11_plg_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "sp") == 0 ){ + SS_ref_db = G_SS_sb11_sp_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "ol") == 0 ){ + SS_ref_db = G_SS_sb11_ol_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "wa") == 0 ){ + SS_ref_db = G_SS_sb11_wa_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "ri") == 0 ){ + SS_ref_db = G_SS_sb11_ri_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "opx") == 0 ){ + SS_ref_db = G_SS_sb11_opx_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "cpx") == 0 ){ + SS_ref_db = G_SS_sb11_cpx_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "hpcpx") == 0 ){ + SS_ref_db = G_SS_sb11_hpcpx_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "ak") == 0 ){ + SS_ref_db = G_SS_sb11_ak_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "gtmj") == 0 ){ + SS_ref_db = G_SS_sb11_gtmj_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "pv") == 0 ){ + SS_ref_db = G_SS_sb11_pv_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "ppv") == 0 ){ + SS_ref_db = G_SS_sb11_ppv_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "mw") == 0 ){ + SS_ref_db = G_SS_sb11_mw_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else if (strcmp( name, "cf") == 0 ){ + SS_ref_db = G_SS_sb11_cf_function(SS_ref_db, gv.research_group, EM_dataset, gv.len_ox, z_b, eps); } + else{ + printf("\nsolid solution '%s' is not in the database\n",name); } + + for (int j = 0; j < SS_ref_db.n_em; j++){ + SS_ref_db.mu_array[FD][j] = SS_ref_db.gbase[j]; + } + } + + for (int j = 0; j < SS_ref_db.n_em; j++){ + SS_ref_db.bounds[j][0] = SS_ref_db.bounds_ref[j][0]; + SS_ref_db.bounds[j][1] = SS_ref_db.bounds_ref[j][1]; + } + + /* Calculate the number of atoms in the bulk-rock composition */ + double fbc = 0.0; + for (int i = 0; i < gv.len_ox; i++){ + fbc += z_b.bulk_rock[i]*z_b.apo[i]; + } + + for (int i = 0; i < SS_ref_db.n_em; i++){ + SS_ref_db.ape[i] = 0.0; + for (int j = 0; j < gv.len_ox; j++){ + SS_ref_db.ape[i] += SS_ref_db.Comp[i][j]*z_b.apo[j]; + } + } + + SS_ref_db.fbc = z_b.fbc; + + if (gv.verbose == 1){ + printf(" %4s:",name); + for (int j = 0; j < SS_ref_db.n_em; j++){ + printf(" %+12.5f",SS_ref_db.gbase[j]); + } + printf("\n"); + printf(" S C A F M N\n"); + for (int i = 0; i < SS_ref_db.n_em; i++){ + for (int j = 0; j < gv.len_ox; j++){ + printf(" %.1f",SS_ref_db.Comp[i][j]); + } + printf("\n"); + } + printf("\n"); + } + + return SS_ref_db; +}; diff --git a/src/SB_database/sb_gss_function.h b/src/SB_database/sb_gss_function.h new file mode 100644 index 0000000..11e15ca --- /dev/null +++ b/src/SB_database/sb_gss_function.h @@ -0,0 +1,22 @@ +/*@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ** + ** Project : MAGEMin + ** License : GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 + ** Developers : Nicolas Riel, Boris Kaus + ** Contributors : Dominguez, H., Assunção J., Green E., Berlie N., and Rummel L. + ** Organization : Institute of Geosciences, Johannes-Gutenberg University, Mainz + ** Contact : nriel[at]uni-mainz.de, kaus[at]uni-mainz.de + ** + ** ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @*/ +#ifndef __SB_GSS_FUNCTION_H_ +#define __SB_GSS_FUNCTION_H_ + +#include "../MAGEMin.h" + +SS_ref G_SS_sb11_EM_function( global_variable gv, + SS_ref SS_ref_db, + int EM_dataset, + bulk_info z_b, + char *name ); + +#endif diff --git a/src/SB_database/sb_gss_init_function.c b/src/SB_database/sb_gss_init_function.c index de99cb8..f90151e 100644 --- a/src/SB_database/sb_gss_init_function.c +++ b/src/SB_database/sb_gss_init_function.c @@ -221,7 +221,7 @@ SS_ref G_SS_sb11_cf_init_function(SS_ref SS_ref_db, global_variable gv){ return SS_ref_db; } -void SS_init_sb11( SS_init_type *SS_init, +void SB_SS_init_sb11( SS_init_type *SS_init, global_variable gv ){ for (int iss = 0; iss < gv.len_ss; iss++){ @@ -257,4 +257,14 @@ void SS_init_sb11( SS_init_type *SS_init, printf("\nsolid solution '%s' is not in the database, cannot be initiated\n", gv.SS_list[iss]); } } +} + +void SB_SS_init( SS_init_type *SS_init, + global_variable gv ){ + + if (gv.EM_database == 0){ // metapelite database // + SB_SS_init_sb11( SS_init, + gv ); + } + } \ No newline at end of file diff --git a/src/SB_database/sb_gss_init_function.h b/src/SB_database/sb_gss_init_function.h index 382cf29..76893a1 100644 --- a/src/SB_database/sb_gss_init_function.h +++ b/src/SB_database/sb_gss_init_function.h @@ -15,7 +15,7 @@ #include "../initialize.h" -void SS_init_sb11( SS_init_type *SS_init, - global_variable gv ); - +void SB_SS_init( SS_init_type *SS_init, + global_variable gv ); + #endif diff --git a/src/TC_database/TC_init_database.c b/src/TC_database/TC_init_database.c index e84589a..9784018 100644 --- a/src/TC_database/TC_init_database.c +++ b/src/TC_database/TC_init_database.c @@ -20,7 +20,6 @@ #include "../all_endmembers.h" #include "TC_init_database.h" - oxide_data oxide_info = { 15, /* number of endmembers */ {"SiO2" ,"Al2O3","CaO" ,"MgO" ,"FeO" ,"K2O" ,"Na2O" ,"TiO2" ,"O" ,"MnO" ,"Cr2O3","H2O" ,"CO2" ,"S" ,"Cl" }, @@ -69,16 +68,16 @@ metabasite_dataset metabasite_db = { {"SiO2" ,"Al2O3","CaO" ,"MgO" ,"FeO" ,"K2O" ,"Na2O" ,"TiO2" ,"O" ,"H2O" }, {"q" ,"crst" ,"trd" ,"coe" ,"law" ,"ky" ,"sill" ,"and" ,"ru" ,"sph" ,"ab" ,"H2O" , "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" ,"spn" }, + {"sp" ,"opx" ,"fsp" ,"liq" ,"mu" ,"ilmm" ,"ol" ,"hb" ,"ep" ,"g" ,"chl" ,"bi" ,"dio" ,"abc" ,"spn" }, {1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 }, // allow solvus? - {936 ,1731 ,231 ,3505 ,4536 ,420 ,11 ,7664 ,110 ,216 ,3980 ,1097 ,1808 ,21 ,196 }, // # of pseudocompound - {0.09 ,0.19 ,0.049 ,0.199 ,0.19 ,0.049 ,0.098 ,0.249 ,0.049 ,0.19 ,0.19 ,0.149 ,0.16 ,0.049 ,0.09 }, // discretization step + {936 ,1731 ,231 ,3505 ,4536 ,298 ,11 ,7664 ,110 ,216 ,3980 ,1097 ,1808 ,21 ,196 }, // # of pseudocompound + {0.09 ,0.19 ,0.049 ,0.199 ,0.19 ,0.09 ,0.098 ,0.249 ,0.049 ,0.19 ,0.19 ,0.149 ,0.16 ,0.049 ,0.09 }, // discretization step - {"sp" ,"opx" ,"fsp" ,"liq" ,"mu" ,"ilm" ,"ol" ,"hb" ,"ep" ,"g" ,"chl" ,"bi" ,"aug" ,"abc" ,"spn" }, + {"sp" ,"opx" ,"fsp" ,"liq" ,"mu" ,"ilmm" ,"ol" ,"hb" ,"ep" ,"g" ,"chl" ,"bi" ,"aug" ,"abc" ,"spn" }, {1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 }, // allow solvus? - {936 ,1731 ,231 ,3505 ,4536 ,420 ,11 ,7664 ,110 ,216 ,3980 ,1097 ,2395 ,21 ,196 }, // # of pseudocompound - {0.09 ,0.19 ,0.049 ,0.199 ,0.19 ,0.049 ,0.098 ,0.249 ,0.049 ,0.19 ,0.19 ,0.149 ,0.24 ,0.049 ,0.09 }, // discretization step + {936 ,1731 ,231 ,3505 ,4536 ,298 ,11 ,7664 ,110 ,216 ,3980 ,1097 ,2395 ,21 ,196 }, // # of pseudocompound + {0.09 ,0.19 ,0.049 ,0.199 ,0.19 ,0.09 ,0.098 ,0.249 ,0.049 ,0.19 ,0.19 ,0.149 ,0.24 ,0.049 ,0.09 }, // discretization step 4.0, /* max dG under which a phase is considered to be reintroduced */ 473.15, /* max temperature above which PGE solver is active */ @@ -696,6 +695,7 @@ global_variable global_variable_TC_init( global_variable gv, gv.K2O_id = -1; gv.O_id = -1; gv.MnO_id = -1; + oxide_data ox_in = oxide_info; for (i = 0; i < gv.len_ox; i++){ for (j = 0; j < ox_in.n_ox; j++){ diff --git a/src/TC_database/objective_functions.c b/src/TC_database/objective_functions.c index 18c21e0..8e18730 100644 --- a/src/TC_database/objective_functions.c +++ b/src/TC_database/objective_functions.c @@ -1954,6 +1954,7 @@ double obj_mb_ilmm(unsigned n, const double *x, double *grad, void *SS_ref_db){ double *mu_Gex = d->mu_Gex; double *sf = d->sf; double *mu = d->mu; + double *d_em = d->d_em; px_mb_ilmm(SS_ref_db,x); for (int i = 0; i < n_em; i++){ @@ -1978,7 +1979,7 @@ double obj_mb_ilmm(unsigned n, const double *x, double *grad, void *SS_ref_db){ mu[0] = gb[0] + R*T*creal(clog(sf[0]*sf[5])) + mu_Gex[0]; mu[1] = gb[1] + R*T*creal(clog(4.0*csqrt(sf[0])*csqrt(sf[1])*csqrt(sf[4])*csqrt(sf[5]))) + mu_Gex[1]; - mu[2] = gb[2] + R*T*creal(clog(sf[3]*sf[6])) + mu_Gex[2]; + mu[2] = gb[2] + R*T*creal(clog(sf[3]*sf[6] + d_em[2])) + mu_Gex[2]; mu[3] = gb[3] + R*T*creal(clog(sf[2]*sf[5])) + mu_Gex[3]; d->sum_apep = 0.0; diff --git a/src/TC_database/tc_gss_function.c b/src/TC_database/tc_gss_function.c index df212ba..2f884ab 100644 --- a/src/TC_database/tc_gss_function.c +++ b/src/TC_database/tc_gss_function.c @@ -795,13 +795,12 @@ SS_ref G_SS_mb_aug_function(SS_ref SS_ref_db, char* research_group, int EM_datas SS_ref_db.bounds_ref[6][0] = 0.0+eps; SS_ref_db.bounds_ref[6][1] = 1.0-eps; if (z_b.bulk_rock[8] == 0.){ //O - SS_ref_db.z_em[5] = 0.0; + SS_ref_db.z_em[4] = 0.0; SS_ref_db.d_em[4] = 1.0; - SS_ref_db.bounds_ref[3][0] = 0.0; - SS_ref_db.bounds_ref[3][1] = 0.0; + SS_ref_db.bounds_ref[2][0] = 0.0; + SS_ref_db.bounds_ref[2][1] = 0.0; } - return SS_ref_db; } @@ -1827,7 +1826,12 @@ SS_ref G_SS_mb_ilmm_function(SS_ref SS_ref_db, char* research_group, int EM_data SS_ref_db.bounds_ref[1][0] = 0.0+eps; SS_ref_db.bounds_ref[1][1] = 1.0-eps; SS_ref_db.bounds_ref[2][0] = 0.0+eps; SS_ref_db.bounds_ref[2][1] = 1.0-eps; - + if (z_b.bulk_rock[8] == 0.){ //O + SS_ref_db.z_em[2] = 0.0; + SS_ref_db.d_em[2] = 1.0; + SS_ref_db.bounds_ref[0][0] = 1.0; + SS_ref_db.bounds_ref[0][1] = 1.0; + } return SS_ref_db; } @@ -11612,8 +11616,6 @@ SS_ref G_SS_mb_EM_function( global_variable gv, if (gv.verbose == 1){ printf(" %4s:",name); - - /* display Gibbs free energy of reference? */ for (int j = 0; j < SS_ref_db.n_em; j++){ printf(" %+12.5f",SS_ref_db.gbase[j]); } diff --git a/src/all_endmembers.h b/src/all_endmembers.h index a00318b..0fb49cf 100644 --- a/src/all_endmembers.h +++ b/src/all_endmembers.h @@ -16,6 +16,7 @@ #include "./TC_database/TC_gem_function.h" /* This include Stixrude & Lithgow-Bertelloni solution phase models */ #include "./SB_database/SB_endmembers.h" + #include "./SB_database/SB_gem_function.h" //#include "STIX_solution_phases.h" #endif \ No newline at end of file diff --git a/src/all_solution_phases.h b/src/all_solution_phases.h index 8d8ebda..5526f91 100644 --- a/src/all_solution_phases.h +++ b/src/all_solution_phases.h @@ -15,6 +15,6 @@ #include "./TC_database/TC_solution_phases.h" /* This include Stixrude & Lithgow-Bertelloni solution phase models */ - //#include "STIX_solution_phases.h" + #include "./SB_database/SB_solution_phases.h" #endif \ No newline at end of file diff --git a/src/dump_function.c b/src/dump_function.c index 2e7e409..2b065c5 100644 --- a/src/dump_function.c +++ b/src/dump_function.c @@ -87,7 +87,7 @@ void reset_output_struct( global_variable gv, double sum_Molar_mass_bulk; strcpy(sp[0].MAGEMin_ver,gv.version); - + strcpy(sp[0].database,gv.db); if (gv.EM_dataset == 62){ strcpy(sp[0].dataset,"tc_ds62"); } @@ -118,6 +118,15 @@ void reset_output_struct( global_variable gv, sp[0].aAl2O3 = gv.system_aAl2O3; sp[0].aMgO = gv.system_aMgO; sp[0].aFeO = gv.system_aFeO; + strcpy(sp[0].buffer,gv.buffer); + + // if (strcmp(gv.buffer, "NONE") != 0){ + // sp[0].buffer_n = gv.buffer_n; + // } + // else{ + // sp[0].buffer_n = 0.0; + // } + sp[0].buffer_n = gv.buffer_n; sp[0].alpha = gv.system_expansivity; sp[0].V = gv.system_volume*10.0; @@ -460,7 +469,7 @@ void fill_output_struct( global_variable gv, sp[0].ph_frac[n] = cp[i].ss_n_mol; sp[0].ph_frac_wt[n] = cp[i].ss_n_wt; - + sp[0].ph_frac_1at[n] = cp[i].ss_n; sp[0].ph_type[n] = 1; sp[0].ph_id[n] = m; sp[0].n_SS += 1; @@ -610,6 +619,7 @@ void fill_output_struct( global_variable gv, sp[0].ph_frac[n] = gv.pp_n_mol[i]; sp[0].ph_frac_wt[n] = gv.pp_n_wt[i]; + sp[0].ph_frac_1at[n] = gv.pp_n[i]; sp[0].ph_type[n] = 0; sp[0].ph_id[n] = m; sp[0].n_PP += 1; diff --git a/src/gem_function.c b/src/gem_function.c index dafd1a6..f37f68a 100644 --- a/src/gem_function.c +++ b/src/gem_function.c @@ -48,8 +48,8 @@ PP_ref G_EM_function( char *research_group, name, state); } - else{ - PP_ref_db = TC_G_EM_function( EM_dataset, + if (strcmp(research_group, "sb") == 0){ + PP_ref_db = SB_G_EM_function( EM_dataset, len_ox, id, bulk_rock, diff --git a/src/initialize.c b/src/initialize.c index 75fa2a0..07111a6 100644 --- a/src/initialize.c +++ b/src/initialize.c @@ -12,6 +12,7 @@ #include "MAGEMin.h" #include "initialize.h" #include "toolkit.h" +#include /* Function to allocate the memory of the data to be used/saved during PGE iterations */ global_variable global_variable_init( global_variable gv, @@ -28,6 +29,9 @@ global_variable global_variable_init( global_variable gv, gv = global_variable_SB_init( gv, z_b ); } + else{ + printf(" wrong group, fix group name\n"); + } return gv; } @@ -53,6 +57,7 @@ char** get_EM_DB_names_tc(global_variable gv) { } char** get_EM_DB_names_sb(global_variable gv) { + EM_db_sb EM_return; int i, n_em_db; n_em_db = gv.n_em_db; @@ -117,7 +122,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.5.6 [17/10/2024]"); /** MAGEMin version */ + strcpy(gv.version,"1.5.7 [17/10/2024]"); /** MAGEMin version */ /* generate parameters */ strcpy(gv.buffer,"none"); @@ -197,7 +202,7 @@ global_variable global_variable_alloc( bulk_info *z_b ){ gv.tot_time = 0.0; /* set default parameters (overwritten later from args)*/ - gv.EM_database = 2; + gv.EM_database = 0; gv.n_points = 1; gv.solver = 2; /* 0 -> Legacy, 1 = PGE, Hybrid PGE/LP */ gv.leveling_mode = 0; @@ -269,6 +274,7 @@ stb_system SP_INIT_function(stb_system sp, global_variable gv){ sp.MAGEMin_ver = malloc(50 * sizeof(char) ); sp.dataset = malloc(50 * sizeof(char) ); + sp.database = malloc(50 * sizeof(char) ); sp.oxides = malloc(gv.len_ox * sizeof(char*) ); sp.elements = malloc(gv.len_ox * sizeof(char*) ); @@ -276,6 +282,7 @@ stb_system SP_INIT_function(stb_system sp, global_variable gv){ sp.oxides[i] = malloc(20 * sizeof(char)); sp.elements[i] = malloc(20 * sizeof(char)); } + sp.buffer = malloc(50 * sizeof(char) ); sp.bulk = malloc(gv.len_ox * sizeof(double) ); sp.gamma = malloc(gv.len_ox * sizeof(double) ); sp.bulk_S = malloc(gv.len_ox * sizeof(double) ); @@ -289,6 +296,7 @@ stb_system SP_INIT_function(stb_system sp, global_variable gv){ sp.ph = malloc(gv.len_ox * sizeof(char*) ); sp.ph_frac = malloc(gv.len_ox * sizeof(double) ); sp.ph_frac_wt = malloc(gv.len_ox * sizeof(double) ); + sp.ph_frac_1at = malloc(gv.len_ox * sizeof(double) ); sp.ph_frac_vol = malloc(gv.len_ox * sizeof(double) ); for (int i = 0; i < gv.len_ox; i++){ sp.ph[i] = malloc(20 * sizeof(char)); @@ -776,6 +784,7 @@ void reset_sp( global_variable gv, sp[0].ph_id[i] = 0; sp[0].ph_frac[i] = 0.0; sp[0].ph_frac_wt[i] = 0.0; + sp[0].ph_frac_1at[i] = 0.0; sp[0].ph_frac_vol[i] = 0.0; } diff --git a/src/pp_min_function.c b/src/pp_min_function.c index 7130ba7..eb90222 100644 --- a/src/pp_min_function.c +++ b/src/pp_min_function.c @@ -59,6 +59,7 @@ global_variable init_em_db( int EM_database, ){ /* initialize endmember database */ + double buffer_n; char state[] = "equilibrium"; int sum_zel; for (int i = 0; i < gv.len_pp; i++){ @@ -200,13 +201,15 @@ global_variable init_em_db( int EM_database, /* Calculate normalizing factor */ double factor = fbc/ape; if (gv.buffer_n < 1e-8){ - gv.buffer_n = 1e-8; + buffer_n = 1e-8; } else if (gv.buffer_n >= 1.0){ - gv.buffer_n = 1.0-1e-8; + buffer_n = 1.0-1e-8; } - - PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + H2O.gbase; + else{ + buffer_n = gv.buffer_n; + } + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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; @@ -242,14 +245,16 @@ global_variable init_em_db( int EM_database, } /* Calculate normalizing factor */ double factor = fbc/ape; - if (gv.buffer_n <= 1e-60){ - gv.buffer_n = 1e-60; + if (gv.buffer_n <= 1e-8){ + buffer_n = 1e-8; } else if (gv.buffer_n >= 1.0){ - gv.buffer_n = 1.0-1e-8; + buffer_n = 1.0-1e-8; } - - PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + O2.gbase; + else{ + buffer_n = gv.buffer_n; + } + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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; @@ -287,13 +292,15 @@ global_variable init_em_db( int EM_database, /* Calculate normalizing factor */ double factor = fbc/ape; if (gv.buffer_n <= 1e-8){ - gv.buffer_n = 1e-8; + buffer_n = 1e-8; } else if (gv.buffer_n >= 1.0){ - gv.buffer_n = 1.0-1e-8; + buffer_n = 1.0-1e-8; } - - PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + MgO.gbase; + else{ + buffer_n = gv.buffer_n; + } + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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; @@ -330,13 +337,15 @@ global_variable init_em_db( int EM_database, /* Calculate normalizing factor */ double factor = fbc/ape; if (gv.buffer_n <= 1e-8){ - gv.buffer_n = 1e-8; + buffer_n = 1e-8; } else if (gv.buffer_n >= 1.0){ - gv.buffer_n = 1.0-1e-8; + buffer_n = 1.0-1e-8; } - - PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + FeO.gbase; + else{ + buffer_n = gv.buffer_n; + } + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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; @@ -373,13 +382,15 @@ global_variable init_em_db( int EM_database, /* Calculate normalizing factor */ double factor = fbc/ape; if (gv.buffer_n <= 1e-8){ - gv.buffer_n = 1e-8; + buffer_n = 1e-8; } else if (gv.buffer_n >= 1.0){ - gv.buffer_n = 1.0-1e-8; + buffer_n = 1.0-1e-8; } - - PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + Al2O3.gbase; + else{ + buffer_n = gv.buffer_n; + } + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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; @@ -387,7 +398,7 @@ global_variable init_em_db( int EM_database, else if (strcmp( gv.PP_list[i], "aTiO2") == 0){ PP_ref TiO2 = G_EM_function( gv.research_group, - gv.EM_dataset, + gv.EM_dataset, gv.len_ox, z_b.id, z_b.bulk_rock, @@ -416,13 +427,16 @@ global_variable init_em_db( int EM_database, /* Calculate normalizing factor */ double factor = fbc/ape; if (gv.buffer_n <= 1e-8){ - gv.buffer_n = 1e-8; + buffer_n = 1e-8; } else if (gv.buffer_n >= 1.0){ - gv.buffer_n = 1.0-1e-8; + buffer_n = 1.0-1e-8; + } + else{ + buffer_n = gv.buffer_n; } - PP_ref_db[i].gbase = z_b.R * z_b.T*log(gv.buffer_n) + TiO2.gbase; + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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; @@ -722,3 +736,216 @@ global_variable init_em_db( int EM_database, } return gv; }; + + + +/** + initialize pure phase database */ +global_variable init_em_db_sb( int EM_database, + bulk_info z_b, + global_variable gv, + PP_ref *PP_ref_db +){ + + /* initialize endmember database */ + double buffer_n; + char state[] = "equilibrium"; + int sum_zel; + for (int i = 0; i < gv.len_pp; i++){ + + if (strcmp( gv.PP_list[i], "aMgO") == 0){ + + PP_ref MgO = G_EM_function( gv.research_group, + gv.EM_dataset, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "pe", + 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){ + buffer_n = 1e-8; + } + else if (gv.buffer_n >= 1.0){ + buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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( gv.research_group, + gv.EM_dataset, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "wu", + 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){ + buffer_n = 1e-8; + } + else if (gv.buffer_n >= 1.0){ + buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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( gv.research_group, + gv.EM_dataset, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + "co", + 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){ + buffer_n = 1e-8; + } + else if (gv.buffer_n >= 1.0){ + buffer_n = 1.0-1e-8; + } + + PP_ref_db[i].gbase = z_b.R * z_b.T*log(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{ + PP_ref_db[i] = G_EM_function( gv.research_group, + gv.EM_dataset, + gv.len_ox, + z_b.id, + z_b.bulk_rock, + z_b.apo, + z_b.P, + z_b.T, + gv.PP_list[i], + state ); + } + + sum_zel = 0; + for (int j = 0; j < z_b.zEl_val; j++){ + + /* If pure-phase contains an oxide absent in the bulk-rock then do not take it into account */ + if (PP_ref_db[i].Comp[z_b.zEl_array[j]] != 0.0){ + sum_zel += 1; + } + + } + + /* If pure-phase contains an oxide absent in the bulk-rock then do not take it into account */ + if (sum_zel != 0){ + gv.pp_flags[i][0] = 0; + gv.pp_flags[i][1] = 0; + gv.pp_flags[i][2] = 0; + gv.pp_flags[i][3] = 1; + } + else{ + if (gv.pp_flags[i][0] != 0){ //here we check if the pure phase is deactivated from the start (O2 for instance) + gv.pp_flags[i][0] = 1; + gv.pp_flags[i][1] = 0; + gv.pp_flags[i][2] = 1; + gv.pp_flags[i][3] = 0; + } + } + + /* If buffer not active then remove it */ + if (gv.pp_flags[i][4] == 1 && strcmp(gv.buffer,gv.PP_list[i]) != 0){ + gv.pp_flags[i][0] = 0; + gv.pp_flags[i][1] = 0; + gv.pp_flags[i][2] = 0; + gv.pp_flags[i][3] = 1; + } + + if (gv.verbose==1){ + printf("\n %4s: %+10f %+10f\n",gv.PP_list[i],PP_ref_db[i].gbase, PP_ref_db[i].factor); + + /* display molar composition */ + + if (EM_database == 0){ + printf(" S C A F G N\n"); + } + for (int j = 0; j < gv.len_ox; j++){ + printf(" %.1f",PP_ref_db[i].Comp[j]); + } + printf("\n"); + } + + } + if (gv.verbose==1){ + printf("\n"); + } + return gv; +}; diff --git a/src/pp_min_function.h b/src/pp_min_function.h index 44e252c..0b3a63e 100644 --- a/src/pp_min_function.h +++ b/src/pp_min_function.h @@ -23,5 +23,11 @@ global_variable init_em_db( int EM_database, global_variable gv, PP_ref *PP_ref_db ); + +global_variable init_em_db_sb( int EM_database, + bulk_info z_b, + + global_variable gv, + PP_ref *PP_ref_db ); #endif diff --git a/src/toolkit.c b/src/toolkit.c index c263ef0..fc3521a 100644 --- a/src/toolkit.c +++ b/src/toolkit.c @@ -768,6 +768,51 @@ void print_cp( global_variable gv, }; +global_variable get_tests_bulks( global_variable gv +){ + + if ( strcmp(gv.research_group, "tc") == 0 ){ + if (gv.EM_database == 0){ + gv = get_bulk_metapelite( gv ); + } + else if (gv.EM_database == 1){ + gv = get_bulk_metabasite( gv ); + } + else if (gv.EM_database == 2){ + gv = get_bulk_igneous( gv ); + } + else if (gv.EM_database == 3){ + gv = get_bulk_igneous_igad( gv ); + } + else if (gv.EM_database == 4){ + gv = get_bulk_ultramafic( gv ); + } + else if (gv.EM_database == 5){ + gv = get_bulk_ultramafic_ext( gv ); + } + else if (gv.EM_database == 6){ + gv = get_bulk_mantle( gv ); + } + else{ + printf(" Wrong database...\n"); + } + } + else if ( strcmp(gv.research_group, "sb") == 0 ){ + if (gv.EM_database == 0){ + gv = get_bulk_stx11( gv ); + } + else{ + printf(" Wrong database...\n"); + } + } + else{ + printf(" Wrong research group...\n"); + } + + + return gv; +} + /** rotate G-hyperplane using Gamma */ diff --git a/src/toolkit.h b/src/toolkit.h index fc1b966..66ae47f 100644 --- a/src/toolkit.h +++ b/src/toolkit.h @@ -24,7 +24,7 @@ bulk_info retrieve_bulk_PT( global_variable gv, void convert_system_comp( global_variable gv, char *sys_in, bulk_info z_b ); - +global_variable get_tests_bulks( global_variable gv ); void get_act_sf_id(int *result, double *A, int n); void inverseMatrix(int *ipiv, double *A1, int n, double *work, int lwork); void MatMatMul( double **A, int nrowA, double **B, int ncolB, int common, double **C); diff --git a/test/test_diagram_test0_mb.jl b/test/test_diagram_test0_mb.jl index 539d666..357dad7 100644 --- a/test/test_diagram_test0_mb.jl +++ b/test/test_diagram_test0_mb.jl @@ -1,55 +1,55 @@ # This checks the diagram produced by test4 (=Gt-migmatite) # -list = [outP{Float64}(0.01, 600.0, 0, -833.0108842123142, ["ilm", "aug", "fsp", "sp", "ilm", "opx", "fsp", "q", "H2O"], [0.015066457359944928, 0.21330918840680843, 0.37352098881926343, 0.02099544269247504, 0.018613772500099955, 0.16775758248197603, 0.015763356043942824, 0.008306545028817468, 0.16666666666667201]), -outP{Float64}(0.01, 700.0, 0, -846.030076204107, ["opx", "ilm", "fsp", "fsp", "aug", "sp", "q", "H2O"], [0.1708405002181094, 0.0413703417323135, 0.014904668123904551, 0.37221372361337474, 0.21902686309645347, 0.009141779388255328, 0.005835723094348328, 0.1666664007332406]), -outP{Float64}(0.01, 800.0, 0, -859.7443487939465, ["fsp", "opx", "fsp", "sp", "ilm", "aug", "q", "H2O"], [0.37109714409419275, 0.16595381815840637, 0.013093386256276013, 0.010404854342544385, 0.03910801034362012, 0.22678129081409082, 0.0069076054563125934, 0.1666538905345569]), -outP{Float64}(0.01, 900.0, 0, -874.1055573144797, ["ilm", "fsp", "opx", "sp", "liq", "aug", "H2O"], [0.03493310504272607, 0.3664670482139637, 0.15743562570577727, 0.014079206827789415, 0.024727896899986633, 0.23587163965875077, 0.16648547765100596]), -outP{Float64}(0.01, 1000.0, 0, -889.0818629821189, ["ilm", "liq", "fsp", "opx", "sp", "aug", "H2O"], [0.0321675901796047, 0.06569414419793271, 0.33713226178904815, 0.1370353944295203, 0.015541922739967646, 0.2462316129261662, 0.1661970737377603]), -outP{Float64}(0.01, 1100.0, 0, -904.7443004740051, ["ilm", "liq", "fsp", "aug", "opx", "H2O"], [0.03936591095366488, 0.25456610684119546, 0.22792952358033264, 0.24643912281652436, 0.06656101172688678, 0.1651383240813959]), -outP{Float64}(0.01, 1200.0, 0, -921.1825225074177, ["aug", "liq", "ilm", "fsp", "ru", "H2O"], [0.23176720111217272, 0.4565592797555339, 0.03223806232587665, 0.11375715765115108, 0.0016633103144025464, 0.1640149888408629]), -outP{Float64}(4.01, 600.0, 0, -820.2199888253759, ["hb", "fsp", "bi", "q", "sph", "H2O"], [0.5742223855506291, 0.20467993635909937, 0.014794632544438937, 0.05645359212182181, 0.021161223050293474, 0.12868823037371724]), -outP{Float64}(4.01, 700.0, 0, -831.9231107278648, ["liq", "hb", "fsp", "aug", "q", "sph", "H2O"], [0.03073724958102105, 0.5785327991509318, 0.19613409378717317, 0.006684737990748561, 0.04703911651538412, 0.017326884864591035, 0.12354511811015027]), -outP{Float64}(4.01, 800.0, 0, -844.4825511858518, ["fsp", "liq", "hb", "aug", "ilm", "H2O"], [0.13650724208005527, 0.20127589442738109, 0.47970676799648887, 0.08068713490536196, 0.01367875654475093, 0.08814420404596192]), -outP{Float64}(4.01, 900.0, 0, -857.8815734580671, ["liq", "ilm", "fsp", "hb", "aug", "H2O"], [0.3055754381980378, 0.013755253588464082, 0.08006406355313087, 0.38824972166949745, 0.1422621402088177, 0.07009338278205204]), -outP{Float64}(4.01, 1000.0, 0, -872.1437684091422, ["fsp", "hb", "liq", "aug", "ilm", "H2O"], [0.008060181580642502, 0.20124938101313386, 0.5138236781026841, 0.21900021569403705, 0.01994536701915965, 0.03792117659034282]), -outP{Float64}(4.01, 1100.0, 0, -887.3798023847376, ["aug", "liq", "ilm", "H2O"], [0.23886003707111372, 0.7167459523052979, 0.036987756739844654, 0.007406253883743712]), -outP{Float64}(4.01, 1200.0, 0, -903.42133180305, ["aug", "liq", "ilm", "ru"], [0.1396869716659256, 0.8182826865187091, 0.04130218681382433, 0.0007281550015410509]), -outP{Float64}(8.01, 600.0, 0, -812.2017692763923, ["ep", "hb", "fsp", "bi", "aug", "q", "sph", "H2O"], [0.05626264946273827, 0.583322992571966, 0.09055376267854699, 0.013927218050824575, 0.019260308534057653, 0.09236351733755692, 0.020752174849651455, 0.12355737651465827]), -outP{Float64}(8.01, 700.0, 0, -823.7713960670384, ["hb", "aug", "liq", "fsp", "q", "sph", "H2O"], [0.5732341578150781, 0.05234203037733809, 0.09170281963231594, 0.09938930553651588, 0.06727967123325339, 0.0171722436221721, 0.0988797717833264]), -outP{Float64}(8.01, 800.0, 0, -836.2695817050887, ["hb", "fsp", "aug", "liq", "sph", "H2O"], [0.5272905998650925, 0.01359904166697939, 0.09151655596003276, 0.33125632565944846, 0.011910436887120193, 0.024427039961326755]), -outP{Float64}(8.01, 900.0, 0, -849.6221172528484, ["hb", "liq", "ilm", "aug", "H2O"], [0.40047041862703925, 0.42387888936690293, 0.011660339482577149, 0.16203093556467624, 0.00195941695880458]), -outP{Float64}(8.01, 1000.0, 0, -863.8192021406793, ["aug", "hb", "liq", "ilm"], [0.2240264187828929, 0.2210430083167862, 0.5349582488123525, 0.01997232408796836]), -outP{Float64}(8.01, 1100.0, 0, -878.8922861419654, ["liq", "opx", "aug", "ilm", "ru"], [0.6811065291061783, 0.012096857909813108, 0.27415362344188415, 0.03143253270444494, 0.0012104568376796081]), -outP{Float64}(8.01, 1200.0, 0, -894.7552181662521, ["liq", "ilm", "aug", "ru"], [0.7630388750518338, 0.03306778519169642, 0.20166809071924116, 0.002225249037228865]), -outP{Float64}(12.01, 600.0, 0, -804.5540942899241, ["aug", "mu", "hb", "bi", "ep", "q", "sph", "H2O"], [0.029805522752731343, 0.01066501685848931, 0.5948008972118187, 0.007718241589938595, 0.11530350671635815, 0.10381688219541273, 0.020639569654606192, 0.11725036302064495]), -outP{Float64}(12.01, 700.0, 0, -815.9711812535147, ["hb", "ep", "liq", "aug", "q", "sph", "H2O"], [0.5844977976485889, 0.06668992144506473, 0.1457460616286092, 0.04375082144958598, 0.0778248337733556, 0.01716072964028762, 0.06432983441450799]), -outP{Float64}(12.01, 800.0, 0, -828.3529905460357, ["aug", "g", "hb", "liq", "ru", "sph"], [0.1226427340754218, 0.037144224770927145, 0.47022907586490403, 0.36422179752901773, 0.004530063192089444, 0.001232104567639927]), -outP{Float64}(12.01, 900.0, 0, -841.6161360283872, ["aug", "liq", "g", "hb", "ilm", "ru"], [0.17385889264010443, 0.4093641989682698, 0.03555575712481469, 0.37400434109717823, 0.003919301896773425, 0.003297508272859409]), -outP{Float64}(12.01, 1000.0, 0, -855.6903400794105, ["liq", "hb", "aug", "ilm"], [0.5024808519048559, 0.24417675011488144, 0.23479927535537, 0.018543122624892532]), -outP{Float64}(12.01, 1100.0, 0, -870.6002345373444, ["ilm", "liq", "aug", "opx", "hb", "ru"], [0.024731157635702276, 0.6330483406038778, 0.2933223698627382, 0.014881334533104851, 0.03185338126990946, 0.0021634160946672892]), -outP{Float64}(12.01, 1200.0, 0, -886.302160497227, ["ilm", "liq", "aug", "ru"], [0.02607848709509566, 0.7196851152848656, 0.25063599629706473, 0.0036004013229741683]), -outP{Float64}(16.01, 600.0, 0, -797.1277675368204, ["hb", "ep", "mu", "aug", "q", "ru", "H2O"], [0.5883910486462253, 0.18807746840825912, 0.019863003863574633, 0.016040525229556176, 0.06803522535901194, 0.007311340736579696, 0.11228138775679321]), -outP{Float64}(16.01, 700.0, 0, -808.3898627997189, ["hb", "ep", "g", "liq", "aug", "q", "ru", "H2O"], [0.5174504279027126, 0.10022493250019461, 0.038757823216333934, 0.08118784329060034, 0.08025343957443823, 0.08719931976249891, 0.0062809577105231735, 0.08864525604269839]), -outP{Float64}(16.01, 800.0, 0, -820.6592203636367, ["g", "aug", "liq", "hb", "ru"], [0.22924389867933836, 0.11241066521729698, 0.3646907212525616, 0.286761381250162, 0.006893333600641061]), -outP{Float64}(16.01, 900.0, 0, -833.834375931435, ["g", "liq", "hb", "aug", "ru"], [0.2606899168898933, 0.40628921629456216, 0.15643097745103998, 0.16935380839274947, 0.0072360809717552195]), -outP{Float64}(16.01, 1000.0, 0, -847.7800590397177, ["liq", "hb", "aug", "g", "ru"], [0.4686132591631269, 0.03948409513398686, 0.23316589510608643, 0.25053434769596233, 0.00820240290083755]), -outP{Float64}(16.01, 1100.0, 0, -862.5039592177286, ["aug", "liq", "ilm", "g", "ru"], [0.27444368181593776, 0.5696100849167286, 0.0037687884033253414, 0.14430729326947445, 0.007870151594533865]), -outP{Float64}(16.01, 1200.0, 0, -878.0214094713955, ["liq", "aug", "ilm", "ru"], [0.6837965591669669, 0.2913815012714052, 0.019973416641598674, 0.004848522920029092]), -outP{Float64}(20.01, 600.0, 0, -789.9128501808186, ["hb", "ep", "mu", "hb", "aug", "g", "q", "ru", "H2O"], [0.08196164939751585, 0.22680652801614817, 0.024359998183413424, 0.38221072327973393, 0.0019139922045525524, 0.12197501116522941, 0.035668738365448295, 0.008181525464035935, 0.1169218339239223]), -outP{Float64}(20.01, 700.0, 0, -801.0558445259801, ["hb", "g", "mu", "ep", "aug", "q", "ru", "H2O"], [0.34366866194243473, 0.1551257423335992, 0.0222135171831378, 0.10577743322050243, 0.14554069900657224, 0.08689837667652378, 0.007479247306024949, 0.1332963223312049]), -outP{Float64}(20.01, 800.0, 0, -813.1903789052379, ["ep", "aug", "hb", "liq", "g", "q", "ru"], [0.007175584263162252, 0.18693600885351827, 0.08405271063842183, 0.36073186162055954, 0.3517788293628542, 0.0010478530685901287, 0.008277152192893846]), -outP{Float64}(20.01, 900.0, 0, -826.2901698737817, ["g", "liq", "aug", "ru"], [0.3888858103984635, 0.3997264085288565, 0.2026378882028045, 0.00874989286987549]), -outP{Float64}(20.01, 1000.0, 0, -840.102363551254, ["liq", "g", "aug", "ru"], [0.43719185357710105, 0.3504886967581837, 0.20356955477876723, 0.00874989488594796]), -outP{Float64}(20.01, 1100.0, 0, -854.630614944856, ["g", "liq", "aug", "ru"], [0.2549562813153897, 0.5110677553322962, 0.22522607059305552, 0.008749892759258568]), -outP{Float64}(20.01, 1200.0, 0, -869.9168035037434, ["aug", "liq", "g", "ru"], [0.248932736408965, 0.6156298471551765, 0.12668751498727385, 0.008749901448584582]), -outP{Float64}(24.01, 600.0, 0, -782.84742311144, ["hb", "g", "mu", "ep", "aug", "q", "ru", "H2O"], [0.36592957042833746, 0.1872158498451617, 0.024699950118928216, 0.187646113429018, 0.056905863473430775, 0.04353002351030021, 0.008348953540858229, 0.12572367565396525]), -outP{Float64}(24.01, 700.0, 0, -793.9431958387042, ["aug", "g", "ep", "mu", "q", "ky", "ru", "H2O"], [0.41032284853577283, 0.27603206271219344, 0.0105253511940095, 0.02577124564626581, 0.06524083891303127, 0.04148660188051113, 0.008749916321166513, 0.1618711347970496]), -outP{Float64}(24.01, 800.0, 0, -805.9341097048041, ["liq", "g", "aug", "q", "ky", "ru", "H2O"], [0.21440435168980138, 0.3518681928155183, 0.3130924813676721, 0.028834702842417768, 0.021597314893405266, 0.008749871993171142, 0.06145308439801379]), -outP{Float64}(24.01, 900.0, 0, -818.8929677230251, ["g", "liq", "aug", "ru"], [0.39587949658390875, 0.3776737381875919, 0.21769687252273465, 0.008749892705764616]), -outP{Float64}(24.01, 1000.0, 0, -832.5705126811523, ["g", "liq", "aug", "ru"], [0.3761899731157409, 0.4110422405996056, 0.2040178948538423, 0.008749891430811136]), -outP{Float64}(24.01, 1100.0, 0, -846.9302236340743, ["aug", "liq", "g", "ru"], [0.19616098227742723, 0.46609855020829355, 0.32899056693343526, 0.008749900580844044]), -outP{Float64}(24.01, 1200.0, 0, -862.0078022924001, ["liq", "g", "aug", "ru"], [0.5608276061297298, 0.222596459693562, 0.2078260382616352, 0.008749895915072973]) +list = [outP{Float64}(0.01, 600.0, 0, -833.018591503985, ["sp", "ilmm", "fsp", "fsp", "opx", "aug", "ilmm", "q", "H2O"], [0.028556708629588434, 0.02484153450630075, 0.015770725590487672, 0.37335927045971884, 0.16378819895160301, 0.21351350664719881, 0.003260852296832066, 0.01025154979001902, 0.16665765312825145]), +outP{Float64}(0.01, 700.0, 0, -846.0406064189539, ["ilmm", "opx", "sp", "fsp", "fsp", "aug", "q", "H2O"], [0.03708129220778433, 0.16767880828486417, 0.014969371543223252, 0.014898318489368803, 0.3720456884664841, 0.21930629498743612, 0.007362082620733023, 0.1666581434001063]), +outP{Float64}(0.01, 800.0, 0, -859.7582656171719, ["opx", "ilmm", "sp", "fsp", "fsp", "aug", "q", "H2O"], [0.16342885743223898, 0.03549059291592639, 0.015137866284388204, 0.37105522757422493, 0.01308898266751593, 0.22702279394047314, 0.008121452184748955, 0.1666542270004835]), +outP{Float64}(0.01, 900.0, 0, -874.1236894100732, ["fsp", "sp", "opx", "ilmm", "liq", "aug", "H2O"], [0.3659710753192874, 0.01693795423706182, 0.1560291524582112, 0.03243288080729109, 0.026045652172651063, 0.23609392490062148, 0.16648936010487597]), +outP{Float64}(0.01, 1000.0, 0, -889.1058056671181, ["liq", "ilmm", "sp", "aug", "fsp", "opx", "H2O"], [0.06893894229577896, 0.03072631379799995, 0.016167542972787177, 0.24611990007455636, 0.3358207096260302, 0.1360081798231298, 0.16621841140971755]), +outP{Float64}(0.01, 1100.0, 0, -904.7854364309327, ["liq", "ilmm", "opx", "fsp", "aug", "H2O"], [0.2702373823245891, 0.037858535432415054, 0.058905457923985625, 0.21934146100616042, 0.24859251548400715, 0.16506464782884273]), +outP{Float64}(0.01, 1200.0, 0, -921.2541454561142, ["liq", "ilmm", "aug", "fsp", "H2O"], [0.46496624260301306, 0.03450737808217542, 0.22768144046900754, 0.10887939262845818, 0.16396554621734583]), +outP{Float64}(4.01, 600.0, 0, -820.2199888527514, ["hb", "fsp", "bi", "q", "sph", "H2O"], [0.5742497795441798, 0.2046694301441336, 0.014777323363389534, 0.05645353242716479, 0.021161443967668725, 0.1286884905534636]), +outP{Float64}(4.01, 700.0, 0, -831.923110725527, ["fsp", "aug", "hb", "liq", "q", "sph", "H2O"], [0.1961461570717106, 0.006685575916907892, 0.5785323337738232, 0.030711649084593243, 0.04704457540101029, 0.01732699571494817, 0.12355271303700663]), +outP{Float64}(4.01, 800.0, 0, -844.4883857670002, ["ilmm", "hb", "aug", "fsp", "liq", "H2O"], [0.01459620430783435, 0.4753621246832736, 0.08233246118330523, 0.13850254408606566, 0.2008751657923221, 0.08833149994719905]), +outP{Float64}(4.01, 900.0, 0, -857.8908410287308, ["hb", "liq", "aug", "ilmm", "fsp", "H2O"], [0.376874173140531, 0.30795523971089034, 0.14652135649609888, 0.015551542955303496, 0.08304562182680127, 0.0700520658703751]), +outP{Float64}(4.01, 1000.0, 0, -872.167579999219, ["ilmm", "aug", "liq", "fsp", "hb", "H2O"], [0.022138827144959013, 0.23168555090999876, 0.5300887128757286, 0.009481977574912707, 0.17032756367323001, 0.036277367821170935]), +outP{Float64}(4.01, 1100.0, 0, -887.4423558737342, ["liq", "ilmm", "aug", "H2O"], [0.720422846319944, 0.034632345979048004, 0.23733346134947747, 0.007611346351530607]), +outP{Float64}(4.01, 1200.0, 0, -903.5014188793097, ["liq", "aug", "ilmm"], [0.8224430231607827, 0.1357187483310143, 0.04183822850820312]), +outP{Float64}(8.01, 600.0, 0, -812.2017692250909, ["hb", "fsp", "ep", "bi", "aug", "q", "sph", "H2O"], [0.5833233007044378, 0.09055360632131977, 0.05626246324509833, 0.01392719408480675, 0.0192602818039003, 0.09236361999750765, 0.020752162477306567, 0.1235573713656228]), +outP{Float64}(8.01, 700.0, 0, -823.7713957985247, ["fsp", "liq", "hb", "aug", "q", "sph", "H2O"], [0.0994333104205156, 0.09160336513389841, 0.5732461070199772, 0.05233773518160196, 0.06729283514133143, 0.01717209810728685, 0.09891454899538855]), +outP{Float64}(8.01, 800.0, 0, -836.2695617502595, ["hb", "fsp", "liq", "aug", "sph", "H2O"], [0.5273133107054441, 0.013510164213048199, 0.33134693296502127, 0.09152785747971635, 0.0119117467966185, 0.024389987840151575]), +outP{Float64}(8.01, 900.0, 0, -849.6305928640046, ["ilmm", "liq", "aug", "hb", "H2O"], [0.013301466966644615, 0.4291522426911397, 0.16625116138574408, 0.39071563847380025, 0.0005794904826712127]), +outP{Float64}(8.01, 1000.0, 0, -863.8404659821086, ["aug", "hb", "liq", "ilmm"], [0.23469587880430898, 0.19869983758351675, 0.5452114717549424, 0.02139281185723177]), +outP{Float64}(8.01, 1100.0, 0, -878.950796834842, ["opx", "liq", "aug", "ilmm"], [0.0036391717331083753, 0.6883693848492245, 0.2752830443273452, 0.03270839909032193]), +outP{Float64}(8.01, 1200.0, 0, -894.8303527407421, ["liq", "ilmm", "aug"], [0.7672642760468534, 0.037406534108737365, 0.19532918984440917]), +outP{Float64}(12.01, 600.0, 0, -804.5540942528808, ["aug", "hb", "bi", "ep", "mu", "q", "sph", "H2O"], [0.029803289731446077, 0.5947723215264056, 0.007752169836009209, 0.11532147572015948, 0.010637015388421663, 0.10382406390905685, 0.02063846179342471, 0.11725120209507647]), +outP{Float64}(12.01, 700.0, 0, -815.971181250509, ["ep", "liq", "hb", "aug", "q", "sph", "H2O"], [0.06669091849180528, 0.14574760437219084, 0.5844970552419615, 0.043750390540647534, 0.07782404282609207, 0.017160812131234132, 0.06432917639606878]), +outP{Float64}(12.01, 800.0, 0, -828.3529905689918, ["hb", "liq", "g", "aug", "ru", "sph"], [0.4702262277894174, 0.3642220378384084, 0.03714520689865603, 0.12264535276043545, 0.004530587837280299, 0.0012305868758022511]), +outP{Float64}(12.01, 900.0, 0, -841.6235501762127, ["hb", "aug", "g", "liq", "ilmm"], [0.36554837252118566, 0.17834226931595867, 0.025735554642195804, 0.4162138041489833, 0.014159999371676547]), +outP{Float64}(12.01, 1000.0, 0, -855.7084890378607, ["liq", "hb", "g", "aug", "ilmm"], [0.5092857497183438, 0.21233708961308478, 0.014777108409380587, 0.2434684725551036, 0.02013157970408718]), +outP{Float64}(12.01, 1100.0, 0, -870.6508151696057, ["ilmm", "liq", "aug", "opx"], [0.03084052909940444, 0.6481666462580145, 0.30233324510548754, 0.018659579537093626]), +outP{Float64}(12.01, 1200.0, 0, -886.3698889829215, ["ilmm", "liq", "aug"], [0.03452176178051468, 0.7240213060122548, 0.24145693220723055]), +outP{Float64}(16.01, 600.0, 0, -797.1277772705217, ["ep", "mu", "hb", "aug", "q", "ru", "H2O"], [0.18823515061144205, 0.019877128981641266, 0.5881311157028573, 0.01616303259043474, 0.06798853264118487, 0.007314741689813592, 0.11229029778262627]), +outP{Float64}(16.01, 700.0, 0, -808.3898630481184, ["hb", "g", "ep", "liq", "aug", "q", "ru", "H2O"], [0.5174513022647408, 0.03875691853234428, 0.10022543413803552, 0.08118807037074996, 0.0802530804691145, 0.08719914473681188, 0.006280960577713748, 0.08864508891048926]), +outP{Float64}(16.01, 800.0, 0, -820.6592205848705, ["g", "liq", "hb", "aug", "ru"], [0.2292222840841752, 0.36469581453317773, 0.2868841929986683, 0.11230419835471624, 0.006893510029262677]), +outP{Float64}(16.01, 900.0, 0, -833.8343753044085, ["liq", "g", "hb", "aug", "ru"], [0.40628953998788614, 0.2606885685475112, 0.15643183147160564, 0.16935399319867386, 0.007236066794323157]), +outP{Float64}(16.01, 1000.0, 0, -847.7842083555013, ["aug", "liq", "g", "ilmm", "hb", "ru"], [0.24266307464499678, 0.47739194452679345, 0.22105678082771288, 0.018737539887073478, 0.038335377289833536, 0.001815282823589933]), +outP{Float64}(16.01, 1100.0, 0, -862.5372488911568, ["aug", "g", "liq", "ilmm"], [0.28204089088990003, 0.10989585050138502, 0.5808599276180125, 0.027203330990702414]), +outP{Float64}(16.01, 1200.0, 0, -878.079824611628, ["ilmm", "aug", "liq"], [0.032558406701846586, 0.27904829299060974, 0.6883933003075438]), +outP{Float64}(20.01, 600.0, 0, -789.912850641355, ["g", "ep", "mu", "aug", "hb", "hb", "q", "ru", "H2O"], [0.12199036528361312, 0.2267987816429559, 0.024362976736896014, 0.0019307004420007902, 0.08237069406061129, 0.3817676499134884, 0.03567144968286964, 0.008181717400929513, 0.11692566483663527]), +outP{Float64}(20.01, 700.0, 0, -801.0558507332044, ["ep", "g", "hb", "mu", "aug", "q", "ru", "H2O"], [0.10581850482058788, 0.1550183035110909, 0.34380647647325774, 0.022211770186985703, 0.1454863307713818, 0.08689486392751937, 0.007479046701772682, 0.13328470360740396]), +outP{Float64}(20.01, 800.0, 0, -813.1903790226486, ["hb", "liq", "ep", "aug", "g", "q", "ru"], [0.08432020977215587, 0.3606941781077958, 0.0071419997423541924, 0.1867873305612864, 0.3517001676178172, 0.0010804707548621997, 0.008275643443728404]), +outP{Float64}(20.01, 900.0, 0, -826.2901692828718, ["g", "aug", "liq", "ru"], [0.3889594862712292, 0.2026317888286557, 0.3996587645752772, 0.008749960324838]), +outP{Float64}(20.01, 1000.0, 0, -840.102362969159, ["liq", "g", "aug", "ru"], [0.4371992530912563, 0.350438828531181, 0.20361193414920958, 0.00874998422835318]), +outP{Float64}(20.01, 1100.0, 0, -854.6323978193435, ["aug", "ilmm", "g", "liq", "ru"], [0.23294040015641737, 0.01238426751962214, 0.23417453285188256, 0.5159201508727466, 0.004580648599331417]), +outP{Float64}(20.01, 1200.0, 0, -869.9489254487681, ["liq", "ilmm", "aug", "g"], [0.6320359496660592, 0.028517920187198058, 0.2640745162050315, 0.07537161394171138]), +outP{Float64}(24.01, 600.0, 0, -782.8474293239113, ["g", "hb", "mu", "ep", "aug", "q", "ru", "H2O"], [0.18716943005526812, 0.3659982603871683, 0.024699776429178744, 0.18767325630787654, 0.05687189664291661, 0.043521429051465284, 0.008348986443059508, 0.1257169646830671]), +outP{Float64}(24.01, 700.0, 0, -793.9431957184438, ["mu", "g", "ep", "aug", "q", "ky", "ru", "H2O"], [0.025768493765234043, 0.2758241975602092, 0.010345809833136665, 0.41059605234018576, 0.06520057901762943, 0.04162439800236827, 0.008750241124358168, 0.1618902283568785]), +outP{Float64}(24.01, 800.0, 0, -805.9341092756522, ["liq", "g", "aug", "q", "ky", "ru", "H2O"], [0.21465509300363386, 0.35185238657639234, 0.3130326976693037, 0.02878042593880082, 0.021596911347035296, 0.008750066536698966, 0.0613324189281349]), +outP{Float64}(24.01, 900.0, 0, -818.8929677258574, ["liq", "g", "aug", "ru"], [0.3776729591288991, 0.3958801989756918, 0.21769693229048578, 0.008749909604923382]), +outP{Float64}(24.01, 1000.0, 0, -832.570512185278, ["liq", "g", "aug", "ru"], [0.4110452580869383, 0.3761989107678498, 0.2040058608679105, 0.008749970277301192]), +outP{Float64}(24.01, 1100.0, 0, -846.930223403727, ["g", "aug", "liq", "ru"], [0.3290033048392981, 0.19614114621343942, 0.4661055216293581, 0.008750027317904429]), +outP{Float64}(24.01, 1200.0, 0, -862.0085981911471, ["g", "liq", "ilmm", "aug", "ru"], [0.20743915896078716, 0.5641411968059892, 0.008168985046454117, 0.21413373412197087, 0.006116925064798702]) ]