Posts for the month of October 2011

Meeting Update 11/01/2011 - Eddie

  • I finished my first hydro code last week. Just as a reminder, it is one dimensional, uses an exact Riemann solver, and uses a first order Godunov upwind scheme for updating.
  • I ran the same tests that Toro runs, and fixed a couple of errors. All the plots for these tests can be found on my build code page.
  • Also, this page is getting more and more detail on it every week. There's still quite a bit of information that I want to put on it that I've already learned. Just haven't gotten around to it yet.
  • The next step is to put some approximate Riemann solvers into my code. I've already begun reading about them, and will probably have this step done by the end of the week. Should I systematically go through Toro and put all approximate solvers in my code, or are there just a couple specific ones that I should use?

Meeting update -- Rotation

The product of initial central angular velocity with the free fall time, as measured from the central density of the clump — omega*tff, can be used as a parameter for the rotation problem. Authors Matsumoto & Hanawa ('03) have found that when this product is low enough to allow collapse, but larger than ~ 0.05 (when a disk can form before first core formation), the clump fragments into a particular morphology. When the BE sphere is critical, with a small density "enhancement", SLOW rotation induces collapse.

These authors initialize angular velocity with the following:

This perturbs the amplitude of the global rotation,

with varying amplitudes of Omega2 and Omega3, the m=2 and m=3 velocity perturbations. The parameter C determines the angular velocity's dependance on r. For C=0, rotation is rigid-body, for C larger, Omega decreases more rapidly with increasing r. Both C and beta, the ratio of rotational to gravitational energy strongly determine the resulting morphology of collapse — the fragmentation into bars, discs, etc.

The authors Banerjee, Pudritz, Holmes on the other hand study the collapse of both high mass and low mass BE spheres. There initial conditions prescribe a rigid body rotation that is also determined by the product Omega*tff. They use this parameter = 0.1, 0.2, 0.3. They go through the calculation of both alpha (thermal/gravitational energy) and beta (rotational/gravitational energy) of BE spheres and explain that beta is maximum at the OUTER edge of the cloud.

My questions are —

1) Rigid body vs. keplerian rotation? What conditions for each rotation? If we start out with rigid body, do we expect the collapse to evolve to a keplerian rotation?

2) Why does slow rotation trigger collapse rather than expansion or nothing? How does alpha and beta change with r and contribute to the total stability of the clump? How does the speed of rotation effect stability?

3) A good way to visualize the simulation output?

4) Initial conditions for the BE sphere — should the set-up be one that is JUST critical, but doesn't collapse? Or should it be a set-up that is already super-critical (thermal energy reduced in grid by 10%)?


Current Run:

Is in the queue in BH

It has:

-An omega*tff = 0.1 (rigid body rotation)

-A critical BE sphere (stable at the critical radius)

-Grid is shifted in x,y by ½ dx

-43+1level cells — ~ 20 zones p/clump radius

Created a page about dissipation methods in the hydro solver

Hydro code finished!

I have finished building my first 1D Hydro code. It uses an exact Riemann solver, and a first order Godunov upwind scheme for updating. At first it was not giving the correct solution, and after staring at it for a while Jonathan found a typo in one line in which the wrong gamma constant was being used. Now it seems to work great. I will run the 5 different tests in Toro, and make plots to show in next week's meeting. Also, I'm working on a separate program to compute the exact solution to directly compare it to the numerical one.

I have definitely learned a lot through this process, and my understanding of the numerical methods is much clearer now. I'm currently working on my buildcode page which documents the important aspects of a 1D Hydro code.

Results of recent sink tests

A shift in the xy plane by 1/10th dx DID NOT form a sink (not in AMR, not in fixed). Evidently, in order to form a sink the potential minima in the grid must be well defined — slightly shifting the grid by 1/10 dx is NOT large enough to force the minima to be in a single center cell.

In contrast, shifting the grid by ½ dx in x and y successfully restricted the minima to reside in a single cell, at the center of the grid (in x and y). Though, it was along a face in the z plane. Since the potential was symmetric along x and y, a sink was able to form. However, since the potential was asymmetric along z, a kick happened.

The above tests were done using an even number of cells — 86 in each direction. They were prompted by the fact that a cell formed a sink with 85 cubed fixed, but not when using AMR with a similar eff. resolution. This indicated that sink formation rested on whether the potential was centered in the middle of the grid, located in a single cell, as it was in 853.

Meeting Update 10/26/2011 - Eddie

Following the examples in Toro, I finished writing an exact Riemann solver and a 1D code that solves the inviscid Burger's equation. It uses a 1st order Godunov upwind scheme. I'm still trying to test it and also write a code that can solve more general problems, but I think it would be good to start looking at approximate Riemann solvers now.

Meeting Update 10.25

The following plots show the AstroBEAR test of various settings in the current sheet outflow problem.

http://www.pas.rochester.edu/~shuleli/CS_outflow.png

The Workman setting has central symmetry while the reconnection sites locate at the two edges. The positive radius of curvature of field lines at the reconnection site results in an inflow to the X point.
The first panel shows its kinetic energy distribution.
The Semenov setting is achieved by applying a resistivity function that varies with position. This resistivity gradient results in bending of magnetic field lines with negative radius of curvature. A vertical outflow from the X point is thus created. The resulted Sweet-Parker region is expanding throughout the simulation till the flow speed at its boundary reaches the Alfven speed. Panels 2 and 3 shows the kinetic energy and vertical flow momentum. Panel 4 shows the kinetic energy of Semenov outflow with an enhanced resistivity contrary.

10/25 Update

Got Shape working/figured out so we can get emission maps of CRL618 done. Mostly have been playing with the Shape interface.

Currently working on making maps of t=50yr,100yr of the 4 runs, and will post them as they are done.

Once these are done, I will need to do maps of

  1. Emission-measure images at a viewing angles of ~60deg and 90deg from the symmetry axis. EM = integral(n2 dl).
  2. Position-Doppler-shift diagrams for each of these four frames

First EM result, 100 years, 90degs: EM on left, PV on right

Clump with Rings
http://www.pas.rochester.edu/~blin/crl618/clumpring/90degclumpring.jpg http://www.pas.rochester.edu/~blin/crl618/clumpring/90degpvclumpring.jpg

Jet with Rings
http://www.pas.rochester.edu/~blin/crl618/jetring/90degjetrings.jpghttp://www.pas.rochester.edu/~blin/crl618/jetring/90degpvjetrings.jpg

Jet without Rings
http://www.pas.rochester.edu/~blin/crl618/jetnoring/90degjetring.jpghttp://www.pas.rochester.edu/~blin/crl618/jetnoring/90degpvjetnr.jpg

Clump without rings will be done as soon as the rerun is done.

Reconstruction Limiters...

Running some 2D colliding flow stuff and noticed that gudonov and ppm give very different results initially. Gudonov method looks diffusive, but correct while PPM looks funny. Not sure if this is just unresolved cooling or a potential bug in PPM or just the way the multi-D limiters are being used…

1st order reconstruction 2nd order 3rd order

I double checked that the ppm and plm methods revert to gudonov type methods when the additional shock flatteners are triggered and found one minor 'bug'. I also compared the results with no limiters at all and they are worse.

Here they are with the corrected form of the current limiters

And here is a map showing the actual value of the slope limiters used

The horizontal striping is a result of the stencil used for calculating the multidimensional limiter. First limiters are calculated using slopes in x and y. Then the multidimensional limiter for a given cell is computed as the minimum of the x and y slopes in the surrounding cells as below

y
x xy x
y

To reduce the striping I modified the multidimensional limiter to instead be the minimum of

xy xy xy
xy xy xy
xy xy xy

Here is the resulting plot of the multidimensional limiter

Finally I also added an extra dissipation that is usually used with PPM schemes that is proportional to the divergence (and is active for compressible flow).

(Left panel has the additional dissipation)

And here is a movie

And finally I did a convergence study at 5 resolutions (left to right) of the 1st order reconstruction (bottom), ppm with diffusion (middle), and ppm without diffusion (top)

I also added some smoothing to the velocity interface to ease the shock in so that the initial velocity gradient was not a function of the resolution.

Note that the lower res PPM schemes are actually better able to resolve some features then the next higher res gudonov scheme (albeit more pixelated).

Here is the movie

To see the modifications in action with mixing see the Colliding Flows Page

Hierarchical Gravitational Collapse

So I've started playing around with gravo-turbulence in 2D. I've tried using linear perturbations to the density within a clump, however the entire clump will collapse in the same time it takes for the local perturbations to grow. So we need a clump with small Jeans-unstable clumplets inside that have a much shorter free fall time. So I tried a density distribution where log(rho) had a flat power spectra. .

The clump density was 400 and the temperature was 10K giving an overall Jeans length of 1.1 pc and free fall time of 1.7 Myr.

Since the perturbations kept the pressure constant, the Temperature was proportional to . As a result, the free-fall time scales like and the Jeans length goes like

The perturbations had a flat spectra from k_min=4 to kmax=48 or from 6 pc to .5 pc. As a result, the local density enhancements had typical size scales of .5 pc (8 cells). Most of these had densities of at least 800 (twice the mean) with several having densities of 1600 and a few with densities of 3200. The ones with a density of 800 would be stable to collapse and those with a density of 1600 have Jeans lengths of .26 pc and could potentially collapse - although the timescales would be similar to the timescales for the global collapse. This is in fact what we see… Two particles form before the entire thing collapses.

So i modified the mean density to be 10, increased the power spectra slope to -2 and strengthened the perturbation. (Plot is now in log scale)

And here is the movie

Everything looks great until the particles form - then the whole thing expands. After some pondering - I realized that this was because the potential of the sink particles needs to be modified in 2D to represent the potential of sink cylinders. The potential goes like and the gravitational force goes like .

Meeting Update Week 10.18

Updated the Magnetized Cloud Project page.
The parallel AMR works now after fixing the convergence issue. Passed the interface test. The Harris sheet test requires more than 8 processors so I am planning to do it on bh or bg.
Also there is a table on the project page listing the runs needed for magnetized cloud simulation. It seems we need to reduce the number of runs: there are totally 45 runs.

Meeting 10/18

Main achievements: -Better understanding of the code -Got a non-zero set of chombo files!!

Main Issues: -Do the data make sense? -Need to talk to Martin about specific parts I still don't understand (e.g. What values pointgravityobj takes care of, why ang. momentum was included in Travis's module, rs/cs values)

Progress:

Last week I managed to get the module to compile properly, and thought everything was more or less in order. Unfortunately, I kept running into "Nans found in ProblemBeforeStep restart requested". I started fiddling around with the mod ules thinking it may have to do with the conditional on lRestart. However, I learned at the testing meeting that it is an entirely different function/feature. Up until then, I had simply been editing Travis's previous module on the Bondi problem, but over the weekend I decided to write it from (almost) scratch, using Travis's module and Eddie's RT Instability module as references. I was careful not to copy/paste anything unless I at least understood what it was doing. I was finally able to get it working earlier today, meaning it compiled, ran without problems and gave me nonzero chombo files. (Hooray!)

These pictures seem like they look right:

http://i.imgur.com/xEomj.jpg http://i.imgur.com/Wqbc2.jpg

Questions/Issues for this week:

These images seem less right:

http://i.imgur.com/UXqxe.jpg http://i.imgur.com/yAZOD.jpg

Basically, the simulation starts out doing exactly what I think it will, then it starts looking kind of funny. Goals for this week then is to diagnose/fix this, and maybe then running postprocessing (bear2fix?). Also, I need to get my webspace set up for at least hosting the images (had to upload to imgur for this week) and I need to figure out how to make those nice .gifs in Visit.

Update 10/18

Radiative instability has been added to the repository. I've been mainly trying to make sure that the code works properly and getting things to compile.

Meeting Update 10/18/2011 - Eddie

  • put together files for restart test
  • read Toro Ch.6 - The Method of Godunov for Non-linear Systems
  • solved Riemann problem analytically
  • looking at exact solver in astrobear, will start coding tomorrow

Meeting Update

I have largely been focusing on course work this week as I am getting ready for midterms. For the group — sink runs are in bluehive for shifting the grid in y and z by 1/10 * dx to check if the sink still moves. I'll post the results when they are up.

Problem Submission Script

I created a problem submission script based in part on the old buildtests.s script. I have debugged it enough to the point where it seems usable, please let me know if there are any issues with it.

The script is called buildproblem and is located at /home/noyesma/buildproblem on grass.

Pre-Conditions

  1. The astrobear code has been checked out to some location
  2. The module directory exists and contains properly formatted *.data files.
  3. A Makefile.inc.[hostname] file exists in the root code directory with [hostname] being the name of the machine being compiled for

Features

  1. Compiles AstroBEAR based on the machine being run on
  2. Populates TABLES and out directories
  3. Copies and runs astrobear to the run directory (currently the run directory is the module directory)
  4. Outputs run log to a (date).log file

Usage

buildproblem [-d code_dir] [-m module] [-np processors] [-t]
buildproblem -d /grassdata/noyesma/astrobear -m RTInstability -np 2
cd /grassdata/noyesma/astrobear
buildproblem -t -np 8
Command line Argument Description Default Value
-d Code Root Directory Current Directory
-m Module to Compile None; will prompt if not provided
-np Number of Processors to Run On 0 (i.e. will not run astrobear. Provide a value > 0 to run astrobear)
-t Specifies testing mode. Will compile all well-formed modules in the code/modules/tests directory Off
-u Upload to clover:/data Off

Tips and Tricks

1) Running

buildproblem [arguments] &> /dev/null &

will compile and run astrobear in the background and pipe all non-error output to null.

Velocity/Density Perturbations and Sink Particles in 2D

I added density and velocity perturbations to the Ambient object and wrote routines for initializing density and velocity spectra with a given range in wavenumbers and slope. Currenty only solenoidal velocity fields are supported. Also ran some 2D gravo-turbulent and TI-turbulent sims and fixed a few bugs with Sink Particles in 2D.

field loop restart test

In addition to the normal tests that we are going to run in the test suite, Martin thought it would be a good idea to have one restart test. The 2D field loop advection problem was chosen for this. I completed all the runs to make sure it worked, copied the necessary files over to the test directory in grass, and updated the wiki page: TestSuite/FieldLoop2D. This test will be run just like all the other tests in the test suite. The only difference is that this test comes with an out directory that already has a chombo, so that astrobear knows when/where to start from.

Meeting Update for 10/11/2011

  • Checked in threaded version of the code along with some minor bug fixes that allowed me to close tickets
  • Posted the scaling results and closed ticket #95 (Scrambler Scaling Tests)
  • Closed ticket #148 as the errors were related to stack issues on Kure.
  • I also closed ticket #98 and instead created the milestone:"Implementing Multi-Physics" as well as individual tickets for various desired multiphysics capabilities
    • #150 (Implementing Self Gravity in Cylindrical Coordinates)
    • #151 (Thermal Diffusion)
    • #152 (Implmenting Magnetic Resistivity)
    • #153 (Implementing Viscosity)
    • #154 (Sink Particles in 2D and Cylindrical)
  • Added testbuilds.s script to repo that compiles every problem module (although it could also be modified to start each problem. This is useful if their are modifications to the way problem modules are interfaced or if stuff is removed from data files…)
    #!/bin/bash
    cd modules
    for i in $(ls -d */ | grep -vE 'objects|tests|Problem') ; do
      rm Problem
      ln -s $i Problem
      echo Testing $i
      touch Problem/problem.f90
      cd ../
      make
      cd modules
     done
    cd ../
    
    for i in $(ls -d tests/*/) ; do
      rm Problem
      ln -s $i Problem
      echo Testing $i
      touch Problem/problem.f90
      cd ../
      make
      cd modules
     done
    cd ../
    
    

Meeting Update 10.10

So I finished the resistive/viscous implementation. It can run the Run Test (tests on the Project page) without problem. Works in parallel and AMR. And it's compatible both in 2 and 3-D.

Only problem is that the original implementation needs to write a whole set of communication routines which seem to need some generalization/organization. So this working version is using an asymmetric solver as a work around, and it can not do subcycling (we may not want it for the next paper anyway).

I think we are at a point to discuss the detailed runs we want to take.

Meeting Update 10/11

Edit: bear2fix error has been fixed. Martin has been making progress in Shape, and should have emission maps within a week or two.

I had trouble logging into the wiki last week, but Rich fixed the issue over the week.

For testing, the Radiative Shocks page is updated and ready. All the test and data files have been copied over to the Martin's testing directory. https://clover.pas.rochester.edu/trac/astrobear/wiki/TestSuite/RadiativeShocks

For CRL 618, http://www.pas.rochester.edu/~blin/crl618/results.html

Bruce noted that clump-with-no-ring has very poor spatial resolution except near the clump itself. This makes a comparison of this model with the others very awkward. We had to lower the resolution significantly due to sizing. The clump chombo sizes are very large, and resolution had to be compensated. Perhaps there is a way to overcome this issue and still have resolved results?

We need:

  • Emission-measure images at a viewing angles of ~60deg and 90deg from the symmetry axis. EM = integral(n2 dl).
  • Position-Doppler-shift diagrams for each of these four frames

However, Martin and I are having trouble trying to do ascii files from 3D hydro chombos using bear2fix.

I am trying restarting runs for Clump-with-Rings since I only ran it to 117 years and we need up to 150 years. However, Bluehive is complaining:

HDF5 error: function h5fopen_f failed with error code -1 .
HDF5 error: function h5fopen_f failed with error code -1 .
HDF5 error: function h5fopen_f failed with error code -1 .

Meeting Update 10/11/2011 - Eddie

Agenda:

  • Test suite preparations done ehansen10062011
  • Read Pat's review article on measurements of physical conditions in jets, sent him an email…hopefully start learning specifics about his cooling table soon
  • Read Toro Ch. 5 - Notions on Numerical Methods, gathering thoughts/ideas for Wednesday's code meeting
  • Galactic Dynamics course project ideas? ehansen10102011b
  • Group organization ideas? ehansen10102011a

Galactic Dynamics Course Project

I have an opportunity to do some research/learning that is relevant to both AstroBEAR and a course that I am taking. Eric Blackman is teaching Galactic Dynamics this semester and I need to do a presentation and write a paper. The topic is supposed to be "controversial" and probably somehow related to galactic dynamics. As an example, Blackman mentioned the Dark Matter vs. MOND debate. I'd like to do something that will help me understand some physical process that is either already in AstroBEAR or maybe something that might go into AstroBEAR in the distant future. My favorite idea so far is to do something on AGN jets. Does anyone have any suggestions for a more specific AGN jet topic? Or any different paper ideas on a different topic related to galactic dynamics?

Group Organization

I think it's important that we clarify what google groups we need and what teams we should have. We also need to discuss who should lead each team, how the team should function, etc. In particular I'm thinking of the Development, Debug, Testing, Documentation, and any other teams. We'll talk more about this in tomorrow's group meeting, but this blog post is intended to generate a discussion so please comment!

Meeting update -- week of 10/10

Test pages updated

I've finished updating some test pages and getting the bear2fix.data files ready. The following test pages have been updated, complete with images from the reference chombo: TestSuite/RayleighTaylor2D,TestSuite/UniformCollapse,TestSuite/FieldLoop2D. All bear2fix.data files necessary to reproduce images from a new chombo, and those necessary for chombo comparison have been copied to Martin's test directory in grass.

The only slightly tricky one was UniformCollapse since it is a 3D problem. I had to make a very minor change in bear2fix that makes a density plot in x and y by averaging over z. So to reproduce the image for that problem, you'll need the current version of bear2fix which I already pushed to the devel repo.

Meeting update

Here is a page of the tests I re-ran and ran on the sink particle formation:

http://www.pas.rochester.edu/~erica/SINKtests.html.

Note, a sink only formed on a fixed grid, of an ODD number of cells. I was unable to form a sink with an even number of cells (and thus also AMR). The complaints were "system unbound".

Also, this sink DID NOT WANDER.

This may be due to it being in the 'true' center of the grid.

In other news, here is the paper by Bannerjee and Pudritz on rotating BE spheres:

http://adsabs.harvard.edu/abs/2004MNRAS.355..248B

First Blog Post 10/04/11

I was assigned the Bondi Accretion problem with which to do my first module. After consulting with Martin extensively I've been reading Bondi's original paper, as well as:

Galactic Black Hole - Falcke & Hehl High Energy Astrophysics - Melia Physics of Astrophysics Vol. 2 - Shu

My initial issue was looking through the equations and finding the simpler solution for particles that have already crossed the sonic barrier for speed, rather than the more complicated full solution to the problem. Eventually, the simple solution of v = sqrt(2GM/r) made sense, as they have just enough energy to get back out to infinity, but since they are being accreted, have a velocity going the opposite direction.

This last week has been a tutorial in transferring physics into code. I managed to compile astrobear for my problem.f90 file, but got many errors while trying to run the program. I spent a good part of the weekend trying to figure out if there was anything in my problem.f90 file that was causing these errors, but eventually realized that it was because all of my .data files were only in the Problem directory, and not the astrobear directory as well. After fixing this, that error went away, but now astrobear complains of an error with my global.data & physics.data.

I feel as if once I have successfully run this simulation, future ones will be a lot simpler to do!

Tweaking the wiki

I modified the old astrobear page to automatically redirect to the wiki. Is this what we want? Or do we want to maintain a separate website?

Also I had Rich add a development team. Previously development tickets were being marked as for the debug team or the parallelization team. I also created two new reports which are useful for getting a handle on open tickets.

Resistive MHD test

So I've decided to use Harris current sheet test as a test problem. It's essentially an instability with easy to identify flow pattern. And it is also a good research problem if we want to dig into it. I will talk about this test in tomorrow's meeting.

Meeting Update 10/03/2011 - Eddie

Activity last week:

  • Read Toro Ch. 4 - The Riemann Problem for the Euler Equations
  • Read Dopita Ch. 3 - Collisional Excitations
  • Created Development Procedure page
  • Updated test suite pages
  • Gathered .data files for RT and UniformCollapse tests