There is now a function that can apply a second function to the 'grid'. For example:

  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
    CALL ApplyFunction(myclump, Info)
  END SUBROUTINE

  SUBROUTINE myclump(pos, t, dt, q, Ghost)
    REAL(KIND=qPREC), DIMENSION(:) :: pos
    REAL(KIND=qPREC), DIMENSION(:) :: q
    REAL(KIND=qPREC) :: t, dt
    LOGICAL :: Ghost
    dq=0
    IF (sum(pos**2) <= 1d0) THEN
       q(1)=10d0
    END IF
  END SUBROUTINE myclump

If you want to smooth your solution into a background state - then provided that 'q' comes in with a background state, you can just modify the function

  f = smoothing_function(pos)
  q(1)=(1-f)*q(1)+f*10d0

Existing objects can use the same routine - as long as the object type can be sent to the function. This is tricky to do in fortran - but is possible by casting the 'Outflow' or 'Clump' object as an 'Object'

  SUBROUTINE OutflowBeforeStep(Info, Outflow)
    TYPE(InfoDef) :: Info
    TYPE(OutflowDef), POINTER :: Outflow
    TYPE(ObjectDef), POINTER :: Object
    CALL TransferOutflowToObject(Outflow, Object)
    CALL ApplyFunction(OutflowSource, Info, IEVERYWHERE, lHydroPeriodic, Outflow%shape, Object)
  END SUBROUTINE OutflowBeforeStep

  SUBROUTINE OutflowSource(pos, t, dt, q, Ghost, object)
    REAL(KIND=qPrec), DIMENSION(:) :: pos,q, dq
    REAL(KIND=qPREC) :: t, dt
    LOGICAL :: Ghost
    TYPE(ObjectDef), POINTER :: Object
    TYPE(OutflowDef), POINTER :: Outflow
    REAL(KIND=qPREC) :: density, velocity, energy

    CALL TransferObjectToOutflow(Object, outflow)
    q(1)=Outflow%density
   ...

Now most objects have shapes that have routines that help transfer coordinates into the shapes orientation etc… But vectors also need to be rotated into and out of the shapes orientation. The momentum (and Magnetic) vectors need all 3 components to properly do the transforms - even if nDim==2. Since q does not always have all 3 components of velocity, this is a minor problem. One solution is to temporarily insert additional components of momentum into the q array. The routines that set the values for 'q' would need to explicitly use q(/ivx,ivy,ivz/) to access the momentum - and can no longer assume that the components of momentum are contiguous. The other option is to keep the routines that transform objects into a shapes reference inside the final function.

The other issue is what form of the fluid variables are stored in 'q'. Primitive is the easiest for the user to code.

• ticket #431 (tracking density and pressure protections) is ready to merge.
  • I tested it on one of my pulsed jet simulations and you can see it working.

movie​

• ticket #433 (tracking S species) is complete, and it will be ready to merge later tonight.
  • I thought there might be an issue because I saw some strange behavior in the emission maps from my 3-D 2-clump Mach stem runs. However, I think the cooling length was just poorly resolved, because increasing the resolution got rid of the strangeness.

movie​

• ticket #436 (tracking max speeds) still needs some work.
  • Work will be slow for the next few days because of the holiday, but I think I can finish this before next Monday.

Here: ​http://www.pas.rochester.edu/~erica/Talks/Fourier_Presentation.pdf

*I ran the simulation again with four times the box size and 30 days. The sonic surface can be seen moving outwards now and there is no inflow from the boundaries.

​http://www.pas.rochester.edu/~mehr/wikistuff/symbox.gif

Initial separation = 4 AU

Temperature = 1500 K

Primary mass = 1.5 solar mass but behave like 0.1379 solar mass star when interact with gas

Secondary mass = 0.5 solar mass

Mass conserved, momentum conserved.

By the way, at this resolution, I did not see the disk tilting. Maybe tilted disk is a resolution issue or illusion?

alpha = 0.12 case

temperature = 2000 K

*No temperature variation wrt angle. Temperature stays constant at the base of the wind.

​http://www.pas.rochester.edu/~mehr/wikistuff/symmetrical.gif

*Still failing to get supersonic outflow at the edges, and therefore the inflow at the end.

*Temperature profile

​http://www.pas.rochester.edu/~mehr/wikistuff/symmtemp.gif

*Working on including radial ambient outflow to match with the initial conditions.

• Divergence Problem with 2.5D MHD in the pnStudy Tried to test Bruce's idea of 45 deg taper:

   Martin has the option to launch collimated flows with a Gaussian taper  To run this you configure problem.data with an opening angle=90 (spherical flow) which Martin modulates with a gaussian of the form exp-(latitude/user-specified-gaussian angle)^2.  You find that modulation function somewhere.  I should think that its form can be changed from
   exp-((phi)/W)^2,                   where phi is the latitude angle and
                                            W is the user-specified gaussian width called "tf" in problem.data
to this:
   exp-((phi-PA)/W)^2,                where PA is the user-specified flow angle

PA might as well default to 45 deg for cones and gaussian taper sims (that is, Outflowtype=2).

But found this divergence Mag issue on bluehive:

​Movie

Need more time to confirm the results and dig out why….

• 2.5D Result from alfalfa

​2.5D Movie

• Test for 45 deg taper with 2D runs (taper15n4e2v200namb4e4)

Fixing a bug just found. Will update results soon…

Binary separation = 7 AU

Temperature = 2000 K

Primary mass = 0.8 AU with alpha = 0.1723 that is, primary star behave like an 0.1379 solar mass star. This is presumably due to radiation pressure on dust grains (so the temperature should be like 1500 K)

I am not sure about the mass loss rate of the primary star.

Secondary accrete gas with initial mass of 0.5 solar mass

First picture shows the mass gain of the secondary and the separation decrease of the two stars

​midplane

This gif file is the mid plane density plot of the whole simulation.

People should not take above result too seriously. Many things, for example, mass conserved primary star, self-consistent boundary condition of the primary star, tidal effect and tidal force, non-isothermal situation, should be taken into consideration.

The following are simulations with mass conserved primary.

Red line is the theoretical angular momentum computed from the instantaneous separation (that is, the blue line) of the two stars and their corresponding mass. It is calculated by:

L=m_1 m_2 \sqrt{\frac{G d}{m_1+m_2}}

Green line is the result of angular momentum from AstroBEAR.

Initial separation is 4 AU, primary mass = 0.8 solar mass with alpha = 0.1723, secondary is 0.5 solar mass. Gas temperature is 1500 K.

The orbit changed to elliptical one since I turn on accretion and mass conservation mechanism suddenly, so it serves as a perturbation. Distance, theoretical angular momentum and computational angular momentum synchronizes afterward. We can see that the mean slope of green line is lower than the red line. That is to say, part of the orbital energy is changed to some other kind of energy. Since the simulation is isothermal, it can only be transformed into gas's kinetic energy or its gravitational potential energy.

I forgot to deduct momentum of the primary in simulation, that is, its momentum is not conserved, that is why angular momentum is increasing.

This is a revisit for the new refinement for pnStudy module (see blog:bliu11182014). This is Bruce's idea:

"My sense is that the wiggles develop only after the rim emerges from the fixed high-res region of the sim at 0 < x < 0.5 and 0 < y < 0.5. To check this out, is there a way to reconfigure the high-res zone so that it is three times thinner and taller? "

Below is the results — the mesh is kept to show the rectangle refinement works. The results show clearly the difference for areas along x and y axis:

With Mesh	No Mesh
​movie with mesh	​movie without mesh

• As discussed earlier, I removed the temperature asymmetry from the problem by commenting out the function we wrote to vary it as a function of angle.

Problem As the image from the initial cycle shows, the right side is still colder than the left.

*Also, after cycle 10 the color bars at the base of the wind do not change, which means the temperature is being maintained as we need, but the contrast between the day and the night is unclear. I am using the unedited version of outflows.f90 for this and I still get the contrast.

Subsonic flow problem *We were not getting supersonic outflow at the boundary. I tried various fixes including increasing the box-size, ambient temperature and density but the mach number continues to stay below one for one side and supersonic for the other. Below curve A is for the left side.

Later on we do see supersonic speeds but it's mostly rightward inflow as shown by the contours below.

Initial Conditions I ran these simulations for lambda=5.3 as done in the paper. The only part where my initial conditions differ from the paper are the direction of the ambient velocity. I gave it a y-directed initial sound speed. I am looking into how to make it radial by modifying the code.

Movie ​http://www.pas.rochester.edu/~mehr/wikistuff/movie.gif

Shear Angle	Final Frame/2 (137)*	Final Frame (273)*
0
15 *
30
60

Table. Magnetic Field vs. Density plots for all shear cases of the "weak-field" simulations (beta 10). The sloped straight line is the 10K temperature, the top horizontal line is the ram pressure, the lower the magnetic pressure. The curved function is the density verses pressure line. Here we have also included a "grey" area, which is the mass weighted density versus pressure post-processing.

* The Beta 10 Shear 15 case does not have a final frame of 273 yet. We are still waiting for the untarring of the rest of the chombos from Stampede. Thus for now I've included 217 and 109 for the middle and end visuals. Will update the table once we get up to 273.

I implement the new refinement idea we discussed yesterday: refine only the outflow launch region to the highest level and see if it helps for the weird wiggles close to y-axis. This run has 6 level AMR refinement (equivalent to 7 level as the box is ½ of runs before) only on the outflow object and maximum 3 level AMR outside the outflow. It also uses Rho as refinement variable instead of px and py… Just as Bruce said, the resolution doesn't help much

​Movie

• Got rid of the inflow from the left and the unexpected bow-shock.

Movie ​http://www.pas.rochester.edu/~mehr/wikistuff/movieb.gif

Density

Mach number -Still seeing higher values on the wrong side

• Finished 2.5-D pulsed jet paper revisions. The referee had asked about how a convolution with the instrumental psf would change my emission maps, so I wrote a program to calculate a convolution. Here are images of the beta = 1 model:

• The 3-D Mach stem simulation with cooling finished. Below is the density movie.

movie​

• Need to work on HEDLA paper revisions.

• Also working on setting up camera to make angled emission maps of my 3-D Mach stem runs.

Ran the original simulation for a longer time (25 days) and wider bounds. Inflow from the boundaries persists.

​http://www.pas.rochester.edu/~mehr/wikistuff/movietime.gif

Gravity Only

Testing whether the point gravity is working correctly.

​http://www.pas.rochester.edu/~mehr/wikistuff/movie_go.gif

New Initial Condition

Ambient has a fixed outward radial velocity. No difference from the original.

One problem is that the outflow speed does not seem to be exceeding the escape speed.

Wind Only

I keep getting the following error if I comment out the ambient to look at a strong wind only.

Fixed grid planet simulation also experienced explosions. Was going to do FFT self-gravity, but that is not fully implemented.

And MATLAB is taking 24 hours to generate new hydrostatic profiles for some unknown reason. Need to reboot my system.

Also, in order to verify the projection routines were working - I made a movie of a simulation of a wire box. Click the image to see a movie.

Here is my Grid Initialization routine…

  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
    INTEGER :: i, j, k, l, mb(3,2)
    DO i=1,3
       DO j=1,3
          IF (i == j) CYCLE
          DO k=1,2
             DO l=1,2
                IF ((Info%mGlobal(i,k) == GmGlobal(i,k)) .AND. Info%mGlobal(j,l)==GmGlob\
al(j,l)) THEN
                   mb(6-i-j,:)=(/1,Info%mX(6-i-j)/)
                   mb(i,:)=(GmGlobal(i,k)-Info%mGlobal(i,1)+1)
                   mb(j,:)=(GmGlobal(j,l)-Info%mGlobal(j,1)+1)
                   Info%q(mb(1,1):mb(1,2), mb(2,1):mb(2,2), mb(3,1):mb(3,2),1)=100d0
                END IF
             END DO
          END DO
       END DO
    END DO
  END SUBROUTINE

November 17-21:

• Write a script to automate jobs.
  • Update (11/18): Completed.
• Run astrobear in post-processing mode on all of the Beta 10 runs (Shear 0, 15, 30, 60), as well as the convergence tests for Beta 10 Shear 15 (levels 2, 4, 6) for the B vs. n plots. Done twice for the mid-way and end of the simulation. Exception: Beta 10, Shear 60, where we have frames out to 328, and the convergence tests which are at frame 50. This yields 12 simulations in total
  • Update (11/18): Finished shear 0, and 30, and 2/3 of Shear 60. Need to get other chombos for Shear 15.
• Make B vs. n plots for the different runs (convergence tests not necessary), publish on wiki.
  • Update (11/18): See my blog post for the first set: ​madams11182014
• Make B vs. n figures (for the paper).
• Run code for beta maps with projected streamlines, 3 times (10.1, 20.1, 27.2 Myr) for each production run, the final frame of Beta 10 Shear 60 (i.e. frame 328), and all of the convergence tests for Shear 15 at frame 15.
  • Update (11/18): Reservation for tomorrow.
  • Update (11/19): Submitted jobs, waiting in queue for noon to roll around.
• Make plots of the beta maps and the projected streamlines, publish on wiki.
• Make figures for the beta maps.

Waiting on Erica before doing the following:

• Run doe for energy spectra on all the production runs.
• Make plots of spectra, publish on wiki.
• Make figures for paper out of best spectra.

November 24-28:

• Run code for density and velocity spectra on the different runs.
• Make plots of spectra, publish on wiki.
• Make figures for paper out of best spectra.
• Run filament analysis on the 3 or 4 best filaments of each of the runs at end time

Fly Through Movies to Make:

• Movie of camera circling around the collision region for one of the colliding flows cases.
• Make a movie that travels down the barrel of the flow, reaches collision region and circles around.
• Make a movie that travels down the barrel of the flow, reaches the collision region, circles around and then zooms into a smaller region of interest to observe its evolution
• Make a movie that travels down the barrel of the flow, reaches the collision region, circles around, zooms into a smaller region of interest and circles arond to watch its evolution, then pans out to circle around the collision region. Do this for each run.

Fig1 shows the roche-lobe and 5 Lagrangian points. The red ring is the boundary of the primary star. Gif2, gif3 and gif4 are animations. Fig5 is the mass gain of the secondary.

I use thermal-driven wind model, which is parker's solution of isothermal wind. I have checked the wind speed for single star case, it is about 30% higher than theoretical solution at this resolution. So the mass loss rate will not deviate too much.

The mass loss rate of a single primary star (without the secondary) is 2\times10^{-9} solar mass per year. The mass of the primary star is 0.1379 solar mass, secondary is 0.1 solar mass. Isothermal temperature is 2000K, sonic radius is at 7.5 AU, also the wind should reach its escape velocity at about 9 AU according to Parker's solution. The separation of the two stars is 4 AU. The scale in Roche-Lobe picture is 1AU per unit length, that is, red ring has radius of 1 AU.

But when I add the secondary there, the accretion rate is 1\times10^{-8} solar mass per year. I think this is reasonable because the secondary is in the subsonic region and its existence is affecting the structure of the wind, or the stellar structure of the primary, which leads to much higher mass loss rate.​ The reason that the mass loss can behave very differently to the solely thermal-driven case may be the secondary is attracting the gas therefore removing the gas from the boundary much faster than thermal-driven wind. Then the more gas will be pushed from the boundary due to pressure (Ideal gas).

rochelobe.pdf​

45degreeview.gif​

equatorialview.gif​

z0.gif​

massgain.pdf​

Trying to figure out how the grid size affect the wiggles developing near the x and y axes for the Bruce's run with the spherical winds into a spherical AGB environment. From the image below it seems the AMR works as expected. While the finest refinement is clearly around the spherical wind area, things seem better for larger Rjet than Rjet=0.5.

• Reproduce Bruce's results with Rjet=0.5 (sphAGBn10v200namb4e4)

	Zoom-in image

• Rjet=5000AU(or 10 cu) and other parameters are same

t=0
t=100y

​Mesh movie for Rjet=5 or 10cu and ​Zomed-in on the edge

Here is a first cut of the Bvn plots. Y axis = Magnetic Pressure, X axis = number density:

B1 60

Density weighted:

Volume weighted:

B10 60

Density weighted:

Volume weighted:

B10 0

Density weighted:

Volume weighted:

Trying to figure out the ambient temperature where close to the jet drops for the pnStudy module (blog:bliu11102014). This hasn't been found before (From Keira's runs). Here I tried to reproduce one of Bruce's runs and also a run of Keira. For Bruce's run, the ambient temperature is 1000K, and the temperatures close to the jet drops to 100K around 660y. For Keira's run, the ambient temperature is 100K. And the temperature of similar area doesn't drop — as listed in 1 and 2 below.

A. Cooling Floor Temperature

Jonathan helped me figure out the cooling played a role here. The floor temperature for the cooling is 100K and when it reaches 100K it shuts off the cooling. That's why in Keira's case, it won't see the drop. But whenever the ambient temperature is higher than 100K, there is a drop

B. Misleading comments in problem.data

It was also found the comments for the temperature parameters in the problem.data was misleading. Here's the current problem.data

!      BACKGROUND or “AMBIENT” SECTION. Values apply to origin
tamb = 1d3           ! ambient temp, 1cu = 0.1K (100K=1000cu)     

...

tjet  = 1000d0      ! flow temp, 1cu = 0.1K (100K=1000cu)

Both of these should be in Kelvin and 1 cu = 10K

!      BACKGROUND or “AMBIENT” SECTION. Values apply to origin
tamb = 1d3           ! ambient temp in Kelvin ( 1cu = 10K)

...

tjet  = 1000d0      ! flow temp in Kelvin ( 1cu = 10K)

C. Results

1. Reproduce Keira's run with current code

Main Parameters	Baowei's Run	Keira's Run
tamb = 1d2; namb = 4e2; outflowType = 1; njet = 1d2; Rjet = 1d0; vjet = 2e7; tjet = 1d1;

2. Reproduce Bruce's run

Main Parameters	Baowei's Run	Bruce's Run
tamb = 1d3; namb = 4e4; outflowType = 1; njet = 4d4; Rjet = 2d0; vjet = 2e7; tjet = 1000d0;

3. Test for different ambient Temperatures — all other parameter are same

Tamb =1d2
Tamb =2d2
Tamb =3d2
Tamb =4d2
Tamb =5d2

3d.gif​

Wind accretion. Boundary set to be outflow only.

Accretion has been turned on.

Feedback from gas to particle is turned off.

• Had some issues with running 3-D Mach stem problem with cooling.
  • It seems that the code can only use the exact Riemann solver with the ideal equation of state (and not a multi-species EOS which is what I use for Z cooling). The code immediately gets restarts due to nan in flux.
  • There are 11 tracer fields: one for each clump, one for the wind, HI, HII, HeI, HeII, HeIII, SII, SIII, SIV. It's likely that the nan is being caused by one or more of these.
  • For now, I have resubmitted the run to the queue using the default (HLLC) solver.

• On the afrank nodes, my 3-D Mach stem with cooling run failed after 7 frames. I got srun errors that I didn't really understand. I'm thinking it's a memory issue, so I switched back to the standard queue. Here is what I have so far:

• Altered my problem module to use 3 clumps now, and also implemented a camera object so that I can make emission maps based on inclination angle. This will be useful for comparing with observations. This just needs to be tested now.

Here's a review of everything that has been accomplished up till now on the S&P problem.

*Introduced asymmetry wrt initial temperature in the outflow module.

*Initialized the ambient at low density and sound speed to allow for a thermally driven wind.

*Ran the simulations for a period of 10 days, temperature values of the order of 10,000K and planetary sizes of 1-2 Jupiter masses. *All having to do with different values of the hydrodynamic escape parameter

*Best results seen in the range of lambda 4-6. Robustness issues.

Following are some plots and a comparison with the S&P results.

lambda = 5.67 Density

*In steady state the pattern matches well with S&P results; purely radial flow at large distances and circulation at small distances.

*Major difference is the position of the shock. I do not see a discontinuity in the theta direction. Instead it's visible between the ambient and the planetary outflow. S&P get a shock discontinuity at about theta = 3pi/4

*Cuts at theta = 0 (K) , pi/2 (L) and pi (M). As expected.

Temperature *Day-night temperature ratio 20,000:200K. S&P kept a ratio of 100:1.

*The position of sonic surface is well in agreement with the S&P result occuring at almost 2*planetary_radius.

Main issue with code *Only gives sensible results for carefully chosen parameter values. Compare the movies below for lambda=5.67 and lambda=6.11 ​http://www.pas.rochester.edu/~mehr/wikistuff/movienew.gif

​http://www.pas.rochester.edu/~mehr/wikistuff/bad.gif

Next Tasks *Finish the write-up. *Owens Adams Paper *Move to 3D regime

Ran 64³+2 simulations of planetary atmospheres and looked at

• fast spinning (10x )
• lower resolution (64³+0)
• higher tolerance (1e-16 instead of 1e-8)

And yes - they all blow up

These are not all at the same time.

Fiducial	Fast Spinning
Low res	Higher Tolerance
1e-6 density contrast	1e-9 density contrast

Ambient temperature is set 1000K initially. As jet moves, ambient temperature at the bottom drops to 100K. Results are similar for both the new code (vJet not depends on Rjet) or the old code (v=vJet/Rjet…) — the jet velocity for the old code is half of the new code due to the Rjet dependence.

	New code(5AMR)	Old code (3AMR)
t=0
t=600y
movie	​New code:jet&Abient Temperature	​Old code:jet&Abient Temperature

• Write script to automate analysis of the high resolution colliding flows runs, which will be useful for spectral analysis (see: ​Erica's blog post)
• Run astrobear in post-processing mode on all of the runs (B10 S60, 30, 15, 0) for the B vs. n plots.
• Make B vs. n figures for the paper of the different runs.
• Make column density figures for the paper of each of the different runs stated above including B1S60 and convergence tests.
• Make a movie that travels down the barrel of the flow, reaches the collision region and circles around (i.e. implement the fly through cameras for our colliding flows data)

Took the wind out of the problem module and added "res" and "temperature," which can be edited for each simulation in the problem.data. Temperature is added to the clump and the res is added to the "Projection Movie." In both simulations they have a res=300 even though they have 1 level of amr. You can see you can see the affect by adding mesh in visit (the image with the mesh is from the second visualization):

​movie

problem.data file for this run:

&ProblemData
rho=10
radius=1
temperature=.1
res=300
ncameras=5
/
&CameraData
pos=-10,4,-10
focus=2,0,0
upvector=1,0,0
time=0d0
/
&CameraData
pos=-5,4,-10
focus=2,0,0
upvector=0,1,0
time=.25d0
/
&CameraData
pos=0,4,-10
focus=2,0,0
upvector=0,0,1
time=.5d0
/
&CameraDatasee by adding the mesh to the visualization
pos=5,4,-10
focus=2,0,0
upvector=0,1,0
time=.75d0
/
&CameraData
pos=10,4,-10
focus=2,0,0
upvector=1,0,0
time=1d0
/

Here we have 5 cameras with a starting position of pos=-10,4,-10 that is viewing the clump in the box with an up vector in the y-direction until frame 25 where it switches perspective and moves along the x direction toward positive 10. So the simulation looks like it is moving around while we move in the x-direction over time.

​movie

problem.data for this run:

&ProblemData
rho=10
radius=1
temperature=.1
res=300
ncameras=4
/
&CameraData
pos=-10,8,-10
focus=2,0,0
upvector=0,1,0
time=0d0
/
&CameraData
pos=-10,-8,14
focus=2,0,0
upvector=0,0,1
time=.33d0
/
&CameraData
pos=14,-8,-10
focus=2,0,0
upvector=1,0,0
time=.66d0
/
&CameraData
pos=-10,8,-10
focus=2,0,0
upvector=0,1,0
time=1d0
/

Here I am using 4 cameras. In this simulation I attempted to rotate the clump while zooming toward and away from it. The visualization is also continuous as the first and last camera have the same parameters. I think this is quite nice and I prefer this visualization to the first one.

I completed the same run as before, but with the exact Riemann solver:

movie​

B v. n plot-

The B v. n code is done - might just need to be tweaked for the best min/max combination for the magnetic pressure.

Here is a page on this plot:

​https://astrobear.pas.rochester.edu/trac/wiki/u/erica/CF_bvn_plot

Marissa, you can find the problem module here:

/grassdata/erica/CollidingFlows/scrambler_3.0/modules/Problem/problem_Bvn_pdf.f90

You will need to copy it to BS and compile it, and run it with a reservation. I set up a reservation for you to start next week, Wednesday. There are *2* pdfs of B v. n — one weighted by volume and one by mass. We can see which one we like best. Each job shouldn't need more than 256 cores I would say. Start by making quick and dirty plots of these for each of the runs including the convergence tests - at frame 100 & 200 please (for the convergence tests at frame 50). With the reservation you should be able to have them all run in parallel easily.

Also Marissa, start early in week on script - I will be back on Wed to look at the B v n plots with you.

B10_60 case -

Ran out to about 35 Myr and got a single sink to form - you can find an updated movie down the barrel for this case on this page - ​https://astrobear.pas.rochester.edu/trac/wiki/u/erica/B10_S60#no1 under attachments, named 'test.mpeg'.

Self-gravity ticket -

Made some progress on my ticket — will post to it when I return

I ran the simulation again for lambda=5.67 for a longer time (10 days, turning off outflow at day 5). The strange box-like artifact seems to have gone away and the steady state looks very close to the original paper.

​http://www.pas.rochester.edu/~mehr/wikistuff/movienew.gif

tracer.gif​

2.5D breakup.

3D boundary wind object is not good.

Week of Nov 3-7th
-Get working code for B v.n maps for marissa
-Get projection of field lines code
-Get beta map code for Marissa

Week of Nov 10-14th
-Get working code for energy spectra for Marissa
-Get working code for velocity and density spectra
-Get code for velocity dispersion
-Email Amira folks for trial license
Week of Nov 17-21
-Get Fabian’s code working and code up different radial analysis we want
-Work on characterizing filaments mathematically for use with Amira

Week of Nov 24-28
-Work with Amira

Test if set jet velocity not depending on the Rjet (blog:bliu11032014)

    v=vjet/velscale*(/0d0,1d0,0d0/)/fact*timef 
    ! v=vjet/velscale*(/0d0,1d0,0d0/)/Rjet*fact*timef !10 mar 2014

	New Results (3 levels AMR)	Old Results (3 levels AMR)
Rjet=2, vJet=200	;
Rjet=2, vJet=100	;	;

​No-Rjet_dependence Code:Rjet=2,Vjet=200;

​No-Rjet_dependence Code:Rjet=2,Vjet=100

​Rjet_dependence Code:Rjet=2,Vjet=100 effectively

z0.gif​

Today I ran seven simulations that play with the CameraData in the problem.data file. I focused on the up vector as well as the focus of the camera. For two cameras I played with a focus {1, 0, 0} and {3, 0, 0} as the previously done simulations have a focus of {2, 0, 0} for all cameras. The simulations I ran with three cameras have the second camera switch up vectors, as well as have cameras switch into ascending or descending order based on the three focuses we know already (1, 2 and 3).

Two Cameras

Run description	Image Final Frame	No. of Frames	CameraData File	Link to movie
2 cameras, both up vectors in y, with focus of 1 in x		50	​problem.data	​camera movie
2 cameras, both up vectors in y, with focus of 3 in x		50	​problem.data	​camera movie

Three Cameras

Run description	Image Final Frame	No. of Frames	CameraData File	Link to movie
3 cameras; first and last camera have up vectors in y; second camera has an up vector in x		100	​problem.data	​camera movie
3 cameras; all cameras have up vectors in y		100	​problem.data	​camera movie
3 cameras; first and last camera have up vectors in y; second camera has an up vector in z		100	​problem.data	​camera movie
3 cameras; all up vectors in y with ascending focus from 1-3 in x		100	​problem.data	​camera movie
3 cameras; all up vectors in y with descending focus from 3-1 in x		100	​problem.data	​camera movie

Runs are ¾ done.

Just moving data to local machines. Should have more movies up soon.

Here is the B10-60 case:

​https://astrobear.pas.rochester.edu/trac/wiki/u/erica/B10_S60

I am thinking we will want to do the bulk analysis on just the Beta = 10 runs, that is shear angle = 0, 15, 30, & 60

For the B1 run- maybe we can just have a small section in the paper showing morphology (column density maps)? This run is only out to t=20 Myr as of now. Don't think it is necessary to run out longer either.

For the resolution study also was just going to show column density map??

For the B v. n maps - thinking of plotting pdfs of the magnetic energy vs number density ? What units should this be in?

Lastly - thoughts on restarting B10 60 case for another 10 - 20 Myr?

• 3-D co-moving clumps, no cooling, gamma = 1.01, nclump/namb = 500, 40 cells/rclump, M = 5 shock for moving clump, difference in relative clump velocities = 10 km/s. Took 16.7 hrs on 96 cores on bluehive. Used default solver (HLLC).

movie​

• For 2.5-D pulsed jets, I am just going to regenerate figures with frame 96 and finish the paper revisions.

• Will also do paper revisions for HEDLA Mach stems paper this week.

• This is a busy week for PHY 101, but I should still have time to do development if we want to start having our weekly development sessions.

In the pnStudy, the velocity of the jet is set according to the radius:

!======== J E T=========:
IF (outflowType == collimated) then
  q(i,j,k,itracer4)=1d0
  fact=1d0 !10 mar 2014
  !fact=exp(-(x**2+z**2)/jet_width**2) ! b 4 10 mar 2014
  qjet(1)=njet/nScale*fact
  v=vjet/velscale*(/0d0,1d0,0d0/)/Rjet*fact*timef !10 mar 2014
  !v=vjet/velscale*(/0d0,y,0d0/)/Rjet*fact*timef !b 4 10 mar 2014
  qjet(imom(1:nDim))=v(1:nDim)*qjet(1)*& !ramp up velocity
  mom_flow !5 may 2014, time dependent mom flux requested by bbalick. 
  qjet(iE)=qjet(1)*tjet/TempScale*gamma7

This makes the Jet velocity is not vJet (which is given in problem.data) when Rjet !=0.

1. Rjet=2d0, vJet-2e7 as in Problem.data

outflowType  = 1    ! TYPE OF FLOW    1 cyl jet, 2 conical wind, 3 is clump
njet  = 4d4         ! flow density at launch zone, 1cu = 1cm^-3
Rjet  = 2d0         ! flow radius at launch zone, 1cu = 500AU
vjet  = 2e7         ! flow velocity , 1cu = cm/s (100km/s=1e7cu)
tjet  = 1000d0      ! flow temp, 1cu = 0.1K (100K=1000cu)
tt    = 0.0d0       ! flow accel time, 1cu = 8250y (0.02 = 165y)
open_angle = 00d0   ! 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

Here's the file of Scales.data

 TIMESCALE       =   260347122628.507     ,
 LSCALE  =  7.479899800000000E+015,
 MSCALE  =  2.099937121547526E+026,
 RSCALE  =  5.017864740000001E-022,
 VELSCALE        =   28730.4876830661     ,
 PSCALE  =  4.141950900000000E-013,
 NSCALE  =   300.000000000000     ,
 BSCALE  =  2.281431350619136E-006,
 TEMPSCALE       =   10.0000000000000     ,
 SCALEGRAV       =  2.269614763656989E-006

Velocity plot shows the jet velocity is 100 km/s instead of 200 km/s. And this can also be calculated from the distance the jet travels during 600 yrs:

14*1.5*E11 km / (660*365*24*3600) = 100 km/s

t=0	t=660 y

2. Rjet=1d0, vJet-2e7 as in Problem.data

outflowType  = 1    ! TYPE OF FLOW    1 cyl jet, 2 conical wind, 3 is clump
njet  = 4d4         ! flow density at launch zone, 1cu = 1cm^-3
Rjet  = 1d0         ! flow radius at launch zone, 1cu = 500AU
vjet  = 2e7         ! flow velocity , 1cu = cm/s (100km/s=1e7cu)
tjet  = 1000d0      ! flow temp, 1cu = 0.1K (100K=1000cu)
tt    = 0.0d0       ! flow accel time, 1cu = 8250y (0.02 = 165y)
open_angle = 00d0   ! 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

Scales.data are same

Velocity plot shows the jet velocity is 200 km/s. And this can also be calculated from the distance the jet travels during 600 yrs:

28*1.5*E11 km / (660*365*24*3600) = 200 km/s

t=0	t=660 y

3. Rjet=4d0, vJet-2e7 as in Problem.data

outflowType  = 1    ! TYPE OF FLOW    1 cyl jet, 2 conical wind, 3 is clump
njet  = 4d4         ! flow density at launch zone, 1cu = 1cm^-3
Rjet  = 4d0         ! flow radius at launch zone, 1cu = 500AU
vjet  = 2e7         ! flow velocity , 1cu = cm/s (100km/s=1e7cu)
tjet  = 1000d0      ! flow temp, 1cu = 0.1K (100K=1000cu)
tt    = 0.0d0       ! flow accel time, 1cu = 8250y (0.02 = 165y)
open_angle = 00d0   ! 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

Scales.data are same

Velocity plot shows the jet velocity is 50 km/s instead of 200 km/s. And this can also be calculated from the distance the jet travels during 600 yrs:

7*1.5*E11 km / (660*365*24*3600) = 50 km/s

t=0	t=660 y

3D Binary simulation can not restart.

z0.gif​

-Here are some scaled plots for values of lambda between 3 and 10. I initialized the ambient outflow as done in the paper and turned it off after 0.5 days, so that the dominant thermally driven results are visible. In the paper, they refer to it as the steady state.

The plots here are the density profile, temperature and density with the velocity/soundspeed contours. The transition from subsonic speeds at small distances to supersonic speeds can be seen.

Lambda = 3.67

Lambda = 5.5

​http://www.pas.rochester.edu/~mehr/wikistuff/movie5.gif

Lambda = 8

Will add cuts at different values of theta to determine the exact position of the shock discontinuity

Below I have a table of visualized simulations I ran of a clump and a wind from last week's blog posts:

​movie

For each of these little tests I changed the final time and frame number in the global.data file. There is a time for the camera in the problem.data file, so I am also interested in whether or not that time has to be equal to t-final or not and how that affects the simulation if it is changed. No matter the number of frames, the simulation will play out to the same point. Essentially the t-final in the global.data acts as a frame rate, or temporal resolution. You can see that 0.1 final time simulations have a wind that hasn't even hit the clump yet by the end of the run. For the 0.25 simulations, the wind hits it by the end, and in 0.5 the wind envelopes the clump. The number of files, bovs or chombos do not bear on what is "physically happening."

The "camera movies" are made in VisIt like how we would make a visualized bov movie. The camera movies are projections, however they are not necessarily "down the barrel," like what we are use to.

In the chombo simulations I took a slice down the y-axis in VisIt. We know for the camera movies the camera movies along the x-axis, which is does also in Visit. It is a projection of a 3D simulation so we do not need a "z" axis.

In each simulation there are two cameras. One is positioned at pos=-10,2,-10 and the other at pos=8,2,-10. So it is changing in the x-direction, but keeps focused on the same object as it moves away from it. Both cameras have a focus of 2 in the x-direction.

So the camera is positioned outside of the box and integrates the density through the domain - in essence it is a column density map.

Why I did this is because there is a time for each camera in the problem.data file. For the nth camera, should we have its time equal to the final computational time in the global.data file?

	0.1 Final Time	0.25 Final Time	0.50 Final Time
10 Frames
	​camera movie, 10 bovs	​camera movie, 10 bovs	​camera movie, 10 bovs
	​chombo movie, 10 chombos	​chombo movie, 10 chombos	​chombo movie, 10 chombos
25 Frames
	​camera movie, 25 bovs	​camera movie, 25 bovs	​camera movie, 50 bovs
	​chombo movie, 25 chombos	​chombo movie, 25 chombos	​chombo movie, 25 chombos
50 Frames
	​camera movie, 50 bovs	​camera movie, 50 bovs	​camera movie, 50 bovs
	​chombo movie, 50 chombos	​chombo movie, 50 chombos	​chombo movie, 50 chombos