Skip to content

Commit

Permalink
MAGEMin_C v1.6.6 (#64)
Browse files Browse the repository at this point in the history 
- Added “name_solvus” option in order to differentiate solvus minerals
- Added bulk conversion functions mol2wt and wt2mol
  • Loading branch information
@NicolasRiel
NicolasRiel authored Dec 20, 2024
1 parent 8210991 commit 28c782e
Show file tree
Hide file tree
Showing 2 changed files with 204 additions and 16 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "MAGEMin_C"
uuid = "e5d170eb-415a-4524-987b-12f1bce1ddab"
authors = ["Nicolas Riel <[email protected]> & Boris Kaus <[email protected]>"]
version = "1.6.5"
version = "1.6.6"

[deps]
CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
Expand Down
218 changes: 203 additions & 15 deletions julia/MAGEMin_wrappers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,15 @@ using DataFrames, Dates, CSV

const VecOrMat = Union{Nothing, AbstractVector{Float64}, AbstractVector{<:AbstractVector{Float64}}}

export anhydrous_renormalization, retrieve_solution_phase_information, remove_phases,
export anhydrous_renormalization, retrieve_solution_phase_information, remove_phases, get_ss_from_mineral,
init_MAGEMin, finalize_MAGEMin, point_wise_minimization, convertBulk4MAGEMin, use_predefined_bulk_rock, define_bulk_rock, create_output,
print_info, create_gmin_struct, pwm_init, pwm_run,
single_point_minimization, multi_point_minimization, AMR_minimization, MAGEMin_Data, W_Data,
MAGEMin_data2dataframe, MAGEMin_dataTE2dataframe,
Initialize_MAGEMin, Finalize_MAGEMin

export wt2mol, mol2wt

export adjust_chemical_system, TE_prediction, get_OL_KDs_database, adjust_bulk_4_zircon

export initialize_AMR, split_and_keep, AMR
Expand Down Expand Up @@ -415,6 +417,7 @@ function single_point_minimization( P :: T1,
T :: T1,
MAGEMin_db :: MAGEMin_Data;
light :: Bool = false,
name_solvus :: Bool = false,
test :: Int64 = 0, # if using a build-in test case,
X :: VecOrMat = nothing,
B :: Union{Nothing, T1, Vector{T1}} = nothing,
Expand Down Expand Up @@ -442,6 +445,7 @@ function single_point_minimization( P :: T1,
T,
MAGEMin_db;
light = light,
name_solvus = name_solvus,
test = test,
X = X,
B = B,
Expand Down Expand Up @@ -525,6 +529,7 @@ function multi_point_minimization(P :: T2,
T :: T2,
MAGEMin_db :: MAGEMin_Data;
light :: Bool = false,
name_solvus :: Bool = false,
test :: Int64 = 0, # if using a build-in test case,
X :: VecOrMat = nothing,
B :: Union{Nothing, T1, Vector{T1}} = nothing,
Expand Down Expand Up @@ -599,9 +604,9 @@ function multi_point_minimization(P :: T2,
end
else
if isnothing(B)
out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data; scp, rm_list)
out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data; scp, rm_list, name_solvus=name_solvus)
else
out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data; buffer_n = B[i], W = W, scp, rm_list)
out = point_wise_minimization(P[i], T[i], gv, z_b, DB, splx_data; buffer_n = B[i], W = W, scp, rm_list, name_solvus=name_solvus)
end
end
elseif light == true
Expand Down Expand Up @@ -783,13 +788,54 @@ function normalize(vector::Vector{Float64})
end


function wt2mol( bulk_wt :: Vector{Float64},
bulk_ox :: Vector{String})

ref_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "O"; "Cr2O3"; "MnO"; "H2O"; "CO2"; "S"];
ref_MolarMass = [60.08; 101.96; 56.08; 40.30; 71.85; 159.69; 94.2; 61.98; 79.88; 16.0; 151.99; 70.937; 18.015; 44.01; 32.06]; #Molar mass of oxides

bulk_mol = zeros(length(bulk_ox));
bulk_wt = normalize(bulk_wt)

for i = 1:length(bulk_ox)
id = findfirst(ref_ox .== bulk_ox[i]);
bulk_mol[i] = bulk_wt[i]/ref_MolarMass[id];
end

bulk_mol .= bulk_mol ./sum(bulk_mol) .* 100.0

return bulk_mol
end

function mol2wt( bulk_mol :: Vector{Float64},
bulk_ox :: Vector{String})

ref_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "O"; "Cr2O3"; "MnO"; "H2O"; "CO2"; "S"];
ref_MolarMass = [60.08; 101.96; 56.08; 40.30; 71.85; 159.69; 94.2; 61.98; 79.88; 16.0; 151.99; 70.937; 18.015; 44.01; 32.06]; #Molar mass of oxides

bulk_wt = zeros(length(bulk_ox));
bulk_mol = normalize(bulk_mol)

for i = 1:length(bulk_ox)
id = findfirst(ref_ox .== bulk_ox[i]);
bulk_wt[i] = bulk_mol[i]*ref_MolarMass[id];
end

bulk_wt .= bulk_wt ./sum(bulk_wt) .* 100.0

return bulk_wt
end

"""
convertBulk4MAGEMin( bulk_in, bulk_in_ox, sys_in)
MAGEMin_bulk, MAGEMin_ox; = convertBulk4MAGEMin(bulk_in::T1,bulk_in_ox::Vector{String},sys_in::String,db::String) where {T1 <: AbstractVector{Float64}}
receives bulk-rock composition in [mol,wt] fraction and associated oxide list and sends back bulk-rock composition converted for MAGEMin use
"""
function convertBulk4MAGEMin(bulk_in::T1,bulk_in_ox::Vector{String},sys_in::String,db::String) where {T1 <: AbstractVector{Float64}}
function convertBulk4MAGEMin( bulk_in :: T1,
bulk_in_ox :: Vector{String},
sys_in :: String,
db :: String ) where {T1 <: AbstractVector{Float64}}

ref_ox = ["SiO2"; "Al2O3"; "CaO"; "MgO"; "FeO"; "Fe2O3"; "K2O"; "Na2O"; "TiO2"; "O"; "Cr2O3"; "MnO"; "H2O"; "CO2"; "S"];
ref_MolarMass = [60.08; 101.96; 56.08; 40.30; 71.85; 159.69; 94.2; 61.98; 79.88; 16.0; 151.99; 70.937; 18.015; 44.01; 32.06]; #Molar mass of oxides
Expand Down Expand Up @@ -985,6 +1031,7 @@ function point_wise_minimization( P ::Float64,
DB,
splx_data;
light = false,
name_solvus = false,
buffer_n = 0.0,
scp = 0,
rm_list = nothing,
Expand Down Expand Up @@ -1057,7 +1104,7 @@ function point_wise_minimization( P ::Float64,
if light == true
out = deepcopy(create_light_gmin_struct(DB));
else
out = deepcopy(create_gmin_struct(DB, gv, time));
out = deepcopy(create_gmin_struct(DB, gv, time; name_solvus = name_solvus));
end
# here we compute specific heat capacity using reactions
if (scp == 1)
Expand Down Expand Up @@ -1305,9 +1352,10 @@ point_wise_minimization(P :: Number,
buffer_n:: Float64 = 0.0,
scp :: Int64 = 0,
rm_list :: Union{Nothing, Vector{Int64}} = nothing,
name_solvus::Bool = false,
data_in :: Union{Nothing, gmin_struct{Float64, Int64}, Vector{gmin_struct{Float64, Int64}}} = nothing,
W :: Union{Nothing, W_Data} = nothing) =
point_wise_minimization(Float64(P),Float64(T), gv, z_b, DB, splx_data; buffer_n, scp, rm_list, data_in, W)
point_wise_minimization(Float64(P),Float64(T), gv, z_b, DB, splx_data; buffer_n, scp, rm_list, name_solvus, data_in, W)

point_wise_minimization(P :: Number,
T :: Number,
Expand All @@ -1319,19 +1367,21 @@ point_wise_minimization(P :: Number,
buffer_n:: Float64 = 0.0,
scp :: Int64 = 0,
rm_list :: Union{Nothing, Vector{Int64}} = nothing,
name_solvus::Bool = false,
data_in :: Union{Nothing, gmin_struct{Float64, Int64}, Vector{gmin_struct{Float64, Int64}}} = nothing,
W :: Union{Nothing, W_Data} = nothing) =
point_wise_minimization(Float64(P),Float64(T), gv, z_b, DB, splx_data; buffer_n, scp, rm_list, data_in, W)
point_wise_minimization(Float64(P),Float64(T), gv, z_b, DB, splx_data; buffer_n, scp, rm_list, name_solvus, data_in, W)

point_wise_minimization(P :: Number,
T :: Number,
data :: MAGEMin_Data;
buffer_n:: Float64 = 0.0,
scp :: Int64 = 0,
rm_list :: Union{Nothing, Vector{Int64}} = nothing,
name_solvus::Bool = false,
data_in :: Union{Nothing, gmin_struct{Float64, Int64}, Vector{gmin_struct{Float64, Int64}}} = nothing,
W :: Union{Nothing, W_Data} = nothing) =
point_wise_minimization(Float64(P),Float64(T), data.gv[1], data.z_b[1], data.DB[1], data.splx_data[1]; buffer_n, scp, rm_list, data_in, W)
point_wise_minimization(Float64(P),Float64(T), data.gv[1], data.z_b[1], data.DB[1], data.splx_data[1]; buffer_n, scp, rm_list, name_solvus, data_in, W)


"""
Expand Down Expand Up @@ -1390,12 +1440,145 @@ end

# pwm_run(gv, z_b, DB, splx_data) = pwm_run( gv, z_b, DB, splx_data)


function get_mineral_name(db, ss, SS_vec)

mineral_name = ss

if db == "ig" || db == "igad"
x = SS_vec.compVariables
if ss == "spn"
if x[3] - 0.5 > 0.0; mineral_name = "cm";
elseif x[2] - 0.5 > 0.0; mineral_name = "mt";
else mineral_name = "sp"; end
elseif ss == "fsp"
if x[2] - 0.5 > 0.0; mineral_name = "ksp";
else mineral_name = "pl"; end
elseif ss == "mu"
if x[4] - 0.5 > 0.0; mineral_name = "pat";
else mineral_name = "mu"; end
elseif ss == "hb"
if x[3] - 0.5 > 0.0; mineral_name = "gl";
elseif -x[3] -x[4] + 0.2 > 0.0; mineral_name = "act";
else
if x[6] < 0.1; mineral_name = "cum";
elseif -1/2*x[4]+x[6]-x[7]-x[8]-x[2]+x[3]>0.5; mineral_name = "tr";
else mineral_name = "hb"; end
end
elseif ss == "ilm"
if -x[1] + 0.5 > 0.0; mineral_name = "hem";
else mineral_name = "ilm"; end
elseif ss == "ness"
if x[2] - 0.5 > 0.0; mineral_name = "nesK";
else mineral_name = "ness"; end
elseif ss == "cpx"
if x[3] - 0.6 > 0.0; mineral_name = "pig";
elseif x[4] - 0.5 > 0.0; mineral_name = "omph";
else mineral_name = "cpx"; end
end

elseif db == "mp" || db == "mpe" || db == "mb" || db == "ume"
x = SS_vec.compVariables
if ss == "sp"
if x[2] - 0.5 > 0.0; mineral_name = "mt";
else mineral_name = "sp"; end
elseif ss == "fsp"
if x[2] - 0.5 > 0.0; mineral_name = "ksp";
else mineral_name = "pl"; end
elseif ss == "mu"
if x[4] - 0.5 > 0.0; mineral_name = "pat";
else mineral_name = "mu"; end
elseif ss == "hb"
if x[3] - 0.5 > 0.0; mineral_name = "gl";
elseif -x[3]-x[4]+0.2>0.0; mineral_name = "act";
else
if x[6] < 0.1; mineral_name = "cum";
elseif -1/2*x[4]+x[6]-x[7]-x[8]-x[2]+x[3]>0.5; mineral_name = "tr";
else mineral_name = "hb"; end
end
elseif ss == "ilmm"
if x[1] - 0.5 > 0.0; mineral_name = "ilmm";
else mineral_name = "hemm"; end
elseif ss == "ilm"
if 1.0 - x[1] > 0.5; mineral_name = "hem";
else mineral_name = "ilm"; end
elseif ss == "dio"
if 2*x[4] > 0.5; mineral_name = "omph";
else mineral_name = "dio"; end
elseif ss == "occm"
if x[2] > 0.5; mineral_name = "sid";
elseif x[3] > 0.5; mineral_name = "ank";
elseif x[1] > 0.25 && x[3] < 0.01; mineral_name = "mag";
else mineral_name = "cc"; end

end

end

return mineral_name
end

"""
This function returns the solution phase name given the mineral name (handling solvus -> solution phase)
"""
function get_ss_from_mineral(db, mrl, mbCpx)

ss = mrl

if db =="ig" || db == "igad"

if mrl == "cm" || mrl == "mt" || mrl == "sp"
ss = "spn"
elseif mrl == "pat" || mrl == "mu"
ss = "mu"
elseif mrl == "ksp" || mrl == "pl"
ss = "fsp"
elseif mrl == "gl" || mrl == "act" || mrl == "hb" || mrl == "cum" || mrl == "tr"
ss = "hb"
elseif mrl == "hem" || mrl == "ilm"
ss = "ilm"
elseif mrl == "pig" || mrl == "omph"
ss = "cpx"
elseif mrl == "nesK"
ss = "ness"
end

elseif db == "mp" || db == "mpe" || db == "mb" || db == "ume"

if mrl == "mt" || mrl == "sp"
ss = "sp"
elseif mrl == "ksp" || mrl == "pl"
ss = "fsp"
elseif mrl == "pat" || mrl == "mu"
ss = "mu"
elseif mrl == "gl" || mrl == "act" || mrl == "hb" || mrl == "cum" || mrl == "tr"
ss = "hb"
elseif mrl == "hem" || mrl == "ilm"
ss = "ilm"
elseif mrl == "hemm" || mrl == "ilmm"
ss = "ilmm"
elseif mrl == "omph" || mrl == "dio"
if mbCpx == 0
ss = "dio"
else
ss = "aug"
end
elseif mrl == "sid" || mrl == "mag" || mrl == "ank" || mrl == "cc"
ss = "occm"
end


end

return ss
end

"""
out = create_gmin_struct(gv, z_b)
This extracts the output of a pointwise MAGEMin optimization and adds it into a julia structure
"""
function create_gmin_struct(DB, gv, time)
function create_gmin_struct(DB, gv, time; name_solvus = false)

stb = unsafe_load(DB.sp)

Expand Down Expand Up @@ -1469,18 +1652,23 @@ function create_gmin_struct(DB, gv, time)
n_SS = stb.n_SS # number of solid solutions
n_mSS = stb.n_mSS # number of solid solutions

ph_frac = unsafe_wrap(Vector{Cdouble},stb.ph_frac, n_ph)
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)
ph_id_db = unsafe_wrap(Vector{Cint}, stb.ph_id_db , n_ph)
ph = unsafe_string.(unsafe_wrap(Vector{Ptr{Int8}}, stb.ph, n_ph)) # stable phases
ph_type = unsafe_wrap(Vector{Cint}, stb.ph_type, n_ph)
ph_id = unsafe_wrap(Vector{Cint}, stb.ph_id , n_ph)
ph_id_db = unsafe_wrap(Vector{Cint}, stb.ph_id_db , n_ph)
ph = unsafe_string.(unsafe_wrap(Vector{Ptr{Int8}}, stb.ph, n_ph)) # stable phases

# extract info about compositional variables of the solution models:
SS_vec = convert.(LibMAGEMin.SS_data, unsafe_wrap(Vector{LibMAGEMin.stb_SS_phase},stb.SS,n_SS))

if name_solvus == true
for i=1:n_SS
ph[i] = get_mineral_name(database, ph[i], SS_vec[i])
end
end
# extract information about metastable solution phases
mSS_vec = convert.(LibMAGEMin.mSS_data, unsafe_wrap(Vector{LibMAGEMin.mstb_SS_phase},stb.mSS,n_mSS))

Expand Down

0 comments on commit 28c782e

Please sign in to comment.