{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3D Phantom - Reduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We reduce (or resample) the Magritte model of the 3D Phantom snapshot that was created in the [previous example](3_create_Phantom_3D.ipynb) using the methods explained in [Ceulemans et al. (2023) in prep]. This is an improvement over the previous implementation of [De Ceuster et al. (2020)](https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.5194D/abstract)."
]
},
{
"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 magrittetorch.tools.mesher as mesher\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 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, 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). We assume here that the scripts of the [previous example](3_create_Phantom_3D.ipynb) have already been executed and go back to that working directory."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"wdir = \"/lhome/thomasc/Magrittetorch-examples/Phantom_3D/\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define file names."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"model_file = os.path.join(wdir, 'model_Phantom_3D.hdf5') # Original Magritte model\n",
"lamda_file = os.path.join(wdir, 'co.txt' ) # Line data file\n",
"redux_file = os.path.join(wdir, 'model_Phantom_3D_red' ) # Reduced Magritte model (no extension!)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Load the original model."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Reading model from: /lhome/thomasc/Magrittetorch-examples/Phantom_3D/model_Phantom_3D.hdf5\n",
"Reading Magrittetorch model\n"
]
}
],
"source": [
"original_model = Model(model_file)\n",
"original_model.read()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Extract the data from the model by casting it into numpy arrays."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"position = original_model.geometry.points.position.get_astropy()\n",
"velocity = original_model.geometry.points.velocity.get_astropy()\n",
"boundary = original_model.geometry.boundary.boundary2point.get()#indices, so need no units\n",
"nCO = original_model.chemistry.species.abundance.get_astropy()[:,0]\n",
"nH2 = original_model.chemistry.species.abundance.get_astropy()[:,1]\n",
"tmp = original_model.thermodynamics.temperature.gas.get_astropy()\n",
"trb = original_model.thermodynamics.turbulence.vturb.get_astropy()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(Optional) Add a spherical outer boundary by setting all values to 0 outside of a certain radius (for example 90% of the model)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# radius = constants.au*180\n",
"\n",
"# indices = np.array([index for index, value in enumerate(position) if np.linalg.norm(value) >= radius])\n",
"\n",
"# velocity[indices] = 0.0 * units.km/units.s\n",
"# nCO[indices] = 1.0 / units.m**3 #CO cannot be set to 0, or the outer points will still be refined\n",
"# nH2[indices] = 0.0 / units.m**3\n",
"# tmp[indices] = 1.0 * units.K #Temperature cannot be set to 0\n",
"# trb[indices] = 0.0 * units.km/units.s"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the (new) built-in remeshing procedure for point clouds"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"new interior points: 79156\n",
"number boundary points: 386\n"
]
}
],
"source": [
"radius = constants.au*180\n",
"positions_reduced, nb_boundary = mesher.remesh_point_cloud(position, nCO, max_depth = 12, threshold = 4e-1, hullorder = 3)\n",
"positions_reduced, nb_boundary = mesher.point_cloud_add_spherical_outer_boundary(positions_reduced, nb_boundary, radius)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compare the number of points in the original and the reduced mesh."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"npoints = original_model.parameters.npoints.get()\n",
"npoints_reduced = len(positions_reduced)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"npoints original = 1034671\n",
"npoints reduced = 72950\n",
"reduction factor = 14.183289924605894\n"
]
}
],
"source": [
"print('npoints original =', npoints)\n",
"print('npoints reduced =', npoints_reduced)\n",
"print('reduction factor =', npoints/npoints_reduced)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Interpolate the data to the reduced mesh. Effectively we find out which points in the original mesh are closes to which point in the reduced mesh and we simpliy map the data between the correpsonding points."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# Find closest points\n",
"corresp_points = cKDTree(position).query(positions_reduced)[1]\n",
"\n",
"# Map data\n",
"position_reduced = positions_reduced\n",
"velocity_reduced = velocity[corresp_points]\n",
"nCO_reduced = nCO [corresp_points]\n",
"nH2_reduced = nH2 [corresp_points]\n",
"tmp_reduced = tmp [corresp_points]\n",
"trb_reduced = trb [corresp_points]\n",
"\n",
"# Extract Delaunay vertices (= Voronoi neighbors)\n",
"delaunay = Delaunay(position_reduced)\n",
"(indptr, indices) = delaunay.vertex_neighbor_vertices\n",
"neighbors = [indices[indptr[k]:indptr[k+1]] for k in range(npoints_reduced)]\n",
"nbs = [n for sublist in neighbors for n in sublist]\n",
"n_nbs = [len(sublist) for sublist in neighbors]\n",
"\n",
"# map to torch\n",
"# The reduction places boundary indices at start of list\n",
"boundary_torch = torch.arange((nb_boundary), dtype=torch.int64)\n",
"nbs_torch = torch.tensor(nbs, dtype=torch.int64)\n",
"n_nbs_torch = torch.tensor(n_nbs, dtype=torch.int64)\n",
"\n",
"# # Convenience arrays\n",
"# zeros = np.zeros(npoints_reduced)\n",
"# ones = np.ones (npoints_reduced)"
]
},
{
"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",
"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": 11,
"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_red.hdf5\n"
]
}
],
"source": [
"reduced_model = Model (f'{redux_file}.hdf5') # Create model object\n",
"reduced_model.geometry.geometryType.set(GeometryType.General3D) # This is a 3D model\n",
"reduced_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",
"reduced_model.geometry.points.position.set_astropy(position_reduced) # Set point positions\n",
"reduced_model.geometry.points.velocity.set_astropy(velocity_reduced) # Set point velocities\n",
"reduced_model.chemistry.species.abundance.set_astropy(np.stack([nCO_reduced, nH2_reduced, np.zeros(npoints_reduced)/units.m**3], axis=1))# Set species number densities\n",
"reduced_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",
"reduced_model.thermodynamics.temperature.gas.set_astropy(tmp_reduced) # Set gas temperature\n",
"reduced_model.thermodynamics.turbulence.vturb.set_astropy(trb_reduced) # Set turbulence velocity\n",
"\n",
"# Set the neighbors\n",
"reduced_model.geometry.points.neighbors.set(nbs_torch) # Set neighbors\n",
"reduced_model.geometry.points.n_neighbors.set(n_nbs_torch) # Set number of neighbors\n",
"reduced_model.geometry.boundary.boundary2point.set(boundary_torch) # Set which points are boundary points\n",
"\n",
"reduced_model = setup.set_boundary_condition_CMB (reduced_model) # Set CMB as boundary condition\n",
"reduced_model = setup.set_uniform_rays (reduced_model, 48) # 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",
"reduced_model = setup.set_linedata_from_LAMDA_file(reduced_model, lamda_file, {'considered transitions': [i for i in range(10)]})\n",
"# reduced_model = setup.set_linedata_from_LAMDA_file(reduced_model, lamda_file) # Consider all transitions\n",
"reduced_model = setup.set_quadrature (reduced_model, 7) # Set number of frequency quadrature points for NLTE radiative transfer\n",
"\n",
"reduced_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": 12,
"metadata": {},
"outputs": [],
"source": [
"ds = yt.load_unstructured_mesh(\n",
" connectivity = delaunay.simplices.astype(np.int64),\n",
" coordinates = position_reduced.to_value(units.cm), # yt expects cm not m\n",
" node_data = {('connect1', 'n'): nCO_reduced[delaunay.simplices].to_value(1/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": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))\n",
"sl.set_zlim( ('connect1', 'n'), 10**5, 10**(11)) #fix colorbar\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": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))\n",
"sl.set_zlim (('connect1', 'n'), 10**5, 10**(11)) #fix colorbar\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": "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
}