|
| 1 | +""" |
| 2 | +This script generates a bifurcation diagram and spatial solutions for a |
| 3 | +planar bifurcation problem. It reads data from specified files, processes |
| 4 | +the data, and creates contour plots to visualize the results. The script |
| 5 | +uses Matplotlib for plotting and NumPy for numerical operations. |
| 6 | +
|
| 7 | +To run this script, ensure you have the required libraries installed: |
| 8 | +- numpy |
| 9 | +- matplotlib |
| 10 | +You can install them using pip: |
| 11 | +pip install numpy matplotlib |
| 12 | +
|
| 13 | +To execute the script, run the following command in your terminal: |
| 14 | +python plot_solutions.py |
| 15 | +""" |
| 16 | + |
| 17 | +import numpy as np |
| 18 | +import matplotlib.pyplot as plt |
| 19 | +from matplotlib.gridspec import GridSpec |
| 20 | +from matplotlib.colors import ListedColormap, BoundaryNorm |
| 21 | + |
| 22 | +# Uncomment the following lines to use LaTeX for text rendering |
| 23 | +plt.rcParams.update({ |
| 24 | + "text.usetex": True, |
| 25 | + "font.family": "sans-serif", |
| 26 | + "font.sans-serif": "Helvetica", |
| 27 | + 'text.latex.preamble': r'\usepackage{amsfonts}' |
| 28 | +}) |
| 29 | + |
| 30 | + |
| 31 | +def cmap(Z, RES=12): |
| 32 | + """ |
| 33 | + Generate a custom colormap and normalization for visualizing data. |
| 34 | +
|
| 35 | + This function creates a colormap based on the 'RdBu_r' colormap, with a |
| 36 | + specified number of levels. The color corresponding to the value 0 is |
| 37 | + replaced with white to emphasize the zero level. The function also |
| 38 | + returns the levels and a normalization object for mapping data values |
| 39 | + to the colormap. |
| 40 | +
|
| 41 | + Parameters: |
| 42 | + ----------- |
| 43 | + Z : numpy.ndarray |
| 44 | + A 2D array of data values for which the colormap will be created. |
| 45 | + The maximum absolute value in `Z` determines the range of the colormap. |
| 46 | + RES : int, optional |
| 47 | + The number of levels in the colormap. Default is 12. It must be even. |
| 48 | + This parameter is used to create evenly spaced levels for the colormap. |
| 49 | + Returns: |
| 50 | + -------- |
| 51 | + levels : numpy.ndarray |
| 52 | + An array of evenly spaced levels for the colormap, symmetric around 0. |
| 53 | + custom_cmap : matplotlib.colors.ListedColormap |
| 54 | + A custom colormap with the color corresponding to 0 replaced by white. |
| 55 | + norm : matplotlib.colors.BoundaryNorm |
| 56 | + A normalization object for mapping data values to the colormap. |
| 57 | +
|
| 58 | + Notes: |
| 59 | + ------ |
| 60 | + - The function assumes that the input array `Z` is 2D and contains numeric |
| 61 | + values. |
| 62 | +
|
| 63 | + Example: |
| 64 | + -------- |
| 65 | + >>> import numpy as np |
| 66 | + >>> Z = np.random.randn(10, 10) |
| 67 | + >>> levels, custom_cmap, norm = cmap(Z) |
| 68 | + """ |
| 69 | + |
| 70 | + print(np.asarray([min(Z.flatten()), max(Z.flatten())])) |
| 71 | + Max = np.max(abs(Z)) |
| 72 | + |
| 73 | + assert RES % 2 == 0 # ensure RES is even |
| 74 | + |
| 75 | + # Define levels (must include 0) |
| 76 | + levels = np.linspace(-Max, Max, RES) |
| 77 | + |
| 78 | + # Create RdBu_r colormap with 12 colors |
| 79 | + original_cmap = plt.get_cmap('RdBu_r', len(levels) - 1) |
| 80 | + colors = original_cmap(np.arange(original_cmap.N)) |
| 81 | + |
| 82 | + # Find the index of the bin that contains 0 |
| 83 | + # Since levels are symmetric, 0 will be in the middle bin |
| 84 | + zero_bin_index = (len(levels) - 1) // 2 # e.g., 6 if 12 bins |
| 85 | + |
| 86 | + # Replace that color with white |
| 87 | + colors[zero_bin_index] = [1, 1, 1, 1] # RGBA for white |
| 88 | + custom_cmap = ListedColormap(colors) |
| 89 | + |
| 90 | + # Create normalization |
| 91 | + norm = BoundaryNorm(levels, custom_cmap.N) |
| 92 | + |
| 93 | + return levels, custom_cmap, norm |
| 94 | + |
| 95 | + |
| 96 | +def transform(file): |
| 97 | + """ |
| 98 | + Transforms data from a file into structured arrays for plotting. |
| 99 | +
|
| 100 | + This function reads data from a file, extracts specific columns, and reshapes |
| 101 | + the data into 2D arrays suitable for plotting. It also removes an offset from |
| 102 | + the computed values. |
| 103 | +
|
| 104 | + Args: |
| 105 | + file (str): Path to the file containing the data. The file should be in a |
| 106 | + format readable by `numpy.loadtxt` and contain at least 5 columns. |
| 107 | +
|
| 108 | + Returns: |
| 109 | + tuple: A tuple containing three 2D numpy arrays: |
| 110 | + - X (numpy.ndarray): 2D array of x-coordinates (reshaped from column 1). |
| 111 | + - Z (numpy.ndarray): 2D array of z-coordinates (reshaped from column 2). |
| 112 | + - T (numpy.ndarray): 2D array of transformed values (computed from column 5). |
| 113 | + """ |
| 114 | + # Load data from the file |
| 115 | + data = np.loadtxt(file) |
| 116 | + |
| 117 | + # Extract columns for plotting (column 1, column 2, and computed value for color) |
| 118 | + x = data[:, 0] # Column 1 - but python is zero-indexed |
| 119 | + z = data[:, 1] # Column 2 |
| 120 | + t = data[:, 4] # Column 5 |
| 121 | + |
| 122 | + Nz = 17 |
| 123 | + Nx = 500 |
| 124 | + X = np.zeros((Nz, Nx)) |
| 125 | + Z = np.zeros((Nz, Nx)) |
| 126 | + T = np.zeros((Nz, Nx)) |
| 127 | + for i in range(Nx): |
| 128 | + X[:, i] = x[i*Nz:(i+1)*Nz] |
| 129 | + Z[:, i] = z[i*Nz:(i+1)*Nz] |
| 130 | + T[:, i] = t[i*Nz:(i+1)*Nz] |
| 131 | + |
| 132 | + T += (Z - 1/2) # Remove the offset |
| 133 | + |
| 134 | + # Print file and max/min |
| 135 | + print(file, np.asarray([min(T.flatten()), max(T.flatten())]), '\n') |
| 136 | + |
| 137 | + return X, Z, T |
| 138 | + |
| 139 | + |
| 140 | +print('Load above') |
| 141 | +dir = '/home/pmannix/Spatial_Localisation/SpectralDoubleDiffusiveConvection/Paper_Data/Beamue2011_Data/' |
| 142 | + |
| 143 | +# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~ |
| 144 | + |
| 145 | +fig = plt.figure(figsize=(16, 6), layout='constrained') |
| 146 | + |
| 147 | +# Create a GridSpec with custom width ratios |
| 148 | +gs = GridSpec(6, 2, figure=fig, width_ratios=[0.5, 1.5]) # Adjust width_ratios as needed |
| 149 | + |
| 150 | +# Create the left column as a single subplot |
| 151 | +ax = fig.add_subplot(gs[:, 0]) |
| 152 | + |
| 153 | +# Create the right column subplots |
| 154 | +axs = np.empty((6, 1), dtype=object) |
| 155 | +for i in range(6): |
| 156 | + axs[i, 0] = fig.add_subplot(gs[i, 1], aspect=.75) |
| 157 | + |
| 158 | + |
| 159 | +# A) Bifurcation diagram |
| 160 | +# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~ |
| 161 | + |
| 162 | +# Load data from bifdiag1.dat and bifdiag2.dat |
| 163 | +data1 = np.loadtxt(dir + 'Planarfigure1_2/bifdiag1.dat') |
| 164 | +data2 = np.loadtxt(dir + 'Planarfigure1_2/bifdiag2.dat') |
| 165 | + |
| 166 | +# Extract columns for plotting (column 2 and column 7, zero-indexed) |
| 167 | +x1, y1 = data1[:, 1], data1[:, 6] |
| 168 | +x2, y2 = data2[:, 1], data2[:, 6] |
| 169 | + |
| 170 | +# Create the plot |
| 171 | +ax.plot(x1, y1, 'k-', linewidth=2, label=r'$L^{1-}_{10}$') |
| 172 | +ax.plot(x2, y2, 'k:', linewidth=2, label=r'$L^{1+}_{10}$') |
| 173 | + |
| 174 | +# Add labels, legend, and title |
| 175 | +ax.set_xlabel(r'$Ra_T$', fontsize=30) |
| 176 | +ax.set_ylabel(r'$\mathcal{E}$', fontsize=30) |
| 177 | + |
| 178 | +# Show the plot |
| 179 | +ax.set_xlim(2660, 2800) |
| 180 | +ax.set_ylim(0, 290) |
| 181 | +ax.tick_params(axis='both', which='major', labelsize=20) |
| 182 | + |
| 183 | +E = np.array([35-2, 65-4, 90-3, 25-3, 50-3, 75]) |
| 184 | +Ra= np.array([2700-4, 2700-7, 2700-7, 2710-6, 2700-6, 2700-6]) |
| 185 | +letters = [r'1',r'2',r'3', r'a',r'b',r'c'] |
| 186 | +for i, xy in enumerate(zip(Ra, E)): |
| 187 | + ax.annotate(letters[i], xy=xy, textcoords='data', fontsize=15) |
| 188 | + |
| 189 | +# B) Spatial solutions |
| 190 | +# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~ |
| 191 | + |
| 192 | +files = ['Planarfigure2_2/SOL_sn1.plt','Planarfigure2_2/SOL_sn2.plt', |
| 193 | + 'Planarfigure2_2/SOL_sn3.plt','Planarfigure2_2/SOL_sna.plt', |
| 194 | + 'Planarfigure2_2/SOL_snb.plt','Planarfigure2_2/SOL_snc.plt'] |
| 195 | + |
| 196 | +# Create a contour plot |
| 197 | +for i, file in enumerate(files): |
| 198 | + X, Z, T = transform(dir + file) |
| 199 | + levels, custom_cmap, norm = cmap(T) |
| 200 | + axs[::-1, 0][i].contour(X, Z, T, levels=levels, cmap=custom_cmap, norm=norm) |
| 201 | + axs[::-1, 0][i].contourf(X, Z, T, levels=levels, cmap=custom_cmap, norm=norm) |
| 202 | + axs[::-1, 0][i].set_xticks([]) |
| 203 | + axs[::-1, 0][i].set_yticks([]) |
| 204 | + |
| 205 | + |
| 206 | +letters = [r'1',r'2',r'3', r'a',r'b',r'c'] |
| 207 | +for i, letter in enumerate(letters): |
| 208 | + axs[::-1, 0][i].annotate(letter, xy=(-0.025, 0.5), xycoords='axes fraction', fontsize=20, color='k') |
| 209 | + |
| 210 | +# ~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~ |
| 211 | + |
| 212 | +# Save the figure |
| 213 | +plt.savefig('planar_bifurcation_diagram.png', dpi=100) |
| 214 | +plt.show() |
0 commit comments