Skip to content

Commit

Permalink
Merge branch 'main' of github.com:lsiess/SPH_plotting_pipeline
Browse files Browse the repository at this point in the history
  • Loading branch information
@MatsEsseldeurs
MatsEsseldeurs committed Oct 17, 2024
2 parents e5a1791 + 93bc0af commit f3ed027
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 166 deletions.
158 changes: 26 additions & 132 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,151 +9,45 @@
![Build status](https://readthedocs.org/projects/plons/badge/?version=latest)


This is the README for the PLONS PLOtting tool for Nice Simulations. It can be used to
read in and analyse data of hydrodynamical simulations with the SPH code [PHANTOM](https://phantomsph.bitbucket.io/)
([Price et al. 2018](https://ui.adsabs.harvard.edu/abs/2018PASA...35...31P/abstract)). PHANTOM returns
(different types of useful data files: 'wind_xxxxx'-files and '.ev'-files. The former can
also be visualised by the SPH visualisation tool [SPLASH](https://users.monash.edu.au/~dprice/splash/) ([Price 2007](https://adsabs.harvard.edu/abs/2007PASA...24..159P)), and can be analysed more toroughly with PLONS.
This is the README for the PLONS PLOtting tool for Nice Simulations.
PLONS can be used to read in and analyse data of hydrodynamical simulations with the SPH code
[PHANTOM](https://phantomsph.bitbucket.io/) ([Price et al. 2018](https://ui.adsabs.harvard.edu/abs/2018PASA...35...31P/abstract)).
This python package is currently tailored to single, binary and triple AGB wind models.
For a more detailed explanation of how to use PLONS, see ([readthedocs](https://plons.readthedocs.io/en/latest/)).

This tool is complementary to the SPH visualisation tool [SPLASH](https://users.monash.edu.au/~dprice/splash/) ([Price 2007](https://adsabs.harvard.edu/abs/2007PASA...24..159P)), but can be analysed more toroughly with PLONS.

PHANTOM
-------

PHANTOM returns (full) dump files every certain timestep during the evolution of the
model. These are the wind_xxxxx-files and contain the relevant data (position, velocity,
mass, density, energy, temperature) of the SPH particles in the model. The last x rows of
the output are the x sink particles in the model. More detailed information of the sink
particles (position, velocity, mass, accreted mass, spin, angular momentum) as a function of
evolution time of the model, can be retrieved in the .ev files.

PLONS - General
PLONS pipeline
------------------

This pipeline is suited for single, binary and triple AGB wind models.

The following info can be attained:

1. 2D slice plots of the global structure of the last dump full of the model.

2. 1D line plots (radial structure) of the global structure of the last dump of the model along the x-, y- and z-axes.

3. Information about the velocity related quantities of the model + Quantitative measurement of the degree of aspherical morphology: morphological parameters eta, Qp and epsilon.

4. Cummulative mass fraction as a function of the polar coordinate theta.

5. Information of the orbital evolution.

6. 1D spherical profiles for single star models



Pipeline - More detailed information
------------------------------------

This pipeline loads the wind.in and wind.setup files first. This is the general setup information
of the model (configuration of the system, information of the thermodynamics, evolution time, wind
outflow setup,...). Next the last full dump is loaded, and using this data, some other useful
quantities are calculated (gas pressure/temperature, speed, sperical coordinates, sound speed of
the gas,...). Lastly, the data of the sink particle(s) is loaded. From this file, also the period
and orbital velocity of the sink particles is calculated, if the system is a multiple, as a function of
time. The last entry in the sink data, corresponds to the full dump loaded. Therefore the
position/velocity/mass/... of this last entry are also saved in the full dump data for easy usage.


<b>(1) 2D slice plots of the global structure of the last dump full of the model.</b>

The 2D slices in the plots are infinitely thin and calculated by using the smoothing kernel from
PHANTOM. The density, speed and temperature are calculated by this method on a grid of n*n
pixels, using the data of the m nearest neighbours (SPH particles) of those grid points. These
plots give a visual representation of the morphology present in the AGB wind of the simulation.


<b>(2) 1D line plots (radial structure) of the global structure of the last dump of the model along the x-, y- and z-axes.</b>
PLONS contains a pipeline that produces useful output files with information about your simulation, and currently contains the following functionalities:

In 1D, the data along the x-, y- and z-axes is constructed in the same way as explained in (1),
but for a 1D grid. This results in plots of the radial structure, which is an important addition
to the 2D slice plots.
1. 2D slice plots showing the density, temperature and velocity distribution of the last full dump of the model.​

2. 1D line plots (radial structure) of the global structure of the last dump of the model along the x-, y- and z-axes.​

<b>(3) Information about the velocity related quantities of the model + Quantitative measurement of the degree of aspherical morphology: morphological parameters eta, Qp and epsilon.</b>
3. Information and plots of the orbital evolution​

The terminal velocity in PHANTOM SPH code is not an input parameter, so it needs to be calculated.
Due to the morphology in the outflow of the AGB star, at a certain radius a wide range of speeds
will be obtained. Therefore calculating only one value for the terminal velocity is not possible.
Instead, three values are obtained, a minimum, mean and maximum terminal velocity. The method used
here is binning the speed as a function of the radius of the SHP particle to the AGB star. Bins of
1AU are constructed and per bin, the minimum, mean and maximum speed is calculated. Using speed of
the outer 20% of the data, three values for the terminal velocity are obtained.
4. 1D spherical profiles, which are useful for single star models

NOTE: For some models (most often with slow initial wind velocity, low mass loss rate and/or low-mass
companion), the interaction of the companion has not reached the outer boundary of the model. Therefore,
this method will not lead to a correct terminal velocity value and thus the unphysical part of the data
will be cut off. Visually, the speed profile as a function of radius shows a strong, unphysical linear
increase, starting with a clear twist.
To run this pipeline, use the script 'main.py' in plons/scripts/ .
To select a certain model and option, you can execute this script as python main.py -m ... -o ... (and filling in the model and option number instead of the ...).

Source file functions and example notebooks
----------
Next to the functionalities included in the pipeline, you can use Plons for more detailed analysises using the functions defined in the source files (in plons/src/plons/).
Examples of how to use these functions can be found in some example notebooks (in docs/src/1_examples/), describing a more thorough analysis that makes use of line slices and sliceplots (/0_creating_figures_package), and an analysis of accretion disks around the companion (/2_accretion_disk).

In order to get some quantitative indication about the morphology of the models, several morphology
parameters are currently in use:

- $\eta = v/v_{orb}$ (see [Saladino et al. 2018](https://ui.adsabs.harvard.edu/abs/2018A%26A...618A..50S/abstract); [El Mellah et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...637A..91E/abstract))

- $Q_p = p_{comp}/p_{wind} = (M_{comp} v_{orb}) / (M_{wind} v_{wind})$
(see [Decin, L. et al. 2020](https://ui.adsabs.harvard.edu/abs/2020Sci...369.1497D/abstract))
- $\epsilon = e_{kin}/e_{grav} = (v_{wind}^2 a)/(24 G^3 M_{comp}^2 M_{AGB})^{1/3}$ (see [Maes et al. 2021](https://ui.adsabs.harvard.edu/abs/2021A%26A...652A..51M/abstract); [Malfait et al. 2021](https://ui.adsabs.harvard.edu/abs/2021A%26A...653A..25M/abstract))
- $R_{capt}/a = (2 G M_{comp})/(v_{wind}^2)/a$

Depending on the value of these parameters, the model is expected to show radial/EDE/complex morphology.

$$R_{hill} = a(1-e)(M_{comp}/ 3M_{AGB})^{1/3}$$
Therefore, Rhill/Rcapt = epsilon, if for both radii the same v_wind is used.

VERY IMPORTANT NOTE:
The parameter '$v_{wind}$' is not unambiguously defined. Even more, if only a binary models is used,
this parameter is most likely not to be constrained properly. For the different morphological parameters,
multiple values for $v_{wind}$ are used.
For eta, different velocities/speed $v$ are used:
- terminal velocity
- speed of the wind at the location of the companion

For Qp and epsilon, two different values are used as v_wind:
- the mean velocity of the wind at the location is used as calculated from the binning
- the average of the min and max velocity at the location of the companion is used

For $R_{capt}/a$ the initial wind velocity is used to get a rough indication.

Better options for $v_{wind}$ would be to use the velocity of the corresponding single model as follows:
- speed of the wind at the location of the companion
- speed resulting from the vector sum of the speed of the wind at the location of the companion
and the orbital velocity of the AGB star.

These can be calculated using the output of this pipeline.


<b>(4) Cummulative mass fraction as a function of the polar coordinate theta.</b>

The cummulative mass fraction (CMF) is calculated as a function of the polar angle theta (theta = 0.5pi
is the orbital plane, theta = 0 the north pole), again by using the method of binning. Only the northern
hemisphere is used, because of symmetry. Because of asymmetry along the x-axis, the positive (right)
and negative (left) part according to x are also calculated seperately. Certainly for eccentric models
this can differ. To calculate the CMF, only the selfsimilar part of the outflow it used, therefore a
factor (called factor) is given. The inner factor * semi-major axis AU is left out in this calculation.

From this CMF, the angle is calculated where 25, 50 and 75% of the mass of the wind is present, called
theta25/50/75 respectively. From this, the parameter delta is calculated, defined as
(theta25-theta50)/(theta50-theta75). When this parameter is normalised to the value for the corresponding
single model, it gives a measure for the degree of EDE.

Also the mean density is calculated using a similar approach as the CMF, using the method of binning.
This mean density can easily be normalised to compare the morphology in different models.

<b>(5) Information of the orbital evolution.</b>

In order to get information about the orbital evolution, the sink files are used (.ev), which gives
relevant quantities of the sink particles as a function of time. This file returns plots of the orbit,
evolution of orbital separation, orbital velocity and accreted mass of the companion.

Reading in PHANTOM data
-------

<b>(6) 1D spherical profiles for single star models</b>
PHANTOM returns different types of useful data files.
Every certain timestep during the evolution of the model, dump files are output, typically named e.g. as 'wind_xxxxx'.
These files contain relevant data (such as position, velocity,
mass, density, energy, temperature) of the SPH particles in the model at a certain timestep.
The last 'x' rows in these datafiles correspond to the 'x' sink particles in the model.
Further, for every run, evolutionary '*.ev' are produced that contain information about the time evolution of usefull quantities such as the mass of a sink particle for example.

In order to check your simulations and compare the 3D simulations to a 1D solution, this routine reads wind_1D.dat to plot the 1D solution, as well as the 3D sph values.

71 changes: 37 additions & 34 deletions scripts/main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,11 +23,11 @@

options = { '1': '(1) 2D slice plots',
'2': '(2) 1D line plots & tube plots',
'3': '(3) velocity related quantities and Morphological parameters',
'4': '(4) Cummulative mass fraction',
'5': '(5) Orbital evolution',
'6': '(6) 1D spherical model',
'7': '(7) Accretion disk'
# '3': '(3) velocity related quantities and Morphological parameters',
# '4': '(4) Cummulative mass fraction',
'3': '(3) Orbital evolution',
'4': '(4) 1D spherical model'
# '7': '(7) Accretion disk'
# '7': '(7) ArchimedianSpiral'
}

Expand Down Expand Up @@ -94,42 +94,45 @@ def runPart(part, run, saveloc, dumpData, setup, sinkData, outerData):
rs.radialStructPlots(run, saveloc, dumpData, setup)
rs.plotRho_r2(run,saveloc, dumpData, setup)

if part == '3':
print('')
print('(3) Start calculations for terminal velocity...')
tmv.main_terminalVelocity(setup, dumpData, sinkData, saveloc, run)
# if part == '3':
# print('')
# print('(3) Start calculations for terminal velocity...')
# tmv.main_terminalVelocity(setup, dumpData, sinkData, saveloc, run)

if part == '4':
print('')
print('(4) Start calculations for the cummulative mass fraction and mean density plots...')
if setup['single_star'] == True:
cmf.CMF_meanRho(run, saveloc, dumpData, setup, factor)
else:
if factor>0:
cmf.CMF_meanRho(run, saveloc, outerData, setup, factor)
else:
cmf.CMF_meanRho(run, saveloc, dumpData, setup, factor)
# if part == '4':
# print('')
# print('(4) Start calculations for the cummulative mass fraction and mean density plots...')
# if setup['single_star'] == True:
# cmf.CMF_meanRho(run, saveloc, dumpData, setup, factor)
# else:
# if factor>0:
# cmf.CMF_meanRho(run, saveloc, outerData, setup, factor)
# else:
# cmf.CMF_meanRho(run, saveloc, dumpData, setup, factor)

if part == '5':
if part == '3':
print('')
if setup['single_star']:
print('(5) A single model has no orbit, and thereby no orbital evolution.')
print('(3) A single model has no orbit, and thereby no orbital evolution.')
print(' The orbital evolution part is therefore skipped.')
else:
print('(5) Start calculations for orbital evolution...')
print('(3) Start calculations for orbital evolution...')

# ov.orbEv_main(run, saveloc, sinkData, setup)
ov.orbEv_main(run, saveloc, sinkData, setup, dumpData)

if part == '6':
if part == '4':
print('')
print('(6) Start calculating for the 1D spherical plots')

print('(4) Start calculating for the 1D spherical plots')
dp.profiles_main(run, loc, saveloc, dumpData, setup)

if part == '7':
print('')
print('(7) Start calculations for accretion disk analysis')
acd.accrDiskAnalysis(run, saveloc, dumpData, setup)


# if part == '7':
# print('')
# print('(7) Start calculations for accretion disk analysis')
# acd.accrDiskAnalysis(run, saveloc, dumpData, setup)

# if part == '7':
# print('')
Expand Down Expand Up @@ -228,13 +231,13 @@ def LoadData_cgs(run, loc, userSettings, bound = None, factor = -1, number = -1,
print()
print(' (1) 2D slice plots of the global structure of the last dump full of the model.')
print(' (2) 1D line plots (radial structure) of the global structure of the last dump of the model along the x-, y- and z-axes.')
print(' (3) Information about the velocity related quantities of the model + Quantitative measurement of the degree of aspherical morphology: morphological parameters eta, Qp and epsilon.')
print(' (4) Cummulative mass fraction as a function of the polar coordinate theta.')
print(' (5) Information of the orbital evolution.')
print(' (6) 1D spherical profiles for single star models.')
print(' (7) Analysis of accretion disk')
# print(' (3) Information about the velocity related quantities of the model + Quantitative measurement of the degree of aspherical morphology: morphological parameters eta, Qp and epsilon.')
# print(' (4) Cummulative mass fraction as a function of the polar coordinate theta.')
print(' (3) Information of the orbital evolution.')
print(' (4) 1D spherical profiles for single star models.')
# print(' (7) Analysis of accretion disk')
print()
print('Choose from 0 to 7, where 0 means \'all\', split multiple components by a space (q or exit to quit):')
print('Choose from 0 to 4, where 0 means \'all\', split multiple components by a space (q or exit to quit):')
runParts = input(' >>> ')
runParts = runParts.split()
if any(part in ('q', 'exit') for part in runParts):
Expand Down

0 comments on commit f3ed027

Please sign in to comment.