Skip to content

Instantly share code, notes, and snippets.

@dutta-alankar

dutta-alankar/cloud-props_paraview-analysis.py

Last active November 14, 2024 13:46
Show Gist options
  • Save dutta-alankar/6d2484212a4ca19e2814e19c5f7102ec to your computer and use it in GitHub Desktop.
Save dutta-alankar/6d2484212a4ca19e2814e19c5f7102ec to your computer and use it in GitHub Desktop.
Download ZIP
cloud-props_paraview-analysis
Raw
cloud-props_paraview-analysis.py
# -*- coding: utf-8 -*-
"""
Created on Mon Oct 28 17:39:48 2024
@author: alankar.
Usage: time python area-analysis.py
"""
import numpy as np
import scipy
import subprocess as sp
from scipy.interpolate import interp1d
import sys
import os
import pickle
import matplotlib
import matplotlib.pyplot as plt
dark = False
## Plot Styling
matplotlib.rcParams["xtick.direction"] = "in"
matplotlib.rcParams["ytick.direction"] = "in"
matplotlib.rcParams["xtick.top"] = False
matplotlib.rcParams["ytick.right"] = True
matplotlib.rcParams["xtick.minor.visible"] = True
matplotlib.rcParams["ytick.minor.visible"] = True
matplotlib.rcParams["axes.grid"] = True
matplotlib.rcParams["grid.linestyle"] = ":"
matplotlib.rcParams["grid.linewidth"] = 0.8
matplotlib.rcParams["grid.color"] = "gray"
matplotlib.rcParams["grid.alpha"] = 0.3
matplotlib.rcParams["lines.dash_capstyle"] = "round"
matplotlib.rcParams["lines.solid_capstyle"] = "round"
matplotlib.rcParams["legend.handletextpad"] = 0.4
matplotlib.rcParams["axes.linewidth"] = 1.0
matplotlib.rcParams["lines.linewidth"] = 3.5
matplotlib.rcParams["ytick.major.width"] = 1.2
matplotlib.rcParams["xtick.major.width"] = 1.2
matplotlib.rcParams["ytick.minor.width"] = 1.0
matplotlib.rcParams["xtick.minor.width"] = 1.0
matplotlib.rcParams["ytick.major.size"] = 11.0
matplotlib.rcParams["xtick.major.size"] = 11.0
matplotlib.rcParams["ytick.minor.size"] = 5.0
matplotlib.rcParams["xtick.minor.size"] = 5.0
matplotlib.rcParams["xtick.major.pad"] = 10.0
matplotlib.rcParams["xtick.minor.pad"] = 10.0
matplotlib.rcParams["ytick.major.pad"] = 6.0
matplotlib.rcParams["ytick.minor.pad"] = 6.0
matplotlib.rcParams["xtick.labelsize"] = 26.0
matplotlib.rcParams["ytick.labelsize"] = 26.0
matplotlib.rcParams["axes.titlesize"] = 24.0
matplotlib.rcParams["axes.labelsize"] = 28.0
matplotlib.rcParams["axes.labelpad"] = 8.0
plt.rcParams["font.size"] = 28
matplotlib.rcParams["legend.handlelength"] = 2
# matplotlib.rcParams["figure.dpi"] = 200
matplotlib.rcParams["axes.axisbelow"] = True
matplotlib.rcParams["figure.figsize"] = (13,10)
if dark:
plt.style.use('dark_background')
tcool_mix_B_tcc = [0.08, 0.10, 0.20, 0.50, 0.80, 1.00, 1.40, 2.50, 8.00,]
rini_B_rcl = [28.268, 35.335, 70.671, 176.677, 282.684, 353.355, 494.697, 883.387, 2826.838,]
tcool_mix_B_tcc = tcool_mix_B_tcc[1:]
rini_B_rcl = rini_B_rcl[1:]
mach = 1.496
Tcl = 4.0e+04
chi = 100
file_ext = "flt.h5"
tcc = np.sqrt(chi)
gamma = 5/3.
till = 100
cs = 1.0/mach # in unit velocity
wind = np.loadtxt('../../CC85_steady-prof_gamma_1.667.txt', skiprows=1)
mach_data = wind[:,3]/np.sqrt(gamma*wind[:,2]/wind[:,1])
rnorm = wind[:,0]
relpos = interp1d(mach_data, rnorm) #inverting the Mach relation
diniBdinj = relpos(mach)
CC85windrho = interp1d(rnorm, wind[:,1])
CC85windprs = interp1d(rnorm, wind[:,2])
CC85windvel = interp1d(rnorm, wind[:,3])
mu = 0.60917
mp = 1.6726e-24
kB = 1.3807e-16
Myr = 1.0e+06 * 365*24*60*60
MSun = 1.99e+33
Tmix = np.sqrt(chi)*Tcl
Tcutoff = 9.0e+04
Temperature_identify = 2.0*Tcl
print("Tmix: %.2e K"%Tmix)
print("Tcutoff: %.2e K"%Tcutoff)
print("Temperature_identify: %.2e K"%Temperature_identify)
root = "../../output"
till = 100 # tcc
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 100
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# t/tcc Area
area_data = []
for key in list(cloud_data.keys()):
area_data.append([float(key), cloud_data[key]['cloud_tot_surface_area']])
area_data = np.array(area_data)
line, = plt.plot(area_data[:till+1,0], area_data[:till+1,1]/area_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,0], (area_data[:till+1,2]/area_data[0,2])**0.845, linestyle ='--', color = line.get_color())
plt.legend(loc="upper left", title=r"$t_{\rm cool, mix}/t_{\rm cc}|_{\rm ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
plt.ylim(ymin=5.0e-02)
plt.xlim(xmin=0., xmax = till)
plt.yscale('log')
plt.xlabel(r"time [$t_{\rm cc,ini}$]")
plt.ylabel(r"Surface area $A_{\rm cold}$ ($T<8\times 10^4$ K) [$A_{\rm cold, ini}$]")
plt.savefig(f"area-cold-trunc_time{'-dark' if dark else ''}.svg", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
ndens_w_ini = UNIT_DENSITY/(mu*mp) # cgs
ndens_cl_ini = chi * ndens_w_ini
prs_w_ini = ndens_w_ini*chi*Tcl # p/kB
v_w_ini = UNIT_VELOCITY # cgs
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 100
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
mass_data = []
for key in list(cloud_data.keys()):
cloud_distance = cloud_data[key]['cloud_distance']/distance_ini
ndens_w = (CC85windrho(cloud_distance * diniBdinj)/CC85windrho(diniBdinj)) * UNIT_DENSITY/(mu*mp) # cgs
cloud_ndens = cloud_data[key]['cloud_density'] * UNIT_DENSITY/(mu*mp) # cgs
cloud_volume = cloud_data[key]['cloud_volume_elems'] * UNIT_LENGTH**3 # cgs
cloud_com = np.sum(cloud_ndens*cloud_volume*cloud_distance)/np.sum(cloud_ndens*cloud_volume)
condition = cloud_ndens > (10. * ndens_w)
cloud_mass = np.sum(cloud_ndens[condition]*cloud_volume[condition])*(mu*mp)
mass_data.append([float(key), cloud_mass, cloud_com])
mass_data = np.array(mass_data)
line, = plt.plot(mass_data[:till+1,0], mass_data[:till+1,1]/mass_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(mass_data[:till+1,0], (mass_data[:till+1,2]/mass_data[0,2])**0.845, linestyle ='--', color = line.get_color())
plt.legend(loc="upper left", title=r"$t_{\rm cool, mix}/t_{\rm cc}|_{\rm ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
plt.ylim(ymin=5.0e-02)
plt.xlim(xmin=0., xmax = till)
plt.yscale('log')
plt.xlabel(r"time [$t_{\rm cc,ini}$]")
plt.ylabel(r"Cloud mass $A_{\rm cold}$ ($T<8\times 10^4$ K, $\rho > 10 \rho _{\rm wind}$) [$M_{\rm cold, ini}$]")
plt.savefig(f"mass-cold+dense-trunc_time{'-dark' if dark else ''}.svg", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 100
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# t/tcc Area Volume
area_data = []
for key in list(cloud_data.keys()):
area_data.append([float(key), cloud_data[key]['cloud_tot_surface_area'], cloud_data[key]['cloud_tot_vol']])
area_data = np.array(area_data)
plt.plot(area_data[:till+1,2]/area_data[0,2], area_data[:till+1,1]/area_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}')
if tcool_mix_B_tcc[select] == 0.2:
index = 2/3.
K1, K2 = 6./10.**index, 40./10.**index
plt.fill_between(area_data[:till+1,2]/area_data[0,2], K1*(area_data[:till+1,2]/area_data[0,2])**index, K2*(area_data[:till+1,2]/area_data[0,2])**index,
color="lemonchiffon", alpha=0.8, zorder=-100)
'''
K1, K2 = 6./10.**(2/3.), 40./10.**(2/3.)
plt.fill_between(area_data[:till+1,2]/area_data[0,2], K1*(area_data[:till+1,2]/area_data[0,2])**(2/3.), K2*(area_data[:till+1,2]/area_data[0,2])**(2/3.),
color="gray", alpha=0.3, zorder=-100)
'''
plt.legend(loc="upper left", title=r"$t_{\rm cool, mix}/t_{\rm cc}|_{\rm ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
plt.ylim(ymin=0.7)
plt.xlim(xmin=0.7)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r"Volume $V_{\rm cold}$ ($T<8\times 10^4$ K) [$V_{\rm cold, ini}$]")
plt.ylabel(r"Surface area $A_{\rm cold}$ ($T<8\times 10^4$ K) [$A_{\rm cold, ini}$]")
plt.savefig(f"area-cold-trunc_volume{'-dark' if dark else ''}.svg", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# t/tcc Mass Area Velocity Position Mdot
area_data = []
for key in list(cloud_data.keys()):
cloud_velocity_r = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_r'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity_th = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_th'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity_ph = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_ph'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity = np.sqrt( cloud_velocity_r**2 + cloud_velocity_th**2 + cloud_velocity_ph**2 )
cloud_com = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_distance'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_mdot = np.sum(-cloud_data[key]['cloud_surface_density']*cloud_data[key]['cloud_tot_surface_area']*cloud_data[key]['cloud_surface_vin_elems'])
area_data.append([float(key),
cloud_data[key]['cloud_tot_mass'],
cloud_data[key]['cloud_tot_surface_area'],
cloud_velocity,
cloud_com,
np.abs(cloud_mdot),
np.average(cloud_data[key]['cloud_surface_vin_elems']),
np.average(cloud_data[key]['cloud_surface_density']),
])
area_data = np.array(area_data)
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 100
wind_vel = CC85windvel(area_data[:till+1,4]/distance_ini * diniBdinj)/CC85windvel(diniBdinj) # code units
mdot_from_mass = np.gradient(area_data[:till+1,0]*tcc, area_data[:till+1,1])/(area_data[0,1]/tcc)
mass_smoothed = scipy.signal.savgol_filter(area_data[:till+1,1], window_length=6, polyorder=4, deriv=0)
mdot_from_mass = scipy.signal.savgol_filter(area_data[:till+1,1], window_length=6, polyorder=4, deriv=1)/area_data[0,1]
# smooth mass test
# line, = plt.plot(area_data[:till+1,0], area_data[:till+1,1]/area_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,0], mass_smoothed/mass_smoothed[0], color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,0], mdot_from_mass, label=f'{tcool_mix_B_tcc[select]:.2f}') # dM/dt [M0/tcc]
# mdot_from_mass = np.gradient(area_data[:till+1,0]*tcc, area_data[:till+1,1])/(area_data[0,1]/tcc)
# plt.plot(area_data[:till+1,0], mdot_from_mass, color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,4]/area_data[0,4], mdot_from_mass, label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,5]/(area_data[0,1]/tcc), color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,4]/area_data[0,4], -mdot_from_mass/(area_data[:till+1,7]*area_data[:till+1,2]), label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,6], color=line.get_color(), linestyle="--")
# velocity
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,3]/cs, label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], wind_vel/cs, color='k', linestyle="--")
# plt.plot(area_data[:till+1,4]/area_data[0,4], np.abs(area_data[:till+1,3]-wind_vel)/wind_vel, label=f'{tcool_mix_B_tcc[select]:.2f}') # del_v
# mass-distance
plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,1]/area_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}', zorder=100+int(8%(select+1)))
# plt.plot(area_data[:till+1,4]/area_data[0,4], K0*(area_data[:till+1,4]/area_data[0,4]), color='k', linestyle="--")
if tcool_mix_B_tcc[select] == 0.2:
K0, K1, K2 = 1.71/1.46, 1.08/1.75, 2.02/1.29
plt.fill_between(area_data[:till+1,4]/area_data[0,4], K1*(area_data[:till+1,4]/area_data[0,4]), K2*(area_data[:till+1,4]/area_data[0,4]),
color="lemonchiffon", alpha=0.8, zorder=-100)
K0, K1, K2 = 1.71/1.46**3, 1.51/1.46**3, 1.71/1.14**3
plt.fill_between(area_data[:till+1,4]/area_data[0,4], K1*(area_data[:till+1,4]/area_data[0,4])**3, K2*(area_data[:till+1,4]/area_data[0,4])**3,
color="gray", alpha=0.3, zorder=-100)
plt.legend(loc="best", title=r"$t_{\rm cool, mix}/t_{\rm cc}|_{\rm ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
plt.ylim(ymin=0.6, ymax=22.0)
plt.xlim(xmin=0.98, xmax=9.9)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r"distance traveled [$d_{\rm ini}$]")
plt.ylabel(r"Cold mass $M_{\rm cold}$ ($T<8\times 10^4$ K) [$M_{\rm cold, ini}$]")
plt.savefig(f"mass-cold-trunc_distance{'-dark' if dark else ''}.svg", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# t/tcc Mass Area Velocity Position Mdot
area_data = []
for key in list(cloud_data.keys()):
cloud_velocity_r = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_r'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity_th = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_th'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity_ph = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_ph'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity = np.sqrt( cloud_velocity_r**2 + cloud_velocity_th**2 + cloud_velocity_ph**2 )
cloud_com = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_distance'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_mdot = np.sum(-cloud_data[key]['cloud_surface_density']*cloud_data[key]['cloud_tot_surface_area']*cloud_data[key]['cloud_surface_vin_elems'])
area_data.append([float(key),
cloud_data[key]['cloud_tot_mass'],
cloud_data[key]['cloud_tot_surface_area'],
cloud_velocity,
cloud_com,
np.abs(cloud_mdot),
np.average(cloud_data[key]['cloud_surface_vin_elems']),
np.average(cloud_data[key]['cloud_surface_density']),
])
area_data = np.array(area_data)
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 100
wind_vel = CC85windvel(area_data[:till+1,4]/distance_ini * diniBdinj)/CC85windvel(diniBdinj) # code units
mdot_from_mass = np.gradient(area_data[:till+1,0]*tcc, area_data[:till+1,1])/(area_data[0,1]/tcc)
mass_smoothed = scipy.signal.savgol_filter(area_data[:till+1,1], window_length=6, polyorder=4, deriv=0)
mdot_from_mass = scipy.signal.savgol_filter(area_data[:till+1,1], window_length=6, polyorder=4, deriv=1)/area_data[0,1]
# smooth mass test
# line, = plt.plot(area_data[:till+1,0], area_data[:till+1,1]/area_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,0], mass_smoothed/mass_smoothed[0], color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,0], mdot_from_mass, label=f'{tcool_mix_B_tcc[select]:.2f}') # dM/dt [M0/tcc]
# mdot_from_mass = np.gradient(area_data[:till+1,0]*tcc, area_data[:till+1,1])/(area_data[0,1]/tcc)
# plt.plot(area_data[:till+1,0], mdot_from_mass, color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,4]/area_data[0,4], mdot_from_mass, label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,5]/(area_data[0,1]/tcc), color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,4]/area_data[0,4], -mdot_from_mass/(area_data[:till+1,7]*area_data[:till+1,2]), label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,6], color=line.get_color(), linestyle="--")
# velocity
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,3]/cs, label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], wind_vel/cs, color='k', linestyle="--")
plt.plot(area_data[:till+1,4]/area_data[0,4], np.abs(area_data[:till+1,3]-wind_vel)/wind_vel, label=f'{tcool_mix_B_tcc[select]:.2f}',
zorder=100+int(8%(select+1))) # del_v
# mass-distance
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,1]/area_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], K0*(area_data[:till+1,4]/area_data[0,4]), color='k', linestyle="--")
plt.legend(loc="best", title=r"$t_{\rm cool, mix}/t_{\rm cc}|_{\rm ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
plt.ylim(ymin=0.06)#, ymax=22.0)
plt.xlim(xmin=0.98, xmax=9.9)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r"distance traveled [$d_{\rm ini}$]")
plt.ylabel(r"Cold mass velocity $|v_{\rm cold}-v_{\rm wind}|$ ($T<8\times 10^4$ K) [$v_{\rm wind}$]")
plt.savefig(f"vel_diff-cold-trunc_distance{'-dark' if dark else ''}.svg", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# t/tcc Mass Area Velocity Position Mdot
area_data = []
for key in list(cloud_data.keys()):
cloud_velocity_r = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_r'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity_th = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_th'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity_ph = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_velocity_ph'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_velocity = np.sqrt( cloud_velocity_r**2 + cloud_velocity_th**2 + cloud_velocity_ph**2 )
cloud_com = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_distance'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
cloud_mdot = np.sum(-cloud_data[key]['cloud_surface_density']*cloud_data[key]['cloud_tot_surface_area']*cloud_data[key]['cloud_surface_vin_elems'])
area_data.append([float(key),
cloud_data[key]['cloud_tot_mass'],
cloud_data[key]['cloud_tot_surface_area'],
cloud_velocity,
cloud_com,
np.abs(cloud_mdot),
np.average(cloud_data[key]['cloud_surface_vin_elems']),
np.average(cloud_data[key]['cloud_surface_density']),
])
area_data = np.array(area_data)
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 100
wind_vel = CC85windvel(area_data[:till+1,4]/distance_ini * diniBdinj)/CC85windvel(diniBdinj) # code units
mdot_from_mass = np.gradient(area_data[:till+1,0]*tcc, area_data[:till+1,1])/(area_data[0,1]/tcc)
mass_smoothed = scipy.signal.savgol_filter(area_data[:till+1,1], window_length=6, polyorder=4, deriv=0)
mdot_from_mass = scipy.signal.savgol_filter(area_data[:till+1,1], window_length=6, polyorder=4, deriv=1)/area_data[0,1]
# smooth mass test
# line, = plt.plot(area_data[:till+1,0], area_data[:till+1,1]/area_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,0], mass_smoothed/mass_smoothed[0], color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,0], mdot_from_mass, label=f'{tcool_mix_B_tcc[select]:.2f}') # dM/dt [M0/tcc]
# mdot_from_mass = np.gradient(area_data[:till+1,0]*tcc, area_data[:till+1,1])/(area_data[0,1]/tcc)
# plt.plot(area_data[:till+1,0], mdot_from_mass, color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,4]/area_data[0,4], mdot_from_mass, label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,5]/(area_data[0,1]/tcc), color=line.get_color(), linestyle="--")
# line, = plt.plot(area_data[:till+1,4]/area_data[0,4], -mdot_from_mass/(area_data[:till+1,7]*area_data[:till+1,2]), label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,6], color=line.get_color(), linestyle="--")
# velocity
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,3]/cs, label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], wind_vel/cs, color='k', linestyle="--")
plt.plot(area_data[:till+1,0], np.abs(area_data[:till+1,3]-wind_vel)/wind_vel, label=f'{tcool_mix_B_tcc[select]:.2f}',
zorder=100+int(8%(select+1))) # del_v
# mass-distance
# plt.plot(area_data[:till+1,4]/area_data[0,4], area_data[:till+1,1]/area_data[0,1], label=f'{tcool_mix_B_tcc[select]:.2f}')
# plt.plot(area_data[:till+1,4]/area_data[0,4], K0*(area_data[:till+1,4]/area_data[0,4]), color='k', linestyle="--")
plt.legend(loc="best", title=r"$t_{\rm cool, mix}/t_{\rm cc}|_{\rm ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
plt.ylim(ymin=0.06)#, ymax=22.0)
plt.xlim(xmin=0., xmax = till)
plt.yscale('log')
plt.xlabel(r"time [$t_{\rm cc,ini}$]")
plt.ylabel(r"Cold mass velocity $|v_{\rm cold}-v_{\rm wind}|$ ($T<8\times 10^4$ K) [$v_{\rm wind}$]")
plt.savefig(f"vel_diff-cold-trunc_time{'-dark' if dark else ''}.svg", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# t/tcc vin distance
area_data = []
for key in list(cloud_data.keys()):
cloud_com = np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems']*cloud_data[key]['cloud_distance'])/np.sum(cloud_data[key]['cloud_density']*cloud_data[key]['cloud_volume_elems'])
area_data.append([float(key),
-cloud_data[key]['vin_avg'],
cloud_com,
])
area_data = np.array(area_data)
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 100
wind_vel = CC85windvel(area_data[:till+1,2]/distance_ini * diniBdinj)/CC85windvel(diniBdinj) # code units
plt.plot(area_data[:till+1,2]/area_data[0,2], np.abs(area_data[:till+1,1])/wind_vel, label=f'{tcool_mix_B_tcc[select]:.2f}',
zorder=100+int(8%(select+1))) # v_in
if tcool_mix_B_tcc[select] == 0.2:
index = -gamma/2
K1, K2 = 1.0e-02/2.0**index, 7.0e-03/1.0**index
plt.fill_between(area_data[:till+1,2]/area_data[0,2], K1*(area_data[:till+1,2]/area_data[0,2])**index, K2*(area_data[:till+1,2]/area_data[0,2])**index,
color="lemonchiffon", alpha=0.8, zorder=-100)
plt.legend(loc="best", title=r"$t_{\rm cool, mix}/t_{\rm cc}|_{\rm ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
# plt.ylim(ymin=0.06)#, ymax=22.0)
plt.xlim(xmin=0.98, xmax=9.9)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r"distance traveled [$d_{\rm ini}$]")
plt.ylabel(r"Inflow velocity $v_{\rm in}$ ($T<8\times 10^4$ K) [$v_{\rm wind}$]")
plt.savefig(f"vel_inf-cold-trunc_distance{'-dark' if dark else ''}.svg", transparent=False)
# plt.show()
plt.close()
os.makedirs(f"profiles-paraview", exist_ok = True)
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
ndens_w_ini = UNIT_DENSITY/(mu*mp)
ndens_cl_ini = chi * ndens_w_ini # cgs
prs_w_ini = ndens_w_ini*chi*Tcl # p/kB
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# print(f"./paraview-cloud-analysis_data-dump/{label}.pickle")
# print(list(cloud_data.keys()))
count = 0
distance_min, distance_max = None, None
for key in list(cloud_data.keys()):
if int(key)%8 > 1.0e-06:
continue
cloud_density = cloud_data[key]['cloud_density']* UNIT_DENSITY/(mu*mp)/ndens_w_ini
cloud_distance = cloud_data[key]['cloud_distance']/distance_ini
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 64
if float(key)>till:
continue
if count==0:
distance_min = np.min(cloud_distance)
distance_max = np.max(cloud_distance)
else:
distance_min = distance_min if (_tmp:=np.min(cloud_distance))>distance_min else _tmp
distance_max = distance_max if (_tmp:=np.max(cloud_distance))<distance_max else _tmp
count += 1
line, = plt.plot(cloud_distance, cloud_density, marker='.', markersize=4, linestyle='None')
plt.plot([], [], color=line.get_color(), label=f"{key}")
'''
if distance_min<0.5*diniBdinj:
distance_min = 0.5*diniBdinj
if distance_max<1.5*diniBdinj:
distance_max = 1.5*diniBdinj
'''
ndens_w = (CC85windrho(distance_w:=np.linspace(0.98*distance_min, distance_max, 1000) * diniBdinj)/CC85windrho(diniBdinj)) * UNIT_DENSITY/(mu*mp) # cgs
prs_w = (CC85windprs(distance_w*diniBdinj)/CC85windprs(diniBdinj)) * prs_w_ini
plt.plot(distance_w, ndens_w/ndens_w_ini, linestyle="--", color="black", linewidth=2.0)
plt.plot(distance_w, prs_w/prs_w_ini*chi, linestyle=":", color="gray", linewidth=2.0)
plt.legend(loc="best", title=r"$t/t_{\rm cc,ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
# plt.ylim(ymin=0.06)#, ymax=22.0)
# plt.xlim(xmin=distance_min, xmax=distance_max)
plt.xlim(xmin=0.8, xmax=1.01*distance_max)
# plt.xscale('log')
plt.yscale('log')
plt.xlabel(r"distance traveled [$d_{\rm ini}$]")
plt.ylabel(r"Density [initial wind density]")
plt.savefig(f"profiles-paraview/density-profile_{tcool_mix_B_tcc[select]:.2f}{'-dark' if dark else ''}.png", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
ndens_w_ini = UNIT_DENSITY/(mu*mp) # cgs
ndens_cl_ini = chi * ndens_w_ini
prs_w_ini = ndens_w_ini*chi*Tcl # p/kB
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# print(f"./paraview-cloud-analysis_data-dump/{label}.pickle")
# print(list(cloud_data.keys()))
count = 0
distance_min, distance_max = None, None
for key in list(cloud_data.keys()):
if int(key)%8 > 1.0e-06:
continue
cloud_pressure = cloud_data[key]['cloud_pressure']* UNIT_DENSITY*UNIT_VELOCITY**2/kB/prs_w_ini
cloud_distance = cloud_data[key]['cloud_distance']/distance_ini
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 64
if float(key)>till:
continue
if count==0:
distance_min = np.min(cloud_distance)
distance_max = np.max(cloud_distance)
else:
distance_min = distance_min if (_tmp:=np.min(cloud_distance))>distance_min else _tmp
distance_max = distance_max if (_tmp:=np.max(cloud_distance))<distance_max else _tmp
count += 1
line, = plt.plot(cloud_distance, cloud_pressure, marker='.', markersize=4, linestyle='None')
plt.plot([], [], color=line.get_color(), label=f"{key}")
'''
if distance_min<0.5*diniBdinj:
distance_min = 0.5*diniBdinj
if distance_max<1.5*diniBdinj:
distance_max = 1.5*diniBdinj
'''
ndens_w = (CC85windrho(distance_w:=np.linspace(0.98*distance_min, distance_max, 1000) * diniBdinj)/CC85windrho(diniBdinj)) * UNIT_DENSITY/(mu*mp) # cgs
prs_w = (CC85windprs(distance_w*diniBdinj)/CC85windprs(diniBdinj)) * prs_w_ini
plt.plot(distance_w, prs_w/prs_w_ini, linestyle="--", color="black", linewidth=2.0)
plt.legend(loc="best", title=r"$t/t_{\rm cc,ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
# plt.ylim(ymin=0.06)#, ymax=22.0)
# plt.xlim(xmin=distance_min, xmax=distance_max)
plt.xlim(xmin=0.8, xmax=1.01*distance_max)
# plt.xscale('log')
plt.yscale('log')
plt.xlabel(r"distance traveled [$d_{\rm ini}$]")
plt.ylabel(r"Pressure [initial wind pressure]")
plt.savefig(f"profiles-paraview/pressure-profile_{tcool_mix_B_tcc[select]:.2f}{'-dark' if dark else ''}.png", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
ndens_w_ini = UNIT_DENSITY/(mu*mp) # cgs
ndens_cl_ini = chi * ndens_w_ini
prs_w_ini = ndens_w_ini*chi*Tcl # p/kB
v_w_ini = UNIT_VELOCITY # cgs
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# print(f"./paraview-cloud-analysis_data-dump/{label}.pickle")
# print(list(cloud_data.keys()))
count = 0
distance_min, distance_max = None, None
for key in list(cloud_data.keys()):
if int(key)%8 > 1.0e-06:
continue
cloud_velocity = np.sqrt(cloud_data[key]['cloud_velocity_r']**2 + cloud_data[key]['cloud_velocity_th']**2 + cloud_data[key]['cloud_velocity_ph']**2) * UNIT_VELOCITY
cloud_distance = cloud_data[key]['cloud_distance']/distance_ini
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 64
if float(key)>till:
continue
if count==0:
distance_min = np.min(cloud_distance)
distance_max = np.max(cloud_distance)
else:
distance_min = distance_min if (_tmp:=np.min(cloud_distance))>distance_min else _tmp
distance_max = distance_max if (_tmp:=np.max(cloud_distance))<distance_max else _tmp
count += 1
line, = plt.plot(cloud_distance, cloud_velocity/v_w_ini, marker='.', markersize=4, linestyle='None')
plt.plot([], [], color=line.get_color(), label=f"{key}")
'''
if distance_min<0.5*diniBdinj:
distance_min = 0.5*diniBdinj
if distance_max<1.5*diniBdinj:
distance_max = 1.5*diniBdinj
'''
v_w = (CC85windvel(distance_w:=np.linspace(0.98*distance_min, distance_max, 1000)*diniBdinj)/CC85windvel(diniBdinj)) * v_w_ini
plt.plot(distance_w, v_w/v_w_ini, linestyle="--", color="black", linewidth=2.0)
plt.legend(loc="best", title=r"$t/t_{\rm cc,ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
# plt.ylim(ymin=0.06)#, ymax=22.0)
# plt.xlim(xmin=distance_min, xmax=distance_max)
plt.xlim(xmin=0.8, xmax=1.01*distance_max)
# plt.xscale('log')
# plt.yscale('log')
plt.xlabel(r"distance traveled [$d_{\rm ini}$]")
plt.ylabel(r"Velocity [initial wind velocity]")
plt.savefig(f"profiles-paraview/velocity-profile_{tcool_mix_B_tcc[select]:.2f}{'-dark' if dark else ''}.png", transparent=False)
# plt.show()
plt.close()
for select in range(len(tcool_mix_B_tcc)):
label = f"c{chi:d},m{mach:.3f},T4e4,t{tcool_mix_B_tcc[select]:.2f},r{rini_B_rcl[select]:.3f}"
directory = f"{root}-{label}"
UNIT_LENGTH = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_LENGTH").split()[-1])
UNIT_DENSITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_DENSITY").split()[-1])
UNIT_VELOCITY = float(sp.getoutput(f"cat {directory}/definitions.h | grep UNIT_VELOCITY").split()[-1])
distance_ini = float(sp.getoutput(f"cat {directory}/pluto.ini | grep RINI").split()[-1])
ndens_w_ini = UNIT_DENSITY/(mu*mp) # cgs
ndens_cl_ini = chi * ndens_w_ini
prs_w_ini = ndens_w_ini*chi*Tcl # p/kB
v_w_ini = UNIT_VELOCITY # cgs
with open(f"./paraview-cloud-analysis_data-dump/{label}.pickle", "rb") as handle:
cloud_data = pickle.load(handle)
# print(f"./paraview-cloud-analysis_data-dump/{label}.pickle")
# print(list(cloud_data.keys()))
count = 0
distance_min, distance_max = None, None
for key in list(cloud_data.keys()):
if int(key)%8 > 1.0e-06:
continue
cloud_velocity = np.sqrt(cloud_data[key]['cloud_velocity_r']**2 + cloud_data[key]['cloud_velocity_th']**2 + cloud_data[key]['cloud_velocity_ph']**2) * UNIT_VELOCITY
cloud_distance = cloud_data[key]['cloud_distance']/distance_ini
if tcool_mix_B_tcc[select]==0.1:
till = 25
elif tcool_mix_B_tcc[select]==0.2:
till = 62
else:
till = 64
if float(key)>till:
continue
if count==0:
distance_min = np.min(cloud_distance)
distance_max = np.max(cloud_distance)
else:
distance_min = distance_min if (_tmp:=np.min(cloud_distance))>distance_min else _tmp
distance_max = distance_max if (_tmp:=np.max(cloud_distance))<distance_max else _tmp
count += 1
v_w = (CC85windvel(cloud_distance*diniBdinj)/CC85windvel(diniBdinj)) * v_w_ini
line, = plt.plot(cloud_distance, (cloud_velocity-v_w)/v_w, marker='.', markersize=4, linestyle='None')
plt.plot([], [], color=line.get_color(), label=f"{key}")
'''
if distance_min<0.5*diniBdinj:
distance_min = 0.5*diniBdinj
if distance_max<1.5*diniBdinj:
distance_max = 1.5*diniBdinj
'''
v_w = (CC85windvel(distance_w:=np.linspace(0.98*distance_min, distance_max, 1000)*diniBdinj)/CC85windvel(diniBdinj)) * v_w_ini
# plt.plot(distance_w, v_w/v_w_ini, linestyle="--", color="black", linewidth=2.0)
plt.legend(loc="best", title=r"$t/t_{\rm cc,ini}$", ncols=3,
prop = { "size": 20 }, title_fontsize=22, fancybox=True)
# plt.ylim(ymin=0.06)#, ymax=22.0)
# plt.xlim(xmin=distance_min, xmax=distance_max)
plt.xlim(xmin=0.8, xmax=1.01*distance_max)
# plt.xscale('log')
# plt.yscale('log')
plt.xlabel(r"distance traveled [$d_{\rm ini}$]")
plt.ylabel(r"$\Delta v = v_{\rm cl} - v_{\rm wind}$ [wind velocity]")
plt.savefig(f"profiles-paraview/velocity_delta-profile_{tcool_mix_B_tcc[select]:.2f}{'-dark' if dark else ''}.png", transparent=False)
# plt.show()
plt.close()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment