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.
[17]:
import magritte.setup as setup # Model setup
import magritte.core as magritte # Core functionality
import magritte.mesher as mesher # Mesher
import numpy as np # Data structures
import vtk # Reading the model
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 vtk.util.numpy_support import vtk_to_numpy # Converting data
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.
[18]:
wdir = "/lhome/thomasc/Magritte-examples/Phantom_3D/"
Define file names.
[19]:
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.
[20]:
model = magritte.Model(model_file)
-------------------------------------------
Reading Model...
-------------------------------------------
model file = /home/owenv/Magritte/docs/src/1_examples/1_post-processing/model_Phantom_3D.hdf5
-------------------------------------------
Reading parameters...
Reading points...
Reading rays...
Reading boundary...
Reading chemistry...
Reading species...
Reading thermodynamics...
Reading temperature...
Reading turbulence...
Reading lines...
Reading lineProducingSpecies...
Reading linedata...
read num 0
read sym CO
nlev = 41
nrad = 1
Reading collisionPartner...
Reading collisionPartner...
Reading quadrature...
Reading radiation...
Reading frequencies...
Not using scattering!
-------------------------------------------
Model read, parameters:
-------------------------------------------
npoints = 1034671
nrays = 2
nboundary = 12634
nfreqs = 31
nspecs = 3
nlspecs = 1
nlines = 1
nquads = 31
-------------------------------------------
Extract the data from the model by casting it into numpy arrays.
[21]:
position = np.array(model.geometry.points.position)
velocity = np.array(model.geometry.points.velocity)
boundary = np.array(model.geometry.boundary.boundary2point)
nCO = np.array(model.chemistry.species.abundance)[:,0]
nH2 = np.array(model.chemistry.species.abundance)[:,1]
tmp = np.array(model.thermodynamics.temperature.gas)
trb = np.array(model.thermodynamics.turbulence.vturb2)
(Optional) Add a spherical outer boundary by setting all values to 0 outside of a certain radius (for example 90% of the model)
[22]:
radius = constants.au.si.value*180
indices = np.array([index for index, value in enumerate(position) if np.linalg.norm(value) >= radius])
velocity[indices] = 0
nCO[indices] = 1 #CO cannot be set to 0, or the outer points will still be refined
nH2[indices] = 0
tmp[indices] = 1 #Temperature cannot be set to 0
trb[indices] = 0
Use the (new) built-in remeshing procedure for point clouds
[23]:
positions_reduced, nb_boundary = mesher.remesh_point_cloud(position, nCO, max_depth = 12, threshold = 4e-1, hullorder = 3)
new interior points: 99773
number boundary points: 386
Compare the number of points in the original and the reduced mesh.
[24]:
npoints = model.parameters.npoints()
npoints_reduced = len(positions_reduced)
[25]:
print('npoints original =', npoints)
print('npoints reduced =', npoints_reduced)
print('reduction factor =', npoints/npoints_reduced)
npoints original = 1034671
npoints reduced = 100159
reduction factor = 10.330284847093122
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.
[26]:
# 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]
# Convenience arrays
zeros = np.zeros(npoints_reduced)
ones = np.ones (npoints_reduced)
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.
[27]:
model = magritte.Model () # Create model object
model.parameters.set_model_name (f'{redux_file}.hdf5') # Magritte model file
model.parameters.set_dimension (3) # This is a 3D model
model.parameters.set_npoints (npoints_reduced) # Number of points
model.parameters.set_nrays (2) # Number of rays
model.parameters.set_nspecs (3) # Number of species (min. 5)
model.parameters.set_nlspecs (1) # Number of line species
model.parameters.set_nquads (31) # Number of quadrature points
model.geometry.points.position.set(position_reduced)
model.geometry.points.velocity.set(velocity_reduced)
model.geometry.points. neighbors.set( nbs)
model.geometry.points.n_neighbors.set(n_nbs)
model.chemistry.species.abundance = np.array((nCO_reduced, nH2_reduced, zeros)).T
model.chemistry.species.symbol = ['CO', 'H2', 'e-']
model.thermodynamics.temperature.gas .set(tmp_reduced)
model.thermodynamics.turbulence.vturb2.set(trb_reduced)
model.parameters.set_nboundary(nb_boundary)
model.geometry.boundary.boundary2point.set(range(nb_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.
[28]:
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_reduced[delaunay.simplices].astype(np.float64)}
)
Plot a slice through the mesh orthogonal to the z-axis.
[29]:
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)
[29]:
Show mesh on the plot.
[30]:
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]})
[30]:
