Kernel: Python 3 (CoCalc)
In [17]:
import numpy as np import matplotlib.pyplot as plt import glob import re import os from collections import OrderedDict import matplotlib.gridspec as gridspec import pickle
In [18]:
from Dabade16 import deltaC_dabade16, calculate_eccentricity from Einarsson15 import deltaC_einarsson15
In [19]:
###load the experiments fittings of C,dC with open('C_dC_dictionary.pkl', 'rb') as f: exp = pickle.load(f) exp = dict(sorted(exp.items(), key=lambda item: item[0])) exp[0.43] = exp[0.47438694107364204] del exp[0.47438694107364204] exp = dict(sorted(exp.items(), key=lambda item: item[0]))
In [20]:
#aspect ratio r=0.56 #Reynolds number Rep=0.47 #Initial condition IC = np.array([0.0,0.9,0.05]) IC = IC / np.linalg.norm(IC) #time-line time = np.linspace(0,500.0,200000) #shear-rate gammadot=3.1
In [21]:
###calculate the theory of dabade dabade16 = deltaC_dabade16(r=r) e0 = calculate_eccentricity(r) ###scaling factor for the oblate ellipsoid ftr = e0**2
In [22]:
einarsson15 = deltaC_einarsson15(r=r,Rep=Rep,dir='',IC=IC,time=time,gammadot=gammadot)
Out[22]:
/Users/ziqiwang/PD/MP_PD/2025-JFM/Jupyter_notebook_figures/figure_4/Einarsson15.py:89: RuntimeWarning: invalid value encountered in arccos
theta = np.arccos(n3[orbit_inds])
Analyse the experiments
In [23]:
ftr = np.loadtxt('C_dC_map/epsilon_zero_ELL06.txt') _,einars15 = np.loadtxt('C_t_kappa_theory.txt',skiprows=1).T
In [24]:
markers = ['o', 's', 'D', '^', 'v', '<', '>', 'P'] green_colors = ['#F7FCF5', '#DDEEDB','#BFE7B7','#8FD08A','#66BD63','#31A354','#006D2C','#00441B'] green_colors = green_colors[::-1] blue_colors = ['#E6F1FA', '#DBE9F6', '#BAD6EB', '#89BEDC', '#539ECD', '#2B7BBA', '#0B559F', '#08306B'] blue_colors = blue_colors[::-1]
In [25]:
all_C = [] all_dC = [] all_reps = [] all_kappas = []
Panel (a)
In [26]:
# Path to your data files file_pattern = "Fig_4_a_data/C_t_Rep_*_kappa_*.dat" # Dictionary to store data data_dict = {} # Regex to extract Rep and kappa pattern = re.compile(r"C_t_Rep_([0-9.]+)_kappa_([0-9.]+)\.dat") # Load all files for filename in glob.glob(file_pattern): match = pattern.search(filename) if match: Rep = float(match.group(1)) kappa = float(match.group(2)) # Load data (adjust depending on your file format) data = np.loadtxt(filename) # Store in dictionary data_dict[(Rep, kappa)] = data Rep_kappa_dict = {} for (rep, kappa), data in data_dict.items(): if rep not in Rep_kappa_dict: Rep_kappa_dict[rep] = {} Rep_kappa_dict[rep][kappa] = data # Sort both Rep and kappa keys Rep_kappa_dict = OrderedDict( (rep, OrderedDict(sorted(kappa_dict.items()))) for rep, kappa_dict in sorted(Rep_kappa_dict.items()) )
In [27]:
def calculate_C_dC(C_in): """calculate C and dC for plotting from the 1D C_in array of orbit constants C""" C = C_in[1:-1] dC = C_in[2:] - C_in[:-2] return C,dC
In [28]:
fig,axs = plt.subplots(1,3,figsize=(10,3),sharey=True) for j,rep in enumerate(Rep_kappa_dict.keys()): for i, kappa in enumerate(Rep_kappa_dict[rep].keys()): _,Cs = Rep_kappa_dict[rep][kappa].T C, dC = calculate_C_dC(Cs) axs[j].plot(C/(C+1),dC/(C**2+1)*ftr/rep,c=green_colors[i+1],marker=markers[i],label=rf'$\kappa={kappa}$', markevery=5,markersize=3) all_C.append(C) all_dC.append(dC) all_reps.append(rep) all_kappas.append(kappa) C, dC = calculate_C_dC(einars15) #axs[j].plot(C/(C+1),dC/(C**2+1)*ftr/0.47,ls='--',c='red',label=r'$Einarsson15$') C=dabade16[:,0] dC = dabade16[:,1]+dabade16[:,2] #axs[j].plot(C/(C+1),dC/(C**2+1)*ftr,ls=':',c='k',label=r'$Dabade16$') axs[j].legend(ncol=2,loc='upper center',bbox_to_anchor=(0.5,1.55)) axs[j].set_title(rf'$Re_p = {rep}$') axs[j].set_xlabel(r'$C/(C+1)$',fontsize=10,labelpad=0.5); axs[0].set_ylabel(r'$\Delta C/(C^2+1)\times Re_p^{-1} \times {\xi_0}^2$',fontsize=10); #plt.savefig('ELL06_phase_space_varying_kappa_3_Reps.png',format='png',dpi=300, bbox_inches='tight');
Out[28]:
Panel b
In [29]:
# Path to your data files file_pattern = "Fig_4_b_data/C_t_Rep_*_kappa_*.dat" # Dictionary to store data data_dict = {} # Regex to extract Rep and kappa pattern = re.compile(r"C_t_Rep_([0-9.]+)_kappa_([0-9.]+)\.dat") # Load all files for filename in glob.glob(file_pattern): match = pattern.search(filename) if match: Rep = float(match.group(1)) kappa = float(match.group(2)) # Load data (adjust depending on your file format) data = np.loadtxt(filename) # Store in dictionary data_dict[(Rep, kappa)] = data kappa_Rep_dict = {} for (rep, kappa), data in data_dict.items(): if kappa not in kappa_Rep_dict: kappa_Rep_dict[kappa] = {} kappa_Rep_dict[kappa][rep] = data # Sort both Rep and kappa keys kappa_Rep_dict = OrderedDict( (kappa, OrderedDict(sorted(rep_dict.items()))) for kappa, rep_dict in sorted(kappa_Rep_dict.items()) )
In [30]:
fig,axs = plt.subplots(1,3,figsize=(10,3),sharey=True) for j,kappa in enumerate(kappa_Rep_dict.keys()): for i, rep in enumerate(kappa_Rep_dict[kappa].keys()): _,Cs = kappa_Rep_dict[kappa][rep].T C, dC = calculate_C_dC(Cs) axs[j].plot(C/(C+1),dC/(C**2+1)*ftr/rep,c=blue_colors[i+1],marker=markers[i],label=rf'$Re_p={rep}$', markevery=5,markersize=3) all_C.append(C) all_dC.append(dC) all_reps.append(rep) all_kappas.append(kappa) C, dC = calculate_C_dC(einars15) #axs[j].plot(C/(C+1),dC/(C**2+1)*ftr/0.47,ls='--',c='red',label=r'$Einarsson15$') C=dabade16[:,0] dC = dabade16[:,1]+dabade16[:,2] #axs[j].plot(C/(C+1),dC/(C**2+1)*ftr,ls=':',c='k',label=r'$Dabade16$') ###legend axs[j].legend(ncol=2,loc='upper center',bbox_to_anchor=(0.5,1.6)) axs[j].set_title(rf'$\kappa = {kappa}$') axs[j].set_xlabel(r'$C/(C+1)$',fontsize=10,labelpad=0.5); axs[0].set_ylabel(r'$\Delta C/(C^2+1)\times Re_p^{-1} \times {\xi_0}^2$',fontsize=10); #plt.savefig('ELL06_phase_space_varying_Rep_3_kappas.png',format='png',dpi=300, bbox_inches='tight');
Out[30]:
Plot the joint figure
In [31]:
label_fontsize=9 marker_size=3 title_fontsize=10 letters_fontsize=8 xticks_fontsize=8 yticks_fontsize=8 exp_colors = ["#FFD700", "#FF8C00", "#DAA520"] exp_makers = ['p','h','8']
In [32]:
scaling_factor = 1.0 fig = plt.figure(figsize=(7.5*scaling_factor, 7*scaling_factor)) n_panels = 2 plt.rcParams.update({'xtick.labelsize': xticks_fontsize, 'ytick.labelsize': yticks_fontsize}) #grid of 2 panels stacked vertically outer_gs = gridspec.GridSpec(n_panels, 1, height_ratios=[1]*n_panels,hspace=0.9) #first panel inner_gs = gridspec.GridSpecFromSubplotSpec(1, 3,subplot_spec=outer_gs[0],wspace=0.1, hspace=0.1) # Create subplots for this panel axes = [] for j,rep in enumerate(Rep_kappa_dict.keys()): ax = fig.add_subplot(inner_gs[0, j]) for i, kappa in enumerate(Rep_kappa_dict[rep].keys()): _,Cs = Rep_kappa_dict[rep][kappa].T C, dC = calculate_C_dC(Cs) ax.plot(C/(C+1),dC/(C**2+1)*ftr/rep,c=green_colors[i+1],marker=markers[i],label=rf'$\kappa={kappa}$', markevery=5,markersize=3) ###theory of Einarsson C, dC = calculate_C_dC(einarsson15[0,:]) # ax.plot(C/(C+1),dC/(C**2+1)*ftr/0.47,ls='--',c='red',label=r'$Einarsson \ et \ al., 2015$') # ax.plot(C/(C+1),dC/(C**2+1)*ftr/0.47,ls='--',c='red',label=r'Theory') ###theory of Dabade C=dabade16[:,0] dC = dabade16[:,1]+dabade16[:,2] # ax.plot(C/(C+1),dC/(C**2+1)*ftr,ls=':',c='k',label=r'$Dabade \ et \ al., 2016$') ax.plot(C/(C+1),dC/(C**2+1)*ftr,ls='--',c='red',label=r'Theory') if j==1: ax.legend(ncol=3,loc='upper center',bbox_to_anchor=(0.5,1.6)) if j==0: ax.set_ylabel(r'$\Delta C/(C^2+1)\times \text{Re}_p^{-1} \times {\xi_0}^2$',fontsize=label_fontsize); if j>0: ax.set_yticklabels([]) ax.set_title(rf'$\text{{Re}}_p = {rep}$',fontsize=title_fontsize) ax.set_xlabel(r'$C/(C+1)$',fontsize=label_fontsize,labelpad=0.5); ax.set_xlim(0,1) ax.set_ylim(-0.5,0.05) ax.axhline(0.0,c='grey',alpha=0.3,zorder=0) if j==0: ax.text(-0.15, 1.15, r'$(a)$',transform=ax.transAxes,fontsize=letters_fontsize, fontweight='bold', va='top', ha='left') #second panel inner_gs = gridspec.GridSpecFromSubplotSpec(1, 3,subplot_spec=outer_gs[1],wspace=0.1, hspace=0.1) # Create subplots for this panel axes = [] for j,kappa in enumerate(kappa_Rep_dict.keys()): ax = fig.add_subplot(inner_gs[0, j]) for i, rep in enumerate(kappa_Rep_dict[kappa].keys()): _,Cs = kappa_Rep_dict[kappa][rep].T C, dC = calculate_C_dC(Cs) ax.plot(C/(C+1),dC/(C**2+1)*ftr/rep,c=blue_colors[i+1],marker=markers[i],label=rf'$\text{{Re}}_p={rep}$', markevery=5,markersize=3) ###experiments # if j == 0: # exp_handles = [] # for i,Rep in enumerate(exp.keys()): # if Rep > 0.1 and Rep < 1.5: # for k in range(len(exp[Rep])): # C,dC = exp[Rep][k] # exp_hand = ax.scatter(C/(C+1),dC/(C**2+1)*ftr/Rep,c=exp_colors[i],marker=exp_makers[i],label=rf'$Exp.-Re_p={Rep:.2f}$') # exp_handles.append(exp_hand) ###theory of Einarsson C, dC = calculate_C_dC(einarsson15[0,:]) # ax.plot(C/(C+1),dC/(C**2+1)*ftr/0.47,ls='--',c='red',label=r'$Einarsson \ et \ al., 2015$') # ax.plot(C/(C+1),dC/(C**2+1)*ftr/0.47,ls='--',c='red',label=r'Theory') ###theory of Dabade C=dabade16[:,0] dC = dabade16[:,1]+dabade16[:,2] # ax.plot(C/(C+1),dC/(C**2+1)*ftr,ls=':',c='k',label=r'$Dabade \ et \ al., 2016$') ax.plot(C/(C+1),dC/(C**2+1)*ftr,ls='--',c='red',label=r'Theory') ###legend if j==1: # handles = exp_handles + ax.get_legend_handles_labels()[0] handles = ax.get_legend_handles_labels()[0] ax.legend(ncol=2,loc='upper center',bbox_to_anchor=(0.4,1.65),handles=handles) if j==0: ax.set_ylabel(r'$\Delta C/(C^2+1)\times \text{Re}_p^{-1} \times {\xi_0}^2$',fontsize=label_fontsize); if j>0: ax.set_yticklabels([]) ax.set_title(rf'$\kappa = {kappa}$',fontsize=title_fontsize) ax.set_xlabel(r'$C/(C+1)$',fontsize=label_fontsize,labelpad=0.5); ax.set_xlim(0,1) ax.set_ylim(-0.5,0.05) ax.axhline(0.0,c='grey',alpha=0.3,zorder=0) if j==0: ax.text(-0.15, 1.15, r'$(b)$',transform=ax.transAxes,fontsize=letters_fontsize, fontweight='bold', va='top', ha='left') #plt.tight_layout(rect=[0, 0.05, 1, 1]) #plt.show() plt.savefig('FIG4.pdf', format='pdf', bbox_inches='tight', dpi=600)
Out[32]:
In [0]: