{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Analytic disk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create a Magritte model from an analytic description of a Keplerian disk (see e.g. [Homan et al. 2018](https://doi.org/10.1051/0004-6361/201732246); [Booth et al. 2019](https://doi.org/10.1051/0004-6361/201834388))." ] }, { "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 magritte.setup as setup # Model setup\n", "import magritte.core as magritte # Core functionality\n", "import magritte.mesher as mesher # Mesher\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 units, constants # Unit conversions\n", "from scipy.spatial import Delaunay, cKDTree # 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/Analytic_disk/\"" ] }, { "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": [ "model_file = os.path.join(wdir, 'model_analytic_disk.hdf5') # Resulting Magritte model\n", "lamda_file = os.path.join(wdir, 'co.txt' ) # Line data file\n", "bmesh_name = os.path.join(wdir, 'analytic_disk' ) # bachground mesh name (no extension!)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use a data file that can be downloaded with the following links." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "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 $lamda_link --output-document $lamda_file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The functions below describe the disk structure, based on the Magritte application presented in [De Ceuster et al. (2019)](https://doi.org/10.1093/mnras/stz3557)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "G = constants.G.si.value\n", "kb = constants.k_B.si.value\n", "m_H2 = 2.01588 * constants.u.si.value\n", "\n", "XCO = 6.0e-4 # [.]\n", "vturb = 1.5e+3 # [m/s]\n", "\n", "M_star = 2.0 * constants.M_sun.si.value\n", "r_star = 2.0 * constants.au.si.value\n", "T_star = 2500 # [K]\n", "\n", "rho_in = 5.0e-12 # [kg/m^3]\n", "r_out = 600.0 * constants.au.si.value\n", "r_in = 10.0 * constants.au.si.value\n", "T_in = 1500 # [L]\n", "\n", "p = -2.125 # [.]\n", "h = 1.125 # [.]\n", "q = -0.5 # [.]\n", "\n", "\n", "def cylindrical(x, y, z):\n", " \"\"\"\n", " Convert cartesian to cylindrical coordinates.\n", " \"\"\"\n", " rxy = np.sqrt (x**2 + y**2)\n", " phi = np.arctan2(y, x)\n", " return (rxy, phi, z)\n", "\n", "\n", "def H(rxy):\n", " \"\"\"\n", " Vertical Gaussian scale height.\n", " \"\"\"\n", " # Compute prefactor\n", " prefactor = r_in * np.sqrt(kb * T_in / m_H2 * r_in / (G * M_star))\n", " # Return power law result\n", " return prefactor * np.power(rxy/r_in, h)\n", "\n", "\n", "def density(rr):\n", " \"\"\"\n", " Keplerian disk density in cylindrical coordinates.\n", " \"\"\"\n", " # Convert carthesian to cylindrical coords\n", " (rxy,phi,z) = cylindrical(rr[0],rr[1],rr[2])\n", " # Compute density\n", " rho = rho_in * np.power(rxy/r_in, p) * np.exp(-0.5 * (z/H(rxy))**2)\n", " # Set radii below r_in to zero (in a way that also work for np arrays)\n", " if hasattr(rxy, \"__len__\"):\n", " rho[rxy\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": 14, "metadata": { "tags": [] }, "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 (min. 5)\n", "model.parameters.set_nlspecs (1) # Number of line species\n", "model.parameters.set_nquads (51) # 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(nb_boundary)\n", "model.geometry.boundary.boundary2point.set(np.arange(nb_boundary))\n", "\n", "direction = np.array([[+1,0,0], [-1,0,0]]) # 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": 15, "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 and x-axis." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))\n", "sl.set_cmap (('connect1', 'n'), 'magma')\n", "sl.zoom (1.2)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sl = yt.SlicePlot (ds, 'x', ('connect1', 'n'))\n", "sl.set_cmap (('connect1', 'n'), 'magma')\n", "sl.zoom (1.2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show meshes on the plots." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))\n", "sl.set_cmap (('connect1', 'n'), 'magma')\n", "sl.zoom (1.2)\n", "sl.annotate_mesh_lines (plot_args={'color':'lightgrey', 'linewidths':[.25]})" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sl = yt.SlicePlot (ds, 'x', ('connect1', 'n'))\n", "sl.set_cmap (('connect1', 'n'), 'magma')\n", "sl.zoom (1.2)\n", "sl.annotate_mesh_lines (plot_args={'color':'lightgrey', 'linewidths':[.25]})" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.13" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }