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.
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:
- 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.
- 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.
- 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)
- Outflow Object (part of object 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
- 3D view
- volume rendering for results in blog:bliu08182015
- 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.
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:
L=200 movie |
- Restart Issues with High CFL
- 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.
- Currently uses an older version of the code and latest module problem which works fine..
pnStudy: 3D Results from IAC 08192015
3D Data with Nlevel=4 and DMcooling of up to 493 y
- Ran for 2.5 hrs on 160 on TeideHPC
- 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 | 3D movie | |
same run in 2D | 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 | 2D slice movie;Sonic Surface movie;Vol movie | |
No Outflow_only for xhigh and ylow | 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.
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.
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
- 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
- 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 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).
(on how to determine the numerical viscosity I refer to one of my other posts:…….
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:
- 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).- 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:
- Why is the collapse faster for higher magnetic field strengths? Magnetic fields seem to accelerate the angular momentum extraction.
- 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?
- and fit quite well for the r=1pc simulations, but why is there such a huge difference for the r=2pc simulations?
- 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
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.