Column Density Maps (CDM) of Wire Simulations
This is the first completed run. Hydro case, where I switched from MaxLevel = 2 to 1 at frame 62.
Here is a basic model I made, to make the axes more clear:
down x
down z
The integrated mass for the CDM down-x is predominantly between 37 and 1000. However looking at x_0=1.5 to x=3 in the CDM down-z, I suppose it makes sense how we see down-x being very high in mass.
For the CDM down-x, I replicated the BOV down each axis three times. For down-z, I replicated it 3 times in the Y-direction to be consistent with a 3x3 wire mesh.
Things to do
- Read Winter 2002 A&A
- Make SED of the radiative transfer simulation.
- Nail down a proper luminosity for the AGB star. Mass loss rate is not clear now, we can treat it as a free parameter but just need to be reasonable.
- Try more species of dust.
Brainstorming for October shots - multiple wires?
The upcoming Astroshock campaign shot day on October 27 will feature driver tests; we also anticipate a few magnetized wire shots (possibly done with new drivers being tested on the same day).
Last year we tried to reproduce magnetospheric structures in an experiment at LLE by having a blast wave overrun a single wire of ~600 um thickness (conductor ~400 um thickness). Due to various real-world complications the magnetic influence on the flow was weak, at best.
With the upcoming shots, we would like to modify both the upstream (i.e. driver) and the downstream (wires) to make the magnetic effects much more obvious, one proposal involves using multiple wires.
Running with this, I was inspired by Eddie's work on Mach stems: Mach stems can be an obvious feature due to brightness, and they may be robust due to hysteresis:
Mach stem evolution can be controlled by adjusting the separation of the obstacles/clumps/wires.
The separation of the obstacles can be controlled by the state of the magnetic field.
Meeting Update
Reviewed Matsakos initial conditions
- Initial magnetic field is a dipole outside that transitions to a uniform Z field inside of a sphere of radius R/2. No details of smoothing distance in paper.
- Initial density field is uniform in side a sphere - and is parker wind solution outside
- Velocity is zero inside sphere and parker wind solution outside
- Pressure is isothermal outside (from parker wind), and hydrostatic inside.
- Note they calculate P from radial momentum equation, isothermal EOS, and solution for parker wind velocity. I calculate radial velocity and use density equation to get rho, and isothermal EOS to get pressure.
- Density is constant inside and parker solution outside.
- Gravity is GM/R2 outside, and uniform sphere gravity inside. Note that M used for outside gravity does not equal 4/3 pi r3 rho inside
- Density, velocity, and pressure are stepped on inside of r=1.5 R (however, we cannot do this without analytic solution for non-isotropic parker wind).
- Magnetic field is free to evolve inside and outside.
- They discuss 4 layer model of Matt and Pudritz
- v_poloidal parallel to B_poloidal
- rho and P held constant
- v set to zero ( v_phi is corotating)
- B_poloidal held constant and B_phi is extrapolated inward to not have any poloidal current (curl of B_phi = 0)
- Each of these layers is approximately only 1 zone thick (the surface of the star is 30 zones in radius)
Parallel HDF5 writes
- The development branch of the code now supports parallel HDF5.
- On bluehive need to use astrobear/b3 module
- Need to add -D_PHDF5 to CPPFLAGS in Makefile.inc
- Set lParallelHDF5=.true. in GlobalData namelist.
- Haven't tested this on BGQ.
- If restarting from older files, need to open and rewrite chombo files without lParallelHDF5=.true.
- There is a lReOutput flag that can be set along with lPostProcess to output the HDF5 files in a format compatible with ParallelHDF5=.true. restarts.
New routines for setting up initial/boundary conditions
- Developed these for potential use by Farrukh on Magnetized Kelvin Helmholtz
Paper draft in pdf form
Working on draft. Attached.
Update 07/28/2015 - Eddie
- new movies: MachStems
- best ones to look at are 2-clump comparisons
- the column density movies are useful for watching the transmitted shocks move through the clumps, and the movies also support the idea that the secondary "messy" emission is due to material shedding off of the clumps
- working on the paper; mainly the theory section to include discussion of various instabilities in these sims.
- Adam, I'll send you an up-to-date version of everything I've written so far on Friday?
- I could continue making emission map and column density movies, but I don't know how much more information we'll get. I have generated 55 movies already, out of a potential total of 92. I'll do some more, but I might not do all 92. I think my focus should shift to reading and writing now. I could also start some new 3-D pulsed jet runs.
- I'm also starting to pull out some specific details from Pat's 2011 third epoch paper. One thing that stuck out was this image:
In the paper, Pat writes "One can also imagine the structure in the froth as arising from the irregularities in the shape of the wing of the bow shock as it sweeps up material". The bow shocks in our simulations show irregularities due to the strong cooling, and I believe this is the Vishniac instability (a thin layer bounded by shock and contact).
PPN
I used radmc3d to produce a radiative transfer simulation based on the result of my hydro simulation. The radiative transfer code use Monte Carlo method. Below are the wavelength picture correspondingly from observation and my simulation.
This is a synthetic radiative transfer simulation. It corresponds to the dense shell modeldense shell hydro simulation. I also analyzed the density and velocity properties of constant wind and dense shell model here. What defines a disk or a torus?
The viewing angle is the same as L2 pup, that is 83 degree. In this radiative transfer simulation, I assumed a 1 AU radius, 3000 K black body in the center of the simulation zone. I also assumed the dust/(dust+gas) mass ratio is 0.05. This is a bit high dust gas ratio but since the actual binary separation is 2 AU (in the simulation it is 4 AU), the optical depth may go higher if the separation shrink to 2 AU.
I do not know the particular species of the dust (probably silicates). I adopted the default dust species setting in radmc3d, I make the scattering to be isotropic even though the dust scattering is likely to be forward.
Further research can be done on:
- Redo 2 AU separation simulation.
- Figure out the dust species and their scattering, absorption properties.
- If there is enough computational resources, I want to do a larger simulation zone. Currently, for 2 AU separation, I can only go 40AU*40AU*40AU with 100*100*100 box plus 3 level AMR and it takes 120CPU*7days.
OutflowWind: Parker wind & Paper Figures
1. low density stellar wind
StellarEnvelope%omega=0d0 ! Assume planet is tidally locked ... lStellarGravity=false
movie of stellar density 10^-5;
Other densities and details can be found here —part 4.
2. Paper Figures
Will upload to the high res pictures to this page
Update (07/28/2015) -- Marissa
Shape
- Shape is now installed on all local machines! (Grass, Clover, Bamboo). It is installed in /opt/Shape/. In order for it to work, now all the machines have Java as well.
- Shape has a difficult time rendering large ASCII files.
Recall how Martin's data looks
3D Module | Render Module |
Martin's ASCII file has 4,803 lines. Note that the 3D Module illustrates if you have properly imported the simulation into Shape, where as the Render module is the actual emission map of the simulation. In the latter module, one has to go through a rendering process. This generally took a few seconds for Martin's data.
Baowei & Bruce's Data
Here is how it looks Visit:
In the 3D Module in Shape:
Now rotated to get another perspective:
Notes:
- The Right view in Screen Shot 1 looks like what might be the face of a box.
- However from the other views in the 3D Module, we can see that the box is not filled in. Wouldn't we expect to see something that does not look so "flat."
- For instance, Martin's data in the 3D Module sort of looks like the general shape that we would expect it to look like in emission map form.
- Perhaps the ambient is too strong?
- Perhaps Shape has a hard time loading in all of the rows?
The fact that we see *something* in the 3D Module that is not a mesh object that you import the data on is a good sign that I have imported the data correctly into Shape.
So at this point I would go on to the Render module to get the emission map. However when I fix my parameters and I am ready to visualize the simulation, Shape has a hard time rendering such a large ASCII file. So I cannot verify if there is a comparable view of the simulation to how it looks in Visit. If Shape were able to render this simulation, I would rotate it to see if there is a view that yields comparable results. However I cannot do that.
The ASCII file that Baowei gave me is in the format px,py,pz,n,vx,vy,vz,T
, i.e. position components, density, velocity components, and temperature. In Martin's simulation I did not use temperature. Although when I import the file with or without the temperature column, it still looks the same. So I do not think there is anything wrong with how I am importing it into Shape.
So now let's take a look at the parameters of the simulation (thanks to Baowei):
nDim = 3 GmX = 64,64,64 MaxLevel = 0 LastStaticLevel = -1 GxBounds = -32d0,0d0,-32d0,32d0,64d0,32d0
So we have a 3-dimensional simulation, which is good. MaxLevel=0 will help. However it seems that the GmX and the GxBounds are kind of big. We'd have 643 = 262144, which is close to the number of lines that we have in the file: 254727. Keeping the right proportions with a smaller GmX, say 163 = 4,096. Seems more reasonable.
Wires
- Visualizing on BH2 for the wire simulations seems to be a roadblock.
- Jonathan suggested using
vnc_start -t 480 --mem 60gb
however I wait for my job to start for a long time. I have never successfully gotten it to get me onto the node. Maybe other people are using it? - debug queue on BH2 (
'interactive -p debug -N 32 -t 60'
): Cannot do contour plots of any frames. It'll register the object you want to visualize in red text. However for that same frame, I can take a slice and visualize it fine. However say I want to visualize the next frame as a slice, then that frame also has the dreaded red text of death. However if I submit a new job and visualize that frame from the get-go, it works. However there is a similar story for the next frame, etc. - Visit with Clover, Grass, Bamboo: with any of the local machines seems to be no good.
- The maximum space for an AMR=2 chombo that I am trying to visualize is ~50 GB, for one with AMR=1, it is 13GB. The same story of visualization woes holds for both levels of AMR…
- I don't think the simulation is corrupt as if it was, then I couldn't visualize any frame after I submit a new visualization job.
- Jonathan suggested using
mass 2 flux ratio
Vazquez-Semadeni 2011 say for a uniform cylinder with a mean magnetic field
, the criticality condition (i.e. the critical mass to flux ratio) in terms of surface density and the field strength is given by:
where sigma is the surface mass density, i.e.
and L is the length of the cylinder.
This can be rearranged for the critical length of a cylinder of mean magnetic field strength B0 and number density n,
Given my initial field is
, and n = 1 cm-3, this critical length is:
This means that for the smaller box simulations, the colliding flows are sub-critical. However, for the larger box simulations they are super-critical.
Now, we do not see global gravitational collapse, so this measure of super-criticality doesn't seem to mean much. It doesn't take into account the kinetic energy in the flows, and the turbulent and thermal pressures that would also be counter-acting gravity. So then, since we do not see collapse even though the purely magneto-gravity critical condition is satisfied means there are other stronger forces in the cylinder preventing collapse.
If instead of putting the initial density in the above equation and we instead looked at the length scale of a cylinder with density = perturbation density, what do we get. First, let's take the average cloud density to be 100x the initial density. Then, we get,
This seems to say, if we assume the collision region is shocked to a mean density of n = 100 n0, and the strength of the field stayed the same as the initial value (which IS an assumption given the field also is shocked), then the length of an unstable collision region cylinder would be 1.5 pc. This is about 10x smaller then the collision region seems to be from eyeballing the cdm.
For the case of subcritical cylinders (i.e. smaller box simulations), how long does it take to acquire enough mass for the cylinder to become supercritical?
Can rearrange the above criticality equation to instead put in terms of the column number density, N:
Setting this equal to a function for N(t),
and solving for t gives:
where,
and
Solving for time gives,
This says the field is dynamically weak and SHOULD not be able to prevent global collapse by the time t = 2.6 Myr at the mass influx of the simulation. However, our simulation time is 30 Myr. So, something is supporting the gas against collapse.
Next, compare turbulent jeans timescales? If that is much longer than this timescale (and the simulation time), might be able to argue that the reason we're not seeing collapse is becuase of turbulent support?
meeting update
I've been working on my paper and am thinking of goals for my upcoming trip to MPI.
What defines a disk or torus?
I picked out the x-axis data (density and velocity in y-direction, which is a measurable data in Doppler shift) then I averaged over time.
The secondary locates at 2.67 AU away from origin and the primary locates at 1.33 AU away from origin. So the density spike within 2.7 AU is natural. The velocity in y direction is very spiky. I do not know why.
I do not think this is a disk or torus, people can be more clear when I compare this data to the dense shell ejection model.
Constant wind result, interpolated, averaged.
Dense shell result, interpolated, averaged.
Update 07/20/2015 - Eddie
Emission Maps
- Run C is the only one left that needs emission maps.
- More movies posted: MachStems
Clump Instabilities
- 4 processes: transmitted shock compresses clump, RT, KH, radiative cooling
- cloud-crushing time: where is the clump/ambient density ratio, is the clump radius, and is the incident shock velocity
- RT time: , and the dominant mode is of order the clump radius, so the RT instability will grow on time scale equal to
- KH time: , and again the dominant mode is of order the clump radius so the KH instability will grow on same time scale
- Radiative cooling is not occurring within the clump because the clump is already very cold (2.5 K) to satisfy pressure equilibrium with the ambient. Even with the transmitted shock, the clump will not reach a temperature large enough to trigger radiative cooling (T > 1000 K).
M | (yrs) |
---|---|
7 | 40.66 |
10 | 28.46 |
15 | 18.98 |
The simulation time is 75 years, and references state that a clump is typically compressed in 1-2
Bow Shock Instabilities
- Vishniac instability: typically cited for blast waves. Occurs when a slab bounded by a shock and a contact discontinuity is thin enough.
- NTSI: usually referred to in colliding flows. A cold slab bounded by two shocks will undergo NTSI if the slab is thin enough.
Tasks going forward (07-20-2015) -- Marissa
Wire Simulations
- Transfer files
- 2D Simulations on BH2 —> grassdata
- Level 1 Hydro Run —> BH2
- Visualize Level 1 Hydro Run
- Keep submitting jobs for Level 2 Hydro Run
- Use code made with different solver for level 2/1 MHD Run
- Blog post comparing the different solvers
Shape
- Visualize Baowei's output from his code
- Spherical shell problem.f90
- Testing with Baowei's code
- Visualize output in Shape
- Documentation of Shape
- Making movies
- HDf5 —> ascii (pipeline of Baowei's code)
Constant AGB wind
Isothermal Temperature: 400 K
Primary mass: 1 SM with alpha=0.134 (behave like a 0.134 SM star to gas)
Secondary mass: 0.5 SM
Escape velocity to the primary star: 14.7 km/s
Separation: 4 AU
AGB wind:
24 km/sIntro 1st draft
The interstellar medium (ISM) is a dynamic and ever changing environment. It exists in various thermal states, appearing to cycle through them continuously - from a hot, sparse phase to a cold, dense phase and back again. The evolution of the gas through these phases paints a picture of star formation, one in which the molecular clouds from which stars form are not long lived steady state structures as once believed, but rather the transient consequence of various instabilities that arise in the ISM (cf. Heitsch et al 2005 for a quick discussion of the various instabilities).
The picture of transient molecular clouds (born out of the cold, dense phase of the ISM) is the foundation of today’s theory of star formation, and is supported by many well known lines of evidence, such as: the vast majority of (if not all) molecular clouds are in the process of forming stars, these stars are young (less than 5 Myr), and molecular clouds contain a wealth of hierarchical structure that would not last long due to star-star scattering or tidal interactions (Elmegreen 2000, Hartmann 2001, Ballesteros-Paredes and Hartmann 2007, and references therein). Indeed, this idea of short cloud lifetimes and its connection to fast star formation has its roots in work that now dates back at least many decades (Hunter 1979, Larson 1981, Hunter et al 1986). As computational methods have become increasingly more powerful, these theories have been tested in simulations designed to follow the highly dynamical evolution of molecular cloud formation.
One such model that has gained widespread recognition is the ‘colliding flows’ model. This model provides a mechanism for generating clouds with the above mentioned characteristics. Importantly, the collision between 2 supersonic streams of gas has been shown to: produce non-linear density structures with ease, the density of which can grow large enough, for long enough time to allow for H2 formation, after which the molecular gas undergoes localized gravitational collapse (i.e. numerical star/cluster formation), and disperses within roughly a dynamical time (Heitsch et al. 2006, and others).
While these models are highly idealized, they are not without motivation. Coherent large-scale streams of gas are plausible in many situations in the ISM: expanding bubbles driven outward from energetic OB associations and/or supernovae, turbulent motions, density waves in the spiral arms of galaxies, and cloud-cloud collisions have been proposed (Hartmann, Ballesteros, Bergin (2001), PPVI chapter 1, Inutsuka 2015). Further, observational evidence is beginning to provide support for such mechanisms - Looney 2006 has demonstrated cloud-cloud collisions, atomic inflows surrounding molecular gas has been observed in Taurus (Ballesteros-Paredes, Hartmann, & Vazquez-Semadeni 1999) and other molecular clouds (Brunt 2003), and molecular clouds have been found at the edges of supershells (Dawson 2011, 2013).
With the large body of data supporting large scale magnetic fields in the ISM (cite), the addition of magnetic fields to these models was necessary (cite sims). Such models have identified additional instabilities formed in the flows (cite), and an overall reduction of SF is found (cite). As a natural next step, models are beginning to study how oblique shocks at the collision interface affect cloud and star formation (cite).
As the flows themselves are not likely to be exactly parallel, the question of oblique shocks is an important one to address. Such oblique flows can arise in X. As flow enters the oblique shocks, a shear velocity field is formed giving rise to the question - how does added shear affect cloud and star formation? Recently In Ostriker et al. …. In the present paper we wish to explore this parameter regime, i.e. of a magnetized oblique colliding flows.
Magnetic ring evolution
Consider the magnetic field to be similar to an atmosphere where we balance magnetic pressure gradients with magnetic tension forces? Than we have
which gives us
We can then integrate this to get the net tension force
and dividing by the inner area gives
which we can combine with flux freezing
and
Which can be reduced to
And then can be combined to
which eventually gives (if I haven't made any mistakes…)
although clearly I have since
Update 07/13/2015 - Eddie
Not much to add since late last week, just a few things:
- progress on generating the rest of the emission maps is steady and on-going (~1 week until they are all finished): MachStems
- updated the page with a couple of helpful diagrams
- I'll put up some more movies tomorrow.
- will also start generating some column density movies this week.
- re-reading Pat's 2011 3rd epoch paper in closer detail
Weekly Update (07-13-2015) -- Marissa
Mainly finished this page in my Shape guide : Importing External Data to Shape
Still working on other parts of my guide. Will post things for Jonathan's stuff tomorrow. Got my soft limit bumped to my hard limit as well.
Meeting update
- OutflowWind Parker solution:
- The result of Parker wind of the star looks stable as shown in the movies of blog:bliu07092015.
- The planet still have problems as the radii is too small/resolution too low and cannot see the temperature profile.. Running with larger radii
- OutflowWind: high-res 3D co-rotating
- In the development branch, Jonathan updated the 3D temperature profile in outflow object by moving the sun/day side was moved to +z direction comparing +x direction. But there's a problem in the temperature profile currently as shown below… Working on debugging it..
- high res runs with smaller omega on BG/Q: fixed some compiling bugs on bg/q. Running the old version of code while debugging 1…
- 3D pnStudy
progress
I have written a 1D hydro code using Fortran and python together. Python is used to call the excutable and fortran provide the numerically compiled exe.
This is a 1D parker wind.
My goal is to write a 1D radiative hydro code with dust formation and cooling.
I did a constant wind binary simulation last week. The escape velocity for AGB star is approximately 15 km/s. The wind is 24 km/s.
We see gas is not escaping from the system due to the deepened gravitational well. Gas material is more inclined to stay in equatorial plane and it is infalling. It does not have polar outward flow anymore but the density in polar direction is low.
meeting update
I got an estimate for the ring last week but it seems to be off by a couple orders Will share at meeting tomorrow
Working on text trimming and intro writing as well.
Some more thoughts on magnetic ring evolution
I will start with a definition of ram pressure for a cylindrical flow
where r_0 is the initial radius of the cylinder, \rho_0 is the post-shock density in the cylinder and c is the post-shock sound speed.
To get the field that is being stretched by the flow we use flux freezing
where S refers to the area the flux threads.
This gives us
If we now use the tension "pressure" going as
then we have ram pressure = magnetic tension leading to
or the balance radius r_b
OutflowWind: Redo stellar wind with new Park solution
This is to redo the stellar wind Park solution (as in blog:bliu07012015) with Jonathan's updated park wind data (blog:johannjc06252015) and compare with the old one…Other than the initial data, the new run has larger box (larger than Sonic surface) while the boundary for old one is inside the supersonic area…
density | ;density plot | ;density plot | |
Vx | ;Vx plot | ;Vx plot |
2D results:
pnStudy: clump bubble
Modified the density for area r<Rjet from
if (outflowType == clump .AND. r2.lt.Rjet) then q(i,j,k,1) = namb/nScale+& njet/nscale*(1d0-(r2/Rjet)**2) ![cu] ... end if
to
if (outflowType == clump .AND. r2.lt.Rjet) then q(i,j,k,1) = njet/nScale ... end if
to make the density inside equal to nJet while it was a value of stratified…
And the results look like this
time | rho_scaled | |
t=0 | ||
t= 0.003 C.U. |
Adam's Edits of Erica's paper
My edits
Dense Spherical Wind followed by AGB wind
I used 400 K isothermal gas temperature. In The dust disk and companion of the nearby AGB star L2 Puppis, they infer that the inner rim of the circumbinary disk is about 1000 K - 1200 K at 6 AU away from the AGB star. So 400 K is approximately the temperature at 24 AU - 40 AU. Since my simulation extend to 80 AU, I think 400 K is somewhat reasonable.
The gas is assumed 100% hydrogen atom. The wind is all spherical this time. The slow dense shell bear initial velocity of 16 km/s and mass loss rate of 5E-5 SM/yr while the AGB wind has initial velocity of 32 km/s and mass loss rate of 1E-7 SM/yr. The dense shell ejection last for 11.6 yrs and the AGB wind last for the rest of time so the total mass loss in this simulation is approximately 6.5E-4 SM. In the L2 Pup paper, the author estimated the dust weigh 1E-7 SM, if the gas dust mass ratio is 100, it implies that the mass of the disk is about 1E-5 SM. Usually, most of ejected material will leave the system and only a fraction stay in the circumbinary disk, I hope this fraction is 10%.
There is no dust and radiation transfer in my simulation so this initial velocity compensate some absence of the accelerate mechanism.
The AGB star is 1 SM and has radiative pressure on gas (I know the physics is much more complicate, e.g. radiative pressure on dust and dust and gas momentum coupling). The secondary is a compact accreting object which is different from the paper. I adopt the accreting object to make a hole in the center of the circumbinary disk, as I can keep removing gas from the disk. In reality, high temperature near the star can expel gas and create this hole. So temperature gradient is needed if we want to see a self consistent hole in the center.
Isothermal Temperature: 400 K
Primary mass: 1 SM with alpha=0.134 (behave like a 0.134 SM star to gas)
Secondary mass: 0.5 SM
Escape velocity to the primary star: 14.7 km/s
Separation: 4 AU
AGB wind:
32 km/sDense shell:
16 km/sWind ejection radius: 1 AU
Dense shell ejection duration: 11.6 yrs
compared to:
pnStudy: test different cooling after merging
This is to test the pnStudy module in 3D after merging with the latest development branch. The data files are modified from a run of 2.5D (located at /home/balick/tap30/t30AGBn4e2v200namb4e3/). Not sure if it's physically correct or not but just test if the code can run when setting different parameters for cooliing in physics.data…
The 2.5D result is for density and temperature while the 3D results are middle section of density only..
2.5D result from Bruce | |
3D; no cooling | |
3D; analytic cooling | |
3D; DM cooling | |
3D; II cooling |
OutflowWind: park solution for stellar wind
Currently the star is not stable
density | density line plot movie; | |
velocity | Vx along x movie | |
temperature | temperature movie |