noe/pipeline
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

[fragdisk] [galactica] fix units, add velocity plots

Browse Source
This commit is contained in:
Noe Brucy
2020-12-16 09:51:23 +01:00
parent b2edc9ba54
commit 854f8da9f3
1 changed files with 149 additions and 105 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fragdisk_galactica.py
+149 -105
View File
@@ -28,15 +28,17 @@
from __future__ import print_function, unicode_literals from __future__ import print_function, unicode_literals
import os import os
import argparse
import h5py import h5py
import numpy as np import numpy as np
import matplotlib as mpl
from astrophysix import units as U from astrophysix import units as U
from astrophysix.simdm import Project, ProjectCategory, SimulationStudy from astrophysix.simdm import Project, ProjectCategory, SimulationStudy
from astrophysix.simdm.datafiles import PlotInfo, PlotType from astrophysix.simdm.datafiles import PlotInfo, PlotType
from astrophysix.utils.file import FileType from astrophysix.utils.file import FileType
from plotter import Plotter, default_params from plotter import Plotter
from pp_params import default_params
from ramses_astrophysix import ramses from ramses_astrophysix import ramses
pp_params = default_params() pp_params = default_params()
@@ -53,6 +55,9 @@ pp_params.pdf.nb_bin = 100
pp_params.process.num_process = 10 pp_params.process.num_process = 10
in_dir = "/drf/projets/alfven-data/nbrucy/simus/fragdisk" in_dir = "/drf/projets/alfven-data/nbrucy/simus/fragdisk"
out_dir = "/dsm/anais/storageA/nbrucy/visus/fragdisk/mnras" out_dir = "/dsm/anais/storageA/nbrucy/visus/fragdisk/mnras"
in_dir_conv = "/drf/projets/alfven-data/nbrucy/simus/conv_disk"
out_dir_conv = "/dsm/anais/storageA/nbrucy/visus/conv_disk"
nml_key = "cloud_params/beta_cool" nml_key = "cloud_params/beta_cool"
@@ -177,32 +182,36 @@ pl_units = Plotter(
info = pl_units.comp.info info = pl_units.comp.info
rd = U.Unit.create_unit( rd = U.Unit.create_unit(
"r_d", "r_d",
latex="r_\\mathrm{d}",
base_unit=0.25 * info["unit_length"] / info["boxlen"], base_unit=0.25 * info["unit_length"] / info["boxlen"],
) )
rhod = U.Unit.create_unit( rhod = U.Unit.create_unit(
"rho_d", "rho_d",
latex="\rho_d", latex="\\rho_\\mathrm{d}",
base_unit=6.36 * info["unit_density"], base_unit=6.36 * info["unit_density"],
) )
rho_u = U.Unit.create_unit( rho_u = U.Unit.create_unit(
"M_*/r_d^3", "rho_u",
latex="M_\\star.r_d^{-3}", latex="\mathrm{M}_\\star.r_\\mathrm{d}^{-3}",
base_unit=info["unit_mass"] * rd ** (-3), base_unit=info["unit_mass"] * rd ** (-3),
) )
orp = U.Unit.create_unit(
"ORP",
base_unit=2 * np.pi * np.sqrt(0.25 ** 3) * info["unit_time"],
)
vkd = U.Unit.create_unit(
"vkd",
latex="v_\\mathrm{k,d}",
base_unit=2 * np.pi * rd / orp,
)
Pd = U.Unit.create_unit( Pd = U.Unit.create_unit(
"P_d", "P_d",
base_unit=(0.2 * info["unit_velocity"]) ** 2 * rhod, base_unit=(0.2 * info["unit_velocity"]) ** 2 * rhod,
) )
orp = U.Unit.create_unit(
"ORP",
base_unit=2 * np.pi * 0.25 ** -3 * info["unit_time"],
)
vkd = rd * orp
P_u = U.Unit.create_unit( P_u = U.Unit.create_unit(
"(k_kd * M_*)/r_d^3", "P_u",
latex="v_{k,d}^2.M_\\star.r_d^{-3}", latex="v_\\mathrm{k,d}^2.\\mathrm{M}_\\star.r_\\mathrm{d}^{-3}",
base_unit=vkd * rho_u, base_unit=vkd ** 2 * rho_u,
) )
Sigmad = U.Unit.create_unit( Sigmad = U.Unit.create_unit(
@@ -211,7 +220,7 @@ Sigmad = U.Unit.create_unit(
base_unit=0.25 * info["unit_density"] * info["unit_length"] / info["boxlen"], base_unit=0.25 * info["unit_density"] * info["unit_length"] / info["boxlen"],
) )
Sigma_u = U.Unit.create_unit( Sigma_u = U.Unit.create_unit(
"M*/r_d^2", "Sigma_u",
latex="M_\\star.r_d^{-2}", latex="M_\\star.r_d^{-2}",
base_unit=rd ** (-2) * info["unit_mass"], base_unit=rd ** (-2) * info["unit_mass"],
) )
@@ -219,82 +228,95 @@ Sigma_u = U.Unit.create_unit(
# ---- Runs ----- # ---- Runs -----
# JR13_TIC groups = ["jr11", "jr12", "jr12_tic", "jr13_tic"]
# groups = ["jr13_tic"]
pls = []
pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}" for group in groups:
pp_params.astrophysix.descr_fmt = """ if group == "jr13_tic":
<p>Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.</p> # JR13_TIC
"""
runs = "*_jr13" pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}"
pl_jr13 = Plotter( pp_params.astrophysix.descr_fmt = """
in_dir, <p>Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.</p>
filter_name=runs, """
in_nums="all",
sort_run_by=nml_key,
path_out=out_dir,
tag="jr13_tic_mnras",
pp_params=pp_params,
unit_time=orp,
)
print("JR13_TIC defined")
# JR12_TIC runs = "*_jr13"
runs_12 = "0[0-9][0-9]_beta*_jr12" pl = Plotter(
pl_jr12_tic = Plotter( in_dir,
in_dir, filter_name=runs,
filter_name=runs_12, in_nums="all",
in_nums="all", sort_run_by=nml_key,
sort_run_by=nml_key, path_out=out_dir,
filter_nml=("cloud_params", "!=", 7), tag="jr13_tic_mnras",
path_out=out_dir, pp_params=pp_params,
tag="jr12_tic_mnras", unit_time=orp,
pp_params=pp_params, )
unit_time=orp,
)
pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}" if group == "jr12_tic":
pp_params.astrophysix.descr_fmt = """ # JR12_TIC
<p>Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}</p> pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}"
""" pp_params.astrophysix.descr_fmt = """
print("JR12_TIC defined") <p>Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.</p>
"""
runs_12 = "0[0-9][0-9]_beta*_jr12"
pl = Plotter(
in_dir,
filter_name=runs_12,
in_nums="all",
sort_run_by=nml_key,
filter_nml=("cloud_params", "!=", 7),
path_out=out_dir,
tag="jr12_tic_mnras",
pp_params=pp_params,
unit_time=orp,
)
# JR12 if group == "jr12":
in_dir_conv = "/drf/projets/alfven-data/nbrucy/simus/conv_disk"
out_dir_conv = "/dsm/anais/storageA/nbrucy/visus/conv_disk"
runs = "[7-8][0-9]_beta*_j*"
pl_jr12 = Plotter(
in_dir_conv,
filter_name=runs,
in_nums="all",
sort_run_by=nml_key,
path_out=out_dir,
tag="jr12_mnras",
pp_params=pp_params,
unit_time=orp,
)
print("JR12 defined")
# JR11 pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}"
runs = "*beta*_jr11" pp_params.astrophysix.descr_fmt = """
pl_l11 = Plotter( <p>Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}</p>
in_dir_conv, """
filter_name=runs,
sort_run_by=nml_key,
filter_nml=("cloud_params/beta_cool", ">", 3),
path_out=out_dir_conv,
tag="jr11_mnras",
pp_params=pp_params,
unit_time=orp,
)
print("JR11 defined") # JR12
runs = "[7-8][0-9]_beta*_j*"
pl = Plotter(
in_dir_conv,
filter_name=runs,
in_nums="all",
sort_run_by=nml_key,
path_out=out_dir,
tag="jr12_mnras",
pp_params=pp_params,
unit_time=orp,
)
if group == "jr11":
pls = [pl_l11, pl_jr12, pl_jr12_tic, pl_jr13] pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}"
pp_params.astrophysix.descr_fmt = """
<p>Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}</p>
"""
# JR11
runs = "*beta*_jr11"
pl = Plotter(
in_dir_conv,
filter_name=runs,
sort_run_by=nml_key,
filter_nml=("cloud_params/beta_cool", ">", 3),
path_out=out_dir_conv,
tag="jr11_mnras",
pp_params=pp_params,
unit_time=orp,
)
pls.append(pl)
print("{} defined".format(group.upper()))
# -------------------------------------------------------------------------------------------------------------------- # # -------------------------------------------------------------------------------------------------------------------- #
redo = True redo = True
for pl in pls: for pl in pls:
pl.pp_params.process.verbose = True pl.pp_params.process.verbose = True
@@ -307,13 +329,17 @@ for pl in pls:
# -------------------------------------------------------------------------------------------------------------------- # # -------------------------------------------------------------------------------------------------------------------- #
select = {
"time": 4.5,
# "filter_nml" : (nml_key, "=", 8),
}
for pl in pls: for pl in pls:
select = {"time": 4.5}
# Edit descriptions # Edit descriptions
pl.rules["slice_rho"].description = "Density slice" pl.rules["slice_rho"].description = "Density slice"
pl.rules["slice_P"].description = "Pressure slice" pl.rules["slice_P"].description = "Pressure slice"
pl.rules["slice_velr"].description = "Radial velocity slice"
pl.rules["slice_velphi"].description = "Orthoradial velocity slice"
pl.rules[ pl.rules[
"pdf_coldens" "pdf_coldens"
].description = """ ].description = """
@@ -321,42 +347,49 @@ for pl in pls:
$\\sigma = \\Sigma/\\overline{\\Sigma}$ with respect to its azimuthal average $\\sigma = \\Sigma/\\overline{\\Sigma}$ with respect to its azimuthal average
""" """
coldens_kwargs = { # define plot parameters
map_kwargs = {
"overwrite": redo, "overwrite": redo,
"overwrite_dep": False, "overwrite_dep": False,
"unit_space": rd, "unit_space": rd,
"center_space": True, "center_space": True,
"unit_time": orp, "unit_time": orp,
"unit": Sigma_u,
"nml_key": "cloud_params/beta_cool", "nml_key": "cloud_params/beta_cool",
"select": select, "select": select,
"label": r"$\Sigma$", }
coldens_kwargs = {
**map_kwargs,
**{
"unit": Sigma_u,
"label": r"$\Sigma$",
},
} }
rho_kwargs = { rho_kwargs = {
"overwrite": redo, **map_kwargs,
"overwrite_dep": False, **{
"unit_space": rd, "unit": rho_u,
"center_space": True, "label": r"$\rho$",
"unit_time": orp, },
"unit": rho_u,
"nml_key": "cloud_params/beta_cool",
"select": select,
"label": r"$\rho$",
} }
P_kwargs = { P_kwargs = {
"overwrite": redo, **map_kwargs,
"overwrite_dep": False, **{
"unit_space": rd, "unit": P_u,
"center_space": True, "label": r"$P$",
"unit_time": orp, },
"unit": P_u,
"nml_key": "cloud_params/beta_cool",
"select": select,
"label": r"$P$",
} }
vel_kwargs = {
**map_kwargs,
**{
"unit": vkd,
},
}
# do plots
pl.coldens("z", **coldens_kwargs) pl.coldens("z", **coldens_kwargs)
pl.coldens("y", **coldens_kwargs) pl.coldens("y", **coldens_kwargs)
@@ -365,6 +398,15 @@ for pl in pls:
pl.slice_P("z", **P_kwargs) pl.slice_P("z", **P_kwargs)
pl.slice_velr(
"z",
label=r"$v_r$",
norm=mpl.colors.SymLogNorm(0.1, vmin=-2, vmax=2, base=10),
autoscale=False,
**vel_kwargs,
)
pl.slice_velphi("z", label=r"$v_\varphi$", **vel_kwargs)
pl.pdf_coldens( pl.pdf_coldens(
"z", "z",
overwrite=redo, overwrite=redo,
@@ -400,15 +442,17 @@ for simu in proj.simulations:
# Change unit to things Galactica understands # Change unit to things Galactica understands
snap.time = (snap.time[0], U.none) snap.time = (snap.time[0], U.none)
snap.description = """ snap.description = """
<p>Snapshot at {time} ORPs.</p> <p>Snapshot at {time:.3g} ORPs.</p>
<h5 class="sn_panel">Notations:</h5>
<h5> Notations </h5>
<p> <p>
$r_d$ : radius of the disk <br /> $r_d$ : radius of the disk <br />
ORP : outer rotation period, that is period of the gas at $r_d$ <br /> ORP : outer rotation period, that is period of the gas at $r_d$ <br />
$M_\\star$ : Mass of the central object <br /> $M_\\star$ : Mass of the central object <br />
$v_{kepl}$ : Keplerian speed at $r_d$ $v_{kepl}$ : Keplerian speed at $r_d$
</p> </p>
<p>All slices are taken at $z=0$ or $y=0$.
""".format( """.format(
time=snap.time[0], kepl="{k,d}" time=snap.time[0], kepl="{k,d}"
) )