{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Analytic disk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create a Magrittetorch 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 magrittetorch.model.model as magritte\n", "import magrittetorch.tools.mesher as mesher \n", "import magrittetorch.tools.setup as setup\n", "import numpy as np # Data structures\n", "from astropy import units, constants # Unit conversions\n", "from astropy.units import Quantity # For function arguments\n", "import os\n", "import torch\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/Magrittetorch-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" ] }, { "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\n", "kb = constants.k_B\n", "m_H2 = 2.01588 * constants.u\n", "\n", "XCO = 6.0e-4 # [.]\n", "vturb = 1.5e+3 * units.m/units.s # [m/s]\n", "\n", "M_star = 2.0 * constants.M_sun\n", "r_star = 2.0 * constants.au\n", "T_star = 2500 * units.K # [K]\n", "\n", "rho_in = 5.0e-12 * units.kg/units.m**3 # [kg/m^3]\n", "r_out = 600.0 * constants.au\n", "r_in = 10.0 * constants.au\n", "T_in = 1500 * units.K # [K]\n", "\n", "p = -2.125 # [.]\n", "h = 1.125 # [.]\n", "q = -0.5 # [.]\n", "\n", "\n", "def cylindrical(x: Quantity[units.m], y: Quantity[units.m], z: Quantity[units.m])-> tuple[Quantity[units.m], Quantity[units.dimensionless_unscaled], Quantity[units.m]]:\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: Quantity[units.m]) -> Quantity[units.m]:\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: [Quantity[units.m], Quantity[units.m], Quantity[units.m]]) -> Quantity[units.kg/units.m**3]:\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 Quantity[units.m**-3]:\n", " \"\"\"\n", " H2 number density function.\n", " \"\"\"\n", " return density(rr) / (2.01588 * constants.u)\n", "\n", "\n", "def abn_nCO(rr: [Quantity[units.m], Quantity[units.m], Quantity[units.m]]) -> Quantity[units.m**-3]:\n", " \"\"\"\n", " CO number density function.\n", " \"\"\"\n", " return XCO * abn_nH2(rr)\n", "\n", "\n", "def temperature(rr: [Quantity[units.m], Quantity[units.m], Quantity[units.m]]) -> Quantity[units.K]:\n", " \"\"\"\n", " Temperature structure.\n", " \"\"\"\n", " # Convert carthesian to cylindrical coords\n", " (rxy,phi,z) = cylindrical(rr[:,0],rr[:,1],rr[:,2])\n", " # Get spherical radial coordinate\n", " r = np.sqrt(rxy**2 + z**2)\n", " # Compute temperature\n", " return T_star * np.power(r/r_star, q)\n", "\n", " \n", "def velocity_f(rr: [Quantity[units.m], Quantity[units.m], Quantity[units.m]]) -> Quantity[units.m/units.s]:\n", " \"\"\"\n", " Velocity structure.\n", " \"\"\"\n", " # Convert carthesian to cylindrical coords\n", " (rxy,phi,z) = cylindrical(rr[:,0],rr[:,1],rr[:,2])\n", " # Get spherical radial coordinate\n", " r = np.sqrt(rxy**2 + z**2)\n", " # Compute angle\n", " d = phi.si.value + 0.5 * np.pi\n", " # Compute velocity (in carthesian)\n", " return np.sqrt(G * M_star/r)[:, np.newaxis] * np.stack([np.cos(d), np.sin(d), np.zeros_like(d)], axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create background mesh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the desired background mesh, a ($75 \\times 75 \\times 75$) cube." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "resolution = 75\n", "\n", "xs = np.linspace(-r_out, +r_out, resolution, endpoint=True)\n", "ys = np.linspace(-r_out, +r_out, resolution, endpoint=True)\n", "zs = np.linspace(-r_out, +r_out, resolution, endpoint=True)\n", "\n", "(Xs, Ys, Zs) = np.meshgrid(xs, ys, zs)\n", "\n", "# Extract positions of points in background mesh\n", "position = np.stack((Xs.ravel(), Ys.ravel(), Zs.ravel()), axis=1)\n", "\n", "# Evaluate the density on the cube\n", "rhos = density(np.stack([Xs, Ys, Zs], axis=1))\n", "rhos_min = np.min(rhos[rhos!=0.0])\n", "rhos += rhos_min" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we remesh this model using the new remesher" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "new interior points: 198688\n", "number boundary points: 1538\n" ] } ], "source": [ "positions_reduced, nb_boundary = mesher.remesh_point_cloud(position, rhos.si.value.ravel(), max_depth=10, threshold= 2e-1, hullorder=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We add a spherical inner boundary at 0.01*r_out" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "number of points in reduced grid: 200525\n", "number of boundary points: 1838\n" ] } ], "source": [ "healpy_order = 5 #using healpy to determine where to place the 12*healpy_order**2 boundary points on the sphere\n", "origin = np.array([0.0,0.0,0.0]).T * units.m #the origin of the inner boundary sphere\n", "positions_reduced, nb_boundary = mesher.point_cloud_add_spherical_inner_boundary(positions_reduced, nb_boundary, 0.01*r_out, healpy_order=healpy_order, origin=origin)\n", "print(\"number of points in reduced grid: \", len(positions_reduced))\n", "print(\"number of boundary points: \", nb_boundary)\n", "#Magritte-torch does not yet support having an inner boundary, so we just treat them as normal points\n", "nb_boundary = 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We add a spherical outer boundary at r_out" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "number of points in reduced grid: 114051\n" ] } ], "source": [ "healpy_order = 15 #using healpy to determine where to place the 12*healpy_order**2 boundary points on the sphere\n", "origin = np.array([0.0,0.0,0.0]).T * units.m #the origin of the inner boundary sphere\n", "positions_reduced, nb_boundary = mesher.point_cloud_add_spherical_outer_boundary(positions_reduced, nb_boundary, r_out, healpy_order=healpy_order, origin=origin)\n", "print(\"number of points in reduced grid: \", len(positions_reduced))" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "npoints = len(positions_reduced)\n", "\n", "# Extract Delaunay vertices (= Voronoi neighbors)\n", "delaunay = Delaunay(positions_reduced)\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", "# Convenience arrays\n", "zeros = np.zeros(npoints)\n", "ones = np.ones (npoints)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convert model functions to arrays based the model mesh." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "position = positions_reduced\n", "velocity = velocity_f (positions_reduced)\n", "nH2 = abn_nH2 (positions_reduced)\n", "nCO = abn_nCO (positions_reduced)\n", "tmp = temperature(positions_reduced)\n", "trb = np.ones(len(tmp))*vturb" ] }, { "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": 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 model to: /lhome/thomasc/Magrittetorch-examples/Analytic_disk/model_analytic_disk.hdf5\n" ] } ], "source": [ "from magrittetorch.model.geometry.geometry import GeometryType\n", "from magrittetorch.model.geometry.boundary import BoundaryType\n", "\n", "model = magritte.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, zeros/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) # Set turbulence velocity\n", "\n", "model = setup.set_Delaunay_neighbor_lists (model) # Automatically computes and sets neighbors for each point, using a Delaunay triangulation\n", "# For unitless quantities, we can also directly set the torch tensors\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", "# 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 quadrature points for NLTE radiative transfer\n", "\n", "# And write the model to disk\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]})" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }