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

[galactica] Improve project description + units

Browse Source
This commit is contained in:
Noe Brucy
2020-12-15 15:45:17 +01:00
parent 0b180c0ac6
commit 9864ba97b8
1 changed files with 232 additions and 170 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fragdisk_galactica.py
+232 -170
View File
@@ -27,42 +27,20 @@
# #
from __future__ import print_function, unicode_literals from __future__ import print_function, unicode_literals
import datetime
import os import os
import h5py import h5py
import numpy as N import numpy as np
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 Datafile, PlotInfo, PlotType from astrophysix.simdm.datafiles import PlotInfo, PlotType
from astrophysix.simdm.experiment import (
AppliedAlgorithm,
ParameterSetting,
ParameterVisibility,
ResolvedPhysicalProcess,
Simulation,
)
from astrophysix.simdm.protocol import (
Algorithm,
AlgoType,
InputParameter,
PhysicalProcess,
Physics,
SimulationCode,
)
from astrophysix.simdm.results import GenericResult, Snapshot
from astrophysix.utils.file import FileType from astrophysix.utils.file import FileType
from matplotlib import gridspec
from mpl_toolkits.axes_grid1 import AxesGrid, Grid
from plotter import *
P.rcParams["text.usetex"] = False
from plotter import Plotter, default_params
from ramses_astrophysix import ramses
pp_params = default_params() pp_params = default_params()
pp_params.input.nml_filename = "disk.nml" pp_params.input.nml_filename = "disk.nml"
# pp_params.out.interactive = True
pp_params.pymses.map_size = 2048 pp_params.pymses.map_size = 2048
pp_params.pymses.zoom = 4 pp_params.pymses.zoom = 4
pp_params.pymses.filter = False pp_params.pymses.filter = False
@@ -77,29 +55,177 @@ 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"
nml_key = "cloud_params/beta_cool" nml_key = "cloud_params/beta_cool"
# --- Runs -----
pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}" # ---------------------------------------------- Project creation ------------------------------------------ #
pp_params.astrophysix.descr_fmt = """ # Available project categories are :
<p>Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.</p> # - ProjectCategory.SolarMHD
# - ProjectCategory.PlanetaryAtmospheres
# - ProjectCategory.StarPlanetInteractions
# - ProjectCategory.StarFormation
# - ProjectCategory.Supernovae
# - ProjectCategory.GalaxyFormation
# - ProjectCategory.GalaxyMergers
# - ProjectCategory.Cosmology
proj = Project(
category=ProjectCategory.StarPlanetInteractions,
project_title="Fragdisk",
alias="FRAGDISK",
short_description="Fragmentation of self-gravitating disks",
general_description="""
<h4> Study of the fragmentation of self-gravitating disks</h4>
<p>
See Brucy & Hennebelle 2021 (submitted) for more details.
This database is currently being completed.
</p>
<h4>Abstract:</h4>
<p> <p>
The number behind jr gives the maximal level of refinement. Self-gravitating disks are believed to play an important role in astrophysics in particular regarding the
An ORP is an Outer Rotation Period and corresponds ro the time took by the gas at the border of the disk to rotates around the central object. star and planet formation process.
All quantities are in code units. In this context, disks subject to an idealized cooling process, characterized by a cooling timescale
$\\beta$ expressed in unit of orbital timescale, have been extensively studied.
We take advantage of the Riemann solver and the 3D Godunov scheme implemented in the code Ramses to
perform high resolution simulations, complementing previous studies that have used Smoothed Particle
Hydrodynamics (SPH) or 2D grid codes.
</p> </p>
"""
pl_orp = Plotter( <h4> Description of the simulations: </h4>
<p>
We simulate a disk of gas undergoing purely hydrodynamics forces,
its own gravity and the $\\beta$-cooling.
The simulation is ran with the 3D-grid code Ramses (Teyssier 2002) which uses a Godunov scheme.
The flux between each cell is computed with the HLLC Riemann solver.
The gravity potential is updated at each timestep with a Poisson solver,
and a source term is added to the energy equation to implement the $\\beta$-cooling.
</p>
<p>
The $\\beta$-cooling consists in removing internal energy from the gas with a cooling time:
$t_\\text{cool} = \\beta \\Omega^{-1}$ with $\\Omega$ the rotation frequency.
</p>
<p>
We use the same initial conditions as in Meru & Bate (2012) to allow comparison.
The specific disk setup for Ramses was inspired by Hennebelle et al. (2017).
The disk is initially close to equilibrium with an initial column density profile
$\\Sigma \\propto r^{-1}$ and a temperature profile $T \\propto r^{-1/2}$
where $r$ is the cylindrical radius.
The disk has a radius $r_d = 0.25$ (code units), after which the density is divided by 100.
The density and temperature at the disk radius $r_d$ are chosen so that the mass
of the disk $M_d = 0.1 M_\\star$, where $M_\\star$ is the mass of the central object,
and the initial value of the Toomre parameter at the disk radius is $Q_{0,d} = 2$.
The adiabatic index of the gas is $\\gamma = 5/3$.
</p>
<p>
The simulation is run within a cube of size $L=2$. Although the problem has a cylindrical symmetry,
we use Cartesian coordinates.
This prevents having a singularity at the centre of the box.
One caveat is the poor resolution on the centre of the cube but this is mitigated
by the use of the adaptive mesh refinement (AMR).
Another caveat is that having a cubic box may introduce spurious
reflection at the border of the simulation.
To avoid this, we maintain a dead zone over a radius of $0.875$ (in code units)
where all variables are replaced by their initial value at each timestep.
This method has been used in Hennebelle et al. (2017) and has proven to be efficient.
</p>
<p>
The simulations presented here were run for several values of $\beta$ and several resolutions.
To reduce the computation time, we use the Ramses's Adaptative Mesh Refinement (AMR).
Only the parts of the simulation which are prone to form fragments are simulated with full resolution.
Each cell is refined until the Jeans's length is covered by at least 20 cells or it reaches the maximum
level of refinement $l_{\\max}$.
The level of refinement of a cell is the number of times the simulation box must be divided in eight
equal part to get the cell.
Thus, the resolution of a simulation is given by the value of $l_{\\max}$.
A first set of simulations with $l_{\\max} = 11$ to $l_{\\max} = 12$ are run until
about 5 Outer Rotation Periods (ORP),
that is that the gas at the border of the disk had 5 orbits around the star.
A second set of simulations, labelled tic, for Turbulent Initial Condition,
were run from relaxed initial conditions for $l_{\\max} = 12$ and $l_{\\max} = 13$.
More precisely, they were restarted from a simulation at $\beta = 20$ and $l_{\\max} = 12$
for which the whole disk reached a gravito-turbulent state (after 2 ORPs).
According to Paardekooper et al. (2011) and Clarke et al.(2007), departing
from such turbulent condition should reduce spurious fragmentation.
</p>
""",
data_description="""
<p>The data available for this project is the underlying data of the article Brucy & Hennebelle 2021.</p>
<p>The data is not already fully uploaded. 3D datacube extraction on demand is planned</p>
""",
directory_path="~nbrucy/simus/fragdisk",
)
print(proj)
# ------------------------------------------------------------------------------------------------------ #
# ---- Units ----
pl_units = Plotter(
in_dir, in_dir,
filter_name="104_beta4_jr13", filter_name="104_beta4_jr13",
in_nums="last", in_nums="last",
path_out=out_dir, path_out=out_dir,
pp_params=pp_params, pp_params=pp_params,
) )
orp = cst.Unit.create_unit("ORP", base_unit=pl_orp.comp.info["unit_time"] * 0.79) info = pl_units.comp.info
rd = U.Unit.create_unit(
"r_d",
base_unit=0.25 * info["unit_length"] / info["boxlen"],
)
rhod = U.Unit.create_unit(
"rho_d",
latex="\rho_d",
base_unit=6.36 * info["unit_density"],
)
rho_u = U.Unit.create_unit(
"M_*/r_d^3",
latex="M_\\star.r_d^{-3}",
base_unit=info["unit_mass"] * rd ** (-3),
)
Pd = U.Unit.create_unit(
"P_d",
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(
"(k_kd * M_*)/r_d^3",
latex="v_{k,d}^2.M_\\star.r_d^{-3}",
base_unit=vkd * rho_u,
)
Sigmad = U.Unit.create_unit(
"Sigma_d",
latex="\\Sigma_d",
base_unit=0.25 * info["unit_density"] * info["unit_length"] / info["boxlen"],
)
Sigma_u = U.Unit.create_unit(
"M*/r_d^2",
latex="M_\\star.r_d^{-2}",
base_unit=rd ** (-2) * info["unit_mass"],
)
# ---- Runs -----
# JR13_TIC # JR13_TIC
pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}"
pp_params.astrophysix.descr_fmt = """
<p>Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.</p>
"""
runs = "*_jr13" runs = "*_jr13"
pl_jr13 = Plotter( pl_jr13 = Plotter(
in_dir, in_dir,
@@ -130,12 +256,6 @@ pl_jr12_tic = Plotter(
pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}" pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}"
pp_params.astrophysix.descr_fmt = """ pp_params.astrophysix.descr_fmt = """
<p>Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}</p> <p>Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}</p>
<p>
The number behind jr gives the maximal level of refinement.
An ORP is an Outer Rotation Period and corresponds ro the time took by the gas at the border of the disk to rotates around the central object.
All quantities are in code units.
</p>
""" """
print("JR12_TIC defined") print("JR12_TIC defined")
@@ -174,62 +294,6 @@ print("JR11 defined")
pls = [pl_l11, pl_jr12, pl_jr12_tic, pl_jr13] pls = [pl_l11, pl_jr12, pl_jr12_tic, pl_jr13]
# ----------------------------------------------- Project creation --------------------------------------------------- #
# Available project categories are :
# - ProjectCategory.SolarMHD
# - ProjectCategory.PlanetaryAtmospheres
# - ProjectCategory.StarPlanetInteractions
# - ProjectCategory.StarFormation
# - ProjectCategory.Supernovae
# - ProjectCategory.GalaxyFormation
# - ProjectCategory.GalaxyMergers
# - ProjectCategory.Cosmology
proj = Project(
category=ProjectCategory.StarPlanetInteractions,
project_title="Fragdisk",
alias="FRAGDISK",
short_description="Fragmentation of self-gravitating disks",
general_description="""
<p>Study of the fragmentation of self-gravitating disks. See Brucy & Hennebelle 2021 (submitted) for more details.
This database is currently being completed.
</p>
<h4>Abstract:</h4>
<p>
Self-gravitating disks are believed to play an important role in astrophysics in particular regarding the star and planet formation process. In this context, disks subject to an idealized cooling process, characterized by a cooling timescale $\\beta$ expressed in unit of orbital timescale, have been extensively studied. We take advantage of the Riemann solver and the 3D Godunov scheme implemented in the code Ramses to perform high resolution simulations, complementing previous studies that have used Smoothed Particle Hydrodynamics (SPH) or 2D grid codes.
</p>
<h4> Description of the simulations: </h4>
<p>
We simulate a disk of gas undergoing purely hydrodynamics forces, its own gravity and the $\\beta$-cooling. The simulation is ran with the 3D-grid code Ramses (Teyssier 2002) which uses a Godunov scheme. The flux between each cell is computed with the HLLC Riemann solver. The gravity potential is updated at each timestep with a Poisson solver, and a source term is added to the energy equation to implement the $\\beta$-cooling.
</p>
<p>
The $\\beta$-cooling consists in removing internal energy from the gas with a cooling time:
$t_\\text{cool} = \\beta \Omega^{-1}$ with $\Omega$ the rotation frequency.
</p>
<p>
We use the same initial conditions as in Meru & Bate (2012) to allow comparison.
The specific disk setup for Ramses was inspired by Hennebelle et al. (2017).
The disk is initially close to equilibrium with an initial column density profile $\Sigma \propto r^{-1}$ and a temperature profile $T \propto r^{-1/2}$ where $r$ is the cylindrical radius. The disk has a radius $r_d = 0.25$ (code units), after which the density is divided by 100. The density and temperature at the disk radius $r_d$ are chosen so that the mass of the disk $M_d = 0.1 M_\star$, where $M_\star$ is the mass of the central object, and the initial value of the Toomre parameter at the disk radius is $Q_{0,d} = 2$. The adiabatic index of the gas is $\gamma = 5/3$.
</p>
<p>
The simulation is run within a cube of size $L=2$. Although the problem has a cylindrical symmetry, we use Cartesian coordinates. This prevents having a singularity at the centre of the box. One caveat is the poor resolution on the centre of the cube but this is mitigated by the use of the adaptive mesh refinement (AMR).
Another caveat is that having a cubic box may introduce spurious reflection at the border of the simulation.
To avoid this, we maintain a dead zone over a radius of $0.875$ (in code units) where all variables are replaced by their initial value at each timestep. This method has been used in Hennebelle et al. (2017) and has proven to be efficient.
</p>
""",
data_description="""The data available for this project is the underlying data of the article Brucy & Hennebelle 2021.
The data is not already fully uploaded. 3D datacube extraction on demand is planned""",
directory_path="~nbrucy/simus/fragdisk",
)
print(proj)
# -------------------------------------------------------------------------------------------------------------------- #
# -------------------------------------------------------------------------------------------------------------------- # # -------------------------------------------------------------------------------------------------------------------- #
redo = True redo = True
for pl in pls: for pl in pls:
@@ -252,71 +316,54 @@ for pl in pls:
pl.rules["slice_P"].description = "Pressure slice" pl.rules["slice_P"].description = "Pressure slice"
pl.rules[ pl.rules[
"pdf_coldens" "pdf_coldens"
].description = "Probability function of the logarithm of the column density fluctuations $\sigma = \Sigma/\overline{\Sigma}$ with respect to its azimuthal average" ].description = """
Probability function of the logarithm of the column density fluctuations
$\\sigma = \\Sigma/\\overline{\\Sigma}$ with respect to its azimuthal average
"""
pl.coldens( coldens_kwargs = {
"z", "overwrite": redo,
overwrite=redo, "overwrite_dep": False,
overwrite_dep=False, "unit_space": rd,
unit_space=cst.cm, "center_space": True,
unit_time=orp, "unit_time": orp,
nml_key="cloud_params/beta_cool", "unit": Sigma_u,
vmin=1e-2, "nml_key": "cloud_params/beta_cool",
vmax=1e2, "select": select,
put_units=False, "label": r"$\Sigma$",
select=select, }
label=r"$\Sigma$ [code_units]",
)
pl.coldens( rho_kwargs = {
"y", "overwrite": redo,
overwrite=redo, "overwrite_dep": False,
overwrite_dep=False, "unit_space": rd,
unit_space=cst.cm, "center_space": True,
unit_time=orp, "unit_time": orp,
nml_key="cloud_params/beta_cool", "unit": rho_u,
vmin=1e-2, "nml_key": "cloud_params/beta_cool",
vmax=1e2, "select": select,
put_units=False, "label": r"$\rho$",
select=select, }
label=r"$\Sigma$ [code_units]",
)
pl.slice_rho( P_kwargs = {
"z", "overwrite": redo,
overwrite=redo, "overwrite_dep": False,
overwrite_dep=False, "unit_space": rd,
unit_space=cst.cm, "center_space": True,
unit_time=orp, "unit_time": orp,
nml_key="cloud_params/beta_cool", "unit": P_u,
put_units=False, "nml_key": "cloud_params/beta_cool",
select=select, "select": select,
label=r"$\rho$ [code_units]", "label": r"$P$",
) }
pl.slice_rho( pl.coldens("z", **coldens_kwargs)
"y", pl.coldens("y", **coldens_kwargs)
overwrite=redo,
overwrite_dep=False,
unit_space=cst.cm,
unit_time=orp,
nml_key="cloud_params/beta_cool",
put_units=False,
select=select,
label=r"$\rho$ [code_units]",
)
pl.slice_P( pl.slice_rho("z", **rho_kwargs)
"z", pl.slice_rho("y", **rho_kwargs)
overwrite=redo,
overwrite_dep=False, pl.slice_P("z", **P_kwargs)
unit_space=cst.cm,
unit_time=orp,
nml_key="cloud_params/beta_cool",
put_units=False,
select=select,
label=r"$P$ [code_units]",
)
pl.pdf_coldens( pl.pdf_coldens(
"z", "z",
@@ -331,13 +378,13 @@ for pl in pls:
) )
# -------------------------------------------------------------------------------------------------------------------- # # ----------------------------------------------------------------------------------------------------------------- #
# Light fake PlotInfo to replace real one to reduce size of the hdf5 file
pi = PlotInfo( pi = PlotInfo(
plot_type=PlotType.LINE_PLOT, plot_type=PlotType.LINE_PLOT,
xaxis_values=N.array([10.0, 20.0, 30.0, 40.0, 50.0]), xaxis_values=np.array([10.0, 20.0, 30.0, 40.0, 50.0]),
yaxis_values=N.array([1.256, 2.456, 3.921, 4.327, 5.159]), yaxis_values=np.array([1.256, 2.456, 3.921, 4.327, 5.159]),
xaxis_log_scale=False, xaxis_log_scale=False,
yaxis_log_scale=False, yaxis_log_scale=False,
xaxis_label="Mass", xaxis_label="Mass",
@@ -347,9 +394,26 @@ pi = PlotInfo(
yaxis_unit=U.Mpc, yaxis_unit=U.Mpc,
) )
# Convert plot_info to HDF5 files
for simu in proj.simulations: for simu in proj.simulations:
for snap in simu.snapshots: for snap in simu.snapshots:
# Change unit to things Galactica understands
snap.time = (snap.time[0], U.none)
snap.description = """
<p>Snapshot at {time} ORPs.</p>
<h5> Notations </h5>
<p>
$r_d$ : radius of the disk <br />
ORP : outer rotation period, that is period of the gas at $r_d$ <br />
$M_\\star$ : Mass of the central object <br />
$v_{kepl}$ : Keplerian speed at $r_d$
</p>
""".format(
time=snap.time[0], kepl="{k,d}"
)
# Convert plotinfo into HDF5 file
for df in snap.datafiles: for df in snap.datafiles:
name = df[FileType.JPEG_FILE].filename name = df[FileType.JPEG_FILE].filename
name = os.path.splitext(name)[0] + ".h5" name = os.path.splitext(name)[0] + ".h5"
@@ -360,15 +424,13 @@ for simu in proj.simulations:
df[FileType.HDF5_FILE] = out_dir + "/" + name df[FileType.HDF5_FILE] = out_dir + "/" + name
df.plot_info = pi df.plot_info = pi
# Trim parameters name
for param in ramses.input_parameters: for param in ramses.input_parameters:
param.key = os.path.basename(param.key) param.key = os.path.basename(param.key)
study = SimulationStudy(project=proj) # Validity check
for sim in study.project.simulations:
for snap in sim.snapshots:
snap.time = (snap.time[0], cst.year)
proj.galactica_validity_check() proj.galactica_validity_check()
# Save
study = SimulationStudy(project=proj)
study.save_HDF5(out_dir + "/fragdisk_study.h5") study.save_HDF5(out_dir + "/fragdisk_study.h5")