Develop

.gitignore

@@ -1,3 +1,24 @@
+# Virtual environment and pip install requirements
+venv
+requirements.txt
+
+# Vscode
+.vscode/settings.json
+
+#Pytest
+.pytest_cache
+
+# Results directory
+results/run1
+results/
+
+# GRT directory (used for reference)
+source/codeYoung/GRT_1300_6000
+
+# temporary progress files
+progress.png 
+progress.txt
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

.vscode/settings.json

@@ -0,0 +1,7 @@
+{
+    "python.testing.pytestArgs": [
+        "tests"
+    ],
+    "python.testing.unittestEnabled": false,
+    "python.testing.pytestEnabled": true
+}

requirements.txt

@@ -0,0 +1,23 @@
+contourpy==1.3.0
+corner==2.2.2
+cycler==0.12.1
+emcee==3.1.6
+fonttools==4.54.1
+importlib_metadata==8.5.0
+importlib_resources==6.4.5
+jax==0.4.30
+jaxlib==0.4.30
+kiwisolver==1.4.7
+llvmlite==0.43.0
+matplotlib==3.9.2
+ml_dtypes==0.5.0
+numba==0.60.0
+numpy==2.0.2
+opt_einsum==3.4.0
+packaging==24.1
+pillow==10.4.0
+pyparsing==3.1.4
+python-dateutil==2.9.0.post0
+scipy==1.13.1
+six==1.16.0
+zipp==3.20.2

source/activity.py

@@ -0,0 +1,60 @@
+from __future__ import annotations
+import numpy as np
+from numpy import log as ln
+
+def get_activity(D, config, get_all=False):
+    """Computes activity coefficient for model species. For ideal mixing, all ln(g) terms equal
+       zero such that there is no deviation from ideal behaviour. Non-ideal mixing is captured by the
+       computed lng (ln gamma) terms below for Si, O, H2, H2O in melt phase, H in metal phase and xB.
+       Optional boolean parameter get_all controls whether all activities are computed and returned,
+       or if only the activities of Si, O and xB are returned.
+       Parameters:
+       D (dict)                 : Dictionary representing variable input for model; 
+                                retrieves mole fractions of species used in calculating coefficients.
+       config (dict)            : Dictionary of global model input, retrieves parameters that are not 
+                                varied in the model: T_surface, T_eq, bool_ideal_mixing. The latter
+                                controls non-ideality of activity coefficients. If False, lng terms will
+                                be computed via empirical formulae. If True, lng terms are zero.
+       get_all (Bool)           : Boolean controlling whether terms for H2, H2O_melt and H_metal are 
+                                computed and returned. These equations are implemented but their 
+                                validity is still under review. Default is False.
+
+       Returns:
+       lng (list)               : List of lng terms for Si, O, H2, H2O (melt), H (metal) and xB."""
+
+    if config["bool_ideal_mixing"] is False:
+        lngSi = -6.65*1873.0/config["T_eq"]-(12.41*1873.0/config["T_eq"])*ln(1.0-D["Si_metal"]) - \
+        ((-5.0*1873.0/config["T_eq"])*D["O_metal"]*(1.0+ln(1-D["O_metal"])/D["O_metal"]-1.0/(1.0-D["Si_metal"]))) + \
+        (-5.0*1873.0/config["T_eq"])*D["O_metal"]**2.0*D["Si_metal"]*(1.0/(1.0-D["Si_metal"])+1.0/(1.0-D["O_metal"]) + \
+                                                             D["Si_metal"]/(2.0*(1.0-D["Si_metal"])**2.0)-1.0)
+        lngO = (4.29-16500.0/config["T_eq"])-(-1.0*1873.0/config["T_eq"])*ln(1.0-D["O_metal"]) - \
+        ((-5.0*1873.0/config["T_eq"])*D["Si_metal"]*(1.0+ln(1-D["Si_metal"])/D["Si_metal"]-1.0/(1.0-D["O_metal"]))) + \
+        (-5.0*1873.0/config["T_eq"])*D["Si_metal"]**2.0*D["O_metal"]*(1.0/(1.0-D["O_metal"])+1.0/(1.0-D["Si_metal"])+ \
+                                                             D["O_metal"]/(2.0*(1.0-D["O_metal"])**2.0)-1.0)
+
+        xB = (D["H2O_melt"]/(1-D["H2O_melt"])*(1/3)) / (1.0 + (D["H2O_melt"]/(1-D["H2O_melt"])*(1/3)))
+
+        if get_all is True:
+            lngH2 = (74829.6/(8.3144*config["T_eq"]))*(1.0 - D["H2_melt"])**2.0
+
+            lngH2O_melt = (74826.0/(8.3144*config["T_eq"]))*(1.0 - D["H2O_melt"])**2.0
+
+            lngH_metal_T_indep = -1.9*(1.0 - D["H_metal"])**2.0
+
+            x_light = D["Si_metal"] + D["O_metal"]
+            lngH_metal = -3.8*x_light*(1.0 + ln(1.0 - x_light)/x_light - 1.0/(1.0-D["H_metal"])) + 3.8*(x_light**2.0)*D["H_metal"]* \
+                (1.0/(1.0 - D["H_metal"]) + 1.0/(1.0 - x_light) + D["H_metal"]/(2.0*(1.0-D["H_metal"])**2.0) - 1.0)
+        else:
+            lngH2 = 0.0
+            lngH2O_melt = 0.0
+            lngH_metal = 0.0
+
+    elif config['bool_ideal_mixing'] is True:
+        lngSi = 0.0
+        lngO = 0.0
+        xB = 0.0
+        lngH2 = 0.0
+        lngH2O_melt = 0.0
+        lngH_metal = 0.0
+
+    return [lngSi, lngO, lngH2, lngH2O_melt, lngH_metal, xB]