|
6 | 6 | from perplex_utils import run_vertex, run_pssect |
7 | 7 | from perplex_utils import create_perplex_table |
8 | 8 |
|
| 9 | + |
| 10 | +def create_table( |
| 11 | + perplex_bindir, |
| 12 | + project_name, |
| 13 | + database, |
| 14 | + perplex_option_file, |
| 15 | + composition, |
| 16 | + pressure_range, |
| 17 | + temperature_range, |
| 18 | + n_pressures, |
| 19 | + n_temperatures, |
| 20 | + outfile, |
| 21 | +): |
| 22 | + # Make the build file |
| 23 | + # This file contains the information about the project, such as the |
| 24 | + # composition, the database, and the options for Perple_X. |
| 25 | + print("* Making build file...") |
| 26 | + make_build_file( |
| 27 | + perplex_bindir, |
| 28 | + project_name, |
| 29 | + database, |
| 30 | + perplex_option_file, |
| 31 | + composition, |
| 32 | + pressure_range, |
| 33 | + temperature_range, |
| 34 | + verbose=False, |
| 35 | + ) |
| 36 | + |
| 37 | + # Run the Perple_X vertex command to do the Gibbs minimization |
| 38 | + # and create output files containing the equilibrium assemblages. |
| 39 | + print("* Running vertex...") |
| 40 | + run_vertex(perplex_bindir, project_name, verbose=False) |
| 41 | + |
| 42 | + # Run the Perple_X pssect command to create a postscript |
| 43 | + # file containing the phase diagram. |
| 44 | + print("* Running pssect...") |
| 45 | + run_pssect(perplex_bindir, project_name, convert_to_pdf=False, verbose=False) |
| 46 | + |
| 47 | + # Create the BurnMan-readable table from the vertex output |
| 48 | + print("* Creating BurnMan-readable table...") |
| 49 | + create_perplex_table( |
| 50 | + perplex_bindir, |
| 51 | + project_name, |
| 52 | + outfile, |
| 53 | + n_pressures, |
| 54 | + n_temperatures, |
| 55 | + pressure_range, |
| 56 | + temperature_range, |
| 57 | + ) |
| 58 | + |
| 59 | + |
9 | 60 | if __name__ == "__main__": |
10 | 61 | # Define the project parameters |
11 | 62 | project_name = "iron_olivine_lo_res" |
12 | 63 | database = databases["stx24"] |
13 | 64 | composition = burnman.Composition( |
14 | 65 | {"MgO": 1.8, "FeO": 0.2, "SiO2": 1.0, "Fe": 0.1}, "molar" |
15 | 66 | ) |
16 | | - pressure_range = [1.0e5, 10.0e9] |
17 | | - temperature_range = [200.0, 3000.0] |
18 | | - n_pressures = 11 |
19 | | - n_temperatures = 15 |
20 | | - outfile = "iron_olivine_lo_res_table.dat" |
| 67 | + pressure_range_total = [1.0e5, 10.0e9] |
| 68 | + temperature_range_total = [200.0, 3000.0] |
| 69 | + n_pressures_per_split = 6 |
| 70 | + n_temperatures_per_split = 8 |
| 71 | + n_splits_pressure = 2 |
| 72 | + n_splits_temperature = 2 |
| 73 | + |
| 74 | + outfile_base = "iron_olivine_lo_res_table" |
21 | 75 |
|
22 | 76 | perplex_dir = os.path.join(os.getcwd(), "perplex-installer/Perple_X") |
23 | 77 | perplex_bindir = os.path.join(os.getcwd(), "perplex-installer/Perple_X/bin") |
|
50 | 104 | # In the refinement stage, the number of nodes is increased. |
51 | 105 | # Any solutions not present on any of the bounding nodes will be excluded |
52 | 106 | # from the calculations. |
53 | | - n_P_exploratory = int((n_pressures + 1) / 2) |
54 | | - n_T_exploratory = int((n_temperatures + 1) / 2) |
| 107 | + n_P_exploratory = int((n_pressures_per_split + 1) / 2) |
| 108 | + n_T_exploratory = int((n_temperatures_per_split + 1) / 2) |
55 | 109 |
|
56 | 110 | # Create the Perple_X options file |
57 | 111 | with open(perplex_option_file, "w") as f: |
|
61 | 115 | f.write("auto_refine auto\n") # Automatically refine the grid |
62 | 116 | f.write("grid_levels 1 1\n") # Do not use adaptive grid refinement |
63 | 117 | f.write( |
64 | | - f"x_nodes {n_P_exploratory} {n_pressures}\n" |
| 118 | + f"x_nodes {n_P_exploratory} {n_pressures_per_split}\n" |
65 | 119 | ) # Exploratory and final number of pressure nodes |
66 | 120 | f.write( |
67 | | - f"y_nodes {n_T_exploratory} {n_temperatures}\n" |
| 121 | + f"y_nodes {n_T_exploratory} {n_temperatures_per_split}\n" |
68 | 122 | ) # Exploratory and final number of temperature nodes |
69 | 123 |
|
70 | | - # Make the build file |
71 | | - # This file contains the information about the project, such as the |
72 | | - # composition, the database, and the options for Perple_X. |
73 | | - print("Making build file...") |
74 | | - make_build_file( |
75 | | - perplex_bindir, |
76 | | - project_name, |
77 | | - database, |
78 | | - perplex_option_file, |
79 | | - composition, |
80 | | - pressure_range, |
81 | | - temperature_range, |
82 | | - verbose=False, |
83 | | - ) |
| 124 | + # Create the table(s) using the defined parameters |
| 125 | + delta_P = (pressure_range_total[1] - pressure_range_total[0]) / n_splits_pressure |
| 126 | + delta_T = ( |
| 127 | + temperature_range_total[1] - temperature_range_total[0] |
| 128 | + ) / n_splits_temperature |
84 | 129 |
|
85 | | - # Run the Perple_X vertex command to do the Gibbs minimization |
86 | | - # and create output files containing the equilibrium assemblages. |
87 | | - print("Running vertex...") |
88 | | - run_vertex(perplex_bindir, project_name, verbose=False) |
| 130 | + for i_P in range(n_splits_pressure): |
| 131 | + for i_T in range(n_splits_temperature): |
| 132 | + pressure_range = [ |
| 133 | + pressure_range_total[0] + i_P * delta_P, |
| 134 | + pressure_range_total[0] + (i_P + 1) * delta_P, |
| 135 | + ] |
| 136 | + temperature_range = [ |
| 137 | + temperature_range_total[0] + i_T * delta_T, |
| 138 | + temperature_range_total[0] + (i_T + 1) * delta_T, |
| 139 | + ] |
| 140 | + if n_splits_pressure > 1 or n_splits_temperature > 1: |
| 141 | + outfile = f"{outfile_base}_P{i_P + 1:02d}_T{i_T + 1:02d}.dat" |
| 142 | + else: |
| 143 | + outfile = f"{outfile_base}.dat" |
89 | 144 |
|
90 | | - # Run the Perple_X pssect command to create a postscript |
91 | | - # file containing the phase diagram. |
92 | | - print("Running pssect...") |
93 | | - run_pssect(perplex_bindir, project_name, convert_to_pdf=False, verbose=False) |
| 145 | + print( |
| 146 | + f"Creating table [{i_P + 1}/{n_splits_pressure}, {i_T + 1}/{n_splits_temperature}]" |
| 147 | + ) |
| 148 | + # Create the table for the current split |
| 149 | + create_table( |
| 150 | + perplex_bindir, |
| 151 | + project_name, |
| 152 | + database, |
| 153 | + perplex_option_file, |
| 154 | + composition, |
| 155 | + pressure_range, |
| 156 | + temperature_range, |
| 157 | + n_pressures_per_split, |
| 158 | + n_temperatures_per_split, |
| 159 | + outfile, |
| 160 | + ) |
| 161 | + |
| 162 | + # If there are splits, merge the output files into one |
| 163 | + if n_splits_pressure > 1 or n_splits_temperature > 1: |
| 164 | + print("Merging output files...") |
| 165 | + |
| 166 | + # Read the header lines from the first split file |
| 167 | + first_split_outfile = f"{outfile_base}_P01_T01.dat" |
| 168 | + with open(first_split_outfile, "r") as infile: |
| 169 | + header_lines = [next(infile) for _ in range(13)] |
| 170 | + |
| 171 | + total_n_pressures = (n_pressures_per_split - 1) * n_splits_pressure + 1 |
| 172 | + total_n_temperatures = (n_temperatures_per_split - 1) * n_splits_temperature + 1 |
| 173 | + header_lines[6] = f" {total_n_pressures}\n" |
| 174 | + header_lines[10] = f" {total_n_temperatures}\n" |
| 175 | + |
| 176 | + # Load all other lines from the split files into a single list of lists |
| 177 | + all_lines = [] |
| 178 | + for i_P in range(n_splits_pressure): |
| 179 | + for i_T in range(n_splits_temperature): |
| 180 | + split_outfile = f"{outfile_base}_P{i_P + 1:02d}_T{i_T + 1:02d}.dat" |
| 181 | + with open(split_outfile, "r") as infile: |
| 182 | + lines = infile.readlines()[13:] |
| 183 | + all_lines.extend(lines) |
| 184 | + |
| 185 | + # Sort the lines by pressure and temperature, looping over pressures first |
| 186 | + all_lines.sort(key=lambda x: (float(x.split()[1]), float(x.split()[0]))) |
| 187 | + |
| 188 | + # Remove duplicate pressures and temperatures from the sorted list |
| 189 | + unique_lines = [] |
| 190 | + seen = set() |
| 191 | + for line in all_lines: |
| 192 | + pressure_temp = (float(line.split()[0]), float(line.split()[1])) |
| 193 | + if pressure_temp not in seen: |
| 194 | + seen.add(pressure_temp) |
| 195 | + unique_lines.append(line) |
| 196 | + |
| 197 | + # Check the total number of unique lines is as expected |
| 198 | + if len(unique_lines) != total_n_pressures * total_n_temperatures: |
| 199 | + raise ValueError( |
| 200 | + f"Expected {total_n_pressures * total_n_temperatures} unique lines, " |
| 201 | + f"but found {len(unique_lines)}." |
| 202 | + ) |
| 203 | + # Write the merged output file |
| 204 | + merged_outfile = f"{outfile_base}.dat" |
| 205 | + with open(merged_outfile, "w") as outfile: |
| 206 | + outfile.writelines(header_lines) |
| 207 | + outfile.writelines(unique_lines) |
| 208 | + print(f"Merged output file created: {merged_outfile}") |
94 | 209 |
|
95 | | - # Create the BurnMan-readable table from the vertex output |
96 | | - print("Creating BurnMan-readable table...") |
97 | | - create_perplex_table( |
98 | | - perplex_bindir, |
99 | | - project_name, |
100 | | - outfile, |
101 | | - n_pressures, |
102 | | - n_temperatures, |
103 | | - pressure_range, |
104 | | - temperature_range, |
105 | | - ) |
106 | 210 | print("Processing complete.") |
0 commit comments