Skip to content

Commit

Permalink
MAGEMin_C v1.7.1
Browse files Browse the repository at this point in the history 
- Corrected issue that was making the ultramafic extended database crash
- Corrected issue when computing graphite mol fraction for metapelite extended database
  • Loading branch information
@NicolasRiel
NicolasRiel authored Jan 14, 2025
1 parent 6dbe161 commit 9e3c552
Show file tree
Hide file tree
Showing 12 changed files with 1,490 additions and 1,260 deletions.
4 changes: 3 additions & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
.vscode
Manifest.toml
/src/*.o
/src/*.o
/src/TC_database/*.o
/src/SB_database/*.o
6 changes: 3 additions & 3 deletions 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.7.0"
version = "1.7.1"

[deps]
CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
Expand All @@ -18,7 +18,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
CEnum = "0.4"
CSV = "0.10.14"
DataFrames = "1.6.1"
MAGEMin_jll = "1.6.3 - 1.6.3"
NLopt_jll = "2.8.0 - 2.8.0"
MAGEMin_jll = "1.6.4 - 1.6.4"
NLopt_jll = "2.9.0 - 2.9.0"
ProgressMeter = "1"
julia = "1.6"
69 changes: 34 additions & 35 deletions gen/magemin_library.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,41 +168,40 @@ const nlopt_precond = Ptr{Cvoid}
NLOPT_GN_ORIG_DIRECT_L = 7
NLOPT_GD_STOGO = 8
NLOPT_GD_STOGO_RAND = 9
NLOPT_LD_LBFGS_NOCEDAL = 10
NLOPT_LD_LBFGS = 11
NLOPT_LN_PRAXIS = 12
NLOPT_LD_VAR1 = 13
NLOPT_LD_VAR2 = 14
NLOPT_LD_TNEWTON = 15
NLOPT_LD_TNEWTON_RESTART = 16
NLOPT_LD_TNEWTON_PRECOND = 17
NLOPT_LD_TNEWTON_PRECOND_RESTART = 18
NLOPT_GN_CRS2_LM = 19
NLOPT_GN_MLSL = 20
NLOPT_GD_MLSL = 21
NLOPT_GN_MLSL_LDS = 22
NLOPT_GD_MLSL_LDS = 23
NLOPT_LD_MMA = 24
NLOPT_LN_COBYLA = 25
NLOPT_LN_NEWUOA = 26
NLOPT_LN_NEWUOA_BOUND = 27
NLOPT_LN_NELDERMEAD = 28
NLOPT_LN_SBPLX = 29
NLOPT_LN_AUGLAG = 30
NLOPT_LD_AUGLAG = 31
NLOPT_LN_AUGLAG_EQ = 32
NLOPT_LD_AUGLAG_EQ = 33
NLOPT_LN_BOBYQA = 34
NLOPT_GN_ISRES = 35
NLOPT_AUGLAG = 36
NLOPT_AUGLAG_EQ = 37
NLOPT_G_MLSL = 38
NLOPT_G_MLSL_LDS = 39
NLOPT_LD_SLSQP = 40
NLOPT_LD_CCSAQ = 41
NLOPT_GN_ESCH = 42
NLOPT_GN_AGS = 43
NLOPT_NUM_ALGORITHMS = 44
NLOPT_LD_LBFGS = 10
NLOPT_LN_PRAXIS = 11
NLOPT_LD_VAR1 = 12
NLOPT_LD_VAR2 = 13
NLOPT_LD_TNEWTON = 14
NLOPT_LD_TNEWTON_RESTART = 15
NLOPT_LD_TNEWTON_PRECOND = 16
NLOPT_LD_TNEWTON_PRECOND_RESTART = 17
NLOPT_GN_CRS2_LM = 18
NLOPT_GN_MLSL = 19
NLOPT_GD_MLSL = 20
NLOPT_GN_MLSL_LDS = 21
NLOPT_GD_MLSL_LDS = 22
NLOPT_LD_MMA = 23
NLOPT_LN_COBYLA = 24
NLOPT_LN_NEWUOA = 25
NLOPT_LN_NEWUOA_BOUND = 26
NLOPT_LN_NELDERMEAD = 27
NLOPT_LN_SBPLX = 28
NLOPT_LN_AUGLAG = 29
NLOPT_LD_AUGLAG = 30
NLOPT_LN_AUGLAG_EQ = 31
NLOPT_LD_AUGLAG_EQ = 32
NLOPT_LN_BOBYQA = 33
NLOPT_GN_ISRES = 34
NLOPT_AUGLAG = 35
NLOPT_AUGLAG_EQ = 36
NLOPT_G_MLSL = 37
NLOPT_G_MLSL_LDS = 38
NLOPT_LD_SLSQP = 39
NLOPT_LD_CCSAQ = 40
NLOPT_GN_ESCH = 41
NLOPT_GN_AGS = 42
NLOPT_NUM_ALGORITHMS = 43
end

function nlopt_algorithm_name(a)
Expand Down
3 changes: 2 additions & 1 deletion julia/MAGEMin_wrappers.jl
View file

Large diffs are not rendered by default.

2,530 changes: 1,339 additions & 1,191 deletions src/TC_database/TC_init_database.c
View file

Large diffs are not rendered by default.

77 changes: 75 additions & 2 deletions src/TC_database/objective_functions.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4085,6 +4085,73 @@ void p2x_um_po(void *SS_ref_db, double eps){
}
}

/**
Endmember to xeos for fsp_mp
*/
void p2x_ume_pl4tr(void *SS_ref_db, double eps){
SS_ref *d = (SS_ref *) SS_ref_db;

d->iguess[1] = d->p[2];
d->iguess[0] = d->p[1];

for (int i = 0; i < d->n_xeos; i++){
if (d->iguess[i] < d->bounds[i][0]){
d->iguess[i] = d->bounds[i][0];
}
if (d->iguess[i] > d->bounds[i][1]){
d->iguess[i] = d->bounds[i][1];
}
}
}
/**
Endmember to xeos for amp
*/
void p2x_ume_amp(void *SS_ref_db, double eps){
SS_ref *d = (SS_ref *) SS_ref_db;

d->iguess[7] = d->p[10];
d->iguess[6] = d->p[8];
d->iguess[2] = d->iguess[6] + d->p[3];
d->iguess[3] = d->p[2] + d->p[9];
d->iguess[4] = d->p[9]/(d->p[2] + d->p[9]);
d->iguess[5] = d->iguess[3] + d->p[0] + d->p[1] + d->iguess[7];
d->iguess[1] = -0.5*d->iguess[3] + d->iguess[5] -d->iguess[6] -d->p[0] -d->iguess[7] + d->iguess[2];
d->iguess[0] = (5.0*d->iguess[5] + 5.0*d->p[4] - 2.0*d->p[5] + d->p[6] + 5.0*d->iguess[2] - 5.0)/(2.0*d->iguess[5] + 2.0*d->iguess[6] + 2.0*d->iguess[7] + 2.0*d->iguess[1] + 2.0*d->iguess[2] - 7.0);
d->iguess[8] = -0.4*d->iguess[5]*d->iguess[0] + 2.0*d->iguess[5] - 0.4*d->iguess[6]*d->iguess[0] + 2.0*d->p[4] - 0.4*d->p[5] + 1.2*d->p[6] - 0.4*d->iguess[7]*d->iguess[0] - 0.4*d->iguess[0]*d->iguess[1] - 0.4*d->iguess[0]*d->iguess[2] + 2.4*d->iguess[0] + 2.0*d->iguess[2] - 2.0;
d->iguess[9] = (-2.0*d->iguess[5]*d->iguess[0] + 5.0*d->iguess[5] + 5.0*d->p[4] + 3.0*d->p[6] - 2.0*d->iguess[0]*d->iguess[2] + 5.0*d->iguess[0] + 5.0*d->iguess[2] - 5.0)/(2.0*d->iguess[6] + 2.0*d->iguess[7] + 2.0*d->iguess[1] - 2.0);

for (int i = 0; i < d->n_xeos; i++){
if (d->iguess[i] < d->bounds[i][0]){
d->iguess[i] = d->bounds[i][0];
}
if (d->iguess[i] > d->bounds[i][1]){
d->iguess[i] = d->bounds[i][1];
}
}
}
/**
Endmember to xeos for aug
*/
void p2x_ume_aug(void *SS_ref_db, double eps){
SS_ref *d = (SS_ref *) SS_ref_db;

d->iguess[6] = d->p[5];
d->iguess[1] = d->p[6] + d->iguess[6];
d->iguess[2] = d->p[4];
d->iguess[4] = d->iguess[2] + d->p[3];
d->iguess[3] = d->p[0] + d->iguess[1];
d->iguess[0] = (2.0*d->iguess[4] + 2.0*d->p[1] + d->p[7] + 2.0*d->iguess[3] - 2.0)/(2.0*d->iguess[4] + d->iguess[1] + d->iguess[3] - 2.0);
d->iguess[5] = (4.0*d->iguess[4]*d->p[1] + 4.0*d->iguess[4]*d->p[2] + 2.0*d->iguess[4]*d->p[7] + 4.0*d->iguess[4]*d->iguess[1] + 4.0*d->iguess[4]*d->iguess[3] - 8.0*d->iguess[4] + 4.0*d->iguess[4]*d->iguess[4] + 4.0*d->p[1]*d->iguess[1] - 4.0*d->p[1] + 2.0*d->p[2]*d->iguess[1] + 2.0*d->p[2]*d->iguess[3] - 4.0*d->p[2] + 2.0*d->p[7]*d->iguess[1] - 2.0*d->p[7] + 4.0*d->iguess[1]*d->iguess[3] - 4.0*d->iguess[1] - 4.0*d->iguess[3] + 4.0)/(d->iguess[4]*d->iguess[1] + 3.0*d->iguess[4]*d->iguess[3] - 4.0*d->iguess[4] + 2.0*d->iguess[4]*d->iguess[4] + d->iguess[1]*d->iguess[3] -d->iguess[1] - 3.0*d->iguess[3] + d->iguess[3]*d->iguess[3] + 2.0);

for (int i = 0; i < d->n_xeos; i++){
if (d->iguess[i] < d->bounds[i][0]){
d->iguess[i] = d->bounds[i][0];
}
if (d->iguess[i] > d->bounds[i][1]){
d->iguess[i] = d->bounds[i][1];
}
}
}

/**************************************************************************************/
/**************************************************************************************/
Expand Down Expand Up @@ -14870,6 +14937,12 @@ void TC_um_P2X_init( P2X_type *P2X_read,
P2X_read[iss] = p2x_um_opx; }
else if (strcmp( gv.SS_list[iss], "po") == 0){
P2X_read[iss] = p2x_um_po; }
else if (strcmp( gv.SS_list[iss], "pl4tr") == 0){
P2X_read[iss] = p2x_ume_pl4tr; }
else if (strcmp( gv.SS_list[iss], "amp") == 0){
P2X_read[iss] = p2x_ume_amp; }
else if (strcmp( gv.SS_list[iss], "aug") == 0){
P2X_read[iss] = p2x_ume_aug; }
else{
printf("\nsolid solution '%s' is not in the database, cannot be initiated\n", gv.SS_list[iss]);
}
Expand All @@ -14896,11 +14969,11 @@ void TC_P2X_init( P2X_type *P2X_read,
TC_igad_P2X_init( P2X_read,
gv );
}
else if (gv.EM_database == 4){ // ultramafic database //
else if (gv.EM_database == 4 || gv.EM_database == 5){ // ultramafic database //
TC_um_P2X_init( P2X_read,
gv );
}
else if (gv.EM_database == 7){ // ultramafic database //
else if (gv.EM_database == 7){ // metapelite ext database //
TC_mpe_P2X_init( P2X_read,
gv );
}
Expand Down
2 changes: 1 addition & 1 deletion src/initialize.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,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.6.3 [05/01/2025]"); /** MAGEMin version */
strcpy(gv.version,"1.6.4 [13/01/2025]"); /** MAGEMin version */

/* generate parameters */
strcpy(gv.buffer,"none");
Expand Down
5 changes: 2 additions & 3 deletions src/pp_min_function.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -703,7 +703,6 @@ global_variable init_em_db( int EM_database,
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 A C M F K N T O Mn H\n");
}
Expand All @@ -723,10 +722,10 @@ global_variable init_em_db( int EM_database,
printf(" S A M F O H S C N\n");
}
else if (EM_database == 6){
printf(" S A M F O H S C N\n");
printf(" S A M F O H S C N \n");
}
else if (EM_database == 7){
printf(" S A C M F N \n");
printf(" S A C M F K N T O Mn H CO2 S\n");
}
for (int j = 0; j < gv.len_ox; j++){
printf(" %.1f",PP_ref_db[i].Comp[j]);
Expand Down
16 changes: 10 additions & 6 deletions src/simplex_levelling.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -780,10 +780,12 @@ global_variable update_global_info( bulk_info z_b,
get initial conditions for active phases
*/
for (i = 0; i < d->n_Ox; i++){

ph_id = d->ph_id_A[i][1];

/* if phase is a pure species */
if (d->ph_id_A[i][0] == 1 ){

gv.pp_flags[ph_id][1] = 1;
gv.pp_flags[ph_id][2] = 0;
gv.pp_n[ph_id] = d->n_vec[i];
Expand All @@ -792,6 +794,7 @@ global_variable update_global_info( bulk_info z_b,
}
/* pure endmembers as solution phase */
if (d->ph_id_A[i][0] == 2){

phase_on[ph_id] = 1;
em_id = d->ph_id_A[i][3];

Expand All @@ -808,7 +811,7 @@ global_variable update_global_info( bulk_info z_b,
SS_ref_db[ph_id],
z_b,
ph_id );

strcpy(cp[id_cp].name,gv.SS_list[ph_id]); /* get phase name */

cp[id_cp].split = 0;
Expand Down Expand Up @@ -872,6 +875,7 @@ global_variable update_global_info( bulk_info z_b,

/* solution phase */
if (d->ph_id_A[i][0] == 3){

pc_id = d->ph_id_A[i][3];
phase_on[ph_id] = 1;

Expand Down Expand Up @@ -1437,17 +1441,17 @@ global_variable run_initial_guess_function( bulk_info z_b,
gv,
PP_ref_db,
SS_ref_db );

/* update global variable gamma */
update_global_gamma_LU( z_b,
splx_data );


/* remove solution from consideration when min driving force is > gv.bnd_filter_pc */
reduce_ss_list( SS_ref_db,
gv );


/* function to send back the updated initial guess, and phases flags */
gv = update_global_info( z_b,
splx_data,
Expand Down Expand Up @@ -1577,7 +1581,7 @@ global_variable run_levelling_function( bulk_info z_b,
reduce_ss_list( SS_ref_db,
gv );


/* function to send back the updated initial guess, and phases flags */
gv = update_global_info( z_b,
splx_data,
Expand Down
16 changes: 9 additions & 7 deletions src/toolkit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1375,7 +1375,9 @@ global_variable compute_phase_mol_fraction( global_variable gv,

sum = 0.0;
for (int j = 0; j < gv.len_ox; j++){
sum += PP_ref_db[i].Comp[j];
if (PP_ref_db[i].Comp[j] > 0.0 ){
sum += PP_ref_db[i].Comp[j];
}
}

n_at_ph = 0.0;
Expand Down Expand Up @@ -1500,7 +1502,7 @@ global_variable compute_density_volume_modulus( int EM_database,
// cp[i].volume_P0 += (dGdP_P0)*cp[i].p_em[j];

/* entropy */
cp[i].phase_entropy += -(dGdTMP)*cp[i].p_em[j];
cp[i].phase_entropy += -(dGdTMP)*cp[i].p_em[j]*cp[i].factor;

/* expansivity */
cp[i].phase_expansivity += (1.0/(dGdP)*((dGdTPP-dGdTMP)/(gv.gb_P_eps)))*cp[i].p_em[j];
Expand Down Expand Up @@ -1528,7 +1530,7 @@ global_variable compute_density_volume_modulus( int EM_database,
}

/* enthalpy */
cp[i].phase_enthalpy = cp[i].phase_entropy*T + G;
cp[i].phase_enthalpy = cp[i].phase_entropy*T + G*cp[i].factor;

/** calculate density from volume */
cp[i].phase_density = (cp[i].mass*1000.0)/(cp[i].volume*10.0);
Expand Down Expand Up @@ -1622,10 +1624,10 @@ global_variable compute_density_volume_modulus( int EM_database,
PP_ref_db[i].phase_expansivity = 1.0/(dGdP)*((dGdTPP-dGdTMP)/(gv.gb_P_eps));

/* entropy */
PP_ref_db[i].phase_entropy = -dGdTMP;
PP_ref_db[i].phase_entropy = -dGdTMP*PP_ref_db[i].factor;

/* enthalpy */
PP_ref_db[i].phase_enthalpy = PP_ref_db[i].phase_entropy*T + PP_ref_db[i].gbase;
PP_ref_db[i].phase_enthalpy = PP_ref_db[i].phase_entropy*T + PP_ref_db[i].gbase*PP_ref_db[i].factor;;

/* bulk modulus */
if ( strcmp(gv.research_group, "sb") == 0 ){
Expand Down Expand Up @@ -1704,7 +1706,7 @@ global_variable compute_density_volume_modulus( int EM_database,
for (int i = 0; i < gv.len_cp; i++){
if (cp[i].ss_flags[1] == 1){
gv.system_density += cp[i].phase_density*cp[i].ss_n_wt;
gv.system_entropy += cp[i].phase_entropy*cp[i].ss_n_mol*cp[i].factor;
gv.system_entropy += cp[i].phase_entropy*cp[i].ss_n_mol;//*cp[i].factor;
gv.system_cp += cp[i].phase_cp*cp[i].ss_n_mol;
gv.system_expansivity += cp[i].phase_expansivity*cp[i].ss_n_mol;
if (strcmp( cp[i].name, "liq") != 0 && strcmp( cp[i].name, "fl") != 0){
Expand All @@ -1715,7 +1717,7 @@ global_variable compute_density_volume_modulus( int EM_database,
for (int i = 0; i < gv.len_pp; i++){
if (gv.pp_flags[i][1] == 1 && gv.pp_flags[i][4] == 0){
gv.system_density += PP_ref_db[i].phase_density*gv.pp_n_wt[i];
gv.system_entropy += PP_ref_db[i].phase_entropy*gv.pp_n_mol[i]*PP_ref_db[i].factor;
gv.system_entropy += PP_ref_db[i].phase_entropy*gv.pp_n_mol[i];//*PP_ref_db[i].factor;
gv.system_cp += PP_ref_db[i].phase_cp*gv.pp_n_mol[i];
gv.system_expansivity += PP_ref_db[i].phase_expansivity*gv.pp_n_mol[i];
if (strcmp( gv.PP_list[i], "H2O") != 0){
Expand Down
18 changes: 9 additions & 9 deletions test/tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,15 @@ using Test
using MAGEMin_C # load MAGEMin (needs to be loaded from main directory to pick up correct library in case it is locally compiled)


# generic test for sb11 database
data = Initialize_MAGEMin("sb11", verbose=true);
test = 1 #KLB1
data = use_predefined_bulk_rock(data, test);
P = 80.0
T = 800.0
out = point_wise_minimization(P,T, data);
@test sort(out.ph) == sort(["gtmj", "hpcpx", "ol" ,"cpx"])
Finalize_MAGEMin(data)
# # generic test for sb11 database
# data = Initialize_MAGEMin("sb11", verbose=true);
# test = 1 #KLB1
# data = use_predefined_bulk_rock(data, test);
# P = 80.0
# T = 800.0
# out = point_wise_minimization(P,T, data);
# @test sort(out.ph) == sort(["gtmj", "hpcpx", "ol" ,"cpx"])
# Finalize_MAGEMin(data)

# generic test for thermocalc database
data = Initialize_MAGEMin("ig", verbose=true);
Expand Down
Loading

0 comments on commit 9e3c552

Please sign in to comment.