Skip to content

Commit b4facd9

Browse files
author
Paul Mannix
committed
Updated figures
1 parent d12270e commit b4facd9

File tree

4 files changed

+155
-31
lines changed

4 files changed

+155
-31
lines changed

‎Gap_Continuation.py‎

Lines changed: 7 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -34,9 +34,9 @@ def Gap_Vary(Ra_s_new, open_filename, frame):
3434
print("\n Loading Ra = %e, Ra_s=%e, Pr=%2.5f, Tau=%2.5f, d=%2.5f and resolution N_fm, N_r = %d,%d \n"%(Ra,Ra_s,Pr,Tau,d,N_fm,N_r))
3535

3636
sign = -1
37-
N_steps = 250
37+
N_steps = 500
3838
Y = np.hstack((X, Ra))
39-
kwargs = {"Ra":Ra,"Ra_s":Ra_s_new,"Tau":Tau,"Pr":Pr,"d":d,"N_fm":N_fm_n,"N_r":N_r_n, "symmetric":False}
39+
kwargs = {"Ra":Ra,"Ra_s":Ra_s_new,"Tau":Tau,"Pr":Pr,"d":d,"N_fm":N_fm_n,"N_r":N_r_n, "symmetric":True}
4040

4141
# Generate a new path
4242
filename, extension = os.path.splitext(open_filename)
@@ -53,13 +53,12 @@ def main():
5353

5454
print('Creating a test directory .... \n')
5555

56-
filenames = ["NewtonSolveRas400d0.345_0.h5", "TimeStepd0.34_0.h5", "TimeStepd0.335_0.h5", "TimeStepd0.33_0.h5"]
57-
frame = -1
58-
Ra_s = 400
56+
open_filename = "ConvectonL10PlusRas400Ras300_1.h5"
57+
frame = 0
5958

60-
for open_filename in filenames:
61-
print("\n file = ",open_filename)
62-
Gap_Vary(Ra_s, open_filename, frame)
59+
for Ras_i in [306.125]:
60+
print("\n Ra_s = ",Ras_i)
61+
Gap_Vary(Ras_i, open_filename, frame)
6362

6463
return None
6564

‎Main.py‎

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -427,7 +427,7 @@ def Time_Step(open_filename=None, frame=-1, d_new=0.31325, Ra_new = 3750.00):
427427
return filename;
428428

429429

430-
def _Newton(X, Ra, Ra_s, Tau, Pr, d, N_fm, N_r, symmetric, dt=.1, tol_newton=1e-08, tol_gmres=1e-04, Krylov_Space_Size=300):
430+
def _Newton(X, Ra, Ra_s, Tau, Pr, d, N_fm, N_r, symmetric, dt=10**4, tol_newton=1e-08, tol_gmres=1e-04, Krylov_Space_Size=300):
431431

432432
"""
433433
Given a starting point and a control parameter compute a new steady-state using Newton iteration
@@ -738,7 +738,7 @@ def _NewtonC(Y, sign, ds, **kwargs_f):
738738
return Y ,sign,ds, Norm,KE,NuT,NuS, exitcode;
739739

740740

741-
def _ContinC(Y_0_dot, Y_0, sign, ds, Ra, Ra_s, Tau, Pr, d, N_fm, N_r, symmetric, dt=.1, tol_newton=1e-08, tol_gmres=1e-04, Krylov_Space_Size=300):
741+
def _ContinC(Y_0_dot, Y_0, sign, ds, Ra, Ra_s, Tau, Pr, d, N_fm, N_r, symmetric, dt=10**4, tol_newton=1e-08, tol_gmres=1e-04, Krylov_Space_Size=300):
742742

743743
"""
744744
Given a starting point and a control parameter compute a new steady-state using Newton iteration
@@ -966,7 +966,7 @@ def _Continuation(filename, N_steps, sign, Y, **kwargs):
966966
# Default parameters
967967
ds =0.01; # The starting step size
968968
ds_min=1.0; #Threshold to switching between Newton & Psuedo
969-
ds_max=10.0; # Max Newton step
969+
ds_max=1.5; # Max Newton step
970970

971971
Result = result()
972972

@@ -1090,12 +1090,12 @@ def Continuation(open_filename, frame=-1):
10901090

10911091
# ~~~~~~~~~ Interpolate ~~~~~~~~~~~~~~~~~~~
10921092
from Matrix_Operators import INTERP_RADIAL,INTERP_THETAS
1093-
N_r_n = N_r; X = INTERP_RADIAL(N_r_n,N_r,X,d);
1094-
N_fm_n = 64; X = INTERP_THETAS(N_fm_n,N_fm,X);
1093+
N_r_n = 25; X = INTERP_RADIAL(N_r_n,N_r,X,d);
1094+
N_fm_n = 128; X = INTERP_THETAS(N_fm_n,N_fm,X);
10951095
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
10961096

1097-
sign = 1;
1098-
N_steps= 250;
1097+
sign = -1;
1098+
N_steps= 500;
10991099
Y = np.hstack( (X,Ra) );
11001100
kwargs = {"Ra":Ra,"Ra_s":Ra_s,"Tau":Tau,"Pr":Pr,"d":d,"N_fm":N_fm_n,"N_r":N_r_n, "symmetric":True}
11011101

@@ -1198,8 +1198,8 @@ def _plot_bif(filename, point = -1):
11981198
print("Initialising the code for running...")
11991199

12001200
# %%
1201-
#file = "AntiConvectonL10PlusRas400_3.h5"
1202-
#Continuation(open_filename=file,frame=-3)
1201+
#file = "ConvectonL10PlusRas303_6.h5"
1202+
#Continuation(open_filename=file,frame=18)
12031203
#trim(filename='Continuationl11Ras150_0.h5',point=-4)
12041204
#_plot_bif(filename='Continuationl11Ras150_0.h5',point=-4) # Good start point
12051205

@@ -1209,9 +1209,9 @@ def _plot_bif(filename, point = -1):
12091209

12101210
# %%
12111211
from Plot_Tools import Cartesian_Plot, Energy,Uradial_plot
1212-
filename = "AntiConvectonL10PlusRas400_4.h5"
1213-
_plot_bif(filename,point=-3)
1214-
Cartesian_Plot(filename,frame=1,Include_Base_State=False)
1212+
filename = "ConvectonL10PlusRas303_7.h5"
1213+
_plot_bif(filename,point=-1)
1214+
#Cartesian_Plot(filename,frame=-1,Include_Base_State=False)
12151215
#Energy(filename,frame=-1)
12161216
# %%
12171217

‎Paper_Figures/Plot_figures_L10_Plus_compare.py‎

Lines changed: 58 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@ def add_to_fig(obj):
6262
N_fm_fold = []
6363
N_r_fold = []
6464
for filename in glob.glob(folder + '/*.h5'):
65-
65+
6666
obj = result()
6767
with h5py.File(filename, 'r') as f:
6868
ff = f["Bifurcation"]
@@ -130,6 +130,32 @@ def Add_Label(X_folds, Nr_folds, Nfm_folds, Ra_folds, ax):
130130

131131
# %%
132132

133+
134+
# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~
135+
# L = 10 Large
136+
# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~
137+
138+
# fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(16, 6), layout='constrained', sharey=True)
139+
140+
141+
# # A) Plot the bifurcation diagram
142+
143+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras303_Convectons/Upper/', ax, line='c-')
144+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras303_Convectons/Lower/', ax, line='k-')
145+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras303_AntiConvectons/', ax, line='k:')
146+
# ax.set_title(r'$Ra_S = 303$', fontsize=25)
147+
148+
# ax.set_xlabel(r'$Ra_T$', fontsize=25)
149+
# ax.tick_params(axis='both', labelsize=25)
150+
# ax.set_xlim([2700, 3750])
151+
# ax.set_ylim([0, 7])
152+
153+
# #plt.savefig('Bifurcation_L10Plus_Ras_Compare_zoom.png', format='png', dpi=100)
154+
# plt.show()
155+
156+
157+
158+
# %%
133159
# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~
134160
# L = 10 Large
135161
# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~
@@ -153,28 +179,49 @@ def Add_Label(X_folds, Nr_folds, Nfm_folds, Ra_folds, ax):
153179
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras200_AntiConvectons/', axs[i], line='k:') # Low Res
154180
axs[i].set_title(r'$Ra_S = 200$', fontsize=25)
155181

156-
i = 2
157-
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras250_Convectons/Upper/', axs[i], line='c-')
158-
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras250_Convectons/Lower/', axs[i], line='k-')
159-
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras250_AntiConvectons/', axs[i], line='k:')
160-
axs[i].set_title(r'$Ra_S = 250$', fontsize=25)
182+
# i = 2
183+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras250_Convectons/Upper/', axs[i], line='c-')
184+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras250_Convectons/Lower/', axs[i], line='k-')
185+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras250_AntiConvectons/', axs[i], line='k:')
186+
# axs[i].set_title(r'$Ra_S = 250$', fontsize=25)
161187

162-
i = 3
188+
i = 2
163189
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras300_Convectons/Upper/', axs[i], line='c-')
164190
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras300_Convectons/Lower/', axs[i], line='k-')
165191
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras300_AntiConvectons/', axs[i], line='k:')
166192
axs[i].set_title(r'$Ra_S = 300$', fontsize=25)
167193

194+
195+
i = 3
196+
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras303_Convectons/Upper/', axs[i], line='c-')
197+
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras303_Convectons/Lower/', axs[i], line='k-')
198+
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras303_AntiConvectons/', axs[i], line='k:')
199+
axs[i].set_title(r'$Ra_S = 303$', fontsize=25)
200+
201+
168202
i = 4
169-
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras400_Convectons/', axs[i], line='k-')
170-
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras400_AntiConvectons/', axs[i], line='k:')
171-
axs[i].set_title(r'$Ra_S = 400$', fontsize=25)
203+
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras306.125_Convectons/Upper/', axs[i], line='c-')
204+
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras306.125_Convectons/Lower/', axs[i], line='k-')
205+
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras306.125_AntiConvectons/', axs[i], line='k:')
206+
axs[i].set_title(r'$Ra_S = 306.125$', fontsize=25)
207+
208+
# i = 5
209+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras325_Convectons/Upper/', axs[i], line='c-')
210+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras325_Convectons/Lower/', axs[i], line='k-')
211+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras325_AntiConvectons/', axs[i], line='k:')
212+
# axs[i].set_title(r'$Ra_S = 325$', fontsize=25)
213+
214+
# i = 6
215+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras400_Convectons/Upper/', axs[i], line='c-')
216+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras400_Convectons/Lower/', axs[i], line='k-')
217+
# X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Ras400_AntiConvectons/', axs[i], line='k:')
218+
# axs[i].set_title(r'$Ra_S = 400$', fontsize=25)
172219

173220
for ax in axs:
174221
ax.set_xlabel(r'$Ra_T$', fontsize=25)
175222
ax.tick_params(axis='both', labelsize=25)
176223
ax.set_xlim([2700, 3750])
177-
ax.set_ylim([0, 7])
224+
ax.set_ylim([0, 6])
178225

179226
plt.savefig('Bifurcation_L10Plus_Ras_Compare.png', format='png', dpi=100)
180227
plt.show()

‎Paper_Figures/Plot_figures_l10_snaking.py‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -483,3 +483,81 @@ def Add_Label(X_folds, Nr_folds, Nfm_folds, Ra_folds, ax):
483483

484484
plt.savefig('L10_Plus_AntiConvectons_Ras175.png', format='png', dpi=100)
485485
plt.show()
486+
487+
488+
489+
490+
# %%
491+
# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~
492+
# (L^{A+}_10) L = 10 Plus Convectons Ra_s = 175 (Blue branch)
493+
# ~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~
494+
495+
fig, axs = plt.subplots(nrows=7, ncols=2, figsize=(16, 6), layout='constrained')
496+
497+
# remove the underlying Axes in the right column
498+
gs = axs[1, 0].get_gridspec()
499+
for axl in axs[:, 0]:
500+
axl.remove()
501+
ax = fig.add_subplot(gs[:, 0])
502+
503+
# A) Plot the bifurcation diagram
504+
X_folds, Nr_folds, Nfm_folds, Ra_folds = Plot_full_bif(dir + 'Convectons_Plus_Ras175/', ax, line='c-')
505+
506+
ax.set_ylabel(r'$\mathcal{E}$', fontsize=25)
507+
ax.set_xlabel(r'$Ra_T$', fontsize=25)
508+
ax.tick_params(axis='both', labelsize=25)
509+
ax.set_ylim([1, 4])
510+
ax.set_xlim([2940, 3040])
511+
ax.tick_params(axis='both', labelsize=25)
512+
513+
514+
# L=10 Minus eigenvector from KE = 0
515+
obj = result()
516+
with h5py.File(dir + 'Convectons_Plus_Ras175/' + "ConvectonL10PlusRas175_3.h5", 'r') as f:
517+
ff = f["Bifurcation"]
518+
for key in ff.keys():
519+
setattr(obj, key, ff[key][()])
520+
N_fm = f['Parameters']["N_fm"][()]
521+
N_r = f['Parameters']["N_r"][()]
522+
523+
D, R = cheb_radial(N_r, d)
524+
525+
point = 25
526+
KE = Kinetic_Energy(obj.X_DATA[point], R, D, N_fm, N_r-1, symmetric=False)
527+
ax.plot(obj.Ra_DATA[point], KE, 'bs')
528+
529+
X_folds.insert(5, obj.X_DATA[point])
530+
Ra_folds.insert(5, obj.Ra_DATA[point])
531+
Nr_folds.insert(5, N_r)
532+
Nfm_folds.insert(5, N_fm)
533+
534+
obj = result()
535+
with h5py.File(dir + 'Convectons_Plus_Ras175/' + "ConvectonL10PlusRas175_1.h5", 'r') as f:
536+
ff = f["Bifurcation"]
537+
for key in ff.keys():
538+
setattr(obj, key, ff[key][()])
539+
N_fm = f['Parameters']["N_fm"][()]
540+
N_r = f['Parameters']["N_r"][()]
541+
542+
D, R = cheb_radial(N_r, d)
543+
544+
point = 37
545+
KE = Kinetic_Energy(obj.X_DATA[point], R, D, N_fm, N_r-1, symmetric=False)
546+
ax.plot(obj.Ra_DATA[point], KE, 'bs')
547+
548+
X_folds.insert(0, obj.X_DATA[point])
549+
Ra_folds.insert(0, obj.Ra_DATA[point])
550+
Nr_folds.insert(0, N_r)
551+
Nfm_folds.insert(0, N_fm)
552+
553+
# D) Plot the points out in terms of their stream-function
554+
Psi_Plot(X_folds, Nr_folds, Nfm_folds, axs=axs[:, 1])
555+
556+
for count in range(6+1):
557+
axs[::-1, 1][count].annotate(r'%d' % count, xy=(-0.05, 0.5), xycoords='axes fraction', fontsize=20)
558+
559+
# E) Add labels
560+
Add_Label(X_folds[::-1], Nr_folds[::-1], Nfm_folds[::-1], Ra_folds[::-1], ax)
561+
562+
plt.savefig('L10_Plus_Convectons_Ras175.png', format='png', dpi=100)
563+
plt.show()

0 commit comments

Comments
 (0)