Posts for the month of August 2015

OutflowWind: planetary wind only (large window and low-res)

in C.U. or is infinite as there's no stellar wind. always. Type III in Matsakos.

http://www.pas.rochester.edu/~bliu/OutflowWind/3D_PlanetWind_Corot/rhoV_lambeda2.5_Omega1_8zones_L600_Outwind_only.png 8-zone-per-radii movie

SED and pictures at specific wavelengths

I use Mie model to calculate the dust optical properties as follow.

Then I put a mixture of 0.1 um olivine and 0.3 um pyroxene in the hydro simulation. The spatial distribution of the dust is a superposition of the following three distribution. Let me start from the output of hydro simulation:

  1. Create a 3 AU radius sphere that has 0 dust since this region is probably too hot for dust. You can imagine we are in a spherical coordinate.
  1. Make the distribution along z-axis as a gaussian and the spread increases as radius increase. Here we are in cylindrical coordinate. The functional form is . The purpose of this layer is to create a thin disk of dust.
  1. Add a dust distribution that is proportional to the gas density. This dust density at equatorial plane is low compared to the dust density in step 2 but higher in polar or other regions.

In this simulation, the total dust mass is only , it is 20 times less than the dust mass in Kervella's paper.

You can compare the result with Kervella's result in http://www.eso.org/public/archives/releases/sciencepapers/eso1523/eso1523a.pdf at page 9.

I put a 3500 K (effective) AGB star at the center of the hydro simulation result. Its radius is 123 solar radii, equivalent to 2000 solar luminosity. This type of AGB star fits in HR diagram with ZAMS evolution. The SED is from 0.5 um to 100 um.

In Kervella's paper, they take a reddened version of AGB spectrum as its initial SED, which I did not do. If I go through that procedure, I think the SED should looks more similar.

Below are four pictures at 0.65 um, 1.24 um, 2.17 um and 4.05 um respectively, they are in the same logarithmic intensity scale.

How to handle particles, outflows, and point gravity objects

Frame Restarts

In general, when a problem module creates an object - it is easiest on a restart, to let the same objects just be recreated. Then the problem module does not need to be modified - and will still have pointers to the objects. The exception, however, is sink particles - since they can be created by self-gravity, they can wander due to gravitational interactions, accrete, etc…

Now a problem module on a restart, has no way to identify particles as their position, mass, velocity, etc… may have changed. (although we could require problem modules to name the particles). Then problem modules would need to have something like the following:

IF (lRestart) THEN
  CALL FindSinkParticle("Star 1", Particle)
ELSE
  CALL CreateParticle("Star 1", Particle)
  Particle%xloc=...
  Particle%mass=...
  ...
END IF

Then the problem module could continue to setup their outflow objects, point gravity objects, etc…

CALL CreatePointGravityObj(Particle%PointGravityObj)
CALL CreateOutflowObj(Particle%OutflowObj)
CALL UpdateParticle(Particle)

The alternative is to store their outflow objects, point gravity objects, etc… in the chombo file. Then a problem module would look like

IF (.NOT. lRestart) THEN
  CALL CreateParticle("Star 1", Particle)
  Particle%xloc...
  CALL CreatePointGravityObj(Particle%PointGravityObj)
  CALL CreateOutflowObj(Particle%OutflowObj)
  CALL UpdateParticle(Particle)
END IF

But now on a restart the user would not easily have pointers to the various outflows and pointgravity objects…

Another complication, is that if a tracer is added to an outflow object that is associated with a particle.

The solution to that is to create the tracers whether there is a restart or not…

Step Restarts

When a restart is triggered, we need a way to roll back any particles to their previous state. It is also possible that particles would have been created or destroyed during the last step, so those would need to be recreated (or deleted).

A particle has additional sub pointers in the derived types that need to be allocated and copied

  • Particle
    • Outflow Object (part of object list
      • Shape
      • Profiles
        • data
        • fields
    • Point Gravity Object (linked list)

For now, let's assume that there is no merging of particles as this simplifies everything.

Before each step:

  • Make copies of particles

On a restart:

  • Find backed up particles and copy contents
  • Remove any new particles and their objets
  • Update particle objects with copied values for position, mass, etc…

When we do implement merging, we could move particles (and their associated objects) out of their respective lists. Then on a restart, we could put the particle back without creating new pointers etc…

3D Clumpy Paper Figures

I had some issues regenerating the projections at higher resolutions, so I just used the old projections. I think this is fine since you can't see the high resolution in a paper anyways. I did get it to work for the one zoom-in image which is arguably the only one that this is important for. Below are the 11 figures going in the results section; the other 4 figures which I did not post here are in the theory and methods sections.

OutflowWind: Co-rot planetary wind

  • Larger box For the setup in blog:bliu08182015, the radius of star is big (~78 CU). Large box will get the star in as shown in the following picture.
http://www.pas.rochester.edu/~bliu/OutflowWind/3D_PlanetWind_Corot/PW_corot_largerBox.png

Not sure if this is OK or should use different parameters since the left and top have persistInBoundaries. Here's the result ( 32 zones per radii and 8-zones per radii for the first 80 frames) of larger box but not including the star in:

http://www.pas.rochester.edu/~bliu/OutflowWind/3D_PlanetWind_Corot/rhoV_stellarRho1E-11_lambeda2.5_Omega1_32zones_L200.png L=200 movie

  • Restart Issues with High CFL
  1. This happens when restarting from a multi-core chombo files. chombo files from Single-core work OK.. I suspect there's a bug in the new parallel hdf writing code.
  2. Currently uses an older version of the code and latest module problem which works fine..

Dust Optical Properties

Using Mie Model.

pnStudy: 3D Results from IAC 08192015

3D Data with Nlevel=4 and DMcooling of up to 493 y

  1. Ran for 2.5 hrs on 160 on TeideHPC
  2. Problem.data

tamb = 1d3           ! ambient temp, 1cu = 0.1K (100K=1000cu)
namb = 4e3           ! ambient density in cell above lunch surface. 1cu = 1cm^-3
stratified = t       ! 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
StraTorus  = 2       ! 1 for Martin's way used before Jan 6th 2015, tracer might not work correctly!!
                     ! 2 for Baowei's way, updated from Martin's code according to Bruce's request
!
!     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  = 1d0         ! flow radius at launch zone, 1cu = 500AU (or clump radius)
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)  (outflowType=2 cones only)
tf    = 30d0        ! conical flow Gaussian taper (deg) for njet and vjet; 0= disable
sigma = 0d0         ! !toroidal.magnetic.energy / kinetic.energy, example 0.6


3D from IAC https://www.pas.rochester.edu/~bliu/pnStudy/IAC_Data/08192015/rho_iac08092015_493y.png 3D movie
same run in 2D https://www.pas.rochester.edu/~bliu/pnStudy/IAC_Data/08192015/rho_iac08092015_493y_2Dcompare.png 2D movie

fixed resolution issue for figures

I fixed the resolution issue. Turns out that when you use a camera object, the code automatically sets the resolution to be equal to the maximum number of cells in the x direction. Thus, if your y or z directions have a larger extent, your data will be coarsened.

Some of my 3D clump simulations have a y domain that is almost 3 times the size of the x domain. This is why I said that the resolution should be 80 cells per clump radius, but the images show approximately 30 cells per clump radius.

We should probably change this in the code so that the camera resolution matches the actual maximum resolution of your data. The old image is on the left, new one on the right.

OutflowWind: Corot PW results with or without Outflow-only for xhigh and ylow

Planetary wind only in Co-rotating frame results of setting/not setting the xhigh and ylow boundaries as Outflow-only

Outflow_only for xhigh and ylow http://www.pas.rochester.edu/~bliu/OutflowWind/3D_PlanetWind_Corot/rhoV_Outflow_only.png 2D slice movie;Sonic Surface movie;Vol movie
No Outflow_only for xhigh and ylow http://www.pas.rochester.edu/~bliu/OutflowWind/3D_PlanetWind_Corot/rhoV_Extrapolate.png movie

Update 08/18/2015 - Eddie

  • applied for guest account on Rice machines
    • we will be using DAVinCI: 2304 cores on 192 nodes (12 cores per node)
    • scratch directory for running jobs, work directory quota is 2TB per group, can use Globus for file transfer
    • can use 48 nodes (576 cores) at a time, but only for 8 hours at a time :(
    • there might also be some fees, but I'm not sure yet. From what I found online, it might be a flat annual fee of 100.00 + 340.00 per million SUs
  • made some figures for the 3-D clump paper
    • Adam, I'm taking all of your suggestions from your latest email, and adding one more figure. In addition to "orientation comparison", I want a "viewing angle" comparison figure.
    • of the figures below, the enlarged one requires the most discussion…

i = reddened, post-shock [S II], cooling regions

ii = bow shock intersection point, possibly Mach stem

iii = post-intersection point "cone"

iv = unstable, fragmenting bow shock

v = the "froth"

Planet Wind Orbital Trajectories

I calculated trajectories for planetary wind particles leaving at the escape speed from the planet surface in a corotating frame.

Here are the trajectories for a separation of 400 AU

and here are the trajectories for a separation of 200 AU

… and 100 AU

We have

and

We can calculate an coriolis radius where

which gives

spectrum

I am still using the default silicate data in RADMC3D.

I find a program that can calculate dust absorption and scattering coefficients for irregular shape dust (ofc regular works) from n-k data.

http://www.ddscat.org/

Gravitational Bounded or Not Bounded

q Measures the ratio of kinetic energy to gravitational potential. In stellar evolution theory, I should take enthalpy instead of just kinetic energy but since I use isothermal temperature, the enthalpy is not conserved so I did not take it.

By Virial Theorem, if h=enthalpy

then the material is stable on its orbit, if it is larger than, it will go to higher orbit, if it is less than, it will go to lower orbit.

is the condition that material will go to infinity (unbounded)

So in my case, q=1 absolutely means unbounded, but gas is not bounded even q is close to 1 because I did not take temperature into account.

The figure shows that the gas in polar direction is not gravitational bounded while it is in equatorial plane. (q<0.5)

At this moment, the circumbinary disk has not reached steady state yet.

Setting Initial Conditions for Planet Wind Calculations

Here are some relevant figures and equations.

First the different regimes

Here is the radius of the wind standoff.

Here is the Hill radius .

Here are parameters for the stellar wind.

Update 08/13/2015 - Eddie

  • waiting to hear back from Andy about getting on Rice machines
  • need to start generating figures for the results section of the clump/bow shock paper
  • need to start some higher resolution 3D pulsed jet runs
  • worked on some 1-D radiative shock models for Pat since he supplied us with a new cooling table
    • The new z cooling is very very close to the old z cooling in terms of total energy losses. It appears to be slightly stronger or slightly weaker depending on the exact parameters.
    • DM cooling is clearly much weaker in the range of the z cooling table, but stronger at high temperatures.

  • I did a convergence study to see how resolution affects the initial post-shock temperature of one of my models. The model appears to begin converging for resolutions > 30 cells/AU.

On the research

The first image shows a streamline plot of constant wind. There is no polar outflow, however, there is a circumbinary disk. Fluid parcels start in polar region just fall onto the disk. The disk looks stable as the red streamline means the fluid parcel position in later time.

The second image shows a streamline plot of dense shell wind. The disk is not stable as there is no limit circle around it. There is no red lines in the streamline indicating that the fluid parcel just run away.

Below shows a polytropic wind with gamma=1.35 . The wind start at 0.5 AU, 1500 K with 20km/s. It drops to about 300 K at 10 AU. This simulation is done with my 1D spherical Riemann Solver.

Heat can be transported by conduction, convection and radiation. High order Riemann solver can handle convection very well. Conduction is not important in the exterior region of star but important in interior (very deep, can be estimated). Radiation is the most important heat transfer mechanism in low density region, for Local Thermal Equilibrium (LTE) situation, Monte Carlo (MC) method is the best both for multi-dimension and 1D. In non-LTE situation, people sometimes use particular cooling function, for example, hydrogen cooling or any other molecular cooling, there are also MC model for non-LTE but I need to research more.

In general, MC method is universal in simulating radiative transfer.

OutflowWind: Planetary Wind only in Co-rot frame

  1. Latest high-res results:
Lambda=5 Omga=0.5 movie 32 zones per radii
Lambda=5 Omga=1 movie 32 zones per radii
Lambda=2.5 Omga=1 movie 32 zones per radii

details can be found in the last part of this page

  1. A bug found when doing restart.. Haven't found how to fix it…

Location of wire simulations

With the coming of my starting graduate school, I moved all of my simulations off of bluehive2 and bluestreak to /grassdata/. You can now find them at:

/grassdata/madams/TurbulentWires

where I changed the group setting to orda. Hopefully any one in the group can have access to them.

Under ProductionRuns/, I have things split between the two machines. Under bh2/ you can find the level = 1 hydro simulation, along with CDMs. Then Under bs/ you'll find the continuation of level = 2 hydro simulation, along with the mhd simulations and some code Jonathan was trying to test solvers with.

Clumpy Simulation Figure Set

The big question now is what figures to make/show. Below are my suggestions.

1) We need an initial figure set for a single 2 clump run which illustrates the basic kinds of phenomena we are seeing. I suggest this should include

a) A column density "movie" figure showing the simulation at 3 (or 4) times. b) An [SII]+Ha intensity "movie" figure with same frame times as a) c) A single (blow up) frame of the [SII]+Ha with arrows/labels for features we want to explore in text

2) An initial figure set for a single 3 clump run which builds off 1)

a) A column density "movie" figure showing the simulation at 3 (or 4) times. b) An [SII]+Ha intensity "movie" figure with same frame times as a)

3) Seperation Comparison

a) no need for "movies" here just the best example(s) illustrating how seperation effects things. Only need [SII]+Ha.

4) Velocty Comparision

a) (same as above) Perhaps we make sure we are using different simulations for each of these comparisons so that in the end we are not repeating showing the same sim.

5) Orientation comparison

a) (same as above)

6) 10 clump run.

a) A single column density frame b) A movie of 3 (or 4) [SII]+Ha intensity images at a single orientation/inclination

Simulations of the circumnuclear disk with a larger inner cavity, Part II -- Marvin

In this post I am presenting 10 simulations of the circumnuclear disk, 5 with an inner cavity of r=1pc (see first figure and first animation below) and 5 with an inner cavity of r=2pc (see second figure and second animation below). The time is given in units of yrs. The r=1pc simulations cover a time of years, the r=2pc simulations years. The following table shows the details of these simulations. (Note: B=1mG corresponds to a plasma-beta of about 0.1, B=0.1mG to a plasma-beta of about 10.)

Number properties linestyle
1. without magn. fields and with outflow blue, dotted
2. with magn. fields (B=1mG) and with outflow red, dashed
3. without magnetic fields and without outflow turquoise, solid
4. like 2., but with B=0.1mG green, dashdotted
5. with magn. fields and without outflow violet, solid

First figure: CND with an inner cavity of r=1pc

Second figure: CND with an inner cavity of r=2pc

First animation: CND with an inner cavity of r=1pc

Second animation: CND with an inner cavity of r=2pc

In all these simulations the inner rim moves inwards in one way or another. Different physical processes can be responsible for that, but they all work on different timescales: Firstly, the inner rim moves inwards due to numerical viscosity on a viscous timescale (on how to determine the numerical viscosity I refer to one of my other posts: https://astrobear.pas.rochester.edu/trac/blog/mblank05302015). Secondly, gas pressure can cause the inner rim to move inwards, the corresponding timescale is with the speed of sound . Then, magnetic pressure can also cause the inner rim to migrate inwards on a timescale , where is the Alfven speed. Finally, the inner rim moves inwards due to the angular momentum extraction caused by the interaction with the wind, the corresponding formula can be found in our paper draft. All these timescales are listed in the following table in units of yrs. is the timescale for the angular momentum extraction according to the formula presented in our paper draft, and is the time the inner rim actually needs to move inwards (as seen in the simulations).

…….

timescale r=1pc r=2pc
4.5 12.7
2400 21700
37.8 113.5
10 30
100 300
49 490
8 42
20 65

Let's first look at the simulations without outflow (number 3 and 5): There is no big difference between these two models. Thus in the absence of the central outflow magnetic fields do not have a large effect on the evolution of the CND, maybe the reason for this is that initially there is no magnetic field present inside the inner cavity. Only numerical viscosity causes an inwards migration of matter, but as the corresponding timescale is much larger than the simulation time not much matter is actually moving inwards.

Examining the simulations without magnetic fields (model 1) and with a low initial magnetic field strength (model 4) shows that initially they do not deviate much from each other (maybe the magnetic field strength that I've chosen for model 4 is a little bit too low). In both simulations the inner rim takes about yrs to reach the outflow object ( yrs for the r=2pc simulation). The first yrs ( yrs for r=2pc) they look more or less the same, after this time the surface density is larger when no magnetic fields are present.

In the simulation with magnetic fields and with outflow (model 2) the inner rim moves inwards very quickly (faster than in the simulations without magnetic fields and with a low magnetic field strength), it reaches the outflow object after about yrs ( yrs for the r=2 simulation). After that the surface density decreases about a factor of two, thus matter seems to be removed from the inner cavity. Then the simulation reaches a "steady state", where the surface density experiences no significant changes.

Conclusions:

  1. The outflow is causing the inner rim to movs inwards, the higher the magnetic field strength, the faster it is moving.

One may be tempted to think that this is due to the magnetic pressure that drives the material inside, as for the simulations with magnetic fields and with outflow (model 2) the timescales correspond to the times the inner rim needs to move inwards. But even without magnetic fields the inner rim moves inwards, one again could think that this is caused by the gas pressure, because roughly corresponds to the times the inner rim needs to move inwards. However, the simulations without outflow do not show such an inwards movement of the inner rim whatsoever, thus it is the outflow that causes the inner rim to shrink (at least it is the main contributor).

  1. After the inner rim has reached the outflow object, the surface density is lower for higher magnetic field strengths.

But as always, answering some questions leads to a larger number of new questions:

  1. Why is the collapse faster for higher magnetic field strengths? Magnetic fields seem to accelerate the angular momentum extraction.
  1. Why is the surface density lower for higher magnetic field strengths? Do magnetic fields prevent the inflow of matter or is there a kind of magnetically driven outflow?
  1. and fit quite well for the r=1pc simulations, but why is there such a huge difference for the r=2pc simulations?
  1. Do we have to answer all these questions in our paper? Maybe the first two can be shifted to our future work, but I think the third needs some additional pondering.

Spreading ring calculations, Part II -- Marvin

In one of my last blog posts (https://astrobear.pas.rochester.edu/trac/blog/mblank05302015) I presented some "spreading ring calculations" that allow to estimate the magnitude of numerical viscosity in disk simulations with astroBEAR. There I introduced the parameter that gives the numerical viscosity in units of the maximum -viscosity .

In the following I show a parameter study of these calculations by varying the resolution of the simulations. However, I have performed this study in 2D to save computational time. The simulations from my previous post have a resolution (cell size) of 0.04 pc, here I additionally show simulations with 0.08, 0.02, and 0.01 pc.

The evolution of the spreading ring for all resolution levels is shown in the Figure below (upper left: 0.08 pc, upper right: 0.04 pc, lower left: 0.02 pc, lower right: 0.01 pc).

The qualiative evolution is quite similar for all the simulations, but as expected the spreading of the ring is slower with higher resolution. I also estimated the parameter following the procedure described in my previous post, and I additionally give the runtime of the simulations, all have been calculated using 64 cores:

resolution runtime
0.08 pc 11.3 min
0.04 pc 1 h
0.02 pc 5.3 h
0.01 pc 31.1 h

The upper right figure with a resolution of 0.04 pc can be compared with the corresponding 3D simulation of my previous post, again the qualitative evolution is similar, but the 2D simulation has a higher numerical viscosity than the 3D simulation, which has .

The numerical viscosity decreases with increasing resolution, as expected. Furthermore, it seems to be a linear function of the resolution, which makes sense: Let us consider the following simplified transport equation: For solving this equation numerically we want to replace the spatial derivative by a difference quotient. Therefore we make a taylor expansion of the function at the point :

Rearranging gives:

And inserting this into the first equation gives:

The term on the right side is responsible for numerical viscosity, and is a linear function of the resolution .

OutflowWind: Co-rotating Planetary Wind & Parameters

1. Testing results about Co-rotating frame

http://www.pas.rochester.edu/~bliu/OutflowWind/3DrotDebug/rhoV_stellarRho2E-5_Omega0.5_32zones_2.png movie

Details in this page

2. Planetary wind only in Co-rotating frame

Only low-res results so far.. Got a lot of High CFL restarts for high-res runs — either by restarting from a low-res frame or starting from beginning directly..

low-res Movie for Lambda=5.3 Omega=0.5 ; low-res Movie for Lambda=10.6 Omega=0.5

Details in this page

3. Compare Parameters with Matsakos

check this page — not quite finished yet..

SED analysis and Hydro binary simulation diagnosis

L2pup is the a binary system that is located only 64 parsec away from earth. It has an AGB star and secondary that maybe a RG. Assume a 3000 K AGB star (blackbody) with 1140 solar luminosity (SL), we can calculate its radius to be 126 solar radius (SR) by standard procedures. The Spectrum Energy Distribution (SED) at 64 parsec away is shown in the picture below. It peaks at 1 um with intensity.

Red line is result from Radmc3d with the same condition of pure blackbody radiation in vacuum. Blue line is SED of blackbody radiation. Simulation is "Redder".

Compare to the results from other simulations by eye, the difference in magnitude is not very significant. Thus much of the light reaches earth.

Where can I get the data in Kervella's 2015 paper?

I find the accreting secondary has accreted too much gas material. In other papers talking about L2pup, they usually guestimate the mass of the dust is on order of 1E-7 solar mass with 1% dust gas ratio. In my simulation, the total mass of ejected gas is 1E-5 solar mass but the secondary accreted 50% of the gas as time goes on and another about 40% escape so the total mass in the disk is only about 1E-6 solar mass which make the simulation not realistic.