3D Phantom - Reduction

We reduce (or resample) the Magritte model of the 3D Phantom snapshot that was created in the previous example 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).

Setup

Import the required functionalty.

import magrittetorch.tools.radiativetransferutils as rtutils
import magrittetorch.tools.setup as setup
import magrittetorch.tools.mesher as mesher
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 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, cKDTree  # 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). We assume here that the scripts of the previous example have already been executed and go back to that working directory.

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

Define file names.

model_file = os.path.join(wdir, 'model_Phantom_3D.hdf5')   # Original Magritte model
lamda_file = os.path.join(wdir, 'co.txt'               )   # Line data file
redux_file = os.path.join(wdir, 'model_Phantom_3D_red' )   # Reduced Magritte model (no extension!)

Load the original model.

original_model = Model(model_file)
original_model.read()

Reading model from:  /lhome/thomasc/Magrittetorch-examples/Phantom_3D/model_Phantom_3D.hdf5
Reading Magrittetorch model

Extract the data from the model by casting it into numpy arrays.

position = original_model.geometry.points.position.get_astropy()
velocity = original_model.geometry.points.velocity.get_astropy()
boundary = original_model.geometry.boundary.boundary2point.get()#indices, so need no units
nCO      = original_model.chemistry.species.abundance.get_astropy()[:,0]
nH2      = original_model.chemistry.species.abundance.get_astropy()[:,1]
tmp      = original_model.thermodynamics.temperature.gas.get_astropy()
trb      = original_model.thermodynamics.turbulence.vturb.get_astropy()

(Optional) Add a spherical outer boundary by setting all values to 0 outside of a certain radius (for example 90% of the model)

# radius = constants.au*180

# indices = np.array([index for index, value in enumerate(position) if np.linalg.norm(value) >= radius])

# velocity[indices] = 0.0 * units.km/units.s
# nCO[indices] = 1.0 / units.m**3                #CO cannot be set to 0, or the outer points will still be refined
# nH2[indices] = 0.0 / units.m**3
# tmp[indices] = 1.0 * units.K                #Temperature cannot be set to 0
# trb[indices] = 0.0 * units.km/units.s

Use the (new) built-in remeshing procedure for point clouds

radius = constants.au*180
positions_reduced, nb_boundary = mesher.remesh_point_cloud(position, nCO, max_depth = 12, threshold = 4e-1, hullorder = 3)
positions_reduced, nb_boundary = mesher.point_cloud_add_spherical_outer_boundary(positions_reduced, nb_boundary, radius)

new interior points:  79156
number boundary points:  386

Compare the number of points in the original and the reduced mesh.

npoints          = original_model.parameters.npoints.get()
npoints_reduced  = len(positions_reduced)

print('npoints original =', npoints)
print('npoints reduced  =', npoints_reduced)
print('reduction factor =', npoints/npoints_reduced)

npoints original = 1034671
npoints reduced  = 72950
reduction factor = 14.183289924605894

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.

# Find closest points
corresp_points = cKDTree(position).query(positions_reduced)[1]

# Map data
position_reduced = positions_reduced
velocity_reduced = velocity[corresp_points]
nCO_reduced      = nCO     [corresp_points]
nH2_reduced      = nH2     [corresp_points]
tmp_reduced      = tmp     [corresp_points]
trb_reduced      = trb     [corresp_points]

# Extract Delaunay vertices (= Voronoi neighbors)
delaunay = Delaunay(position_reduced)
(indptr, indices) = delaunay.vertex_neighbor_vertices
neighbors = [indices[indptr[k]:indptr[k+1]] for k in range(npoints_reduced)]
nbs       = [n for sublist in neighbors for n in sublist]
n_nbs     = [len(sublist) for sublist in neighbors]

# map to torch
# The reduction places boundary indices at start of list
boundary_torch = torch.arange((nb_boundary), dtype=torch.int64)
nbs_torch      = torch.tensor(nbs,      dtype=torch.int64)
n_nbs_torch    = torch.tensor(n_nbs,    dtype=torch.int64)

# # Convenience arrays
# zeros = np.zeros(npoints_reduced)
# ones  = np.ones (npoints_reduced)

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.

reduced_model = Model (f'{redux_file}.hdf5')                     # Create model object
reduced_model.geometry.geometryType.set(GeometryType.General3D) # This is a 3D model
reduced_model.geometry.boundary.boundaryType.set(BoundaryType.Sphere3D) # With a spherical boundary

# In order to make unit conversions trivial, we use astropy quantities as input
reduced_model.geometry.points.position.set_astropy(position_reduced) # Set point positions
reduced_model.geometry.points.velocity.set_astropy(velocity_reduced) # Set point velocities
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
reduced_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

reduced_model.thermodynamics.temperature.gas.set_astropy(tmp_reduced) # Set gas temperature
reduced_model.thermodynamics.turbulence.vturb.set_astropy(trb_reduced) # Set turbulence velocity

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

reduced_model = setup.set_boundary_condition_CMB  (reduced_model) # Set CMB as boundary condition
reduced_model = setup.set_uniform_rays            (reduced_model, 48) # 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
reduced_model = setup.set_linedata_from_LAMDA_file(reduced_model, lamda_file, {'considered transitions': [i for i in range(10)]})
# reduced_model = setup.set_linedata_from_LAMDA_file(reduced_model, lamda_file)   # Consider all transitions
reduced_model = setup.set_quadrature              (reduced_model, 7) # Set number of frequency quadrature points for NLTE radiative transfer

reduced_model.write()

Not considering all radiative transitions on the data file but only the specified ones!
Writing model to:  /lhome/thomasc/Magrittetorch-examples/Phantom_3D/model_Phantom_3D_red.hdf5

Plot model

Load the data in a yt unstructured mesh.

ds = yt.load_unstructured_mesh(
         connectivity = delaunay.simplices.astype(np.int64),
         coordinates  = position_reduced.to_value(units.cm), # yt expects cm not m
         node_data    = {('connect1', 'n'): nCO_reduced[delaunay.simplices].to_value(1/units.m**3)}
)

Plot a slice through the mesh orthogonal to the z-axis.

sl = yt.SlicePlot (ds, 'z', ('connect1', 'n'))
sl.set_zlim( ('connect1', 'n'), 10**5, 10**(11))   #fix colorbar
sl.set_cmap       (('connect1', 'n'), 'magma')
sl.zoom           (1.1)

Show mesh on the plot.

sl = yt.SlicePlot      (ds, 'z', ('connect1', 'n'))
sl.set_zlim            (('connect1', 'n'), 10**5, 10**(11))   #fix colorbar
sl.set_cmap            (('connect1', 'n'), 'magma')
sl.zoom                (1.1)
sl.annotate_mesh_lines (plot_args={'color':'lightgrey', 'linewidths':[.25]})