{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", " This is an example file to run the model with standard set of parameters.\n", " A small visualization routine is also provided to plot and check the outputs.\n", "\n", " Currently we are developing this model and more functionalities will be added in due course.\n", "\n", " Author: Soumyaranjan Dash\n", " Email: dash.soumya922@gmail.com\n", " Date: 15th Jan 2025\n", "\n", "\"\"\"\n", "# Import basic packages for simulation and visualization\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from tqdm import tqdm\n", "\n", "from matplotlib import rc\n", "import matplotlib.style\n", "## Plotting canvas properties.\n", "params = {'legend.fontsize': 12,\n", " 'axes.labelsize': 10,\n", " 'axes.titlesize': 10,\n", " 'xtick.labelsize' :10,\n", " 'ytick.labelsize': 10,\n", " 'grid.color': 'k',\n", " 'grid.linestyle': ':',\n", " 'grid.linewidth': 0.5,\n", " 'mathtext.fontset' : 'stix',\n", " 'mathtext.rm' : 'DejaVu serif',\n", " 'font.family' : 'DejaVu serif',\n", " 'font.serif' : \"Times New Roman\", # or \"Times\" \n", " }\n", "matplotlib.rcParams.update(params)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "os.chdir('../../')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import packagess from the SFT model\n", "from src.grid import create_grid\n", "from src.initial_conditions import initialize_field\n", "from src.transport_profiles import meridional_flow, differential_rotation\n", "from src.time_step import calculate_time_step\n", "from src.diffusion import calculate_diffusion\n", "from src.advection import calculate_advection" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define and initialize some essential model input parameters.\n", "\n", "## Grid creation (180x360) Uniform in latitude-longitude\n", "grid_sft = create_grid(180,360)\n", "## Meridional flow profile\n", "mf_ = meridional_flow(grid_sft.copy())\n", "## Differential rotation profile\n", "dr_ = differential_rotation(grid_sft.copy())\n", "## Initial magnetic field (dipole)\n", "## There is a choice between Dipole or HMI CR map for initial condition\n", "## For HMI map, set 'read' in place of 'dipole'.\n", "field = initialize_field(grid_sft.copy(), 'dipole')\n", "## Magnetic diffusivity\n", "diffusivity = 2.5 * 10**8 # cm^2/s\n", "## Time step calculation\n", "ts, ndt = calculate_time_step(grid_sft.copy(), diffusivity)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Create basic grid for visualization\n", "colatitude = grid_sft['colatitude']\n", "longitude = grid_sft['longitude']\n", "delta_theta = grid_sft['dtheta']\n", "delta_phi = grid_sft['dphi']\n", "solar_radius = 6.955 * 10**8 # Solar radius in meters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Visualize the flow profiles (set 1 to show the plot)\n", "if 0:\n", " plt.figure(figsize=[9,3])\n", " ax1 = plt.subplot(121)\n", " ax1.plot(np.rad2deg(grid_sft['colatitude']-np.pi/2),mf_[:,180],color='k',label='Meridional flow',lw=2)\n", " ax1.set_xlabel('Latitude [deg]')\n", " ax1.set_ylabel('Velocity [m/s]')\n", " ax1.set_title('Meridional flow')\n", " ax1.set_xlim([-90,90])\n", " ax1.set_xticks(np.arange(-90,91,30))\n", "\n", " ax2 = plt.subplot(122)\n", " ax2.plot(np.rad2deg(grid_sft['colatitude']-np.pi/2),dr_[:,180]*np.sin(grid_sft['colatitude'])*solar_radius,c='k',label='Differential rotation',lw=2)\n", " ax2.set_xlabel('Latitude [deg]')\n", " ax2.set_ylabel('Velocity [m/s]')\n", " ax2.set_title('Differential rotation')\n", " ax2.set_xlim([-90,90])\n", " ax2.set_xticks(np.arange(-90,91,30))\n", "\n", " plt.subplots_adjust(wspace=0.3)\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Run the SFT evolution.\n", "\n", "## Total number of days to run the simulation\n", "## Time-step is rounded off to fit exactly into one day\n", "num_days = 15*365 # Example total time [days]\n", "\n", "## Grid properties\n", "num_theta = grid_sft['colatitude'].size\n", "num_phi = grid_sft['longitude'].size\n", "\n", "## Initialize arrays to store the data\n", "bfly_data = np.zeros((num_days+1,num_theta))\n", "all_br_data = np.zeros((num_days+1,num_theta,num_phi))\n", "all_br_data[0,:,:] = field.copy()\n", "bfly_data[0,:] = np.mean(field,axis=1)\n", "\n", "## Define temporary arrays for updating the field\n", "B_temp = field.copy()\n", "B_temp_update = np.zeros_like(field)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Time loop for evolution\n", "for t in tqdm(range(1,num_days+1),desc='Simulation days: '):\n", " delta_t = ts\n", " for steps in range(ndt):\n", " # Calculate the diffusion term\n", " B_temp_diff = calculate_diffusion(B_temp, diffusivity, grid_sft.copy())\n", " # Calculate the advection term\n", " B_temp_adv = calculate_advection(B_temp, dr_, mf_, grid_sft.copy())\n", " \n", " # Update magnetic field using all terms\n", " B_temp_update = B_temp + delta_t * (1.0*B_temp_diff - 1.0*B_temp_adv)\n", "\n", " # Apply periodic boundary conditions in the phi direction\n", " B_temp_update[:, 0] = B_temp_update[:, -2] # First column matches second-to-last column\n", " B_temp_update[:, -1] = B_temp_update[:, 1] # Last column matches second column\n", "\n", " # Apply open boundary conditions in the theta direction\n", " B_temp_update[0, :] = B_temp_update[1, :] # Northern boundary (pole)\n", " B_temp_update[-1, :] = B_temp_update[-2, :] # Southern boundary (pole)\n", " B_temp = B_temp_update.copy()\n", " # Save the butterfly diagram\n", " bfly_data[t,:] = np.mean(B_temp,axis=1)\n", " all_br_data[t,:,:] = B_temp.copy()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Visualize the butterfly diagram (set 1 to run the following script)\n", "if 0:\n", " plt.figure(figsize=[7,3])\n", " ax1 = plt.subplot(111)\n", " bmax = 5\n", " pm1 = ax1.pcolormesh(np.arange(num_days+1),np.rad2deg(grid_sft['colatitude']-np.pi/2),bfly_data.T,cmap='bwr',vmax=bmax,vmin=-bmax)\n", " ax1.set_xlabel('Time [days]')\n", " ax1.set_ylabel('Latitude [deg]')\n", " plt.colorbar(pm1)\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Calculate the total unsigned flux from the magnetic field\n", "def calc_usflx(all_br_data):\n", " temp_bsinth = np.zeros_like(all_br_data[0,:,:])\n", " usflx = np.zeros(num_days+1)\n", " for i in range(num_days+1):\n", " for i1 in range(longitude.shape[0]):\n", " temp_bsinth[:,i1] = np.abs(all_br_data[i,:,i1]) * np.sin(colatitude)\n", " usf_1d = np.sum(temp_bsinth) * delta_theta * delta_phi * (solar_radius*1e2)**2\n", " usflx[i] = usf_1d\n", " return usflx" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Calculate the total unsigned flux\n", "usflx_diff = calc_usflx(all_br_data)\n", "\n", "## Plot the USFLUX variation (set 1 to run the following script)\n", "if 0:\n", " plt.plot(np.arange(num_days+1),usflx_diff)\n", " plt.xlabel('Time [days]')\n", " plt.ylabel('Total unsigned flux [Mx]')\n", " plt.title('Total unsigned flux')\n", " plt.show()" ] } ], "metadata": { "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 2 }