Cross-comparison/Fig_resolution_marlpde.py at main · MindTheGap-ERC/Cross-comparison · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#!/usr/bin/env python3

import matplotlib.pyplot as plt
import numpy as np
import glob
import h5py

folders = [
    'data/resolution/marlpde/Tstar_10_200nnx',
    'data/resolution/marlpde/Tstar_10_300nnx',
    'data/resolution/marlpde/Tstar_10_400nnx',
]

labels = ["Aragonite", "Calcite", "Porosity", r"$Ca^{2+}$", r"$CO_3^{2-}$"]
line_styles = ['-', '--', '-.', ':', (0, (3, 1, 1, 1))]
colors = ['#E69F00', 'red', '#000000', '#009E73', '#0072B2']

# Depth index for time series (-1 = bottom of the system)
depth = -1

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

Ts = 131.9 / 0.1**2

all_y_values = []

for idx, folder in enumerate(folders):
    h5_files = glob.glob(f'{folder}/*.hdf5')
    if not h5_files:
        print(f"Warning: No .hdf5 file found in {folder}")
        continue
    filename = h5_files[0]
    ax = axes[idx]

    sol = h5py.File(filename, 'r')
    df = np.array(sol['data'])
    dt = np.array(sol['times'])
    sol.close()

    t_plot = dt * Ts / 1000

    AR  = df[:, 0, depth]
    CA  = df[:, 1, depth]
    phi = df[:, 4, depth]
    Ca  = df[:, 2, depth]
    CO  = df[:, 3, depth]

    ax.plot(t_plot, AR,  label=labels[0], color=colors[0], linestyle=line_styles[0])
    ax.plot(t_plot, CA,  label=labels[1], color=colors[1], linestyle=line_styles[1])
    ax.plot(t_plot, phi, label=labels[2], color=colors[2], linestyle=line_styles[2])
    ax.plot(t_plot, Ca,  label=labels[3], color=colors[3], linestyle=line_styles[3])
    ax.plot(t_plot, CO,  label=labels[4], color=colors[4], linestyle=line_styles[4])

    all_y_values.extend(AR)
    all_y_values.extend(CA)
    all_y_values.extend(phi)
    all_y_values.extend(Ca)
    all_y_values.extend(CO)

    folder_name = folder.split('/')[-1]
    parts = folder_name.split('_')
    nnx = parts[2].replace('nnx', '')

    ax.set_title(f'nnx = {nnx}', fontsize=10)
    ax.set_xlabel('Time [ky]')
    ax.set_ylabel('Concentration/Porosity')

    panel_label = chr(97 + idx)
    ax.text(0.02, 0.98, f'({panel_label})', transform=ax.transAxes,
            fontsize=12, fontweight='bold', va='top', ha='left')

    if idx == 2:
        ax.legend(loc='lower right', fontsize=8)

y_min = min(all_y_values)
y_max = max(all_y_values)
y_margin = (y_max - y_min) * 0.05
for ax in axes:
    ax.set_ylim(y_min - y_margin, y_max + y_margin)

plt.tight_layout()
plt.savefig('figs/Fig_resolution_marlpde.png',
            format='png', dpi=300, bbox_inches='tight')

plt.close()