Posts for the month of October 2015

Decomposing gradients into wave strengths

Consider the 2D hydro linearization of the wave equation in each dimension

where

The matrix has eigenvalues , , , corresponding to left propagating sound wave, contact, shear discontinuity, and right moving sound wave. And we can construct the matrix of left eigenvalues (where each row is a left eigenvector)

and the matrix of right eigenvectors (where each column is a right eigenvector)

and the diagonal matrix of eigenvalues where each diagonal term is an eigenvalue

The eigenvectors are orthogonal

and we also have

which combining gives

Now since the eigenvectors are orthonormal, we can multiply any left eigenvector by an arbitrary normalization constant, provided we divide the corresponding right eigenvector by the same amount.

Rewriting the equation

Now if we renormalize the left eigenvectors

and the right eigenvectors

then

where

has units of

This gives us the relative strength of each wave at each cell. This looks at gradients in the fluid variables and how they correlate with the various types of waves.

For example, if we look at the Sod Shock Tube test problem whose solution entails a left rarefaction, contact, and right moving shock, we see that the 3 different waves are correlated with the amplitude of the projections of the gradients onto the left eigenvectors associated with each wave.

And if we look at the Brio Wu shock tube which consists of a fast rarefaction, a slow compound wave, a contact, a slow shock and a fast rarefaction, we see the following wave strengths

Meeting update

Restart of HD2794b with additional refinement

MHD run consistently stalls while outputting frame 109. Output file appears to be 3x larger than previous frames. Must be a bug in the io routines.

meeting update

Heavily updated the M2FR and the protocluster formation sections of the paper. Also cleaned up the intro. Am now working the rest of the way through the paper.

HD209458b: PlanetaryWind and Sonic Surface

When trying to produce a high-res picture of HD planetary wind with larger planet, I found the different-looking of sonic surface for L=60 and L40 as shown in blog:bliu10192015 .

I tended to think the L=60 is more correct(?) as

  • The SF of L=60 is more like the 2D PW simulations I did before (comparing with 2D sonic surface ).
  • In early frames of L=40 has similar sonic surface as shown in this picture

http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/SonicSurface/L40_4AMR_frame115.png

Thought this is resolution related. I did several restarts with higher AMR levels. Here's some of the results

L=40 4-AMR outside the planet http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/SonicSurface/L40_4AMR_amb.png
L=40 5-AMR outside the planet http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/SonicSurface/L40_5AMR_amb.png

movie

Will try higher-res for L=60 to see if the results could be consistent..

Higher-res results for L=60 shows similar sonic surface as those of L=40. So we conclude the sonic surface of higher resolution is more correct although it looks slightly different from our 2.5D results

L=40 4-AMR outside the planethttp://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/SonicSurface/sonicSurface_L60_4AMR.png movie

Update on HD sim

movie

And WireTurbulence run is at frame 106/200

meeting update 10/19/15

High-res PW Picture for HD209458b:

L=60 http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/MassLoss/HD209_L60_high.png movie 1.56E+09 gram per sec
L=40 http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/MassLoss/HD209_L40_high.png movie 1.42E+09 gram per sec

Update 10/19/2015 - Eddie

Below is image/movie showing the hydro 3-D pulsed jet: slice of x-z plane, slice of y-z plane, and column density.

movie

It might be helpful to review how the outflow object is set up: OutflowObjects

My jet uses a radius = Rjet/2, thickness = Rjet/4, base = Rjet, and open_angle = 0. The outflow object steps on some cells, but the rest of the lower boundary is extrapolated; should probably change this to reflecting.


Things To Try

  • use different Riemann solver
  • adjust CFL
  • adjust boundary conditions of outflow object (radius, thickness, fade, etc.)
  • turn off pulsations
  • have precession fade from 0 to Rjet


UPDATE

Things i have tried thus far:

  • For some reason, the exact Riemann solver does not work well. Code stops early due to nan in flux.
  • Added velocity fade from 0 to Rjet. This is not a precession fade, just a velocity magnitude fade. Still see waves.
  • Turned off characteristic limiters, use primitive limiters instead. Still see waves.
  • Turned on ApplyDiffusion.
  • Lowered target CFL from 0.3 to 0.1.
  • Turned off pulsations.

M2-9: 3D pnStudy with rotating/spinning conical wind

  • Test the 3D pnStudy module by adding a rotating velocity to the conical wind:
Non-rotate http://www.pas.rochester.edu/~bliu/pnStudy/M2-9/rhoV_M29_3D_noRot.png 4AMR non Rotate movie; 4AMR non Rotate 2D slice movie;
Rotate with 150y period http://www.pas.rochester.edu/~bliu/pnStudy/M2-9/M29_3D_rho_rot_f44.png 4AMR Rotate movie; 4AMR Rotate 2D slice movie;
  • The data file for this run:
tamb = 1d3           ! ambient temp, 1cu = 0.1K (100K=1000cu)     
namb = 4e4           ! ambient central density cm^-3. Usually 400 for 1/r^2 or torus.
stratified = f       ! true = add a 1/r^2 background 'AGB stellar wind'
torus      = f       ! true - add torus to the background
torusalpha = 0.7     ! alpha and beta specify the geometry
torusbeta  = 10d0    ! see Frank & Mellema, 1994ApJ...430..800F
rings      = f       ! true - add radial density modulations to AGB wind
!
!     FLOW  DESCRIPTION SECTION, values apply at origin at t=0
outflowType  = 2    ! TYPE OF FLOW    1 cyl jet, 2 conical wind, 3 is clump 
njet  = 4d2         ! flow density at launch zone, 1cu = 1cm^-3
Rjet  = 2d0         ! flow radius at launch zone, 1cu = 500AU (outflowType=1 only)
vjet  = 2e7         ! flow velocity , 1cu = cm/s (100km/s=1e7cu)
tjet  = 1d3         ! flow temp, 1cu = 0.1K (100K=1000cu)
tt    = 0.0d0       ! flow accel time, 1cu = 8250y (0.02 = 165y)
open_angle = 90d0    ! conical flow open angle (deg)
tf    = 15d0        ! conical flow Gaussian taper (deg) for njet and vjet; 0= disable
sigma = 0d0         ! !toroidal.magnetic.energy / kinetic.energy, example 0.6

HD29:planetary wind and Mass loss

1. Planetary Wind With Temperature Profile

I have been trying to get a high-res picture for planetary wind HD29 with temperature profile. This picture shows there's a shock from the night-side and there's a back-flow.

frame 40 http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/MassLoss/RhoV0040.png
frame 48 http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/MassLoss/RhoV0048.png

The shock disappeared as time goes further, although still need to double-check if this comes form the restarting problem

frame 52 (2AMR) http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/MassLoss/RhoV0052.png movie for this run
frame 76 (3AMR) http://www.pas.rochester.edu/~bliu/OutflowWind/HD209458b/MassLoss/RhoV0052.png movie for this run

2. Estimate of the Mass Loss Rate

  • 1). The mass loss rate is estimated by calculating the flux out of the cubic surfaces as show below:
http://www.pas.rochester.edu/~bliu//OutflowWind/HD209458b/MassLoss/integral.png

Where the size of box is 0.4*Orbital separation between the star and the planet, or 0.019 AU and the planet is sitting in the center of the box.

  • 2). The mass loss rate for no-temperature-profile case is 3.82E+09 gram per sec, see blog:bliu09222015
  • 3). For with-temperature-profile case, I checked the mass loss rate for different frames
frame Picture Mass loss rate
Frame 48 http://www.pas.rochester.edu/~bliu//OutflowWind/HD209458b/MassLoss/RhoV0048.png 1.33E+09 gram per sec
Frame 70 3AMR http://www.pas.rochester.edu/~bliu//OutflowWind/HD209458b/MassLoss/rhoV0070.png 1.55E+09 gram per sec
Frame 76 3AMR http://www.pas.rochester.edu/~bliu//OutflowWind/HD209458b/MassLoss/rhoV0076.png 1.56E+09 gram per sec

Mach stem exploration

Running Eddie's setup,

time = 39.7 tau; M = 30.02; d = distance between clumps centers;

Summary: Simulation Extrapolation/Interpolation

d [R_clump]\gamma 5/3 1.4 1.2
3 TMS, SB TMS, SB TMS, SB
3.25 TMS, SB TMS, SB SMS
3.5 TMS, SB TMS, SB RR
4 TMS, SB TMS, SB RR
4.5 TMS, SB TMS, SB RR
5 TMS, SB SMS RR
5.5 TMS, SB RR RR
6 TMS, SB RR RR
6.5 TMS, SB RR RR
7 SMS RR RR
7.5 SMS RR RR
8 SMS RR RR
8.5 SMS RR RR
9 RR RR RR

Snapshots:

gamma = 1.4; d = 6, 5, 4 clump radii:

Update 10/13/2015 - Eddie

2D Mach stem results with low M

low M = 1.55

TMS, SB = Transient Mach Stem, develops into a Single Bow shock

SMS = Stable Mach Stem

RR = Regular Reflection

gamma d Result
5/3 24 RR
5/3 23 SMS
5/3 22 SMS
5/3 21 SMS
5/3 18 SMS
5/3 12 SMS
5/3 11 SMS
5/3 10 SMS
5/3 9 TMS, SB
1.4 21 RR
1.4 20 SMS
1.4 19 SMS
1.4 18 SMS
1.4 17 SMS
1.4 15 SMS
1.4 12 SMS
1.4 11 SMS
1.4 10 SMS
1.4 9 SMS
1.4 8 TMS, SB
1.4 6.5 TMS, SB
1.4 6 TMS, SB
1.4 5 TMS, SB
1.4 4.5 TMS, SB
1,2 18 RR
1,2 17 SMS
1,2 16 SMS
1,2 15 SMS
1,2 13 SMS
1.2 12 SMS
1.2 9 SMS
1.2 8 SMS
1.2 7 TMS, SB
1.2 6 TMS, SB
1.2 5 TMS, SB
1.2 4 TMS, SB


3D Pulsed Jets

Ran the MHD, beta = 1 model farther at 32 cells/Rjet.

movie

Started a production-quality hydro run at 64 cells/Rjet.

movie


3D Clump/Bow Interaction Paper

  • Having some trouble with ds9 and the .fits images. I just want images that look like the ones from Pat's 2011 paper, but I can't seem to get the scaling quite right.

Meeting update

  • Continue to work on implementing ThermalConduction solvers.
  • Plugging away at MHD Wire Turbulence sims using ½ of BGQ. Getting a 26GB frame every 8 hours…
  • Ran HD… sim out longer - still hasn't reached a steady state…

On the research

This model is still the dense shell + constant stellar wind model.

AGB mass: 2 solar mass with alpha=0.134

Secondary mass: 0.5 solar mass

Mass loss rate in the burst is: 9.3E-6 solar mass / year last for 11.6 years

Mass loss rate of stellar wind is: 1.65E-7 solar mass / year

The speed of burst is between the escape velocity of the AGB star and the binary system. The speed of stellar wind is above the escape velocity of the binary system.

Simulation is in isothermal state at 400 K.

0.65 micron and 0.55 micron respectively. They looks identical…

SED looks more reasonable. The red line is the SED I got in current model. Note the flat spectrum energy signature from 1 micron to 4 micron. This meets the observation in Kervella et al (2015). Also, the 10 micron bump is there.

Dust composition is 0.8 % of olivine and 0.2 % of pyroxene.

topview

sideview

Anisotropic Thermal Diffusion

Working on implementing Anisotropic Thermal Diffusion

The 4 columns correspond to Isotropic diffusion, Anisotropic but with , parallel diffusion, and perpendicular diffusion. And the 3 rows correspond to Backward Euler, Crank Nicholson, and Forward Euler (explicit).

The magnetic field is uniform in the +x+y direction and In each case the effective diffusion constant is 1, , each box is 2 units across, and the total time = 1. Hydro is turned off and these are run using a single core in fixed grid.

The pseudocolor is showing Temperature and the initial circle has a temperature ½ that of the surrounding.

movie

Newton iteration for solving implicit non-linear PDE's

HD...

Figures

Meeting update

2D Mach stem results for low M, gamma = 1.4

M = 1.55, gamma = 1.4

d = 8, 9, 10, 11, 12

All appear to show a Mach stem at the end of the simulation (triple points still visible) except for d = 8 which forms a single bow.

Now increasing d to find regular reflection