Posts for the month of February 2015

OutflowWind: sonic surfaces check

  1. Fixed bugs causing the subsonic speed instead of supersonic issue for the stellar wind. Here's the low res (for the frames after stellar wind kicked in) results
http://www.pas.rochester.edu/~bliu/OutflowWind/bugFix_022815/rhoT0200_sonicCheck.png movie
  1. 2.5D High Res (3 levels AMR) and longer runtime

movie; Zoomed in movie

  1. 3D Low res run

density; temperature with sonic surface

Meeting Update 02/28/15

Tried to solve the issues about small standoff radius and turbulence with 2 levels of AMR as shown in blog:bliu02232015. Found and fixed bugs in the OutflowWinds problem module. Here's the first cut results with 2 levels of AMR after the bug-fixing: parameters are same as the Stone&Proga paper. More results are coming. Will redo the case of moving ambient also…

http://www.pas.rochester.edu/~bliu/OutflowWind/bugFix_022815/rhoT0200.png movie
http://www.pas.rochester.edu/~bliu/OutflowWind/bugFix_022815/rhoT_0172.png

Meeting Update 02/24/2015 - Eddie

3-D Mach Stems

All 5 runs that I wanted to do for the 3-D Mach stem study are basically done. I still need to make a few more emission maps from different angles.

One issue that I'm running into is plotting the different Mach number runs with the same scale. The peak H-alpha emission is approximately 3 orders of magnitude stronger in the M = 10 cases compared to the M = 5 cases. See image/movie below.

movie

I think I need to plot the two models with different scales, or maybe we should consider running different models…

One idea that I had was perhaps testing something other than velocity effects. We could leave M = 10, use the two different separations (d = 3, 6), and implement an initially constant magnetic field perpendicular to the flow. This would become a 3-D study on the effects of magnetic fields on Mach stems and intersecting bow shocks.

Perhaps this is too much, since we did not do any 2-D runs with magnetic fields. We could stick to simply studying more separation distances while leaving M constant at 10.

Below are emission and column density maps with 0 deg inclination for all M = 10 models. From left to right: single, d = 6, d = 3.

movie

movie


Other Stuff

  • Last week we started discussing how to implement magnetic fields into the outflow module. Perhaps it would be good to have a separate meeting for this this week.
  • I/We need to write an abstract for the AstroNum talk.

Meeting update

Update e v. t plot

Last time, we were thinking this plot of total energy in the collision region over time looked a bit funny with the leveling off of the gravitational energy (yellow) in the collision region,

since the mass in the collision region was thought to be a strictly increasing function. Recall, was made by taking a section of the grid, a cylinder of radius 20 pc and length 20 pc, and then doing a weighted variable sum of the various quantities, see this page for reference.

Gravitational energy was calculated with the expression,

<gasphi*rho>,

which only accounts for the *gas* gravitational energy. This doesn't matter before 10 Myr, before which any sinks form. But after 10 Myr, it was curious how that curve just leveled off..

Here is a recent plot we made showing the particle mass compared to 1) the whole simulation box, and 2) the same collision region we used for the above Evt. plot. We see that the particles themselves add to about 10% of the total mass of the collision region. This then might imply that they aren't negligible, and should be added to the total gravitational energy expression. Further, it is interesting that the total mass of the collision region seens to level off (note this IS a log plot, making differences appear less than they are). This implies the mass flux in = mass flux out in this region.

Talk from CIRC symposium

If you go to slide 10, you will see side-on views of the 30 and 60 case. It was interesting to see them again like this.. It seems like the 60 case straightens out a lot more than the 30 case

https://docs.google.com/presentation/d/1JJ4nTgvgm9IvAxEc38g9LNKCDBRxzqYlPBCf-kr3nIc/pub?start=false&loop=false&delayms=3000

Curious what the field is doing here, wanted to see the streamlines at the final time for the different runs. Here is a projection of the streamlines through the domain at 27 Myr,

and just another pretty shot of the mesh and the flow,

Isotropic Turbulent Sim

Got a version of this module up and running for 2D to generate filaments for Amira to analyze,

movie

Looking at the V2 spectra, we see the characteristic turbulent decay,

Amira output -

x,y positions of segments along filaments:

Nodes, segments, points:

Comparing the filaments Amira finds of a certain threshold to a pseudocolor plot of density,

3D single star with pulsation

This post is aimed to reproduce G. H. Bowen 1988

profile is fixed.

The density at the base which is at 1AU radii is also fixed as 7.52E-11 (g/cc).

Initial velocity varies in time. Its velocity in quiet phase is almost zero, the pulsed velocity is much below the escape velocity - 29.7 km/s.

The flux varies from to

flux.mov

rho.mov

OutflowWind: moving ambient and standoff radius

1. Moving Ambient

Check the idea of setting ambient moving starting from the beginning at the same speed as the stellar wind. Stellar wind will kick in after some time t0, just as the case of ambient not moving… The following results are for ambient density 1e-4 g/cc with speed 1.638e3 km/s — the calculated stellar wind speed from the stone&proga paper. It shows the planet wind couldn't expand from the top…

http://www.pas.rochester.edu/~bliu/OutflowWind/movingAmbient/zoomedRho0073.png; http://www.pas.rochester.edu/~bliu/OutflowWind/movingAmbient/zoomedT0073.png density movie; temperature movie

2. Standoff Radius

check the standoff radius for different wind velocity. The ambient has initial velocity zero. The results show the standoff radius is around 7cu for subsonic speed wind and 6 cu for supersonic speed wind.

v=8.92e2 km/s, 50X100, 1AMR http://www.pas.rochester.edu/~bliu/OutflowWind/standOffcheck/Rho_50X100_1AMR_v8.92d2.png; http://www.pas.rochester.edu/~bliu/OutflowWind/standOffcheck/T_50X100_1AMR_v8.92d2.png density;temperature
v=1e2 km/s, 25X100, 0AMR http://www.pas.rochester.edu/~bliu/OutflowWind/standOffcheck/Rho_25X100_0AMR_v1d2.png; http://www.pas.rochester.edu/~bliu/OutflowWind/standOffcheck/T_25X100_0AMR_v1d2.png density;temperature
v=1e0 km/s, 25X100, 0AMR http://www.pas.rochester.edu/~bliu/OutflowWind/standOffcheck/Rho_25X100_0AMR_1d0.png; http://www.pas.rochester.edu/~bliu/OutflowWind/standOffcheck/T_25X100_0AMR_1d0.png density;temperature
v=1e0 km/s, 25X100, 2AMR http://www.pas.rochester.edu/~bliu/OutflowWind/standOffcheck/Rho_25X100_2AMR_v1d2.png; http://www.pas.rochester.edu/~bliu/OutflowWind/standOffcheck/T_25X100_2AMR_v1d2.png density;temperature

OutflowWind: 3D

Higher interpolation order + H_viscocity; Wind_velocity=8.19e2 km/s

density without wind

temperature without wind

density with wind

temperature with wind

OutflowWind: 2.5D Higer oder Interpolation

Linear Interpolation + H_Viscocity

With Wind (vel 8.19e2km/s)

Zoomed In density

Zoomed in Temperature --fixed colorbar

Density

Temperature

With No Wind

1. symmetric profile

density temperature

2. asymmetric profile

density temperature

Completed mass query for theta = 0 case


Figure. A mass versus time plot for the four sinks that form in the case of CollidingFlows, along with M_Simulation (the total mass in the simulation box), and M_Collision (the total mass in the collision region within the box).

  • Time interval is for ~13 Myr.
  • Mass axes is logarithmic in order to better illustrate the trend. A linear axes lets the sink masses look close to zero in comparison to the mass in the box and cylindrical collision region.

Followed a similar procedure as to what we did for the E vs. t line plots, however did a query for rho instead.

I did the query using both VNC and the interactive queue on both Erica's and my own account on BH2. For the first 200 frames of the M_simulation query, Clover collected most of the data based on the chombos stored locally. Surprisingly, Clover was faster than both VNC and the interactive queue. Ultimately both of these types of remote visualization settings are really unreliable and result in a lot of babysitting. It took me about 3 days to collect all the data for the box and cylindrical collision region in Visit. Here are some issues & interesting things I encountered:

  • Want to use a GPU? Need to use VNC which has a GUI. GUIs are super flaky and prone to time out errors it seems. Here is the command to run visit with a GPU though:
    module load virtualgl visit/2.7.3 && vglrun visit
  • Attached is the script I used for use VNC. CIRC Website also has some stuff on remote visualization (using VNC). It might be faster and nicer to use on data that isn't too demanding on memory.
  • Wanted to query using the interactive -p standard for say -t 300. Whenever I tried to do this, after Visit collected data for a few frames, it would time out. Seems like there were some memory issues. So I just stuck to an hour in the debug queue and monitored the memory/cpu percentage on the node I was on. Here is a website explaining how to do that. This implies I can only query for approximately 10 frames per interactive job.
  • Apparently using query on data that utilizes the cylindrical clip operator requires more memory than just query for the total mass in the box. Visit is probably doing extra work. Just an FYI.

Jonathan suggested making a post processing element in astrobear that'll just spit out all the data into a curve file during a batch job. Think if we want these for the three other runs I will just do that…

Moral of the story: Using query in visit for large data sets is finicky, be careful!

VNC Script

Notes:

  • Sign into your machine/local machine.
  • emacs vnc_connect_linux.sh -nw
  • Paste script in.
  • Make sure it is executable and run it.
  • It should prompt you to sign into BH2, and for how long you want your session to be, along with the size of the GUI, etc.
  • Hit enter when it says a password is found. (FYI) You'll have to make an extra password for your account too when it establishes the tunnel.
#!/bin/bash -i                                                                                                                              
via=bluehive2.circ.rochester.edu
#TurboVNCDir="/opt/TurboVNC/"                                                                                                               
#vncviewer=$TurboVNCDir/bin/vncviewer                                                                                                       
vncviewer=vncviewer

read -p "Please enter your netid: " user
read -p "Do you need to start a VNC server? [y]:" vnc_start
read -p "Set a timeout for your VNC server [60]:" vnc_timeout
read -p "Choose a resolution [1280x1024]:" vnc_resolution

if [[ -z "$vnc_timeout" ]]; then
  vnc_timeout=60
fi

if [[ -z "$vnc_resolution" ]]; then
  vnc_resolution="1280x1024"
fi

if [[ $vnc_start =~ ^[Yy]$ ]] || [[ -z "$vnc_start" ]]; then
  echo
  echo "Now connecting to bluehive and starting the "
  echo "VNC server."
  ssh $user@$via "vnc_start -t $vnc_timeout -g $vnc_resolution" # | grep "vncserver running on " |  cut -d " " -f 4                         
fi

read -p "Please enter server (ie bhx0101:1) " server
host=`echo $server | awk -F ':' '{print $1}'`
display=`echo $server | awk -F ':' '{print $2}'`
port=$(expr 5900 + $display)

echo "Establishing ssh tunnel"
TMPSOCK=`mktemp -u XXXXXX.control`
ssh -fN -o ExitOnForwardFailure=yes -M -S /tmp/$TMPSOCK -L $port:$host:$port $user@$via
echo "Launching VNC viewer"
$vncviewer localhost:$port
echo "Disconnecting ssh tunnel"
ssh -S /tmp/$TMPSOCK -O exit $user@$via

Slab symmetric sims of planet atmospheres

2D version with stepped on boundaries

2D version with outflow only wind at top (semi-extrapolating)

meeting update

Not much to report, working on the paper (will have a draft of the sections I am working on in a couple of days on the wiki), and going to put together a short talk for the CIRC symposium Friday

Meeting Update 02/16/2014 - Eddie

Fixed an error in the wind boundary in my Mach stem module, and reran one of the simulations. Images/movies below are for the d = 3 rclump, M = 10 model. The d = 6, M = 10 model is running right now and should finish sometime today.

movie

movie

previous movie

Planetary Atmospheres

Ran simulation using hydrostatic profile from r_inner to the boundary, but began to see unusual behavior. Very high temperatures caused the timestepping to slow to a crawl.

And this shows lineouts along the x axis of

  • Temperature
  • density
  • mach number
  • |vx|
  • Pressure
  • |g|
  • |grad(p)/rho|
  • |grad(p)/rho + g|
  • |v dot grad(v)|

I then reran this with the outer boundary being stepped on and saw the following:

Summary of things we've tried

  • Pressure Equilibrium - Led to temperature jump which broke radiation solver
  • Isothermal atmosphere - Self-gravity solver couldn't reliably converge
  • Point Gravity with both inner and outer boundaries - heating of lowest density material
  • Point Gravity with just inner boundary - very high temperatures and tiny time steps.
  • Also - low res simulations cause density protections regardless of time step due to 'weakness' with 6 solve CTU.

Meeting Update 02/09/2015

  • 2D OutflowWind module Fixed a bug in wind object about the velocity direction when the wind applied on the +y direction. The following is the new results with wind velocity . The sonic surface looks very different from the Stone&Proga.
Fig 6 reproduced http://www.pas.rochester.edu/~bliu/OutflowWind/2DwithWind/zoomedRhoV_lambda5_wind1.623e3kmps.png;http://www.pas.rochester.edu/~bliu/OutflowWind/2DwithWind/zoomedTV_lambda5_wind1.623e3kmps.png density;temperature
Stone&Proga http://www.pas.rochester.edu/~bliu/OutflowWind/2DwithWind/Fig6.png
  • Ablative RT Growth rate result from Ruil:

http://www.pas.rochester.edu/~bliu/AblativeRT/3Dcase/ThickTarget/astrobear-bubble_thickT.jpg

Re-running the job to longer time.

Core mass script documentation

Toward the end of last week I worked on writing a script that would open, read and write the mass data in the sinks_*.okc files to a new .csv file (In particular for the CollidingFlows problem, see our data here: CollidingFlowsFigures). The purpose of this was to gain all of the sink data over time, so that we could visualize it, as you can see by the results I have posted below. These charts will allow us to observe how the sinks accumulate mass over the time of the simulation. Here I will document the development of my code, and discuss what else I would like to do with it.

Objectives

Editing code:

  • Write information of the density/mass of the cores to a .csv file
  • Convert the information from the .csv into solar masses and Myr.
  • Take converted .csv file and make graphs using matplotlib, or excel. Excel is quick, but a code that can do all of this and generate plots in one swoop would be super efficient for later uses.
  • Document and post concise and general form under the VisIt page as a script tab.

So far I have completed the first bullet~ (02-09-2015)

Editing charts:

  • Convert mass to solar masses.
  • Convert time to Myr.
  • Crop the x-axes to when the first sink forms. Get rid of the 0 values.
  • Fix x and y-axes to be the same for all the runs.

okc_to_csv.py v. 1

Screen capture of code Screen capture of .okc file
okc to csv v. 1 screencap screen cap of okc

The code reads the first line of the .okc file, splits the numbers into a list of ints. It then uses those numbers in order to access the data below that starts at the 34th line. The only parts that are hard coded are the headers for the columns of the csv file (L17) and the number of lines of data it has to read into the .csv (L27). Essentially you change this by the number of sinks that form by the end of the simulation. You can check this by counting the number of lines of data that are spit out by the final frame of the simulation.

Results 02-09-2015

Shear 0 chart of unconverted sink masses (s0)
Shear 15 chart of unconverted sink masses (s15)
Shear 30 chart of unconverted sink masses (s30)
Shear 60 chart of unconverted sink masses (s60)

Update 02/03/2015 - Eddie

2D Mach stems

gamma = 1.4, d = 5 rclump results in wide Mach stem

movie


3D Mach stems

Now have column density plots for d = 3, M = 10 case. I also changed how the angled emission maps are produced; inclination angles of 0 deg and 90 deg do not use a camera object. It will be tricky to look at the different angles side-by-side, but at least we can now be sure that any blacked-out regions are actually in the data and not being caused by a camera object.

Inclination Density | Emission
0
movie
30
movie
60
movie
90
movie


Other Stuff

meeting update

Working on the results section of the paper to get out to folks.

Worked on some new streamline plots last week in collaboratory on the figures page.

Working on some mass/sink plots with Marissa

OutflowWind: Tests with Wind

Updated the problem module with a flag lWind to turn on/off stellar wind. The stellar wind is applied after t_wind>0 to allow the planet with asymmetric temperature profile to relax to a stable state before the wind kicks in. Also the wind is applied along -y direction for 2.5D and -x direction for…

  • Stellar Wind Speed In Stone&Proga and . Correspondingly in AstroBEAR, , . For case, , so the calculated stellar wind speed . The speed is way too high for the code with the current parameters…
  • 2.5D restuls
density;temperature
http://www.pas.rochester.edu/~bliu/OutflowWind/2DwithWind/zommedRhoV0136.png;http://www.pas.rochester.edu/~bliu/OutflowWind/2DwithWind/zommedTV0136.png density; temperature;zoomed density;zoomed temperature
http://www.pas.rochester.edu/~bliu/OutflowWind/2DwithWind/zoomedRhoV0200_windSpeed16.3kmps.png;http://www.pas.rochester.edu/~bliu/OutflowWind/2DwithWind/zoomedTV0200_windSpeed16.3kmps.png density; temperature
Stone & Proga http://www.pas.rochester.edu/~bliu/OutflowWind/2DwithWind/Fig6.png
  • Interface growth rate for Ablative RT Regenerated 50 frames txt files for Rui as 200 & 100 frames takes too long to read in…
  • New user Helped Karan set up on local machine and worked him through using the code.

on fall back disk

Circumbinary Disk

Typical circumbinary disks are mainly made of dusts and accompanied with some gas (Jura & Kahane 1999). The disk is at low temperature, usually 100 K to several hundreds K (Bujarrabal et al. 2005). Our module is purely hydro, and we can not simulate dust formation. However, if the disk is purely gas, it will disperse or be reaccreted onto the stars.

The inner radii of a disk is 2-20 AU, the outer radii of a disk is 100-2000 AU. (I am looking for specific datas)

People has not found a stable circumbinary disk in simulation yet. It is possible that certain physics is not included in all the simulations (M. Akashi & N. Soker 2007).

Tidal Force

Strong tidal force will make the boundary of the primary star invalid. If we set the upper limit ratio of tidal force to gravitational force to be , the radius of the primary is , separation is , primary mass is , secondary mass is , then:

\[\frac{2 r G m_2}{d3}∧frac{m_2 G}{d2}<\gamma\]

That means:

\[\frac{2r}{d}<\gamma\]

We set our , if let , then . This is an example of lower limit of binary separation.

Separation, Mass Ratio and Accretion Mode

Different binary separation and mass ratio have impact on shape of the Roche Lobe. If the wind emitted from the primary star is subsonic at L1 point, it will become wind-RLOF, if the wind is supersonic at L1 point, it is more like wind accretion picture, if the envelop its self is filling the Roche Lobe of the primary, it is RLOF.

In simulations, we want to fix the initial mass of the primary and the secondary and to vary the effective mass of the primary to adjust the Roche Lobe while keeping the orbital period the same. In the Temperature section we can see how the sonic radii behaves when we change the effective gravity coefficient - the greater the effective gravity, the farther the sonic radii. Therefore, to get a RLOF, we need a high , to get a wind accretion mode, we need a low .

Radiation Pressure

The effective gravity of the primary exerted on the gas is:

\[f_g=\frac{\alpha G m_1}{r2}\]

Where

Prof. Morris pointed out that we can make alpha to be a function of radius, I believe this is true because in reality the dust formation and condensation will change the opacity and momentum coupling.

Another fact is that will decrease to zero during the helium flash since the luminosity of AGB star is exceeding the Eddington Luminosity but then increase slowly during the decaying phase and reach a certain value at quiet phase, thus making time dependent, too. If we have enough computational resource to simulation the whole period of a pulsation, I think we can see the envelop pushed out and then fall back.

Feed Back from the Secondary

As more material accrete onto the secondary, more gravitational energy will be released from the star. It is a good reason for us to let the high accretion rate to create a strong photon flux and offset some secondary's gravity (similar to primary's ). This will make set an upper limit of the accretion and slow down the reaccretion onto the secondary.

Temperature

Wind-RLOF require the wind speed to be subsonic at L1 point, this set a lower limit on temperature. According to Parker's solution:

\[r_s=\frac{\alpha G m_1}{c_{s}2}=\frac{\alpha G m_1}{k T/m}\]

Where is the mean particle mass. At temperature between to atoms start to condensate into molecules and dusts. So we can take .

The temperature in the photonsphere of an AGB star is about .

On the other hand, we do not want the dispersion to be too large, we can incorporate some cooling, especially low temperature cooling.

Mass Loss Rate

Typical AGB mass loss rate is between . We can make the quiet phase mass loss rate to be and He-flash mass loss to be . The duration of He-flash is several years.

Wind Speed

Quiet phase wind speed is greater than the combined escape velocity of the binary stars. Pulsation wind speed is less than the primary's escape velocity.

Pulsation

Some reference say the typical duration of pulsation is 200 days to 2000 days (Kyung-Won Suh 2014), which is more friendly for simulation. So, I assume the pulsation duration to be 1000 days at this moment and decays one E fold in 10000 days. The density during pulsation is 500 times higher than the quiet phase, so after 50000 days, the density will drop to quiet phase density.

From my experience, we can choose following parameters:

parameterswind-RLOFwind accretion
separation 20502050
primary mass 1.51.51.51.5
secondary mass 0.50.50.50.5
effective gravity 0.120.30.060.16
sonic radii (if isothermal) 12.631.46.316.8
temperature 2000200020002000
quiet phase wind speed 20.030.013.017.0
pulsation wind speed 17.024.011.015.0
pulsation duration (days)1000100010001000
pulsation decay time (days)10000100001000010000
period (days)23100913002310091300

A typical simulation box for separation binary will be which will be divided by base cells. Therefore a base cell has physical dimension of . Using 2 levels of AMR will make the finest grid be , thus the primary star will have a radii of 8 cells. This is actually 8 times larger than my previous simulations.

If the simulation time range for this simulation is 1000 yrs. I guess it could take 200*120 CPU hrs or more.