Synth. Obs.: Create and image 3D Phantom without saving

Although it’s recommended to initially create and save a model for the sake of reproducibility, Magrittetorch does not strictly demand this step before applying any postprocessing. In this example, we create a model starting from a 3D Phantom hydrodynamics model (as in this model creation example) and post-process it in the same way as this postprocessing example.

Setup

Import the required functionalty.

import magrittetorch.tools.radiativetransferutils as rtutils
import magrittetorch.tools.setup as setup
import magrittetorch.tools.plot as plot
import magrittetorch.algorithms.solvers as solvers
import torch
from magrittetorch.model.model import Model # Model class
from magrittetorch.model.geometry.geometry import GeometryType # GeometryType enum
from magrittetorch.model.geometry.boundary import BoundaryType # BoundaryType enum
# 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, units   # Unit conversions
from scipy.spatial  import Delaunay           # Finding neighbors
from yt.funcs       import mylog              # To avoid yt output
mylog.setLevel(40)                            # as error messages

Define a working directory (you will have to change this; it must be an absolute path).

wdir = "/lhome/thomasc/Magrittetorch-examples/Phantom_3D/"

Create the working directory.

!mkdir -p $wdir

Define file names.

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.

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).

%%capture
!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.

# Loading the data
setupData = plons.LoadSetup(wdir, "wind")
dumpData  = plons.LoadFullDump(dump_file, setupData)

position = dumpData["position"]*units.cm # position vectors [cm]
velocity = dumpData["velocity"]*units.km/units.s # velocity vectors [km/s]
rho      = dumpData["rho"]*units.g/units.cm**3 # density [g/cm^3]
u        = dumpData["u"] # internal energy density [erg/g]
tmp      = dumpData["Tgas"]*units.K # temperature [K]
tmp[tmp<2.725*units.K] = 2.725*units.K # Cut-off temperatures below 2.725 K

# Extract the number of points
npoints = len(rho)

# Convert rho (total density) to abundances
nH2 = rho * constants.N_A / (2.02 * units.g / units.mol)
nCO = nH2 * 1.0e-4

# Define turbulence at 150 m/s
trb = 150.0*units.m/units.s

# 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))])
#convert to torch tensor
boundary_torch = torch.from_numpy(boundary)
nbs_torch = torch.Tensor(nbs).type(torch.int64)
n_nbs_torch = torch.Tensor(n_nbs).type(torch.int64)

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 6960643/6960643 [00:33<00:00, 209299.10it/s]

Create model

Now all data is read, we can use it to construct a Magrittetorch 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.

model = Model("filename.hdf5") # Create model object
model.geometry.geometryType.set(GeometryType.General3D) # This is a 3D model
model.geometry.boundary.boundaryType.set(BoundaryType.Sphere3D) # With a spherical boundary

# In order to make unit conversions trivial, we use astropy quantities as input
model.geometry.points.position.set_astropy(position) # Set point positions
model.geometry.points.velocity.set_astropy(velocity) # Set point velocities
model.chemistry.species.abundance.set_astropy(np.stack([nCO, nH2, np.zeros(npoints)/units.m**3], axis=1))# Set species number densities
model.chemistry.species.symbol.set(np.array(['CO', 'H2', 'e-'], dtype='S')) #Set species symbols; should correspond to the LAMDA file format
#Note: the dtype='S' is necessary to correctly save and read the species symbols to/from the hdf5 file

model.thermodynamics.temperature.gas.set_astropy(tmp) # Set gas temperature
model.thermodynamics.turbulence.vturb.set_astropy(trb*np.ones(npoints)) # Set turbulence velocity

# Set the neighbors
model.geometry.points.neighbors.set(nbs_torch) # Set neighbors
model.geometry.points.n_neighbors.set(n_nbs_torch) # Set number of neighbors
model.geometry.boundary.boundary2point.set(boundary_torch) # Set which points are boundary points

# # For unitless quantities, we can also directly set the torch tensors
# nb_boundary = len(boundary)
# boundary_indices = torch.arange(nb_boundary, dtype = torch.int64)
# model.geometry.boundary.boundary2point.set(boundary_indices) # Set which points are boundary points
# model = setup.set_Delaunay_neighbor_lists (model) # Automatically computes and sets neighbors for each point, using a Delaunay triangulation

# Conveniently, the remeshing function puts the boundary points in front of the positions array
model = setup.set_boundary_condition_CMB  (model) # Set CMB as boundary condition
model = setup.set_uniform_rays            (model, 12) # Number of rays for NLTE raytracing; has be of the form 12*2**n

#As this example does not do NLTE, we might as well only consider the first 10 transitions of CO
model = setup.set_linedata_from_LAMDA_file(model, lamda_file, {'considered transitions': [i for i in range(10)]})
# model = setup.set_linedata_from_LAMDA_file(model, lamda_file)   # Consider all transitions
model = setup.set_quadrature              (model, 7) # Set number of frequency quadrature points for NLTE radiative transfer

# We are not saving the model in this example
# model.write()

Not considering all radiative transitions on the data file but only the specified ones!

Model the medium

Initilize the model

model.dataCollection.infer_data()

In this example we will work with the LTE level populations and do not demand statistical equilibrium.

# If you have a GPU, we can use it to speed up the computations
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")

# Getting the full model on GPU might be a bit too much for its memory, so you might want to use your CPU instead for the computations
# device = torch.device("cpu")

# Iterate level populations until statistical equilibrium
# solvers.compute_level_populations(model, device=device, 20)

Make synthetic observations

Now we can make synthetich observations of the model

# Get the line frequency corresponding to the CO 2-1 line
fcen = model.sources.lines.lineProducingSpecies[0].linedata.frequency.get(device)[1]
vpix = 1500   # velocity pixel size [m/s]
nfreqs = 31   # number of frequency bins
dd   = vpix * (nfreqs-1)/2 / constants.c.value # frequency bin size [Hz]
fmin = fcen - fcen*dd
fmax = fcen + fcen*dd

image_freqs = torch.linspace(fmin, fmax, nfreqs, device=device, dtype=torch.float32)

# Compute image in the specified direction at the given frequencies
solvers.image_model(model, torch.Tensor([0,0,1]).type(torch.float32).to(device), image_freqs, device)

Plot observations

Plot the resulting channel maps with the Magrittetorch matplotlib back end.

plot.image_mpl(
    model,
    image_nr =  -1,
    zoom     = 1.3,
    npix_x   = 256,
    npix_y   = 256,
    x_unit   = units.au,
    v_unit   = units.km / units.s,
    method = 'nearest'
)

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 31/31 [00:15<00:00,  1.98it/s]

<function magrittetorch.tools.plot.image_mpl.<locals>.<lambda>(v)>

(The plot is only interactive in a live notebook.)

Save the image cube in a fits file.

plot.save_fits(model)

WARNING: No regularly spaced frequency bins!
Written file to: /lhome/thomasc/github/Magritte-torch/docs/src/1_examples/2_synthetic_observations/images/image.fits

(Optional: To create your own plots) Overview of data stored in the Image object

xdir = model.images[-1].image_direction_x.get()#directions of the x-and y-vectors of the image
ydir = model.images[-1].image_direction_y.get()
zdir = model.images[-1].image_direction_z.get()#this is direction in which we observe the object
print("image directions: ", xdir, ydir, zdir)
nfreqs = model.images[-1].nfreqs.get() #number of frequency bins
freqs = model.images[-1].freqs.get_astropy() #frequency bins [Hz]
print("# of frequencies: ", nfreqs, " frequencies :", freqs)
ImX = model.images[-1].imX.get_astropy()#X position in image [m]
ImY = model.images[-1].imY.get_astropy()#Y position in image [m]
I = model.images[-1].I.get_astropy()#Intensity at the corresponding ImX, ImY position [W/(m^2*Hz*sr)], at a given frequency bin
# print("Intensities :", I, " ImX:", ImX, "ImY:", ImY) #prints a lot of output

image directions:  tensor([1., 0., 0.]) tensor([0.0000e+00, 1.0000e+00, 4.3711e-08]) tensor([ 0.0000e+00, -4.3711e-08,  1.0000e+00])
# of frequencies:  31  frequencies : [2.3051970e+11 2.3052085e+11 2.3052201e+11 2.3052316e+11 2.3052432e+11
 2.3052547e+11 2.3052662e+11 2.3052778e+11 2.3052893e+11 2.3053009e+11
 2.3053124e+11 2.3053238e+11 2.3053355e+11 2.3053469e+11 2.3053586e+11
 2.3053700e+11 2.3053815e+11 2.3053931e+11 2.3054046e+11 2.3054162e+11
 2.3054277e+11 2.3054392e+11 2.3054508e+11 2.3054623e+11 2.3054739e+11
 2.3054854e+11 2.3054968e+11 2.3055085e+11 2.3055199e+11 2.3055316e+11
 2.3055430e+11] Hz