Skip to content

Commit

Permalink
radmc-3d: Adding mass conversation check, changing theta grid, adding…
Browse files Browse the repository at this point in the history 
… metadata file.
  • Loading branch information
@stammler
stammler committed Jun 15, 2023
1 parent 486e070 commit 13d9bbc
Showing 1 changed file with 107 additions and 27 deletions.
134 changes: 107 additions & 27 deletions dustpylib/radtrans/radmc3d/radmc3d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@

import dsharp_opac as do
from dustpy import Simulation
from dustpy import constants as c
from dustpylib.grid.refinement import refine_radial_local
import numpy as np
import os
Expand Down Expand Up @@ -92,8 +93,8 @@ def __init__(self, sim, ignore_last=True):
#: Dust scale heights array in cm from ``DustPy``
self.H_dust_ = None

#: Dust midplane density profile in g/cm³ from ``DustPy``
self.rho_dust_ = None
#: Dust surface density profile in g/cm² from ``DustPy``
self.Sigma_dust_ = None

#: Directory to store the ``RADMC-3D input files``
self.datadir = "radmc3d"
Expand Down Expand Up @@ -122,22 +123,23 @@ def __init__(self, sim, ignore_last=True):
self.T_gas_ = self.T_gas_[:-1]
self.a_dust_ = self.a_dust_[:-1, :]
self.H_dust_ = self.H_dust_[:-1, :]
self.rho_dust_ = self.rho_dust_[:-1, :]
self.Sigma_dust_ = self.Sigma_dust_[:-1, :]

#: Radial grid cell interfaces for ``RADMC-3D`` model
self.ri_grid = refine_radial_local(self.ri_grid_, 0., num=3)

lam1 = np.geomspace(0.1e-4, 7.e-4, 20, endpoint=False)
lam2 = np.geomspace(7.e-4, 25.e-4, 100, endpoint=False)
lam3 = np.geomspace(25.e-4, 1., 30, endpoint=True)
self.lam_grid = np.concatenate([lam1, lam2, lam3])

Ntheta = 64
theta_threshold = 0.25*np.pi
theta1 = np.linspace(0., theta_threshold, 5, endpoint=False)
theta2 = np.linspace(theta_threshold,
0.5*np.pi, Ntheta+1-theta1.shape[0])
theta_up = np.concatenate([theta1, theta2])
linthresh = 1.e-3
Ntheta = 128
Ntheta_lin = 32
thetai_lin = np.pi/2 \
- np.linspace(0, linthresh, Ntheta_lin, endpoint=False)
thetai_log = np.pi/2 \
- np.geomspace(linthresh, np.pi/2, Ntheta+1-Ntheta_lin)
theta_up = np.hstack((thetai_lin, thetai_log))[::-1]
theta_down = (np.pi - theta_up[:-1])[::-1]
self.thetai_grid = np.concatenate([theta_up, theta_down])

Expand Down Expand Up @@ -268,7 +270,7 @@ def _init_from_dustpy(self, sim):

self.a_dust_ = sim.dust.a
self.H_dust_ = sim.dust.H
self.rho_dust_ = sim.dust.rho
self.Sigma_dust_ = sim.dust.Sigma

def _init_from_namespace(self, sim):
"""
Expand All @@ -292,7 +294,7 @@ def _init_from_namespace(self, sim):

self.a_dust_ = sim.dust.a
self.H_dust_ = sim.dust.H
self.rho_dust_ = sim.dust.rho
self.Sigma_dust_ = sim.dust.Sigma

def write_files(self, datadir=None, opacity=None, write_opacities=True):
"""
Expand All @@ -319,6 +321,7 @@ def write_files(self, datadir=None, opacity=None, write_opacities=True):
self._write_dust_temperature_dat(datadir=datadir)
if write_opacities:
self.write_opacity_files(datadir=datadir, opacity=opacity)
self._write_metadata(datadir=datadir)

def write_opacity_files(self, datadir=None, opacity=None):
"""
Expand Down Expand Up @@ -478,37 +481,63 @@ def _write_dust_density_inp(self, datadir=None):
r_grid = R_grid * np.sin(theta_grid)
z_grid = R_grid * np.cos(theta_grid)

rho_grid = np.empty(
(self.rc_grid.shape[0],
self.thetac_grid.shape[0],
self.phic_grid.shape[0],
# Binning surface density onto new size bins
Sigma = np.empty(
(self.rc_grid_.shape[0],
self.ac_grid.shape[0])
)
for i in range(Sigma.shape[1]):
mask = (
(self.a_dust_ >= self.ai_grid[i]) &
(self.a_dust_ < self.ai_grid[i+1])
)
Sigma[:, i] = np.where(mask, self.Sigma_dust_, 0.).sum(-1)

# Interpolate dust scale heights on size grid
# Interpolate dust scale heights on RADMC-3D size grid
x = np.repeat(self.rc_grid_, self.a_dust_.shape[1])
y = self.a_dust_.flatten()
z = self.H_dust_.flatten()
xi = self.rc_grid_
yi = self.ac_grid
H_grid_ = griddata(
H = griddata(
(x, y), z, (xi[:, None], yi[None, :]),
method="linear", rescale=True)

rho = Sigma / (np.sqrt(2.*np.pi)*H)
rho_grid = np.empty(
(self.rc_grid.shape[0],
self.thetac_grid.shape[0],
self.phic_grid.shape[0],
self.ac_grid.shape[0])
)

for i in range(self.ac_grid.shape[0]):
rho = np.where(
(self.a_dust_ >= self.ai_grid[i])
& (self.a_dust_ < self.ai_grid[i+1]),
self.rho_dust_,
0.
).sum(-1)
f_H = interp1d(self.rc_grid_, H_grid_[:, i],
f_H = interp1d(self.rc_grid_, H[..., i],
fill_value="extrapolate", kind="linear")
f_rho = interp1d(self.rc_grid_, rho,
f_rho = interp1d(self.rc_grid_, rho[..., i],
fill_value="extrapolate", kind="linear")
rho_grid[..., i] = np.maximum(
1.e-100, f_rho(r_grid) * np.exp(-0.5*(z_grid/f_H(r_grid))**2))

# Check for mass conservation
# Mass from DustPy import
M_dp = (np.pi * (self.ri_grid_[1:]**2-self.ri_grid_[:-1]**2)
* self.Sigma_dust_[..., :].sum(-1)).sum()
# Computing volume elements
Ri_grid, thetai_grid, phii_grid = np.meshgrid(
self.ri_grid, self.thetai_grid, self.phii_grid,
indexing="ij"
)
dR = np.diff(Ri_grid, axis=0)[:, :-1, :-1]
dtheta = np.diff(thetai_grid, axis=1)[:-1, :, :-1]
dphi = np.diff(phii_grid, axis=2)[:-1, :-1, :]
dV = R_grid**2 * np.sin(theta_grid) * dR * dtheta * dphi
# Getting Mass
M = (rho_grid.sum(-1)*dV).sum()
if np.isclose(self.thetai_grid[-1], 0.5*np.pi):
M *= 2.
dM = np.abs(M-M_dp)/M_dp

with open(path, "w") as f:
f.write("1\n")
f.write("{:d}\n".format(rho_grid[..., 0].flatten().shape[0]))
Expand All @@ -517,6 +546,15 @@ def _write_dust_density_inp(self, datadir=None):
f.write("{:.6e}\n".format(val))
print("done.")

if dM > 1.e-3:
print()
print(" Warning: The total dust mass in the DustPy setup")
print(" and the RADMC-3D model differ.")
print(" M_dustpy = {:.3e} M_sun".format(M_dp/c.M_sun))
print(" M_radmc3d = {:.3e} M_sun".format(M/c.M_sun))
print(" If this is not intended, please refine the grids.")
print()

def _write_dust_temperature_dat(self, datadir=None):
"""
Function writes the 'dust_temperature.dat' input file.
Expand Down Expand Up @@ -599,10 +637,13 @@ def _write_dustkapscatmat_inp(self, datadir=None, opacity=None):
datadir : str, optional, default: None
Data directory in which the files are written. None defaults to
the datadir attribute of the parent class.
opacity : str, optional, default: None
Opacity model to be used. Either 'birnstiel2018' or 'ricci2010'.
None defaults to 'birnstiel2018'.
"""

datadir = self.datadir if datadir is None else datadir
opacity = opacity or self.opacity or 'birnstiel2018'
opacity = opacity or self.opacity or "birnstiel2018"
Path(datadir).mkdir(parents=True, exist_ok=True)

Nangle = 181
Expand Down Expand Up @@ -669,6 +710,35 @@ def _write_dustkapscatmat_inp(self, datadir=None, opacity=None):
)
print("done.")

def _write_metadata(self, datadir=None):
"""
Function writes the 'metadata*.npz' file.
It is not required for ``RADMC-3D`` model, but it contains
additional data such as particle size bins used in the setup.
Parameters
----------
datadir : str, optional, default: None
Data directory in which the files are written. None defaults to
the datadir attribute of the parent class.
"""

filename = "metadata.npz"
datadir = self.datadir if datadir is None else datadir
Path(datadir).mkdir(parents=True, exist_ok=True)

path = os.path.join(datadir, filename)

print("Writing {}.....".format(path), end="")

np.savez(
path,
a=self.ac_grid,
ai=self.ai_grid,
)

print("done.")


def read_model(datadir=""):
"""
Expand All @@ -695,6 +765,16 @@ def read_model(datadir=""):
if os.path.isfile(path):
d["T"] = _read_dust_temperature_dat(datadir=datadir)

path = os.path.join(datadir, "metadata.npz")
if os.path.isfile(path):
metadata = np.load(path)
d["grid"].a = metadata["a"]
d["grid"].ai = metadata["ai"]
if d["grid"].a.shape[0] != d["rho"].shape[-1]:
s = "Warning: " \
"The meta data file does not match the RADMC-3D setup."
print(s)

return SimpleNamespace(**d)


Expand Down

0 comments on commit 13d9bbc

Please sign in to comment.