Module diskchef.physics.physmodels

Classes

class ReadAndesPhysData (star_mass: Unit("solMass") = <Quantity 1. solMass>, xray_plasma_temperature: Unit("K") = <Quantity 10000000. K>, xray_luminosity: Unit("erg / s") = <Quantity 1.e+31 erg / s>, cr_padovani_use_l: bool = False, folder: Union[str, os.PathLike] = None, index: int = 0, time_mark: Unit("s") = <Quantity 0. s>, age: Unit("Myr") = <Quantity 1. Myr>, read_envelope: bool = False, dim: Literal['2d', '3d'] = '2d', coord_type: Literal['z', 'zr'] = 'zr')

Class to read a version of ANDES (accounting for jets + no chemistry) output for following usage in diskchef.maps In testing

Expand source code

class ReadAndesPhysData(PhysicsBase):
    """
    Class to read a version of ANDES (accounting for jets + no chemistry) output for following usage in `diskchef.maps`
    In testing
    """
    folder: PathLike = None
    index: int = 0
    time_mark: u.s = 0 * u.s
    age: u.Myr = 1 * u.Myr
    read_envelope: bool = False
    dim: Literal["2d", "3d"] = '2d'
    coord_type: Literal["z", "zr"] = 'zr'

    @property
    def table(self) -> CTable:
        return self._table

    def read(self, index: int = None) -> CTable:
        """Return the associated diskchef.CTable with r, z, and dust and gas properties"""
        if index is None:
            index = self.index
        reader = astropy.io.ascii.get_reader(Reader=astropy.io.ascii.CommentedHeader)
        reader.header.splitter.delimiter = '|'
        # if physics[["ir", "iz"]] != chemistry[["ir", "iz"]]:
        #     raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and chemistry tables do not match!")
        try:
            self._config = self._config_read(os.path.join(self.folder, "../config.ini"))
        except FileNotFoundError:
            self._config = self._config_read(os.path.join(self.folder, "../config_sep.ini"))
        self.age = float(self._config['last_time_mark  [Myr]']) * u.Myr
        if self.dim == '2d':
            try:
                if index == 0:
                    physics = reader.read(os.path.join(self.folder, f"physical_structure_{index:05d}"))
                else:
                    physics = reader.read(os.path.join(self.folder, f"physical_structure_{index - 1:05d}"))
            except FileNotFoundError:
                physics = reader.read(os.path.join(self.folder, f"physical_structure"))
            table = CTable(
                [
                    physics["Radius"] << u.au,
                    physics["Height"] << u.au,
                ],
                names=["Radius", "Height"]
            )
            table["Gas temperature"] = physics["Gas temperature"] << (u.K)
            table["Dust temperature"] = physics["Dust temperature"] << (u.K)
        if self.dim == '3d':
            physics = reader.read(os.path.join(self.folder, f"all_data_{index:05d}"))
            table = CTable(
                [
                    physics["Radius"] << u.au,
                    physics["Height"] << u.au,
                    physics["Angle"] << u.rad,
                ],
                names=["Radius", "Height", "Phi"], dims='3d'
            )
            table["Gas temperature"] = physics["Temperature"] << (u.K)
            table["Dust temperature"] = physics["Temperature"] << (u.K)

        if self.coord_type == 'zr':
            table["Height to radius"] = np.nan * u.dimensionless_unscaled
            for izr in physics["iz"]:
                indices = physics["iz"] == izr
                table["Height to radius"][indices] = np.nanmedian(
                    physics["Height"][indices] / physics["Radius"][indices])
                table["Height"] = table.zr * table.r

        if self.coord_type == 'z':
            table["Height to radius"] = table["Height"]/table.r

        table["Distance to star"] = np.sqrt(table.r ** 2 + table.z ** 2)
        table["Gas density"] = physics["Gas density"] << (u.g / u.cm ** 3)
        table["Dust density"] = physics["Dust density"] << (u.g / u.cm ** 3)


        if self.read_envelope:
            try:
                if self.dim == '2d':
                    envelope_mask_table = astropy.io.ascii.read(os.path.join(self.folder, f"Envelope_test1.dat"))
                    if np.any(envelope_mask_table[["ir", "iz"]] != physics[["ir", "iz"]]):
                        raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and envelope tables do not match!")
                    table["Envelope mask"] = envelope_mask_table['is_envelope']
                elif self.dim == '3d':
                    table["Envelope mask"] = physics["env"]
            except FileNotFoundError:
                self.logger.info("Envelope mask file not found. Assuming no envelope")
                table["Envelope mask"] = np.zeros(len(table['Gas density']))
            dust_types = np.sort(np.unique(table["Envelope mask"]))
            if len(dust_types) > 1:
                for dt in dust_types:
                    dtmask = table["Envelope mask"] == dt
                    table[f"Dust density {dt}"] = table[f"Dust density"]
                    table[f"Dust temperature {dt}"] = table[f"Dust temperature"]
                    table[f"Dust density {dt}"][~dtmask] = 0
                    table[f"Dust temperature {dt}"][~dtmask] = 0
        else:
            table["Envelope mask"] = np.zeros(len(table['Gas density']))

        dust_types = np.sort(np.unique(table["Envelope mask"]))
        if len(dust_types) > 1:
            for dt in dust_types:
                dtmask = table["Envelope mask"] == dt
                table[f"Dust density {dt}"] = table[f"Dust density"]
                table[f"Dust temperature {dt}"] = table[f"Dust temperature"]
                table[f"Dust density {dt}"][~dtmask] = 0
                table[f"Dust temperature {dt}"][~dtmask] = 0

        if self.dim == '2d':
            if self.coord_type == 'zr':
                table.sort(['Height to radius', 'Radius'])
            elif self.coord_type == 'z':
                table.sort(['Height', 'Radius'])
        elif self.dim == '3d':
            if self.coord_type == 'zr':
                table.sort(['Phi', 'Height to radius', 'Radius'])
            elif self.coord_type == 'z':
                table.sort(['Phi', 'Height', 'Radius'])


        try:
            self.physics = diskchef.physics.base.PhysicsBase(
                star_mass=float(self._config['stellar mass [MSun]']) * u.solMass,
            )
        except KeyError:
            self.physics = diskchef.physics.base.PhysicsBase(
                star_mass=(float(self._config['stellar mass 1 [MSun]'])+float(self._config['stellar mass 2 [MSun]'])) * u.solMass,
            )
        self.physics.table = table

        return table

    def __post_init__(self):
        super().__post_init__()
        self._table = self.read()

    def _config_read(self, path: PathLike) -> dict:
        out = {}
        with open(path, 'r') as configfile:
            for line in configfile.readlines():
                if len(line.split('::')) == 1:
                    continue
                value, key = [entry.strip() for entry in line.split('::')]
                out[key] = value
        return out

    def recalc_dust_dens(self):
        dust_types = np.sort(np.unique(self.table["Envelope mask"]))
        if len(dust_types) > 1:
            for dt in dust_types:
                dtmask = self.table["Envelope mask"] == dt
                self.table[f"Dust density {dt}"] = self.table[f"Dust density"]
                self.table[f"Dust temperature {dt}"] = self.table[f"Dust temperature"]
                self.table[f"Dust density {dt}"][~dtmask] = 0
                self.table[f"Dust temperature {dt}"][~dtmask] = 0

Ancestors

• PhysicsBase

Class variables

var age : Unit("Myr")

var coord_type : Literal['z', 'zr']

var dim : Literal['2d', '3d']

var folder : Union[str, os.PathLike]

var index : int

var read_envelope : bool

var time_mark : Unit("s")

Methods

def read(self, index: int = None) ‑> CTable

Return the associated diskchef.CTable with r, z, and dust and gas properties

def recalc_dust_dens(self)

Inherited members

• PhysicsBase:
  • column_density_to
  • cosmic_ray_padovani18
  • ionization
  • plot_column_densities
  • plot_density
  • plot_temperatures
  • table
  • xray_bruderer

class ReadGFARGO2Data (star_mass: Unit("solMass") = <Quantity 1. solMass>, xray_plasma_temperature: Unit("K") = <Quantity 10000000. K>, xray_luminosity: Unit("erg / s") = <Quantity 1.e+31 erg / s>, cr_padovani_use_l: bool = False, folder: Union[str, os.PathLike] = None, index: int = 0, zr_max: float = 1, nzr: int = 10, rho_grain: Unit("g / cm3") = <Quantity 3. g / cm3>)

Class to read Hurakan output for following usage in diskchef.maps Unfinished

Expand source code

class ReadGFARGO2Data(PhysicsBase):
    """
    Class to read Hurakan output for following usage in `diskchef.maps`
    Unfinished
    """
    folder: PathLike = None
    index: int = 0
    zr_max: float = 1
    nzr: int = 10
    rho_grain: u.g/u.cm**3 = 3 * u.g/u.cm**3

    @property
    def table(self) -> CTable:
        return self._table

    def read(self, index: int = None) -> CTable:
        """Return the associated diskchef.CTable with r, z, and dust and gas properties"""

        if index is None:
            index = self.index

        self.config = self._config_read()
        radii_edges = np.array(
            [self.config['Rmin'] * np.exp(i / self.config['Nrad'] * np.log(self.config['Rmax'] / self.config['Rmin'])) for i in range(self.config['Nrad'] + 1)])
        radii_points = np.sqrt(radii_edges[1:]*radii_edges[:-1])
        phi_edges = np.array([2 * np.pi * i / self.config['Nsec'] for i in range(self.config['Nsec'] + 1)])
        phi_points = 0.5*(phi_edges[1:] + phi_edges[:-1])
        zr_edges = np.linspace(0, self.zr_max, self.nzr+1)
        zr_points = 0.5*(zr_edges[1:]+zr_edges[:-1])

        ZR, R, PHI = np.meshgrid(zr_points, radii_points, phi_points, indexing='ij')
        R = R * u.au
        Z = ZR * R
        PHI = PHI * u.rad
        R_2D, PHI_2D = np.meshgrid(radii_points, phi_points, indexing='ij')

        gas_surf_dens = np.fromfile(os.path.join(self.folder, f'gas_dens_{index:06d}.dat'))
        gas_surf_dens_2D = np.reshape(gas_surf_dens, (self.config['Nrad'], self.config['Nsec']))
        self.gas_surf_dens = gas_surf_dens_2D
        gas_surf_dens = np.repeat(gas_surf_dens_2D[np.newaxis, :, :], self.nzr, axis=0) * u.Msun / u.au**2
        ndust = len(self.config['DustSizeBin'])
        dust_surf_dens = []

        nstokes = []
        self.dust_sizes = []
        self.dust_surf_dens = []
        for i in range(ndust):
            temp_data = np.fromfile(os.path.join(self.folder, f'dust_dens_s{i}_{index:06d}.dat'))
            temp_data = np.reshape(temp_data, (self.config['Nrad'], self.config['Nsec']))
            self.dust_surf_dens.append(temp_data)
            temp_data = np.repeat(temp_data[np.newaxis, :, :], self.nzr, axis=0) * u.Msun / u.au**2
            dust_surf_dens.append(temp_data)
            if self.config['DustConstStokes'] == 'YES':
                nstokes.append(np.zeros((self.nzr, self.config['Nrad'], self.config['Nsec'])) + self.config['DustSizeBin'][i])
                dsz_distro = (2/np.pi * self.config['DustSizeBin'][i] * gas_surf_dens / self.rho_grain[i]).to(u.cm)
                self.dust_sizes.append(np.max(dsz_distro))
            else:
                self.dust_sizes.append(self.config['DustSizeBin'][i] * u.cm)
                nst = (np.pi/2 * self.rho_grain[i] * self.config['DustSizeBin'][i] * u.cm / gas_surf_dens).to(u.dimensionless_unscaled).value
                nstokes.append(nst)


        if self.config['Adiabatic'] == 'YES':
            gas_energy = np.fromfile(os.path.join(self.folder, f'gas_energy_{index:06d}.dat'))
            gas_energy_2D = np.reshape(gas_energy, (self.config['Nrad'], self.config['Nsec']))
            gas_energy = np.repeat(gas_energy_2D[np.newaxis, :, :], self.nzr, axis=0) * u.Msun / u.au**3 # check the real units !!


            cs_sq = (self.config['AdiabaticIndex']-1)*gas_energy/gas_surf_dens

            cs_sq_dust = []
            for i in range(ndust):
                cs_sq_dust.append(cs_sq * (self.config['AlphaViscosity']/(self.config['AlphaViscosity']+nstokes[i])))
        if self.config['Adiabatic'] == 'NO':
            raise CHEFNotImplementedError("Not written for now")

        # gas_temperature = 2.35 * c.m_p / c.k_B * cs_sq

        self.gas_aspect_ratio = np.sqrt(cs_sq*R).to(u.dimensionless_unscaled).value
        self.dust_aspect_ratio = [np.sqrt(cs_sq_dust[i]*R).to(u.dimensionless_unscaled).value for i in range(ndust)]

        f_gas = np.exp(-1/cs_sq*(1/R - 1/np.sqrt(R**2 + Z**2)))
        f_int = 2*np.trapz(f_gas, x=Z, axis=0)
        f_int = np.repeat(f_int[np.newaxis, :, :], self.nzr, axis=0)
        gas_dens = gas_surf_dens * f_gas / f_int

        dust_dens = []
        for i in range(ndust):
            f = np.exp(-1 / cs_sq_dust[i] * (1 / R - 1 / np.sqrt(R ** 2 + Z ** 2)))
            f_int = 2*np.trapz(f, x=Z, axis=0)
            f_int = np.repeat(f_int[np.newaxis, :, :], self.nzr, axis=0)
            dust_dens.append(dust_surf_dens[i]*f/f_int)

        table = CTable(
            [
                R.flatten(), Z.flatten(), PHI.flatten(),
            ],
            names=["Radius", "Height", 'Phi'], coord_base='cylinder', dims='3d',
            edges=[radii_edges*u.au, zr_edges, phi_edges*u.rad]
        )

        table['Height to radius'] = ZR.flatten()
        table['Distance to star'] = np.sqrt(table['Radius']**2 + table['Height']**2)
        table['Gas density'] = gas_dens.flatten()


        table.dust_list = [
            DustPopulationOptool(opacity_file=f'{self.folder}/dustkappa_{self.dust_sizes[i].value:.2f}.inp', table=table,
                                 name=f'{self.dust_sizes[i].value:.2f}', size=self.dust_sizes[i]) for i in range(len(self.dust_sizes))]

        for i in range(ndust):
            table[f'{table.dust_list[i].name} density'] = dust_dens[i].flatten()

        # if not self.isotherm:
        #     table['Gas temperature'] = gas_temperature.flatten()
        #     for i in range(ndust):
        #         table[f'Dust temperature {i}'] = gas_temperature.flatten()

        table.sort(['Phi', 'Height to radius', 'Radius'])


        return table

    def _config_read(self):
        cfg = np.loadtxt(os.path.join(self.folder, 'disk.cfg'), dtype='str', usecols=(0, 1))
        cfg = {key: value for key, value in cfg}
        for key in cfg:
            try:
                if ',' in cfg[key]:
                    cfg[key] = cfg[key].split(',')
                if 'N' in key:
                    cfg[key] = int(cfg[key])
                else:
                    cfg[key] = np.float64(cfg[key])
            except ValueError:
                continue

        if type(cfg['DustSizeBin']) == np.float64:
            cfg['DustSizeBin'] = [cfg['DustSizeBin']]
        if type(cfg['DustMassBin']) == np.float64:
            cfg['DustMassBin'] = [cfg['DustMassBin']]

        return cfg


    def __post_init__(self):
        super().__post_init__()
        self._table = self.read()

Ancestors

• PhysicsBase

Class variables

var folder : Union[str, os.PathLike]

var index : int

var nzr : int

var rho_grain : Unit("g / cm3")

var zr_max : float

Methods

def read(self, index: int = None) ‑> CTable

Return the associated diskchef.CTable with r, z, and dust and gas properties

Inherited members

• PhysicsBase:
  • column_density_to
  • cosmic_ray_padovani18
  • ionization
  • plot_column_densities
  • plot_density
  • plot_temperatures
  • table
  • xray_bruderer

class ReadHurakanData (star_mass: Unit("solMass") = <Quantity 1. solMass>, xray_plasma_temperature: Unit("K") = <Quantity 10000000. K>, xray_luminosity: Unit("erg / s") = <Quantity 1.e+31 erg / s>, cr_padovani_use_l: bool = False, folder: Union[str, os.PathLike] = None, index: int = 0, time_mark: Unit("s") = <Quantity 0. s>, age: Unit("Myr") = <Quantity 1. Myr>, read_envelope: bool = False)

Class to read Hurakan output for following usage in diskchef.maps Unfinished

Expand source code

class ReadHurakanData(PhysicsBase):
    """
    Class to read Hurakan output for following usage in `diskchef.maps`
    Unfinished
    """
    folder: PathLike = None
    index: int = 0
    time_mark: u.s = 0 * u.s
    age: u.Myr = 1 * u.Myr
    read_envelope: bool = False

    @property
    def table(self) -> CTable:
        return self._table

    def read(self, index: int = None) -> CTable:
        """Return the associated diskchef.CTable with r, z, and dust and gas properties"""
        if index is None:
            index = self.index
        physics = astropy.io.ascii.read(os.path.join(self.folder, f"result-{index:04d}.dat"), header_start=3, data_start=4)#, data_end=-101)
        theta_zeros = physics["Theta"] == 0
        physics["Theta"][theta_zeros] = np.abs(np.pi - physics["Theta"][-1])
        # physics = physics[~theta_zeros]
        # Will need fixing !!!
        table = CTable(
            [
                physics["R"] << u.cm,
                physics["Theta"] << u.rad,
            ],
            names=["Distance to star", "Theta"], coord_base='sphere', coord_style='edge'
        )
        self._config = astropy.io.ascii.read(os.path.join(self.folder, f"result-{index:04d}.dat"), data_end=3)
        self.time_mark = self._config['col2'][1] * u.s
        table["Height"] = table.d * np.cos(table.theta)
        table["Radius"] = table.d * np.sin(table.theta)
        table["Height to radius"] = np.tan(np.pi/2*u.rad - table.theta)

        table["Gas density"] = physics["Rho"] << (u.g / u.cm ** 3)
        table["Dust density"] = table["Gas density"]/100
        table["Gas temperature"] = physics["T"] << (u.K)
        table["Dust temperature"] = table["Gas temperature"]

        if self.read_envelope: #to be done
            try:
                envelope_mask_table = astropy.io.ascii.read(os.path.join(self.folder, f"envelope_{0:05d}"))
                if np.any(envelope_mask_table[["ir", "iz"]] != physics[["ir", "iz"]]):
                    raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and envelope tables do not match!")
                table["Envelope mask"] = envelope_mask_table['is_envelope']
            except FileNotFoundError:
                self.logger.info("Envelope mask file not found. Assuming no envelope")
                table["Envelope mask"] = np.zeros(len(table['Gas density']))
        else:
            table["Envelope mask"] = np.zeros(len(table['Gas density']))

        table.sort(['Height to radius', 'Radius'])

        self.physics = diskchef.physics.base.PhysicsBase(
            star_mass=float(self._config['col2'][2]) * u.g,
        )
        self.physics.table = table

        return table

    def __post_init__(self):
        super().__post_init__()
        self._table = self.read()

Ancestors

• PhysicsBase

Class variables

var age : Unit("Myr")

var folder : Union[str, os.PathLike]

var index : int

var read_envelope : bool

var time_mark : Unit("s")

Methods

def read(self, index: int = None) ‑> CTable

Return the associated diskchef.CTable with r, z, and dust and gas properties

Inherited members

• PhysicsBase:
  • column_density_to
  • cosmic_ray_padovani18
  • ionization
  • plot_column_densities
  • plot_density
  • plot_temperatures
  • table
  • xray_bruderer