{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3D AMRVAC"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We create a Magritte model from a snapshot of 3D AMRVAC hydrodynamics simulation.\n",
"The hydro model was kindly provided by Jan Bolte.\n",
"Currently, the AMRVAC binary files can not yet be used directly to extract the snapshot data.\n",
"Hence we use the corresponding `.vtu` files."
]
},
{
"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\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; it must be an **absolute path**)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"wdir = \"/lhome/thomasc/Magritte-examples/AMRVAC_3D/\""
]
},
{
"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": [
"input_file = os.path.join(wdir, 'model_AMRVAC_3D.vtu' ) # AMRVAC snapshot\n",
"model_file = os.path.join(wdir, 'model_AMRVAC_3D.hdf5') # Resulting Magritte model\n",
"lamda_file = os.path.join(wdir, 'co.txt' ) # Line data file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use a snapshot and data file that can be downloaded with the following links."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"input_link = \"https://owncloud.ster.kuleuven.be/index.php/s/mqtDyDSMPm2TjmG/download\"\n",
"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 $input_link --output-document $input_file\n",
"!wget $lamda_link --output-document $lamda_file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extract data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The script below extracts the required data from the snapshot `.vtu` file."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████████| 89/89 [00:00<00:00, 307237.08it/s]\n",
"100%|████████████████████████████████| 297984/297984 [00:04<00:00, 73535.32it/s]\n",
"100%|██████████████████████████████| 2141089/2141089 [00:21<00:00, 99629.04it/s]\n"
]
}
],
"source": [
"# Create a vtk reader to read the AMRVAC file and extract its contents.\n",
"reader = vtk.vtkXMLUnstructuredGridReader()\n",
"reader.SetFileName(input_file)\n",
"reader.Update()\n",
"\n",
"# Extract the grid output\n",
"grid = reader.GetOutput()\n",
"\n",
"# Extract the number of cells\n",
"ncells = grid.GetNumberOfCells() \n",
"\n",
"# Extract cell data\n",
"cellData = grid.GetCellData()\n",
"for i in tqdm(range(cellData.GetNumberOfArrays())):\n",
" array = cellData.GetArray(i)\n",
" if (array.GetName() == 'rho'):\n",
" rho = vtk_to_numpy(array)\n",
" if (array.GetName() == 'temperature'):\n",
" tmp = vtk_to_numpy(array)\n",
" if (array.GetName() == 'v1'):\n",
" v_x = vtk_to_numpy(array)\n",
" if (array.GetName() == 'v2'):\n",
" v_y = vtk_to_numpy(array)\n",
" if (array.GetName() == 'v3'):\n",
" v_z = vtk_to_numpy(array)\n",
"\n",
"# Convert rho (total density) to abundances\n",
"nH2 = rho * 1.0e+6 * constants.N_A.si.value / 2.02\n",
"nCO = nH2 * 1.0e-4\n",
"\n",
"# Convenience arrays\n",
"zeros = np.zeros(ncells)\n",
"ones = np.ones (ncells)\n",
"\n",
"# Convert to fractions of the speed of light\n",
"velocity = np.array((v_x, v_y, v_z)).transpose() / constants.c.cgs.value\n",
"\n",
"# Define turbulence at 150 m/s\n",
"trb = (150.0/constants.c.si.value)**2 * ones\n",
"\n",
"# Extract cell centres as Magritte points\n",
"centres = []\n",
"for c in tqdm(range(ncells)):\n",
" cell = grid.GetCell(c)\n",
" centre = np.zeros(3)\n",
" for i in range(8):\n",
" centre = centre + np.array(cell.GetPoints().GetPoint(i))\n",
" centre = 0.125 * centre\n",
" centres.append(centre)\n",
"centres = np.array(centres)\n",
"centres *= 1.0e-2 # convert [cm] to [m]\n",
"position = centres\n",
"\n",
"# We assume that the geometry to be rectangular. We still need to define the boundary however.\n",
"# For adaptive mesh refinement models, we put a boundary box around it (to ensure that we have a convex boundary around the model; this is mainly important for imaging the model).\n",
"\n",
"delta_xyz = (np.max(centres, axis=0) - np.min(centres, axis=0))/2.0#half the width in each direction\n",
"center_xyz = (np.max(centres, axis=0) + np.min(centres, axis=0))/2.0\n",
"xyz_max = center_xyz + 1.001*delta_xyz\n",
"xyz_min = center_xyz - 1.001*delta_xyz\n",
"\n",
"# Add a boundary box around the model, putting the boundary points on a regular grid\n",
"box = mesher.create_cubic_uniform_hull(xyz_min, xyz_max, order=5)#the boundary box contains 2**order+1 points in each direction.\n",
"nbox = len(box)\n",
"\n",
"# Find closest points, in order to extrapolate the values towards the boundary box.\n",
"corresp_points = cKDTree(position).query(box)[1]\n",
"\n",
"# Map data of model to box\n",
"position_box = box\n",
"velocity_box = velocity[corresp_points]\n",
"nCO_box = nCO [corresp_points]\n",
"nH2_box = nH2 [corresp_points]\n",
"tmp_box = tmp [corresp_points]\n",
"trb_box = trb [corresp_points]\n",
"\n",
"#and extend data arrays\n",
"nCO = np.concatenate((nCO_box, nCO))\n",
"nH2 = np.concatenate((nH2_box, nH2))\n",
"zeros = np.concatenate((np.zeros(nbox), zeros))\n",
"ones = np.concatenate((np.ones(nbox), ones))\n",
"velocity = np.concatenate((velocity_box, velocity), axis=0)\n",
"tmp = np.concatenate((tmp_box, tmp))\n",
"trb = np.concatenate((trb_box, trb))\n",
"position = np.concatenate((box, position), axis=0)\n",
"\n",
"\n",
"# Extract Delaunay vertices (= Voronoi neighbors)\n",
"delaunay = Delaunay(position)\n",
"(indptr, indices) = delaunay.vertex_neighbor_vertices\n",
"neighbors = [indices[indptr[k]:indptr[k+1]] for k in range(len(position))]\n",
"nbs = [n for sublist in neighbors for n in sublist]\n",
"n_nbs = [len(sublist) for sublist in neighbors]\n",
"\n",
"# Compute the indices of the boundary particles of the mesh, extracted from the Delaunay vertices\n",
"# Note: technically, we already know the boundary (the constructed box), but as we compute a Delaunay triangulation either way to find the neighbors, we might as well use it to compute the boundary points.\n",
"boundary = set([])\n",
"for i in tqdm(range(delaunay.neighbors.shape[0])):\n",
" for k in range(4):\n",
" if (delaunay.neighbors[i][k] == -1):\n",
" nk1,nk2,nk3 = (k+1)%4, (k+2)%4, (k+3)%4 \n",
" boundary.add(delaunay.simplices[i][nk1])\n",
" boundary.add(delaunay.simplices[i][nk2])\n",
" boundary.add(delaunay.simplices[i][nk3])\n",
"\n",
"boundary = list(boundary)\n",
"boundary = np.array(boundary)"
]
},
{
"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",
"
"
]
},
{
"cell_type": "code",
"execution_count": 8,
"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 (model_file) # Magritte model file\n",
"model.parameters.set_dimension (3) # This is a 3D model\n",
"model.parameters.set_npoints (len(position)) # Number of points\n",
"model.parameters.set_nrays (12) # 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)\n",
"model.geometry.points.velocity.set(velocity)\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, nH2, zeros)).T\n",
"model.chemistry.species.symbol = ['CO', 'H2', 'e-']\n",
"\n",
"model.thermodynamics.temperature.gas .set(tmp)\n",
"model.thermodynamics.turbulence.vturb2.set(trb)\n",
"\n",
"model.parameters.set_nboundary(boundary.shape[0])\n",
"# model.geometry.boundary.boundary2point.set(np.arange(boundary.shape[0]))\n",
"model.geometry.boundary.boundary2point.set(boundary)\n",
"\n",
"setup.set_boundary_condition_CMB (model)\n",
"setup.set_uniform_rays (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": 9,
"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."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 10,
"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": "markdown",
"metadata": {},
"source": [
"Show mesh on the plot."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" "
],
"text/plain": [
""
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))\n",
"sl.set_cmap (('connect1', 'n'), 'magma')\n",
"sl.zoom (1.0)\n",
"sl.annotate_mesh_lines (plot_args={'color':'lightgrey', 'linewidths':[.25]})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the [next example](2_reduce_AMRVAC_3D.ipynb) we demonstrate how to reduce this model as in [De Ceuster et al. (2020)](https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.5194D/abstract)."
]
}
],
"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
}