Computing seismic kernels ========================= This tutorial explains how to use ``persephone`` to compute seismic rotational kernels of evolutionary models obtained through the interface. The user will find an extensive description and discussion of the rotational kernel formalism presented in this tutorial in Chapter 8 of the book by J. Christensen-Dalsgaard, `Lecture Notes on Stellar Oscillations `__. .. code:: ipython3 import persephone as pph import persephone.kernels.rotation as pph_rot import pygyre import os pph.__version__ .. parsed-literal:: '0.1.1' Computing the stellar model --------------------------- In this tutorial, we will consider a :math:`1.3~\rm M_\odot` stellar model along its evolution. .. code:: ipython3 mesa_version = "r15140" parameters = {"controls":{"initial_mass":"1.3", "initial_z":"0.0134"}} model_dir = "kernel_example" The ``create_mesa_model`` function provides the wrapper necessary to run a single model in its working directory. .. code:: ipython3 pph.grids.create_mesa_model (model_dir=model_dir, mesa_version=mesa_version, parameters=parameters, run_mesa=True, verbose=False) Low-frequency modes computation with GYRE ----------------------------------------- We want to follow the evolution of the properties of low-frequencies gravity modes (g modes) of the model along stellar evolution. To do this, we need to use a dedicated GYRE input template file, as the standard input file included in ``persephone`` is dedicated to compute the properties of high-frequency acoustic modes (p modes). The customed template file to use is shown below. **Note:** for the sake of consistency between filenames templates and further computation of 2D-kernels, we use GYRE to compute :math:`m=1` modes, even if we do not actually include rotation in GYRE inputs (see the ``rot`` namelist below). Therefore, here, :math:`m=0` and :math:`m=1` mode computed by GYRE should have the exact same properties. .. code:: ipython3 pph.grids.show_kernel_example_gyre_input () .. parsed-literal:: &constants / &model model_type = 'EVOL' ! Obtain stellar structure from an evolutionary model file_format = 'MESA' ! File format of the evolutionary model add_center = .FALSE. / &mode l = 1 m = 1 n_pg_min = -5 n_pg_max = 5 / &osc inner_bound = 'REGULAR' outer_bound = 'VACUUM' ! Use a zero-pressure outer mechanical boundary condition / &rot Omega_rot_source = 'UNIFORM' Omega_rot = 0. Omega_rot_units = 'UHZ' / &num diff_scheme = 'COLLOC_GL4' ! 4th-order collocation scheme for difference equations / &scan grid_type = 'INVERSE' freq_min = 10 freq_max = 300 n_freq = 500 freq_units = 'UHZ' / &grid w_osc = 10 ! Oscillatory region weight parameter w_exp = 2 ! Exponential region weight parameter w_ctr = 10 ! Central region weight parameter / &ad_output summary_item_list = 'l,m,n_pg,n_p,n_g,freq,freq_units, E, E_norm, E_p, E_g, R_star, M_star, Delta_p, Delta_g' detail_item_list = 'l,m,n_pg,n_p,n_g,omega,freq,x,xi_r, xi_h,c_1,As,V_2,Gamma_1,P,rho, eul_P,prop_type, E, E_norm, E_p, E_g, R_star, M_star, Delta_p, Delta_g, M_r' freq_units = 'UHZ' freq_frame = 'INERTIAL' / &nad_output / By default, in order to save disk memory, ``persephone`` uses a GYRE template that does not save mode details file. However, here, they are required to compute the seismic kernel, the procedure to save the detail files can be simply included in the interface by setting ``add_details=True``. In this case, ``persephone`` will add to the input namelists above a ``detail_template`` argument and will take care of creating the directories where the detail files are stored. .. code:: ipython3 _gyre_in = pph.grids.get_kernel_example_gyre_input () with _gyre_in as template_file : pph.grids.analyse_gyre (model_dir, run=True, template_file=template_file, add_details=True, verbose=False) Kernels along evolution ----------------------- This list of details directory can then be accessed through the ``get_listdir_gyre_details`` function. Selecting a details directory, the details file can be obtained using ``get_detail_file``. .. code:: ipython3 list_details = pph.grids.get_listdir_gyre_details (model_dir, sort=True) list_details .. parsed-literal:: ['kernel_example/LOGS/profile1_details', 'kernel_example/LOGS/profile2_details', 'kernel_example/LOGS/profile3_details', 'kernel_example/LOGS/profile4_details', 'kernel_example/LOGS/profile5_details', 'kernel_example/LOGS/profile6_details', 'kernel_example/LOGS/profile7_details', 'kernel_example/LOGS/profile8_details', 'kernel_example/LOGS/profile9_details'] Making the approximation of shellular rotation, it is possible to compute one-dimensional kernels that do not depend on the azimuthal number :math:`m`. In this example, we focus on the :math:`n=1` g mode of the models (if computed, see Note below). **Note:** Depending on the stellar structure, not all the models we computed have the :math:`n=1`, :math:`\ell=1` g mode in the considered frequency range, or g modes at all (if they are fully convective). For this reason, the ``file_exists`` boolean status is used to avoid ``FileNotFoundError`` exceptions. .. code:: ipython3 for details_dir in list_details : filename, file_exists = pph.grids.get_detail_file (details_dir, n=-1, l=1, m=1) if file_exists : d = pygyre.read_output(filename) beta_nl, K_nl = pph_rot.compute_shellular_kernel (d["x"], d["rho"], d["xi_r"], d["xi_h"], l=d.meta["l"]) title = r"{} ($\nu={:.1f}$ $\mu$Hz)".format (os.path.basename (details_dir), d.meta["freq"].real) fig = pph_rot.plot_shellular_kernel (d["x"], K_nl, color="red", lw=2, title=title, figsize=(6,4)) .. image:: kernel_tutorial_files/kernel_tutorial_17_0.png .. image:: kernel_tutorial_files/kernel_tutorial_17_1.png .. image:: kernel_tutorial_files/kernel_tutorial_17_2.png .. image:: kernel_tutorial_files/kernel_tutorial_17_3.png .. image:: kernel_tutorial_files/kernel_tutorial_17_4.png .. image:: kernel_tutorial_files/kernel_tutorial_17_5.png .. image:: kernel_tutorial_files/kernel_tutorial_17_6.png Finally, we can compute two-dimensional rotational kernels as a function of radius and latitude. .. code:: ipython3 for details_dir in list_details : filename, file_exists = pph.grids.get_detail_file (details_dir, n=-1, l=1, m=1) if file_exists : d = pygyre.read_output(filename) (r, theta, beta_nlm, K_nlm) = pph_rot.compute_general_kernel (d["x"], d["rho"], d["xi_r"], d["xi_h"], l=d.meta["l"], m=d.meta["m"], theta_res=180, theta_0=0, theta_1=180, ) title = r"{} ($\nu={:.1f}$ $\mu$Hz)".format (os.path.basename (details_dir), d.meta["freq"].real) fig = pph_rot.plot_2D_kernel (r, theta, K_nlm, cmap="magma", title=title) .. image:: kernel_tutorial_files/kernel_tutorial_19_1.png .. image:: kernel_tutorial_files/kernel_tutorial_19_2.png .. image:: kernel_tutorial_files/kernel_tutorial_19_3.png .. image:: kernel_tutorial_files/kernel_tutorial_19_4.png .. image:: kernel_tutorial_files/kernel_tutorial_19_5.png .. image:: kernel_tutorial_files/kernel_tutorial_19_6.png .. image:: kernel_tutorial_files/kernel_tutorial_19_7.png