TASK 2 implementation¶
In [ ]:
# importing the required modules
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
fig_path = 'figures/task2/'
fig_save = True
Required functions¶
In [ ]:
# defining spatial deritivative functions
def spatial_derivative_B(B, dx, order = 6, ghost_zone_type = 'relative anti-symmetric'):
'''
This function calculates the spatial derivatives of the magnetic field B.
Parameters:
B: 2D array
The magnetic field B.
dx: float
The spatial step size.
order: int
The order of the finite difference scheme. It can be 2, 4, 6, 8 or 10.
Default is 6.
ghost_zone_type: str
The type of the ghost zone. It can be 'anti-symmetric', 'symmetric', 'relative anti-symmetric', 'smooth symmetric', or None.
Default is 'relative anti-symmetric'.
Returns:
dBdx: 2D array
The spatial-derivative of the magnetic field B.
ddBddx: 2D array
The second z-derivative of the magnetic field B.
'''
dBdx = np.zeros(B.shape)
ddBddx = np.zeros(B.shape)
# print(f"The shape of B is {B.shape}")
# print(f"The shape of dBdx is {dBdx.shape}")
# print(f"The shape of ddBddx is {ddBddx.shape}")
B = np.copy(B) # to avoid modifying the original array
if order not in ['np',2, 4, 6, 8, 10]:
raise ValueError("Order should be 'np',2, 4, 6, 8 or 10")
if ghost_zone_type not in ['anti-symmetric', 'symmetric', 'relative anti-symmetric', 'smooth',None]:
raise ValueError('ghost_zone_type should be anti-symmetric, symmetric, relative anti-symmetric or None')
if order == 'np':
ghost_zone = 3 # 3 cells on each side
else:
ghost_zone = int(order/2)
if ghost_zone_type == 'relative anti-symmetric':
B = np.pad(B, ((0, 0), (ghost_zone, ghost_zone)), mode='constant')
for i in range(ghost_zone):
B[:,i] = (2 * B[:,ghost_zone]) - B[:,ghost_zone + ghost_zone - i] # correcting the start of the array
B[:, -1 - i] = 2 * B[:,-ghost_zone - 1] - B[:, (- ghost_zone - 1) - ghost_zone + i] # correcting the end of the array
elif ghost_zone_type == 'anti-symmetric':
B = np.pad(B, ((0, 0), (ghost_zone, ghost_zone)), mode='reflect') # makes symmetric padding
for i in range(ghost_zone):
B[:,i] = - B[:, i] # making anti-symmetric at the start of the array
B[:, -1 - i] = - B[:, -1 - i] # making anti-symmetric at the end of the array
elif ghost_zone_type == 'symmetric':
B = np.pad(B, ((0, 0), (ghost_zone, ghost_zone)), mode='reflect')
elif ghost_zone_type == 'smooth':
B = np.pad(B, ((0, 0), (ghost_zone, ghost_zone)), mode='reflect', reflect_type='odd')
else:
B = np.pad(B, ((0, 0), (ghost_zone, ghost_zone)), mode='constant') # zero padding
# print(f"The shape of B after padding is {B.shape}")
# print(f"The shape of dBdx after padding is {dBdx.shape}")
# print(f"The shape of ddBddx after padding is {ddBddx.shape}")
if order == 6:
for i in range(ghost_zone, B.shape[1] - ghost_zone):
dBdx[:,i-ghost_zone] = ((- B[:,i - 3]) + (9 * B[:,i - 2]) - (45 * B[:,i - 1]) + (45 * B[:,i + 1]) - (9 * B[:,i + 2]) + (B[:,i + 3])) / (60 * dx)
ddBddx[:,i-ghost_zone] = ((2 * B[:,i - 3]) - (27 * B[:,i - 2]) + (270 * B[:,i - 1]) - (490 * B[:, i]) + (270 * B[:,i + 1]) - (27 * B[:,i + 2]) + (2 * B[:,i + 3])) / (180 * (dx ** 2))
elif order == 2:
for i in range(ghost_zone, B.shape[1] - ghost_zone):
dBdx[:,i-ghost_zone] = ((- B[:,i - 1]) + (B[:,i + 1])) / (2 * dx)
ddBddx[:,i-ghost_zone] = (B[:,i - 1] - (2 * B[:, i]) + B[:,i + 1]) / (dx ** 2)
elif order == 4:
for i in range(ghost_zone, B.shape[1] - ghost_zone):
dBdx[:,i-ghost_zone] = ((B[:,i - 2]) - (8 * B[:,i - 1]) + (8 * B[:,i + 1]) - (B[:,i + 2])) / (12 * dx)
ddBddx[:,i-ghost_zone] = ((- B[:,i - 2]) + (16 * B[:,i - 1]) - (30 * B[:,i]) + (16 * B[:,i + 1]) - (B[:,i + 2])) / (12 * (dx ** 2))
elif order == 8:
for i in range(ghost_zone, B.shape[1] - ghost_zone):
dBdx[:,i-ghost_zone] = ((3*B[:,i - 4]) - (32 * B[:,i - 3]) + (168 * B[:,i - 2]) - (672 * B[:,i - 1]) + (672 * B[:,i + 1]) - (168 * B[:,i + 2]) + (32 * B[:,i + 3]) - (3 * B[:,i + 4])) / (840 * dx)
ddBddx[:,i-ghost_zone] = ((-9 * B[:,i - 4]) + (128 * B[:,i - 3]) - (1008 * B[:,i - 2]) + (8064 * B[:,i - 1]) - (14350 * B[:,i]) + (8064 * B[:,i + 1]) - (1008 * B[:,i + 2]) + (128 * B[:,i + 3]) - (9 * B[:,i + 4])) / (5040 * (dx ** 2))
elif order == 10:
for i in range(ghost_zone, B.shape[1] - ghost_zone):
dBdx[:,i-ghost_zone] = ((-2 * B[:,i - 5]) + (25 * B[:,i - 4]) - (150 * B[:,i - 3]) + (600 * B[:,i - 2]) - (2100 * B[:,i - 1]) + (2100 * B[:,i + 1]) - (600 * B[:,i + 2]) + (150 * B[:,i + 3]) - (25 * B[:,i + 4]) + (2 * B[:,i + 5])) / (2520 * dx)
ddBddx[:,i-ghost_zone] = ((8 * B[:,i - 5]) - (125 * B[:,i - 4]) + (1000 * B[:,i - 3]) - (6000 * B[:,i - 2]) + (42000 * B[:,i - 1]) - (73766 * B[:, i]) + (42000 * B[:,i + 1]) - (6000 * B[:,i + 2]) + (1000 * B[:,i + 3]) - (125 * B[:,i + 4]) + (8 * B[:,i + 5])) / (25200 * (dx ** 2))
else:
dBdx = np.gradient(B, dx, axis=1)
ddBddx = np.gradient(dBdx, dx, axis=1)
# removing the ghost zones
dBdx = dBdx[:,ghost_zone:-ghost_zone]
ddBddx = ddBddx[:,ghost_zone:-ghost_zone]
# update with other orders
return dBdx, ddBddx
In [ ]:
def first_derivative_dYdx_FD(Y, dx, order, ghost_zone_type):
dYdx = np.zeros(Y.shape)
Y = np.copy(Y) # to avoid modifying the original array
if order not in ['np',2, 4, 6, 8, 10]:
raise ValueError("Order should be 'np',2, 4, 6, 8 or 10")
if ghost_zone_type not in ['anti-symmetric', 'symmetric', 'relative anti-symmetric','smooth', None]:
raise ValueError('ghost_zone_type should be anti-symmetric, symmetric, relative anti-symmetric or None')
if order == 'np':
ghost_zone = 3 # 3 cells on each side
else:
ghost_zone = int(order/2)
if ghost_zone_type == 'relative anti-symmetric':
Y = np.pad(Y, (ghost_zone, ghost_zone), mode='constant')
for i in range(ghost_zone):
Y[i] = (2 * Y[ghost_zone]) - Y[ghost_zone + ghost_zone - i] # correcting the start of the array
Y[ -1 - i] = 2 * Y[-ghost_zone - 1] - Y[ (- ghost_zone - 1) - ghost_zone + i] # correcting the end of the array
elif ghost_zone_type == 'symmetric':
Y = np.pad(Y, (ghost_zone, ghost_zone), mode='reflect')
elif ghost_zone_type == 'anti-symmetric':
Y = np.pad(Y, (ghost_zone, ghost_zone), mode='reflect') # makes symmetric padding
for i in range(ghost_zone):
Y[i] = - Y[ i] # making anti-symmetric at the start of the array
Y[ -1 - i] = - Y[ -1 - i] # making anti-symmetric at the end of the array
elif ghost_zone_type == 'smooth':
Y = np.pad(Y, (ghost_zone, ghost_zone), mode='reflect', reflect_type='odd')
else:
Y = np.pad(Y, (ghost_zone, ghost_zone), mode='constant') # zero padding
# print(f"shape of Y is {Y.shape}")
if order == 6:
for i in range(ghost_zone, len(Y) - ghost_zone):
dYdx[i-ghost_zone] = ((- Y[i - 3]) + (9 * Y[i - 2]) - (45 * Y[i - 1]) + (45 * Y[i + 1]) - (9 * Y[i + 2]) + (Y[i + 3])) / (60 * dx)
elif order == 2:
for i in range(ghost_zone, len(Y) - ghost_zone):
dYdx[i-ghost_zone] = ((- Y[i - 1]) + (Y[i + 1])) / (2 * dx)
elif order == 4:
for i in range(ghost_zone, len(Y) - ghost_zone):
dYdx[i-ghost_zone] = ((Y[i - 2]) - (8 * Y[i - 1]) + (8 * Y[i + 1]) - (Y[i + 2])) / (12 * dx)
elif order == 8:
for i in range(ghost_zone, len(Y) - ghost_zone):
dYdx[i-ghost_zone] = ((3*Y[i - 4]) - (32 * Y[i - 3]) + (168 * Y[i - 2]) - (672 * Y[i - 1]) + (672 * Y[i + 1]) - (168 * Y[i + 2]) + (32 * Y[i + 3]) - (3 * Y[i + 4])) / (840 * dx)
elif order == 10:
for i in range(ghost_zone, len(Y) - ghost_zone):
dYdx[i-ghost_zone] = ((-2 * Y[i - 5]) + (25 * Y[i - 4]) - (150 * Y[i - 3]) + (600 * Y[i - 2]) - (2100 * Y[i - 1]) + (2100 * Y[i + 1]) - (600 * Y[i + 2]) + (150 * Y[i + 3]) - (25 * Y[i + 4]) + (2 * Y[i + 5])) / (2520 * dx)
else:
dYdx = np.gradient(Y, dx, axis=1)
# removing the ghost zones
dYdx = dYdx[ghost_zone:-ghost_zone]
# update with other orders
return dYdx
In [ ]:
def simulate_B_field(r_i, r_f, Nr, T, Nt, B_r0, B_phi0, dBrdt, dBphidt, order = 6, ghost_zone_type = 'relative anti-symmetric', iterations_to_save = None):
'''
This function simulates the time evolution of the magnetic field B.
Parameters:
r_i: float
The initial value of the radial coordinate.
r_f: float
The final value of the radial coordinate.
Nr: int
The number of spatial points.
T: float
The final time.
Nt: int
The number of time points.
B_r0: function
The initial condition for the radial component of the magnetic field B.
B_phi0: function
The initial condition for the azimuthal component of the magnetic field B.
dBrdt: function
The function defining the time derivative of the radial component of the magnetic field B.
dBphidt: function
The function defining the time derivative of the azimuthal component of the magnetic field B.
order: int
The order of the finite difference scheme. It can be 2, 4, 6, 8, 10 or 'np'.
Default is 6.
ghost_zone_type: str
The type of the ghost zone. It can be 'anti-symmetric', 'symmetric', 'relative anti-symmetric' or None.
Default is 'relative anti-symmetric'.
iterations_to_save: list
The list of time iterations to save. If None, only the initial and final states are saved.
If 'all', all iterations are saved.
Default is None.
Returns:
B_list: 3D array
The solution of the magnetic field B in time.
shape: (len(iterations_to_save), 2, Nr)
'''
dr = (r_f - r_i) / (Nr - 1) # Spatial step size
dt = T / Nt # Temporal step size
# Create arrays to store the solution in time
B_list = []
# creating initial conditions from z, B_r0 and B_phi0
B = np.zeros((2, Nr))
r = np.linspace(r_i, r_f, Nr)
B[0, :] = B_r0(r, 0, 0) # first row is B_r
B[1, :] = B_phi0(r, 0, 0) # second row is B_phi
if iterations_to_save is None:
iterations_to_save = [0, Nt]
elif iterations_to_save == 'all':
iterations_to_save = list(range(Nt + 1))
if 0 in iterations_to_save:
B_list.append(np.copy(B))
#--------- RK4 time-stepping block ---------#
for n in tqdm(range(1, Nt + 1)):
# Compute spatial derivatives
dBdr, ddBddr = spatial_derivative_B(B, dr, order, ghost_zone_type)
k1_r = dBrdt(B[0], B[1], dBdr[0], ddBddr[0], (n - 1) * dt)
k1_phi = dBphidt(B[0], B[1], dBdr[1], ddBddr[1], (n - 1) * dt)
k2_r = dBrdt(B[0] + (0.5 * dt * k1_r), B[1] + (0.5 * dt * k1_phi), dBdr[0], ddBddr[0], ((n-0.5) * dt))
k2_phi = dBphidt(B[0] + (0.5 * dt * k1_r), B[1] + (0.5 * dt * k1_phi), dBdr[1], ddBddr[1], ((n-0.5) * dt))
k3_r = dBrdt(B[0] + (0.5 * dt * k2_r), B[1] + (0.5 * dt * k2_phi), dBdr[0], ddBddr[0], ((n-0.5) * dt))
k3_phi = dBphidt(B[0] + (0.5 * dt * k2_r), B[1] + (0.5 * dt * k2_phi), dBdr[1], ddBddr[1], ((n-0.5) * dt))
k4_r = dBrdt(B[0] + (dt * k3_r), B[1] + (dt * k3_phi), dBdr[0], ddBddr[0], n * dt)
k4_phi = dBphidt(B[0] + (dt * k3_r), B[1] + (dt * k3_phi), dBdr[1], ddBddr[1], n * dt)
B[0, :] += (dt / 6) * (k1_r + (2 * k2_r) + (2 * k3_r) + k4_r)
B[1, :] += (dt / 6) * (k1_phi + (2 * k2_phi) + (2 * k3_phi) + k4_phi)
# print(B)
if n in iterations_to_save:
B_list.append(np.copy(B))
#------------------------------------------#
B_list = np.array(B_list)
return B_list
Primary implementation and evolution of components¶
In [ ]:
# defining constants and parameters
eta = 1
r_i = 0.01 # start of spatial domain
r_f = 8.01 # end of spatial domain
T = 1.0 # Total simulation time
Nr = 100 # Number of spatial grid points
Nt = 4000 # Number of time steps
h = 1
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
iterations_to_save = 'all' # which time steps to save
iterations_to_save_plot = [int(i * Nt/8) for i in range(8 + 1)] # which time steps to plot
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
omega_0 = 50
r_omega = 1
h = 1
alpha = 0.5
alpha_cap = 1 #UPDATE
# defining the spatial dependent omega
omega = omega_0 / (np.sqrt(1 + (r / r_omega) ** 2))
# defining q
neg_q_omega = r * first_derivative_dYdx_FD(omega, dr, order = 6, ghost_zone_type='smooth')
R_alpha = alpha * h / eta
R_omega = (neg_q_omega * (h ** 2)) / eta
D = R_alpha * R_omega
# defining initial conditions
def B_phi0(r, phi, z):
return (3 * np.sin(((r - r_i)/(r_f - r_i))*np.pi)) + (1 * np.sin(((r - r_i)/(r_f - r_i))*2*np.pi)) + (2 * np.sin(((r - r_i)/(r_f - r_i))*5*np.pi))
def B_r0(r, phi, z):
return (3 * np.sin(((r - r_i)/(r_f - r_i))*4*np.pi)) + (1 * np.sin(((r - r_i)/(r_f - r_i))*2*np.pi)) + (2 * np.sin(((r - r_i)/(r_f - r_i))*3*np.pi))
# definiing the time evolution of the magnetic field Br
def dBrdt(Br, Bphi, dBrdr, ddBrdr, t):
return (- (2/(np.pi*h)) * R_alpha * alpha_cap* Bphi) + (eta * (ddBrdr + (dBrdr / r) - (Br/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Br)))
# definiing the time evolution of the magnetic field Bphi
def dBphidt(Br, Bphi, dBphidr, ddBphidr, t):
return (R_omega * Br) + (eta * (ddBphidr + (dBphidr / r) - (Bphi/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Bphi)))
In [ ]:
# plotting the fynamo number vs r
plt.figure(figsize=(12, 6))
plt.plot(r, D)
plt.xlabel('r')
plt.ylabel('Dynamo number')
plt.title('Dynamo number vs r')
if fig_save:
plt.savefig(fig_path + 'Dynamo_number_vs_r.png')
plt.show()
# plotting omega vs r
plt.figure(figsize=(12, 6))
plt.plot(r, omega)
plt.xlabel('r')
plt.ylabel(r'$\omega$')
plt.title(r'$\omega$ vs r')
if fig_save:
plt.savefig(fig_path + 'omega_vs_r.png')
plt.show()
In [ ]:
# running the simulations
B_list = simulate_B_field(r_i, r_f, Nr, T, Nt, B_r0, B_phi0, dBrdt, dBphidt, order = order, ghost_zone_type = ghost_zone_type, iterations_to_save = iterations_to_save)
100%|██████████| 4000/4000 [00:15<00:00, 254.90it/s]
In [ ]:
# plotting Br and Bphi time evolition in two plots side by side
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
for i, ax in enumerate(axs):
for j in range(len(iterations_to_save_plot)):
ax.plot(r, B_list[iterations_to_save_plot[j], i], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
ax.set_xlabel('r')
if i == 0:
name = r'$B_r$'
else:
name = r'$B_{\phi}$'
ax.set_ylabel(name)
ax.set_title(f'Time evolution of {name}')
ax.legend()
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = {ghost_zone_type}\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\alpha = {alpha}, \\hat\\alpha = {alpha_cap}$\n'
figure_text3 = 'seed fields: $B_r = 3 \\sin(\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(1.5\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
if fig_save:
plt.savefig(fig_path + 'Br_Bphi_time_evolution.png')
plt.show()
In [ ]:
# plotting magnetic field strength time evolution
B_strength = np.sqrt(B_list[:, 0] ** 2 + B_list[:, 1] ** 2)
plt.figure(figsize=(12, 6))
for j in range(len(iterations_to_save_plot)):
plt.plot(r, B_strength[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
plt.xlabel('r')
plt.ylabel(r'$|B|$')
plt.title('Time evolution of $|B|$')
plt.legend()
if fig_save:
plt.savefig(fig_path + 'B_strength_time_evolution.png')
plt.show()
Study of Magnetic strength decay¶
In [ ]:
B_strength2 = B_strength.T
# plotting magnetic field strength at 7 different spatial points
cutoff_index = Nt // 2
plt.figure(figsize=(12, 6))
for i in np.linspace(0, Nr - 1, 9, dtype=int)[1:-1]:
plt.plot(np.linspace(0, T, Nt + 1)[cutoff_index:], B_strength2[i][cutoff_index:], label=f'r = {r[i]:.2f}')
plt.yscale('log')
plt.xlabel('t')
plt.ylabel(r'$|B|$')
plt.title('Time evolution of $|B|$ at different spatial points')
plt.grid()
plt.legend()
if fig_save:
plt.savefig(fig_path + 'B_strength_time_evolution_at_different_r.png')
plt.show()
In [ ]:
# plotting the magnetic field strength decay at r = 4
cutoff_index = 3 * Nt // 4
plt.figure(figsize=(12, 6))
plt.plot(np.linspace(0, T, Nt + 1)[cutoff_index:], B_strength2[Nr // 2][cutoff_index:], label = f'r = {r[Nr // 2]:.2f}')
plt.yscale('log')
plt.xlabel('t')
plt.ylabel(r'$|B|$')
# finding the decay rate
t = np.linspace(0, T, Nt + 1)[cutoff_index:]
B = B_strength2[Nr // 2][cutoff_index:]
slope, intercept = np.polyfit(t, np.log10(B), 1)
growth_rate = 10**(slope) - 1
# plotting the fit
plt.plot(t, 10**(slope * t + intercept),'--k', label = f'Fit, Growth rate = {growth_rate:.4f}')
plt.grid()
plt.title('Time evolution of $|B|$ at middle of r')
plt.legend()
if fig_save:
plt.savefig(fig_path + 'B_strength_time_evolution_at_middle_r.png')
plt.show()
Study of pitch angle variation over space¶
In [ ]:
pitchangles = np.degrees(np.arctan2(B_list[:, 1], B_list[:, 0]))
# plotting pitch angle time evolution
plt.figure(figsize=(12, 6))
for j in range(len(iterations_to_save_plot)):
plt.plot(r, pitchangles[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
plt.xlabel('r')
plt.ylabel('pitch angle (degrees)')
plt.title('Time evolution of pitch angle')
plt.legend()
if fig_save:
plt.savefig(fig_path + 'pitch_angle_time_evolution.png')
plt.show()
Study of different seed fields¶
In [ ]:
# study of different seed fields
# defining constants and parameters
eta = 1
r_i = 0.01 # start of spatial domain
r_f = 8.01 # end of spatial domain
T = 1.0 # Total simulation time
Nr = 100 # Number of spatial grid points
Nt = 4000 # Number of time steps
h = 1
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
iterations_to_save = 'all' # which time steps to save
iterations_to_save_plot = [int(i * Nt/8) for i in range(8 + 1)] # which time steps to plot
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
omega_0 = 50
r_omega = 4
h = 1
alpha = 0.5
alpha_cap = 1 #UPDATE
# defining the spatial dependent omega
omega = omega_0 / (np.sqrt(1 + (r / r_omega) ** 2))
# defining q
neg_q_omega = r * first_derivative_dYdx_FD(omega, dr, order = 6, ghost_zone_type='smooth')
R_alpha = alpha * h / eta
R_omega = (neg_q_omega * (h ** 2)) / eta
D = R_alpha * R_omega
# defining initial conditions
def B_phi01(r, phi, z):
return (3 * np.sin(((r - r_i)/(r_f - r_i))*np.pi)) + (1 * np.sin(((r - r_i)/(r_f - r_i))*2*np.pi)) + (2 * np.sin(((r - r_i)/(r_f - r_i))*5*np.pi))
def B_r01(r, phi, z):
return (3 * np.sin(((r - r_i)/(r_f - r_i))*4*np.pi)) + (1 * np.sin(((r - r_i)/(r_f - r_i))*2*np.pi)) + (2 * np.sin(((r - r_i)/(r_f - r_i))*3*np.pi))
def B_phi02(r, phi, z):
return (3 * np.sin(((r - r_i)/(r_f - r_i))*3*np.pi)) + (1 * np.sin(((r - r_i)/(r_f - r_i))*3*np.pi)) + (2 * np.sin(((r - r_i)/(r_f - r_i))*5*np.pi))
def B_r02(r, phi, z):
return (3 * np.sin(((r - r_i)/(r_f - r_i))*1*np.pi)) + (1 * np.sin(((r - r_i)/(r_f - r_i))*6*np.pi)) + (2 * np.sin(((r - r_i)/(r_f - r_i))*3*np.pi))
# definiing the time evolution of the magnetic field Br
def dBrdt(Br, Bphi, dBrdr, ddBrdr, t):
return (- (2/(np.pi*h)) * R_alpha * alpha_cap* Bphi) + (eta * (ddBrdr + (dBrdr / r) - (Br/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Br)))
# definiing the time evolution of the magnetic field Bphi
def dBphidt(Br, Bphi, dBphidr, ddBphidr, t):
return (R_omega * Br) + (eta * (ddBphidr + (dBphidr / r) - (Bphi/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Bphi)))
In [ ]:
# running 2 simulations with different seed fields
B_list1 = simulate_B_field(r_i, r_f, Nr, T, Nt, B_r01, B_phi01, dBrdt, dBphidt, order = order, ghost_zone_type = ghost_zone_type, iterations_to_save = iterations_to_save)
B_list2 = simulate_B_field(r_i, r_f, Nr, T, Nt, B_r02, B_phi02, dBrdt, dBphidt, order = order, ghost_zone_type = ghost_zone_type, iterations_to_save = iterations_to_save)
100%|██████████| 4000/4000 [00:14<00:00, 271.64it/s] 100%|██████████| 4000/4000 [00:14<00:00, 272.70it/s]
In [ ]:
# plotting Br and Bphi time evolition in two plots side by side
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
axs[0].set_title('Evolution of $B_r$')
axs[1].set_title('Evolution of $B_{\phi}')
axs[0].set_xlabel('r')
axs[1].set_xlabel('r')
axs[0].set_ylabel(r'$B_r$')
axs[1].set_ylabel(r'$B_{\phi}$')
for i, ax in enumerate(axs):
for j in range(len(iterations_to_save_plot)):
ax.plot(r, B_list1[iterations_to_save_plot[j], i], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
ax.legend()
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = {ghost_zone_type}\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\alpha = {alpha}, \\hat\\alpha = {alpha_cap}$\n'
figure_text3 = 'seed field 1: $B_r = 3 \\sin(\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(4\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
if fig_save:
plt.savefig(fig_path + 'Br_Bphi_time_evolution_seed_field1.png')
plt.show()
# plotting Br and Bphi time evolition in two plots side by side in seed field 2
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
axs[0].set_title('Evolution of $B_r$')
axs[1].set_title('Evolution of $B_{\phi}')
axs[0].set_xlabel('r')
axs[1].set_xlabel('r')
axs[0].set_ylabel(r'$B_r$')
axs[1].set_ylabel(r'$B_{\phi}$')
for i, ax in enumerate(axs):
for j in range(len(iterations_to_save_plot)):
ax.plot(r, B_list2[iterations_to_save_plot[j], i], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
ax.legend()
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = {ghost_zone_type}\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\alpha = {alpha}, \\hat\\alpha = {alpha_cap}$\n'
figure_text3 = 'seed field 2: $B_r = 3 \\sin(3\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(3\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(1 \\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(6\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
if fig_save:
plt.savefig(fig_path + 'Br_Bphi_time_evolution_seed_field2.png')
plt.show()
In [ ]:
# evolution of their magnetic field strength
B_strength1 = np.sqrt(B_list1[:, 0] ** 2 + B_list1[:, 1] ** 2)
B_strength2 = np.sqrt(B_list2[:, 0] ** 2 + B_list2[:, 1] ** 2)
plt.figure(figsize=(12, 6))
for j in range(len(iterations_to_save_plot)):
plt.plot(r, B_strength1[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
plt.xlabel('r')
plt.ylabel(r'$|B|$')
plt.title('Time evolution of $|B|$ with seed field 1')
plt.legend()
if fig_save:
plt.savefig(fig_path + 'B_strength_time_evolution_seed_field1.png')
plt.show()
plt.figure(figsize=(12, 6))
for j in range(len(iterations_to_save_plot)):
plt.plot(r, B_strength2[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
plt.xlabel('r')
plt.ylabel(r'$|B|$')
plt.title('Time evolution of $|B|$ with seed field 2')
plt.legend()
if fig_save:
plt.savefig(fig_path + 'B_strength_time_evolution_seed_field2.png')
plt.show()
Study of different boundary conditions¶
In [ ]:
# simulating the time evolution of the magnetic field with different ghost zone types
B_list1 = simulate_B_field(r_i, r_f, Nr, T, Nt, B_r01, B_phi01, dBrdt, dBphidt, order = order, ghost_zone_type = 'anti-symmetric', iterations_to_save = iterations_to_save)
B_list2 = simulate_B_field(r_i, r_f, Nr, T, Nt, B_r01, B_phi01, dBrdt, dBphidt, order = order, ghost_zone_type = 'symmetric', iterations_to_save = iterations_to_save)
B_list3 = simulate_B_field(r_i, r_f, Nr, T, Nt, B_r01, B_phi01, dBrdt, dBphidt, order = order, ghost_zone_type = 'relative anti-symmetric', iterations_to_save = iterations_to_save)
100%|██████████| 4000/4000 [00:14<00:00, 271.67it/s] 100%|██████████| 4000/4000 [00:14<00:00, 268.68it/s] 100%|██████████| 4000/4000 [00:14<00:00, 268.38it/s]
In [ ]:
# plotting Br and Bphi time evolition in (3,2) grid
fig, axs = plt.subplots(1, 3, figsize=(18, 6))
axs[0].set_title(r'Evolution of $B_r$')
axs[1].set_title(r'Evolution of $B_{\phi}$')
axs[2].set_title(r'Evolution of $|B|$')
axs[0].set_ylabel(r'$B_r$')
axs[1].set_ylabel(r'$B_{\phi}$')
axs[2].set_ylabel(r'$|B|$')
for i, ax in enumerate(axs):
ax.set_xlabel('r')
for j in range(len(iterations_to_save_plot)):
if i < 2:
ax.plot(r, B_list1[iterations_to_save_plot[j], i], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
else:
B_strength = np.sqrt(B_list1[:, 0] ** 2 + B_list1[:, 1] ** 2)
ax.plot(r, B_strength[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
ax.legend()
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = anti-symmetric\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\alpha = {alpha}, \\hat\\alpha = {alpha_cap}$\n'
figure_text3 = 'seed field 1: $B_r = 3 \\sin(\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(4\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
if fig_save:
plt.savefig(fig_path + 'Br_Bphi_time_evolution_ghost_zone_anti-symmetric.png')
plt.show()
fig, axs = plt.subplots(1, 3, figsize=(18, 6))
axs[0].set_title(r'Evolution of $B_r$')
axs[1].set_title(r'Evolution of $B_{\phi}$')
axs[2].set_title(r'Evolution of $|B|$')
axs[0].set_ylabel(r'$B_r$')
axs[1].set_ylabel(r'$B_{\phi}$')
axs[2].set_ylabel(r'$|B|$')
for i, ax in enumerate(axs):
ax.set_xlabel('r')
for j in range(len(iterations_to_save_plot)):
if i < 2:
ax.plot(r, B_list2[iterations_to_save_plot[j], i], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
else:
B_strength = np.sqrt(B_list2[:, 0] ** 2 + B_list2[:, 1] ** 2)
ax.plot(r, B_strength[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
ax.legend()
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = symmetric\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\alpha = {alpha}, \\hat\\alpha = {alpha_cap}$\n'
figure_text3 = 'seed field 1: $B_r = 3 \\sin(\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(4\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
if fig_save:
plt.savefig(fig_path + 'Br_Bphi_time_evolution_ghost_zone_symmetric.png')
plt.show()
fig, axs = plt.subplots(1, 3, figsize=(18, 6))
axs[0].set_title(r'Evolution of $B_r$')
axs[1].set_title(r'Evolution of $B_{\phi}$')
axs[2].set_title(r'Evolution of $|B|$')
axs[0].set_ylabel(r'$B_r$')
axs[1].set_ylabel(r'$B_{\phi}$')
axs[2].set_ylabel(r'$|B|$')
for i, ax in enumerate(axs):
ax.set_xlabel('r')
for j in range(len(iterations_to_save_plot)):
if i < 2:
ax.plot(r, B_list3[iterations_to_save_plot[j], i], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
else:
B_strength = np.sqrt(B_list3[:, 0] ** 2 + B_list3[:, 1] ** 2)
ax.plot(r, B_strength[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
ax.legend()
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = relative anti-symmetric\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\alpha = {alpha}, \\hat\\alpha = {alpha_cap}$\n'
figure_text3 = 'seed field 1: $B_r = 3 \\sin(\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(4\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
if fig_save:
plt.savefig(fig_path + 'Br_Bphi_time_evolution_ghost_zone_relative_anti-symmetric.png')
plt.show()
In [ ]:
# magnetic field strength evolution with different ghost zone types
B_strength1 = np.sqrt(B_list1[:, 0] ** 2 + B_list1[:, 1] ** 2)
B_strength2 = np.sqrt(B_list2[:, 0] ** 2 + B_list2[:, 1] ** 2)
B_strength3 = np.sqrt(B_list3[:, 0] ** 2 + B_list3[:, 1] ** 2)
plt.figure(figsize=(12, 6))
for j in range(len(iterations_to_save_plot)):
plt.plot(r, B_strength1[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
plt.xlabel('r')
plt.ylabel(r'$|B|$')
plt.title('Time evolution of $|B|$ with anti-symmetric ghost zone')
plt.legend()
if fig_save:
plt.savefig(fig_path + 'B_strength_time_evolution_ghost_zone_anti-symmetric.png')
plt.show()
plt.figure(figsize=(12, 6))
for j in range(len(iterations_to_save_plot)):
plt.plot(r, B_strength2[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
plt.xlabel('r')
plt.ylabel(r'$|B|$')
plt.title('Time evolution of $|B|$ with symmetric ghost zone')
plt.legend()
if fig_save:
plt.savefig(fig_path + 'B_strength_time_evolution_ghost_zone_symmetric.png')
plt.show()
plt.figure(figsize=(12, 6))
for j in range(len(iterations_to_save_plot)):
plt.plot(r, B_strength3[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
plt.xlabel('r')
plt.ylabel(r'$|B|$')
plt.title('Time evolution of $|B|$ with relative anti-symmetric ghost zone')
plt.legend()
if fig_save:
plt.savefig(fig_path + 'B_strength_time_evolution_ghost_zone_relative_anti-symmetric.png')
plt.show()
Study of Magnetic field strength evolution close to critical Dynamo number¶
In [ ]:
# defining constants and parameters
eta = 1
r_i = 0.01 # start of spatial domain
r_f = 10.01 # end of spatial domain
T = 10.0 # Total simulation time
Nr = 200 # Number of spatial grid points
Nt = int(T * 6000) # Number of time steps
h = 1
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
iterations_to_save = 'all' # which time steps to save
iterations_to_save_plot = [int(i * Nt/8) for i in range(8 + 1)] # which time steps to plot
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
trial = 12
omega_0 = 50
r_omega = 5
h = 1
alpha_cap = 1 #UPDATE
# defining the spatial dependent omega
omega = omega_0 / (np.sqrt(1 + (r / r_omega) ** 2))
# defining q
neg_q_omega = r * first_derivative_dYdx_FD(omega, dr, order = 6, ghost_zone_type='smooth')
R_alpha = 1.71
R_omega = (neg_q_omega * (h ** 2)) / eta
D = R_alpha * R_omega
# defining initial conditions
def B_phi0(r, phi, z):
return (3 * np.sin(((r - r_i)/(r_f - r_i))*np.pi)) + (1 * np.sin(((r - r_i)/(r_f - r_i))*2*np.pi)) + (2 * np.sin(((r - r_i)/(r_f - r_i))*5*np.pi))
def B_r0(r, phi, z):
return (3 * np.sin(((r - r_i)/(r_f - r_i))*4*np.pi)) + (1 * np.sin(((r - r_i)/(r_f - r_i))*2*np.pi)) + (2 * np.sin(((r - r_i)/(r_f - r_i))*3*np.pi))
# definiing the time evolution of the magnetic field Br
def dBrdt(Br, Bphi, dBrdr, ddBrdr, t):
return (- (2/(np.pi*h)) * R_alpha * alpha_cap* Bphi) + (eta * (ddBrdr + (dBrdr / r) - (Br/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Br)))
# definiing the time evolution of the magnetic field Bphi
def dBphidt(Br, Bphi, dBphidr, ddBphidr, t):
return (R_omega * Br) + (eta * (ddBphidr + (dBphidr / r) - (Bphi/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Bphi)))
In [ ]:
# plotting the omega values and dynamo number values across the spatial domain side by side
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
axs[0].plot(r, R_omega)
axs[0].set_xlabel('r')
axs[0].set_ylabel(r'$R_{\omega}$')
axs[0].set_title(r'$R_{\omega}$ vs r')
axs[1].plot(r, D)
axs[1].set_xlabel('r')
axs[1].set_ylabel('Dynamo number')
axs[1].set_title('Dynamo number vs r')
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = {ghost_zone_type}\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\hat\\alpha = {alpha_cap}, R_\\alpha = {R_alpha}$\n'
figure_text3 = 'seed fields: $B_r = 3 \\sin(\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(1.5\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
if fig_save:
plt.savefig(fig_path + f'omega_Dynamo_number_vs_r_trial{trial}.png')
plt.show()
In [ ]:
# running the simulations
B_list = simulate_B_field(r_i, r_f, Nr, T, Nt, B_r0, B_phi0, dBrdt, dBphidt, order = order, ghost_zone_type = ghost_zone_type, iterations_to_save = iterations_to_save)
0%| | 0/60000 [00:00<?, ?it/s]100%|██████████| 60000/60000 [07:26<00:00, 134.29it/s]
In [ ]:
# plotting Br and Bphi time evolition in two plots side by side
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
for i, ax in enumerate(axs):
for j in range(len(iterations_to_save_plot)):
ax.plot(r, B_list[iterations_to_save_plot[j], i], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
ax.set_xlabel('r')
if i == 0:
name = r'$B_r$'
else:
name = r'$B_{\phi}$'
ax.set_ylabel(name)
ax.set_title(f'Time evolution of {name}')
ax.grid()
ax.legend()
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = {ghost_zone_type}\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\hat\\alpha = {alpha_cap}, R_\\alpha = {R_alpha}$\n'
figure_text3 = 'seed fields: $B_r = 3 \\sin(\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(1.5\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
if fig_save:
plt.savefig(fig_path + f'Br_Bphi_time_evolution_trial{trial}.png')
plt.show()
print(f"saved figure at {fig_path + f'Br_Bphi_time_evolution_trial{trial}.png'}")
saved figure at figures/task2/Br_Bphi_time_evolution_trial12.png
In [ ]:
# plotting magnetic field strength time evolution
B_strength = np.sqrt(B_list[:, 0] ** 2 + B_list[:, 1] ** 2)
plt.figure(figsize=(12, 6))
for j in range(len(iterations_to_save_plot)):
plt.plot(r, B_strength[iterations_to_save_plot[j]], label=f't = {iterations_to_save_plot[j] * T/Nt:.4f}')
plt.xlabel('r')
plt.ylabel(r'$|B|$')
plt.title('Time evolution of $|B|$')
plt.legend()
figure_text1 = f'$r_i = {r_i}$, $r_f = {r_f}$, $T = {T}$, $N_r = {Nr}$, $N_t = {Nt}$, order = {order}, ghost_zone_type = {ghost_zone_type}\n'
figure_text2 = f'$\\omega_0 = {omega_0}, r_\\omega = {r_omega}, h = {h}, \\hat\\alpha = {alpha_cap}, R_\\alpha = {R_alpha}$\n'
figure_text3 = 'seed fields: $B_r = 3 \\sin(\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(5\\pi \\frac{r - r_i}{r_f - r_i})$\n'
figure_text4 = '$B_{\\phi} = 3 \\sin(1.5\\pi \\frac{r - r_i}{r_f - r_i}) + \\sin(2\\pi \\frac{r - r_i}{r_f - r_i}) + 2 \\sin(3\\pi \\frac{r - r_i}{r_F - r_i})$'
plt.figtext(0.5, 0.1, figure_text1 + figure_text2 + figure_text3 + figure_text4, horizontalalignment='center', fontsize=12)
plt.gcf().subplots_adjust(bottom=0.37)
plt.grid()
if fig_save:
plt.savefig(fig_path + f'B_strength_time_evolution_trial{trial}.png')
plt.show()
In [ ]:
B_strength = B_strength.T
In [ ]:
cutoff_index1 = int(Nt * 0.05)
# cutoff_index2 = int(Nt)
# cutoff_index1 = int(Nt * 0.9)
cutoff_index2 = int(1 * Nt // 10)
slopes = []
for i in range(Nr):
slope, intercept = np.polyfit(np.linspace(0,T,Nt+1)[cutoff_index1:cutoff_index2], np.log(B_strength[i][cutoff_index1:cutoff_index2]), 1)
slopes.append(slope)
slopes = np.array(slopes)
growth_rates = 10**(slopes) - 1
decay_rates = 1 - 10**(slopes)
gamma = slopes
# points where gamma crosses 0
r_zeros = r[np.where(np.diff(np.sign(gamma)))[0]]
print(f'Points where gamma crosses 0: {r_zeros}')
print(f'Dynamo numbers at these points: {D[np.where(np.diff(np.sign(gamma)))[0]]}')
# plotting gamma vs r
plt.figure(figsize=(12, 6))
plt.plot(r, gamma)
plt.xlabel('r')
plt.ylabel(r'$\gamma$')
plt.title('Growth rate vs r')
plt.grid()
plt.plot(r_zeros, np.zeros(len(r_zeros)), 'ro', label = 'gamma = 0')
for i in range(len(r_zeros)):
plt.annotate(f'D = {D[np.where(np.diff(np.sign(gamma)))[0][i]]:.2f}', (r_zeros[i], 0), textcoords="offset points", xytext=(0,10), ha='center')
plt.annotate(f'r = {r_zeros[i]:.2f}', (r_zeros[i], 0), textcoords="offset points", xytext=(0,-20), ha='center')
plt.legend()
if fig_save:
plt.savefig(fig_path + f'gamma_vs_r_trial{trial}_withDc.png')
plt.show()
Points where gamma crosses 0: [1.66829146] Dynamo numbers at these points: [-8.12468735]
In [ ]:
# cutoff_index1 = 0
cutoff_index2 = int(Nt)
cutoff_index1 = int(Nt * 0.9)
# cutoff_index2 = int(1 * Nt // 3)
slopes = []
for i in range(Nr):
slope, intercept = np.polyfit(np.linspace(0,T,Nt+1)[cutoff_index1:cutoff_index2], np.log(B_strength[i][cutoff_index1:cutoff_index2]), 1)
slopes.append(slope)
slopes = np.array(slopes)
growth_rates = 10**(slopes) - 1
decay_rates = 1 - 10**(slopes)
gamma = slopes
### Theoretical values (sir's paper)
t_d = 1
Dc = - (np.pi ** 2) / 32
local_gamma_theoretical = (np.sqrt(2 / np.pi) / t_d) * (np.sqrt(-D) - np.sqrt(-Dc))
# plotting gamma vs r
plt.figure(figsize=(12, 6))
plt.plot(r, gamma, label = f'Numerical global growth rate at end of simulation = {np.mean(gamma):.4f}')
plt.plot(r, local_gamma_theoretical, label = f'Theoretical local growth rate')
plt.xlabel('r')
plt.ylabel(r'$\gamma$')
plt.title('Growth rate vs r')
plt.grid()
plt.legend()
if fig_save:
plt.savefig(fig_path + f'gamma_vs_r_trial{trial}_global.png')
plt.show()
Alternate rate definition trials¶
In [ ]:
decay_rates = []
# cutoff_index1 = 0
# cutoff_index2 = int(Nt)
cutoff_index1 = int(Nt * 0.2)
cutoff_index2 = int(1 * Nt // 2)
slopes = []
for i in range(Nr):
slope, intercept = np.polyfit(np.linspace(0,T,Nt+1)[cutoff_index1:cutoff_index2], np.log10(B_strength[i][cutoff_index1:cutoff_index2]), 1)
slopes.append(slope)
slopes = np.array(slopes)
growth_rates = 10**(slopes) - 1
decay_rates = 1 - 10**(slopes)
# plotting decay rates of the magnetic field strength with respect to r
plt.figure(figsize=(12, 6))
plt.plot(r, growth_rates)
plt.xlabel('r')
plt.ylabel('Growth rate')
plt.title(f'Growth rates of the magnetic field strength, t = [{cutoff_index1 * T/Nt:.4f}, {cutoff_index2 * T/Nt:.4f}]')
# finding the indices at which the growth rate crosses zero
zero_indices = np.where(np.diff(np.sign(growth_rates)))[0]
print(f'Growth rate crosses zero at r = {r[zero_indices]}')
print(f"Dynamo number at r = {r[zero_indices]} is {D[zero_indices]}")
# plotting the growth rate and the dynamo number at the points where the growth rate crosses zero
plt.plot(r[zero_indices], growth_rates[zero_indices], 'ro', label='Growth rate crosses zero')
plt.annotate(f'({r[zero_indices][0]:.2f}, {growth_rates[zero_indices][0]:.4f})', (r[zero_indices][0], growth_rates[zero_indices][0]), textcoords="offset points", xytext=(0,10), ha='center')
plt.annotate(f'({r[zero_indices][1]:.2f}, {growth_rates[zero_indices][1]:.4f})', (r[zero_indices][1], growth_rates[zero_indices][1]), textcoords="offset points", xytext=(0,10), ha='center')
# annotating the dynamo number at the points where the growth rate crosses zero
plt.annotate(f'D = {D[zero_indices][0]:.2f}', (r[zero_indices][0], growth_rates[zero_indices][0]), textcoords="offset points", xytext=(0,-10), ha='center')
plt.annotate(f'D = {D[zero_indices][1]:.2f}', (r[zero_indices][1], growth_rates[zero_indices][1]), textcoords="offset points", xytext=(0,-10), ha='center')
plt.gcf().subplots_adjust(bottom=0.15)
plt.grid()
if fig_save:
plt.savefig(fig_path + f'decay_rate_across_space_trial{trial}.png')
plt.show()
Growth rate crosses zero at r = [1.81904523 8.35170854] Dynamo number at r = [1.81904523 8.35170854] is [-6.99766896 -6.57962909]
In [ ]:
# plotting growth rates of the magnetic field strength with respect to r
plt.figure(figsize=(12, 6))
plt.plot(r, growth_rates)
plt.xlabel('r')
plt.ylabel('Growth rate')
plt.title(f'Growth rates of the magnetic field strength, t = [{cutoff_index1 * T/Nt:.4f}, {cutoff_index2 * T/Nt:.4f}]')
# putting a growth rate cutoff and finding the indices of each peak
growth_rate_cutoff = 6
peak_indices = np.where(growth_rates > growth_rate_cutoff)[0]
# separating the peaks
peak_indices = np.split(peak_indices, np.where(np.diff(peak_indices) != 1)[0] + 1)
# plotting the peaks and annotating the dynamo number at the peak's center
# for peak in peak_indices:
# peak_center = np.mean(r[peak])
# plt.plot(r[peak], growth_rates[peak], label=f'Peak at r = {peak_center:.2f}')
# plt.plot(r[np.argmin(np.abs(r - peak_center))] , growth_rates[np.argmin(np.abs(r - peak_center))], 'ro')
# plt.annotate(f'D = {D[np.argmin(np.abs(r - peak_center))]:.2f}', (peak_center, growth_rates[np.argmin(np.abs(r - peak_center))]), textcoords="offset points", xytext=(0,-10), ha='center')
Dc_values = []
for peak in peak_indices:
peak_center = np.argmax(growth_rates[peak])
plt.plot(r[peak], growth_rates[peak], label=f'Peak at r = {peak_center:.2f}')
plt.plot(r[peak][peak_center], growth_rates[peak][peak_center], 'ro')
plt.annotate(f'D = {D[peak][peak_center]:.2f}', (r[peak][peak_center], growth_rates[peak][peak_center]), textcoords="offset points", xytext=(0,-10), ha='center')
Dc_values.append(D[peak][peak_center])
# plt.gcf().subplots_adjust(bottom=0.15)
plt.legend()
plt.grid()
if fig_save:
plt.savefig(fig_path + f'growth_rate_across_space_trial{trial}_alternative.png')
plt.show()
In [ ]:
# plotting critical dynamo nymbers in the dynamo number vs r plot
plt.figure(figsize=(12, 6))
plt.plot(r, D)
plt.xlabel('r')
plt.ylabel('Dynamo number')
plt.title('Dynamo number vs r')
for Dc in Dc_values:
plt.plot(r[np.argmin(np.abs(D - Dc)), ], Dc, 'ro')
plt.annotate(f'D = {Dc:.2f}', (r[np.argmin(np.abs(D - Dc)), ], Dc), textcoords="offset points", xytext=(0,-10), ha='center')
plt.gcf().subplots_adjust(bottom=0.15)
plt.grid()
if fig_save:
plt.savefig(fig_path + f'Dynamo_number_vs_r_trial{trial})with_dc.png')
plt.show()
