figure 11.ipynb

Kernel: Python 3 (ipykernel)

import numpy as np
import h5py
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import cmocean

# Aesthetics
sns.set_theme()
sns.set_style('ticks')
sns.set_context('paper')
plt.rc('mathtext', fontset='stix')
plt.rc('font', family='serif')

Define consistent colours and line styles for variations in $Re$, $Gr$, and $Pr$

camp = sns.color_palette('flare', as_cmap=True)
cspeed = cmocean.tools.crop_by_percent(cmocean.cm.tempo, 30, which='both', N=None)

def Re_col(lR):
    return camp((lR - 2.5)/1.5)

def Gr_col(lG):
    return cspeed((lG - 6)/2)

def Pr_stl(Pr):
    if Pr==1:
        stl = '-'
    elif Pr==4:
        stl = '--'
    elif Pr==10:
        stl = ':'
    else:
        stl = '-.'
    return stl

Average profiles using symmetry or anti-symmetry across the midplane

def sym_prof(v):
    n = v.size
    return 0.5*(v[:n//2] + v[-1:-n//2-1:-1])
def asym_prof(v):
    n = v.size
    return 0.5*(v[:n//2] - v[-1:-n//2-1:-1])

fig, axs = plt.subplots(2,2, figsize=(5.2,3.2), layout='constrained', dpi=200, sharex='col', sharey='col')

with h5py.File('../data/profile_record.h5','r') as fp:

    for grp in fp.__iter__():
        # Read mean profiles
        x = fp[grp+'/xm'][:]
        # xs = x[:x.size//2]
        xr = fp[grp+'/xmr'][:]
        xsr = xr[:xr.size//2]
        Tbar = -asym_prof(fp[grp+'/Sbar'][:])
        vybar = asym_prof(fp[grp+'/vybar'][:])
        vzbar = sym_prof(fp[grp+'/vzbar'][:])

        # Read control parameters (logarithms)
        lGr = int(grp[2])
        lRe = float(grp[-4:])
        lRi = lGr - 2*lRe

        # Exclude cases with no pressure gradient
        if lRe > 0:
            # Compute free-fall and bulk velocity scales
            Uf = 10**np.min([lGr/2 - lRe, 0])
            Ub = 10**np.min([lRe - lGr/2, 0])
    
            Gr = 10**lGr
            Re = 10**lRe
            Pr = int(grp[6:-7])

            # Compute dimensionless friction velocities and total shear Reynolds number
            Wtau = (vzbar[0]/x[0] / np.max([Gr**0.5,Re]) )**0.5
            Retauz2 = Wtau*np.max([Gr**0.5,Re])            
            Vtau = (np.abs(vybar)[0]/x[0] / np.max([Gr**0.5,Re]) )**0.5
            Retauy2 = Vtau*np.max([Gr**0.5,Re])
            Retau = (Retauz2**4 + Retauy2**4)**0.25

            # Compute Nusselt number from boundary gradient to give friction temperature scale Thet
            Nu = (0.5 - Tbar[0])/xr[0]
            Thet = Nu/Retau
    
            n = xr.size
            if Gr*Pr==1e7:
                axs[0,0].plot(xsr, Tbar, color=Re_col(lRe), linestyle=Pr_stl(Pr))
                axs[0,1].plot(xsr*Retau, (0.5 - Tbar)/Thet, color=Re_col(lRe), linestyle=Pr_stl(Pr))
            if lRi==0:
                axs[1,0].plot(xsr, Tbar, color=Gr_col(lGr), linestyle=Pr_stl(Pr))
                axs[1,1].plot(xsr*Retau, (0.5 - Tbar)/Thet, color=Gr_col(lGr), linestyle=Pr_stl(Pr))
        
# Add theoretical comparison
xl = 10**np.linspace(-2,1.2)
for ax in axs[:,1]:
    ax.plot(xl, xl, 'k--')
    ax.semilogx()
axs[1,0].set_xlabel('$y/H$')
axs[1,0].set_xlim([0,0.5])
axs[1,0].set_ylim([0,0.5])

# Add labels
for ax in axs[:,0]:
    ax.set_ylabel('$\overline{\\theta}/\Delta$')
for ax in axs[:,1]:
    ax.set_ylabel(r'$Pr\ \overline{\theta}^+$')
axs[1,1].set(
    xlabel='$y^+$',
    xlim=[1e-1,400],
    ylim=[0,20]
)

# Add standalone colorbar for Re
norm = mpl.colors.Normalize(vmin=2.5, vmax=4)
cbR = fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=camp),
             ax=axs[0,:], label='$Re\ (Ra=10^7)$')
Rlst = [2.5, 3, 3.5, 4]
cbR.set_ticks(Rlst, labels=['$10^{'+str(R)+'}$' for R in Rlst])

# Add standalone colorbar for Gr
norm = mpl.colors.Normalize(vmin=6, vmax=8)
cbG = fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cspeed),
             ax=axs[1,:], label='$Gr\ (Ri=1)$')
Glst = [6, 7, 8]
cbG.set_ticks(Glst, labels=['$10^{'+str(G)+'}$' for G in Glst])

# Add legend
ax = axs[1,0]
Prs = [1, 4, 10]
for P in Prs:
    ax.plot(0,0, color=Gr_col(6), linestyle=Pr_stl(P), label='$Pr=%i$' % P)
ax.legend()

for ax in axs.flatten():
    ax.grid(True)

lbls = ['$(a)$','$(b)$','$(c)$','$(d)$']
for i, lbl in enumerate(lbls):
    axs[i//2,i%2].annotate(lbl, (0.02, 0.97), xycoords='axes fraction', ha='left', va='top')

ax = axs[1,1]
ax.annotate('$Pr\ \\theta^+ = y^+$', (1.0,3.5),ha='center', va='center', rotation=25)

# fig.savefig('T_profiles.pdf')

plt.show()