{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3D Phantom"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We create a Magrittetorch 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": [],
"source": [
"import magrittetorch.tools.radiativetransferutils as rtutils\n",
"import magrittetorch.tools.setup as setup\n",
"import torch\n",
"from magrittetorch.model.model import Model # Model class\n",
"from magrittetorch.model.geometry.geometry import GeometryType # GeometryType enum\n",
"from magrittetorch.model.geometry.boundary import BoundaryType # BoundaryType enum\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, units # 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/Magrittetorch-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",
"!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": 7,
"metadata": {},
"outputs": [],
"source": [
"# Loading the data\n",
"setupData = plons.LoadSetup(wdir, \"wind\")\n",
"dumpData = plons.LoadFullDump(dump_file, setupData)\n",
"\n",
"position = dumpData[\"position\"]*units.cm # position vectors [cm]\n",
"velocity = dumpData[\"velocity\"]*units.km/units.s # velocity vectors [km/s]\n",
"rho = dumpData[\"rho\"]*units.g/units.cm**3 # density [g/cm^3]\n",
"u = dumpData[\"u\"] # internal energy density [erg/g]\n",
"tmp = dumpData[\"Tgas\"]*units.K # temperature [K]\n",
"tmp[tmp<2.725*units.K] = 2.725*units.K # Cut-off temperatures below 2.725 K\n",
"\n",
"# Extract the number of points\n",
"npoints = len(rho)\n",
"\n",
"# Convert rho (total density) to abundances\n",
"nH2 = rho * constants.N_A / (2.02 * units.g / units.mol)\n",
"nCO = nH2 * 1.0e-4\n",
"\n",
"# Define turbulence at 150 m/s\n",
"trb = 150.0*units.m/units.s"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 6960643/6960643 [00:32<00:00, 213188.65it/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))])\n",
"#convert to torch tensor\n",
"boundary_torch = torch.from_numpy(boundary)\n",
"nbs_torch = torch.Tensor(nbs).type(torch.int64)\n",
"n_nbs_torch = torch.Tensor(n_nbs).type(torch.int64)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now all data is read, we can use it to construct a Magrittetorch model.\n",
"\n",
"
\n",
"\n",
"Warning\n",
"\n",
"Including all radiative transitions can be computationally expensive (both in time and memory cost) for self-consistent NLTE radiative transfer. For LTE radiative transfer, this is not the case, altough if one wants to image a specific line, that line must be in the list of considered transitions. For these examples, we include the first 10 radiative transitions of CO (J=1-0 to J=10-9). To consider all transitions, use `setup.set_linedata_from_LAMDA_file` as in the commented line below it.\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 model to: /lhome/thomasc/Magrittetorch-examples/Phantom_3D/model_Phantom_3D.hdf5\n"
]
}
],
"source": [
"model = Model(model_file) # Create model object\n",
"model.geometry.geometryType.set(GeometryType.General3D) # This is a 3D model\n",
"model.geometry.boundary.boundaryType.set(BoundaryType.Sphere3D) # With a spherical boundary\n",
"\n",
"# In order to make unit conversions trivial, we use astropy quantities as input\n",
"model.geometry.points.position.set_astropy(position) # Set point positions\n",
"model.geometry.points.velocity.set_astropy(velocity) # Set point velocities\n",
"model.chemistry.species.abundance.set_astropy(np.stack([nCO, nH2, np.zeros(npoints)/units.m**3], axis=1))# Set species number densities\n",
"model.chemistry.species.symbol.set(np.array(['CO', 'H2', 'e-'], dtype='S')) #Set species symbols; should correspond to the LAMDA file format\n",
"#Note: the dtype='S' is necessary to correctly save and read the species symbols to/from the hdf5 file\n",
"\n",
"model.thermodynamics.temperature.gas.set_astropy(tmp) # Set gas temperature\n",
"model.thermodynamics.turbulence.vturb.set_astropy(trb*np.ones(npoints)) # Set turbulence velocity\n",
"\n",
"# Set the neighbors\n",
"model.geometry.points.neighbors.set(nbs_torch) # Set neighbors\n",
"model.geometry.points.n_neighbors.set(n_nbs_torch) # Set number of neighbors\n",
"model.geometry.boundary.boundary2point.set(boundary_torch) # Set which points are boundary points\n",
"\n",
"# # For unitless quantities, we can also directly set the torch tensors\n",
"# nb_boundary = len(boundary)\n",
"# boundary_indices = torch.arange(nb_boundary, dtype = torch.int64)\n",
"# model.geometry.boundary.boundary2point.set(boundary_indices) # Set which points are boundary points\n",
"# model = setup.set_Delaunay_neighbor_lists (model) # Automatically computes and sets neighbors for each point, using a Delaunay triangulation\n",
"\n",
"# Conveniently, the remeshing function puts the boundary points in front of the positions array\n",
"model = setup.set_boundary_condition_CMB (model) # Set CMB as boundary condition\n",
"model = setup.set_uniform_rays (model, 12) # Number of rays for NLTE raytracing; has be of the form 12*2**n\n",
"\n",
"#As this example does not do NLTE, we might as well only consider the first 10 transitions of CO\n",
"model = setup.set_linedata_from_LAMDA_file(model, lamda_file, {'considered transitions': [i for i in range(10)]})\n",
"# model = setup.set_linedata_from_LAMDA_file(model, lamda_file) # Consider all transitions\n",
"model = setup.set_quadrature (model, 7) # Set number of frequency quadrature points for NLTE radiative transfer\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 = position.to_value(units.cm), # yt expects cm not m \n",
" node_data = {('connect1', 'n'): nCO[delaunay.simplices].to_value(units.m**-3)}\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.11.4"
},
"widgets": {
"application/vnd.jupyter.widget-state+json": {
"state": {},
"version_major": 2,
"version_minor": 0
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}