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
ComputationalThermodynamics · Nov 7, 2024 · 7aea204 · 7aea204
1 parent ebd989f
commit 7aea204
Show file tree

Hide file tree

Showing 27 changed files with 1,652 additions and 183 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "MAGEMin_C"
 uuid = "e5d170eb-415a-4524-987b-12f1bce1ddab"
-authors = ["Boris Kaus <kaus@uni-mainz.de> & Nicolas Riel <riel@uni-mainz.de>"]
-version = "1.6.0"
+authors = ["Nicolas Riel <riel@uni-mainz.de> & Boris Kaus <kaus@uni-mainz.de>"]
+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
@@ -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,22 +2303,20 @@ 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}
     longidx::Cint
     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
diff --git 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
@@ -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
@@ -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
@@ -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); 
 	// }