Skip to content

Commit

Permalink
Added 3D stellar wind example to docs
Browse files Browse the repository at this point in the history
  • Loading branch information
@FredDeCeuster
FredDeCeuster committed Aug 8, 2024
1 parent 8f8946a commit a37a8db
Show file tree
Hide file tree
Showing 13 changed files with 11,856 additions and 20 deletions.
54 changes: 37 additions & 17 deletions docs/src/examples/1D_stellar_wind.ipynb
View file

Large diffs are not rendered by default.

2,476 changes: 2,476 additions & 0 deletions docs/src/examples/3D_stellar_wind.dump
View file

Large diffs are not rendered by default.

1,824 changes: 1,824 additions & 0 deletions docs/src/examples/3D_stellar_wind.in
View file

Large diffs are not rendered by default.

388 changes: 388 additions & 0 deletions docs/src/examples/3D_stellar_wind.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,388 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3D Stellar wind\n",
"\n",
"---\n",
"In this example, we consider an smoothed-particle hydrodynamics (SPH) model of a companion-perturbed stellar wind."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import os\n",
"from astropy import units, constants"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"---\n",
"\n",
"## Model setup\n",
"First we download a snapshot of this SPH model."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"setup_file = '3D_stellar_wind_data/wind.setup'\n",
"input_file = '3D_stellar_wind_data/wind.in'\n",
"dump_file = '3D_stellar_wind_data/wind.dump'"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--2024-08-08 16:03:19-- https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2024a/v10e00/wind.setup\n",
"Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.109.133, 185.199.108.133, 185.199.110.133, ...\n",
"Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.109.133|:443... connected.\n",
"HTTP request sent, awaiting response... 200 OK\n",
"Length: 1209 (1.2K) [text/plain]\n",
"Saving to: ‘3D_stellar_wind_data/wind.setup’\n",
"\n",
"3D_stellar_wind_dat 100%[===================>] 1.18K --.-KB/s in 0s \n",
"\n",
"2024-08-08 16:03:19 (30.7 MB/s) - ‘3D_stellar_wind_data/wind.setup’ saved [1209/1209]\n",
"\n",
"--2024-08-08 16:03:20-- https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2024a/v10e00/wind.in\n",
"Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.111.133, 185.199.110.133, 185.199.108.133, ...\n",
"Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.111.133|:443... connected.\n",
"HTTP request sent, awaiting response... 200 OK\n",
"Length: 5579 (5.4K) [text/plain]\n",
"Saving to: ‘3D_stellar_wind_data/wind.in’\n",
"\n",
"3D_stellar_wind_dat 100%[===================>] 5.45K --.-KB/s in 0s \n",
"\n",
"2024-08-08 16:03:20 (33.2 MB/s) - ‘3D_stellar_wind_data/wind.in’ saved [5579/5579]\n",
"\n",
"--2024-08-08 16:03:21-- https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2024a/v10e00/wind_v10e00\n",
"Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.111.133, 185.199.110.133, 185.199.108.133, ...\n",
"Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.111.133|:443... connected.\n",
"HTTP request sent, awaiting response... 200 OK\n",
"Length: 70236960 (67M) [application/octet-stream]\n",
"Saving to: ‘3D_stellar_wind_data/wind.dump’\n",
"\n",
"3D_stellar_wind_dat 100%[===================>] 66.98M 37.1MB/s in 1.8s \n",
"\n",
"2024-08-08 16:03:27 (37.1 MB/s) - ‘3D_stellar_wind_data/wind.dump’ saved [70236960/70236960]\n",
"\n"
]
}
],
"source": [
"!wget 'https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2024a/v10e00/wind.setup' --output-document $setup_file\n",
"!wget 'https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2024a/v10e00/wind.in' --output-document $input_file\n",
"!wget 'https://raw.githubusercontent.com/Ensor-code/phantom-models/main/Malfait%2B2024a/v10e00/wind_v10e00' --output-document $dump_file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use [plons](https://github.com/Ensor-code/plons) to open the data."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"import plons\n",
"\n",
"setupData = plons.LoadSetup(f'{os.getcwd()}/3D_stellar_wind_data', 'wind')\n",
"dumpData = plons.LoadFullDump(f'{os.getcwd()}/{dump_file}', setupData)\n",
"\n",
"position = dumpData[\"position\"]*1e-2 # position vectors [cm -> m]\n",
"velocity = dumpData[\"velocity\"]*1e3 # velocity vectors [km/s -> m/s]\n",
"rho = dumpData[\"rho\"] # density [g/cm^3]\n",
"tmp = dumpData[\"Tgas\"] # temperature [K]\n",
"tmp[tmp<2.7] = 2.7 # Cut-off temperatures below 2.7 K\n",
"\n",
"# Unpack velocities\n",
"v_x, v_y, v_z = velocity.T\n",
"\n",
"# Convert rho (total density) to H2 abundance\n",
"nH2 = rho * 1.0e+6 * constants.N_A.si.value / 2.02\n",
"\n",
"# Define turbulence at 150 m/s\n",
"trb = 150.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we map the particle data to a regular Cartesian mesh."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"from pomme.haar import Haar\n",
"\n",
"# Map point data to a regular grid\n",
"haar = Haar(position, q=8)\n",
"# Zoom in on the centre region to avoid edge effects\n",
"imin = 2**(haar.q-3)\n",
"imax = 3*imin\n",
"# Map data to a regular grid\n",
"nH2_dat = haar.map_data(nH2, interpolate=True)[-1][imin:imax,imin:imax,imin:imax]\n",
"tmp_dat = haar.map_data(tmp, interpolate=True)[-1][imin:imax,imin:imax,imin:imax]\n",
"v_x_dat = haar.map_data(v_x, interpolate=True)[-1][imin:imax,imin:imax,imin:imax]\n",
"v_y_dat = haar.map_data(v_y, interpolate=True)[-1][imin:imax,imin:imax,imin:imax]\n",
"v_z_dat = haar.map_data(v_z, interpolate=True)[-1][imin:imax,imin:imax,imin:imax]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"### TensorModel\n",
"With all the data in place, we can start building a pomme model.\n",
"First, we store all model parameters as a TensorModel object and store this in an HDF5 file.\n",
"We will use this later as the ground truth to verify our reconstructions against."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pomme.model import TensorModel\n",
"\n",
"model = TensorModel(shape=nH2_dat.shape, sizes=haar.xyz_L)\n",
"model['log_H2' ] = np.log(nH2_dat).astype(np.float64)\n",
"model['log_temperature' ] = np.log(tmp_dat).astype(np.float64)\n",
"model['velocity_x' ] = v_x_dat .astype(np.float64)\n",
"model['velocity_y' ] = v_y_dat .astype(np.float64)\n",
"model['velocity_z' ] = v_z_dat .astype(np.float64)\n",
"model['log_v_turbulence'] = np.log(trb)*np.ones(model.shape, dtype=np.float64)\n",
"model.save('3D_stellar_wind_truth.h5')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"### GeneralModel\n",
"First, we define the functions that can generate the model distributions from the model parameters."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import torch\n",
"from pomme.utils import planck, T_CMB\n",
"\n",
"def get_velocity(model):\n",
" \"\"\"\n",
" Get the velocity from the TensorModel.\n",
" \"\"\"\n",
" return model['velocity_z']\n",
"\n",
"def get_temperature(model):\n",
" \"\"\"\n",
" Get the temperature from the TensorModel.\n",
" \"\"\"\n",
" return torch.exp(model['log_temperature'])\n",
"\n",
"def get_abundance(model, l):\n",
" \"\"\"\n",
" Get the abundance from the TensorModel.\n",
" \"\"\"\n",
" # Define the assumed molecular fractions w.r.t. H2\n",
" X_mol = [3.0e-4, 5.0e-6]\n",
" # Return the abundance\n",
" return torch.exp(model['log_H2']) * X_mol[l]\n",
"\n",
"def get_turbulence(model):\n",
" \"\"\"\n",
" Get the turbulence from the TensorModel.\n",
" \"\"\"\n",
" return torch.exp(model['log_v_turbulence'])\n",
"\n",
"def get_boundary_condition(model, frequency):\n",
" \"\"\"\n",
" Get the boundary condition from the TensorModel.\n",
" model: TensorModel\n",
" The TensorModel object containing the model.\n",
" frequency: float\n",
" Frequency at which to evaluate the boundary condition.\n",
" \"\"\"\n",
" # Compute the incoming boundary intensity\n",
" Ibdy = torch.ones((model.shape[0], model.shape[1], len(frequency)), dtype=torch.float64)\n",
" Ibdy *= planck(temperature=T_CMB, frequency=frequency)\n",
" # Return the result\n",
" return Ibdy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using these functions, we can build a GeneralModel object that can be used to generate synthetic observations or reconstruct the required parameters. (Cfr. the SphericalModel class for spherically symmetric models.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pomme.model import GeneralModel\n",
"\n",
"gmodel = GeneralModel(model=model)\n",
"gmodel.get_velocity = get_velocity\n",
"gmodel.get_abundance = get_abundance\n",
"gmodel.get_turbulence = get_turbulence\n",
"gmodel.get_temperature = get_temperature\n",
"gmodel.get_boundary_condition = get_boundary_condition"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"### Spectral lines\n",
"We base our reconstructions on synthetic observations of two commonly observed rotational lines CO $J=4-3$ and SiO $J=3-2$. We explicitly provide the molar mass for SiO, since this is not extracted correctly from the line data file."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"You have selected line:\n",
" CO(J=4-3)\n",
"Please check the properties that were inferred:\n",
" Frequency 4.610407682e+11 Hz\n",
" Einstein A coeff 6.126000000e-06 1/s\n",
" Molar mass 28.0101 g/mol\n",
"You have selected line:\n",
" sio-h2(J=03-02)\n",
"Please check the properties that were inferred:\n",
" Frequency 1.302686830e+11 Hz\n",
" Einstein A coeff 1.058000000e-04 1/s\n",
" Molar mass 44.0849 g/mol\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/frederikd/.local/lib/python3.9/site-packages/astroquery/lamda/core.py:145: UserWarning: The first time a LAMDA function is called, it must assemble a list of valid molecules and URLs. This list will be cached so future operations will be faster.\n",
" warnings.warn(\"The first time a LAMDA function is called, it must \"\n"
]
}
],
"source": [
"from pomme.lines import Line\n",
"from pomme.utils import get_molar_mass\n",
"\n",
"lines = [\n",
" Line(species_name='CO', transition=3),\n",
" Line(species_name='sio-h2', transition=2, molar_mass=get_molar_mass('SiO'))\n",
"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"### Frequencies\n",
"Next, we define the velocity/frequency range.\n",
"We observe the lines in 100 frequency bins, centred around the lines, with a spacing of 120 m/s."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"vdiff = 120 # velocity increment size [m/s]\n",
"nfreq = 100 # number of frequencies\n",
"\n",
"velocities = nfreq * vdiff * torch.linspace(-1, +1, nfreq, dtype=torch.float64)\n",
"frequencies = [(1.0 + velocities / constants.c.si.value) * line.frequency for line in lines]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"## Synthetic observations\n",
"We can now generate synthetic observations, directly from the Model object.\n",
"We will use these later to derive our reconstructions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "magritte",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Loading

0 comments on commit a37a8db

Please sign in to comment.