{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3D AMRVAC - Reduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We reduce (or resample) the Magritte model of the 3D AMRVAC snapshot that was created in the [previous example](1_create_AMRVAC_3D.ipynb) using the methods explained in [Ceulemans et al. (2024) 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](1_create_AMRVAC_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/AMRVAC_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_AMRVAC_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_AMRVAC_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/AMRVAC_3D/model_AMRVAC_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": [
"Use the (new) built-in remeshing procedure for point clouds"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"new interior points: 30635\n",
"number boundary points: 1538\n"
]
}
],
"source": [
"positions_reduced, nb_boundary = mesher.remesh_point_cloud(position, nCO, max_depth=9, threshold= 4e-1, hullorder=4)\n",
"boundary_reduced = positions_reduced[:nb_boundary,:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compare the number of points in the original and the reduced mesh."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"npoints = original_model.parameters.npoints.get()\n",
"npoints_reduced = len(positions_reduced)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"npoints original = 304130\n",
"npoints reduced = 32173\n",
"reduction factor = 9.452957448792466\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": 9,
"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)\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)"
]
},
{
"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": 10,
"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/AMRVAC_3D/model_AMRVAC_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.AxisAlignedCube) # With a cubic 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, 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",
"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": 11,
"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": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))\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": 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_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.11.4"
},
"widgets": {
"application/vnd.jupyter.widget-state+json": {
"state": {},
"version_major": 2,
"version_minor": 0
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}