{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analytic spiral"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We create a Magritte model from an analytic description of an Archimedean spiral outflow.\n",
"The description of the model is based on [Homan et al. (2015)](https://www.aanda.org/articles/aa/full_html/2015/07/aa25933-15/aa25933-15.html)."
]
},
{
"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_spiral/\""
]
},
{
"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_spiral.hdf5') # Resulting Magritte model\n",
"lamda_file = os.path.join(wdir, 'co.txt' ) # Line data file\n",
"bmesh_name = os.path.join(wdir, 'analytic_spiral' ) # 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": [
"## Spiral parameters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The functions below describe the spiral structure, based on the parameters in [Homan et al. (2015)](https://www.aanda.org/articles/aa/full_html/2015/07/aa25933-15/aa25933-15.html)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"phi_0 = 0 / 180.0 * np.pi # [radian]\n",
"alpha = 15 / 180.0 * np.pi # [radian]\n",
"sigma = 20 * constants.au.si.value # [m]\n",
"b = 270 / (2.0*np.pi) * constants.au.si.value # [m]\n",
"\n",
"\n",
"def r_0(r, theta, phi):\n",
" return b*(phi - phi_0)\n",
"\n",
"\n",
"def sigma_r(r, theta, phi):\n",
" return sigma\n",
"\n",
"\n",
"def theta_0(r, theta, phi):\n",
" return 0.5*np.pi\n",
"\n",
"\n",
"def sigma_theta(r, theta, phi):\n",
" return alpha\n",
"\n",
"\n",
"def S(r, theta, phi):\n",
" \"\"\"\n",
" Spiral-shaped distribution function.\n",
" \"\"\"\n",
" rd = (r - r_0 (r,theta,phi)) / sigma_r (r,theta,phi)\n",
" td = (theta - theta_0(r,theta,phi)) / sigma_theta(r,theta,phi)\n",
" return np.exp(-0.5*(rd**2 + td**2))\n",
"\n",
"\n",
"def SS(r, theta, phi):\n",
" \"\"\"\n",
" Spiral-shaped distribution function.\n",
" (Including 10 windings.)\n",
" \"\"\"\n",
" result = 0.0\n",
" for i in range(10):\n",
" result += S(r, theta, phi+2.0*np.pi*i)\n",
" return result"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Physical parameters"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"R_max = 600.0 * constants.au.si.value # [m]\n",
"R_star = 1.0 * constants.R_sun.si.value # [m]\n",
"v_inf = 14.5e+3 # [m/s]\n",
"beta = 0.4 # [.]\n",
"Mdot = 1.5e-6 * (constants.M_sun/units.year).si.value # [kg/s]\n",
"R_0 = 2.0 * R_star # [m]\n",
"epsilon = 0.5 # [.]\n",
"T_star = 2330.0 # [K]\n",
"T_0 = 60.0 # [K]\n",
"rho_max = 1.0e-11 # [kg/m^3]\n",
"\n",
"\n",
"def v_r(r):\n",
" \"\"\"\n",
" Beta-law velocity profile.\n",
" \"\"\"\n",
" return v_inf * (1.0 - R_0/r)**beta\n",
"\n",
"\n",
"def rho_regular(r):\n",
" \"\"\"\n",
" Density profile assuming a regular, spherically symmetric outflow.\n",
" \"\"\"\n",
" return Mdot / (4.0*np.pi*r**2*v_r(r))\n",
"\n",
"\n",
"def rho_spiral(r, theta, phi):\n",
" \"\"\"\n",
" Density profile assuming a spiral outflow.\n",
" \"\"\"\n",
" return rho_max * (2.0*np.pi*b/r)**2 * SS(r, theta, phi)\n",
"\n",
"\n",
"def rho_wind(r, theta, phi):\n",
" \"\"\"\n",
" Wind density, combined regular and spiral outflow.\n",
" \"\"\"\n",
" return rho_regular(r) + rho_spiral(r, theta, phi)\n",
"\n",
"\n",
"def T_regular(r):\n",
" \"\"\"\n",
" Temperature profile assuming a regular, spherically symmetric outflow.\n",
" \"\"\"\n",
" return T_star * (R_star/r)**epsilon\n",
"\n",
"\n",
"def T_spiral(r, theta, phi):\n",
" \"\"\"\n",
" Temperature profile assuming a spiral outflow.\n",
" \"\"\"\n",
" return T_0 * SS(r, theta, phi)\n",
"\n",
"\n",
"def T_wind(r, theta, phi):\n",
" \"\"\"\n",
" Wind temperature, combined regular and spiral outflow.\n",
" \"\"\"\n",
" return T_regular(r) + T_spiral(r, theta, phi)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def spherical(x,y,z):\n",
" \"\"\"\n",
" Convert cartesian to spherical coordinates.\n",
" \"\"\"\n",
" r = np.sqrt (x**2 + y**2 + z**2)\n",
" t = np.arccos (z/r)\n",
" p = np.arctan2(y,x)\n",
" return (r,t,p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create background mesh (with desired mesh element sizes)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define the desired mesh element size in a ($100 \\times 100 \\times 100$) cube."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"resolution = 100\n",
"\n",
"xs = np.linspace(-R_max, +R_max, resolution, endpoint=True)\n",
"ys = np.linspace(-R_max, +R_max, resolution, endpoint=True)\n",
"zs = np.linspace(-R_max, +R_max, resolution, endpoint=True)\n",
"\n",
"(Xs, Ys, Zs) = np.meshgrid(xs, ys, zs)\n",
"\n",
"position = np.array((Xs.ravel(), Ys.ravel(), Zs.ravel())).T\n",
"\n",
"# Compute the corresponding spherical coordinates\n",
"Rs = np.sqrt (Xs**2 + Ys**2 + Zs**2)\n",
"Ts = np.arccos (Zs/Rs)\n",
"Ps = np.arctan2(Ys,Xs)\n",
"\n",
"# Evaluate the density on the cube\n",
"rhos = rho_wind(Rs, Ts, Ps)\n",
"rhos = np.nan_to_num(rhos, nan=1.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we remesh the grid using the built-in remesher"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"new interior points: 583111\n",
"number boundary points: 1538\n"
]
}
],
"source": [
"positions_reduced, nb_boundary = mesher.remesh_point_cloud(position, rhos.ravel(), max_depth=10, threshold= 2e-1, hullorder=4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We add a spherical inner boundary at 0.04*R_max"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"number of points in reduced grid: 584918\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 #the origin of the inner boundary sphere\n",
"positions_reduced, nb_boundary = mesher.point_cloud_add_spherical_inner_boundary(positions_reduced, nb_boundary, 0.04*R_max, healpy_order=healpy_order, origin=origin)\n",
"print(\"number of points in reduced grid: \", len(positions_reduced))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We add a spherical outer boundary at R_max"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"number of points in reduced grid: 290761\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 #the origin of the inner boundary sphere\n",
"positions_reduced, nb_boundary = mesher.point_cloud_add_spherical_outer_boundary(positions_reduced, nb_boundary, R_max, healpy_order=healpy_order, origin=origin)\n",
"print(\"number of points in reduced grid: \", len(positions_reduced))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"npoints = len(positions_reduced)\n",
"\n",
"XCO = 6.0e-4 # [.]\n",
"vturb = 1.5e+3 # [m/s]\n",
"\n",
"\n",
"def density(rr):\n",
" \"\"\"\n",
" Density function.\n",
" \"\"\"\n",
" (r, t, p) = spherical(rr[0], rr[1], rr[2])\n",
" if (r < R_0):\n",
" r = R_0\n",
" return rho_wind(r, t, p)\n",
"\n",
"\n",
"def abn_nH2(rr):\n",
" \"\"\"\n",
" H2 number density function.\n",
" \"\"\"\n",
" return density(rr) / (2.01588 * constants.u.si.value)\n",
"\n",
"\n",
"def abn_nCO(rr):\n",
" \"\"\"\n",
" CO number density function.\n",
" \"\"\"\n",
" return XCO * abn_nH2(rr)\n",
"\n",
"\n",
"def temp(rr):\n",
" \"\"\"\n",
" Temperature function.\n",
" \"\"\"\n",
" (r, t, p) = spherical(rr[0], rr[1], rr[2])\n",
" if (r < R_0):\n",
" r = R_0\n",
" return T_wind(r, t, p)\n",
"\n",
" \n",
"def turb(rr):\n",
" \"\"\"\n",
" !!! Peculiar Magritte thing...\n",
" Square turbulent speed as fraction of the speed of light.\n",
" \"\"\"\n",
" return (vturb/constants.c.si.value)**2\n",
"\n",
"\n",
"def velo(rr):\n",
" \"\"\"\n",
" Velocity function.\n",
" \"\"\"\n",
" (r, t, p) = spherical(rr[0], rr[1], rr[2])\n",
" if (r < R_0):\n",
" r = R_0\n",
" return v_r(r) / constants.c.si.value * np.array([rr[0], rr[1], rr[2]]) / r\n",
"\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": 15,
"metadata": {},
"outputs": [],
"source": [
"position = positions_reduced\n",
"velocity = np.array([velo (rr) for rr in positions_reduced])\n",
"nH2 = np.array([abn_nH2(rr) for rr in positions_reduced])\n",
"nCO = np.array([abn_nCO(rr) for rr in positions_reduced])\n",
"tmp = np.array([temp (rr) for rr in positions_reduced])\n",
"trb = np.array([turb (rr) for rr in positions_reduced])"
]
},
{
"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",
"\n",
"