Merge pull request #663 from bobmyhill/add_phase_diagram_mask

add optional phase diagram mask

Author: bobmyhill

contrib/perplex/create_phase_diagram.py

@@ -1,13 +1,14 @@
 import os
 import shutil
-
+import numpy as np
 import burnman
 from perplex_utils import databases, make_build_file
 from perplex_utils import run_vertex, run_pssect
 from perplex_utils import create_mode_table
 from perplex_utils import pretty_plot_phase_diagram
 
 import matplotlib.pyplot as plt
+import pickle
 
 
 def run_perplex(
@@ -71,25 +72,11 @@ def run_perplex(
         "molar",
     ),
     "harzburgite": burnman.Composition(
-        {
-            "SiO2": 36.07,
-            "MgO": 56.51,
-            "FeO": 6.07,
-            "CaO": 0.81,
-            "Al2O3": 0.53,
-            "Na2O": 0.001,
-        },
+        {"SiO2": 36.07, "MgO": 56.51, "FeO": 6.07, "CaO": 0.81, "Al2O3": 0.53},
         "molar",
     ),
     "modified_harzburgite": burnman.Composition(
-        {
-            "SiO2": 36.04,
-            "MgO": 56.54,
-            "FeO": 5.97,
-            "CaO": 0.79,
-            "Al2O3": 0.65,
-            "Na2O": 0.001,
-        },
+        {"SiO2": 36.04, "MgO": 56.54, "FeO": 5.97, "CaO": 0.79, "Al2O3": 0.65},
         "molar",
     ),
 }
@@ -100,15 +87,15 @@ def run_perplex(
     project_name = "basalt"
     database = databases["stx21"]
     composition = compositions[project_name]
-    pressure_range_total = [1.0e5, 140.0e9]
-    temperature_range_total = [200.0, 3000.0]
+    pressure_range_total = [1.0, 140.0e9]
+    temperature_range_total = [200.0, 4000.0]
 
     # Split pressure and temperature so that PerpleX
     # no temperature splits seems to make diagram with
     # less prominent discontinuities
-    n_pressures_per_split = 121
-    n_temperatures_per_split = 601
-    n_splits_pressure = 14
+    n_pressures_per_split = 101
+    n_temperatures_per_split = 761
+    n_splits_pressure = 28
     n_splits_temperature = 1
 
     # If this script has already been run, and you just want to
@@ -258,6 +245,10 @@ def run_perplex(
 
             modes_filenames.append(split_project_name + "_modes.dat")
 
+    mask_polygon = np.array(
+        [[0.0, 2150.0], [12.0e9, 4000.0001], [0.0, 4000.001], [0.0, 2150.0]]
+    )
+
     fig = plt.figure(figsize=(12, 8))
     ax = fig.add_subplot(1, 1, 1)
 
@@ -273,15 +264,23 @@ def run_perplex(
         label_scaling=3.0,
         label_clearance=0.01,
         number_small_fields=True,
+        mask_polygon=mask_polygon,
     )
 
-    with open(f"{project_name}_field_ids.txt", "w") as f:
+    # Save the assemblage data for all the small fields
+    with open(f"{project_name}_phase_diagram_field_ids.txt", "w") as f:
         for i, small_field in enumerate(small_fields):
             line = f"{i+1}: {small_field}"
             print(line)
             f.write(f"{line}\n")
 
-    fig.savefig(f"{project_name}.pdf")
+    # Save figure as a pdf
+    fig.savefig(f"{project_name}_phase_diagram.pdf")
+
+    # Save figure as a pickle file to reload later
+    with open(f"{project_name}_phase_diagram.pkl", "wb") as f:
+        pickle.dump(fig, f)
+
     plt.show()
 
     print("Processing complete.")

contrib/perplex/perplex_utils.py

@@ -809,7 +809,11 @@ def _smooth_polygon_edges_between_intersections(
 
 
 def get_fields_assemblages_and_bounds(
-    filename, phase_name_replacements, smoothing_window, smoothing_order
+    filename,
+    phase_name_replacements,
+    smoothing_window,
+    smoothing_order,
+    mask_polygon=None,
 ):
     """
     Extract phase assemblage fields and plot bounds from a mode tab file.
@@ -823,6 +827,9 @@ def get_fields_assemblages_and_bounds(
     :type smoothing_window: int
     :param smoothing_order: Polynomial order for smoothing polygon edges.
     :type smoothing_order: int
+    :param mask_polygon: 2D array of points (P in Pa, T in K),
+    defining a polygon to mask out.
+    :type mask_polygon: 2D numpy array
     :return: Tuple containing list of polygon data dicts and bounds
     [[P_min, P_max], [T_min, T_max]].
     :rtype: (list[dict], list[list[float]])
@@ -832,8 +839,21 @@ def get_fields_assemblages_and_bounds(
     P, T = data[:2]
     pressures = P[:, 0]
     temperatures = T[0]
-
     phase_modes = data[2:]
+
+    if mask_polygon is not None:
+
+        mask_polygon_bar = mask_polygon.copy()
+        mask_polygon_bar[:, 0] = mask_polygon_bar[:, 0] / 1.0e5
+
+        poly = Polygon(mask_polygon_bar)
+        points = np.vstack((P.ravel(), T.ravel())).T
+        mask = np.array(
+            [poly.contains(Point(pt)) or poly.touches(Point(pt)) for pt in points]
+        )
+        mask = mask.reshape(P.shape)
+        phase_modes[:, mask] = 0.0
+
     num_phases = phase_modes.shape[0]
     phase_names = [phase_name_replacements.get(name, name) for name in column_names[2:]]
 
@@ -1067,6 +1087,7 @@ def pretty_plot_phase_diagram(
     label_scaling=3.0,
     label_clearance=0.01,
     number_small_fields=True,
+    mask_polygon=None,
 ):
     """
     Plot the Perple_X calculated phase diagram on a matplotlib axis.
@@ -1097,6 +1118,9 @@ def pretty_plot_phase_diagram(
     :param number_small_fields: Replace assemblage labels with numbers if
     there is not enough clearance around the label.
     :type number_small_fields: bool
+    :param mask_polygon: 2D array of points (P in Pa, T in K),
+    defining a polygon to mask out.
+    :type mask_polygon: 2D numpy array
     :return: List of assemblages corresponding to numbered small fields.
     :rtype: [str]
     """
@@ -1110,6 +1134,7 @@ def pretty_plot_phase_diagram(
             phase_name_replacements,
             smoothing_window,
             smoothing_order,
+            mask_polygon,
         )
         polygon_data.extend(polygon_data_i)
         bounds.append(bounds_i)