
import pandas as pd
import numpy as np
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rc
import scipy
# from scipy import interpolate
import sys
import os 
sys.path.append(os.environ['HOME'] + '/ns_plot/scripts')
from plot_ascii import combine_files, plot_0D_data
from plot_hdf5 import *
from helper import interpol, find_correct_path, default_rc_params, find_tan_using_ellipse
from configparser import ConfigParser, ExtendedInterpolation
import json
import glob
sys.path.append(os.environ['HOME'] + "/ns_plot/scripts/inverse_distance_weighting")
# import idw


default_rc_params()


sim_log = ConfigParser()
sim_log.optionxform = str
config = ConfigParser()
config.optionxform = str
sim_log.read(os.environ["SIM_LOG_PATH"])
# sim_log.read("/glade/u/home/yufengl/ns_plot/ucar/sim_log.txt")
# this parameter will be set by the master run scripts
run_name = os.environ["RUN_NAME"]
print(os.environ["SIM_LOG_PATH"], os.environ["RUN_NAME"])
root_dir = os.path.realpath(os.path.join("/glade/scratch/yufengl", "simulations", run_name))
# num_dir = int(sim_log[run_name]['num_runs'])  # -1 case
num_dir = int(os.environ["NUM_DIR_RUNNING"])
perturb_type = sim_log[run_name]["perturb"]
sim = sim_log[run_name]['parfile']
machine = sim_log["system"]["machine"]
rns_model = sim_log[run_name]['rns_model']
        
config.read(f"/glade/u/home/yufengl/ns_plot/params/{rns_model}_params.txt")
Req = Rx = float(config['RNS_param']["Rx"])
Ry = float(config['RNS_param']["Ry"])
RHO0 = float(config['RNS_param']["RHO0"])
PRESS_C = float(config['RNS_param']["PRESS_C"])
OMEGA_rot = float(config['RNS_param']["OMEGA_rot"])
M_ADM = float(config['RNS_param']["M_ADM"])
M_0 = float(config['RNS_param']["M_0"])
J_cocal = float(config['RNS_param']["J_cocal"])
t_P = 2 * numpy.pi / OMEGA_rot
t_dyn = 1/ np.sqrt(RHO0)
num_grid = float(sim_log[run_name]['num_grid'])
# t_P = t_dyn
dx = dy = dz = Req / num_grid * (2**4)
r_gw = json.loads(config['grid']['r_gw'])

output_num = num_dir
# if plot_perturb:
#   sim = sim + f'_{perturb_type}'

/glade/u/home/yufengl/ns_plot/ucar/sim_log.txt C01s17_switch_rxry_0419


from pathlib import Path


# plot_output_dir = Path("./", run_name, f"output-{output_num:04}")
plot_output_dir = Path("./", f"output-{output_num:04}")


if not os.path.exists(plot_output_dir):
    os.makedirs(plot_output_dir)
    os.mkdir(Path(plot_output_dir, 'scalar'))
    os.mkdir(Path(plot_output_dir, '1D'))
    os.mkdir(Path(plot_output_dir, '2D'))
    print("created all directories for", plot_output_dir)
print(plot_output_dir)

created all directories for output-0000
output-0000


print(root_dir, num_dir)

/glade/scratch/yufengl/simulations/C01s17_switch_rxry_0419 0


path_template = find_correct_path(root_dir, 0, sim)

Didn't find path with this one: /glade/scratch/yufengl/simulations/C01s17_switch_rxry_0419/output-{num_dir:04}/{filename} 
trying new one
found path /glade/scratch/yufengl/simulations/C01s17_switch_rxry_0419/output-{num_dir:04}/C01s17_vel_switch_rxry/{filename}


constraints = combine_files("baikal-aux_variables.norm2.asc", path_template, num_dir, sim_name=sim)


# Hamiltonian constraint
print("plotting Hamiltonian constraint")
ham_norm = np.array(constraints[:, :3], copy=True)
ham_norm[:, 2] *= dx**2
mom_norm = np.zeros_like(ham_norm)
mom_norm = np.array(constraints[:, :3], copy=True)
mom_norm[:, 2] = np.linalg.norm(constraints[:, 3:]*dx**3, axis=1)


plt.figure(figsize=(8, 8))
plt.semilogy(ham_norm[:, 1] / t_P, ham_norm[:, 2], label=sim)
plt.xlabel(r'$t / t_P$')
plt.ylabel(r'log $|H|$')
plt.ylim([1E-5, 1])
plt.legend()
# plt.tight_layout()
plt.savefig(Path(plot_output_dir, 'scalar', sim +f'_{perturb_type}_output_{output_num}'  + '_ml_ham_norm2.pdf'), bbox_inches="tight")

# Momentum constraint
print("plotting Momentum constraint")
plt.figure(figsize=(8, 8))
plt.semilogy(mom_norm[:, 1] / t_P, mom_norm[:, 2],  label=sim)
plt.xlabel(r'$t / t_P$')
plt.ylabel(r'log $|M|$')
plt.ylim([1E-9, 1])
plt.legend()
# plt.tight_layout()
# plt.show()
plt.savefig(Path(plot_output_dir, 'scalar', sim +f'_{perturb_type}_output_{output_num}'  + '_ml_mom_norm2.pdf'), bbox_inches="tight")

plotting Hamiltonian constraint
plotting Momentum constraint


# # COM movement
# print("plotting center of mass movement")
# if not plot_perturb:
#     center_loc = combine_files('hydro_analysis-hydro_analysis_rho_max_loc.maximum.asc', root_dir, num_dir)
# else:
#     center_loc = combine_files('hydro_analysis-hydro_analysis_rho_max_loc.maximum.asc', root_dir_perturb, num_dir_perturb)


primitives = combine_files("illinoisgrmhd-grmhd_primitives_allbutbi.maximum.asc", path_template, num_dir, sim_name=sim)


rho_max = primitives[:, :3]
plt.figure(figsize=(10,10))
plt.plot(rho_max[:, 1] / t_P, rho_max[:, -1]/RHO0)
# plt.plot(rho_max_low[:, 1] * np.sqrt(RHO0), rho_max_low[:, -1]/RHO0)
# plt.plot(rho_un_press[:, 1] * np.sqrt(RHO0), rho_un_press[:, -1]/RHO0)
# plt.scatter(rho_max_without[:, 1]/(2*np.pi/0.03), rho_max_without[:, -1]/RHO0)
# plt.xlim([0, 0.02]) # this is happening since B = 0.5
# plt.ylim([0.99, 1.01])
plt.xlabel(r't/$t_{P}$')
plt.ylabel(r'$\rho_c / \rho_{c,t=0}$')
plt.savefig(Path(plot_output_dir, 'scalar', sim + '_rho_c.pdf'))


# max pressure vs. time
press_max = primitives[:, [0, 1, 3]]
print("plotting maximum pressure vs. time")
plt.figure(figsize=(8, 8))
plt.plot(press_max[:, 1] / t_P, press_max[:, -1] / PRESS_C, label=sim)
# plt.ylim([0, 1.1])
plt.xlabel(r'time $t/t_P$')
plt.ylabel(r'$p_c / p_{c,t=0}$')
plt.legend()
# plt.ylim([0.5, 1.5])
# plt.tight_layout()
plt.savefig(Path(plot_output_dir, 'scalar', sim +f'_{perturb_type}_output_{output_num}'  + '_press_c.pdf'), bbox_inches="tight")

plotting maximum pressure vs. time


cons_sum = combine_files("illinoisgrmhd-grmhd_conservatives.sum.asc", path_template, num_dir)


rho_sum = cons_sum[:, [0, 1, 2]]
print("plotting maximum pressure vs. time")
plt.figure(figsize=(8, 8))
plt.plot(rho_sum[:, 1] / t_P, rho_sum[:, 2] 
             * dx * dy * dz / M_0, label=sim)
# plt.ylim([0, 1.1])
plt.xlabel(r"$t / t_P$")
plt.ylabel(r"$M_0 / M_{0, cocal}$")
plt.legend()
plt.ylim([-1.5, 1.5])
# plt.tight_layout()
plt.savefig(Path(plot_output_dir, 'scalar', sim +f'_{perturb_type}_output_{output_num}'  + '_rest_mass.pdf'), bbox_inches = "tight")

plotting maximum pressure vs. time


plt.figure(figsize=(8, 8))
adm_mass = combine_files('admmass-admmass_masses.maximum.asc', path_template, num_dir)
plt.plot(adm_mass[:, 1] / t_P, adm_mass[:, -1]
                    / M_ADM, label=sim)
plt.xlabel(r"$t / t_P$")
plt.ylabel(r"$M_{ADM} / M_{ADM, cocal}$")
plt.legend()
plt.ylim([0.8, 1.2])
# plt.tight_layout()
plt.savefig(Path(plot_output_dir, 'scalar', sim +f'_{perturb_type}_output_{output_num}'  + '_adm_mass.pdf'), bbox_inches="tight")


print("plotting stellar modes")
# this should be from sim_log to compare perturbed sim with other sims
plot_perturb = False
if not plot_perturb:
    r_modes = combine_files('stellarmodes-cm.maximum.asc', path_template, num_dir)
else:
    r_modes = combine_files('stellarmodes-cm.maximum.asc', path_template_perturb, num_dir_perturb)
modes = np.zeros_like(r_modes)
# imodes = np.genfromtxt(path_template + 'stellarmodes-cm.inorm2.asc')
for i in range(2, 13, 2):
    modes[:, i] = np.sqrt(r_modes[:, i]**2 + r_modes[:, (i+1)]**2)
modes[:, 1] = r_modes[:, 1]

# |Cm|
plt.figure(figsize=(8, 8))
plt.semilogy(modes[:, 1] / t_P, modes[:, 4] / modes[:, 2], '-.', label='m=1')
plt.semilogy(modes[:, 1] / t_P, modes[:, 6] / modes[:, 2], ':', label='m=2')
plt.semilogy(modes[:, 1] / t_P, modes[:, 8] / modes[:, 2], '--', label='m=3')
plt.semilogy(modes[:, 1] / t_P, modes[:, 10] / modes[:, 2], '-', label='m=4')
# plt.semilogy(modes[:, 1], modes[:, 12] / modes[:,2], label='m=5')
plt.legend()
plt.xlabel(r"$t / t_P$")
plt.ylabel(r"$|C_m| / C_0$")
plt.ylim([1E-5, 1])
# plt.tight_layout()
plt.savefig(Path(plot_output_dir, 'scalar', sim +f'_{perturb_type}_output_{output_num}'  + '_stellar_modes_magnitude.pdf'), bbox_inches="tight")

plotting stellar modes


# angular momentum
print("plotting angular momentum")
plt.figure(figsize=(8, 8))
angularmom = combine_files(
    'angularmom-jphi.maximum.asc', path_template, num_dir)
plt.plot(angularmom[:, 1] / t_P, angularmom[:, 2]
                            / J_cocal, label='Jphi int')
plt.legend()
plt.xlabel(r"$t / t_P$")
plt.ylabel(r'$J / J_{cocal}$')
# plt.ylim([0.8, 1.2])
# plt.tight_layout()
plt.savefig(Path(plot_output_dir, 'scalar', sim +f'_{perturb_type}_output_{output_num}'  + '_angular_mom.pdf'), bbox_inches="tight")

plotting angular momentum


# print('plotting the gw signal')
# # radi = [1.99, 3.09, 4.42, 5.53, 6.63]
# l = 2
# m = 2
# # m = 0, l
# plt.figure(figsize=(8, 8))
# # for r in radi:
# # r = 8.84 is the best res
# # r = r_gw
# for r in r_gw:
#     if not plot_perturb:
#         psi4 = combine_files('gw/' + "mp_Psi4_l"+str(l) +
#                             "_m"+str(m)+f"_r{r:.2f}.asc", path_template, num_dir)
#     else:
#         psi4 = combine_files('gw/' + "mp_Psi4_l"+str(l) +
#                         "_m"+str(m)+f"_r{r:.2f}.asc", path_template_perturb, num_dir_perturb)
#     plt.plot((psi4[:, 0] - r) / t_P, psi4[:, 1] * r, label=f'r={r}')
# plt.xlabel(r'$(t-r)/t_P$')
# plt.ylabel(r'$r \Psi_4 (l=2, m=2)$')
# plt.legend()
# # plt.tight_layout()
# plt.savefig(Path(plot_output_dir, sim +f'_{perturb_type}_output_{output_num}'  + '_gw.pdf', bbox_inches="tight")


Iij = combine_files(
    "quadmoment-iij.maximum.asc", path_template, num_dir)
com = combine_files(
     "quadmoment-com.maximum.asc", path_template, num_dir)


plt.figure(figsize=(10, 10))
t_dyn = 1 / np.sqrt(RHO0)
plt.plot(Iij[:, 1] / t_dyn, Iij[:, 2], label="Ixx")
plt.plot(Iij[:, 1]/ t_dyn, Iij[:, 3], label="Ixy")
plt.plot(Iij[:, 1]/ t_dyn, Iij[:, 4], label="Ixz")
plt.plot(Iij[:, 1]/ t_dyn, Iij[:, 5], label="Iyy")
plt.plot(Iij[:, 1]/ t_dyn, Iij[:, 6], label="Iyz")
plt.plot(Iij[:, 1]/ t_dyn, Iij[:, 7], label="Izz")
# plt.plot(Iij[:, 1], Iij[:, 2], label="Ixx")
plt.legend()

<matplotlib.legend.Legend at 0x7f9b09a97ac0>


times = cactus_hdf5_get_timeinfo(path_template.format(num_dir=num_dir, filename='illinoisgrmhd-grmhd_primitives_allbutbi.x.h5'))


output_step = 512
init_it = list(times.keys())[0]
final_it = list(times.keys())[-1]
init_it_step = int(init_it / output_step)
final_it_step = int(final_it / output_step)
output_time_steps = output_step * np.arange(init_it_step, final_it_step, int((final_it_step - init_it_step)/5.))
print(output_time_steps)
iterations = output_time_steps
iteration_init = iterations[0]
iteration_end = iterations[-1]

# iterations = [0]
# iteration_init = 0

[    0 19968 39936 59904 79872]


principle_axes = {}
for i, it in enumerate(iterations):
    Qmat = np.array([[Iij[i, 2], Iij[i, 3]], 
                    [Iij[i, 3], Iij[i, 5]]])
    eigenvalues, eigenvectors = np.linalg.eig(Qmat)
    idx = eigenvalues.argsort()[::-1]
    #     print(eigenvalues[idx], eigenvectors[:, idx])
    eigenvectors = eigenvectors[:, idx]
    principle_axes[it] = eigenvectors


print("plotting density x axis")
fig, ax = plt.subplots(1, 2, figsize=(20.0, 8.0))
interpolate_1d=False
for it in iterations:
    if interpolate_1d:
        rho_xy_1d, x1d, y1d = cactus_hdf5_get_level_data(
                path_template.format(num_dir=num_dir, filename='illinoisgrmhd-grmhd_primitives_allbutbi.xy.h5'), 'ILLINOISGRMHD::rho_b', it, 0)
        inter_rho_xy_plane = scipy.interpolate.interp2d(x1d, y1d, rho_xy_1d, kind='cubic')
        lin = np.linspace(-1.2 * Req, 1.2 * Req, 1000)
        inter_rho_x_axis = inter_rho_xy_plane(principle_axes[it][1, 0] * lin, principle_axes[it][1, 1] * lin)
        ax[0].plot(lin / Req, inter_rho_x_axis[int(inter_rho_x_axis.shape[0]/2), :] / RHO0,
                label=fr't={times[it]/t_P:.3g} $t_P$')
        inter_rho_y_axis = inter_rho_xy_plane(principle_axes[it][0, 0] * lin, principle_axes[it][0, 1] * lin)
        ax[1].plot(lin / Req, inter_rho_y_axis[:, int(inter_rho_y_axis.shape[1]/2)] / RHO0,
                label=fr't={times[it]/t_P:.3g} $t_P$')
    else:
        rho_x_1d = cactus_hdf5_get_1d_data(
            path_template.format(num_dir=num_dir, filename='illinoisgrmhd-grmhd_primitives_allbutbi.x.h5'),  'ILLINOISGRMHD::rho_b', it)
        rho_y_1d = cactus_hdf5_get_1d_data(
            path_template.format(num_dir=num_dir, filename='illinoisgrmhd-grmhd_primitives_allbutbi.y.h5'),  'ILLINOISGRMHD::rho_b', it)
        ax[0].plot(rho_x_1d[:, 0] / Req, rho_x_1d[:, 1] / RHO0,
                label=fr't={times[it]/t_P:.3g} $t_P$')
        ax[1].plot(rho_y_1d[:, 0] / Req, rho_y_1d[:, 1] / RHO0,
                label=fr't={times[it]/t_P:.3g} $t_P$')
ax[0].set_xlabel(r"$x /R_x$")
ax[0].set_ylabel(r"$\rho_x / \rho_{c,t=0}$")
ax[0].set_xlim([-1.2, 1.2])
ax[0].set_ylim([0, 1.2])
ax[0].legend(loc='upper right')
# plt.savefig(Path(plot_output_dir, sim+ f'_output_{num_dir}' +  '_rho_x' + '_plot.pdf', bbox_inches="tight")

ax[1].set_xlabel(r"$y /R_x$")
ax[1].set_ylabel(r"$\rho_y / \rho_{c,t=0}$")
ax[1].set_xlim([-1.2, 1.2])
ax[1].set_ylim([0, 1.2])
ax[1].legend(loc='upper right')
plt.savefig(Path(plot_output_dir, "1D", sim+ f'_output_{num_dir}' +  '_rho_x_y' + '_plot.pdf'), bbox_inches="tight")

plotting density x axis


path_template.format(num_dir=num_dir, filename='admbase-lapse.x.h5')

'/glade/scratch/yufengl/simulations/C01s17_switch_rxry_0419/output-0000/C01s17_vel_switch_rxry/admbase-lapse.x.h5'


print("plotting rotation curve across x axis")
plt.figure(figsize=(8., 8.))
for it in iterations:
    lapse_on_xaxis = cactus_hdf5_get_1d_data(
        path_template.format(num_dir=num_dir, filename='admbase-lapse.x.h5'), 'ADMBASE::alp', it)
    shift_y_on_xaxis = cactus_hdf5_get_1d_data(
        path_template.format(num_dir=num_dir, filename='admbase-shift.x.h5'), 'ADMBASE::betay', it)
    vel_y_on_xaxis = cactus_hdf5_get_1d_data(
        path_template.format(num_dir=num_dir, filename='hydrobase-vel.x.h5'), 'HYDROBASE::vel[1]', it)
    rest_vel_x = lapse_on_xaxis[:, 1] * vel_y_on_xaxis[:, 1] - shift_y_on_xaxis[:, 1]
    omega_x = rest_vel_x / vel_y_on_xaxis[:, 0]
    plt.plot(vel_y_on_xaxis[:, 0] / Req, omega_x / OMEGA_rot, label=fr't={times[it]/t_P:.3g} $t_P$')
plt.xlim([-1.2, 1.2])
plt.ylim([-1.5, 1.5])
plt.xlabel(r'$x / R_x$')
plt.ylabel(r'$\Omega_x / \Omega_{cocal}$')
plt.legend(loc='lower right')
plt.savefig(Path(plot_output_dir, "1D", sim+ f'_output_{num_dir}' +  '_omega_x_axis' + '_plot.pdf'), bbox_inches='tight')

plotting rotation curve across x axis


plt.figure(figsize=(8., 8.))
for it in iterations:
    vel_y_on_xaxis = cactus_hdf5_get_1d_data(
        path_template.format(num_dir=num_dir, filename='illinoisgrmhd-grmhd_primitives_allbutbi.x.h5'), 'ILLINOISGRMHD::vy', it)
    rest_vel_x = vel_y_on_xaxis[:, 1]
    omega_x = rest_vel_x / vel_y_on_xaxis[:, 0]
    plt.plot(vel_y_on_xaxis[:, 0] / Req, omega_x / OMEGA_rot, label=fr't={times[it]/t_P:.3g} $t_P$')
plt.xlim([-1.2, 1.2])
plt.ylim([-1.5, 1.5])
plt.xlabel(r'$x / R_x$')
plt.ylabel(r'$\Omega_x / \Omega_{cocal}$')
plt.legend(loc='lower right')
plt.savefig(Path(plot_output_dir, "1D", sim+ f'_output_{num_dir}' +  '_omega_x_axis_ILGRMHD' + '_plot.pdf'), bbox_inches='tight')


print("plotting lapse across x")
fig, ax = plt.subplots(1, 1, figsize=(18, 8.))
for it in iterations:
    lapse_on_xaxis = cactus_hdf5_get_1d_data(
        path_template.format(num_dir=num_dir, filename='admbase-lapse.x.h5'), 'ADMBASE::alp', it)
    ax.plot(lapse_on_xaxis[:, 0] / Req, lapse_on_xaxis[:, 1], label=fr't={times[it]/t_P:.3g} $t_P$')
    ax.set_xlabel(r'$x / R_x$')
    ax.set_ylabel(r'$\alpha$')
    ax.set_xlim([-1.2, 1.2])
    ax.legend(loc='upper right')
    # lapse_on_yaxis = cactus_hdf5_get_1d_data(
    #     f'{root_dir}/output-000{num_dir}/admbase-lapse.y.h5', 'ADMBASE::alp', it)
    # ax[1].plot(lapse_on_yaxis[:, 0] / Req, lapse_on_yaxis[:, 1])
    # ax[1].xlabel(r'$y / R_x$')
    # ax[1].ylabel(r'\alpha')
    # ax[1].set_xlim([-1.2 * Req, 1.2 * Req])
    # ax[1].legend(loc='upper right')
plt.show()
plt.savefig(Path(plot_output_dir, "1D", sim+ f'_output_{num_dir}' +  '_lapse_x' + '_line_plot.pdf'), bbox_inches='tight')
# plt.close()

plotting lapse across x

<Figure size 432x288 with 0 Axes>


xy_2d_path = path_template.format(num_dir=num_dir, filename='illinoisgrmhd-grmhd_primitives_allbutbi.xy.h5')
rho_in_xy_plane_init = cactus_hdf5_get_2d_data(xy_2d_path, 'ILLINOISGRMHD::rho_b', iteration_init)
rho_in_xy_plane_final = cactus_hdf5_get_2d_data(xy_2d_path, 'ILLINOISGRMHD::rho_b', iteration_end)


print(xy_2d_path)

/glade/scratch/yufengl/simulations/C01s17_switch_rxry_0419/output-0000/C01s17_vel_switch_rxry/illinoisgrmhd-grmhd_primitives_allbutbi.xy.h5


plt.figure(figsize=(16, 8))
lin = np.linspace(-Req * 2, Req * 2, 100)
plt.subplot(1, 2, 1)
plt.scatter(rho_in_xy_plane_init[:, 0], rho_in_xy_plane_init[:, 1], c=np.log(rho_in_xy_plane_init[:, 2]/RHO0))
plt.plot(principle_axes[iteration_init][1, 0] * lin, principle_axes[iteration_init][1, 1] * lin, label='x axis')
plt.plot(principle_axes[iteration_init][0, 0] * lin, principle_axes[iteration_init][0, 1] * lin, label='y axis')
plt.legend()
plt.xlim([-2 * Req, 2 * Req])
plt.ylim([-2 * Req, 2 * Req])

plt.subplot(1, 2, 2)
plt.scatter(rho_in_xy_plane_final[:, 0], rho_in_xy_plane_final[:, 1], c=np.log(rho_in_xy_plane_final[:, 2]/RHO0))
plt.plot(principle_axes[iteration_end][1, 0] * lin, principle_axes[iteration_end][1, 1] * lin, label='x axis')
plt.plot(principle_axes[iteration_end][0, 0] * lin, principle_axes[iteration_end][0, 1] * lin, label='y axis')
plt.legend()
plt.xlim([-2 * Req, 2 * Req])
plt.ylim([-2 * Req, 2 * Req])

(-0.8842785408075154, 0.8842785408075154)


print("plotting density xy")
# plt.figure(figsize=(24, 10))
# plt.subplot(1, 2, 1)
fig, ax = plt.subplots(1, 2, figsize=(24, 10))
cactus_hdf5_plot2d_pseudocolor(
    xy_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_init, -1.5 * Req, 1.5*Req, -1.5 * Req, 1.5*Req, const_normalize=RHO0, 
    fig_handle=fig, ax_handle=ax[0])

plt.subplot(1, 2, 2)
cactus_hdf5_plot2d_pseudocolor(
    xy_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_end, -1.5 * Req, 1.5*Req, -1.5 * Req, 1.5*Req, const_normalize=RHO0,
    fig_handle=fig, ax_handle=ax[1])
plt.savefig(Path(plot_output_dir, "2D", sim+ f'_output_{num_dir}' +  '_rho_xy' + '_pseudocolor_plot.pdf'), bbox_inches="tight")
plt.savefig(Path(plot_output_dir, "2D", sim+ f'_output_{num_dir}' +  '_rho_xy' + '_pseudocolor_plot.png'), bbox_inches="tight")

plotting density xy


x0 = rho_in_xy_plane_init[:, 0]
y0 = rho_in_xy_plane_init[:, 1]
rho0 = rho_in_xy_plane_init[:, 2]

xf = rho_in_xy_plane_final[:, 0]
yf = rho_in_xy_plane_final[:, 1]
rhof = rho_in_xy_plane_final[:, 2]

N = 3000
x0i, y0i, rho0i = interpol(x0, y0, rho0, N)
xfi, yfi, rhofi = interpol(xf, yf, rhof, N)


fig = plt.figure(figsize=(10.0,8.0))
ax1 = fig.add_subplot(1,1,1)
ax1.grid(color='black', linestyle=':', linewidth=0.5)
levels = [0.4, 0.5, 0.7, 0.9, 0.95]
levcol = [ 'green', 'cyan', 'brown', 'orange', 'magenta' ]
contour1 = ax1.contour(x0i/Req, y0i/Req, rho0i/RHO0, levels=levels, colors=levcol, linestyles='dotted')
contour2 = ax1.contour(xfi/Req, yfi/Req, rhofi/RHO0, levels=levels, colors=levcol, linestyles='-')
# ax1.clabel(contour1, levels, inline=True, inline_spacing=0, fontsize=14)
# ax1.clabel(contour2, levels, inline=True, inline_spacing=1, fontsize=14)
ax1.set_xlim([-1, 1])
ax1.set_ylim([-1, 1])
# ax1.legend()

(-1.0, 1.0)


xz_2d_path = path_template.format(num_dir=num_dir, filename='illinoisgrmhd-grmhd_primitives_allbutbi.xz.h5')
rho_in_xz_plane_init = cactus_hdf5_get_2d_data(
    xz_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_init)

rho_in_xz_plane_final = cactus_hdf5_get_2d_data(
    xz_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_end)
# fig, ax = plt.subplots(1, 2, figsize=(25.0, 8.0))
# dens_plot = ax[0].scatter(rho_in_xz_plane_init[:, 0] / Req, rho_in_xz_plane_init[:, 1] / Req,
#                             c=numpy.log(rho_in_xz_plane_init[:, 2] / RHO0), 
#                             label=fr"$\rho_0/\rho_0(c,t=0)$ at time={(times[iteration_init] / t_P):.3g} t_P",
#                             vmin=-4)
# ax[0].set_aspect('equal')
# ax[0].set_xlim(-1.5, 1.5)
# ax[0].set_ylim(-1.5, 1.5)
# ax[0].set_title("Scatter plot using xz data")
# ax[0].set_xlabel("x [M]")
# ax[0].set_ylabel("z [M]")
# ax[0].legend(loc='upper right')
# fig.colorbar(dens_plot, ax=ax[0])

# dens_plot = ax[1].scatter(rho_in_xz_plane_final[:, 0] / Req, rho_in_xz_plane_final[:, 1] / Req,
#                             c=numpy.log(rho_in_xz_plane_final[:, 2] / RHO0), 
#                             label=fr"$\rho_0/\rho_0(c,t=0)$ at time={(times[iteration_end] / t_P):.3g} t_P",
#                             vmin=-4)
# ax[1].set_aspect('equal')
# ax[1].set_xlim(-1.5, 1.5)
# ax[1].set_ylim(-1.5, 1.5)
# ax[1].set_title("Scatter plot using xz data")
# ax[1].set_xlabel("x [M]")
# ax[1].set_ylabel("z [M]")
# ax[1].legend(loc='upper right')
# fig.colorbar(dens_plot, ax=ax[1])
# plt.savefig(Path(plot_output_dir, sim+ f'_output_{num_dir}' +  '_rho_xz' + '_scatter_plot.pdf', bbox_inches="tight")
# plt.savefig(Path(plot_output_dir, sim+ f'_output_{num_dir}' +  '_rho_xz' + '_scatter_plot.png', bbox_inches="tight")


plt.figure(figsize=(24, 10))
fig, ax = plt.subplots(1, 2, figsize=(24, 10))
cactus_hdf5_plot2d_pseudocolor(
    xz_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_init, -1.5 * Req, 1.5*Req, -1.5 * Req, 1.5*Req, const_normalize=RHO0, 
    fig_handle=fig, ax_handle=ax[0])

# plt.subplot(1, 2, 2)
cactus_hdf5_plot2d_pseudocolor(
    xz_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_end, -1.5 * Req, 1.5*Req, -1.5 * Req, 1.5*Req, const_normalize=RHO0,
    fig_handle=fig, ax_handle=ax[1])
plt.savefig(Path(plot_output_dir, "2D", sim+ f'_output_{num_dir}' +  '_rho_xz' + '_pseudocolor_plot.pdf'), bbox_inches="tight")
plt.savefig(Path(plot_output_dir, "2D", sim+ f'_output_{num_dir}' +  '_rho_xz' + '_pseudocolor_plot.png'), bbox_inches="tight")

<Figure size 1728x720 with 0 Axes>


# xz plane 
x20 = rho_in_xz_plane_init[:, 0]
z20 = rho_in_xz_plane_init[:, 1]
rho20 = rho_in_xz_plane_init[:, 2]

x2f = rho_in_xz_plane_final[:, 0]
z2f = rho_in_xz_plane_final[:, 1]
rho2f = rho_in_xz_plane_final[:, 2]

N = 3000
x20i, z20i, rho20i = interpol(x20, z20, rho20, N)
x2fi, z2fi, rho2fi = interpol(x2f, z2f, rho2f, N)

# fig = plt.figure(figsize=(10.0,8.0))
# ax1 = fig.add_subplot(1,1,1)
# ax1.grid(color='black', linestyle=':', linewidth=0.5)
# levels = [0.4, 0.5, 0.7, 0.9, 0.95]
# levcol = [ 'green', 'cyan', 'brown', 'orange', 'magenta' ]
# contour1 = ax1.contour(x20i/Req, z20i/Req, rho20i/RHO0, levels=levels, colors=levcol, linestyles='dotted')
# contour2 = ax1.contour(x2fi/Req, z2fi/Req, rho2fi/RHO0, levels=levels, colors=levcol, linestyles='-')
# ax1.set_xlim([-1, 1])
# ax1.set_ylim([-1, 1])


# get quiver plot using unstructured interpolation of all data (slow)
print("plotting velocity xy")
fig, ax = plt.subplots(1, 2, figsize=(20.0, 8.0))
velx_in_xy_plane_init = cactus_hdf5_get_2d_data(
    xy_2d_path, 
    'ILLINOISGRMHD::vx', iteration_init)
vely_in_xy_plane_init = cactus_hdf5_get_2d_data(
    xy_2d_path, 
    'ILLINOISGRMHD::vy', iteration_init)
# get a uniform grid
x1d = numpy.linspace(-1.2, 1.2, 20, endpoint=True)
y1d = numpy.linspace(-1.2, 1.2, 20, endpoint=True)
x, y = numpy.meshgrid(x1d, y1d, indexing='xy')
# interpolate onto grid
u = scipy.interpolate.griddata(
    (velx_in_xy_plane_init[:, 0]/Req, velx_in_xy_plane_init[:, 1]/Req), velx_in_xy_plane_init[:, 2], (x, y))
v = scipy.interpolate.griddata(
    (vely_in_xy_plane_init[:, 0]/Req, vely_in_xy_plane_init[:, 1]/Req), vely_in_xy_plane_init[:, 2], (x, y))
# ax[0].scatter(rho_in_xy_plane_init[:, 0] / Req, rho_in_xy_plane_init[:, 1] / Req,
#                             c=numpy.log(rho_in_xy_plane_init[:, 2] / RHO0), 
#                             label=fr"$\rho_0/\rho_0(c,t=0)$ at time={times[iteration_init] / t_P:.3g} t_P",
#                             vmin=-4)
cactus_hdf5_plot2d_pseudocolor(
    xy_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_init, -1.2, 1.2, -1.2, 1.2, const_normalize=RHO0,
    xy_normalize=Req, ax_handle=ax[0], fig_handle=fig, plot_reflevel=False)

ax[0].quiver(x, y, u, v, label="t=%.3g" % (times[iteration_init] / t_P), angles='xy',scale=6,width=0.005,color='black', zorder=100)
ax[0].set_aspect('equal')
# ax[0].set_xlim(-1.2, 1.2)
# ax[0].set_ylim(-1.2, 1.2)
ax[0].set_title("3-velocity in xy plane")
ax[0].set_xlabel("x [M]")
ax[0].set_ylabel("y [M]")
ax[0].legend(loc='upper right')

velx_in_xy_plane_final = cactus_hdf5_get_2d_data(
    xy_2d_path, 
    'ILLINOISGRMHD::vx', iteration_end)
vely_in_xy_plane_final = cactus_hdf5_get_2d_data(
    xy_2d_path, 
    'ILLINOISGRMHD::vy', iteration_end)

# get a uniform grid
x1d = numpy.linspace(-1.2, 1.2, 20, endpoint=True)
y1d = numpy.linspace(-1.2, 1.2, 20, endpoint=True)
x, y = numpy.meshgrid(x1d, y1d, indexing='xy')
# interpolate onto grid
u = scipy.interpolate.griddata(
    (velx_in_xy_plane_final[:, 0]/Req, velx_in_xy_plane_final[:, 1]/Req), velx_in_xy_plane_final[:, 2], (x, y))
v = scipy.interpolate.griddata(
    (vely_in_xy_plane_final[:, 0]/Req, vely_in_xy_plane_final[:, 1]/Req), vely_in_xy_plane_final[:, 2], (x, y))
# ax[1].scatter(rho_in_xy_plane_final[:, 0] / Req, rho_in_xy_plane_final[:, 1] / Req,
#                             c=numpy.log(rho_in_xy_plane_final[:, 2] / RHO0), 
#                             label=fr"$\rho_0/\rho_0(c,t=0)$ at time={times[iteration_end] / t_P:.3g} t_P",
#                             vmin=-4)
cactus_hdf5_plot2d_pseudocolor(
    xy_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_end, -1.2, 1.2, -1.2, 1.2, const_normalize=RHO0,
    xy_normalize=Req, ax_handle=ax[1], fig_handle=fig, plot_reflevel=False)
ax[1].quiver(x, y, u, v, label="t=%.3g" % (times[iteration_end] / t_P), angles='xy',scale=6,width=0.005,color='black', zorder=100)
ax[1].set_aspect('equal')
# ax[1].set_xlim(-1.2, 1.2)
# ax[1].set_ylim(-1.2, 1.2)
ax[1].set_title("3-velocity in xy plane")
ax[1].set_xlabel("x [M]")
ax[1].set_ylabel("y [M]")
ax[1].legend(loc='upper right')

plt.savefig(Path(plot_output_dir, "2D", sim+ f'_output_{num_dir}' +  '_vel_xy' + '_quiver_plot.pdf'), bbox_inches="tight")
plt.savefig(Path(plot_output_dir, "2D", sim+ f'_output_{num_dir}' +  '_vel_xy' + '_quiver_plot.png'), bbox_inches="tight")

plotting velocity xy


# sample_velxi = idw.tree(velx_in_xy_plane_init[:, :2], velx_in_xy_plane_init[:, 2])
# sample_velyi = idw.tree(vely_in_xy_plane_init[:, :2], vely_in_xy_plane_init[:, 2])

# sample_velxf = idw.tree(velx_in_xy_plane_final[:, :2], velx_in_xy_plane_final[:, 2])
# sample_velyf = idw.tree(vely_in_xy_plane_final[:, :2], vely_in_xy_plane_final[:, 2])

# fig = plt.figure(figsize=(10.0, 10.0))
# ax1 = fig.add_subplot(1,1,1)
# ax1.grid(color='black', linestyle=':', linewidth=0.5)
# levels = [0.4, 0.5, 0.7, 0.9, 0.95]
# levcol = [ 'green', 'cyan', 'brown', 'orange', 'magenta' ]
# contour1 = ax1.contour(x0i/Req, y0i/Req, rho0i/RHO0, levels=levels, colors=levcol, linestyles='dotted', label=fr"$\rho_0$ at time={times[iteration_init] / t_P:.3g} t_P")
# contour2 = ax1.contour(xfi/Req, yfi/Req, rhofi/RHO0, levels=levels, colors=levcol, linestyles='-', label=fr"$\rho_0$ at time={times[iteration_end] / t_P:.3g} t_P")

# step = 20
# for cont in contour1.allsegs:
#     coord = cont[0]
#     velxi = sample_velxi(coord)
#     velyi = sample_velyi(coord)
#     ax1.quiver(coord[::step, 0], coord[::step, 1], velxi[::step], velyi[::step], angles='xy',scale=3,width=0.005,color='red')

# for cont in contour2.allsegs:
#     coord = cont[0]
#     velxf = sample_velxf(coord)
#     velyf = sample_velyf(coord)
#     ax1.quiver(coord[::step, 0], coord[::step, 1], velxf[::step], velyf[::step], angles='xy',scale=3,width=0.005,color='blue')

# ax1.set_xlim([-1, 1])
# ax1.set_ylim([-1, 1])
# ax1.set_aspect("equal")


x0 = rho_in_xy_plane_init[:, 0]
y0 = rho_in_xy_plane_init[:, 1]
rho0 = rho_in_xy_plane_init[:, 2]
N = 3000
x0i, y0i, rho0i = interpol(x0, y0, rho0, N)
# xfi, yfi, rhofi = interpol(xf, yf, rhof, N)

fig = plt.figure(figsize=(10.0, 10.0))
ax1 = fig.add_subplot(1,1,1)
ax1.grid(color='black', linestyle=':', linewidth=0.5)
levels = [0.5, 0.7, 0.8, 0.9, 0.95]
levcol = [ 'green', 'cyan', 'brown', 'orange', 'magenta' ]
contour1 = ax1.contour(x0i/Req, y0i/Req, rho0i/RHO0, levels=levels, colors=levcol, linestyles='dotted', label=fr"$\rho_0$ at time={times[iteration_init] / t_P:.3g} t_P")
# levels = [0.85, 0.87, 0.9]
# levcol = [ 'green', 'cyan', 'brown']
# contour1 = ax1.contour(x0i/Req, y0i/Req, alp0i, levels=levels, colors=levcol, linestyles='-')
# contour2 = ax1.contour(xfi/Req, yfi/Req, rhofi/RHO0, levels=levels, colors=levcol, linestyles='-', label=fr"$\rho_0$ at time={times[iteration_end] / t_P:.3g} t_P")
rho_mask = [((0.94 < rho0 / RHO0) & (rho0 / RHO0 < 0.96)), ((0.88 < rho0 / RHO0) & (rho0 / RHO0 < 0.92)),
            ((0.78 < rho0 / RHO0) & (rho0 / RHO0 < 0.82)), ((0.68 < rho0 / RHO0) & (rho0 / RHO0 < 0.72)),
           ((0.48 < rho0 / RHO0) & (rho0 / RHO0 < 0.52))]
steps = 20
for i in range(len(rho_mask)):
    ax1.quiver(velx_in_xy_plane_init[rho_mask[i], 0][::steps]/Req, velx_in_xy_plane_init[rho_mask[i], 1][::steps]/Req, 
               velx_in_xy_plane_init[rho_mask[i], 2][::steps], vely_in_xy_plane_init[rho_mask[i], 2][::steps], angles='xy',scale=4,width=0.005,color='red')
    
#     ax1.quiver(velx_in_xy_plane_init_non[rho_mask[i], 0][::steps]/Req, velx_in_xy_plane_init_non[rho_mask[i], 1][::steps]/Req, 
#                velx_in_xy_plane_init_non[rho_mask[i], 2][::steps], vely_in_xy_plane_init_non[rho_mask[i], 2][::steps], angles='xy',scale=4,width=0.005,color='green')
    
#     ax1.quiver(velx_in_xy_plane_init_vel_cor[rho_mask[i], 0][::steps]/Req, velx_in_xy_plane_init_vel_cor[rho_mask[i], 1][::steps]/Req, 
#                velx_in_xy_plane_init_vel_cor[rho_mask[i], 2][::steps], vely_in_xy_plane_init_vel_cor[rho_mask[i], 2][::steps], angles='xy',scale=4,width=0.005,color='blue')
#     tan_vec = find_tan_using_ellipse(velx_in_xy_plane_init[rho_mask[i], 0]/Req, velx_in_xy_plane_init[rho_mask[i], 1]/Req,
#                                         Rx, Ry)
#     ax1.quiver(velx_in_xy_plane_init[rho_mask[i], 0]/Req, velx_in_xy_plane_init[rho_mask[i], 1]/Req,
#                tan_vec[:, 0], tan_vec[:, 1], angles='xy',scale=4,width=0.005,color='black')
# step = 20
# for cont in contour1.allsegs:
#     coord = cont[0]
#     velxi = sample_velxi(coord)
#     velyi = sample_velyi(coord)
#     ax1.quiver(coord[::step, 0], coord[::step, 1], velxi[::step], velyi[::step], angles='xy',scale=2,width=0.005,color='red')
#     tan_vec = np.zeros_like(coord)
#     for i in range(len(coord)):
#         # if i - 1 < 0:
#         #     p1 = coord[-1]
#         # else:
#         #     p1 = coord[i-1]
#         # p2 = coord[i]
#         # if i + 1 == len(coord):
#         #     p3 = coord[0]
#         # else:
#         #     p3 = coord[i+1]
#         # tan_vec[i] = find_unit_tangential_vec(p1, p2, p3)
#         tan_vec[i] = find_tan_using_ellipse(coord[i, 0], coord[i, 1], Rx, Ry)
#     ax1.quiver(coord[::step, 0], coord[::step, 1], tan_vec[::step, 0], tan_vec[::step, 1], angles='xy',scale=15,width=0.005,color='black')
ax1.set_xlim([-1, 1])
ax1.set_ylim([-1, 1])
ax1.set_aspect("equal")

/glade/scratch/yufengl/ipykernel_63186/2543727547.py:13: UserWarning: The following kwargs were not used by contour: 'label'
  contour1 = ax1.contour(x0i/Req, y0i/Req, rho0i/RHO0, levels=levels, colors=levcol, linestyles='dotted', label=fr"$\rho_0$ at time={times[iteration_init] / t_P:.3g} t_P")


print("plotting velocity xz")
fig, ax = plt.subplots(1, 2, figsize=(20.0, 8.0))
velx_in_xz_plane_init = cactus_hdf5_get_2d_data(
    xz_2d_path, 
    'ILLINOISGRMHD::vx', iteration_init)
velz_in_xz_plane_init = cactus_hdf5_get_2d_data(
    xz_2d_path, 
    'ILLINOISGRMHD::vz', iteration_init)

# get a uniform grid
x1d = numpy.linspace(-1.2, 1.2, 20, endpoint=True)
z1d = numpy.linspace(-1.2, 1.2, 20, endpoint=True)
x, z = numpy.meshgrid(x1d, z1d, indexing='xy')

# interpolate onto grid
u = scipy.interpolate.griddata(
    (velx_in_xz_plane_init[:, 0]/Req, velx_in_xz_plane_init[:, 1]/Req), velx_in_xz_plane_init[:, 2], (x, z))
v = scipy.interpolate.griddata(
    (velz_in_xz_plane_init[:, 0]/Req, velz_in_xz_plane_init[:, 1]/Req), velz_in_xz_plane_init[:, 2], (x, z))
# ax[0].scatter(rho_in_xz_plane_init[:, 0] / Req, rho_in_xz_plane_init[:, 1] / Req,
#                             c=numpy.log(rho_in_xz_plane_init[:, 2] / RHO0), 
#                             label=fr"$\rho_0/\rho_0(c,t=0)$ at time={(times[iteration_init] / t_P):.3g} t_P",
#                             vmin=-4)
cactus_hdf5_plot2d_pseudocolor(
    xz_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_init, -1.2, 1.2, -1.2, 1.2, const_normalize=RHO0,
    xy_normalize=Req, ax_handle=ax[0], fig_handle=fig, plot_reflevel=False)
ax[0].quiver(x, z, u, v, label="t=%.3g" % (times[iteration_init] / t_P), angles='xy',scale=6,width=0.005,color='black', zorder=100)
ax[0].set_aspect('equal')
ax[0].set_xlim(-1.2, 1.2)
ax[0].set_ylim(-1.2, 1.2)
ax[0].set_title("3-velocity in xz plane")
ax[0].set_xlabel("x [M]")
ax[0].set_ylabel("z [M]")
ax[0].legend(loc='upper right')
# http://einsteintoolkit.org/thornguide/EinsteinBase/HydroBase/documentation.html (1) -> v^i
velx_in_xz_plane_final = cactus_hdf5_get_2d_data(
    xz_2d_path, 
    'ILLINOISGRMHD::vx', iteration_end)
velz_in_xz_plane_final = cactus_hdf5_get_2d_data(
    xz_2d_path, 
    'ILLINOISGRMHD::vz', iteration_end)

# get a uniform grid
x1d = numpy.linspace(-1.2, 1.2, 20, endpoint=True)
z1d = numpy.linspace(-1.2, 1.2, 20, endpoint=True)
x, z = numpy.meshgrid(x1d, z1d, indexing='xy')

# interpolate onto grid
u = scipy.interpolate.griddata(
    (velx_in_xz_plane_final[:, 0]/Req, velx_in_xz_plane_final[:, 1]/Req), velx_in_xz_plane_final[:, 2], (x, z))
v = scipy.interpolate.griddata(
    (velz_in_xz_plane_final[:, 0]/Req, velz_in_xz_plane_final[:, 1]/Req), velz_in_xz_plane_final[:, 2], (x, z))
# ax[1].scatter(rho_in_xz_plane_final[:, 0] / Req, rho_in_xz_plane_final[:, 1] / Req,
#                             c=numpy.log(rho_in_xz_plane_final[:, 2] / RHO0), 
#                             label=fr"$\rho_0/\rho_0(c,t=0)$ at time={(times[iteration_end] / t_P):.3g} t_P",
#                             vmin=-4)
cactus_hdf5_plot2d_pseudocolor(
    xz_2d_path, 
    'ILLINOISGRMHD::rho_b', iteration_end, -1.2, 1.2, -1.2, 1.2, const_normalize=RHO0,
    xy_normalize=Req, ax_handle=ax[1], fig_handle=fig, plot_reflevel=False)
ax[1].quiver(x, z, u, v, label="t=%.3g" % (times[iteration_end] / t_P), angles='xy',scale=6,width=0.005,color='black', zorder=100)
ax[1].set_aspect('equal')
ax[1].set_xlim(-1.2, 1.2)
ax[1].set_ylim(-1.2, 1.2)
ax[1].set_title("3-velocity in xz plane")
ax[1].set_xlabel("x [M]")
ax[1].set_ylabel("z [M]")
ax[1].legend(loc='upper right')
# plt.show()
plt.savefig(Path(plot_output_dir, "2D", sim+ f'_output_{num_dir}' +  '_vel_xz' + '_quiver_plot.pdf'), bbox_inches="tight")
plt.savefig(Path(plot_output_dir, "2D", sim+ f'_output_{num_dir}' +  '_vel_xz' + '_quiver_plot.png'), bbox_inches="tight")

plotting velocity xz


# sample_velxi = idw.tree(velx_in_xz_plane_init[:, :2], velx_in_xz_plane_init[:, 2])
# sample_velzi = idw.tree(velz_in_xz_plane_init[:, :2], velz_in_xz_plane_init[:, 2])

# sample_velxf = idw.tree(velx_in_xz_plane_final[:, :2], velx_in_xz_plane_final[:, 2])
# sample_velzf = idw.tree(velz_in_xz_plane_final[:, :2], velz_in_xz_plane_final[:, 2])

# fig = plt.figure(figsize=(10.0, 10.0))
# ax1 = fig.add_subplot(1,1,1)
# ax1.grid(color='black', linestyle=':', linewidth=0.5)
# levels = [0.4, 0.5, 0.7, 0.9, 0.95]
# levcol = [ 'green', 'cyan', 'brown', 'orange', 'magenta' ]
# contour1 = ax1.contour(x20i/Req, z20i/Req, rho20i/RHO0, levels=levels, colors=levcol, linestyles='dotted')
# contour2 = ax1.contour(x2fi/Req, z2fi/Req, rho2fi/RHO0, levels=levels, colors=levcol, linestyles='-')

# for cont in contour1.allsegs:
#     coord = cont[0]
#     velxi = sample_velxi(coord)
#     velzi = sample_velzi(coord)
#     ax1.quiver(coord[:, 0], coord[:, 1], velxi, velzi, angles='xy',scale=3,width=0.005,color='red')

# for cont in contour2.allsegs:
#     coord = cont[0]
#     velxf = sample_velxf(coord)
#     velzf = sample_velzf(coord)
#     ax1.quiver(coord[:, 0], coord[:, 1], velxf, velzf, angles='xy',scale=3,width=0.005,color='blue')

# ax1.set_xlim([-1, 1])
# ax1.set_ylim([-1, 1])
# ax1.set_aspect("equal")

Scalar diagnostics¶

Central density¶

Rest mass and ADM mass¶

Stellar modes¶

Angular momentum $J_\phi$¶

GW $\psi_4$¶

QuadMoments¶

1D output¶

density change along x-axis¶

$\Omega_x$ Rotation curve on x-axis¶

Lapse on x axis $\alpha_x$¶

2D output¶

Density in xy plane $\rho_{xy}$¶

Density in xz plane $\rho_{xz}$¶

Velocity on xy plane $v_{xy}$¶

Velocity on xz plane $v_{xz}$¶