{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3D Phantom"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We create a Magritte model from a snapshot of 3D Phantom hydrodynamics simulation.\n",
"The hydro model was kindly provided by Jolien Malfait.\n",
"The Phantom binary files are used directly, and are loaded using the [plons](https://plons.readthedocs.io) package."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import the required functionalty."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"hwloc/linux: Ignoring PCI device with non-16bit domain.\n",
"Pass --enable-32bits-pci-domain to configure to support such devices\n",
"(warning: it would break the library ABI, don't enable unless really needed).\n"
]
}
],
"source": [
"import magritte.setup as setup # Model setup\n",
"import magritte.core as magritte # Core functionality\n",
"import plons # Loading phantom data\n",
"import numpy as np # Data structures\n",
"import warnings # Hide warnings\n",
"warnings.filterwarnings('ignore') # especially for yt\n",
"import yt # 3D plotting\n",
"import os\n",
"\n",
"from tqdm import tqdm # Progress bars\n",
"from astropy import constants # Unit conversions\n",
"from scipy.spatial import Delaunay # Finding neighbors\n",
"from yt.funcs import mylog # To avoid yt output \n",
"mylog.setLevel(40) # as error messages"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define a working directory (you will have to change this; it must be an **absolute path**)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"wdir = \"/lhome/thomasc/Magritte-examples/Phantom_3D/\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create the working directory."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"!mkdir -p $wdir"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define file names."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"dump_file = os.path.join(wdir, 'model_Phantom_3D' ) # Phantom full dump (snapshot)\n",
"setup_file = os.path.join(wdir, 'wind.setup' ) # Phantom setup file\n",
"input_file = os.path.join(wdir, 'wind.in' ) # Phantom input file\n",
"model_file = os.path.join(wdir, 'model_Phantom_3D.hdf5') # Resulting Magritte model\n",
"lamda_file = os.path.join(wdir, 'co.txt' ) # Line data file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use a snapshot and data file that can be downloaded with the following links."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"dump_link = \"https://github.com/Ensor-code/phantom-models/raw/main/Malfait+2021/v05e50/wind_v05e50\"\n",
"setup_link = \"https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2021/v05e50/wind.setup\"\n",
"input_link = \"https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2021/v05e50/wind.in\"\n",
"lamda_link = \"https://home.strw.leidenuniv.nl/~moldata/datafiles/co.dat\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Dowload the snapshot and the linedata (``%%capture`` is just used to suppress the output)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"%%capture\n",
"print(setup_file)\n",
"!wget $dump_link --output-document $dump_file\n",
"!wget $setup_link --output-document $setup_file\n",
"!wget $input_link --output-document $input_file\n",
"!wget $lamda_link --output-document $lamda_file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extract data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The script below extracts the required data from the snapshot `phantom dump` file."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Loading the data\n",
"setupData = plons.LoadSetup(wdir, \"wind\")\n",
"dumpData = plons.LoadFullDump(dump_file, setupData)\n",
"\n",
"pos_x = dumpData[\"x\"]*1e-2 # x position vectors [cm -> m]\n",
"pos_y = dumpData[\"y\"]*1e-2 # y position vectors [cm -> m]\n",
"pos_z = dumpData[\"z\"]*1e-2 # z position vectors [cm -> m]\n",
"position = np.array((pos_x, pos_y, pos_z)).T\n",
"vel_x = dumpData[\"vx\"]*1e3 # x velocity vectors [km/s -> m/s]\n",
"vel_y = dumpData[\"vy\"]*1e3 # y velocity vectors [km/s -> m/s]\n",
"vel_z = dumpData[\"vz\"]*1e3 # z velocity vectors [km/s -> m/s]\n",
"velocity = np.array((vel_x, vel_y, vel_z)).T\n",
"velocity = velocity/constants.c.si.value # velocity vectors [m/s -> 1/c]\n",
"rho = dumpData[\"rho\"] # density [g/cm^3]\n",
"u = dumpData[\"u\"] # internal energy density [erg/g]\n",
"tmp = dumpData[\"temp\"] # temperature [K]\n",
"tmp[tmp<2.725] = 2.725 # Cut-off temperatures below 2.725 K\n",
"\n",
"# Extract the number of points\n",
"npoints = len(rho)\n",
"\n",
"# Convenience arrays\n",
"zeros = np.zeros(npoints)\n",
"ones = np.ones (npoints)\n",
"\n",
"# Convert rho (total density) to abundances\n",
"nH2 = rho * 1.0e+6 * constants.N_A.si.value / 2.02\n",
"nCO = nH2 * 1.0e-4\n",
"\n",
"# Define turbulence at 150 m/s\n",
"trb = (150.0/constants.c.si.value)**2 * ones"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 6960643/6960643 [00:38<00:00, 178604.41it/s]\n"
]
}
],
"source": [
"# Extract Delaunay vertices (= Voronoi neighbors)\n",
"delaunay = Delaunay(position)\n",
"(indptr, indices) = delaunay.vertex_neighbor_vertices\n",
"neighbors = [indices[indptr[k]:indptr[k+1]] for k in range(npoints)]\n",
"nbs = [n for sublist in neighbors for n in sublist]\n",
"n_nbs = [len(sublist) for sublist in neighbors]\n",
"\n",
"# Compute the indices of the boundary particles of the mesh, extracted from the Delaunay vertices\n",
"boundary = set([])\n",
"for i in tqdm(range(delaunay.neighbors.shape[0])):\n",
" for k in range(4):\n",
" if (delaunay.neighbors[i][k] == -1):\n",
" nk1,nk2,nk3 = (k+1)%4, (k+2)%4, (k+3)%4 \n",
" boundary.add(delaunay.simplices[i][nk1])\n",
" boundary.add(delaunay.simplices[i][nk2])\n",
" boundary.add(delaunay.simplices[i][nk3])\n",
" \n",
"boundary = list(boundary)\n",
"boundary = np.array(boundary)\n",
"\n",
"# The above calculation turned out to be unsatisfactory.\n",
"# Since the outer boundary is assumed to be a sphere,\n",
"# we add all points which fall inside the boundary defined above.\n",
"b_nms = np.linalg.norm(position[boundary], axis=1)\n",
"p_nms = np.linalg.norm(position, axis=1)\n",
"boundary = np.array([i[0] for i in np.argwhere(p_nms >= np.min(b_nms))])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now all data is read, we can use it to construct a Magritte model.\n",
"\n",
"
\n",
"\n",
"Warning\n",
"\n",
"Since we do not aim to do self-consistent non-LTE simulations in these examples, we only consider the first radiative transition of CO (J=1-0). To consider all transitions, use `setup.set_linedata_from_LAMDA_file` as in the commented line below it.\n",
"We also only consider 2 rays here (up and down the direction we want to image). To consider all directions, comment out the indecated lines and use `setup.set_uniform_rays` as in the commented line below.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Not considering all radiative transitions on the data file but only the specified ones!\n",
"Writing parameters...\n",
"Writing points...\n",
"Writing rays...\n",
"Writing boundary...\n",
"Writing chemistry...\n",
"Writing species...\n",
"Writing thermodynamics...\n",
"Writing temperature...\n",
"Writing turbulence...\n",
"Writing lines...\n",
"Writing lineProducingSpecies...\n",
"Writing linedata...\n",
"ncolpoar = 2\n",
"--- colpoar = 0\n",
"Writing collisionPartner...\n",
"(l, c) = 0, 0\n",
"--- colpoar = 1\n",
"Writing collisionPartner...\n",
"(l, c) = 0, 1\n",
"Writing quadrature...\n",
"Writing populations...\n",
"Writing radiation...\n",
"Writing frequencies...\n"
]
}
],
"source": [
"model = magritte.Model () # Create model object\n",
"\n",
"model.parameters.set_model_name (model_file) # Magritte model file\n",
"model.parameters.set_dimension (3) # This is a 3D model\n",
"model.parameters.set_npoints (npoints) # Number of points\n",
"model.parameters.set_nrays (2) # Number of rays \n",
"model.parameters.set_nspecs (3) # Number of species\n",
"model.parameters.set_nlspecs (1) # Number of line species\n",
"model.parameters.set_nquads (31) # Number of quadrature points\n",
"\n",
"model.geometry.points.position.set(position)\n",
"model.geometry.points.velocity.set(velocity)\n",
"\n",
"model.geometry.points. neighbors.set( nbs)\n",
"model.geometry.points.n_neighbors.set(n_nbs)\n",
"\n",
"model.chemistry.species.abundance = np.array((nCO, nH2, zeros)).T\n",
"model.chemistry.species.symbol = ['CO', 'H2', 'e-']\n",
"\n",
"model.thermodynamics.temperature.gas .set(tmp)\n",
"model.thermodynamics.turbulence.vturb2.set(trb)\n",
"\n",
"model.parameters.set_nboundary(boundary.shape[0])\n",
"model.geometry.boundary.boundary2point.set(boundary)\n",
"\n",
"direction = np.array([[0,0,+1], [0,0,-1]]) # Comment out to use all directions\n",
"model.geometry.rays.direction.set(direction) # Comment out to use all directions\n",
"model.geometry.rays.weight .set(0.5 * np.ones(2)) # Comment out to use all directions\n",
"\n",
"# setup.set_uniform_rays (model) # Uncomment to use all directions\n",
"setup.set_boundary_condition_CMB (model)\n",
"setup.set_linedata_from_LAMDA_file(model, lamda_file, {'considered transitions': [0]})\n",
"# setup.set_linedata_from_LAMDA_file(model, lamda_file) # Consider all transitions\n",
"setup.set_quadrature (model)\n",
"\n",
"model.write()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Load the data in a `yt` unstructured mesh."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"ds = yt.load_unstructured_mesh(\n",
" connectivity = delaunay.simplices.astype(np.int64),\n",
" coordinates = delaunay.points.astype(np.float64) * 100.0, # yt expects cm not m \n",
" node_data = {('connect1', 'n'): nCO[delaunay.simplices].astype(np.float64)}\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot a slice through the mesh orthogonal to the z-axis."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))\n",
"sl.set_cmap (('connect1', 'n'), 'magma')\n",
"sl.zoom (1.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Show mesh on the plot."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))\n",
"sl.set_cmap (('connect1', 'n'), 'magma')\n",
"sl.zoom (1.1)\n",
"sl.annotate_mesh_lines (plot_args={'color':'lightgrey', 'linewidths':[.25]})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the [next example](4_reduce_Phantom_3D.ipynb) we demonstrate how to reduce this model as in [De Ceuster et al. (2020)](https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.5194D/abstract)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.14"
},
"widgets": {
"application/vnd.jupyter.widget-state+json": {
"state": {},
"version_major": 2,
"version_minor": 0
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}