TASK 3 implementation¶
In [ ]:
# importing the required modules
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
from matplotlib.animation import FuncAnimation
from astropy import units as u
import os
fig_path = 'figures/task3/'
data_path = 'data/task3/'
current_path = os.getcwd()
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
In [ ]:
def simulate_B_field_till_steady(r_i, r_f, Nr, Nt_per_unitT, B_r0, B_phi0, dBrdt, dBphidt, order = 6, ghost_zone_type = 'relative anti-symmetric', iterations_to_save = None, global_growth_rate_tolerance = 1e-2, check_tolerance_after = 10, check_tolerance_every = 1, growth_rate_calc_using = 0.1, patience = 3, max_iterations = 100):
'''
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.
Nt_per_unitT: int
The number of time points per unit time.
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.
global_growth_rate_tolerance: float
The tolerance for the global growth rate to consider the system as steady.
Default is 1e-2.
check_tolerance_after: int
The number of unit times to wait before starting to check the global growth rate.
Default is 10.
check_tolerance_every: int
The number of unit times to wait before checking the global growth rate.
Default is 1.
growth_rate_calc_using: float
The fraction of the last unit time simulated, used to calculate the growth rate.
Default is 0.1.
patience: int
The number of unit times to wait before stopping the simulation.
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 = 1 / Nt_per_unitT # 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_per_unitT]
elif iterations_to_save == 'all':
iterations_to_save = list(range(Nt_per_unitT + 1))
if 0 in iterations_to_save:
B_list.append(np.copy(B))
N_unitT_passed = 0
iterations_below_tolerance = 0
N_datapoints_to_check = int(Nt_per_unitT * growth_rate_calc_using)
flag_separate_growth_save = False
# checking if the magnetic fields needed for growth rate is saved automatically, else saving separately
try:
temp = (iterations_to_save[-N_datapoints_to_check] - 1) == (Nt_per_unitT - N_datapoints_to_check)
if temp and iterations_to_save[-1] == Nt_per_unitT:
flag_separate_growth_save = True
growth_B_list = []
# print('separate growth rate calculation is enabled')
except:
flag_separate_growth_save = True
growth_B_list = []
# print('separate growth rate calculation is enabled')
check_growth_flag = False
growth_interval = 0 # variable to keep track of checking the growth rate every 'check_tolerance_every' unit times
while True:
if (N_unitT_passed >= check_tolerance_after) :
growth_interval += 1
if growth_interval == check_tolerance_every:
check_growth_flag = True
growth_interval = 0
if flag_separate_growth_save:
growth_B_list = [] # reset the growth list for the next check and start saving the growth data
#--------- RK4 time-stepping block ---------#
for n in tqdm(range(1, Nt_per_unitT + 1), desc=f"{N_unitT_passed + 1}th "):
# 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))
if check_growth_flag and flag_separate_growth_save and n > (Nt_per_unitT - N_datapoints_to_check):
growth_B_list.append(np.copy(B))
#------------------------------------------#
try:
growth_B_list = np.array(growth_B_list)
except:
pass
N_unitT_passed += 1
if check_growth_flag:
gammas = []
if flag_separate_growth_save:
B_strength = np.sqrt(growth_B_list[:, 0] ** 2 + growth_B_list[:, 1] ** 2).T
else:
B_strength = np.sqrt(B_list[-N_datapoints_to_check:, 0] ** 2 + B_list[-N_datapoints_to_check:, 1] ** 2).T
curr_time = np.linspace(Nt_per_unitT * (N_unitT_passed - growth_rate_calc_using) + 1, (Nt_per_unitT * N_unitT_passed), N_datapoints_to_check) * (1 / Nt_per_unitT)
# print(f'Shape of B_strength is {B_strength.shape} and shape of curr_time is {curr_time.shape}')
for i in range(Nr):
slope, _ = np.polyfit(curr_time, np.log(B_strength[i]), 1)
gammas.append(slope)
gammas = np.array(gammas)
mean_gamma = np.mean(gammas)
gammas = np.abs(gammas - mean_gamma)
mean_diff = np.mean(gammas)
print(f"Current mean gamma variation: {mean_diff:.5f}")
if (gammas < global_growth_rate_tolerance).all():
if iterations_below_tolerance == 0:
print(f"System has reached steady state after {N_unitT_passed} unit times. Patience interval started..")
iterations_below_tolerance += 1
else:
if iterations_below_tolerance > 0:
print(f"System did not maintain steady state. Patience interval reset..")
iterations_below_tolerance = 0
if iterations_below_tolerance >= patience:
break
check_growth_flag = False
try:
growth_B_list = []
except: pass
if N_unitT_passed == max_iterations:
break
B_list = np.array(B_list)
return B_list, mean_gamma
In [ ]:
def Dynamo_number_steady_search_z_implementntation(omega_0_start: float, omega_0end: float, omega_0step: float, interpolation: bool, *simulation_args):
'''
This function calculates the dynamo number for different values of omega_0.
Parameters:
omega_0_start: float
The starting value of omega_0.
omega_0_end: float
The ending value of omega_0.
omega_0_step: float
The step size of omega_0.
interpolation: bool
If True, if there is a transition from decay to growth, the function will interpolate the critical value of omega_0.
If False, the function will return the value of omega_0 closest to the transition.
*simulation_args: tuple
The arguments to be passed to the simulate_B_field_till_steady function.
Returns:
omegas: 1D array
The values of omega_0.
gammas: 1D array
The values of the corresponding growth rates.
critical_omega(s): 1D array
The values of omega_0 where the system transitions from decay to growth. (if any)
'''
omegas = []
gammas = []
for omega_val in np.arange(omega_0_start, omega_0end, omega_0step):
print(f"Simulating for omega_0 = {i}")
# defining the spatial dependent omega
omega = omega_val / (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_omega = (neg_q_omega * (h ** 2)) / eta_t
D = R_alpha * R_omega
B_list, mean_gamma = simulate_B_field_till_steady(*simulation_args)
omegas.append(omega)
gammas.append(mean_gamma)
omegas = np.array(omegas)
gammas = np.array(gammas)
critical_omega = []
for i in range(1,len(gammas)):
if gammas[i] * gammas[i-1] < 0:
if interpolation:
critical_omega.append(np.interp(0, [gammas[i-1], gammas[i]], [omegas[i-1], omegas[i]]))
else:
if np.abs(gammas[i-1]) < np.abs(gammas[i]):
critical_omega.append(omegas[i-1])
else:
critical_omega.append(omegas[i])
In [ ]:
def B_list_growth_analyzer(times, B_str_list, window = 25, stride = 1):
'''
This function analyzes the growth of the magnetic field B to calculate growth rate across r over time, and return the data.
Parameters:
times: 1D array
The time values corresponding to the B_list.
B_str_list: 3D array
The solution of the magnetic field strength B in time.
shape: (len(times), Nr)
Returns:
all_times: 1D array
The time values corresponding to the growth rates.
all_gammas: 2D array
The growth rates across r over time.
shape: (len(times), Nr)
'''
all_gammas = []
all_times = []
N_data = int((len(times) - window) / stride) + 1
for i in range(N_data):
start = i * stride
end = start + window
B_str = B_str_list[start:end]
gammas = []
for j in range(B_str.shape[1]):
slope, _ = np.polyfit(times[start:end], np.log(B_str[:, j]), 1)
gammas.append(slope)
all_gammas.append(gammas)
all_times.append(times[end-1])
all_times = np.array(all_times)
all_gammas = np.array(all_gammas)
return all_times, all_gammas
In [ ]:
def animate_to_gif(X_list, Y_lists, T_list, file_name, frames_per_second = 25, X_label = 'X', Y_label = 'Y', title = 'Evolution at t =', log_scale = False, x_lim = None, y_lim = None):
'''
This function creates an animation of the given data and saves it as a gif.
Parameters:
X_list: 1D array
The values of the x-axis.
Y_lists: 2D array
The values of the y-axis.
shape: (len(T_list), len(X_list))
T_list: 1D array
The values of time.
file_name: str
The name of the file to save the gif.
frames_per_second: int
The number of frames per second in the gif.
Default is 25.
X_label: str
The label for the x-axis.
Default is 'X'.
Y_label: str
The label for the y-axis.
Default is 'Y'.
title: str
The title of the plot.
Default is 'Evolution at t ='.
log_scale: bool
If True, the y-axis will be in log scale.
x_lim: tuple
The limits of the x-axis.
Default is None.
y_lim: tuple
The limits of the y-axis.
Default is None.
'''
# Creating a plot
figure = plt.figure(figsize=(12, 6))
plt.xlabel(X_label)
plt.ylabel(Y_label)
plt.title(title)
plt.xlim(np.min(X_list), np.max(X_list))
plt.grid()
def update(frame):
plt.cla()
line1 = plt.plot(X_list, Y_lists[frame])
plt.title(f"{title} {T_list[frame]:.4f} unit time")
plt.xlabel(X_label)
plt.ylabel(Y_label)
if x_lim is not None:
plt.xlim(x_lim)
if y_lim is not None:
plt.ylim(y_lim)
if log_scale:
plt.yscale('log')
return line1
anim = FuncAnimation(figure, update, frames = len(Y_lists))
anim.save(fig_path + file_name, writer='pillow', fps = frames_per_second)
In [ ]:
def animate_to_gif_advanced(X_list, Y_listss, T_list, file_name, frames_per_second = 25, X_label = 'X', Y_label = 'Y', title = 'Evolution at t =', labels=None, log_scale = False, x_lim = None, y_lim = None, colors = None, linestyles = None):
'''
This function creates an animation of the given data and saves it as a gif.
Parameters:
X_list: 1D array
The values of the x-axis.
Y_listss: 3D array
The values of the y-axis.
shape: (N_lines, len(T_list), len(X_list))
where N_lines is the number of lines to plot within the same plot.
T_list: 1D array
The values of time.
file_name: str
The name of the file to save the gif.
frames_per_second: int
The number of frames per second in the gif.
Default is 25.
X_label: str
The label for the x-axis.
Default is 'X'.
Y_label: str
The label for the y-axis.
Default is 'Y'.
title: str
The title of the plot.
Default is 'Evolution at t ='.
labels: list
The labels for the different lines.
Default is None.
if labels is 2D array, it should be of shape (N_lines, len(T_list))
log_scale: bool
If True, the y-axis will be in log scale.
Default is False.
x_lim: tuple
The limits of the x-axis.
Default is None.
y_lim: tuple
The limits of the y-axis.
Default is None.
colors: list
The colors of the lines.
Default is None.
shape: (N_lines,)
linestyles: list
The linestyles of the lines.
Default is None.
shape: (N_lines,)
'''
N_lines = len(Y_listss)
if labels is None:
labels = [str(i) for i in range(N_lines)]
labels = np.array(labels)
labels_2D = False
if len(labels.shape) == 2:
labels_2D = True
if colors is None:
color_list = ['b', 'g', 'r', 'c', 'm', 'y', 'k']
colors = []
for i in range(N_lines):
colors.append(color_list[i % len(color_list)])
if linestyles is None:
linestyles = ['-' for i in range(N_lines)]
# Creating a plot
figure = plt.figure(figsize=(12, 6))
plt.xlabel(X_label)
plt.ylabel(Y_label)
plt.title(title)
plt.xlim(np.min(X_list), np.max(X_list))
plt.grid()
def update(frame):
plt.cla()
if labels_2D:
for i in range(N_lines):
line1 = plt.plot(X_list, Y_listss[i][frame], label=labels[i][frame], c=colors[i], ls=linestyles[i])
else:
for i in range(N_lines):
line1 = plt.plot(X_list, Y_listss[i][frame], label=labels[i], c=colors[i], ls=linestyles[i])
plt.xlabel(X_label)
plt.ylabel(Y_label)
plt.legend()
if x_lim is not None:
plt.xlim(x_lim)
if y_lim is not None:
plt.ylim(y_lim)
if log_scale:
plt.yscale('log')
plt.title(f"{title} {T_list[frame]:.4f} unit time")
return line1
anim = FuncAnimation(figure, update, frames=len(Y_listss[0]))
anim.save(fig_path + file_name, writer='pillow', fps = frames_per_second)
In [ ]:
def animate_to_gif_slider(X_list, Y_lists, window, stride, file_name, frames_per_second = 25, X_label = 'X', Y_label = 'Y', title = 'Evolution at t =', log_scale = False, x_lim = None, y_lim = None):
'''
This function creates an animation of the given data in a sliding window and saves it as a gif.
Parameters:
X_list: 1D array
The values of the x-axis.
Y_lists: 1D array
The values of the y-axis.
shape: (X_list)
window: int
The size of the window.
stride: int
The stride of the window.
file_name: str
The name of the file to save the gif.
frames_per_second: int
The number of frames per second in the gif.
Default is 25.
X_label: str
The label for the x-axis.
Default is 'X'.
Y_label: str
The label for the y-axis.
Default is 'Y'.
title: str
The title of the plot.
Default is 'Evolution at t ='.
x_lim: tuple
The limits of the x-axis.
Default is None.
y_lim: tuple
The limits of the y-axis.
Default is None.
'''
# Creating a plot
# Creating a plot
figure = plt.figure(figsize=(12, 6))
plt.xlabel(X_label)
plt.ylabel(Y_label)
plt.title(title)
plt.grid()
N_frames = int((len(Y_lists) - window) / stride) + 1
def update(frame):
start = frame * stride
end = start + window
plt.cla()
line1 = plt.plot(X_list[start:end], Y_lists[start:end])
plt.title(f"{title} {X_list[start]:.4f} to {X_list[end-1]:.4f}")
plt.xlabel(X_label)
plt.ylabel(Y_label)
if x_lim is not None:
plt.xlim(x_lim)
if y_lim is not None:
plt.ylim(y_lim)
if log_scale:
plt.yscale('log')
return line1
anim = FuncAnimation(figure, update, frames = N_frames)
anim.save(fig_path + file_name, writer='pillow', fps = frames_per_second)
In [ ]:
def animate_to_gif_slider_advanced(X_list, Y_listss, window, stride, file_name, frames_per_second = 25, X_label = 'X', Y_label = 'Y', title = 'Evolution at t =', log_scale = False, labels=None, x_lim = None, y_lim = None):
'''
This function creates an animation of the given data in a sliding window and saves it as a gif.
Parameters:
X_list: 1D array
The values of the x-axis.
Y_listss: 3D array
The values of the y-axis.
shape: (N_lines, len(X_list))
where N_lines is the number of lines to plot within the same plot.
file_name: str
The name of the file to save the gif.
frames_per_second: int
The number of frames per second in the gif.
Default is 25.
X_label: str
The label for the x-axis.
Default is 'X'.
Y_label: str
The label for the y-axis.
Default is 'Y'.
title: str
The title of the plot.
Default is 'Evolution at t ='.
log_scale: bool
If True, the y-axis will be in log scale.
labels: list
The labels for the different lines.
Default is None.
if 2D array, the labels for each line at each frame.
shape: (N_lines, len(Y_listss))
x_lim: tuple
The limits of the x-axis.
Default is None.
y_lim: tuple
The limits of the y-axis.
Default is None.
'''
N_frames = int((Y_listss.shape[1] - window) / stride) + 1
N_lines = Y_listss.shape[0]
if labels is None:
labels = [str(i) for i in range(N_lines)]
label_2D_flag = False
if len(labels.shape) == 2:
label_2D_flag = True
# Creating a plot
figure = plt.figure(figsize=(12, 6))
plt.xlabel(X_label)
plt.ylabel(Y_label)
plt.title(title)
plt.grid()
def update(frame):
start = frame * stride
end = start + window
plt.cla()
if label_2D_flag:
for i in range(N_lines):
line1 = plt.plot(X_list[start:end], Y_listss[i][start:end], label=labels[i][frame])
else:
for i in range(N_lines):
line1 = plt.plot(X_list[start:end], Y_listss[i][start:end], label=labels[i])
plt.xlabel(X_label)
plt.ylabel(Y_label)
plt.legend()
if x_lim is not None:
plt.xlim(x_lim)
if y_lim is not None:
plt.ylim(y_lim)
if log_scale:
plt.yscale('log')
plt.title(f"{title} {X_list[start]:.4f} to {X_list[end-1]:.4f}")
return line1
anim = FuncAnimation(figure, update, frames=N_frames)
anim.save(fig_path + file_name, writer='pillow', fps = frames_per_second)
In [ ]:
# equations of Magnetic field evolution
# definiing the time evolution of the magnetic field Br
def dBrdt(Br, Bphi, dBrdr, ddBrdr, t):
return (- (2/(np.pi*h)) * R_alpha * Bphi) + (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) + (ddBphidr + (dBphidr / r) - (Bphi/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Bphi))
# definiing the time evolution of the magnetic field Br
def dBrdt_withVz(Br, Bphi, dBrdr, ddBrdr, t):
return (- (2/(np.pi*h)) * V_z * Br) + (- (2/(np.pi*h)) * R_alpha * Bphi) + (ddBrdr + (dBrdr / r) - (Br/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Br))
# definiing the time evolution of the magnetic field Bphi
def dBphidt_withVz(Br, Bphi, dBphidr, ddBphidr, t):
return (- (2/(np.pi*h)) * V_z * Bphi) + (R_omega * Br) + (ddBphidr + (dBphidr / r) - (Bphi/(r ** 2)) - (((np.pi ** 2) / (4 * (h**2))) * Bphi))
N.B. Even if Pylance may complain about variables not defined in the above section, After setting the grid and initializing the implementation, the variables will be defined and the code will run without any issues.
Our parameters for system¶
In [ ]:
# Defined parameters 3 - eta_t = 1?
r_i_ = 0.01 * u.kpc # kpc
r_f_ = 15.01 * u.kpc # kpc
l_kpc_ = 0.1 * u.kpc # kpc
u0_ = 10 * u.km / u.s # km/s
omega_0_ = 220 * (u.km / u.s) / u.kpc # km/s/kpc
r_omega_ = 2 * u.kpc # kpc
r_h_ = 10 * u.kpc # kpc
alpha_ = 1.5 * u.km / (u.s) # km/s
B_0_ = 1e-5 * u.G # G
h_ = 0.4 * u.kpc # kpc
# t_d_ = 0.733 * u.Gyr # Gyr
# eta_t_ = (h_ ** 2) / t_d_ # kpc^2 / Gyr
eta_t_ = l_kpc_ * u0_ / 3 # kpc * km / s
t_d_ = ((h_ ** 2) / eta_t_).to(u.Gyr) # Gyr
V_z_ = 1 * (u.km / u.s) # km/s
In [ ]:
# dimensionless parameters scaling
length_scale = h_
time_scale = t_d_
r_i = (r_i_ / length_scale).value
r_f = (r_f_ / length_scale).value
l_kpc = (l_kpc_ / length_scale).value
u0 = (u0_.to(u.kpc / u.Gyr) * time_scale / length_scale).value
omega_0 = (omega_0_.to(1 / u.Gyr) * time_scale).value
r_omega = (r_omega_ / length_scale).value
r_h = (r_h_ / length_scale).value
B_0 = (B_0_ / B_0_).value
t_d = (t_d_ / time_scale).value
h = (h_ / length_scale).value
eta_t = (eta_t_.to(u.kpc ** 2 / u.Gyr) * time_scale / (length_scale ** 2)).value
alpha = (alpha_.to(u.kpc / u.Gyr) * (time_scale / length_scale)).value
V_z = (V_z_.to(u.kpc / u.Gyr) * (time_scale / length_scale)).value
In [ ]:
# Printing the practical and dimensionless parameters in a table
params_table = PrettyTable()
params_table.field_names = ['Parameter', 'Practical Value', 'Dimensionless Value']
params_table.add_row(['r_i', f'{r_i_:.4f}', f'{r_i:.4f}'])
params_table.add_row(['r_f', f'{r_f_:.4f}', f'{r_f:.4f}'])
params_table.add_row(['l_kpc', f'{l_kpc_:.4f}', f'{l_kpc:.4f}'])
params_table.add_row(['u0', f'{u0_:.4f}', f'{u0:.4f}'])
params_table.add_row(['omega_0', f'{omega_0_:.4f}', f'{omega_0:.4f}'])
params_table.add_row(['r_omega', f'{r_omega_:.4f}', f'{r_omega:.4f}'])
params_table.add_row(['r_h', f'{r_h_:.4f}', f'{r_h:.4f}'])
params_table.add_row(['alpha', f'{alpha_:.4f}', f'{alpha:.4f}'])
params_table.add_row(['B_0', f'{B_0_:.4f}', f'{B_0:.4f}'])
params_table.add_row(['h', f'{h_:.4f}', f'{h:.4f}'])
params_table.add_row(['t_d', f'{t_d_:.4f}', f'{t_d:.4f}'])
params_table.add_row(['eta_t', f'{eta_t_:.4f}', f'{eta_t:.4f}'])
params_table.add_row(['V_z', f'{V_z_:.4f}', f'{V_z:.4f}'])
print(params_table)
+-----------+-----------------------+---------------------+ | Parameter | Practical Value | Dimensionless Value | +-----------+-----------------------+---------------------+ | r_i | 0.0100 kpc | 0.0250 | | r_f | 15.0100 kpc | 37.5250 | | l_kpc | 0.1000 kpc | 0.2500 | | u0 | 10.0000 km / s | 12.0000 | | omega_0 | 220.0000 km / (kpc s) | 105.6000 | | r_omega | 2.0000 kpc | 5.0000 | | r_h | 10.0000 kpc | 25.0000 | | alpha | 1.5000 km / s | 1.8000 | | B_0 | 0.0000 G | 1.0000 | | h | 0.4000 kpc | 1.0000 | | t_d | 0.4693 Gyr | 1.0000 | | eta_t | 0.3333 km kpc / s | 1.0000 | | V_z | 1.0000 km / s | 1.2000 | +-----------+-----------------------+---------------------+
Checking the R_omega and D profiles¶
In [ ]:
r = np.linspace(r_i, r_f, 100)
dr = ((r_f - r_i) / 99)
# defining the spatial dependent omega
omega = omega_0 / (np.sqrt(1 + (r / r_omega) ** 2))
# defining spatial dependent h
h = h * np.sqrt(1 + (r / r_h) ** 2)
# defining the spatial dependent alpha
alpha = omega * (l_kpc ** 2) / (h)
# defining q
neg_q_omega = r * first_derivative_dYdx_FD(omega, dr, order = 6, ghost_zone_type='smooth')
R_alpha = alpha
R_omega = neg_q_omega # * (h ** 2)) / eta_t # sjnce all the terms are already dimensionless
D = R_alpha * R_omega
# plot alpha and omega
plt.figure(figsize=(12, 6))
plt.plot(r, alpha, label='alpha')
plt.plot(r, omega, label='omega')
plt.xlabel('r')
plt.ylabel('Value')
plt.title('alpha and omega')
plt.legend()
plt.grid()
plt.show()
# plot R_omega and D
plt.figure(figsize=(12, 6))
plt.plot(r*length_scale, R_omega)
plt.xlabel('r [kpc]')
plt.ylabel(r'$R_\omega$')
plt.title(r'Variation of $R_\omega$ with $r$')
plt.savefig(fig_path + 'R_omega_wth_r.png')
plt.show()
plt.figure(figsize=(12, 6))
plt.plot(r*length_scale, D)
plt.xlabel('r [kpc]')
plt.ylabel(r'$D$')
plt.title(r'Variation of Dynamo number with $r$')
plt.savefig(fig_path + 'Dynamo_number_wth_r.png')
plt.show()
# plotting R_alpha
plt.figure(figsize=(12, 6))
plt.plot(r*length_scale, R_alpha)
plt.xlabel('r [kpc]')
plt.ylabel(r'$R_\alpha$')
plt.title(r'Variation of $R_\alpha$ with $r$')
plt.savefig(fig_path + 'R_alpha_wth_r.png')
plt.show()
Implementation¶
In [ ]:
# defining simulation parameters
Nt_per_unitT = 5000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial32_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
dt = t_d / Nt_per_unitT
rand_seed = 0 # seed for Initial magnetic field
# checking the courant condition
C_val = (2 * eta_t * dt) / ((dr) ** 2)
if C_val > 1:
print(f"Courant condition is not satisfied. C = {C_val:.4f} > 1")
print(f"Need eta_t < {eta_t / C_val:.4f} (dimensionless)")
print(f"or Nt_per_unitT > {int(Nt_per_unitT * C_val)}")
print(f'or need Nr > {int(np.sqrt(((Nr - 1)**2) / C_val) + 1)}')
else:
print(f"Courant condition is satisfied. C = {C_val:.4e} < 1")
if C_val > 0.01:
print(f"Warning: C is close to 1. Consider increasing Nt_per_unitT or decreasing Nr")
# defining the outflow
V_z_ = 4 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
r_Vz_= 10 * u.kpc # kpc
r_Vz = (r_Vz_ / h_).value
fraction = 0.5
V_z = V_z * (1 + fraction * ((((r - r_Vz) / r_Vz) ** 2) - 1))
# disk - flaring setting
h_ = 0.4 * u.kpc # kpc
h = (h_ / h_).value
# r_h_ = 10 * u.kpc # kpc
# r_h = (r_h_ / h_).value
# h = h * ((1 + ((r / r_h) ** 2)) ** 0.5)
# defining the spatial dependent omega
omega = omega_0 / (np.sqrt(1 + (r / r_omega) ** 2))
# defining the spatial dependent alpha
alpha = omega * (l_kpc ** 2) / (h)
# defining q
neg_q_omega = r * first_derivative_dYdx_FD(omega, dr, order = 6, ghost_zone_type='smooth')
R_alpha = alpha
R_omega = neg_q_omega # * (h ** 2)) / eta_t # sjnce all the terms are already dimensionless
D = R_alpha * R_omega
# defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
np.random.seed(rand_seed)
return -1.5* 1e-4 * B_0 * np.random.normal(5, 1, r.shape) * np.sin((r - r_i) * np.pi / (r_f - r_i))
def B_r0(r, phi, z):
np.random.seed(rand_seed)
return B_0 * 1e-4 * np.random.normal(5, 1, r.shape) * np.sin((r - r_i) * np.pi / (r_f - r_i))
# return np.zeros(r.shape)
Courant condition is satisfied. C = 2.7878e-03 < 1
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 + trial_name + 'Dynamo_number_vs_r.png')
plt.show()
# plotting omega vs r
plt.figure(figsize=(12, 6))
plt.plot(r*h_, omega/t_d_)
plt.xlabel('r [kpc]')
plt.ylabel(r'$\omega$ [1/Gyr]')
plt.title(r'$\omega$ vs r')
if fig_save:
plt.savefig(fig_path + trial_name + 'omega_vs_r.png')
plt.show()
# plotting R_omega vs r
plt.figure(figsize=(12, 6))
plt.plot(r*h_, R_omega)
plt.xlabel('r [kpc]')
plt.ylabel(r'$R_{\omega}$')
plt.title(r'$R_{\omega}$ vs r')
if fig_save:
plt.savefig(fig_path + trial_name + 'R_omega_vs_r.png')
plt.show()
try:
if len(h) > 1:
# plotting h with r
plt.figure(figsize=(12, 6))
plt.plot(r*h_, h*h_)
plt.xlabel('r [kpc]')
plt.ylabel('h [kpc]')
plt.title('h vs r')
if fig_save:
plt.savefig(fig_path + trial_name + 'h_vs_r.png')
plt.show()
except:
pass
try:
if len(V_z) > 1:
# plotting V_z with r
plt.figure(figsize=(12, 6))
plt.plot(r*h_, V_z)
plt.xlabel('r [kpc]')
plt.ylabel(r'$V_z$ [kpc/Gyr]')
plt.title(r'$V_z$ vs r')
if fig_save:
plt.savefig(fig_path + trial_name + 'V_z_vs_r.png')
plt.show()
except:
pass
In [ ]:
# plotting initial conditions
plt.figure(figsize=(12, 6))
plt.plot(r, B_r0(r, 0, 0), label = r'$B_r$')
plt.plot(r, B_phi0(r, 0, 0), label = r'$B_{\phi}$')
plt.xlabel('r')
plt.ylabel('B')
plt.title('Initial conditions')
plt.legend()
plt.grid()
if fig_save:
plt.savefig(fig_path + trial_name + 'Initial_conditions.png')
plt.show()
In [ ]:
# running the simulations
B_list, mean_gamma = simulate_B_field_till_steady(r_i, r_f, Nr, Nt_per_unitT, B_r0, B_phi0, dBrdt_withVz, dBphidt_withVz, order = order, ghost_zone_type = ghost_zone_type, iterations_to_save = iterations_to_save_plot, global_growth_rate_tolerance = global_growth_rate_tolerance, check_tolerance_after = check_tolerance_after, check_tolerance_every = check_tolerance_every, growth_rate_calc_using = growth_rate_calc_using, patience = patience, max_iterations = max_iterations)
1th : 0%| | 0/5000 [00:00<?, ?it/s]1th : 100%|██████████| 5000/5000 [00:18<00:00, 266.45it/s] 2th : 100%|██████████| 5000/5000 [00:18<00:00, 265.69it/s] 3th : 100%|██████████| 5000/5000 [00:18<00:00, 266.64it/s] 4th : 100%|██████████| 5000/5000 [00:19<00:00, 262.16it/s] 5th : 100%|██████████| 5000/5000 [00:18<00:00, 264.24it/s] 6th : 100%|██████████| 5000/5000 [00:19<00:00, 258.46it/s] 7th : 100%|██████████| 5000/5000 [00:18<00:00, 268.96it/s] 8th : 100%|██████████| 5000/5000 [00:18<00:00, 267.80it/s] 9th : 100%|██████████| 5000/5000 [00:18<00:00, 263.84it/s] 10th : 100%|██████████| 5000/5000 [00:18<00:00, 263.95it/s] 11th : 100%|██████████| 5000/5000 [00:19<00:00, 260.99it/s]
Current mean gamma variation: 0.06186
12th : 100%|██████████| 5000/5000 [00:19<00:00, 262.72it/s]
Current mean gamma variation: 0.01951
13th : 100%|██████████| 5000/5000 [00:18<00:00, 270.42it/s]
Current mean gamma variation: 0.00779
14th : 100%|██████████| 5000/5000 [00:18<00:00, 263.79it/s]
Current mean gamma variation: 0.00337
15th : 100%|██████████| 5000/5000 [00:19<00:00, 262.26it/s]
Current mean gamma variation: 0.00151 System has reached steady state after 15 unit times. Patience interval started..
16th : 100%|██████████| 5000/5000 [00:18<00:00, 264.12it/s]
Current mean gamma variation: 0.00069
17th : 100%|██████████| 5000/5000 [00:18<00:00, 264.53it/s]
Current mean gamma variation: 0.00032
In [ ]:
np.save(data_path + trial_name + 'B_list.npy', B_list)
print(f"Mean growth rate: {mean_gamma:.5f}")
Mean growth rate: 5.12978
Loading the data and processing¶
In [ ]:
trial_name_ = trial_name
# trial+name_ = 'trial1_'
B_list = np.load(data_path + trial_name_ + 'B_list.npy')
print(B_list.shape)
(426, 2, 100)
In [ ]:
# Calculating the corresponding time to the B_list
total_time = (((B_list.shape)[0] - 1) // frames_per_unitT)
times = np.linspace(0, total_time, B_list.shape[0])
print(times.shape)
(426,)
In [ ]:
# plotting final eigen modes
plt.figure(figsize=(12, 6))
plt.plot(r, B_list[-1, 0], label = r'$B_r$')
plt.plot(r, B_list[-1, 1], label = r'$B_{\phi}$')
plt.plot(r, np.sqrt(B_list[-1, 0] ** 2 + B_list[-1, 1] ** 2), label = r'$|B|$')
plt.xlabel('r')
plt.ylabel('B')
plt.title(f'Final eigen modes - {trial_name_}')
plt.legend()
plt.grid()
if fig_save:
plt.savefig(fig_path + trial_name_ + 'Final_eigen_modes.png')
plt.show()
In [ ]:
animate_to_gif_advanced(r, [B_list[:,0], B_list[:,1], np.sqrt((B_list[:,0] ** 2) + (B_list[:,1] ** 2))], times, f"{trial_name}B_all.gif", frames_per_second=frames_per_unitT, X_label=r"r [kpc]", Y_label=r"$B_r, B_{\phi},|B| [G]$", title=r"Magnetic field evolution", labels = [r"$B_r$",r"$B_{\phi}$",r"$|B|$"])
In [ ]:
# plotting the B_r, B_phi and |B| wr.t. t at certain r indices
r_indices = np.linspace(3, 53, 6) # 7 points in focus
r_indices = r_indices.astype(int)
r_vals_at_indices = r[r_indices]
B_strength_list = np.sqrt(B_list[:,0]**2 + B_list[:,0]**2)
In [ ]:
all_times, all_gammas = B_list_growth_analyzer(times, B_strength_list, window = frames_per_unitT, stride = 1)
print(all_times.shape)
print(all_gammas.shape)
# plotting the growth rate at different r values
animate_to_gif(r, all_gammas, all_times, trial_name + 'gamma_evolution.gif', frames_per_second = frames_per_unitT, X_label = 'r [dimensionless]', Y_label = r'$\gamma [dimensionless]$', title = r'$\gamma$ evolution at t =')
C:\Users\adhil\AppData\Local\Temp\ipykernel_15536\1483822731.py:32: RuntimeWarning: divide by zero encountered in log slope, _ = np.polyfit(times[start:end], np.log(B_str[:, j]), 1)
(402,) (402, 100)
In [ ]:
r_labels = [f'r = {r_vals_at_indices[i]:.4f}: ' for i in range(len(r_indices))]
r_labels2 = [[f'r = {r_vals_at_indices[i]:.4f}: slope = {all_gammas[j, i]:.4f}' for i in range(len(r_indices))] for j in range(all_gammas.shape[0])]
r_labels = np.array(r_labels)
r_labels2 = np.array(r_labels2)
print(times.shape)
print(B_strength_list[:, r_indices].shape)
# plotting the growth of the magnetic field at certain r indices in log scale for full time
plt.figure(figsize=(12, 6))
for i in range(len(r_indices)):
plt.plot(times, B_strength_list[:, r_indices[i]], label = r_labels[i])
plt.xlabel('t [dimensionless]')
plt.ylabel(r'$|B| [dimensionless]$')
plt.title(r'$|B|$ growth at certain r indices for full time')
plt.yscale('log')
plt.legend()
plt.grid()
if fig_save:
plt.savefig(fig_path + trial_name + 'B_growth_at_r_indices_all_time.png')
plt.show()
(426,) (426, 6)
In [ ]:
# plotting the growth of the magnetic field at certain r indices in log scale for full time sliding window
animate_to_gif_slider_advanced(times, B_strength_list[:, r_indices].T, window = frames_per_unitT, stride = 1, file_name = trial_name + 'B_strength_evolution_slider_at_diff_r.gif', frames_per_second = frames_per_unitT, X_label = 't [Gyr]', Y_label = '|B| [G]', title = 'Evolution of |B| at t =', log_scale = True, labels = r_labels2.T)
In [ ]:
pitch_angles = np.arctan(B_list[:,0]/B_list[:,1]) * 180 / np.pi
print(pitch_angles.shape)
print(type(pitch_angles[0,0]), pitch_angles[0,0])
# plotting the pitch angle evolution
# animate_to_gif(r, pitch_angles, times, f"{trial_name}pitch_angle_evolution.gif", frames_per_second=frames_per_unitT, X_label=r"r [kpc]", Y_label=r"$\theta$ [degrees]", title=r"Pitch angle evolution")
(426, 100) <class 'numpy.float64'> nan
C:\Users\adhil\AppData\Local\Temp\ipykernel_15536\443814250.py:1: RuntimeWarning: invalid value encountered in divide pitch_angles = np.arctan(B_list[:,0]/B_list[:,1]) * 180 / np.pi
In [ ]:
# plotting the pitch angle vs r at last time
plt.figure(figsize=(12, 6))
plt.plot(r, pitch_angles[-1])
plt.xlabel('r [dimensionless]')
plt.ylabel(r'$\theta$ [degrees]')
plt.title(r'Pitch angle vs r at last time')
if fig_save:
plt.savefig(fig_path + trial_name + 'pitch_angle_vs_r_at_last_time.png')
plt.show()
# plotting the pitch angle vs t for certain r indices
plt.figure(figsize=(12, 6))
for i in range(len(r_indices)):
plt.plot(times, pitch_angles[:, r_indices[i]], label = f'r = {r_vals_at_indices[i]:.4f}')
plt.xlabel('t [dimensionless]')
plt.ylabel(r'$\theta$ [degrees]')
plt.title(r'Pitch angle vs t')
plt.legend()
if fig_save:
plt.savefig(fig_path + trial_name + 'pitch_angle_vs_t_at_diff_r.png')
plt.show()
In [ ]:
# pitch angle at different r values
animate_to_gif_slider_advanced(times, pitch_angles[:, r_indices].T, window = frames_per_unitT, stride = 1, file_name = trial_name + 'pitch_angle_evolution_slider_at_diff_r.gif', frames_per_second = frames_per_unitT, X_label = 't [Gyr]', Y_label = r'$\theta$ [degrees]', title = 'Evolution of pitch angle at t =', log_scale = False, labels = r_labels.T)
Combined analysis (without disk flaring)¶
In [ ]:
# drawing growth rate evolution for differetn Vz
V_z_list = [0, 10, 1 ,2 , 4, 6, 6.7, 8.5, 8]
trial_names = [f'trial{i}_' for i in range(1, 10)]
label_V_z = [f'V_z = {V_z_list[i]:.4f}' for i in range(len(V_z_list))]
# picking only even V_z values
mask_even = [V_z_list[i] % 2 == 0 for i in range(len(V_z_list))]
V_z_list = np.array(V_z_list)
V_z_list = V_z_list[mask_even]
label_V_z = np.array(label_V_z)
label_V_z = label_V_z[mask_even]
trial_names = np.array(trial_names)
trial_names = trial_names[mask_even]
all_gammas_all_Vz = []
for i in tqdm(range(len(V_z_list))):
B_list = np.load(data_path + trial_names[i] + 'B_list.npy')
B_strength_list = np.sqrt(B_list[:,0]**2 + B_list[:,0]**2)
all_times, all_gammas = B_list_growth_analyzer(times, B_strength_list, window = frames_per_unitT, stride = 1)
all_gammas_all_Vz.append(np.copy(all_gammas))
all_gammas_all_Vz = np.array(all_gammas_all_Vz)
all_times = np.array(all_times)
print(all_gammas_all_Vz.shape)
print(all_times.shape)
print(f'animating the evolution of the growth rate for different V_z values')
animate_to_gif_advanced(r, all_gammas_all_Vz, all_times, f"gamma_evolution_all_Vz2.gif", frames_per_second=frames_per_unitT, X_label=r"r [kpc]", Y_label=r"$\gamma [dimensionless]$", title=r"$\gamma$ evolution at t =", labels = label_V_z)
100%|██████████| 6/6 [00:10<00:00, 1.76s/it]
(6, 377, 100) (377,) animating the evolution of the growth rate for different V_z values
Combined analysis (without disk flaring) 2 (new alpha variation)¶
In [ ]:
# drawing growth rate evolution for differetn Vz
V_z_list = [0, 2, 4, 6, 8, 10]
trial_names = [f'trial{i}_' for i in range(19, 25)]
label_V_z = [f'V_z = {V_z_list[i]:.4f}' for i in range(len(V_z_list))]
# picking only even V_z values
mask = [V_z_list[i] % 2 == 0 for i in range(len(V_z_list))]
mask = [i % 1 == 0 for i in range(len(V_z_list))]
V_z_list = np.array(V_z_list)
V_z_list = V_z_list[mask]
label_V_z = np.array(label_V_z)
label_V_z = label_V_z[mask]
trial_names = np.array(trial_names)
trial_names = trial_names[mask]
all_gammas_all_Vz = []
for i in tqdm(range(len(V_z_list))):
B_list = np.load(data_path + trial_names[i] + 'B_list.npy')
B_strength_list = np.sqrt(B_list[:,0]**2 + B_list[:,0]**2)
all_times, all_gammas = B_list_growth_analyzer(times[:len(B_strength_list)], B_strength_list, window = frames_per_unitT, stride = 1)
all_gammas_all_Vz.append(np.copy(all_gammas))
# print(all_gammas_all_Vz[0].shape)
min_length = np.min([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
max_length = np.max([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
# padding the arrays for animation for the ones that converged earlier
convergence_indices = [all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))]
for i in range(len(all_gammas_all_Vz)):
if all_gammas_all_Vz[i].shape[0] < max_length:
all_gammas_all_Vz[i] = np.pad(all_gammas_all_Vz[i], ((0, max_length - all_gammas_all_Vz[i].shape[0]), (0, 0)), 'edge')
all_gammas_all_Vz = np.array(all_gammas_all_Vz)
all_times = np.array(all_times)
print(all_gammas_all_Vz.shape)
print(all_times.shape)
# adding 'converging...' and 'converged' labels
all_labels_Vz = []
for i in range(max_length):
curr_label = []
for j in range(len(label_V_z)):
if i < convergence_indices[j]:
curr_label.append(f'{label_V_z[j]} - converging...')
else:
curr_label.append(f'{label_V_z[j]} - converged')
all_labels_Vz.append(np.copy(curr_label))
all_labels_Vz = np.array(all_labels_Vz)
print(all_labels_Vz.shape)
# print(all_labels_Vz[0])
print(f'animating the evolution of the growth rate for different V_z values')
animate_to_gif_advanced(r, all_gammas_all_Vz, all_times, f"gamma_evolution_all_Vz2_new.gif", frames_per_second=frames_per_unitT, X_label=r"r [kpc]", Y_label=r"$\gamma [dimensionless]$", title=r"$\gamma$ evolution at t =", labels = all_labels_Vz.T)
0%| | 0/6 [00:00<?, ?it/s]C:\Users\adhil\AppData\Local\Temp\ipykernel_15536\1483822731.py:32: RuntimeWarning: divide by zero encountered in log slope, _ = np.polyfit(times[start:end], np.log(B_str[:, j]), 1) 100%|██████████| 6/6 [00:10<00:00, 1.72s/it]
(6, 377, 100) (377,) (377, 6) animating the evolution of the growth rate for different V_z values
Combined analysis (with disk flaring) (mistake)¶
In [ ]:
# drawing growth rate evolution for differetn Vz
V_z_list = [0, 1, 3, 5, 7]
trial_names = [f'trial{i}_' for i in range(10, 15)]
label_V_z = [f'V_z = {V_z_list[i]:.4f}' for i in range(len(V_z_list))]
mask = [V_z_list[i] % 2 == 0 for i in range(len(V_z_list))]
mask = [i % 1 == 0 for i in range(len(V_z_list))]
V_z_list = np.array(V_z_list)
V_z_list = V_z_list[mask]
label_V_z = np.array(label_V_z)
label_V_z = label_V_z[mask]
trial_names = np.array(trial_names)
trial_names = trial_names[mask]
all_gammas_all_Vz = []
for i in tqdm(range(len(V_z_list))):
B_list = np.load(data_path + trial_names[i] + 'B_list.npy')
B_strength_list = np.sqrt(B_list[:,0]**2 + B_list[:,0]**2)
all_times, all_gammas = B_list_growth_analyzer(times[:len(B_strength_list)], B_strength_list, window = frames_per_unitT, stride = 1)
all_gammas_all_Vz.append(np.copy(all_gammas))
# print(all_gammas_all_Vz[0].shape)
min_length = np.min([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
max_length = np.max([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
convergence_indices = [all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))]
for i in range(len(all_gammas_all_Vz)):
if all_gammas_all_Vz[i].shape[0] < max_length:
all_gammas_all_Vz[i] = np.pad(all_gammas_all_Vz[i], ((0, max_length - all_gammas_all_Vz[i].shape[0]), (0, 0)), 'edge')
all_gammas_all_Vz = np.array(all_gammas_all_Vz)
all_times = np.array(all_times)
print(all_gammas_all_Vz.shape)
print(all_times.shape)
all_labels_Vz = []
for i in range(max_length):
curr_label = []
for j in range(len(label_V_z)):
if i < convergence_indices[j]:
curr_label.append(f'{label_V_z[j]} - converging...')
else:
curr_label.append(f'{label_V_z[j]} - converged')
all_labels_Vz.append(np.copy(curr_label))
all_labels_Vz = np.array(all_labels_Vz)
print(all_labels_Vz.shape)
print(all_labels_Vz[0])
print(f'animating the evolution of the growth rate for different V_z values')
animate_to_gif_advanced(r, all_gammas_all_Vz, all_times, f"gamma_evolution_all_Vz_disk_flaring.gif", frames_per_second=frames_per_unitT, X_label=r"r [kpc]", Y_label=r"$\gamma [dimensionless]$", title=r"$\gamma$ evolution at t =", labels = all_labels_Vz.T)
100%|██████████| 5/5 [00:14<00:00, 2.91s/it]
(5, 777, 100) (777,) (777, 5) ['V_z = 0.0000 - converging...' 'V_z = 1.0000 - converging...' 'V_z = 3.0000 - converging...' 'V_z = 5.0000 - converging...' 'V_z = 7.0000 - converging...'] animating the evolution of the growth rate for different V_z values
Combined analysis (with disk flaring) (new- alpha variation and corrected)¶
In [ ]:
# drawing growth rate evolution for differetn Vz
V_z_list = [0, 2, 4, 6, 8, 10]
trial_names = [f'trial{i}_' for i in range(25, 31)]
label_V_z = [f'V_z = {V_z_list[i]:.4f}' for i in range(len(V_z_list))]
mask = [V_z_list[i] % 2 == 0 for i in range(len(V_z_list))]
mask = [i % 1 == 0 for i in range(len(V_z_list))]
V_z_list = np.array(V_z_list)
V_z_list = V_z_list[mask]
label_V_z = np.array(label_V_z)
label_V_z = label_V_z[mask]
trial_names = np.array(trial_names)
trial_names = trial_names[mask]
all_gammas_all_Vz = []
for i in tqdm(range(len(V_z_list))):
B_list = np.load(data_path + trial_names[i] + 'B_list.npy')
B_strength_list = np.sqrt(B_list[:,0]**2 + B_list[:,0]**2)
all_times, all_gammas = B_list_growth_analyzer(times[:len(B_strength_list)], B_strength_list, window = frames_per_unitT, stride = 1)
all_gammas_all_Vz.append(np.copy(all_gammas))
# print(all_gammas_all_Vz[0].shape)
min_length = np.min([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
max_length = np.max([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
convergence_indices = [all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))]
for i in range(len(all_gammas_all_Vz)):
if all_gammas_all_Vz[i].shape[0] < max_length:
all_gammas_all_Vz[i] = np.pad(all_gammas_all_Vz[i], ((0, max_length - all_gammas_all_Vz[i].shape[0]), (0, 0)), 'edge')
all_gammas_all_Vz = np.array(all_gammas_all_Vz)
all_times = np.array(all_times)
print(all_gammas_all_Vz.shape)
print(all_times.shape)
all_labels_Vz = []
for i in range(max_length):
curr_label = []
for j in range(len(label_V_z)):
if i < convergence_indices[j]:
curr_label.append(f'{label_V_z[j]} - converging...')
else:
curr_label.append(f'{label_V_z[j]} - converged')
all_labels_Vz.append(np.copy(curr_label))
all_labels_Vz = np.array(all_labels_Vz)
print(all_labels_Vz.shape)
# print(all_labels_Vz[0])
print(f'animating the evolution of the growth rate for different V_z values')
animate_to_gif_advanced(r, all_gammas_all_Vz, all_times, f"gamma_evolution_all_Vz_disk_flaring_new.gif", frames_per_second=frames_per_unitT, X_label=r"r [kpc]", Y_label=r"$\gamma [dimensionless]$", title=r"$\gamma$ evolution at t =", labels = all_labels_Vz.T)
0%| | 0/6 [00:00<?, ?it/s]C:\Users\adhil\AppData\Local\Temp\ipykernel_15536\1483822731.py:32: RuntimeWarning: divide by zero encountered in log slope, _ = np.polyfit(times[start:end], np.log(B_str[:, j]), 1) 100%|██████████| 6/6 [00:10<00:00, 1.82s/it]
(6, 452, 100) (452,) (452, 6) animating the evolution of the growth rate for different V_z values
Combined analysis (no disk flaring vs disk flaring)¶
In [ ]:
# drawing growth rate evolution for differetn Vz (with disk - flaring)
V_z_list1 = [0, 2, 4, 6, 8, 10]
trial_names1 = [f'trial{i}_' for i in range(25, 31)]
colors = ['r', 'g', 'b', 'c', 'm', 'k']
line_styles = ['-' for i in range(len(V_z_list1))]
V_z_list2 = [0, 2, 4, 6, 8, 10]
trial_names2 = [f'trial{i}_' for i in range(19, 25)]
colors.extend(['r', 'g', 'b', 'c', 'm', 'k'])
line_styles.extend(['--' for i in range(len(V_z_list2))])
label_V_z = [f'V_z = {V_z_list1[i]:.4f}' for i in range(len(V_z_list1))]
label_V_z.extend([f'DF: V_z = {V_z_list2[i]:.4f}' for i in range(len(V_z_list2))])
# print(label_V_z)
all_gammas_all_Vz = []
for i in tqdm(range(len(V_z_list2))):
B_list = np.load(data_path + trial_names2[i] + 'B_list.npy')
B_strength_list = np.sqrt(B_list[:,0]**2 + B_list[:,0]**2)
all_times, all_gammas = B_list_growth_analyzer(times[:len(B_strength_list)], B_strength_list, window = frames_per_unitT, stride = 1)
all_gammas_all_Vz.append(np.copy(all_gammas))
for i in tqdm(range(len(V_z_list1))):
B_list = np.load(data_path + trial_names1[i] + 'B_list.npy')
B_strength_list = np.sqrt(B_list[:,0]**2 + B_list[:,0]**2)
all_times, all_gammas = B_list_growth_analyzer(times[:len(B_strength_list)], B_strength_list, window = frames_per_unitT, stride = 1)
all_gammas_all_Vz.append(np.copy(all_gammas))
# print(all_gammas_all_Vz[0].shape)
min_length = np.min([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
max_length = np.max([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
convergence_indices = [all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))]
for i in range(len(all_gammas_all_Vz)):
if all_gammas_all_Vz[i].shape[0] < max_length:
all_gammas_all_Vz[i] = np.pad(all_gammas_all_Vz[i], ((0, max_length - all_gammas_all_Vz[i].shape[0]), (0, 0)), 'edge')
all_gammas_all_Vz = np.array(all_gammas_all_Vz)
all_times = np.linspace(0, (max_length - 1)// frames_per_unitT, max_length)
print(all_gammas_all_Vz.shape)
print(all_times.shape)
all_labels_Vz = []
for i in range(max_length):
curr_label = []
for j in range(len(label_V_z)):
if i < convergence_indices[j]:
curr_label.append(f'{label_V_z[j]} - converging...')
else:
curr_label.append(f'{label_V_z[j]} - converged')
all_labels_Vz.append(np.copy(curr_label))
all_labels_Vz = np.array(all_labels_Vz)
print(all_labels_Vz.shape)
# print(all_labels_Vz[0])
print(f'animating the evolution of the growth rate for different V_z values')
animate_to_gif_advanced(r, all_gammas_all_Vz, all_times, f"gamma_evolution_all_Vz_disk_flaring_comparison.gif", frames_per_second=frames_per_unitT, X_label=r"r [kpc]", Y_label=r"$\gamma [dimensionless]$", title=r"$\gamma$ evolution at t =", labels = all_labels_Vz.T, colors = colors, linestyles = line_styles)
0%| | 0/6 [00:00<?, ?it/s]C:\Users\adhil\AppData\Local\Temp\ipykernel_15536\1483822731.py:32: RuntimeWarning: divide by zero encountered in log slope, _ = np.polyfit(times[start:end], np.log(B_str[:, j]), 1) 100%|██████████| 6/6 [00:10<00:00, 1.81s/it] 100%|██████████| 6/6 [00:10<00:00, 1.82s/it]
(12, 452, 100) (452,) (452, 12) animating the evolution of the growth rate for different V_z values
Combined analysis (Vz = 4 km / s) constant, decreasing and decreasing then increasing¶
In [ ]:
trial_names = [f'trial{i}_' for i in [21, 31,32]]
V_z_labels = ['constant', 'decreasing till end', 'decreasing till 10 kpc']
all_gammas_all_Vz = []
for i in tqdm(range(len(trial_names))):
B_list = np.load(data_path + trial_names[i] + 'B_list.npy')
B_strength_list = np.sqrt(B_list[:,0]**2 + B_list[:,0]**2)
all_times, all_gammas = B_list_growth_analyzer(times[:len(B_strength_list)], B_strength_list, window = frames_per_unitT, stride = 1)
all_gammas_all_Vz.append(np.copy(all_gammas))
min_length = np.min([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
max_length = np.max([all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))])
convergence_indices = [all_gammas_all_Vz[i].shape[0] for i in range(len(all_gammas_all_Vz))]
for i in range(len(all_gammas_all_Vz)):
if all_gammas_all_Vz[i].shape[0] < max_length:
all_gammas_all_Vz[i] = np.pad(all_gammas_all_Vz[i], ((0, max_length - all_gammas_all_Vz[i].shape[0]), (0, 0)), 'edge')
all_gammas_all_Vz = np.array(all_gammas_all_Vz)
all_times = np.linspace(0, (max_length - 1)// frames_per_unitT, max_length)
print(all_gammas_all_Vz.shape)
print(all_times.shape)
all_labels_Vz = []
for i in range(max_length):
curr_label = []
for j in range(len(V_z_labels)):
if i < convergence_indices[j]:
curr_label.append(f'{V_z_labels[j]} - converging...')
else:
curr_label.append(f'{V_z_labels[j]} - converged')
all_labels_Vz.append(np.copy(curr_label))
all_labels_Vz = np.array(all_labels_Vz)
print(all_labels_Vz.shape)
animate_to_gif_advanced(r, all_gammas_all_Vz, all_times, f"gamma_evolution_diff_Vz_profiles.gif", frames_per_second=frames_per_unitT, X_label=r"r [kpc]", Y_label=r"$\gamma [dimensionless]$", title=r"Different $V_z$ profile's $\gamma$ evolution at t =", labels = all_labels_Vz.T)
0%| | 0/3 [00:00<?, ?it/s]C:\Users\adhil\AppData\Local\Temp\ipykernel_15536\1483822731.py:32: RuntimeWarning: divide by zero encountered in log slope, _ = np.polyfit(times[start:end], np.log(B_str[:, j]), 1) 100%|██████████| 3/3 [00:05<00:00, 1.79s/it]
(3, 402, 100) (402,) (402, 3)
STORAGE OF ALL TRIAL SETTINGS¶
In [ ]:
######################################################
# -------------------- Trial 1 --------------------- #
######################################################
# defining simulation parameters
Nt_per_unitT = 4000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial1_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
V_z_ = 10 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
np.random.seed(0)
# 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
R_omega = (neg_q_omega * (h ** 2)) / eta_t
D = R_alpha * R_omega
# # defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
return 1.5*B_0 * np.random.normal(5, 1, r.shape)
def B_r0(r, phi, z):
return - B_0 * np.random.normal(5, 1, r.shape)
# 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))
# Note: Trial to test the code and add feautures to the code, tested global gowth rate check, animation, pitch angle graphs etc. bonus: got better task 2 results! (without Vz)
######################################################
# -------------------- Trial 2 --------------------- #
######################################################
# same as trial 1 but with Vz = 10 km/s
######################################################
# -------------------- Trial 3 --------------------- #
######################################################
# defining simulation parameters
Nt_per_unitT = 4000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial3_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
V_z_ = 1 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
np.random.seed(0)
# 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
R_omega = (neg_q_omega * (h ** 2)) / eta_t
D = R_alpha * R_omega
# # defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
return 1.5*B_0 * np.random.normal(5, 1, r.shape)
def B_r0(r, phi, z):
return - B_0 * np.random.normal(5, 1, r.shape)
# Note: V_z = 1, same other paramters
######################################################
# -------------------- Trial 4 --------------------- #
######################################################
# defining simulation parameters
Nt_per_unitT = 4000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial4_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
V_z_ = 2 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
np.random.seed(0)
# 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
R_omega = (neg_q_omega * (h ** 2)) / eta_t
D = R_alpha * R_omega
# # defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
return 1.5*B_0 * np.random.normal(5, 1, r.shape)
def B_r0(r, phi, z):
return - B_0 * np.random.normal(5, 1, r.shape)
# Note: V_z = 2, same other paramters
######################################################
# -------------------- Trial 5 --------------------- #
######################################################
# defining simulation parameters
Nt_per_unitT = 4000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial5_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
V_z_ = 4 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
np.random.seed(0)
# 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
R_omega = (neg_q_omega * (h ** 2)) / eta_t
D = R_alpha * R_omega
# # defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
return 1.5*B_0 * np.random.normal(5, 1, r.shape)
def B_r0(r, phi, z):
return - B_0 * np.random.normal(5, 1, r.shape)
######################################################
# -------------------- Trial 6 --------------------- #
######################################################
# note: same as trial 5 but with Vz = 6
######################################################
# -------------------- Trial 7 --------------------- #
######################################################
# note: same as trial 5 but with Vz = 6.7
######################################################
# -------------------- Trial 8 --------------------- #
######################################################
# note: same as trial 5 but with Vz = 8.5
######################################################
# -------------------- Trial 9 --------------------- #
######################################################
# note: same as trial 5 but with Vz = 8
######################################################
# -------------------- Trial 10 -------------------- #
######################################################
# note: adding disk flaring
h_ = 0.4 * u.kpc # kpc
h = (h_ / h_).value
h = h * ((1 + ((r / r_f) ** 2)) ** 0.5)
######################################################
# -------------------- Trial 11 -------------------- #
######################################################
# note: adding disk flaring and Vz = 1 km/s
######################################################
# -------------------- Trial 12 -------------------- #
######################################################
# note: adding disk flaring and Vz = 3 km/s (converged at less than 40)
######################################################
# -------------------- Trial 13 -------------------- #
######################################################
# note: adding disk flaring and Vz = 5 km/s (converged at 40)
######################################################
# -------------------- Trial 14 -------------------- #
######################################################
# note: adding disk flaring and Vz = 7 km/s (did not converge at 40)
######################################################
# -------------------- Trial 15 -------------------- #
######################################################
# note: adding disk flaring with r_h_ = 10 kpc and Vz = 7 km/s (weird eigen function)
h_ = 0.4 * u.kpc # kpc
h = (h_ / h_).value
r_h_ = 10 * u.kpc # kpc
r_h = (r_h_ / h_).value
h = h * ((1 + ((r / r_h) ** 2)) ** 0.5)
######################################################
# -------------------- Trial 16 -------------------- #
######################################################
# note: adding disk flaring with r_h_ = 10 kpc and Vz = 0 km/s
######################################################
# ---------------- Trial 15_new -------------------- #
######################################################
# removed the mistake of multiplying the q with h^2 in R_omega, disk flaring was affecting the eigen function; ran with Vz = 7 km/s, r_h = 10 kpc
R_omega = neg_q_omega # * (h ** 2)) / eta_t # sjnce all the terms are already dimensionless
######################################################
# -------------------- Trial 17 -------------------- #
######################################################
# defining simulation parameters
Nt_per_unitT = 5000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial18_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
dt = t_d / Nt_per_unitT
rand_seed = 0
# checking the courant condition
C_val = (2 * eta_t * dt) / ((dr) ** 2)
if C_val > 1:
print(f"Courant condition is not satisfied. C = {C_val:.4f} > 1")
print(f"Need eta_t < {eta_t / C_val:.4f} (dimensionless)")
print(f"or Nt_per_unitT > {int(Nt_per_unitT * C_val)}")
print(f'or need Nr > {int(np.sqrt(((Nr - 1)**2) / C_val) + 1)}')
else:
print(f"Courant condition is satisfied. C = {C_val:.4e} < 1")
if C_val > 0.01:
print(f"Warning: C is close to 1. Consider increasing Nt_per_unitT or decreasing Nr")
# defining the outflow
V_z_ = 10 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
# disk - flaring
h_ = 0.4 * u.kpc # kpc
h = (h_ / h_).value
r_h_ = 10 * u.kpc # kpc
r_h = (r_h_ / h_).value
h = h * ((1 + ((r / r_h) ** 2)) ** 0.5)
# 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
R_omega = neg_q_omega # * (h ** 2)) / eta_t # sjnce all the terms are already dimensionless
D = R_alpha * R_omega
# defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
np.random.seed(rand_seed)
return -1.5* 1e-4 * B_0 * np.random.normal(5, 1, r.shape)
def B_r0(r, phi, z):
np.random.seed(rand_seed)
return B_0 * 1e-4 * np.random.normal(5, 1, r.shape)
# return np.zeros(r.shape)
# note reduced seed field, same seed field for Br and Bphi, Vz = 0, disk flaring with r_h = 10 kpc
######################################################
# -------------------- Trial 18 -------------------- #
######################################################
# note: V-z = 10 , other parameters same as trial 17 (to compare to sukanta bhaiyya's results)
######################################################
# -------------------- Trial 19 -------------------- #
######################################################
# note: no disk flaring, alpha variation with r, Vz = 0 km/s,
# defining simulation parameters
Nt_per_unitT = 5000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial19_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
dt = t_d / Nt_per_unitT
rand_seed = 0
# checking the courant condition
C_val = (2 * eta_t * dt) / ((dr) ** 2)
if C_val > 1:
print(f"Courant condition is not satisfied. C = {C_val:.4f} > 1")
print(f"Need eta_t < {eta_t / C_val:.4f} (dimensionless)")
print(f"or Nt_per_unitT > {int(Nt_per_unitT * C_val)}")
print(f'or need Nr > {int(np.sqrt(((Nr - 1)**2) / C_val) + 1)}')
else:
print(f"Courant condition is satisfied. C = {C_val:.4e} < 1")
if C_val > 0.01:
print(f"Warning: C is close to 1. Consider increasing Nt_per_unitT or decreasing Nr")
# defining the outflow
V_z_ = 0 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
# disk - flaring setting
h_ = 0.4 * u.kpc # kpc
h = (h_ / h_).value
# r_h_ = 10 * u.kpc # kpc
# r_h = (r_h_ / h_).value
# h = h * ((1 + ((r / r_h) ** 2)) ** 0.5)
# defining the spatial dependent omega
omega = omega_0 / (np.sqrt(1 + (r / r_omega) ** 2))
# defining the spatial dependent alpha
alpha = omega * (l_kpc ** 2) / (h)
# defining q
neg_q_omega = r * first_derivative_dYdx_FD(omega, dr, order = 6, ghost_zone_type='smooth')
R_alpha = alpha
R_omega = neg_q_omega # * (h ** 2)) / eta_t # sjnce all the terms are already dimensionless
D = R_alpha * R_omega
# defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
np.random.seed(rand_seed)
return -1.5* 1e-4 * B_0 * np.random.normal(5, 1, r.shape) * np.sin((r - r_i) * np.pi / (r_f - r_i))
def B_r0(r, phi, z):
np.random.seed(rand_seed)
return B_0 * 1e-4 * np.random.normal(5, 1, r.shape) * np.sin((r - r_i) * np.pi / (r_f - r_i))
# return np.zeros(r.shape)
######################################################
# -------------------- Trial 20 -------------------- #
######################################################
# note: no disk flaring, alpha variation with r, Vz = 2 km/s
######################################################
# -------------------- Trial 21 -------------------- #
######################################################
# note: no disk flaring, alpha variation with r, Vz = 4 km/s
######################################################
# -------------------- Trial 22 -------------------- #
######################################################
# note: no disk flaring, alpha variation with r, Vz = 6 km/s
######################################################
# -------------------- Trial 23 -------------------- #
######################################################
# note: no disk flaring, alpha variation with r, Vz = 8 km/s
######################################################
# -------------------- Trial 24 -------------------- #
######################################################
# note: no disk flaring, alpha variation with r, Vz = 10 km/s
######################################################
# -------------------- Trial 25 -------------------- #
######################################################
# note: disk flaring, alpha variation with r, Vz = 0 km/s
# defining simulation parameters
Nt_per_unitT = 5000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial25_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
dt = t_d / Nt_per_unitT
rand_seed = 0
# checking the courant condition
C_val = (2 * eta_t * dt) / ((dr) ** 2)
if C_val > 1:
print(f"Courant condition is not satisfied. C = {C_val:.4f} > 1")
print(f"Need eta_t < {eta_t / C_val:.4f} (dimensionless)")
print(f"or Nt_per_unitT > {int(Nt_per_unitT * C_val)}")
print(f'or need Nr > {int(np.sqrt(((Nr - 1)**2) / C_val) + 1)}')
else:
print(f"Courant condition is satisfied. C = {C_val:.4e} < 1")
if C_val > 0.01:
print(f"Warning: C is close to 1. Consider increasing Nt_per_unitT or decreasing Nr")
# defining the outflow
V_z_ = 0 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
# disk - flaring setting
h_ = 0.4 * u.kpc # kpc
h = (h_ / h_).value
r_h_ = 10 * u.kpc # kpc
r_h = (r_h_ / h_).value
h = h * ((1 + ((r / r_h) ** 2)) ** 0.5)
# defining the spatial dependent omega
omega = omega_0 / (np.sqrt(1 + (r / r_omega) ** 2))
# defining the spatial dependent alpha
alpha = omega * (l_kpc ** 2) / (h)
# defining q
neg_q_omega = r * first_derivative_dYdx_FD(omega, dr, order = 6, ghost_zone_type='smooth')
R_alpha = alpha
R_omega = neg_q_omega # * (h ** 2)) / eta_t # sjnce all the terms are already dimensionless
D = R_alpha * R_omega
# defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
np.random.seed(rand_seed)
return -1.5* 1e-4 * B_0 * np.random.normal(5, 1, r.shape) * np.sin((r - r_i) * np.pi / (r_f - r_i))
def B_r0(r, phi, z):
np.random.seed(rand_seed)
return B_0 * 1e-4 * np.random.normal(5, 1, r.shape) * np.sin((r - r_i) * np.pi / (r_f - r_i))
# return np.zeros(r.shape)
######################################################
# -------------------- Trial 26 -------------------- #
######################################################
# note: disk flaring, alpha variation with r, Vz = 2 km/s
######################################################
# -------------------- Trial 27 -------------------- #
######################################################
# note: disk flaring, alpha variation with r, Vz = 4 km/s
######################################################
# -------------------- Trial 28 -------------------- #
######################################################
# note: disk flaring, alpha variation with r, Vz = 6 km/s
######################################################
# -------------------- Trial 29 -------------------- #
######################################################
# note: disk flaring, alpha variation with r, Vz = 8 km/s
######################################################
# -------------------- Trial 30 -------------------- #
######################################################
# note: disk flaring, alpha variation with r, Vz = 10 km/s
######################################################
# -------------------- Trial 31 -------------------- #
######################################################
# note: no disk flaring, alpha variation with r, Vz = decreasing form
# defining simulation parameters
Nt_per_unitT = 5000 # Number of time grid points per unit time
Nr = 100 # Number of spatial grid points
order = 6 # order of the finite difference scheme
ghost_zone_type = 'anti-symmetric' # type of the ghost zone
frames_per_unitT = 25 # number of frames per unit time
iterations_to_save_plot = [int(i * (Nt_per_unitT / frames_per_unitT)) for i in range(0, frames_per_unitT + 1)] # saving only for the frames
trial_name = 'trial31_'
global_growth_rate_tolerance = 1e-2
check_tolerance_after = 10
check_tolerance_every = 1
growth_rate_calc_using = 0.1
patience = 3
max_iterations = 40
omega_0_start = omega_0 * 0.5
omega_0end = omega_0 * 1.5
omega_0step = 10
interpolation = True
r = np.linspace(r_i, r_f, Nr)
dr = (r_f - r_i) / (Nr - 1)
dt = t_d / Nt_per_unitT
rand_seed = 0
# checking the courant condition
C_val = (2 * eta_t * dt) / ((dr) ** 2)
if C_val > 1:
print(f"Courant condition is not satisfied. C = {C_val:.4f} > 1")
print(f"Need eta_t < {eta_t / C_val:.4f} (dimensionless)")
print(f"or Nt_per_unitT > {int(Nt_per_unitT * C_val)}")
print(f'or need Nr > {int(np.sqrt(((Nr - 1)**2) / C_val) + 1)}')
else:
print(f"Courant condition is satisfied. C = {C_val:.4e} < 1")
if C_val > 0.01:
print(f"Warning: C is close to 1. Consider increasing Nt_per_unitT or decreasing Nr")
# defining the outflow
V_z_ = 4 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
r_Vz_= 10 * u.kpc # kpc
r_Vz = (r_Vz_ / h_).value
fraction = 0.5
V_z = V_z * (1 + fraction * ((((r - r_Vz) / r_Vz) ** 2) - 1))
# disk - flaring setting
h_ = 0.4 * u.kpc # kpc
h = (h_ / h_).value
# r_h_ = 10 * u.kpc # kpc
# r_h = (r_h_ / h_).value
# h = h * ((1 + ((r / r_h) ** 2)) ** 0.5)
# defining the spatial dependent omega
omega = omega_0 / (np.sqrt(1 + (r / r_omega) ** 2))
# defining the spatial dependent alpha
alpha = omega * (l_kpc ** 2) / (h)
# defining q
neg_q_omega = r * first_derivative_dYdx_FD(omega, dr, order = 6, ghost_zone_type='smooth')
R_alpha = alpha
R_omega = neg_q_omega # * (h ** 2)) / eta_t # sjnce all the terms are already dimensionless
D = R_alpha * R_omega
# defining the initial conditions as random gaussian field multiplied by B0
def B_phi0(r, phi, z):
np.random.seed(rand_seed)
return -1.5* 1e-4 * B_0 * np.random.normal(5, 1, r.shape) * np.sin((r - r_i) * np.pi / (r_f - r_i))
def B_r0(r, phi, z):
np.random.seed(rand_seed)
return B_0 * 1e-4 * np.random.normal(5, 1, r.shape) * np.sin((r - r_i) * np.pi / (r_f - r_i))
# return np.zeros(r.shape)
######################################################
# -------------------- Trial 32 -------------------- #
######################################################
# note: no disk flaring, alpha variation with r, Vz = decreasing then increasing
V_z_ = 4 * (u.km / u.s) # km/s
V_z = (V_z_.to(u.kpc / u.Gyr) * (t_d_ / h_)).value
r_Vz_= 15 * u.kpc # kpc
r_Vz = (r_Vz_ / h_).value
fraction = 0.5
V_z = V_z * (1 + fraction * ((((r - r_Vz) / r_Vz) ** 2) - 1))