Posts for the month of December 2013

NLUF summary

Previous simulation setup:
http://www.pas.rochester.edu/~shuleli/NLUF_sep/setup.png

Temperature setup:
0.026ev in the ambient, 0.0012ev in the target initially, 0.026~0.8ev in the ambient, 0.0012~0.3ev in the target after preheated by radiation). The Spitzer resistivity has a floor temperature: any temperature below 0.26ev (3000K), is treated as 0.26ev when calculating resistivity. This means that only after preheating, the hot end of the target will start to diffuse magnetic field.

http://astrobear.pas.rochester.edu/trac/astrobear/attachment/blog/shuleli09072013/hydronew.png

Comparison animation for beta = 1, 10 and infinite:
http://www.pas.rochester.edu/~shuleli/NLUF_sep/movieall.gif

AstroBEAR has resistivity implemented and calculates the amount of magnetic diffusion at every simulation cell given their density, temperature and ionization fraction. We have since obtained ionization table of Argon and SiO2 from Jaques at LLE (Z(rho,T)), and have implemented code to read the table and compute local resistivity with the correct Z value.
Other design suggestions: eliminate the backwall shock and put in an ablator in front (density wall)

Meeting Update 12/16/2013 - Eddie

  • made all of the images and movies from the 2 clumps simulations ehansen12052013
  • also did some analytic work to show that the equation that I was using for bow shock shape does not do a very good job of matching the bow shocks that I am getting in my simulations ehansen12112013
  • checked the output from my 3D pulsed jet run. The first things I checked were pressure and temperature, and the pressure does something strange at the injection boundary.

movie

I am doing some runs with different parameters to see if I can find or fix the issue.

First cut of AstroBEAR GPU

http://www.pas.rochester.edu/~shuleli/gpupro/magneticdiffusion_0019.png

Diffusion test movie here:
http://www.pas.rochester.edu/~shuleli/gpupro/magneticdiffusion.gif

Features:
User can specify blocks and threads setup in a data file. In the cuda code, we use simple setup as:
dim3 grid(XBLOCK,YBLOCK,ZBLOCK);
dim3 block(mx/XBLOCK,my/YBLOCK,mz/ZBLOCK);

Can use two different setups: one can use cudaMalloc to compile 3D data into an 1D array, and does calculation on the 1D array (more pressure on global memory since its poor data locality) or use cudaMalloc3D (not working yet, I think I did something wrong when dereferencing PitchPtrs). They use different kernels. It would be ideal to work out the Malloc3D setup.
Launching two different kernels:
JKernel<<<grid,block>>> to calculate edge-centered current
auxKernel<<<grid,block>>> to convert the currents into aux fields and fluxes

Launching multiple kernels seems to have a small overhead (3~5 cycles) unless new data is copied, but it would be nice to explore interblock communication options since we do need all the currents being updated before updating the aux fields. I'm not sure if there is hardware supported barrier on Tesla, but one can easily use an atomic+ on a global memory address to implement a rudimentary barrier between blocks (no sense-reversal needed here since when returning from kernel launch, all blocks will be synchronized and metadata destroyed).

Real time image processing: I have a code written (with a C# front end though) that does real time processing and image rendering on screen on the GPU using openGL. It would be nice to hook this up to the GPU kernel and render simulation results in real time (ideally can have a webapp that take simple user inputs and display simulation results in real time). For the above diffusion test, each frame takes 10ms to compute, 10ms to process (total time though is not 20ms since I spawn two threads one does the computation, the other only gets interrupted when new data available). It could be greater than 16ms if hydrodynamics is involved so some fine tuning is necessary.

This week: will be working on finishing the TSF paper this week, and some job applications. After the paper is out, I may go back to explore the GPU stuff a bit more.

Meeting update 12/16/13 - Andy

Unfortunately, due to long queue times I was not able to run until this morning. However, Baowei's new planet_profile.dat file solved the problem, but a more elegant solution will be achieved in the future. For now I just need it to run.

This is what the first frame looks like, I just have the density profile of the planet plotted:

The simulation box size only surrounds the planet. When I increased the box size, I see the star form as well but the planet is no longer formed. I believe that the planet is too small to be resolved. At this point I need to investigate the code more and get a better understanding of what is going on.

Meeting Update 12/16/2013 -- Baowei

  • Tickets
    1. new: none
    2. closed: #328
  • Resources
    1. Skype premium account: need to set up a recursive payment and may need Adam's credit card to do that

Meeting update

  • wiki page updates - updated standard out page, and the job sizes page, although the job sizes page needs a bit more work. Question - how many aux variables in 2.5 D and 2.5 D with angular momentum?
  • also want to update the debugging page to include stack memory issues (definition of stack and the quiet seg fault you see with this type of limit, and also how fortran can run out of stack fast). Is default stack 10 mb on most clusters?, and this is per node right? Want to also write an explanation of the size of astrobear alone, and how you can do a baseline test to see what you are actually running. Also want to document that the size of info allocations in the standard out only includes the actual info structure, not any other large arrays such a post-processing might use. also want to update the ticket for CF problem.
  • I read a ton this past week and need to update my library with papers and such
  • I am working on making wiki page for CFs meeting tomorrow
  • I emailed christina about using my code, no word back from her yet, emailed her again today to check up on her
  • Should I register for hotel accommodations for Santa Barbara conference
  • Simulation died with a seg-fault when running 5 levels, but not 4 levels. The TI instability seems to be resolution dependant,…?
  • I am seeing more memory usage than expected on some simulations, going to make a note of on a ticket and come back to later
  • A good reference for MHD numerical methods?

Vertical Magnetic Fields in the CND - Marvin

I simulated the CND with a purely vertical initial magnetic field with a field strength of 0.01 mG. The first animation shows a face-on view of the disks magnetic energy density with streamlines and shows that the vertical field is converted quickly into a torodial field. In the corresponding animation of the edge-on view we see that the vertical Field completely vanishes after a short period of time. The third animation shows the magnetic streamlines within a quarter of the disk: The initially vertical field gets dragged around and transforms into a torodial field.

Interestingly the late stage of these simulations resembles the simulations that I posted last week, although they start with a completely different initial field configuration. Even the magnetic energy density inside the disk reaches values of approx 10-6 G2 (corresponding to the observed values) although the simulation starts with a uniform energy density of 10-10 G2.

Animation of face-on view of the magn. energy density

Animation of edge-on view of the magn. energy density

Animation of the magnetic stream lines in a quarter of the disk

Magnetic Fields (face-on):

Magnetic Fields (edge-on):

Magnetic Streamlines in one quarter of the disk:

Sims updates and memory finding on BH.log

My 643 + 5 was running okay, until I got some nans, then the memory exceeded that requested, and bluehive.log gave me an odd error log (check attached). run:/scratch/ekamins2/CollidingFlows/TEST_2013NOV21/TestStackScript/5levels_nans I can restart from the frame before the NaNs, and possibly use a debug version of the code to try for helpful messages in debugging this?

I am also re-running with 4 levels to see how that run does. run:/scratch/ekamins2/CollidingFlows/TEST_2013NOV21/TestStackScript/4levels

At least with 2 levels (completed in ~10 hrs), I do not think I saw these NaNs arising.. run:/scratch/ekamins2/CollidingFlows/TEST_2013NOV21/TestStackScript and files:/grassdata/erica/from_bluehive/CollidingFlows/3DMHDUnlStack

Also found that the memory used as reported in Bluehive.log ~1-2 GB higher than standard out is reporting for some low res test jobs, and that this memory report is different depending on whether I run on 1 or 2 nodes. Perhaps this mismatch between standard out and bluehive.log is due to some overhead in the inter-node, inter-core communication? run:/scratch/ekamins2/CollidingFlows/TEST_2013NOV21/TestBluehivelog

For the high res job (643 + 5), the mem used was reported as 80 GB > 40 GB total I requested, and also > 50 GB standard out is reporting I used. This indicates that perhaps there was an overhead of about 30 GBs… ? That is, maybe the simulation was using 50 GB, which already exceeded the memory request, but with the overhead, it was actually using 80 GB… and finally was killed by the queuing system

This makes it slightly unclear as to how much memory we should be requesting on the cluster if the memory standard out reports /= memory used by cluster..

bow shock shape

Now that I did some simulations just looking at the shape and interactions of bow shocks from clumps, I wanted to compare with the formula that I had been using for my mach stem simulations.

It looks like the equation for the shape of the bow shock given by Ostriker (paper) does not work well for different gammas. There is not a big difference in shape between = 5/3 and = 1.01.

From my 2 clumps simulation with separation 40 rclump, you can see that the bow shocks are very different.

Furthermore, the equation does not even come close for my gamma = 5/3 run. I think the Ostriker formula assumes a radiative bow shock, which makes more sense for its shape.

The equation comes closer for gamma = 1.01, but there is still a big discrepancy. Maybe it is accurate for a gamma = 1, but I doubt it comes close enough for the purposes that I need.

This leads me to think that perhaps I have to derive my own formula for bow shock shape from my simulations in order to do a proper mach stem study with astrobear.

Update: Bz component of magnetized disk

As you suggested I plotted the Bz component of the disk's magnetic field in units of 1mG (this is the initial field strength of the toroidal field). In the ambient medium the Bz component can take up to 1% of the initial field, these are probably just random variations/noise. Along the z-axis of the disk the Bz component takes up to 10% of the initial field strength

Animation of the Bz component

Meeting Update 12/09/2013 - Eddie

  • Finished all 12 of the 2 clumps runs. Just working on visualizing right now. The rest of the images and movies will be up later today. ehansen12052013
  • The 3D pulsed jet run finally made it through the queue on bluestreak. There was something wrong with the system, so they gave me a reservation to make it run. It still produced a lot of restarts though, so I will have to debug that this week.
  • Also on my agenda for this week is to make some progress with the non-equilibrium cooling.

Meeting update

  • Making the stack unlimited fixed the seg faults for the CF runs! ! :)
  • Have complete now a 3D, 643 + 2 levels run
  • Am running a 643 + 5 levels run now on AFRANK queue
  • Completed latest draft for BE paper

Updates to this link

Updates to http://bearclaw.pas.rochester.edu/trac/astrobear/blog/martinhe11182013

Also working with A Ciardi on synthetic images.

Meeting Update 12/09/2013 -- Baowei

  • Tickets
    1. new: #328(Seg fault in Planetary Atmospheres module (AstroBEAR 3.0))
    2. closed: none
  • Resources
    1. working with Frank on the Skype account
    2. Asked Dave to check the Laptop and the wireless worked on 4th floor. Will show him again if still not working on 3rd floor.
  • New Users
    1. Andy from Rice (modeling of proposed magnetized shock experiments at LLE)
    2. from Western Kentucky University (modeling plasma jets of blazars)
  • Worked on
    1. #309: got same result with both Betti's data and my data. checking on the BCs:

http://astrobear.pas.rochester.edu/trac/astrobear/ticket/309#comment:37, http://astrobear.pas.rochester.edu/trac/astrobear/ticket/309#comment:39, http://astrobear.pas.rochester.edu/trac/astrobear/ticket/309#comment:41

  1. new youtube movies from Martin: http://www.youtube.com/user/URAstroBEAR

Magnetized disk without outflow object

I simulated the circumnuclear disk without outflow object and extended the simulation time from 105 to 106 yr (corresponds to 3 orbital timescales at the accretion disks outer rim (4 pc) and to 100 orbital timescales at the "inner rim" 0.4 pc).

The first animation shows the 3D mass density. After ca. 40% of the simulation time the disk reaches a more or less settled state (i.e. its behaviour is not dominated by the initial conditions anymore). So I think 106 yr is an adequate simulation time.

The second animation shows egde-on and face-on views of the magnetic energy density with streamlines. We see that the ambient medium is in a quite turbulent state, with field strengths that vary by two orders of magnitude and magnetic loops that move away from the disk.

It seems that even without outflow object the magnetic field at the z-axis is converted into a vertical field. This can be seen best in the third Animation, where I plot the normalized z-component of the magnetic field.

Animation of the 3D mass density

Animation of the magnetic energy density (face-on and edge-on)

Animation of the normalized Bz-component

Mass density:

Magnetic Fields (face-on and edge-on):

Magnetic Fields (normalized Bz-component):

new simulations of bow shocks from 2 clumps

Parameters

Now that I have a clear plan for these simulations, I thought I would give a little more detail about the problem set up. These are all 2D, so the clumps are actually slices of cylinders. For simplicity, I am still referring to them as clumps.

For all simulations, these parameters are the same:

Tamb 1000 K
namb 1000 cm-3
Twind 1000 K
nwind 1000 cm-3
vwind 50 km/s
nclump 107 cm-3
rclump 10 AU
Ly 40 rclump
my 50 cells
AMR levels 4
tfinal 329 yrs

The wind is supersonic with a mach number of approximately 13.5. Based on the wind speed and length of the domain, this run time is about 8.67 wind crossing times. The clump temperature is set to be in pressure equilibrium with the ambient.

Simulations

I have a set of 12 runs that I am working on. They differ by clump separation distance (d) and . The length of the y domain is fixed, but the length of the x domain changes. The resolution stays fixed at 20 cells/rclump. Here is a list of separation distances and grid size for each type of run:

d [rclump] Lx [rclump] mx [cells]
40 48 60
20 28 35
15 24 30
10 20 25
5 12 15

The run with d = 40 rclump was added to see a set up where the bow shocks did not interact much if at all. Each set up will be run with a = 5/3 and a = 1.01 for a total of 10 runs. These runs all have open boundary conditions. The remaining 2 runs will be with d = 20 and 15, = 5/3, and they will have a reflecting boundary in the middle instead of simulating both clumps. We want to make sure that the reflecting boundary simulations behave just like the 2 clump simulations.

I will post images and movies here over the next few days…

Images and Movies

movie

movie

movie

movie

movie

movie

bow shocks from 2-clumps simulations

I ran four simulations with 2 clumps and a wind. Each run has a different clump separation: 20, 15, 10, and 5 clump radii.

Much of the setup is the same as what I used for the mach stem runs with some slight changes to make the simulations run faster: rclump = 10 AU, nclump = 106 cm-3, namb = nwind = 1000 cm-3, Tamb = Twind = 1000 K, and vwind = 50 km/s. This means the wind is supersonic with a mach number of approximately 13.5.

movie

Once you get down to separations of about 5 rclump, the bow shocks combine into one large bow shock. This is why my mach stem runs have not behaved as expected. The separation distances I have been testing are typically on the order of 2 rclump.

I have a couple of ideas on what to do next:

  1. The formula that I am using for the bow shock shape could be incorrect or not applicable in this context. I could develop my own formula from simulations.
  2. It might just be impossible to form mach stems by forming bow shocks right next to each other in this way. A possible way around this is to form the bow shocks away from each other and then bring them closer together.

Meeting update

I am trying to debug a seg-fault on BH that is occuring when I run a 3D, AMR simulation of the MHD colliding flows module in the debug queue,

http://astrobear.pas.rochester.edu/trac/astrobear/ticket/326

I had a bit of trouble running with Valgrind, but figured out why (will add to the wiki page on valgrind). First, needed to load the valgrind module on bluehive, and 2nd, needed to modify my .pbs script from this call,

mpirun -np $NUMNODES valgrind -v --leak-check=full --leak-resolution=high --num-callers=40 --read-var-info=yes --track-origins=yes --log-file="valgrind_dump.txt" -machinefile nodefile.$PBS_JOBID ./$PROG > astrobear.log

to this call,

mpirun -np $NUMNODES valgrind -v --leak-check=full --leak-resolution=high --num-callers=40 --read-var-info=yes --track-origins=yes --log-file="valgrind_dump.txt" ./$PROG

Note, I removed the additional log files (machine log and astrobear log). I would like to have those printed as well, and hope by fidgeting with the command can get it to do so. Now, though, I am running with Valgrind, and will add the output to the ticket shortly.

Another problem I ran into when trying to debug the seg fault was that I was unable to get error logs when running with debug flags. Christina was also not getting traceback information. Jonathan said he suspected this might be due to the stack limit on Bluehive, and that I might

-create a launch.s script in your run directory that looks like this:

#!/bin/bash ulimit -s unlimited ./astrobear

-then make it executable by running

chmod u+x launch.s

-then modify your pbs script to run launch.s instead of astrobear - ie,

mpirun -n 64 launch.s

I am also trying this. Will update by meeting hopefully.

I am in communication with a CIRC tech, and am seeing if he can help me fix a bug in bluehive.log. When you run a sim on 1 node, there are 2 helpful variables that are reported: 1) total memory requested, and 2) memory used. However, when you run on multiple nodes, the 2nd value is reported as N/A. It would be nice to have both of these values however as an additional measure of memory usage. I also learned that you can only request 2 nodes in the debug queue, even though the CIRC wiki states you can request up to 16 processors in the queue. This is interesting given there are more than 2 nodes available for the debug queue.

Something on sink particle

This is an experiment. A drop of mercury is vibrating in water. It resemble the fractal of accreting sink particle very much.

Actually they have physical resemblance, too. In this experiment, the vibration of mercury and water are coupled. In our simulation, high density part and low density part are also coupled. Waves serve as the driving force in the experiment.

vibratingmercuryinwater.flv(you need to download it to watch)

It may imply that our simulation is physically right!

Meeting Update 1203

Been writing the triggering paper mostly, here's a directory for paper figures. Have been making the 3D rendered movie for the rotating disk cases, they are not finished yet.

http://www.pas.rochester.edu/~shuleli/tsf_paper/paper_figs/

Resubmitted the NLUF jobs with Jonathan's extrapolation routine to bluehive. Should get some results this week if things go on OK. This is using the old design parameters Peter sent around a while ago (setup another conference call with the NLUF team?).

Meeting Update 12/03/2013 - Eddie

  • Marvin is starting to use NEQ cooling in his disk module. I will help with any issues that may come up.
  • I need to remake the synthetic emission maps for Andrea. The previous ones show the "bubble", but what he really wants to be able to see is the jet. If there is still emission in the jet then this should just be a matter of scaling differently in visit.

Mach Stems

I have been trying different 2D runs for the mach stems. The first attempt looked the same as the slice from the 3D runs. I thought that the behavior we were seeing might be caused by the wind, so I moved the wind boundary farther away. Here is the result,

movie

The situation here is not as clear cut as a wind of uniform density sweeping over a clump and making a bow shock. I think the fact that the wind creates a forward and reverse shock is changing the behavior from what we might expect.

I then tried setting the ambient and the wind to be the same density. I also ran the simulation for twice as long.

movie

One mighty bow shock is formed, and the right reflective boundary (the one we are supposedly not interested in) seems to be playing a role as well.

I want to turn off the reflective boundaries and just look at the shape of the bow shock that is formed. I could also just keep running the simulation to see if I ever get to a steady state, but I still don't really like all the lateral movement going on. Any other ideas?

Meeting Update 12/03/2013 -- Baowei

  • Tickets
    1. new: #325(Investigate Grapevine Load Balancing), #326(seg fault on BH with colliding flows), #327(time variable)
    2. closed: #327

Magnetized Cirumnuclear Disk

This is the first simulation of the Circumnuclear Disk with magnetic fields. It includes an initial torodial field configuration with a constant field strength of 1 mG. The first animation shows the mass density, it looks more or less the same than the simulations without mgnetic fields. The next two animations show edge-on and face-on views of the magnetic energy density with fieldlines in black. The edge-on view shows a number of magnetic loops that are created and carried outwards, furthermore the central outflow seems to change the magnetic field in the center from a torodial to a vertical configuration. In the face-on animation the magnetic field maintains more or less its torodial configuration, but later the magnetic field seems to vanish from the inner region.

Animation of the 3D mass density

Animation of the magnetic energy density (edge-on)

Animation of the magnetic energy density (face-on)

Mass density:

Magnetic Fields (edge-on):

Magnetic Fields (face-on):

Meeting update 12/3/13 - Andy

My overall goal is to place a Jupiter-like planet around a star which I will eventually put so close to the star so that it gets tidally shredded. For a long time I attempted to use a BE sphere to do this but it seems that the BE sphere gets ram pressure stripped because the density contrast at the boundary is not high enough to withstand the velocities at the shredding radius. Jason suggested that I try a polytrope object or the density profile that Jonathan and collaborators were working on. I decided to go with the density profile which Jonathan supplied me with and had to get the current version of AstroBEAR(3.0) to get it working. With Baowei's help, I was able to compile 3.0. It seems that lines 69 and 70 are the same in the makefile which is for the physics/EOS.f90 file. This created an error with my compiler which was easily fixed by deleting one of the lines. After changing that and a lot of emails with Jonathan I was able to get his Planetary Atmospheres module to compile successfully.

At this point I keep seg faulting when I go to run the code. After running with the error flags on I got the following error message:

 forrtl: severe (408): fort: (2): Subscript #1 of the array DATA has value 1 which is greater than the upper bound of 0

This occurs at line 161 of the problem.f90 file which is where the code tries to assign r_cutoff. It does this by looping through the density profile and finding the point where the density is lower than the ambient and then assigning r_cutoff to the be that radius.

 DO i=1,nEntries-1
   IF (Clump%profile%data(i,2) < rho_amb) EXIT
 END DO
 r_cutoff=Clump%profile%data(i,1)

When I looked at the density profile it seems that the density is never lower than the ambient. I believe this is the source of the problem but when I changed the density profile for the last density entry to be lower than the ambient it still seg faults at the same point. Although, if the density is never lower than the ambient I think it should just assign r_cutoff with the radius at i=nEntries-1…

LineTransfer

So it's not much - but line transfer works for Fixed Grid, Single Processor, +X direction only

And now it works in parallel (here it is running on 13 processors)

And now it works in AMR and parallel (13 processors w/ 2 levels)

movie