Commit 59303a94 authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Added wiring for defining a power-law in dust temperature.

parent f777e9d9

etl/assets/commons/__init__.py

+14 −2
Original line number Diff line number Diff line
@@ -5,6 +5,7 @@ import astropy.units 
import sqlalchemy

import yaml

import numpy as np

import shutil

import os

import urllib.request

import base64
@@ -280,17 +281,23 @@ def get_physical_px_size(grid_metadata: dict) -> List[u.Quantity]: 


def get_moldata(species_names: list,

                logger: logging.Logger,

                path: Union[str, None] = None):

                path: Union[str, None] = None,

                use_cache: bool = True):

    """

    Downloads the molecular data from the Leiden molecular database; check whether the molecule in mapped in

        leiden_url_mapping

    :param species_names: the names of the species for which to download data

    :param logger: logger for output printing

    :param path: path where the files must be saved

    :param use_cache: use molecular input files in cache, if possible

    """

    _path = validate_parameter(path, default=os.path.join('mdl', 'radmc_files'))

    for species in species_names:

        if not os.path.isfile(os.path.join(_path, f'molecule_{species}.inp')):

        molecular_file_name = f'molecule_{species}.inp'

        if use_cache is True:

            shutil.copy(os.path.join('mdl', 'data', molecular_file_name),

                        os.path.join(_path, molecular_file_name))

        if not os.path.isfile(os.path.join(_path, molecular_file_name)):

            logger.info(f'Downloading file molecule_{species}.inp...')

            data = urllib.request.urlopen(leiden_url_mapping[species]).read().decode()

            with open(os.path.join(_path, f'molecule_{species}.inp'), 'w') as outfile:
@@ -362,6 +369,11 @@ def convert_frequency_to_wavelength(frequency: astropy.units.Quantity, 




def get_pg_engine(logger: logging.Logger) -> sqlalchemy.engine.Engine:

    """

    Return the SQLAlchemy engine, given the credentials in the external file

    :param logger: the logger to use

    :return: the SQLAlchemy engine

    """

    credentials = get_credentials(logger=logger,

                                  credentials_filename=os.path.join('credentials', 'db_credentials.yml'))

    url = f"postgresql://{credentials['DB_USER']}:{credentials['DB_PASS']}@{credentials['DB_HOST']}:{credentials['DB_PORT']}/{credentials['DB_NAME']}"

etl/main.py

+6 −3
Original line number Diff line number Diff line
@@ -13,11 +13,11 @@ from prs.prs_compute_integrated_fluxes_and_ratios import main as prs_main 
from prs.prs_inspect_results import main as prs_inspection_main





def compute_grid(tdust, nh2, line, density_keyword):

def compute_grid(tdust, nh2, line, density_keyword, dust_temperature_keyword):

    scratch_dir = os.path.join('mdl', 'scratches', str(uuid.uuid4()))

    overrides = {

        'grid': {

            'dust_temperature': tdust,

            dust_temperature_keyword: tdust,

            density_keyword: nh2,

        }

    }
@@ -43,7 +43,9 @@ def build_model_grid(run_id: str, 
                     cleanup_scratches: bool = True):

    stg_config = load_config_file(os.path.join('stg', 'config', 'config.yml'))

    pl_density_index = float(stg_config['grid']['density_powerlaw_idx'])

    pl_dust_temperature_idx = float(stg_config['grid']['dust_temperature_powerlaw_idx'])

    _model_type = 'spherical' if pl_density_index != 0 else 'homogeneous'

    _tdust_model_type = 'heated' if pl_dust_temperature_idx != 0 else 'isothermal'



    config = load_config_file(os.path.join('config', 'config.yml'))
@@ -56,8 +58,9 @@ def build_model_grid(run_id: str, 
    line_set = set(chain.from_iterable(line_pairs))



    density_keyword = 'central_density' if _model_type == 'homogeneous' else 'density_at_reference'

    dust_temperature_keyword = 'dust_temperature' if _model_type == 'isothermal' else 'dust_temperature_at_reference'



    parallel_args = product(dust_temperatures, densities, line_set, [density_keyword])

    parallel_args = product(dust_temperatures, densities, line_set, [density_keyword], [dust_temperature_keyword])

    with Pool(8) as pool:

        results = pool.starmap(compute_grid, parallel_args)

etl/stg/config/config.yml

+4 −3
Original line number Diff line number Diff line
@@ -9,9 +9,10 @@ grid: 
    maximum_radius_unit: 'pc'

    distance_reference: 0.5

    distance_reference_unit: 'pc'

    dust_temperature: 15

    dust_temperature: 2000

    dust_temperature_at_reference: 30

    dust_temperature_unit: 'K'

    dust_temperature_powerlaw_idx: 0

    dust_temperature_powerlaw_idx: -0.5

    microturbulence: 1.5

    microturbulence_unit: 'km/s'

    dim1: {"size": 2, "size_units": "pc", "shape": 21, "refpix": 10}
@@ -22,7 +23,7 @@ grid: 
lines:

    species_to_include: ['e-ch3oh']

    molecular_abundances: {

                              "e-ch3oh": 1e-10,

                              "e-ch3oh": 1e-8,

                              "p-h2": 1,

    }

    lines_mode: 'lvg'

etl/stg/stg_radmc_input_generator.py

+28 −16
Original line number Diff line number Diff line
@@ -121,19 +121,12 @@ def get_solid_body_rotation_y(grid_metadata): 




def get_profiles(grid_metadata: dict) -> dict:

    try:

        density_ref = {

            'value_at_reference': (

                (float(grid_metadata['density_at_reference']) * u.Unit(grid_metadata['density_unit']))

                .to(u.cm ** -3)).value,

            'distance_reference': (float(grid_metadata['distance_reference']) * u.Unit(

                grid_metadata['distance_reference_unit'])).to(u.cm).value

        }

    except (KeyError, AttributeError):

        density_ref = {

            'value_at_reference': None,

            'distance_reference': 1.0

        }

    density_ref = get_reference_value(grid_metadata=grid_metadata,

                                      quantity_name='density',

                                      desired_unit='cm^-3')

    tdust_ref = get_reference_value(grid_metadata=grid_metadata,

                                    quantity_name='dust_temperature',

                                    desired_unit='K')

    profiles_mapping = {

        'gas_number_density': {

            'central_value': (float(grid_metadata['central_density']) * u.Unit(grid_metadata['density_unit']))
@@ -144,9 +137,9 @@ def get_profiles(grid_metadata: dict) -> dict: 
        },

        'dust_temperature': {

            'central_value': float(grid_metadata['dust_temperature']),

            'power_law_index': 0,

            'value_at_reference': None,

            'distance_reference': 1.0

            'power_law_index': float(grid_metadata['density_powerlaw_idx']),

            'value_at_reference': tdust_ref['value_at_reference'],

            'distance_reference': tdust_ref['distance_reference']

        },

        'microturbulence': {

            'central_value': (float(grid_metadata['microturbulence']) * u.Unit(grid_metadata['microturbulence_unit']))
@@ -203,6 +196,25 @@ def get_profiles(grid_metadata: dict) -> dict: 
    return profiles





def get_reference_value(grid_metadata: dict,

                        quantity_name: str,

                        desired_unit: str) -> dict:

    try:

        reference_value_dict = {

            'value_at_reference': (

                (float(grid_metadata[f'{quantity_name}_at_reference']) * u.Unit(grid_metadata[f'{quantity_name}_unit']))

                .to(u.Unit(desired_unit))).value,

            'distance_reference': (float(grid_metadata['distance_reference']) * u.Unit(

                grid_metadata['distance_reference_unit'])).to(u.cm).value

        }

    except (KeyError, AttributeError):

        reference_value_dict = {

            'value_at_reference': None,

            'distance_reference': 1.0

        }

    return reference_value_dict





def write_grid_input_files(profiles: dict,

                           grid_metadata: dict,

                           path: Union[str, None] = None):

etl/tests/test_commons.py

+34 −1
Original line number Diff line number Diff line
@@ -13,7 +13,8 @@ from assets.commons import (compute_power_law_radial_profile, 
                            validate_parameter,

                            parse_grid_overrides,

                            get_credentials,

                            setup_logger)

                            setup_logger,

                            get_moldata)

from unittest import TestCase
@@ -261,3 +262,35 @@ class Test(TestCase): 
        self.assertDictEqual(credentials, expected_results)

        for key in expected_results:

            self.assertEqual(str(expected_results[key]), os.getenv(key))



    def test_get_moldata_from_cache(self):

        species = 'e-ch3oh'

        os.chdir('..')

        expected_path = os.path.join('tests', 'test_files', f'molecule_{species}.inp')

        try:

            os.remove(expected_path)

        except IOError:

            pass

        get_moldata(species_names=[species],

                    logger=self.logger,

                    path=os.path.join('tests', 'test_files'),

                    use_cache=True)

        self.assertTrue(os.path.isfile(expected_path))

        os.remove(expected_path)

        os.chdir('tests')



    def test_get_moldata_negate_cache(self):

        species = 'e-ch3oh'

        os.chdir('..')

        expected_path = os.path.join('tests', 'test_files', f'molecule_{species}.inp')

        try:

            os.remove(expected_path)

        except IOError:

            pass

        get_moldata(species_names=[species],

                    logger=self.logger,

                    path=os.path.join('tests', 'test_files'),

                    use_cache=False)

        self.assertTrue(os.path.isfile(expected_path))

        os.remove(expected_path)

        os.chdir('tests')