3D Phantom
We create a Magritte model from a snapshot of 3D Phantom hydrodynamics simulation. The hydro model was kindly provided by Jolien Malfait. The Phantom binary files are used directly, and are loaded using the plons package.
Setup
Import the required functionalty.
[1]:
import magritte.setup as setup # Model setup
import magritte.core as magritte # Core functionality
import plons # Loading phantom data
import numpy as np # Data structures
import warnings # Hide warnings
warnings.filterwarnings('ignore') # especially for yt
import yt # 3D plotting
import os
from tqdm import tqdm # Progress bars
from astropy import constants # Unit conversions
from scipy.spatial import Delaunay # Finding neighbors
from yt.funcs import mylog # To avoid yt output
mylog.setLevel(40) # as error messages
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
Define a working directory (you will have to change this; it must be an absolute path).
[2]:
wdir = "/lhome/thomasc/Magritte-examples/Phantom_3D/"
Create the working directory.
[3]:
!mkdir -p $wdir
Define file names.
[4]:
dump_file = os.path.join(wdir, 'model_Phantom_3D' ) # Phantom full dump (snapshot)
setup_file = os.path.join(wdir, 'wind.setup' ) # Phantom setup file
input_file = os.path.join(wdir, 'wind.in' ) # Phantom input file
model_file = os.path.join(wdir, 'model_Phantom_3D.hdf5') # Resulting Magritte model
lamda_file = os.path.join(wdir, 'co.txt' ) # Line data file
We use a snapshot and data file that can be downloaded with the following links.
[5]:
dump_link = "https://github.com/Ensor-code/phantom-models/raw/main/Malfait+2021/v05e50/wind_v05e50"
setup_link = "https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2021/v05e50/wind.setup"
input_link = "https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2021/v05e50/wind.in"
lamda_link = "https://home.strw.leidenuniv.nl/~moldata/datafiles/co.dat"
Dowload the snapshot and the linedata (%%capture is just used to suppress the output).
[6]:
%%capture
print(setup_file)
!wget $dump_link --output-document $dump_file
!wget $setup_link --output-document $setup_file
!wget $input_link --output-document $input_file
!wget $lamda_link --output-document $lamda_file
Extract data
The script below extracts the required data from the snapshot phantom dump file.
[7]:
# Loading the data
setupData = plons.LoadSetup(wdir, "wind")
dumpData = plons.LoadFullDump(dump_file, setupData)
position = dumpData["position"]*1e-2 # position vectors [cm -> m]
velocity = dumpData["velocity"]*1e3 # velocity vectors [km/s -> m/s]
velocity = velocity/constants.c.si.value # velocity vectors [m/s -> 1/c]
rho = dumpData["rho"] # density [g/cm^3]
u = dumpData["u"] # internal energy density [erg/g]
tmp = dumpData["Tgas"] # temperature [K]
tmp[tmp<2.725] = 2.725 # Cut-off temperatures below 2.725 K
# Extract the number of points
npoints = len(rho)
# Convenience arrays
zeros = np.zeros(npoints)
ones = np.ones (npoints)
# Convert rho (total density) to abundances
nH2 = rho * 1.0e+6 * constants.N_A.si.value / 2.02
nCO = nH2 * 1.0e-4
# Define turbulence at 150 m/s
trb = (150.0/constants.c.si.value)**2 * ones
[8]:
# Extract Delaunay vertices (= Voronoi neighbors)
delaunay = Delaunay(position)
(indptr, indices) = delaunay.vertex_neighbor_vertices
neighbors = [indices[indptr[k]:indptr[k+1]] for k in range(npoints)]
nbs = [n for sublist in neighbors for n in sublist]
n_nbs = [len(sublist) for sublist in neighbors]
# Compute the indices of the boundary particles of the mesh, extracted from the Delaunay vertices
boundary = set([])
for i in tqdm(range(delaunay.neighbors.shape[0])):
for k in range(4):
if (delaunay.neighbors[i][k] == -1):
nk1,nk2,nk3 = (k+1)%4, (k+2)%4, (k+3)%4
boundary.add(delaunay.simplices[i][nk1])
boundary.add(delaunay.simplices[i][nk2])
boundary.add(delaunay.simplices[i][nk3])
boundary = list(boundary)
boundary = np.array(boundary)
# The above calculation turned out to be unsatisfactory.
# Since the outer boundary is assumed to be a sphere,
# we add all points which fall inside the boundary defined above.
b_nms = np.linalg.norm(position[boundary], axis=1)
p_nms = np.linalg.norm(position, axis=1)
boundary = np.array([i[0] for i in np.argwhere(p_nms >= np.min(b_nms))])
100%|██████████| 6960643/6960643 [00:38<00:00, 178604.41it/s]
Create model
Now all data is read, we can use it to construct a Magritte model.
Warning
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. 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.
[9]:
model = magritte.Model () # Create model object
model.parameters.set_model_name (model_file) # Magritte model file
model.parameters.set_dimension (3) # This is a 3D model
model.parameters.set_npoints (npoints) # Number of points
model.parameters.set_nrays (2) # Number of rays
model.parameters.set_nspecs (3) # Number of species
model.parameters.set_nlspecs (1) # Number of line species
model.parameters.set_nquads (31) # Number of quadrature points
model.geometry.points.position.set(position)
model.geometry.points.velocity.set(velocity)
model.geometry.points. neighbors.set( nbs)
model.geometry.points.n_neighbors.set(n_nbs)
model.chemistry.species.abundance = np.array((nCO, nH2, zeros)).T
model.chemistry.species.symbol = ['CO', 'H2', 'e-']
model.thermodynamics.temperature.gas .set(tmp)
model.thermodynamics.turbulence.vturb2.set(trb)
model.parameters.set_nboundary(boundary.shape[0])
model.geometry.boundary.boundary2point.set(boundary)
direction = np.array([[0,0,+1], [0,0,-1]]) # Comment out to use all directions
model.geometry.rays.direction.set(direction) # Comment out to use all directions
model.geometry.rays.weight .set(0.5 * np.ones(2)) # Comment out to use all directions
# model = setup.set_uniform_rays (model) # Uncomment to use all directions
model = setup.set_boundary_condition_CMB (model)
model = setup.set_linedata_from_LAMDA_file(model, lamda_file, {'considered transitions': [0]})
# model = setup.set_linedata_from_LAMDA_file(model, lamda_file) # Consider all transitions
model = setup.set_quadrature (model)
model.write()
Not considering all radiative transitions on the data file but only the specified ones!
Writing parameters...
Writing points...
Writing rays...
Writing boundary...
Writing chemistry...
Writing species...
Writing thermodynamics...
Writing temperature...
Writing turbulence...
Writing lines...
Writing lineProducingSpecies...
Writing linedata...
ncolpoar = 2
--- colpoar = 0
Writing collisionPartner...
(l, c) = 0, 0
--- colpoar = 1
Writing collisionPartner...
(l, c) = 0, 1
Writing quadrature...
Writing populations...
Writing radiation...
Writing frequencies...
Plot model
Load the data in a yt unstructured mesh.
[10]:
ds = yt.load_unstructured_mesh(
connectivity = delaunay.simplices.astype(np.int64),
coordinates = delaunay.points.astype(np.float64) * 100.0, # yt expects cm not m
node_data = {('connect1', 'n'): nCO[delaunay.simplices].astype(np.float64)}
)
Plot a slice through the mesh orthogonal to the z-axis.
[11]:
sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))
sl.set_cmap (('connect1', 'n'), 'magma')
sl.zoom (1.1)
[11]:
Show mesh on the plot.
[12]:
sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))
sl.set_cmap (('connect1', 'n'), 'magma')
sl.zoom (1.1)
sl.annotate_mesh_lines (plot_args={'color':'lightgrey', 'linewidths':[.25]})
[12]:
In the next example we demonstrate how to reduce this model as in De Ceuster et al. (2020).