{
"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": 17,
"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 vtk # Reading the model\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 # Unit conversions\n",
"from vtk.util.numpy_support import vtk_to_numpy # Converting data\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": 18,
"metadata": {},
"outputs": [],
"source": [
"wdir = \"/lhome/thomasc/Magritte-examples/Phantom_3D/\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define file names."
]
},
{
"cell_type": "code",
"execution_count": 19,
"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": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
"-------------------------------------------\n",
" Reading Model... \n",
"-------------------------------------------\n",
" model file = /home/owenv/Magritte/docs/src/1_examples/1_post-processing/model_Phantom_3D.hdf5\n",
"-------------------------------------------\n",
"Reading parameters...\n",
"Reading points...\n",
"Reading rays...\n",
"Reading boundary...\n",
"Reading chemistry...\n",
"Reading species...\n",
"Reading thermodynamics...\n",
"Reading temperature...\n",
"Reading turbulence...\n",
"Reading lines...\n",
"Reading lineProducingSpecies...\n",
"Reading linedata...\n",
"read num 0\n",
"read sym CO\n",
"nlev = 41\n",
"nrad = 1\n",
"Reading collisionPartner...\n",
"Reading collisionPartner...\n",
"Reading quadrature...\n",
"Reading radiation...\n",
"Reading frequencies...\n",
"Not using scattering!\n",
" \n",
"-------------------------------------------\n",
" Model read, parameters: \n",
"-------------------------------------------\n",
" npoints = 1034671\n",
" nrays = 2\n",
" nboundary = 12634\n",
" nfreqs = 31\n",
" nspecs = 3\n",
" nlspecs = 1\n",
" nlines = 1\n",
" nquads = 31\n",
"-------------------------------------------\n",
" \n"
]
}
],
"source": [
"model = magritte.Model(model_file)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Extract the data from the model by casting it into numpy arrays."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"position = np.array(model.geometry.points.position)\n",
"velocity = np.array(model.geometry.points.velocity)\n",
"boundary = np.array(model.geometry.boundary.boundary2point)\n",
"nCO = np.array(model.chemistry.species.abundance)[:,0]\n",
"nH2 = np.array(model.chemistry.species.abundance)[:,1]\n",
"tmp = np.array(model.thermodynamics.temperature.gas)\n",
"trb = np.array(model.thermodynamics.turbulence.vturb2)"
]
},
{
"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": 22,
"metadata": {},
"outputs": [],
"source": [
"radius = constants.au.si.value*180\n",
"\n",
"indices = np.array([index for index, value in enumerate(position) if np.linalg.norm(value) >= radius])\n",
"\n",
"velocity[indices] = 0\n",
"nCO[indices] = 1 #CO cannot be set to 0, or the outer points will still be refined\n",
"nH2[indices] = 0\n",
"tmp[indices] = 1 #Temperature cannot be set to 0\n",
"trb[indices] = 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the (new) built-in remeshing procedure for point clouds"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"new interior points: 99773\n",
"number boundary points: 386\n"
]
}
],
"source": [
"positions_reduced, nb_boundary = mesher.remesh_point_cloud(position, nCO, max_depth = 12, threshold = 4e-1, hullorder = 3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compare the number of points in the original and the reduced mesh."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"npoints = model.parameters.npoints()\n",
"npoints_reduced = len(positions_reduced)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"npoints original = 1034671\n",
"npoints reduced = 100159\n",
"reduction factor = 10.330284847093122\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": 26,
"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",
"# 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 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",
"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": 27,
"metadata": {},
"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 (f'{redux_file}.hdf5') # Magritte model file\n",
"model.parameters.set_dimension (3) # This is a 3D model\n",
"model.parameters.set_npoints (npoints_reduced) # 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 (31) # Number of quadrature points\n",
"\n",
"model.geometry.points.position.set(position_reduced)\n",
"model.geometry.points.velocity.set(velocity_reduced)\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_reduced, nH2_reduced, zeros)).T\n",
"model.chemistry.species.symbol = ['CO', 'H2', 'e-']\n",
"\n",
"model.thermodynamics.temperature.gas .set(tmp_reduced)\n",
"model.thermodynamics.turbulence.vturb2.set(trb_reduced)\n",
"\n",
"model.parameters.set_nboundary(nb_boundary)\n",
"model.geometry.boundary.boundary2point.set(range(nb_boundary))\n",
"\n",
"direction = np.array([[0,0,+1], [0,0,-1]]) # 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": 28,
"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_reduced[delaunay.simplices].astype(np.float64)}\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot a slice through the mesh orthogonal to the z-axis."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 29,
"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": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 30,
"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]})"
]
}
],
"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.9.16"
},
"widgets": {
"application/vnd.jupyter.widget-state+json": {
"state": {},
"version_major": 2,
"version_minor": 0
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}