MAGEMin_C 1.6.5 (#63)

- Added frac_M_vol, frac_F_vol and frac_S_vol output for MAGEMin_C - Corrected wrong set of W's for spinel in the metapelite and metapelite extended database - Added molar mass of CO2 for MAGEMin_C (which was leading to failed mpe bulk loading in the App)
ComputationalThermodynamics · Dec 3, 2024 · 8210991 · 8210991
1 parent dac9941
commit 8210991
Show file tree

Hide file tree

Showing 17 changed files with 228 additions and 111 deletions.
diff --git a/Project.toml b/Project.toml
@@ -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.4"
+authors = ["Nicolas Riel <[email protected]>  & Boris Kaus <[email protected]>"]
+version = "1.6.5"
 
 [deps]
 CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
@@ -17,6 +17,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 CEnum = "0.4"
 CSV = "0.10.14"
 DataFrames = "1.6.1"
-MAGEMin_jll = "1.6.0 - 1.6.0"
+MAGEMin_jll = "1.6.1 - 1.6.1"
 ProgressMeter = "1"
 julia = "1.6"
diff --git a/gen/magemin_library.jl b/gen/magemin_library.jl
@@ -1177,10 +1177,13 @@ struct stb_systems
     rho_F::Cdouble
     bulk_S_wt::Ptr{Cdouble}
     frac_S_wt::Cdouble
+    frac_S_vol::Cdouble
     bulk_M_wt::Ptr{Cdouble}
     frac_M_wt::Cdouble
+    frac_M_vol::Cdouble
     bulk_F_wt::Ptr{Cdouble}
     frac_F_wt::Cdouble
+    frac_F_vol::Cdouble
     n_ph::Cint
     n_PP::Cint
     n_SS::Cint
@@ -1323,7 +1326,7 @@ mutable struct metapelite_datasets
     n_pp::Cint
     n_ss::Cint
     ox::NTuple{11, NTuple{20, Cchar}}
-    PP::NTuple{23, NTuple{20, Cchar}}
+    PP::NTuple{24, NTuple{20, Cchar}}
     SS::NTuple{18, NTuple{20, Cchar}}
     verifyPC::NTuple{18, Cint}
     n_SS_PC::NTuple{18, Cint}
@@ -1349,7 +1352,7 @@ mutable struct metabasite_datasets
     n_pp::Cint
     n_ss::Cint
     ox::NTuple{10, NTuple{20, Cchar}}
-    PP::NTuple{23, NTuple{20, Cchar}}
+    PP::NTuple{24, NTuple{20, Cchar}}
     SS::NTuple{17, NTuple{20, Cchar}}
     verifyPC::NTuple{17, Cint}
     n_SS_PC::NTuple{17, Cint}
@@ -1505,7 +1508,7 @@ mutable struct metapelite_datasets_ext
     n_pp::Cint
     n_ss::Cint
     ox::NTuple{13, NTuple{20, Cchar}}
-    PP::NTuple{25, NTuple{20, Cchar}}
+    PP::NTuple{26, NTuple{20, Cchar}}
     SS::NTuple{23, NTuple{20, Cchar}}
     verifyPC::NTuple{23, Cint}
     n_SS_PC::NTuple{23, Cint}

diff --git a/julia/MAGEMin_wrappers.jl b/julia/MAGEMin_wrappers.jl
diff --git a/src/MAGEMin.c b/src/MAGEMin.c
@@ -914,7 +914,7 @@ global_variable SetupDatabase(			global_variable 	 gv,
 
 	}
 
-	if (gv.verbose == 1){	
+	if (gv.verbose == 2){	
 		printf("\n");	
 		printf("--verbose     : verbose              = %i \n", 	 	   		gv.verbose			);
 		printf("--rg          : research group       = %s \n", 	 	  		gv.research_group	);

diff --git a/src/MAGEMin.h b/src/MAGEMin.h
@@ -769,10 +769,10 @@ typedef struct stb_systems {
 	double *bulk_M; double frac_M; double rho_M; 	/* Melt system informations 												*/
 	double *bulk_F; double frac_F; double rho_F; 	/* Fluid system informations 												*/
 
-	double *bulk_S_wt; double frac_S_wt;   			/* Solid system informations 												*/
-	double *bulk_M_wt; double frac_M_wt;  			/* Melt system informations 												*/
-	double *bulk_F_wt; double frac_F_wt;  			/* Fluid system informations 												*/
-	
+	double *bulk_S_wt; double frac_S_wt; double frac_S_vol;  			/* Solid system informations 												*/
+	double *bulk_M_wt; double frac_M_wt; double frac_M_vol;  			/* Melt system informations 												*/
+	double *bulk_F_wt; double frac_F_wt; double frac_F_vol; 			/* Fluid system informations 												*/
+
 	int     n_ph;									/* number of predicted stable phases 										*/
 	int     n_PP;									/* number of predicted stable pure phases 									*/
 	int     n_SS;									/* number of predicted stable solution phases 								*/

diff --git a/src/TC_database/NLopt_opt_function.c b/src/TC_database/NLopt_opt_function.c
@@ -7718,7 +7718,7 @@ SS_ref NLopt_opt_aq17_function(global_variable gv, SS_ref SS_ref_db){
        SS_ref_db.ub[i] = SS_ref_db.bounds[i][1];
     }
 
-    SS_ref_db.opt = nlopt_create(NLOPT_LN_COBYLA, (n)); 
+    SS_ref_db.opt = nlopt_create(NLOPT_LD_SLSQP, (n)); 
     nlopt_set_lower_bounds(SS_ref_db.opt, SS_ref_db.lb);
     nlopt_set_upper_bounds(SS_ref_db.opt, SS_ref_db.ub);
     nlopt_set_min_objective(SS_ref_db.opt, obj_aq17, &SS_ref_db);

diff --git a/src/TC_database/TC_gem_function.c b/src/TC_database/TC_gem_function.c
@@ -1035,11 +1035,11 @@ PP_ref G_FS_function(		int 		 len_ox,
 	PP_ref_db.factor  =  factor;
 	PP_ref_db.phase_shearModulus  = 0.0;
 
-	// printf(" %4s %+10f | factor: %+10f\n",name,PP_ref_db.gbase,PP_ref_db.factor);
+	// printf(" %4s %+10f %+10f | factor: %+10f\n",name,G,PP_ref_db.gbase,PP_ref_db.factor);
 	// for (i = 0; i < len_ox; i++){
 	// 	printf("%+10f",PP_ref_db.Comp[i]*PP_ref_db.factor); 
 	// }
-	// printf("\n");
+	// printf("\n\n");
 
 	return (PP_ref_db);
 }
diff --git a/src/TC_database/TC_init_database.c b/src/TC_database/TC_init_database.c
@@ -35,10 +35,10 @@ oxide_data oxide_info = {
 metapelite_dataset metapelite_db = {
 	62,							/* Endmember default dataset number */
 	11,							/* number of oxides */			
-	23,							/* number of pure phases */
+	24,							/* number of pure phases */
 	17,							/* number of solution phases */
 	{"SiO2"	,"Al2O3","CaO"	,"MgO"	,"FeO"	,"K2O"	,"Na2O"	,"TiO2"	,"O"	,"MnO"	,"H2O"											},
-	{"q"	,"crst"	,"trd"	,"coe"	,"stv"	,"ky"	,"sill"	,"and"	,"ru"	,"sph"	,"O2"  ,"H2O"	,
+	{"q"	,"crst"	,"trd"	,"coe"	,"stv"	,"ky"	,"sill"	,"and"	,"ru"	,"sph"	,"O2"  ,"H2O"	,"zo"	,
 	"qfm"	,"qif"	,"nno"	,"hm"	,"cco"	,"aH2O"	, "aO2"	,"aMgO"	,"aFeO"	,"aAl2O3"		,"aTiO2"								},
 	{"liq"	,"fsp"	,"bi"	,"g"	,"ep"	,"ma"	,"mu"	,"opx"	,"sa"	,"cd"	,"st"	,"chl"	,"ctd"	,"sp"  ,"mt"  ,"ilm"  ,"ilmm"  ,"aq17"	},
 
@@ -63,10 +63,10 @@ metapelite_dataset metapelite_db = {
 metabasite_dataset metabasite_db = {
 	62,							/* Endmember default dataset number */
 	10,							/* number of oxides */			
-	23,							/* number of pure phases */
+	24,							/* number of pure phases */
 	17,							/* number of solution phases */
 	{"SiO2"	,"Al2O3","CaO"	,"MgO"	,"FeO"	,"K2O"	,"Na2O"	,"TiO2"	,"O"	,"H2O"													},
-	{"q"	,"crst"	,"trd"	,"coe"	,"law"	,"ky"	,"sill"	,"and"	,"ru"	,"sph" 	,"ab"	,"H2O"	,
+	{"q"	,"crst"	,"trd"	,"coe"	,"law"	,"ky"	,"sill"	,"and"	,"ru"	,"sph" 	,"ab"	,"H2O"	,"zo"	,
 	"qfm"	,"qif"	,"nno"	,"hm"	,"cco"	,"aH2O"	,"aO2"	,"aMgO"	,"aFeO"	,"aAl2O3"		,"aTiO2"								},
 	{"sp"	,"opx"	,"fsp"	,"liq"	,"mu"	,"ilmm"	,"ilm"	,"ol"	,"hb"	,"ep"	,"g"	,"chl"	,"bi"	,"dio"	,"aug"	,"abc"  ,"spn"			},
 
@@ -232,10 +232,10 @@ mantle_dataset mantle_db = {
 metapelite_dataset_ext metapelite_ext_db = {
 	62,							/* Endmember default dataset number */
 	13,							/* number of oxides */			
-	25,							/* number of pure phases */
+	26,							/* number of pure phases */
 	23,							/* number of solution phases */
 	{"SiO2"	,"Al2O3","CaO"	,"MgO"	,"FeO"	,"K2O"	,"Na2O"	,"TiO2"	,"O"	,"MnO"	,"H2O"	,"CO2"	,"S"										},
-	{"q"	,"crst"	,"trd"	,"coe"	,"stv"	,"ky"	,"sill"	,"and"	,"ru"	,"sph"	,"O2" 	,"pyr"	,"gph"	,"law"	,
+	{"q"	,"crst"	,"trd"	,"coe"	,"stv"	,"ky"	,"sill"	,"and"	,"ru"	,"sph"	,"O2" 	,"pyr"	,"gph"	,"law"	,"zo"	,
 	"qfm"	,"qif"	,"nno"	,"hm"	,"cco"	,"aH2O"	, "aO2"	,"aMgO"	,"aFeO"	,"aAl2O3"		,"aTiO2"											},
 	{"liq"	,"fsp"	,"bi"	,"g"	,"ep"	,"ma"	,"mu"	,"opx"	,"sa"	,"cd"	,"st"	,"chl"	,"ctd"	,"sp"  ,"mt"  ,"ilm"  ,"ilmm"  ,"occm"	,"fl"	,"po"	,"dio"	,"aug"	,"hb"		},
 

diff --git a/src/TC_database/TC_init_database.h b/src/TC_database/TC_init_database.h
@@ -22,7 +22,7 @@
 		int 	n_pp;
 		int 	n_ss;
 		char    ox[11][20];
-		char    PP[23][20];
+		char    PP[24][20];
 		char    SS[18][20];
 
 		int 	verifyPC[18];
@@ -54,7 +54,7 @@
 		int 	n_pp;
 		int 	n_ss;
 		char    ox[10][20];
-		char    PP[23][20];
+		char    PP[24][20];
 
 		char    SS[17][20];
 		int 	verifyPC[17];
@@ -246,7 +246,7 @@
 		int 	n_pp;
 		int 	n_ss;
 		char    ox[13][20];
-		char    PP[25][20];
+		char    PP[26][20];
 		char    SS[23][20];
 
 		int 	verifyPC[23];

diff --git a/src/TC_database/aqueous.md b/src/TC_database/aqueous.md
@@ -0,0 +1,59 @@
+# Aqueous species formulation, using Miron et al., 2017 database.
+
+
+```math
+G_{aq} = \sum_{i=1}^{n^{em}} n_i \mu_i
+```
+
+where $n_i$ the molar fraction and $\mu_i$ is the chemical potential of the species $i$.
+
+```math
+\mu_i  = G_i^0 + RT log(\gamma_i m_i)
+```
+
+where $m_i$ is the molal fraction (!?)
+
+And then develops into:
+
+```math
+\mu_i  = G_i^0 + RT( log(\gamma_i) + log(m_i))
+```
+
+---
+
+For solute (aqueous species, excluding water), and assuming $\gamma_i$ = 1 (for now)
+
+```math
+\mu_i  = G_i^0 + RT log(m_i)
+```
+
+where 
+
+```math
+m_i  = \frac{n_i}{m_{H2O}}
+```
+
+so here $n_i$ is again the molar fraction and $m_{H2O}$ is mass of water? Should it be water density?
+
+I could find the following definition:
+
+```math
+m_i  = \frac{n_i}{n_{H2O} M_{H2O}}
+```
+and here $n_i$ is the molarity and $M_{H2O}$ is the Molar mass of water (18.015). However, here density of water like we discussed is not appearing, is that right?
+
+Should it simply be:
+
+```math
+m_i  = \frac{n_i}{\rho_{H2O}}
+```
+
+---
+
+For solvant (water):
+
+```math
+\mu_i  = G_i^0 + RT log(\frac{1}{M_{H2O}})
+```
+
+!?
diff --git a/src/TC_database/objective_functions.c b/src/TC_database/objective_functions.c
@@ -5645,7 +5645,6 @@ double obj_mp_sp(unsigned n, const double *x, double *grad, void *SS_ref_db){
     sf[3]          = 1.0 - x[0];
     sf[4]          = x[0];
 
-
     mu[0]          = R*T*creal(clog(sf[0]*sf[4])) + gb[0] + mu_Gex[0];
     mu[1]          = R*T*creal(clog(sf[0]*sf[3])) + gb[1] + mu_Gex[1];
     mu[2]          = R*T*creal(clog(sf[4]*sf[1] + d_em[2])) + gb[2] + mu_Gex[2];
@@ -14523,7 +14522,7 @@ double obj_aq17(unsigned n, const double *x, double *grad, void *SS_ref_db) {
                                 d->g,  
                                 d->epsilon,
                                 xiw                );  //fraction of water
-        // cor          = HSC_to_SUPCRT(ElH, Comp[i], len_ox);
+        cor          = HSC_to_SUPCRT(ElH, Comp[i], len_ox);
         mu[i]        = gb[i]  + (log(pow(10.0,loggamma)) + log(1000.0/18.0153) + log(x[i]/Xw) - log(xiw/Xw) - xiw/Xw + 1.0 )/1000.0;
         m_all       += x[i];
         if (charge[i] != 0.0){
@@ -14552,7 +14551,7 @@ double obj_aq17(unsigned n, const double *x, double *grad, void *SS_ref_db) {
                                         m_all       );
 
     /* set chemical potential of water */
-    // cor          = HSC_to_SUPCRT(ElH, Comp[0], len_ox);
+    cor          = HSC_to_SUPCRT(ElH, Comp[0], len_ox);
     mu[0] = gb[0]  + ( log(logawater) + log(xiw/Xw) - Xw/xiw - xiw/Xw + 2.0)/1000.0;
 
     px_aq17(SS_ref_db,x);

diff --git a/src/TC_database/tc_gss_function.c b/src/TC_database/tc_gss_function.c
@@ -4246,12 +4246,12 @@ SS_ref G_SS_mp_sp_function(SS_ref SS_ref_db, char* research_group, int EM_datase
     };
     //int si = sizeof(SF_tmp)/sizeof(SF_tmp[0]); printf("sz_tmp: %d n_sf %d\n",si,SS_ref_db.n_sf);
 
-    SS_ref_db.W[0] = 16.;
-    SS_ref_db.W[1] = 2.;
-    SS_ref_db.W[2] = 20.;
-    SS_ref_db.W[3] = 18.;
-    SS_ref_db.W[4] = 36.;
-    SS_ref_db.W[5] = 30.;
+    SS_ref_db.W[0] = 0.0;
+    SS_ref_db.W[1] = 18.5;
+    SS_ref_db.W[2] = 27.0;
+    SS_ref_db.W[3] = 40.0;
+    SS_ref_db.W[4] = 30.0;
+    SS_ref_db.W[5] = 0.0;
 
 
     em_data herc_eq 	= get_em_data(		research_group, EM_dataset, 
@@ -13401,12 +13401,12 @@ SS_ref G_SS_mpe_sp_function(SS_ref SS_ref_db, char* research_group, int EM_datas
     };
     //int si = sizeof(SF_tmp)/sizeof(SF_tmp[0]); printf("sz_tmp: %d n_sf %d\n",si,SS_ref_db.n_sf);
 
-    SS_ref_db.W[0] = 16.;
-    SS_ref_db.W[1] = 2.;
-    SS_ref_db.W[2] = 20.;
-    SS_ref_db.W[3] = 18.;
-    SS_ref_db.W[4] = 36.;
-    SS_ref_db.W[5] = 30.;
+    SS_ref_db.W[0] = 0.0;
+    SS_ref_db.W[1] = 18.5;
+    SS_ref_db.W[2] = 27.0;
+    SS_ref_db.W[3] = 40.0;
+    SS_ref_db.W[4] = 30.0;
+    SS_ref_db.W[5] = 0.0;
 
 
     em_data herc_eq 	= get_em_data(		research_group, EM_dataset, 

diff --git a/src/TC_database/tc_gss_init_function.c b/src/TC_database/tc_gss_init_function.c
@@ -1730,7 +1730,7 @@ SS_ref G_SS_mpe_opx_init_function(SS_ref SS_ref_db, global_variable gv){
     SS_ref_db.n_cat     = 0;
     SS_ref_db.is_liq    = 0;
     SS_ref_db.symmetry  = 0;
-    SS_ref_db.n_sf      = 10;
+    SS_ref_db.n_sf      = 11;
     SS_ref_db.n_em      = 7;
     SS_ref_db.n_v       = 7;
     SS_ref_db.n_w       = 21;

diff --git a/src/dump_function.c b/src/dump_function.c
@@ -173,6 +173,9 @@ void reset_output_struct(		global_variable 	 gv,
 	sp[0].frac_S_wt				 = 0.0;
 	sp[0].frac_M_wt				 = 0.0;
 	sp[0].frac_F_wt				 = 0.0;
+	sp[0].frac_S_vol			 = 0.0;
+	sp[0].frac_M_vol			 = 0.0;
+	sp[0].frac_F_vol			 = 0.0;
 	sp[0].rho_S				 	 = 0.0;
 	sp[0].rho_M				 	 = 0.0;
 	sp[0].rho_F				 	 = 0.0;
@@ -570,6 +573,7 @@ void fill_output_struct(		global_variable 	 gv,
 					if (gv.n_phase == 1){
 						sp[0].frac_M 				= 1.0;
 						sp[0].frac_M_wt				= 1.0;
+						sp[0].frac_M_vol		    = 1.0;
 					}
 					else{
 						sp[0].frac_M 				= cp[i].ss_n_mol;
@@ -686,6 +690,19 @@ void fill_output_struct(		global_variable 	 gv,
 	for (int i = 0; i < gv.len_cp; i++){
 		if ( cp[i].ss_flags[1] == 1){
 			sp[0].ph_frac_vol[n] = sp[0].ph_frac_wt[n] / sp[0].SS[n].rho;
+
+			if (strcmp( cp[i].name, "liq") == 0 || strcmp( cp[i].name, "fl") == 0 ){
+				if (strcmp( cp[i].name, "liq") == 0){
+						sp[0].frac_M_vol      = sp[0].ph_frac_vol[n];
+				}
+				else{
+						sp[0].frac_F_vol      = sp[0].ph_frac_vol[n];
+				}
+			}
+			else{
+				sp[0].frac_S_vol += sp[0].ph_frac_vol[n];
+			}
+
 			sum_vol += sp[0].ph_frac_vol[n];
 			sum_mol += sp[0].ph_frac[n];
 			sum_wt  += sp[0].ph_frac_wt[n];
@@ -696,6 +713,14 @@ void fill_output_struct(		global_variable 	 gv,
 	for (int i = 0; i < gv.len_pp; i++){
 		if (gv.pp_flags[i][1] == 1 && gv.pp_flags[i][4] == 0){
 			sp[0].ph_frac_vol[n] =  sp[0].ph_frac_wt[n] / sp[0].PP[m].rho;
+
+			if  (strcmp( gv.PP_list[i], "H2O") != 0){
+				sp[0].frac_S_vol += sp[0].ph_frac_vol[n];
+			}
+			if  (strcmp( gv.PP_list[i], "H2O") == 0){
+				sp[0].frac_F_vol = sp[0].ph_frac_vol[n];
+			}
+
 			sum_vol += sp[0].ph_frac_vol[n];
 			sum_mol += sp[0].ph_frac[n];
 			sum_wt  += sp[0].ph_frac_wt[n];
@@ -709,6 +734,7 @@ void fill_output_struct(		global_variable 	 gv,
 		sp[0].ph_frac_wt[i] 	/= sum_wt;
 	}
 
+
 	/* The following section normalizes the entries for S (solid), M (melt) and F (fluid) which are entries useful for geodynamic coupling */
 	// normalize rho_S and bulk_S
 	if (sp[0].frac_S > 0.0){
@@ -763,6 +789,12 @@ void fill_output_struct(		global_variable 	 gv,
 	sp[0].frac_M_wt /= sum;
 	sp[0].frac_S_wt /= sum;
 
+	sum = sp[0].frac_F_vol + sp[0].frac_M_vol + sp[0].frac_S_vol;
+
+	sp[0].frac_F_vol /= sum;
+	sp[0].frac_M_vol /= sum;
+	sp[0].frac_S_vol /= sum;
+
 	/* compute cp as J/K/kg for given bulk-rock composition */
 	sp[0].s_cp 					= sp[0].cp_wt/mass_bulk*1e6;
 	if (strcmp(gv.research_group, "tc") 	== 0 ){