test
test
Spitzer Resistivity
I added spitzer resistivity to AstroBEAR and it does appear to make a significant difference to the formation of the stand-off/stagnation shocks.
For the spitzer resistivity, I am using
And then the magnetic diffusivity (with units of cm^{2}/s) as
HEDP Simulations
Results
There were 4 sequences of runs to compare the effect of
- Cooling vs no cooling on fiducial data
- Effect of magnetic field with cooling
- Effect of changes in magnetic diffusion length (for cooling and 2X fiducial field)
- Effect of changing the cooling length from .1 to 10x for the 2X fiducial field with ½ the fiducial diffusion length
Reference Image |
First comparison is the fiducial setup but with and without cooling
With cooling |
Without cooling |
Next we just used cooling but adjusted the magnetic field from 1.75 T to 14 T in multiples of 2
1.75 T |
3.5 T |
7 T |
14 T |
We then looked at the 7T case and explored the effect of diffusion comparing the case with no diffusion, ½ the fiducial, the fiducial, and 10X the fiducial lengths
0 |
.366 mm (½ fiducial) |
.732 mm (fiducial) |
7.32 mm (10 * fiducial) |
And finally we looked at the effect of changing the cooling in the 7T half diffusion length case
.1 cooling |
½ cooling |
fiducial cooling |
2x cooling |
10x cooling |
Thoughts on time stepping and MHD
Thoughts on Energy Tracking in AstroBEAR
AstroBEAR does not track the thermal, magnetic, and kinetic energy separately (because it is trying to conserve total energy via conservative fluxing). It does however separately track density, momentum, and magnetic fields - which can be used to derive magnetic and kinetic energy - using thermal energy as a reservoir for discretization errors.
When you have flows that are dominated by non-thermal forms of energy (
or ), discretization errors in total energy (while dynamically unimportant) can still lead to significant relative errors in thermal energy - which can be problematic if there are significant temperature-dependent microphysical processes. In those cases it might be better to solve the thermal energy equation independently and not worry about conserving total energy.Thoughts on time-stepping
Simulations typically run for a few dynamical times - and for flows that are kinetic energy dominated (
and ) the computational time is independent of the flow speed (only a function of resolution ). However, for flows that are magnetic or thermally energy dominated - the time stepping is limited by the Alfven or sound speeds respectively - and the computational time goes as (ignoring the extra factor of due to changes in the number of zones with resolution)This can be combined as
This makes simulating the RT instability
relatively computational expensive. Alfven waves can also restrict the time stepping when.
So simulations with modest
but very small will also be relatively computational expensive.High Alfven speeds are also somewhat easier to generate - since the magnetic fields/energy and density are not as tightly coupled as the thermal energy and density - due to material being free to leave along field lines - and the lack of flux freezing when magnetic resistivity is used.
Explicit Magnetic Resistivity
Astrobear currently implements magnetic resistivity explicitly - without subcycling - so time steps are limited by the smaller of the cell diffusion time
and the cell crossing time . The cell diffusion time will be smaller than the cell crossing time when
So the computational time to simulate a crossing time will go as
So - you pay a penalty with explicit time stepping when
that could be avoided with implicit time stepping.Simulations involving advection around and diffusion through an obstacle
Simulations involving advection around and diffusion through an obstacle will need to run for the longer of the crossing time around - and the diffusion time within (which is longer than the crossing time by a factor of
) - so in the explicit case we have
Making the resistive solve implicit would reduce this to
Hall MHD
Added some documentation on MHD in AstroBEAR (including Hall terms)
Changes in climate sensitivity as a function of Temperature
So apparently radiative forcing going like log CO2 has a name … https://en.wikipedia.org/wiki/Svante_Arrhenius#Greenhouse_effect though the explanation is somewhat complicated… See
https://escholarship.mcgill.ca/downloads/db78th05j?locale=en
Also, from Wikipedia, we have
where the climate sensitivity
and we can approximate the logarithmic response of radiative forcing to changes in CO2
Combining these we get
or
or
and
So the climate sensitivity is not constant - but falls off by a factor of
every 4.3K change in temperature.Wind-Clump interactions
Here are lineouts from 4 of the runs with Bz field only. Note
and Also - the time reported should be scaled by 12.3 - so the final time is 426 nsEach picture shows lineouts at three locations shown here
Each plot shows
- ram + thermal pressure
- magnetic pressure
- total pressure
- density
- magnetic field
- velocity (in x)
- Temperature
And the values are scaled as follows
Density | 1e17 cm^{-3 } |
Pressure | Mbar |
B | Tesla |
T | eV |
v | 10 km/s |
Click for movies
No Cooling | Cooling | |
No Diffusion | ||
Diffusion |
Wind-Clump interactions
I ran a series of numerical experiments to get a feel for how things like magnetic diffusion, cooling, field orientation might affect the results of the setup shown here (B-field is parallel to wire obstacles - into the board (z direction)). Cases with Diffusion have a diffusion scale length of .3 cm - and anything that says Tiny Diffusion uses a diffusion length scale of .01 cm. I also ran a few cases where the field is inclined - so By = .1*Bz
Left panel is density more or less scaled to final max density. Right panel is Bz autoscaled Cases with non-zero By field show contours of the vector potential that trace the field.
testing
testing
test
test
abc
abc
Ionization transfer
To solve the ionization equations with scattering, you would want to use a diffusion equation. Something like
Now the equations are non-linear - so one has to iterate towards convergence. I found a good write-up of the method necessary in a paper on Enzo-RT (where they also use FLD for ionizing radiation transport).
Meeting Update
Looking into tighter coupling between hyperbolic fluxes and source terms
https://astrobear.pas.rochester.edu/trac/wiki/u/johannjc/scratchpad24
And ways to improve corotating frame and possibly implement cylindrical geometry
https://astrobear.pas.rochester.edu/trac/blog/johannjc10222013
meeting update
Did some modeling of random particle distributions to understand galactic colonization in a static system.
Here is a plot showing the number of 'isolated' subgroups of systems given a probe range. The number of isolated subgroups (y-axis) is scaled by N and the probe range (x-axis) is scaled by (1/N^{D}) where D is the dimension. This demonstrates the cutoff in propagation for static systems when the density drops to only a few habitable systems within range on average. The two sets of red lines are for 2D and 3D and they were calculated with 1000 and 10,000 systems. (The 2D system has the longer tail)
And here is another plot showing the fraction of the number of systems that can be reached starting at a point chosen at random as a function of the average number of systems within a sphere of radius probe range.
A good read
I updated the developer guide with the various graphs showing the stencil components and their dependencies… See SweepScheme
Meeting update
- considering trajectories in lab frame
- Wire turbulence MHD run is at 157 and running well…
- Working on paper…
New sims for HD209458b 102% done
Here are the sims to frame 50 - though they are completed to frame 100
Neutral hydrogen |
Density |
volume rendering of log(rho) |
Plan is to follow two more orbits with maximally refined around planet out to 7 planet radii
t=10 orbits | t=10.1 orbits | t=10.2 orbits |
Also, new runs fix issues with odd flow patterns near star
link to krumholz paper
Equations 18-22 and surrounding text discuss various opacities…
mean free paths
cross section for hydrogen hydrogen collisions can be estimated as
where is approximately the bohr radius… This gives
The MFP is then
The material leaving the planet is at
so the MFP is then 1e-7 AU << dxLine Emission from HD204598b
Blackbody luminosity
We have previously calculated the blackbody luminosity from the Intensity as
where
which gives
where
and
has units of ergs/s/HzWe also have from Bourrier et al, a theoretical line emission curve for Lyman Alpha
"Theoretical instrinsic stellar emission line profile scaled to the Earth distance"
with units of ergs/s/cm^{2}/A
After multiplying by 4*pi*47.1 pc ^{2} to get Luminosity per angstrom, we need to convert to Luminostiy per Hz
To go from
towe need to multiply by
So after some conversions, I get the following which shows the previously used black body, along with the line emission profile.
Now if I assume that the atmosphere is optically thin, and integrate the emission profile (+ continuum) over the absorption cross section -
where
and
This cutoff frequency is much higher than the Lyman Alpha wings… Which is not terribly surprising… Lyman Alpha ( even thermally broadened) will not ionize any hydrogen… So, I don't need to change the Ionization Routines (except to adjust the temperature if necessary etc…)
And creating synthetic emission lines - is just a matter of multiplying the background brightness at each wavelength by the emission profile…
The first row shows from left to right the data extracted from the plot in Bourrier et al and the spline fit, the same plot with the x-axis in Angstroms instead of km/s, and the Luminosity per frequency from the line emission as well as the black body.
The second row shows the brightness as a function of time and frequency assuming the Line emission profile as well as lineouts at various times and frequencies. (Time is ordered correctly in all of these plots)
Synthetic Absoprtion profiles for HD209458b
Normalized transmission as a function of velocity for gas in front of the stellar surface. The y-axis shows the time and there is a faint band from t=5 to t=7 hours that correspond to the time when the planet itself is in front of the star. The velocity is the projection of the gas velocity towards the camera - so positive velocities should correspond to shorter wavelengths and higher frequencies. The shifting to higher velocities at early times I believe can be attributed to the gas emitted from the day side that gets swept out in front of the planet and then gets accelerated outward by the stellar wind.
(Note the earlier version of this plot had the velocity scale incorrectly going from -300 to 300 km/s)
And here is a movie showing the emission at the center wavelength
Prescription for making line projections from Planetary Sims
The prescription of Bourrier et al is to calculate (at each pixel) the average optical depth for range of velocity bins with spacing (or frequency bins with size , so that
So we have
and
Now to get the line profile in terms of velocities, we can start with the line profile in terms of frequencies:
and realizing that this is a distribution, we can use
where
so
and
to get
Putting this all together we have…
Now if we bin the column densities into velocity bins, we have
and
and we have
or
where
Now
is of order 20 km/s and is 76 m/s… so if , then the integral may as well go to infinity and we have
and for
(wave arms in air for a bit)
so we have
Or at least that's how I get agreement with their equation 11…
So were I to implement this approach, I would create projections (integrations along LOS) of
where
where Saha equation - or balancing radiative recombination with direct ionization… They give very different answers for the ionization fraction at the planet wind temp.
corresponds to the bin that the LOS cell velocity lies within and is the ionization fraction computing using either theRemaining questions
- How to chose the ionization fraction? - Need to balance photo ionization with radiative recombination
- We should use thermal broadening (which dominates over natural broadening)..
- Binning particle velocities before convolution seems numerically expedient - but binning errors are orders larger than convolution corrections…?
- Why not just numerically approximate integral of line profile - and calculate contributions to each absorption bin for each cell.
Thermal broadening
It should be more accurate to calculate the LOS velocity and the thermally broadened line profile, and then calculate the contribution of that profile to each bin in frequency space. The thermal broadening profile is just a Gaussian, so using the error function we can calculate the contribution to each bin…
Now without thermal broadening, we have
however, with thermal broadening, we have
where
which gives us
Now the integral over
is a convolution with a lorentz and gaussian profile - which gives a Voight profile… however, the thermal broadening is much larger than the natural broadening, so we can approximate the natural broadening as a delta function…
which gives
or upon integrating the exponential, we get
So for each frequency bin, we can create a LOS integration of a corresponding integrated opacity
where
is the number density of neutral HydrogenPhotoionization
Also need to update code to include photoionization rates as well as direct ionization and radiative recombination…
Photoionization rate can be calculated from
We can approximate
where
and
and we can calculate the luminosity from Planck's law integrated over the area and solid angle…
where
and
Putting this all together, we get
Now assuming a temperature of
, this gives
where
and a photoionization time scale (at the .047 AU) of .0435 yr
Here is the same plot as before, but included are the equilibrium ionization fraction as a function of temperature for material at the orbital separation and with densities of 1e7 and 1e9 particles/cm^{3 }
Now in Bourrier et al, they leave the photoionization rate a free parameter and adjust it between .5, 1, and 5 times the solar value.
They quote the ionizing flux at the solar minimum from
which cite
http://www.aanda.org/articles/aa/pdf/2008/43/aa8810-07.pdf
which have a photoionization rate at 1 AU of 1e-7/s or 115 days. At .047 AU this would be 6 hours which is close to the 8 hours shown in Bourrier et al.
This would imply a photoionizing rate of
currently I am using 3.5808d19
Decomposing gradients into wave strengths
Consider the 2D hydro linearization of the wave equation in each dimension
where
The matrix
has eigenvalues , , , corresponding to left propagating sound wave, contact, shear discontinuity, and right moving sound wave. And we can construct the matrix of left eigenvalues (where each row is a left eigenvector)and the matrix of right eigenvectors
(where each column is a right eigenvector)
and the diagonal matrix of eigenvalues
where each diagonal term is an eigenvalueThe eigenvectors are orthogonal
and we also have
which combining gives
Now since the eigenvectors are orthonormal, we can multiply any left eigenvector by an arbitrary normalization constant, provided we divide the corresponding right eigenvector by the same amount.
Rewriting the equation
Now if we renormalize the left eigenvectors
and the right eigenvectors
then
where
has units of
This gives us the relative strength of each wave at each cell. This looks at gradients in the fluid variables and how they correlate with the various types of waves.
For example, if we look at the Sod Shock Tube test problem whose solution entails a left rarefaction, contact, and right moving shock, we see that the 3 different waves are correlated with the amplitude of the projections of the gradients onto the left eigenvectors associated with each wave.
And if we look at the Brio Wu shock tube which consists of a fast rarefaction, a slow compound wave, a contact, a slow shock and a fast rarefaction, we see the following wave strengths
Meeting update
- Continue to work on implementing ThermalConduction solvers.
- Investigated Newton's Method for solving non-linear PDE's
- Currently, Anisotropic & Isotropic, Implicit, & Explicit Thermal diffusion seems to work on Guassian test. johannjc10092015
- Works in parallel and in AMR, with Hydro turned off. Still need to
- Resolve issue with ghost zoning
- Implement flux limiting
- Run various tests - MTI, Ablative RT, ???
- Plugging away at MHD Wire Turbulence sims using ½ of BGQ. Getting a 26GB frame every 8 hours…
- Ran HD… sim out longer - still hasn't reached a steady state…
Anisotropic Thermal Diffusion
Working on implementing Anisotropic Thermal Diffusion
The 4 columns correspond to Isotropic diffusion, Anisotropic but with
, parallel diffusion, and perpendicular diffusion. And the 3 rows correspond to Backward Euler, Crank Nicholson, and Forward Euler (explicit).The magnetic field is uniform in the +x+y direction and In each case the effective diffusion constant is 1,
, each box is 2 units across, and the total time = 1. Hydro is turned off and these are run using a single core in fixed grid.The pseudocolor is showing Temperature and the initial circle has a temperature ½ that of the surrounding.
movie |
Fiducial parameters for HD209458b
0.013893 | ratio of planet radius to orbital separation | |
0.117818 | ratio of stellar radius to orbital separation | |
0.057598 | ratio of Hill radius to orbital separation (from planet) | |
0.240468 | ratio of bow shock radius to orbital separation (from planet) | |
0.074592 | ratio of planet sonic radius to orbital separation (from planet) | |
0.287192 | ratio of coriolis radius to orbital separation (from planet) | |
22.853798 | ratio of densities at bow shock | |
22.088202 | stellar lambda | |
10.738334 | planetary lambda | |
0.000000 | ||
0.000000 | ||
Orbital separation | 0.047470 | AU |
Mass of Star | 1.148000 | solar masses |
Mass of Planet | 0.690396 | Jupiter masses |
Radius of Star | 1.203000 | solar radii |
Radius of Planet | 1.380000 | Jupiter radii |
Temperature of Star | 999999.671375 | Kelvin |
Temperature of Planet | 10005.732583 | Kelvin |
Density of Star | 1.000000e-13 | g/cc |
Density of Planet | 1.000000e-15 | g/cc |
Orbital period | 3.523762 | days |
Mass loss from Star | 1.768867e-16 | solar masses per yr |
Mass loss from Planet | 3.114146e+09 | g/s |
Mach number of Stellar wind at shock | 6.926105e-01 | |
Mach number of Planetary wind at shock | 1.448805e+00 | |
lScale for sim | 7.097443e+11 | cm |
TimeScale | 3.044531e+05 | s |
rScale | 1.000000e-15 | g/cc |
Location of planet in units of lscale | 1.000000e+00 | |
predicted bow shock radius in units of lscale | 2.406060e-01 |
And here is schematic showing the star/planet surfaces, their sonic surfaces, and the location of the bow shock, hill radius, and coriolis radius. Note the bow shock will likely be closer to the planet due to the planet wind being sub-parker and that the coriolis radius will probably be larger for the same reason.
And here is the same setup, but with a stellar density 10x higher
0.013893 | ratio of planet radius to orbital separation | |
0.117818 | ratio of stellar radius to orbital separation | |
0.057598 | ratio of Hill radius to orbital separation (from planet) | |
0.240468 | ratio of bow shock radius to orbital separation (from planet) | |
0.074592 | ratio of planet sonic radius to orbital separation (from planet) | |
0.287192 | ratio of coriolis radius to orbital separation (from planet) | |
22.853798 | ratio of densities at bow shock | |
22.088202 | stellar lambda | |
10.738334 | planetary lambda | |
0.000000 | ||
0.000000 | ||
Orbital separation | 0.047470 | AU |
Mass of Star | 1.148000 | solar masses |
Mass of Planet | 0.690396 | Jupiter masses |
Radius of Star | 1.203000 | solar radii |
Radius of Planet | 1.380000 | Jupiter radii |
Temperature of Star | 999999.671375 | Kelvin |
Temperature of Planet | 10005.732583 | Kelvin |
Density of Star | 3.211200e-13 | g/cc |
Density of Planet | 3.211200e-15 | g/cc |
Orbital period | 3.523762 | days |
Mass loss from Star | 5.680185e-16 | solar masses per yr |
Mass loss from Planet | 1.000014e+10 | g/s |
Mach number of Stellar wind at shock | 6.926105e-01 | |
Mach number of Planetary wind at shock | 1.448805e+00 | |
lScale for sim | 7.097443e+11 | cm |
TimeScale | 3.044531e+05 | s |
rScale | 3.211200e-15 | g/cc |
Location of planet in units of lscale | 1.000000e+00 | |
predicted bow shock radius in units of lscale | 2.406060e-01 |
Meeting Update
Thoughts on parker wind
The solution for the parker wind is found by solving the equations below
Note however, that the actual radius of the planet (and lambda) does not really matter…
All that a higher lambda implies, is a smaller planet radius (compared to the sonic radius) - and sampling more and more of a subsonic atmosphere before having a hard boundary.
Also - per our discussion the other day, it seems that
So this effect is more pronounced when the orbital sepration is small and when the planet radius and mass ratio is close to - ie Hot Jupiters
How to handle particles, outflows, and point gravity objects
Frame Restarts
In general, when a problem module creates an object - it is easiest on a restart, to let the same objects just be recreated. Then the problem module does not need to be modified - and will still have pointers to the objects. The exception, however, is sink particles - since they can be created by self-gravity, they can wander due to gravitational interactions, accrete, etc…
Now a problem module on a restart, has no way to identify particles as their position, mass, velocity, etc… may have changed. (although we could require problem modules to name the particles). Then problem modules would need to have something like the following:
IF (lRestart) THEN CALL FindSinkParticle("Star 1", Particle) ELSE CALL CreateParticle("Star 1", Particle) Particle%xloc=... Particle%mass=... ... END IF
Then the problem module could continue to setup their outflow objects, point gravity objects, etc…
CALL CreatePointGravityObj(Particle%PointGravityObj) CALL CreateOutflowObj(Particle%OutflowObj) CALL UpdateParticle(Particle)
The alternative is to store their outflow objects, point gravity objects, etc… in the chombo file. Then a problem module would look like
IF (.NOT. lRestart) THEN CALL CreateParticle("Star 1", Particle) Particle%xloc... CALL CreatePointGravityObj(Particle%PointGravityObj) CALL CreateOutflowObj(Particle%OutflowObj) CALL UpdateParticle(Particle) END IF
But now on a restart the user would not easily have pointers to the various outflows and pointgravity objects…
Another complication, is that if a tracer is added to an outflow object that is associated with a particle.
The solution to that is to create the tracers whether there is a restart or not…
Step Restarts
When a restart is triggered, we need a way to roll back any particles to their previous state. It is also possible that particles would have been created or destroyed during the last step, so those would need to be recreated (or deleted).
A particle has additional sub pointers in the derived types that need to be allocated and copied
- Particle
- Outflow Object (part of object list
- Shape
- Profiles
- data
- fields
- Point Gravity Object (linked list)
- Outflow Object (part of object list
For now, let's assume that there is no merging of particles as this simplifies everything.
Before each step:
- Make copies of particles
On a restart:
- Find backed up particles and copy contents
- Remove any new particles and their objets
- Update particle objects with copied values for position, mass, etc…
When we do implement merging, we could move particles (and their associated objects) out of their respective lists. Then on a restart, we could put the particle back without creating new pointers etc…
Planet Wind Orbital Trajectories
I calculated trajectories for planetary wind particles leaving at the escape speed from the planet surface in a corotating frame.
Here are the trajectories for a separation of 400 AU
and here are the trajectories for a separation of 200 AU
… and 100 AU
We have
and
We can calculate an coriolis radius
where
which gives
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/R^{2 outside, and uniform sphere gravity inside. Note that M used for outside gravity does not equal 4/3 pi r}3^{ 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
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
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
Implementing Dipole fields in Planetary Studies
Planetary magnetic fields are dipole fields of the form
or
and
which (if we assume that the dipole is oriented along z) gives a field of
at the pole and a field of
at the equator.
This also corresponds to a vector potential
Now if we assume that
Then our vector potential is
Also, the dipole field used in Matsakos was a constant field inside of the planetary core (r<r_p/2)
They quote a constant field of
at the center would require a potential of the form
which agrees with the potential outside at
So we don't need any smoothing… just a piecewise definition for the potential that switches at
Now, if we want to truncate the field lines so they do not pass out the inflowing boundaries, we need to have the derivative of the potential go to zero. This is tricky since the curl of the potential needs to go to zero in both the theta and r directions. However, the following potential has the desired properties
For
the potential is just 0.This gives
and
for
and forSmooth | UnSmoothed | |
log | ||
linear |
Wrote matlab code to solve Parker Wind
Here is a matlab routine to solve for the parker wind given lambda x is in units of planet radius, and y is velocity in units of sound speed.
Also notified Titos about typo in arxiv preprint that made initial attempts to solve the equation nonsensical.
Also note that the original version of this post assumed that
when in fact
So here is the correct version of the matlab routine
lambda=11.5; N=81; decs=4; xx=logspace(0,decs,N); yy=zeros(N,1); syms y positive; for i=1:N yy(i)=vpasolve(y - log(y) == -3 - 4*log(lambda/2)+4*log(xx(i))+2*lambda/xx(i),y,(xx(i) > lambda/2) * 100+.001); end fid=fopen(['parker_wind_',num2str(lambda),'.dat'],'w') fprintf(fid, '%i\n',N); fprintf(fid,'radius, density, velocity\n') for i=1:N fprintf(fid,'%e %e %e\n',xx(i),1d0/xx(i)^2/sqrt(yy(i)/yy(1)), sqrt(yy(i))); end fclose(fid);
Angular Momentum conservation
I ran the rotating variant of the TrueLove Collapse problem. Here is an animation with density contours and the particles along with their spin vectors.
Here is the angular momentum of the gas, the particles, and the internal spin of the particles in the X, Y, and Z directions (top to bottom)
And here is a closeup of Jz of the gas, the gas+particle, and the gas+particle+spin
I'm a little surprised by the ratio of spin to angular momentum, but I guess in a binary system, there is a lot more angular momentum in the orbit of the stars, then in the rotation of the stars themselves. It is also surprising that the Z component of angular momentum of the particles becomes positive at times - but perhaps that has more to do with the choice of the center of the origin…
Angular Momentum and Sink Particles
Looked over Federrath Outflow Model and noticed significant issues with how we handle angular momentum in sink particles.
Both of us update the mass, momentum, and center of mass the same
But our treatment of angular momentum has glossed over some subtleties
A better treatment is to conserve total angular momentum.
where
and S is the spin.
which if we approximate to first order in
which implies our method is only accurate if
andProposed fiducial parameters for Wire Turbulence study
Proposed fiducial parameters for Wire Turbulence | |
---|---|
Wire diameter (in units of wire spacing) | .25 |
Mach Number | 20 |
Beta | 1 |
Gamma | 1.1 |
Field orientation | y (flow is in x and mesh is aligned with yz unit vectors) |
wire density (in units of wind density) | 1000 |
ambient density (in units of wind density) | .01 |
distance from wind boundary to wire center (in units of wire spacing) | 1 |
Number of zones per wire spacing | 1024 |
Working on Parallel IO
I've been continuing to streamline the chombo io process and now have parallel write working. And it was 2x as fast using 8 processors on a 200 MB chombo file.
Planetary Atmospheres
Ran simulation using hydrostatic profile from r_inner to the boundary, but began to see unusual behavior. Very high temperatures caused the timestepping to slow to a crawl.
And this shows lineouts along the x axis of
- Temperature
- density
- mach number
- |vx|
- Pressure
- |g|
- |grad(p)/rho|
- |grad(p)/rho + g|
- |v dot grad(v)|
I then reran this with the outer boundary being stepped on and saw the following:
Summary of things we've tried
- Pressure Equilibrium - Led to temperature jump which broke radiation solver
- Isothermal atmosphere - Self-gravity solver couldn't reliably converge
- Point Gravity with both inner and outer boundaries - heating of lowest density material
- Point Gravity with just inner boundary - very high temperatures and tiny time steps.
- Also - low res simulations cause density protections regardless of time step due to 'weakness' with 6 solve CTU.
Planetary Atmospheres
Updated boundary conditions to go from rho = 1e-9 to tau = 1e4
Inner boundary: (tau=1e4)
- r=7.815240e+09
- rho= 4.972392e-03
- T= 2.306330e+03
- M_enclosed=2.173174e+30
- P=4.145923e+08
Point particle
- M_point=2.1632e30 = M_enclosed - V(r)*rho
Outer Boundary (rho < 1e-9)
- r=8.181350e+09
- rho=9.488437e-10
- T=1.115419e+03
- P=3.826197e+01
Also, calculated the effective aspect ratio for simulation region: 'w' = 75.672
which gives N_{cells} = N_{r}^{3}w^{2} = 5726 N_{r}^{3} = (18 N_{r})^{3}
So if we want 100 zones in radius, we would be effectively simulating an 1800^{3} sim.
Here are the previous values I used
Inner boundary: (tau=10)
- r=7.9568e9
- rho=3.8901e-5
- T=1.4682e3
- M_enclosed=2.1733e30
- P=2.0648e6 (2.0648 atm)
Point particle
- M_point=2.1732e30 = M_enclosed - V(r)*rho
Outer Boundary (rho < 1e-9)
- r=8.1814e9 (3% of inner boundary)
- rho=9.4884e-10
- T=1.1154e+03
- P=38.2620
PLoit showing mesh
Plot showing density (click for movie)
Plot showing temperature (click for movie)
Convective Stability
Worked on calculating convective instability and opacities. Found description of stellar convection here
Plot showing
as well asPlot showing optical depth from star
Image showing density in log space as well as optical depths of 1 and 2 in black, and the region of stability is outline in pink
And here are the density, temperature and pressure profiles.
New structure of object routines
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.
Testing projections
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
Planetary Atmospheres
Ran 64^{3}+2 simulations of planetary atmospheres and looked at
- fast spinning (10x )
- lower resolution (64^{3}+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 |
Meeting update for 07/03/2014
- Eddie's images
- Working on Parallel HDF5 routines - and cleaning up io/chombo stuff.
- Photo Ablated Clump post
Hypre tolerances
So I ran hypre with varying tolerances for the radiation solve and found that the solution only converged when I set the tolerance to be 1e-10
Movie |
So I looked more closely at what the 'tolerance' actually refers to in hypre. The hypre user guide is not much help…
Convergence can be controlled by the number of iterations, as well as various tolerances such as relative residual, preconditioned residual, etc. Like all parameters, reasonable defaults are used. Users are free to change these, though care must be taken. For example, if an iterative method is used as a preconditioner for a Krylov method, a constant number of iterations is usually required.
I'm assuming that hypre is trying to solve Ax=b and instead ends up solving Ax'=b' where r stands for the residual (r=b'-b) and x' is the approximate solution to x. I'm not sure what the _C means or what the <C*b,b> represents… presumably this is the inner product of some matrix 'C' times b with b.
Iters ||r||_C conv.rate ||r||_C/||b||_C ----- ------------ --------- ------------ 1 1.667442e-07 265.427358 3.581554e-08 2 2.421703e-07 1.452347 5.201657e-08 3 2.037404e-07 0.841310 4.376208e-08 4 1.176484e-07 0.577442 2.527008e-08 5 7.646604e-08 0.649954 1.642440e-08 6 4.446094e-08 0.581447 9.549914e-09 7 2.173844e-08 0.488933 4.669272e-09 8 1.033716e-08 0.475525 2.220354e-09 9 5.190075e-09 0.502079 1.114794e-09 10 2.514604e-09 0.484502 5.401202e-10 11 1.291694e-09 0.513677 2.774473e-10 12 6.662719e-10 0.515812 1.431108e-10 13 4.688375e-10 0.703673 1.007032e-10 14 2.451621e-10 0.522915 5.265918e-11 15 1.654131e-10 0.674709 3.552962e-11 16 8.812279e-11 0.532744 1.892819e-11 17 5.478943e-11 0.621740 1.176841e-11 18 7.612170e-11 1.389350 1.635043e-11 19 1.416387e-10 1.860687 3.042304e-11 20 2.346052e-10 1.656364 5.039163e-11 21 4.389153e-10 1.870868 9.427609e-11 22 1.170443e-09 2.666672 2.514034e-10 23 3.069051e-09 2.622127 6.592117e-10 24 3.268667e-09 1.065042 7.020879e-10 25 3.349935e-09 1.024863 7.195436e-10 26 1.136404e-09 0.339232 2.440919e-10 27 2.175246e-10 0.191415 4.672284e-11 28 7.067671e-11 0.324914 1.518089e-11 29 2.139085e-11 0.302658 4.594611e-12 30 7.493659e-12 0.350321 1.609588e-12 31 3.051552e-12 0.407218 6.554531e-13 <C*b,b>: 2.167928e+01
So I looked for an explanation of convergence for linear system solves and found this document
The condition number of a matrix is related to how much a relative error in the rhs affects the relative error in the lhs. ie
which would imply that
In general we have
and
indicates how close to singular a matrix is. (ie a singular matrix would have an infinite condition number)Also because of machine precision, there is an unavoidable error related to the condition number
AstroBEAR solveres
Now we have two solvers available in the code… One is GMres and the other is PCG (preconditioned conjugate gradient). GMres does not converge unless the tolerance is set to ~ 1e-3 - but then it does nothing. PCG apparently expects symmetric positive definite matrices - but the matrices are certainly not symmetric. Not sure how much that matters… I'm also not sure whether PCG is actually doing any preconditioning automatically. I am not setting the preconditioner etc…
Also the actual solution vector 'x' is in terms of the Radiation Energy, and not radiation temperature - and since E ~ T^{4}, there and there is a 3 order range of temperatures - there will be a 12 order range in Energies - and the relative error in the vector norm will be dominated by the larger values… A better metric would be the maximum local relative error…
Time steps for radiation transfer in AstroBEAR
See FluxLimitedDiffusion#no10 for a description…
This image shows the absorption time scale, the equilibrium time scale, the diffusion time scale, and the resulting 'safe' time step for the 1D planetary atmosphere setup.
If we use the single temperature solution, then we can switch to just the diffusion time scale which looks like it is 1000x longer? This might work for the two temperature model as well?
Paper accepted to ApJ
- Submitted paper to ApJ
- Found very useful tool for visualizing latex diffs - appropriately named latexdiff
- See marked up pdf
- Plan to look at single temperature radiation solver this week.
1D Planetary Models
I switched to 1D - and set the ambient density to 1e4 what it was. However I still got odd behavior. After looking at the hypre matrices and vectors, I discovered that there was a typo in the code…
MinTemp*TempScale
should have been
MinTemp/TempScale
If TempScale is 1, this does not matter, but when it is 1e6 it creates problems
In 1D there is no 'self-gravity' but the simulations are able to evolve on a dynamical time scale.
2D results
Timescales for Planetary Simulations
Time Scale | Equation | In the planet | In the ambient |
Sound Crossing Time | 1.5 | 1e-3 | |
Orbital Time | 2.2 | ||
Radiation Diffusion Time | 2.3e-8 | 9.2e5 | |
Radiation Equilibrium Time | 3.8e-17 | 5.1e-4 | |
Light crossing time | 1.1e-5 |
Notes from today's meeting with Jim Kasting
Background
- Habitable zone limited near by by runaway greenhouse at 1 AU and limited far away by the Maximum CO2-H2O greenhouse effect at 1.8 AU
- With H2 green house gas you can extend out to 10 AU (3 M_earth with 40 bars of H2 atm)
- Homopause (turbopause) at 100 km - where turbulent mixing ceases and atmosphere is not well mixed. Transport of hydrogen governed by molecular diffusion
- Exopause at 500 km - where gas becomes collisionless and fast moving hydrogen particles can escape. Velocity distribution no longer Boltzmann. Rate governed by Jeans evaporation rate.
- On Earth, exopause is fairly warm - and Jeans evaporation is quick - so H2 loss rate controlled by diffusion through the Homopause
- On early Mars, the exopause is colder and therefore Jeans evaporation is slower and might have been responsible for regulating H2 loss. Going from 1D to 2D reduced jeans evaporation rate by ~4 which could increase the steady state H2 concentration to the 20% needed to create enough of a green house effect to place early Mars in the habitable zone - and explain the fluvial features on Mars.
Goal
Calculate an upper limit on the escape rate from the Exopause and show that it is lower than the diffusion limit and hopefully consistent with the 20% concentration necessary to explain the fluvial features on Mars
Complexities of the physical model
- Geometry - Obliquity, magnetic fields,
- Multi-species
- Chemistry
- diffusion
- non-Lte
- heating and cooling
- multi-line
- collisionless - moment equations from the boltzmann equation
Possible simulations
- 1D - 2 fluid with most of the physics (apply an approximation to model the moments of the velocity distribution where it becomes collisionless) to demonstrate the h2 loss rate < diffusion limit
- Evenly spaced in log altitude (not radius)?
- 3 points per pressure scale height
- multi-fluid?
- diffusion
- line transfer - would need multi-bin line transfer
- cooling
- equation of state
- Redo Stone and Proga in 3D and measure the Jeans evaporation rate.
- MHD
- rotation
Self gravitating equilibrium
Self-gravitating equilibrium profile using density profile and external pressure
Mathematically, we are solving the equation of hydrostatic equilibrium
subject to the constraint that
.Since
we can rewrite the equation for HSE asand we can express the enclosed mass as
Now we can discretize equation 2 as
where
and equation 1 as
but a taylor series expansion shows that the error at each step in the integration goes as
so we need to limit out point spacing - or improve the accuracy of our discretization.
If we apply linear interpolation to our density profile, we can discretize equation 2 as
and equation 4 becomes a bit of a monster - and we have to worry about limits as r → 0.
if (r0 == 0d0) then dp=drho*((pi*drho*dr**4)/4d0 + (4d0*pi*rho*dr**3)/9d0) + (rho*(pi*drho*dr**3 + 2d0*pi*rho*dr**2))/3d0 else epsilon=log(r0/(r0+dr)) dp=((-36d0*drho*r0)*epsilon*M0*dr + (-36d0*drho*r0**2)*epsilon*M0 + (- 12d0*pi*drho**2*r0**5 + 48d0*pi*rho*drho*r0**4)*epsilon*dr + (-12d0*pi*drho**2*r0**6 + 48d0*pi*rho*drho*r0**5)*epsilon + (36d0*rho - 36d0*drho*r0)*M0*dr + (9d0*pi*drho**2*r0)*dr**5 + (17d0*pi*drho**2*r0**2 + 28d0*pi*rho*drho*r0)*dr**4 + (2d0*pi*drho**2*r0**3 + 64d0*pi*drho*r0**2*rho + 24d0*pi*r0*rho**2)*dr**3 + (- 6d0*pi*drho**2*r0**4 + 24d0*pi*drho*r0**3*rho + 72d0*pi*r0**2*rho**2)*dr**2 +(- 12d0*pi*drho**2*r0**5 + 48d0*pi*rho*drho*r0**4)*dr)/(36d0*r0**2 + 36d0*dr*r0) end if data(i,i_Press)=data(i+1,i_Press)+ScaleGrav*dp
Momentum Conserving Self Gravity
Making the momentum be conserved requires recasting the gravitational force as a total differential
so we need to find such that
Since
we can substitute for
and we havewhich is equivalent (in 1D) to
where we can identify the equivalent momentum flux tensor as
- In more than 1D we have
or we can write it as
Momentum Conserving Self Gravity Source Terms in AstroBEAR
- q - initial state
- qLx - left edge of x - interface in predictor
- qLy - lower edge of y - interface in predictor
- qRy - upper edge of y - interface in predictor
- q2 - intermediate cell centered quantities
- q2Lx - updated left edge of x - interface used for final riemann solve
- fx - predictor x flux
- f2x - final x flux
- q3 - new values
Normal CTU:
- qLx = Reconstruct(q), etc… (for all 6 qxx states)
- fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
- q2 = q+fx+fy+fz
- q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
- f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
- q3=q+f2x+f2y+f2z
Strang split self-gravity:
- q = q + sg(q,phi)*hdt
- qLx = Reconstruct(q), etc… (for all 6 qxx states)
- fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
- q2 = q+fx+fy+fz
- q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
- f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
- q3=q+f2x+f2y+f2z
- q3=q3+sq(q3,phi)*hdt
where
sg=sg_x+sg_y+sg+z and sg_x has two non-zero terms…
Momentum conserving self-gravity
- qLx = Reconstruct(q), etc… (for all 6 qxx states)
- qLx(px) =qLx(px)+sg_x(q,phi)*hdt
- qLy(py) = qLy(py)+sg_y(q,phi)*hdt
- qLy(pz) = qLy(pz)+sg_z(q,phi)*hdt
- fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
- q2 = q+fx+fy+fz + sg(q,phi)
- q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
- q2Lx=q2Lx+sg_y(q,phi)+sg_z(q,phi) (for all 6 q2xx states)
- f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
- f2x(p)=f2x(p)+SG_PFlux_x(phi) (same for y and z)
- q3=q+f2x+f2y+f2z+SG_E(f2x(rho),f2y(rho), f2z(rho), phi)
The final updates for momentum include true 'fluxes'
SG_PFlux_x(phi) has only 1 non-zero term
and an energy update term that is derived from mass fluxes.
SG_E(f2x, f2y, f2z, phi)
For a description of an algorithm that also conserves energy see http://adsabs.harvard.edu/abs/2013NewA...19...48J
Meeting Update
- Worked on scaling implementation
- Still testing merged version of code
- Working on a hydrostatic Jupiter… Getting funny behavior in ambient - but density contrast is extremely high…
Setting Scales in AstroBEAR
I added the ability to set the TimeScale in physics.data, however this requires more careful testing of which combination of parameters can be set…
There are some constraints between scales… ie
pScale=nScale*TempScale*Boltzmann
rScale=nScale*Xmu*amu
lScale/TimeScale=sqrt(pScale/rScale)
and so on…
There are currently 7 scales that can be set - and three equations that constrain them.
Xmu, rScale, nScale, TempScale, TimeScale, pScale, and lScale
Therefore there can be 4 parameters set in physics.data - though not certain combinations. To test this we can construct a linear system in the log of each parameter…
log(pScale) - log(nScale) - log(TempScale) = log(Boltzmann)
log(TimeScale) = log(1d10)…
Which would give a 7x7 linear system that can potentially be solved. If the determinant is zero, then the 4 parameters that were set by the user are not independent. If they are - then the matrix can be inverted and solved for the remaining parameters.
Setting Orbital Parameters
I've had to do this a couple of times - and thought it might be good to post this here…
Given 2 bodies in orbit about each other - calculate the initial positions and velocities given the semi-major axis, the masses, and the eccentricity.
If we assume the orbit is about the z-axis, and that aphelion occurs when both bodies are on the x-axis then at aphelion the velocities must be in the y-axis and the locations are on the x-axis. This leaves us with 4 unknowns
.The center of mass residing at the origin implies that
and we have
which can be solved for
We also must have zero net momentum so
which implies that
Which is consistent with them having the same angular velocity at aphelion.
Now
And we know that twice the entire area of the ellipse is swept out in one orbital period
we have
or
which then gives us
which together with
define
and
Equation of State functions
ideal gas equations
and definitions of sound speed and enthalpy
Here is a list of all of the various functions that are equation of state dependent in the EOS module
- P(e)
- P(h)
- e(P)
- h(P)
- T(e)
- T(p)
- e(T)
- p(T)
- c_{s}(e)
- c_{s}(P)
Which can also be obtained from
- T(e) and its inverse
- p(T) and its inverse
- c_{s}(T)
- h(T) and its inverse - though not absolutely necessary
Meeting update 09/03
- Ideas on protection improvements
- Radiation update
- Finished Explosion test (does not test explicit terms)
- Both explosion test and marshak wave test seem to work with AMR
- Implemented advection test (tests both implicit and explicit terms)
PN Stuff
In the PN stuff the wind and jet temperature are << shock temperature - which can cause feedback at the spherical boundary… Here is an image of the same run with an outflow and a jet in which the temperatures were changed from 1000 to 10000K. The shock temperatures are ~ 100000K
And here's a lemon
Radiation interface
To use radiation you need to
- Add configuration options to global.data
!============================================================================================ !Variables used for radiative transfer module !============================================================================================ &RadTransferData solver = 2 ! Which solver to use for linear system solve? StructPCG=1, StructGMRes=2, StructSMG=3 [1] tolerance = 1e-4 ! Solver convergence tolerance [1e-6] printlevel = 0 ! Sets print level for hypre solver. Increasing this will make solvers more verbose. [0] MaxIters = 10000 ! Maximum number of iterations to attempt convergence. [1000] hverbosity = 0 ! Anything > 0 will cause linear systems to be outputted as matrix, rhs, lhs files. [0] mthbc = 6 2 2 1 2 2 ! What type of boundary conditions to use for gravitational potential ! format: (x1, y1, z1, x2, y2, z2) ! 0-Free Streaming, ! 1-Extrapolated, ! 2-Periodic ! 3-Reflecting ! 4-User defined boundary ! 5-User defined radiative flux RadFlux = 824931.169876 ! User defined radiative flux used for mthbc 5 Crankiness = 1 ! parameter that adjusts from backward euler (0) to crank-nicholson(1) time stepping. (2 is forward euler but nothing over 1 should be used - and 0 should be more stable but less accurate) cfl_rad = .01
- Set iRadiation in physics.data
iRadiation = 1 ! 0 no radiation, 1 implicit terms only (heating/cooling), 2 implicit and explicit (includes pressure forces and advection) iLimiter = 1 ! LP_FLUX_LIMITER=0, DIFFUSION_ONLY=1, NO_FLUX_LIMITER=2
- Set opacity related variables
!================================================================================ ! Opacity related section !================================================================================ PlanckOpacityType = 1 ! [1]-Specific Opacity, 2-Power law opacity PlanckSpecificOpacity = 1000 !computational units RoselandOpacityType = 2 ! [1]-Specific Opacity, 2-Power law opacity RoselandSpecificOpacity = 1 !computational units Roseland_dpow=2 Roseland_Tpow=-3
- Make sure to intialize the radiation field. Equilibrium value is given by scalerad*T4. Not all objects will support this.
Jets etc...
Thought about updating the outflow module to be more transparent and to help launch a spherical outflow with an overlaid collimated outflow or conical outflow…. See attached powerpoint slide for discussion…
Also updated Ruka's module to launch jets in 2.5D (not using multiple outflow objects though)
Meeting update for 07/08/2013
So Marshak boundary conditions are mixed
One could instead do
or
so that
In any even, if we discretize this mixed equation we get (E_g is ghost zone, E_i is internal zone at left edge)
which we can solve for
as
and we have
If we then plug this into the time evolution equation ignoring portions that are unchanged we get:
where
we get
where
is the time averaged energy and depends on the time stepping (crank-nicholson or forward euler)and finally we arrive at
This is very similar to what we normally get ie
and noting that
for the limiteror flux is just
Note that
for
and that if
we revert to what we would expect if there was no Radiation energy in the ghost zone.Fixing latex tags
Thought this might be handy if we ever want to do any batch editing of wiki pages.
#!/bin/bash for i in `trac-admin ./ wiki list | tail -n +4 | head -n -1 | awk -F ' ' '{print $1}'` do trac-admin ./ wiki export $i | perl -pe 's/\[\[latex\(([^\$])(.*?)\)]]/\[\[latex\(\$$1$2\$\)]]/g' | trac-admin ./ wiki import $i done
It basically selects all of the wiki pages and then parses them through a perl regex parser that replaces
[[latex($<some stuff>$)]]
with
[[latex($<some stuff>$)]]
but leaves
[[latex($<some_otherstuff>$)]]
alone. This should fix most things - though it may miss things like
[[latex($ $a=b$$)]]
And here is a script for batch-editing blog posts. Since the blog posts are managed through a plugin they have their own table - and trac-admin does not provide any easy front end… So you should only modify the trac database after shutting down the web server.
#!/bin/bash #make sure to shut down wiki server sqlite3 db/trac.db << EOF .out fullblog_dump .dump fullblog_posts .exit EOF perl -p -i -e 's/\[\[latex\(([^\$])(.*?)\)]]/\[\[latex\(\$$1$2\$\)]]/g' fullblog_dump sqlite3 db/trac.db << EOF DROP TABLE fullblog_posts; .read fullblog_dump .exit EOF
This will unfortunately modify this blog post as well
Meeting update for 07/08/2013
- Reformatted poster. Adam - does this look ok?
- Modified open mp calls to launch 1 to split threads once per grid update. Line label refers to the total number of threads in use per node. The solid lines are times for forking the threads at every do loop - while the dotted lines are times for 1 fork at the start of a grid update.
Thinking about evanescent solutions to self-gravity...
So we start with the linearized differential equations
which we Fourier transform in physical space and solve the third equation for
and substitute it into the second.which gives us a matrix equation for the time evolution operator. Now we just need to find eigenvalues and eigenvectors.
The characteristic equation is
and we have
with eigen vectors
However given the presence of self-gravity we see that
can be real as well as imaginary.So for stable waves we have
where is real and positive. And there are two solutions… (a left and right going gravity-modified sound wave. Gravo-accoustic modes)And if
where is real - then we have two solutions as well.these correspond to an exponentially increasing collapse and an exponentially decreasing decay (which is just the time reversed version of the 1st) Note the velocity field is phase shifted.
Radiation Transport Working
The Implicit radiation transport solver is working in parallel for fixed grid.
It takes hypre a long time to converge initially. I'm assuming this is because the opacity is so low yet the radiation field is very far out of equilibrium.
So after some investigation it looks like the matrix tends to be very ill conditioned (condition numbers of 1e9).
The solution is to start with a small time step and change the cfl relaxation parameter to be ~ 1e-3… If the hydro time step (or the radiative cooling time step) is ~10 orders of magnitude longer than the diffusion time scale - and the system is not in thermal equilibrium - then the matrix becomes too ill-conditioned for hypre to solve.
Not that anyone should do this carelessly
So I wanted to modify the list of reports to have the useful ones at the top… And since there was an open slot (report #2) I used sqlite3 to interface with the trac database in /data/trac/astrobear/db (after making a backup etc…) The id needs to be unique - otherwise one could break the trac database….
johannjc@clover db > sqlite3 trac.db sqlite> update report set id=2 where id=27; sqlite> .quit johannjc@clover db >
mathjax working again
I successfully installed a version of mathjax - that supports the inline features.
You can use
\(a=b\)
for inline equations ie. \(a=b\)
or
$$a=b$$
for block equations
note that the old latex macro still renders images on the server
[[latex($a=b$)]]
gives
so the only problem will occur if people have used
[[latex($a=b$)]]
(or removed their dollar signs when we switched from tracmath to tracmathjax)
timinginfo
Added timers to every do loop in hyperbolic update
Here is the break down in terms of percent time and threaded speedup (using 32 threads)
id | serial time | speedup | 32 thread time | percent of total | description |
82 | 3.231047 | 9.128274 | 0.353960343 | 9.89% | predictor riemann solves |
80 | 3.121115 | 8.762267 | 0.356199486 | 9.55% | predictor riemann solves |
84 | 2.997481 | 9.587972 | 0.312629303 | 9.17% | predictor riemann solves |
69 | 2.919109 | 22.452476 | 0.130012788 | 8.93% | eigen system |
213 | 2.752028 | 9.994905 | 0.275343087 | 8.42% | final riemann solves |
210 | 2.63407 | 9.832083 | 0.26790559 | 8.06% | final riemann solves |
216 | 2.617439 | 9.956359 | 0.262891183 | 8.01% | final riemann solves |
23 | 1.504078 | 23.787341 | 0.063230186 | 4.60% | characteristic tracing |
25 | 1.499858 | 23.942636 | 0.062643812 | 4.59% | characteristic tracing |
20 | 1.495662 | 24.030097 | 0.062241197 | 4.58% | characteristic tracing |
19 | 1.38434 | 23.166956 | 0.059754937 | 4.24% | 2 spatial reconstruction |
18 | 0.689309 | 22.116419 | 0.031167297 | 2.11% | spatial reconstruction |
2 | 0.128533 | 19.941572 | 0.00644548 | 0.39% | cons_to_prim |
97 | 0.110934 | 20.652223 | 0.005371528 | 0.34% | upwinded emfs predictor |
96 | 0.099056 | 19.557486 | 0.005064864 | 0.30% | upwinded emfs predictor |
95 | 0.098898 | 19.709309 | 0.005017832 | 0.30% | upwinded emfs predictor |
225 | 0.091149 | 19.994068 | 0.004558802 | 0.28% | upwinded emfs |
159 | 0.090894 | 20.499695 | 0.00443392 | 0.28% | prim_to_cons2 |
147 | 0.090303 | 20.1791 | 0.004475076 | 0.28% | prim_to_cons2 |
And here is a nifty plot showing the time spend on each loop vs threading efficiency.
The lower right cluster is the 6 different riemann solve loops (3 dimensions for predictor and final) which shows the most promise for improving the speedup. The upper right is the calculation of the eigen system which seems to scale quite well. Then in the center is the spatial reconstruction and characteristic tracing
The weighted average efficiency is around 12% - mostly due to the riemann solves.
So the problem was the allocation of temporary variables inside of the riemann solvers. Different threads would allocate these arrays at the same time and they would end up being close in memory which presumably lead to false sharing. In any event, the scaling of the riemann solvers is much better - and is a small portion of the runtime. Now the threaded bits get an average speedup of 21. However the non-threaded parts now make up 90% of the time…
Here is an updated graph showing the scaling performance of the same do loops.
I also looked at reducing the time spent in calc-eigens. I suspected that the manipulation of a MaxVars by MaxVars static array might be responsible so I reduced MaxVars from 20 to 10. The result was a significant speedup
I also added timers for the source routines (even though there was no source terms) and they seemed responsible for the 90% of the time for the threaded runs. If I comment those out I get a 16 fold speedup with threads on a bluegene node for a fixed grid computation.
I next tried an AMR run but that seems to be running slower due to extra refinement. See attached screenshot.
Obviously the threading is producing noise which is triggering refinement.
So it turns out I was attempting to use allocated variables with the local module scope as private within a threaded do loop… But apparently this was not working with openmp. I changed them to be a static size and declared them with the subroutine scope and it seems to work fine now. And on a 16^{2}+2 field loop advection run, the single node threaded version was 20% faster than the single node mpi version
Meeting Update
- Compared astrobear with loop and sweep direction reordering
- 10% speedup on 2D 16^{2}+4
- 20% speedup on 3D 64^{3}+0
- Next going to look at ways of using shared memory efficiently.
Meeting update
- Flux limited diffusion is working in static case… need to still add explicit source terms. movie
- Working with Baowei and optimizing AstroBEAR's hydro update code.
- Confirmed correct growth of RT instability (yet again...)
Findng the RT Instability growth rate
- First load the database into visit
- create a vector expression called gradrho which is 'gradient(rho)'
- create a scalar expression called gradrhoy which is 'gradrho[1]'
- make a pseudocolor plot for gradrhoy
- make slice(s) along y axis. (So for a 3D data set this would be 2 slice operators: orthogonal to x with 50 percent first — yz plane, orthogonal to z with 50 percent then — xy plane, or vice versa. Also make sure unmark 'Project to 2D' )
- Then perform a 'max' query using the 'actual data' (not the 'original data')
- This will spit out the maximum value of the gradient as well as the position.
- You could advance the time slider and repeat the query - but this is tedious so instead you can
- clear the query results
- Open the command dialog (from Control menu)
- copy and paste the following script
for state in range(TimeSliderGetNStates()): SetTimeSliderState(state) Query('max', use_actual_data=1)
- execute the script
- At this point the query dialog should have the location of the peak density gradient for each frame. Copy all the printed lines like this one:
VisIt: Message - gradrhoy -- Max = 21.989 (zone 4584 in patch level3,patch7 at coord <0.151042, 0.5, 0.1>)
Save this to a file and name the file visit_grepped.txt
- Now you can clean up the data using the following bash one-liner
sed 's/</ /g' visit_grepped.txt | sed 's/>)/ /g' | sed 's/,/ /g' | awk '{print .1*(NR-1), $15}' > data.out
Where we have used the awk expression
.1*(NR-1)
to calculate the time for each frame assuming each frame is .1 computational time units apart.
- Then in gnuplot you can load the data and fit it to any function you please… However you may want to also trim the data file so that you only have points during the linear growth phase.
- Then in gnuplot
g(x)=a+b*exp(c*x) fit g(x) './data.out' via a,b,c plot g(x), 'data.out'
Meeting update
- I've been studying the BlueGene Q architecture and the general tuning of codes to run on BlueGene. We should be able to get much better performance with a minor amount of effort. See this wiki page
- Hopefully with a little work, AstroBEAR can run more efficiently (6x) and I can resume the magnetized colliding flows run.
- FLD is on back burner until colliding flow run has resumed.
Meeting Update 04/16
- MHD Mixed Results
- Radiation Transfer coded, ready for testing.
MHD mixed results
The MHD run ran out to frame 99 = 12 Myr before beginning to throw out Nan's in the reconstruction… So I copied the 5 GB file and loaded it into Visit - and looked at the results…
First the Good
—- Here are density isocontours
movie |
—- And here is a close up of a single core with Field lines
movie |
—- And here is the same but with shifted velocity stream lines
movie |
And the bad
—- The velocity shows a red/black pattern (that is brought out by the color table
—- But the density also shows red/black fluctuations of order 1
Glimmer of hope
However given that the field lines look ok
—- and the velocity stream lines show that this material is not reaching the collision region
—- And that it only seems to be on the root level
Perhaps there is a way to salvage the run… Obviously the boundary conditions need to be modified - but perhaps some diffusion/filter could be used to smooth out the 'ambient' to allow the run to continue?
After more analysis I believe what happened was that early on, outward velocities were achieved in the ambient around the incoming flow. Since the boundary conditions were 'outflow only', the boundary conditions were extrapolated for these cells - I don't know why gravity did not reverse this effect, or why the neighboring cells were able to maintain their 'inflow'… But given the red-black pattern - i suspect hypre truncation errors were somehow amplified by the boundary conditions… However, I don't know why this was not seen in the hydro runs..
movie |
—-
movie |
—-
So I did three fixed grid simulations to see what was responsible
movie |
And it looks like magnetic fields make the difference - although self-gravity seems to exacerbate the problem.
So I also tried a reflecting wall boundary condition to see if that would prevent the striations - and that will finish soon. But in the meantime…
Fixing the striations
It is high frequency noise that could potentially be filtered? Though we don't want to introduce any divergence… Another option is to add local diffusion… ie \(q_{i,j,k} = q_{i,j,k}+ G_{i,j-½,k}-G_{i,j+½,k} +H_{i,j,k-½} - H_{i,j,k+½} \) where \(G_{i,j+½,k} = \alpha_y \left ( q_{i,j,k} -q_{i,j+1,k} \right ) \) and so on…
This would give a new value of \(q_{i,j,k}=q_{i,j,k} \times \left (1 - 2\alpha_y-2\alpha_x-2\alpha_z \right ) + \alpha_x \left (q_{i-1,j,k}+q_{i+1,j,k} \right ) + \alpha_y \left (q_{i,j-1,k}+q_{i,j+1,k} \right ) + \alpha_z \left (q_{i,j,k-1}+q_{i,j,k+1} \right ) \)
This can be repeated as needed, but setting \(\alpha_x = \alpha_y = \alpha_z = 1/12\) should give reasonable results…
Now to keep the B-field divergence free, we need to calculate an 'emf' to update the B-field with… \(\mathbf{B}=\mathbf{B}+\nabla \times \mathbf{E}\). Now we want to diffuse the B-field using \(\mathbf{B}=\mathbf{B}+\alpha \nabla^{2 \mathbf{B} = \mathbf{B}-\alpha \nabla \times \left ( \nabla \times \mathbf{B} \right ) \) so we can identify \(E=-\alpha \nabla \times \mathbf{B} \) This can be discretized as \(E}z_{i+½,j+½,k} = -\alpha \left (B^{y_{i+1,j+½,k}-B}y_{i,j+½,k} + B^{x_{i+½,j,k}-B}x_{i+½,j+1,k} \right ) \) and so on…
This would give
Now if we are only worried about Bx, and smoothing in y and z we can simplify this to
This should keep the ratio of \(\rho / B_x\) constant… consistent with flux freezing.
We can also reduce the diffusion where the mixing parameter is high - or where the grid is refined… if we use a diffusion factor \(f = 1 - e^{{1-4 \left ( \frac{\rho_L+\rho_R}{\rho}\right )}2 } \) and then set \( f = 0\) if a cell is refined. In addition we can force diffusion near the inflow boundaries by scaling \(f = \max(f, \left (\frac{r-.75}{.25} \right ) ^{2 ) \) where \(r \) is a parameter that goes from 0 at the midplane - to 1 at the boundaries. }
Bayes theorem...
So, Bayes theorem says that
so as an example, if we knew that a jar had equal odds of having either 3 or 4 numbered jelly beans, and we drew one at random that was numbered 1, we could calculate the odds of there being 3 as opposed to 4 jelly beans.
In this case
andNow the probability of drawing the number 1 jelly bean IF there are 3 jelly beans is
and the odds of there being 3 jelly beans (as opposed to 4) are but what is ? Well, we have to add all the ways we could draw a '1'. where . So then we have.
Now
was drawn from the discrete set of events , but we could imagine a continuous set of events (\x$)]]. Then Bayes theorem reduces toSo let's say someone placed you randomly along a line that could be anywhere from 0 to L meters long with equal probability… And you found yourself at 1 meter from the start. What does the probability distribution for the length of the line look like now?. Well initially
. Now, for and otherwise. So we have
The following image shows the prior probability distribution for line lengths of 10 and 20 as well as the new probability distributions after finding oneself at 1 meter. Also shown are the 50th percentile and the average expected length. For example, if the initial distribution is between 0 and 10 meters, then the modified distribution has an average expected length of 3.90 and the length will be less than 3.1 50% of the time.
The expected value and the 50th percentile are dependent on our initial guess as to the maximum line length. For example if we assume the length could be from 0 to 20 meters instead, then the 50th percentile is at 4.5 and the expected value is at 6.3. In general the expected value will go as
and the 50th percentile will go asThe above analysis is for a distribution in which there are equal probabilities for being found at any point along the line.
, and an equal likelihood for the line being any length less than L , however we can also consider a system in which the probability of being at any given point grows exponentially with distance (or time). while still keepingFor example if you chose a random human in history, the odds of finding a human living at time t, assuming the human civilization grows exponentially for x years before ending, should be something like
In this case we have
and . So we haveThis gives an expected value of
and an . In this case the expected value is fairly insensitive to and quickly asymptotes to 2. So in an exponentially growing system, a random selection will be heavily skewed towards the end and would give the expectation that the system would end within one more timescale.So if civilizations continue to grow indefinitely, then we are likely near the end of ours… however if civilizations cap out at a fixed number of people, then our expected lifetime becomes dependent on the maximum expected civilization lifetime…
So perhaps this is an argument against colonizing other planets? Or curbing the population growth?
Doomsday statistics worked out...
Assume there is some rate for civilizations to be born A, and each civilization has some chance D to die before a time B at which point they go extra-planetary… Then at any given time T there would have been A*T civilizations born, A*(T-B)*(1-D) extra-planetary civilizations, and at most A*B earth-like civilizations.
Extra planetary | earth like | Dead | Total |
A*(T-B)*(1-D) | < A*B | > A*(T-B)*D | A*T |
If we assume a constant birth rate for all civilizations then the odds that you are born into an earth-like civilization instead of an extra planetary civilization is at best p = B/(T-(T-B)*D) or a lower bound on D = (pT - B)/(pT-pB)
We can guess D and calculate how lucky we are to be in a young civilization, or we can guess p and estimate how unlucky we probably will be.
If B = 1e4 yr, and t = 1e10 yr and…
If d = 0 then we are 1 in 1,000,000 and if d = .5 then we are 1 in 500,000
On the other hand if we assume that p = .5 then our odds of surviving would be 1 in 500,000
Of course we have not used the fact that our civilization has already been around for Y years … Or how close are civilization is to going extra-planetary or how close Y is to B - and it does not appear the authors of the paper considered this either - but as it turns out it can be fairly important if the existential threats are similar to planet killing asteroids that are likely to strike at any time.
Let's assume that lifetimes of civilizations follow some exponential distribution (a*exp(-a(t-t0)) for t < t0+B
Furthermore if (1-D) fraction of civilizations survive for B years, and the lifetime of civilizations that live less than B follows an exponential distribution, then (1-exp(-lambda B)) = D which gives lambda = -ln(1-D)/B. So given that our civilization is Y years old, the odds that we are destroyed before B is (exp(-lambda Y)) - (exp(-lambda B)) = (1-D)^{(Y/B)} - (1-D)
When Y = 0, this gives D consistent with the paper… but obviously when Y = B, this gives 0 as we would be in the process of colonizing other planets…
If we expect to be extra planetary in 100 yrs, then Y = .99 B and if we assume p = .5 which gives D = (1-2e-6) then this gives our odds of us not going extra-planetary at only 2e-7! Which is much better odds than (1-2e-6)…
Of course the value used for D assumed that civilizations that die live for B years, and not the average value from the exponential. We could calculate the average as a function of lambda, and then use p to get lambda, and then use lambda to get D… But given that we are close to B, means that it must happen with some frequency - which means that most civilizations are long lived, and we must just be lucky to have found ourselves in a pre-extra planetary civilization
Thinking about jars of numbered jelly beans
Let's assume we have a jar of numbered jelly beans of unknown size… And we randomly draw a jelly bean with a 1 on it… Does this tell us anything about the jar of jelly beans? Well - we now know the jar had at least 1 jelly bean with absolute certainty… But what are the odds the jar only had 2 jelly beans to choose from - or 3, 4, etc… Or if we were someone gave us even odds on their being more than X jelly beans - should we take it?
Intuitively, we might expect that the odds that there are more than 1,000 jelly beans are slim at best… To actually calculate the odds, we need to assume some distribution for the possible number of jelly beans in the jar… Let's assume that any number of jelly beans is equally likely up to some limit L. Then for each L', the odds of us drawing a jelly bean with the number 1, are just 1/L'. Now to determine the new expected distribution of jelly beans, we can weight each expectation by the relative odds of drawing a 1 compared to the overall odds of drawing a 1. The overall odds are just 1 + ½ + 1/3 + ¼ + 1/5 + … 1/L = H_{L} where H is the harmonic series. The odds of there being L' Jelly beans were just 1/L, but the odds of drawing a 1 with L' jelly beans is 1/L', so the new odds of there being L' jelly beans is 1/L' / H_{L} instead of 1/L. So we can ask ourselves what are the odds of there being only 1 jelly bean? If we thought there could be somewhere between 1 and 100 jelly beans, then the odds are just 1/H_{100} ~ 20%. The odds that there are either 1 or 2 jelly beans would then be (1/1 + ½) / H_{100} = H_{2}/H_{100} ~30 and the odds are 50/50 that there are fewer than 7 Jelly beans…
Now if there could be N jelly beans with equal probability, and we drew a jelly bean numbered 1, what would be the upper limit of beans that gave us 50/50 odds? This is just H_{x}/H_{N} = .5 Since the harmonic series is sum(1/n) can be approximated by integral(1/x) = ln(x), the harmonic series asymptotes to ln(N), and we have ln(x)/ln(N)=.5 or x ~ sqrt(N). So for L=100, we would estimate 10 beans would be a sufficient upper limit… (It was 7)… But for 10,000 beans, even odds say there are no more than 100 beans… But we can make no useful predictions without having some upper limit on the possible number of jelly beans… We can also ask what the average number of jelly beans we should expect is - as opposed to the 50/50 cutoff… This would just be sum (L' 1/L')/H_{L} = L/H_{L} which assymptotes to L/ln(L) …
What about exponential distibutions instead of linear
So now if the jelly bean jar is flared, and divided into L vertical sections and there are exp(z) jelly beans per section. and each jelly bean is labeled by it's vertical section. Then the odds of drawing a jelly bean from section z' is just exp(z') …
Thoughts on Scales
Currently in the code we have 4 different types of 'scale' variables:
- Computational scales which obey expected relationships (VelScale=lScale / TimeScale)
- lScale - length scale
- mScale - mass scale
- rScale - density (rho) scale
- VelScale - velocity scale
- TimeScale - time scale
- BScale - magnetic field scale (Gauss)
- pScale - energy density (pressure) scale
- Computationally scaled constants
- ScaleGrav - G in computational units
- ScaleC - speed of light in computational units
- Microphysical scalings
- TempScale = pScale/rScale * k_{B}/ XMU
- nScale = rScale / XMU
- Inverse computational scales used in source routines
- eDotScale = 1d0 / (pScale/TimeScale)
- And then inverse scales that involve microphysical scalings
- ScaleCool = 1d0 / (pScale/TimeScale/nScale^{2})
- ChemScale = 1d0 / (nScale/TimeScale)
A computational length, density, pressure, velocity etc… derive meaning through the euler equations - however computational temperature, and computational number density are meaningless. Only the corresponding physical quantity has any meaning. This is why the microphysical scalings are not consistent with the other scales (nScale /= lScale^{-3}) but instead are designed to give quick ways of getting physical values of Temperature and number density from computational mass density and pressure (provided XMU is the mean molecular weight mmw)
- n_{p} = rho_{c}* nScale → nScale = n_{p}/rho_{c}=n_{p}/rho_{p}*rScale = XMU*rScale
- T_{p} = P_{c}/rho_{c} * TempScale → TempScale = rho_{c} * T_{p}/P_{c}=pScale*n_{P}*XMU/rScale * T_{p}/n_{p}*k_{B}*T_{p} = pScale*XMU/rScale/k_{B} = pScale/nScale/k_{B}
The problem is further complicated when XMU can vary (as occurs with multi-species gas) Then the deviation of the local value (xmu) from some fiducial value (XMU) must be added to calculate the actual number density or Temperature
- n_{p} = rho_{c}*nScale*XMU/xmu
- T_{p} = P_{c}/rho_{c} * TempScale * xmu/XMU = P_{c}/rho_{c} * mTempScale
And now nScale can have meaning if we store the computational number density in tracer fields - though it is still not consistent with lscale^{-3}
- n_{p} = nScale * n_{c}
And if we are tracking computational number density, there is a consistency relation between computational mass density and computational number density:
- rho_{c} = n_{p} / nScale / XMU * xmu = n_{c} * xmu/XMU
Which allows us to write an expression for a modified TempScale
- mTempScale=TempScale*xmu/XMU = TempScale*rho_{c}/n_{c}
In any event I vote we
- leave the true computational scales, as well as the microphysical scales (nScale and TempScale) as is
- rename the inverse scales to have the Scale in front (ScaleCool, !ScaleEDot, ScaleChem)
- and then to avoid confusion change ScaleGrav to ScaledGrav and !ScaleC to !ScaledC
Meeting Update
Finished the documentation on FluxLimitedDiffusion in the Developer's Guide
New movie of Magnetized Colliding Flows
Meeting update 3/17/2013
- Worked with Eddie on the multi-species EOS and the effective gamma. Need to distinguish between tracked and untracked species.
- Looked over Shule's implicit diffusion solver and read up on Mixed-Frame Flux-Limited Diffusion. I should have it working by weeks end. Also began documenting the solver in the developer's guide.
- Also came across some old documents dealing with Radiation Transport I had prepared 'back in the day'.. See the attachments
All about gamma
In a multi-species gas - one can have different gamma's associated with different species. Combining them to give an effective gamma is not entirely straightforward.
One might expect a particle number weighted gamma might be appropriate
but it turns out that the value that should be weighted is the degrees of freedom
There is a relation between the degrees of freedom and the adiabatic index
so an effective gamma can be derived by taking
Meeting Update 2/25/2013
- astrobear alias - currently points to clover and we back up to botwin
- Any reason not to move wiki server and repositories to botwin and backup elsewhere?
- Need to find 'clover' references and rename to 'astrobear'
- Updated bov2jpeg program and documented it on the wiki
- Updated ImageObject and ProjectionObjects pages
- Implemented non equilibirum DM and improved cooling/chemistry/EOS module structure.
- Having trouble with Christina's run on BlueStreak
Out of memory trying to allocate 482088528 bytes
- Tried restarting with 2x processors but same error
- Will try re-restarting to make sure nothing changed on system end
- Will try restarting with 2x nodes but same number of cores (8 cores/node)
- This week:
- Test new cooling routines
- Implement HelmHoltz EOS
- work on FLD
True Color images in visit
I was thinking about Eddie's emission plots and how where there is h-alpha and sulfur 2 emission, the plot is brown instead of bright yellow - and this has to do with the way visit adds transparent pseudocolors… Basically it adds them as though they were pigments - not sources of light. Red+Green dye = brown, Red+Green light = yellow..
I thought about doing the colors in reverse - basically using anti-red (cyan) and anti-green (magenta) to get anti-yellow (blue), and then inverting the final image to get red, green and yellow - however, visit doesn't add transparent plots in a commutative fashion (the order of the semi-transparent images matters). Then I discovered visit's new true color plot which allows you to specify a rgb vector field and then plot the resulting image.
So if you create an expression for H-alpha emission and S-II emission that is scaled between 0 and 255,
H_sc=max(0, min(255, (H-Hmin)/(Hmax-Hmin)*256))) S_sc=max(0, min(255, (S-Smin)/(Smax-Smin)*256)))
then you can create a vector
{Halpha_sc, 0, SII_sc}
that you can then plot with a truecolor plot.
(I'm not sure why green is 3rd and not 2nd like 'rgb' - and I couldn't get blue to plot at all)
The following was produced by having the red value go from 0 to 255 across the left half and then stay at 255 on the right half, and the opposite for green. That way the center is at 255, 0, 255 = bright yellow….
More NEQ stuff...
Splitting DM into Hydrogen-Helium-Metal-Molecular processes
- So I spent a fair amount of time breaking down the curve in DM into 12 various parts that we use in non-equilibrium cooling. Here is a plot showing all of the different processes:
Now defining:
- x = HII/nH
- xHe = HeII/nHe
- yHe = HeIII/nHe
- y = ne/nH
each process can be identified with proper weight factors (before being multiplied by nH^{2})
Process | weight factors * f(Temp) |
e Z excitation and bremsstrahlung | y |
e H excitation | (1-x)*y |
e He excitation | (1-xHe-yHe)*y |
e HeII excitation | xHe*y |
H ionization | (1-x)*y |
HII recombination | x*y |
He ionization | (1-xHe-yHe)*y |
HeII recombination | xHe * y |
HeII ionization | xHe * y |
HeIII recombination | yHe * y |
Molecular cooling | 1 |
H Z excitation | 1 |
And here is the reconstructed cooling function using equilibrium values from the code plotted alongside the total equilibrium curves extracted from the DM paper:
Other than the rise at 10^{4} due to helium excitation, the majority of the cooling comes from metals excitation by electrons (especially Oxygen at 10^{5}). Of course the DM curve assumes some metal fractions (which are different from Pat's) so merging these two might be difficult. Though it would be possible to modify the DM curve to use Pat's metal fractions…
Breaking up the metal cooling
The form of the cooling function is dependent on the abundances of each element, and each of the sources used to construct the various tables use slightly different abundances:
CT | DM | PH | GT | RC | |
H | 12.00 | 12.00 | 12.00 | 12.00 | 12.00 |
He | 11.20 | 10.92 | 10.93 | ||
C | 8.60 | 8.57 | 8.52 | 8.52 | |
N | 8.04 | 7.94 | 7.96 | 7.96 | |
O | 8.95 | 8.64 | 8.82 | 8.82 | |
Ne | 8.70 | 7.41 | 7.92 | ||
Mg | 7.43 | 7.42 | 7.42 | ||
Si | 7.50 | 7.51 | 7.52 | 7.52 | |
S | 7.30 | 7.15 | 7.20 | 7.20 | |
Fe | 7.51 | 7.60 | 7.60 | ||
Ar | 6.80 | ||||
Ca | 6.30 | ||||
Ni | 6.30 |
DM Curve
Low Temperatures
The low temperature DM curve is dominated by carbon, silicon, and iron
The left plot shows the functional form for each of the lines considered from DM, while the right plot shows the total of all of the lines for each element along with data points extracted from a similar curve in the DM paper.
High Temperature
I also looked at how the abundance affects the cooling at high temperature. The DM curves are only for low temperature, although the completed curve uses cooling functions from Cox and Tucker. So, first I Combined the DM curves at low temperature, with the Cox and Tucker metal cooling curves (though using Gould and Thakur's hydrogen and helium excitation curves).
and then compared the resulting cooling function with total DM curve
The largest errors are around T=1e6 where the dm curve is weaker by a factor of ~2. Tucker and Koren give additional cooling functions for iron at high temperatures that I have not included…
Pat's Data
The data from Pat is a 3D-table
while the DM metal curve is only a function of TSo I divided Pat's table by
to get
so that
since Helium is not ionized in the range of Pat's tables…I then wanted to construct curves at constant nH so I took slices through the 3D data set at
where n_{Hi} = [1e2,1e3,1e4]. I then took lineouts at xi = [.01, .1, and 1] to construct 9 curves which I then plotted against the DM metal curve.They all agree fairly well at high temperatures and deviate at low temperatures when the ionization fraction is high (which is significantly out of equilibrium). But I took an average of these 9 curves to get a single metal cooling curve from Pat's table.
I also looked at the DM cooling function at low temperatures but using Pat's abundances.
The cooling strength from Pat's tables is still much less than the DM curve. But the DM curve appears to have several iron transitions - and if all of them are included, it does not agree with Pat's tables. If I play around with the various iron and carbon contributions I can get a curve that agrees at low temperatures with Pat's tables.
I then combined that curve, with Pat's 1D average curve, and the original DM curve to construct a new metal cooling curve that is more consistent with Pat's cooling tables.
Comparing the cooling curves
This plot shows the DM cooling curves built out of Pat's abundances, the DM abundances, as well as the total metal curve extracted from the DM graph. Also are shown metal cooling functions from Pat's tables at ionizations of .01, .1, and 1 and densities of 1e2, 1e3, and 1e4, along with their average value, and a complete cooling curve that blends the DM graph with Pat's data.
And here is a close up of the range covering Pat's tables
And here is a fun image showing the shape of Pat's metal cooling functions for each species.
As well as the total cooling function
and a breakdown of which species dominates the cooling
I also did a breakdown of Pat's tables where I averaged in X and ne to get functions of T for each species
As well as taking the maximum value in X and ne
Since the iron cooling seems to stick out - I looked in more detail at the DM iron cooling processes:
e+Fe^{+}(a^{6}D_{9/2}) → e+Fe^{+}(a^{6}D_{7/2}) | Seaton 1955 |
e+Fe^{+}(a^{6}D_{9/2}) → e+Fe^{+}(a^{6}D_{5/2}) | Seaton 1955 |
e+Fe^{+}(a^{6}D) → e+Fe^{+}(a^{4}F_{7/2}) | Pottasch 1968 |
e+Fe^{+}(a^{6}D) → e+Fe^{+}(a^{4}F_{9/2}) | Pottasch 1968 |
Oxygen and Nitrogen cooling contributions
Since Oxygen and Nitrogen ionization states are tied to hydrogen via charge exchange, I looked at dividing out 'xi' dependence to make all of the cooling graphs 2D instead of 3D.
Here is a comparison of the original cooling functions for OI and OII with the reconstructed version built from a 2D table.
Also talked with Pat about extending tables to higher temperatures and what new species would need to be tracked.
About the attachments
See the attachments for excitation cooling of hydrogen and helium as well as metals. The DMCooling spreadsheet has all of the calculations, while the DMTables has the various temperature dependent functions for the original DM. And the DMZ Blend has the modified metal cooling function that is more consistent with ZCooling (Pat's tables). And the PHCooling.bov/dat files can be loaded into visit to see the function PH'(ne, x, T)
Working on modifying neqcooling stuff to work with equation of state module
I added a new EoS called EOS_MULTISPECIES_GAS which uses an array of species mass, free electron number, and gamma, to calculate pressure, temperature, etc…
And I've also started to look through the neqcooling routines to see how the code handles some chemistry stuff… Apparently the last version of AstroBEAR did not support any molecular hydrogen - just hydrogen and helium ionizations…
I tried looking back in the repository for molecular hydrogen routines and found the 1st revision checked into the repository had
OIcool=OI_cool(nvec(iH2),nvec(iH),Temp) H2cool=H2_cool(nvec(iH2),nvec(iH),Temp) dmcool = DMCoolingRate(Temp)*ne*(nvec(iH)+nvec(iHII)) H2diss = H2_diss(nvec(iH2), nvec(iH), nvec(iHe), ne, Temp) H2term = -H2diss + H2_dust_Recomb_table(nvec(iH)+nvec(iHII),nvec(iH),Temp) hion = H_ioniz_table(nvec(iH),ne,Temp) hrec = H_recomb_table(nvec(iHII),ne,Temp) heion = He_ioniz_table(nvec(iHe),ne,Temp) herec = He_recomb_table(nvec(iHeII),ne,Temp) dqdt(iE)=-(OICool+H2Cool+dmcool+BindH2*H2diss+IonH*hion+IonHe*heion) dqdt(iH2) = H2term dqdt(iHII) = (hion-hrec) dqdt(iH) = (-(hion-hrec) - 2.d0*H2term) dqdt(iHe) = (herec-heion) dqdt(iHeII) = -dqdt(iHe)
And there were two sections commented out. The first one had to do with a critical density for molecular H2 cooling
CDH2=CriticalDenH2(nvec(iH),Temp) ! function used to smothly conect high HI and low HI density cooling rates CD=one/(one+CDH2) IF(CD<0.01) THEN H2cool=H2_critical_cool_table(nvec(iH2),Temp)*(one-CD) ELSE IF(CD<0.99) THEN H2cool=H2_cool(nvec(iH2),nvec(iH),Temp)*CD+H2_critical_cool_table(nvec(iH2),Temp)*(one-CD) ELSE H2cool=H2_cool(nvec(iH2),nvec(iH),Temp) END IF
and the other had to do with dust cooling
dustcool=dust_cool_table(nHneuc-nvec(iH2),Temp) ! dust cooling due to H sticking dqdt(iE)=-(OICool+H2Cool+dmcool+BindH2*H2diss+IonH*hion+IonHe*heion+dustcool)
Then in 2009 Bobby made a few changes
He added Helium double ionization though forgot to modify the source term for HeII to account for ionization / recombination
heIIion = HeII_ioniz_table(nvec(iHeII),ne,Temp) heIIIrec = HeIII_recomb_table(nvec(iHeIII),ne,Temp) dqdt(iHeIII) = (heIIion-heIIIrec) dqdt(iHeII) = -dqdt(iHe)
He also (accidentally?) removed molecular hydrogen reactions
He then modified the dm cooling contribution from
DMCoolingRate(Temp)*ne*(nvec(iH)+nvec(iHII))
to
DmCoolingRate(Temp)*nvec(iHII) * (nvec(iH)+nvec(iHii))
which amounts to assuming that the number of free electrons is equal to the number of ionized hydrogen atoms… Which isn't true if there is helium around to ionize…
He also added hydrogen excitation cooling
Hex_cool=ne*nvec(iH)*HexCoolingRate(Temp)
and recombination energy losses
Rec_cool=(Boltzmann*Temp)*hrec + (Boltzmann*Temp)*herec
While the average energy of the recombining electron is 3/2 k_{B} T, I vaguely remember Pat saying something about integrating over the cross sections and the velocity distribution etc…. and that the energy loss is just k_{B} T
And finally he added BBC Cooling and changed the name of the dmcooling table file to DMcooling.tab instead of cooling.tab
Then 1 month later Jacob turned off HeII ionization and HeIII recombination completely which seems reasonable as it was only half-implemented…
And then Bobby implemented a Van der Waals equation of state…
Cooling Processes
Process | notes | Dalgarno McCray | BBC | II Cool | NeqCooling | ZCooling |
Hydrogen Ionization | X | X | ||||
Hydrogen Recombination | X | X | ||||
Helium Ionization | X | X | ||||
Helium Recombination | X | X | ||||
gas-grain collisions | ||||||
H2 dissociation | ||||||
H2 formation | ||||||
H2 rovibrational | ||||||
OH rovibrational | ||||||
H2O rovibrational | ||||||
CO rovibrational | ||||||
atomic collisions | ||||||
brehsstrahlung (free-free) | ||||||
hydrogen-impact excitation of metals | Important at low x (and low T) | X | X (weighted by x) | X (weighted by x) | ||
electron-impact excitation of metals | = | X | X | |||
electron-impact excitation of hydrogen | = | 2X | 2X | |||
photo-electric heating from small grains and PAH's | ||||||
cosmic rays | ||||||
soft x-rays |
Meeting Update 2/4/2013
- Magnetized Colliding Flow Update
movie | movie |
- Looked at riemann solvers for general EOS johannjc01312013
- Van der Waals EOS implemented and tested johannjc02012013
- New equations of state possible - just need to provide functions (or tables) for:
- Pressure given internal energy and density + tracers
- Sound Speed given internal energy and density + tracers
- Internal energy given pressure and density + tracers
- Temperature given internal energy and density + tracers (technically not needed for euler equations)
- Does this automatically go into 'Golden Version' since collaborators are interested in using it? or is there a centralized development branch we should be pushing/pulling into?
- Organized existing documentation into Developer's Guide and various team pages.
- Updated the code to be doxygen compatible and enabled the doxygen tab
- Plans for this week
- Write and implement Helmholtz EOS
- Sesame tables? - maybe someone else wants to take this on?
- FLD solver
Van der Waals EOS
So for any new equation of state, we need to be able at the very least to calculate:
If we want to know the actual temperature - we also will need a thermal EOS
and if we want to do any roe averaging we will need
As an example, I looked at a Van der Waals gas
So for a Van der Waals gas the pressure takes the form
However to solve the euler equations we need a caloric equation of stateSo we need the internal energy density (which can't be derived from the pressure function alone)
where f is the degrees of freedom.This allows us to write our caloric EOS
as well as solve for the pressure given the internal energy density
so we can calculate the specific enthalpy
and then use the specific enthalpy to calculate the sound speed using various thermodynamic identities…
which when
reverts to
and when
reverts to
which gives damping sound waves when
and when
reverts to
and finally we need to invert the equation for the specific enthalpy
So we ran the basic shock tube tests comparing ideal gas and VanDerWaals gas.
Note that as
then pressure should begin to dramatically increase with density due to limited volume.And if
then attractive Van der Waals forces begin to significantly modify the pressure.The parameters in cgs were
which corresponds to a mean molecular weight of 36 amu |
Note since
And
so molecular attraction is negligible.In what follows, the red line is Van der Waals, and the blue line is ideal gas.
If we match the temperatures in the left and right states we get
and if we match the pressures
and if we match the internal energy
The energy match looks identical to the temperature match because the " a " constant is so small (although the bounds are different on the plots).
Arbitrary EOS update
I have successfully moved nearly all references to gamma to the EOS module - and all uses of pressure get the pressure from the press functions which can be modified to handle arbitrary EOS.
Now I just need to modify the Riemann solvers to handle an arbitrary EOS - or build a new one.
- For hydrodynamics, the riemann problem has three waves
- left shock/rarefaction
- entropy/contact
- right shock/rarefaction
- Neither the pressure nor velocity changes across the entropy/contact wave, so we can define a p* and u* that connects to the left state across a shock/rarefaction and to the right state across a shock/rarefaction. We can calculate dP*/du* so that newton rhapson iteration can be used.
For arbitrary EOS, there are now two difficulties.
- First, since the density can change across the contact, we know have two more unknowns gamma*_{L} and gamma*_{R} and two more equations which must go into the solution.
- And we no longer can easily calculate dP*/du* so we must resort to the secant method.
The two new equations are
where
Unfortunately the above solution is designed around the exact Riemann solver - which does not work with MHD. What about approximate Riemann solvers like the HLL?
The HLL method does not require the solution for the pressure in the star region - but just uses the conservation law and the assumption of a constant uniform state between the left and right waves.
Then integrating around the entire volume we have
And then integrating around the right half we have
Or the left half
which we can solve F*
so the only modification would be in estimating the wave speeds - but this can be done using the sound speed (or fast speed) in the left and right states - so only modification would be to calculation of sound speed. Seems very straight forward.
What about the HLLC? Once we have U* from the HLL we can estimate the speed of the contact
and then the jump conditions can be solved for p* from either the right or left as:
which then give everything else - from jump conditions. No EOS assumed
What about HLLD? Same basic principle - everything derives from jump conditions… Should be trivial to modify all of these HLL type solvers to work with arbitrary EOS.
Meeting Update 1/28/2013
- Magnetized Colliding Flow run on BlueStreak running well after successful restart
- Worked on writing your own problem module and updated information on using magnetic fields
- Consolidating operations dealing with equation of state to EOS module - updated riemann solvers to use same ordering as Info%q (rho, p, E) instead of (rho, E, p) so that many of the same routines don't need to be duplicated.
- Plans for this week
- Scaling tests on Kraken for proposal
- Work on developer's guide as well as other wiki projects
- implement tabulated EOS
- Start on FLD solver
Meeting Update 1/21/2013
Worked with Christine in setting up her problem module.
- Uses AmbientObjects, DiskObjects, RefinementObjects, ProjectionObjects, ProfileObjects, and implements the stream the ol' fashioned way - though we may want to make a StreamObjects or a JetObjects to do 'jets' from a boundary.
- Ran into problems with the boundaries due to cooling time scales being short and DM Cooling has no equilibrium temps.
- Solution was to have very high density contrast between disk (1e-7 g/cc) and outer edge of ambient (1e-20 g/cc).
- Also wanted to have hydrostatic equilibrium without requiring temperature extremes - so set up isothermal solution…
- First we added a constant pressure to p_inf to get the desired pressure at the disk edge - but since the density changes, this makes the ambient no longer isothermal.
- IF we try to solve for to give pressure matching at the disk edge we end up with a transcendental equation
- If instead we fix T or and solve for we get an easier equation to solve.
Meeting update 1/13/2013
PPM now converges to 3rd order using the new less-restrictive limiters.
Still getting oscillations in BrioWu Shock tube
Not sure what is going on - though probably not a bug with PPM - as that would ruin convergence.
Christina's run
Went 20 or so frames during the 24 hour reservation on bluestreak without difficulties. System has been down - so I've been unable to restart and run further.
FLD
Still waiting on working implicit thermal diff solver…
Things to do:
- Try to nail down source of oscillations in BrioWu shock tube
- Continue to work on documentation issues (Props to Erica for New Users Guide)
- Write up basic differences between directional split schemes (ie MUSCL-Hancock) with unsplit schemes (ie CTU) and the use of the Sweep algorithm to do unsplit updates for developer's guide and to go over briefly at next week's group meeting?
- Work with Christine to setup problem module
A few wiki documentation things to sort out
- contact us link on Download Page does not work if not logged in
- Instructions for installing AstroBEAR discuss using mercurial but don't mention where to check out the official 'repo' – which would of course only work for local folks anyways…
- AstroBearStandardOut needs to be updated to describe elliptic time usage. It would be helpful to also show some common errors that are seen.
- ModulesOnAstroBear
does not discuss the ProblemBeforeGlobalStep subroutineSimulation data section should not encourage the use of numerical indices for things like px, or Energy as these will be dependent on equation of state, lMHD, nTracers, etc…Dimensions section should be updated since 1D problems are possibleDimensions section should also encourage the use of CellPos function to get (x,y,z from i,j,k)Units and scaling section should mention scales.data file that is dumped by codeAlso worth mentioning that value of info%q should technically be average of physical solution over cell… Motivation for sub-cycling/smoothing etc…Also perhaps setting up a clump would be a good example - since you could then show how to rewrite the problem module using the ambient and clumps objects.Units and scaling section should not encourage the changing of physical scalesInitializing a Grid says that each infodef has a %dx which is not trueInitializing a grid has an incorrect expression for the number of ghost zonesFlagging cells for refinement implies that the errflag array extends into ghost zones – which it doesn'tAdditionalPhysics should probably be a separate pageAdditionalPhysics says that MaintainAuxArrays needs to be set to true which is no longer correct- AdditionalPhysics should discuss equation of state options as well.
Aux fields only need to be initialized if lMaintainAuxArrays is true (in 1D MHD it is false)Also an example that shows how to calculate the location of aux field entries would be good – as well as an example showing how to initialize a divergence free auxfield using a vector potential.- It would be nice to stop storing the Cooling Objects in the chombo file – which would simplify the problem module.
- Not sure if selfgravity still has a problem with a uniform density field… this should be checked – and it would only be a problem with peridiodic BC's for the potential – also a discussion (or link to a discussion) about what needs to bet set in global.data for self-gravity to work.
- I believe you can have sink particles without self-gravity (ie just point gravity)
- This should also be updated with viscosity/resistivity/thermal conduction stuff…
- AstroBearObjects
- Disks need to be added to Initial conditions documentation
- VectorPerturbations needs to be added to the Sub-Objects and documented
- All of the individual object pages need to be (created) or updated to reflect changes made by Ivan and others
- ScramblerIntro
- Need to decide what to do about refinement objects etc…
- How to setup a properly posed problem should probably be linked from somewhere else…
- DataFileTutorial
- io.data, communication.data, modules.data, and process.data are depracated
Data File Tutorials should be moved inside of ScramblerIntro and the 'Getting Started with AstroBEAR' should be renamed to 'AstroBEAR user's guide' or something… and linked below the Download link on the main wiki page.Then Tutorials will have information on things besides AstroBEAR - though the visit tutorial could be linked from the visualization page.The link to Code on the main wiki page is essentially a link to a developers guideMeeting pages are depracatedCollaborators & Projects, Publications, and Image Gallery probably don't belong under 'NewUsers'
Meeting Update 01/08/2013
- Things I worked on over break…
- MHD convergence - done
- PPM convergence - almost
- Things to work on this week…
- FLD solver
- Documentation ???
- Developer's manual ???
- Christina's MHD run got to frame 3 in 24 hours on 4000 cpus with refining up to 4 levels. Frame 3 was 11 GB! Resubmitted with 3 levels of AMR except where jeans criterion triggers 4 levels - and it has been in queue for 6 days.
Meeting Recap 12/17/2012
- Spent some time working on Christina's setup and improving refinement criteria. 'I' shape due to density gradients produced due to bending magnetic field lines and strong shear. In 3D this I is a cylinder. So I modified module to use three refinement objects.
- First maintain refinement to level 3 within a 2 pc wide region around the initial interface (nested inside of coarser patches)
- Next allow for density gradient based refinement to level 3 within a 10 pc region centered on the initial interface (nested inside of coarser patches)
- And finally allow for jeans length based refinement (64 cells per jeans length) all the way to level 4 so that collapsing structures will be better resolved.
- 3D MHD with self-gravity on bluestreak runs at 1475 cell updates/second/core
- On 32 bluestreak cores, at 32^{3} + 3 it was spending 76% of the time advancing, 23% solving elliptic equations, and < 1% everything else… and was slotted to finish in 8 days (30 Myr).
- In principle if everything scales up - the initial resolution of 128^{3}+4 could require between 256 and 4096 times as many resources depending on the filling fraction of the highest level…
- At 256 this would give 16 days on ¼ of bluestreak (4000) which seems doable? This is 1.5 Mhrs of CPU time on bluestreak - or probably .4 Mhrs on kraken
- Things to work on over break…
- MHD convergence - fast waves now converge - slow waves still do not.
- PPM convergence
- FLD solver
- Documentation ???
- Developer's manual ???
- Bash script for running parameter studies…
interporders=( 1 2 3 ) Gmxs=( 100 200 400 ) Waves=( 1 2 3 4 5 6 7 ) mkdir out mkdir testdata for (( i=0;i<${#interporders[@]};i++)); do interporder=${interporders[i]} ../../scripts/subst.s InterpOrder $interporder solver.data for (( j=0;j<${#Waves[@]};j++)); do wave=${Waves[j]} ../../scripts/subst.s WaveFamily $wave problem.data for (( m=0;m<${#Gmxs[@]};m++)); do Gmx=${Gmxs[m]} ../../scripts/subst.s Gmx $Gmx,1,1 global.data mpirun -n 1 ../../astrobear > output.out grep dx output.out | awk '{print $3, $4}' >> scaling_${interporder}_${j}.curve mv out/data.curve data_${interporder}_${j}_${Gmx}.curve cat output.out done done done
Meeting update 11/26/2012
- Getting ready for trip to Berkeley/Caltech
- Christina's setup takes 2 months on 4096 bluestreak cores. Scaling is ok - so likely a combination of
- higher filling fractions
- MHD
- Bluestreak speed issues.
Can't look at mesh since didn't get first frame. Another job with a higher frame rate has been in the queue since Saturday morning.
Meeting Update 11/19/2012
- Looked at BlueStreak performance
- Made new figures for paper
- Working on Christina's run #267
- Preparing talk for Berkeley on 12/7…
Performance comparison between BlueStreak and Grass
I ran a fixed grid colliding flow setup - with self gravity and cooling etc… for a few time steps to compare the performance of bluestreak when using 1, 8, or 16 cores on a single node as well as a single core on grass. I also compared different levels of optimization (-O3 and -O4). I tried -O5 but the code dies soon after outputting the first frame. I have tried compiling with traceback enabled to find out where the code is dying - but I don't seem able to get any useful information from the core dumps. I'm not sure what is going wrong there.
The figure below shows the cell updates/core/second over the first 10 timesteps for various numbers of cores on bluestreak - and for optimization levels 3 and 4 - as well as a single core run on grass with level 3 optimization (1_3g).
In any event the performance is about what you would expect. The code seems to slow down after a few time steps due to either CPU temperature - or perhaps cooling source terms… Optimization does not make much of a difference < 1% (at least going from O3 to O4). Also in going from 1 to 16 cores per node - the efficiency only drops around 10% which is surprisingly good. Overall the speed per core is about 25% that of grass. Given that grass has a 2.5 GHz processor and that each core on BlueStreak is supposedly running at 1.6 GHz - this is a little surprising. A single bluestreak core appears to be performing at 39% of what you would expect taking into account CPU speed.
Meeting Update 11/5
- Worked on Spectral Prolongation
- Using BlueStreak for post-processing of colliding flow runs
- Created ticket #264 for parallel IO
- Setting Up Christina's MHD Colliding FLow run on BlueStreak
Spectral prolongation
I've written a routine that can take a coarse grid, do a fourier transform, extrapolate in fourier space, and tranform back to real space resulting in a spectrally prolongated fine mesh. This could be useful for avoiding excess power in AMR spectra.
The fourier prolongation in 1D is fairly straightforward. The only tricky part is that every other fine mesh point does not coincide with the coarser mesh - as would be the case for finite element instead of finite volume. As a result each wave needs to be shifted by the offset between the 1st fine cell and the 1st coarse cell (which works out to be .5*dx_fine)
Here is the result without any shifting…
To shift each wave by .5*dx_fine we have to add a phase = k*dx to each wave. For each element in fourier space we have to determine the corresponding wave vector and then multiply by .5*dx*k.
I also added the ability to do linear prolongation as well as spectral prolongation of layouts for use with FFTs. Here are view of what the data looks like going into the FFT. Constant interpolation is lower left, than linear (monotone center) interpolaton lower right, parabolic interpolation is upper left, and spectral interpolation is upper right
And here are the spectra for those four methods
And here are the side by side view of the spectral prolongation and the fixed grid version.
I also looked at the colliding flow runs. Of course without windowing there ends up being ringing in the spectral prolongation due to the non-periodic nature of the region.
Windowing before doing the spectral prolongation reduces this dramatically
And here is the same window for a later frame
And a comparison between the resulting spectra with and without spectral prolongation
Meeting Update
- Found memory leak using valgrind
- Looked at Colliding Flow Spectra https://clover.pas.rochester.edu/trac/astrobear/wiki/CollidingFlows#Spectra
- AstroBEAR II Paper accepted to JComp Phys
- Started what could be user guide? https://clover.pas.rochester.edu/trac/astrobear/wiki/ScramblerIntro2?version=2
meeting summary
- worked on Martin's paper
- Hydrostatic Star Problem
- Applied for Miller fellowship
- Working on Spectra for colliding flow papers
installed mypage plugin from source
export PYTHONPATH=/data/trac/astrobear/plugins:$PYTHONPATH; easy_install --install-dir /data/trac/astrobear/plugins tracmypageplugin/0.11/
New buildproblem script and current tests plugin
I modified the buildproblem script and currenttests plugin to use your username in the various paths. This will prevent multiple users from stepping on each other's toes. To see your current test output you need to create a page (or somewhere on your personal page) that has the following text…
(where johannjc is replaced with your username)
Of course this won't work until you've updated your buildproblem script but that should be in the main repo soon…
meeting update
Worked with Shule and Ivan on code implementation and prepped for journal club. Left for wedding on Friday
Meeting Summary
- Completed test modules and documentation for:
- Updated Matt's tests macro plugin to show links to test descriptions (see AstroBearTesting)
Instructions for taking a problem module and turning it into a test module
- Link and build your problem module.
- Modifying the data files to follow the format of those in Template
- Adjust the number of frames to 1 and the final time to be something small so that the run completes in under 1 or 2 minutes (optimized)
- make a ref directory parallel to the out directory and copy over chombo00001.hdf
- make a logs directory
- make an images directory
- Run bear2fix and select operations to make a nice image that can be used to verify the solution is 'correct'
- run ../../bear2fix/ps2png.pl out/pgout00001.ps
- mv out/pgout00001.png images/ref.png
- Now rerun the same code with a different number of processors and pipe the output to logs/testlog
- Rerun bear2fix in batch mode to produce another image for this run
- run ../../bear2fix/ps2png.pl out/pgout00001.ps
- mv out/pgout00001.png images/sim.png
- mv bear2fix.data bear2fix.data.img
- Run bear2fix and select operation 10
- mv bear2fix.data bear2fix.data.test
- mv out/chombo00001.hdf logs/sim.hdf
- cp max_errors.data to the ref directory and mv max_errors.data mean_errors.data and max_rel_errors.data to the logs directory
- add bear2fix.data.img and bear2fix.data.test to the mercurial repository (hg add bear2fix.data.img bear2fix.data.test)
- Now navigate back up to the code base directory and add your problem module to the list of problems defined near the top of /scrambler/buildproblem. Add it to the front of the list for faster debugging!
- Commit the new file additions and the modification of buildproblem to the repo.
- Now you will want to use svn to check out a new version of the test repo. svn co file:///cloverdata/repositories/tests ~/tests
- Create a directory for your problem module in the test repo parallel to Bondi, etc…
- Navigate into your newly created directory and move the logs, images, and ref directory from your code directory
- Now add your test problem directory to the svn repo using svn add MolecularCloudFormation (or whatever yours is called)
- Now we can commit those additions to the svn repo. svn ci -m 'added my problem module'
- Now you need to let the trac test macro know about the new module directory. So log on to clover and go to /data/tests and run svn update and then run trac-admin /data/trac/astrobear changeset added tests N where N is the revision number you just updated to.
- At this point you should be able to go to the test page on the wiki and see your sample test results https://clover.pas.rochester.edu/trac/astrobear/wiki/AstroBearTesting
- Now go back to your code directory and run buildproblem -np 8 -t to start the problem testing
- If this passes then you are all set to check in your mercurial repo! If not then seek help
- Don't forget to add your problem module to https://clover.pas.rochester.edu/trac/astrobear/wiki/TestSuite
- And then create a page for your module using the TestTemplate template
Thoughts on Self Gravity, Phi, etc...
So…. I was looking at Erica's BE setup and decided to start advecting the stable sphere and ambient - to test the stability/accuracy of self-gravity… First I noticed that occasionally the value of PhiDot would be very large… This was because a root step would finish just before a frame dump - so the code would take a very small time step - this would cause noise or small fluctuations in the potential to generate very large time derivatives… If these large derivatives are then used for a normal size time step - they will produce drastic errors… This is because of the drastic extrapolation in time of the potential. One way to mitigate this is to avoid having to take very small time steps by anticipating such scenarios and taking two half-size steps instead of a full step and a tiny step.
(Errors end up being amplified by the ratio of time steps)
Fortunately we no longer use old time derivatives for subsequent steps (except to set the initial guess for the solution vector which could make converging a bit difficult for visit) - but it does make the value for PhiDot plotted in visit look way off…
The other problem is regridding… As the grid adapts - regions that previously were refined will continue to have a correct fine-grid solution while regions that were not previously refined will have prolongated solutions… If a cell goes from having prolongated to a calculated solution - the time derivative will have a component in it that is just because of changes in the discretization of the grid - which could then be used to set boundary conditions for finer cells etc… Here the solution is to calculate the potential after regridding - before advancing - and then calculate it again after the advance to get a time derivative for the potential that is solely due to gas motions and not regridding effects… It is a little more expensive (though hypre should converge quickly since the initial guess will be fairly accurate) but well worth it in terms of stability etc…
Meeting Update for 8/28/2012
- Worked on static AMR improvements. #242
- Found a bug with self gravity and sink particles.
- Updated my tickets report {7} to show due date.
- Noticed that data files explained section of wiki should be updated with golden release… Created an enhancement ticket #246 - but it is currently unassigned…
- Got PhD
Meeting Summary
- I've been working on getting ready for my defense - Tomorrow at 2:30…
- I've finshed the talk - just need to run through it and time it a couple times
Energy budget for non-accreting case
I looked at the total energies for the non-accreting case and everything looks good.
Total energy seems to flatten out after the first collapse… The 5% error could be due to the softening approximation for the gravitational energy I am using in visit… But the temperature on this run also looks good
Here we see the temperature at the 1st and last frame on a linear scale.
However, with accretion turned on a hot spot develops instantly within the accretion kernel before the disk has had time to collapse or heat… This is certainly a bug that I will look into this weekend.
It looks like the problem module was creating two point gravity objects - and that the additional point gravity object had no softening. This was causing the inside out collapse and the large heating at the center.
HydroStatic Isothermal Disk
First every distance is softened to avoid trouble at origin
Distance from origin:
Distance from rotation axis:
Disk radius:
Then angular velocity is set to be keplerian and balance force along the 's' direction
Then the density within the disk is set to be constant along the midplane
with a decaying exponential to give vertical hydrostatic equilibrium.
where the scale height is given by
At the outer rim of the disk when
Omega goes to zero and the profile changes slightly to
which is just the profile for an envelope about a point mass in hydrostatic equilibrium ie
Then for each point
where is the disk background density.Finally additionally smoothing is done in density and pressure between the disk solution and the original cell values at the outer radius using the disk%thickness to define this smoothing distance - although I did just notice a typo if one uses the hydrostatic disk without using an isothermal EOS…
s(iE)=(gamma7*Disk%pressure-q(iE))*f
should be replaced with
s(iE)=(gamma7*Disk%Pressure/Disk%Density*rho-q(iE))*f
so that if one uses a gamma of 1.0001 the disk is at the same temperature.
Solid type boundaries with accelerating source terms
So Erica was seeing high velocities and low densities near the boundaries even though the ghost zones were being stepped on with zero velocity and moderate densities every time step. So I setup a simple test problem to verify that this behavior was reasonable. I used a uniform gravity setup but with the top half of the grid being stepped on each time step with the initial values of density, pressure and 0 velocity. I also did a comparison with a reflecting boundary instead to see if it was much different - although there are 'wiggles' that appear in the reflecting case perhaps because of density protections? Here everything has a value of 1… rho, g, P, cs, etc…
Attached is a movie showing the behavior of the density and velocity fields for both cases.
Gource rendering of scrambler...
(The video shows all of the files in the scrambler repo - and who was working on which file when… - though drastically sped up
After a while the module directory begins to take on a life of it's own - and later on you can spot bear2fix - and the test modules… Also .f90 files are in green and data files in purple See if you can spot the module-objects directory
Don't forget to blast your favorite heavy metal tune in the background while watching…
And for those who reminisce about AstroBear 1.0
Self-gravity thoughts...
Thoughts about the potential having negative curvature
3D
So in 3D - with radial symmetry - you can have a potential that radially has a negative curvature (ie
) and still satisfy poisson's equation since but ratherSo the potential can go like r^{2} inside a uniform sphere and then switch to -1/r outside.
1D
For an infinite slab of uniform density - the potential inside goes like x^{2} and switches to just x^{1} outside… - so it has zero curvature outside as expected - and never has negative curvature anywhere.
Thoughts about the 1st derivative of the potential - or the Gravitational Field Strength |g|
Uniform density
Imagine a uniform density sphere in a very low density environment. As you move towards the sphere from infinity - the gravitational field increases because you are getting closer to all of the mass. Now as you enter the mass two competing effects happen. On the one hand you are getting closer to most of the mass, but at the same time you effectively stop feeling the rest of the surrounding mass.
Or in equation form…
Here the first term represents the fact that you are getting closer to the remaining interior mass - while the second term reflects the fact that the amount of matter inside is decreasing.
For a uniform density mass distribution,
so we haveSo the strengh of the gravitational field peaks at the outer edge of the cloud. This means that the outer material feels a stronger acceleration and falls in faster… This makes sense since if everything accelerated at the same rate - you would get material piling up at the center - and the uniform collapse is uniform. So the uniform collapse is essentially an outward in collapse
1/r mass distribution
If the mass distribution goes like 1/r something very interesting happens…
and we have
So every point feels the same inward acceleration. Matter will accumulate fastest at the center and you have an inside out collapse..
Point mass
The other extreme is to have a very centrally concentrated mass distribution (ie a point mass). In that case the second term is zero since as you get closer to a point mass - you don't see any less mass - so the strength of the field increases like
- and so the acceleration is fastest near the point particle - and you would get an inside out collapse.BE Sphere
A BE sphere is something in between a uniform density and a 1/r potential… And the inflection point will probably be where the density switches from a constant to 1/r - or it looks like at a non-dimensional radius of ~ 3… Of course we have completely neglected pressure gradients which are extremely important for a BE sphere…
BE sphere in a uniform density background
So now if we imagine a BE sphere surrounded by a uniform density ambient….
First the ambient effectively sees a point mass (provided that the BE sphere is much heavier than the ambient) - so it will want to collapse from the inside out. So a rarefaction should develop in the ambient. Except at the edge of the BE sphere - where the infalling material must pile up (or heat up and re-expand)… In any event this would add additional mass to the outer edge of the BE sphere and send a pressure wave through the cloud. The center of the BE sphere - must be completely ignorant of any additional mass piling up outside - so it will only see pressure waves. If the pressure wave is slow enough - or low enough amplitude - the central region should remain stable while the outside is free to undergo collapse…
Meeting Summary
- Have good initial scaling results on Kure up to 256 cores #202… Moving to Kraken and 10,000 cores…
- Page describing AstroBEAR's output
- Worked on threading development for scaling results and various tickets…
Made a page describing AstroBEAR's super-standard out..
https://clover.pas.rochester.edu/trac/astrobear/wiki/AstroBearStandardOut
Or if you want to navigate there yourself it is wiki⇒Tutorials⇒Getting Started with AstroBEAR⇒AstroBEAR 2.0 standard out
Some exciting new features!!!
Impact parameter for companion disk
If we modify simple BH accretion about a fixed object by applying an acceleration
that is perpendicular to the direction of the flow - then it is reasonable to expect the 'backflow' to be aimed towards a position somewhere in between the object's original position and current position because the material that is captured has been accelerated on average towards a retarted position. If we take the time scale for this capture as R_{Bondi}/v_{wind} we can calculate a distance the object will have moved due to acceleration as . If the acceleration is not perpendicular to the flow but at some angle then it is the projection of the acceleration perpendicular to the flow that matters…If we now imagine that we are in an intertial frame of reference that is comoving with the secondary's velocity at t_{0} - we have the secondary accelerating towards the primary at
and the projection relative to the incoming flow velocity as just . Given that acceleration of captured material occurs over a time scale we get that the offset should go as and should be radially outwards. This would imply that the backflow would return outside of the secondary and would then proceed into a prograde orbit about the secondary.Meeting Update
- Created a work sheet for calculating various run parameters…
X Linear size of entire base grid N Number of cpu's D Physical dimension of problem C Cell updates per cpu per second T Run time x Linear size of base grid on each cpu L number of levels or refinement R Refinement ratio F Filling ratios between AMR levels (assumed to be fixed across all levels) E Ratio of workloads between levels (E=F R^{(D+1)})
- In general for a given linear resolution X, there will be X cell updates per crossing time - and the number of cells to update will be X^{D} so the total work load for a problem goes as X^{(D+1)}. If a single cpu can update C cells/second, then we have C N T = X^{D+1}
- Now if we divide a domain that is X^{D} into N pieces to distribute, then each piece will have X^{D}/N cells and have a linear dimension x = (X^{D}/N)^{1/D} = X/N^{1/D}
- Now with weak scaling the number of cells per cpu 'x' is kept constant - so we have X^{D} ~ N. If we were actually doing simulations then the walltime would go as T=(X^{D+1})/N C ~ X ~ N^{1/D} so the 1,000 core simulation would take 10 times as long (assuming 3D) because there would be 10 times as many time steps though each time step would take the same amount of time since x is constant
- For strong scaling the goal is to keep X constant so x = X/N^{(1/D)} ~ 1/N^{(1/D)}. The number of time steps is unchanged but the time per step goes as x^{D} ~ 1/N so the over all wall time T ~ 1/N. (Of course the memory requirements also go as 1/N so a strong scaling run would take 1000x as much memory on 1 core as 1000 cores.
- For hybrid scaling the goal is to keep the wall time constant since this is what effectively sets the resolutions we run at. So X^{(D+1)}/N is kept constant so X^{D+1} ~ N and x ~ X/N^{1/D} ~ N^{1/D+1}/N^{1/D} ~ N^{1/[D(D+1)]} or x^{D} ~ N^{1/(D+1)} which is a fairly weak scaling - so hybrid scaling is very similar to weak scaling - but there is a slight decrease in the workload per cpu because in general more cpus → more time steps → shorter time per step.
- With hybrid scaling - the invariant is the wall-time which can be chosen intelligently to be 1 day or 1 week … But with strong or weak scaling we have to motivate a choice for x (in the case of weak) or X (in the case of strong). The best way to do this it so choose a target number of cpus and a wall-time. Then you can back out what X and x is for that number of cpus and that wall-time and use those values for the strong and weak scaling respectively.
- Finally with AMR - if we have a base grid with linear size X, then there will be F X^{D} cells marked for refinement - each of which will create R^{D} child cells that will need to take R substeps - so for each coarse level step there will be F X^{D} R^{D+1} level 1 steps and since the whole simulation will consist of X coarse steps, this will be a total of F X^{D+1} R^{D+1} level 1 steps. And for each level 1 step there will be (F X^{D} R^{D+1})^{2} level 2 steps… So for the entire simulation there will be X^{D+1} (1+F R^{D+1} + (F R^{D+1})^{2} + … + (F R^{D+1})^{L}) cell updates = X^{D+1} (1-E^{L})/(1-E) where E = F R^{D+1} So if we want to keep the wall-time constant as we add levels of AMR then we need to keep X^{D+1}(1-E^{L})/(1-E) a constant - so X ~ [(1-E)/(1-E^{L})]^{1/(D+1)}
- And we have the 'master equation' N C T = X^{D+1}(1-E^{L})/(1-E)
- So in summary:
Weak Scaling fixed #cells per processor Strong Scaling fixed #cells Hybrid Scaling fixed wall-time Hybrid-AMR Scaling fixed wall-time ( where is the filling fraction and is the depth )
- There is also the issue of grid sizes causing a slow down in computation. In general the overhead required to update a grid of size X reduces the efficiency to X^{D}/(X+M)^{D} where M is of order 4 or 5… If the hyperbolic advance uses up G zones then the extended ghost zone advance will reduce the efficiency to 2 X^{D}/((X+2G+M)^{D}+(X+M)^{D)… For X of 10 and G of 4 and M of 5 this gives 2000/(23}3^{+15}3)=13% efficiency - without taking into account any latency or communication bottlenecks…. This is just do to having small grids of size 10^{3} - combined with extended ghost zones… If the grids are 20^{3} this becomes 33%… But smaller grids means a lot more overhead - especially as the linear size of a typical grid approaches 10-20 cells
- Fixed a few threading bugs…
- Started working on thesis…
Remapping from Turbulence sims to 3D Cloud collapse sims
So I have two data sets from isotropic turbulence (isothermal and iicool) that consist of the density, velocity, (and energy for the IICool run) on a 512^{3} mesh that I am remapping to a the middle half of a level 3 region of a simulation with a 1024^{3} effective resolution…
Both should have been driven to the same mach number though the isothermal simulation has a higher temperature and XMU (20 K and 1.4 and a gamma of 1.0) where the IICool simulation has an equilibrium temp of 16.98 and an XMU of 1.27 and a gamma of 1.6666
The Isothermal run can be rescaled in density, length, or time - to a different density and/or temperature and/or size as long as the mach number is constant. This leaves us with two free parameters. We want a virialized cloud, and we want the peak density to be within the Jeans length at the 512 resolution…?
The current peak density is 40622 which is to be expected for mach 21.5 turbulence with a background density of 100? 100*21.5^{2}=46225
However at a temperature of 10 this gives a jeans length of only 0.192981859204160 pc. If we want to resolve everything by a jeans length of 32 cells… then this requires each cell be .012 pc or the cloud at 512^{2} resolution - be only 6 pc across.. We could lower the density, or raise the temperature …
So if we have a virialized cloud, the Jeans length at the mean density will be
and density enhancements of M^{2} will have Jeans lengths of order which means we will need a resolution of which allows for Mach 4 turbulence…
Options:
- Generate lower mach number turbulence and make region larger and less dense?
- Begin simulation with higher resolution… But it would be nice to get a chombo file complete with a solution for the potential to start with…
It would be nice to start with turbulence that wouldn't be strongly gravitating - at least locally … but would still be virialized… With IICool this should be possible - since virial equilibrium could be predominantly thermal support with lower mach number flows - that would trigger thermal instabilities and collapse would ensue…
Meeting Update
- Worked on cleaning up threading algorithm and improved the load balancing algorithm. Early results are quite encouraging. See Ticket #191
- Ran 512^{3} turbulence sims on Kraken using IICool and Isothermal. Velocity spectra seem very similar - but still analyzing other statistics…
- Working on mapping data from 512^{3} simulations onto gravitational collapse simulations. See blog:johannjc03062012b for more info..
- Managed to get code to compile and run with gfortran. Had to build hypre, openmpi, and pgplot with gfortran. They have been added to the modules but some of them are the new default on grass, alfalfa, clover… so check your modules and your paths. See ticket #195
- Updated CollidingFlows page
- Added CameraObjects to allow for projections from any vantage point.
- Added IO routines for layouts to make it easier to transfer data from one simulation to another.
Meeting Update
- Created SpectraObject page as well as PFFT and LayoutObjects
- Worked on creating isotropic turbulence of a given mach number for initial conditions for cloud collappsing cloud sims. See this blog
- Modified IICool routine to use analytic formula that can be adjusted at run time. See this blog
- Started isotropic turbulence runs on Kraken
- Added shapes to various ProcessingObjects. See this blog
Histogram funniness
So the smooth colliding flows run showed a peak at .4 as well as 1 in the mass-weighted histograms of the mixing ratio. (Red line in plot). It turns out that the expression for the mixing ratio was to blame. It was set to
to avoid nan's in the region where neither tracer was present instead ofThis means however in regions near the edges of the cylinder, where both flows collide and mix with the ambient, the mixing ratio will be ~ .5 or so… I replaced the expression with
and the second hump went away (Green line)I also remade the histograms using only the data within the colliding region and got the same result regardless of the expression… (dark and light blue lines)
First run on Kraken
Kraken Simulations
res | wall time | advances | restarts | cores | cell updates/s/cpu | notes |
256^{3} | 86410 | 4442 | 124 | 24 | 36934 | Isothermal |
128^{3} | 27475 | 12664 | 194 | 24 | 40893 | Isothermal |
Achieving a target kinetic energy by varying the forcing strength...
Energy dissipates in a crossing time which is
which means that the specific kinetic energy dissipation rate is . If we have a given accleration field , it should inject specific kinetic energy at a rate proportional toSo assume we have a system
where
and are unknown constantsand we want to adjust
so that or so thatwe could start by assuming that
and and set and then run until we reach a steady state energy. At that point we can evaluate and then adjust the accelerationThe kinetic energy could potentially be widely varying on short time scales - which makes it tricky to determine a steady state energy. Additionally any change in the forcing will take of order a dynamical time for the system to relax… So each adjustment to the forcing should be followed by 1 dynamical time of allowing the system to relax - followed by a half dynamical time of averaging. The initial state should probably be allowed to evolve for 2 dynamical times before averaging.
Meeting Summary
Spent most of the past week sorting out sources of asymmetry in the code…
The sources of asymmetry could be broken down into 4 categories…
- Non-symmetric addition [ a+b-c-d /= (a-c)+(b-d) ]
- Diffusion algorithm's calculation of the divergence
- Characteristic limiters that do matrix multiplication from slowest to fastest wave speed
- Calculation of back forces on particles
- Upwind type calculations that expect quantities to be < or > 0…
- Riemann fans that have an even number of states (odd number of waves) such as HLLC, HLLD, ExactRS
- MHD eigenvectors which depend on the sign of Bx or the 'polarization angle' [atan(Bz/By)] which break down when Bx = 0 or By=Bz=0
- Bugs in the code
- Extrapolated MHD physical boundaries not updating cell centered b fields using aux
- Clump object's default behavior to persist in boundaries
- Wind object's failure to zero out other fields - miscellaneous tracers etc…
- Machine optimization
- Apparently optimization can lead to a breakdown in accuracy (perhaps because it reorders additions/subtractions to be faster?). Fortunately this behavior can be corrected by the addition of -fp-model precise to the compiler options.
Developed co-rotating reference frame source terms
or if we assume that
then we havewhere
movie |
movie |
Thoughts on binary sims
So Eric and I were discussing the specific angular momenta in the gas yesterday and you can calculate the velocity of the wind at any given point (assuming the secondary does not influence it's motion) by solving for the characteristic that leaves from the surface of the primary with a velocity vector pointed at the position. It is difficult to solve exactly and some approximations must be made… If we assume that the wind velocity is much greater than the orbital velocity, and that the distance to a given position is much greater than the primary's orbital radius, then we can easily esimate the time delay (t-t_{0}) between the wind leaving the primary at t=t_{0} and arriving at a position x at time t as t-t_{0}=|x|/v_{wind}. We can then determine the location of the primary at this time, and then assuming that the primary's radius is << |x|, we can estimate the direction of the wind as being the same as the vector going from the primary towards |x| or v(x,t) x (x-x_{p}(t_{0})) = 0. Then to calculate the magnitude of the wind we know that it will be constrained between v_{w} +- v_{orbit} depending on the degree to which v_{p}(t_0) is aligned with (x-x_{p}(t_{0}))
See the following illustration
The upshot, is that the secondary tends to be in a region where the specific angular momenta from the primary is positive (which is retrograde)…
Here is the 'corrected' approximate solution plotted (here the orbit is clockwise and the z-axis if pointing out of the plane)
Click here for the original version. Note where |x| = 1 - there is no difference
And here is another more exact solution where the 'time delay' is not calculated based on the distance from the origin, but by the distance to the current location of the primary. The exact solution would require solving a transcendental equation to find the distance to the position of the primary at the retarted time…
and here is the calculations that went into the above approximation
And here is the exact solution calculated from the simulation
Here the yellow corresponds to positive specific angular momenta - and the system as a whole has negative angular momenta - so yellow ⇒ retrograde disk
Here is a few frames later after the secondary has begun to torque the wind…
The secondary should likely give a net torque that is prograde (so more blue then yellow) but positive specific angular momenta does not automatically form a disk since the disk has to be orbiting the system's center of mass as well as orbiting the secondary. So does the infalling material have to have a net specific angular momenta that is greater than the secondary to form a prograde disk?
And here is a zoomed out version of the approximate solution… Note the whole thing rotates with the same pattern speed - so just spin your monitor to see a movie
Also was thinking about what the flow looked like from the point of view of the secondary - I don't see any shear in the flow, but there is of course divergence. But in the co-rotating frame, the secondary sees a nearly uniform flow subject to a Coriolis force as well as a gravitational force… something like
So setting up the initial solution can be accomplished by the following:
First start with the assumption that the retarted time is the current time
DO
Calculate the displacement vector from the primary at the retarded time
Calculate the wind normal so that
It is possible to solve
for the unit vector
Calculate the wind velocity from the primary
Update the retarded time using the new distance and wind speed
END DO
The only problem occurs when there are multiple solutions for the retarded time…
This will occur once we reach distances of order
If we switch to a rotating frame that rotates counter to the orbit so the angular speed is
, then
and
so that
Meeting Update
Mostly working on turbulence stuff and developing routines for creating random fields using FFT's as well as spectra…
Created Layout objects which allow for a distributed fixed subgrid across all of the processors. There are also routines for loading data from the AMR mesh into the layout and extracting data from the layout into the AMR mesh - as well as transferring data from one layout to another
Created PFFTPlan objects which allow for easily performing parallel FFT's of data. It creates nDim layouts where in the 1st layout there are no divisions in x, in the second layout there are no divisions in y and so on…
Creating spectra processing objects which use PFFTPlans to load various fields from the AMR mesh into a local data array - then perform a parallel FFT - and then bin the results into radial bins and output the results in a curve file. Supports scalar fields (density spectra) as well as vector fields(velocity spectra…) Vector fields output the spectra of each component as well as the spectra of the total, the solenoidal, and the diverging components.
Plans can be used to generate initial conditions much faster than summing up sines and cosines… And layouts should allow to take the data in the base grid - and remap it on to a subsection of a finer level (for using forced turbulence results as initial conditions for a clouds velocity field…)…
movie |
Meeting Summary
Hypre
Sorted out several issues with bluegene and hypre - although I am still a little perplexed by hypre's fortran interface… Was able to use it successfully on bluegene - but when I tried the same on alfalfa - I ran into weird mpi errors… I think it has something to with the assumed length of integers passed through the hypre interface??? Our own interface works correctly - so this is probably not worth pursuing… Also discovered that when running with AMR on 100's of processors, that hypre's memory requirement may be quite large unless you use a version of hypre configured with "—with-no-global-partition". See Hypre Memory usage
Restarts
I modified the restart IO to be much faster on bluegene… Restart time was reduced from 3 hours to 3 minutes! See #180. I still need to check this into the repository…
Colliding Flows
These runs are going well on kure. Need to start thinking about analysis… Here's a movie of a clumpy flow. I have a similar setup but with a smooth flow running on kure now.
Gravo Turb
After fixing hypre issues, runs were able to proceed for about .25 free fall times. Then some combination of protections/resteps/etc… caused runs to 'lock up'… Maybe initial setup was too ambitious - or refinement criterion was not sensitive enough to avoid large errors in energy/gravitational potential at coarse fine boundaries… Looking at a similar setup in 2D on alfalfa to try and determine what may have gone wrong… Here's a movie
Collected some code improvement ideas on to DevelopmentProjects
Hypre memory usage
Level | nBoxes | nCells | x_{min} | y_{min} | z_{min} | x_{max} | y_{max} | z_{max} | MinBoxSize | MaxBoxSize | Table_{1} | Table_{2} | Table_{3} | TableSize | TableSize/nBoxes | TableSize/nCells | |
0 | 1024 | 262144 | 1 | 1 | 1 | 64 | 64 | 64 | 256 | 256 | 8 | 8 | 160 | 10,240 | 10 | 0.0390625 | |
1 | 246 | 157536 | 33 | 25 | 33 | 96 | 94 | 96 | 512 | 1344 | 20 | 21 | 8 | 3,360 | 13 | 0.021 | |
2 | 869 | 1132952 | 65 | 53 | 67 | 190 | 188 | 186 | 512 | 3584 | 44 | 46 | 24 | 48,576 | 55 | .042 | |
3 | 4242 | 8477904 | 131 | 107 | 135 | 378 | 376 | 370 | 512 | 20480 | 112 | 116 | 95 | 1,234,240 | 290 | .145 | |
4 | 4201 | 24779672 | 315 | 313 | 321 | 704 | 696 | 704 | 512 | 46080 | 192 | 190 | 189 | 6,894,720 | 1641 | .27 | |
5 | 10653 | 8818504 | 657 | 627 | 659 | 1380 | 1368 | 1362 | 512 | 4608 | 354 | 369 | 337 | 44,020,962 | 4132 | 5 | |
6 | 1902 | 1404112 | 1024 | 1024 | 1024 | 2728 | 2710 | 2642 | 512 | 2744 | 465 | 449 | 444 | 92,700,540 | 48738 | 66 |
The code dies with the following message Out of memory trying to allocate 412094160 bytes
- nBoxes is the number of grids passed into hypre for that level's solve
- nCells it the number of cells inside of those grids
- x_{min} is the lowest index in the x-direction of all the boxes passed into hypre
- x_{max} is the highest index in the x-direction…
- y_{min} is the lowest index in the y-direction etc…
- MinBoxSize is the number of cells in the smallest box
- MaxBoxSize is the number of cells in the largest box
- Table_{1} , Table_{2} , & Table_{3} are the sizes used to calculate the index table
- Table size is the product of those sizes.
So hypre was working as indicated assuming that it was compiled without the '-with-no-global-partition' option. This compiler option is strongly recommended on BlueGene for amr problems using hypre.
Here is a graphical example of hypre's global partition index table… This would results in a 9x6 table with an entry for every black dot. The horizontal lines are along the top and bottom of every grid, and the vertical lines are along the left and right edge.
Weekly Summary
- Patched Back Links Macro
- See ProcessingObjects for an example using
[[BackLinksMenu]]
- While patching figured out how to modify wiki macros and recompile them etc… so I modified BackLinks to automatically be embedded in a collapsible structure instead of appearing next to the
[[PageOutline]]
- After you modify the source code in BackLinksMacro.py you have to:
- compile the source python code to generate a .pyc file.
- restart the web server to refresh the cached versions of the plugin
johannjc@clover plugins > emacs BackLinksMenu.py -nw johannjc@clover plugins > python Python 2.6.6 (r266:84292, Sep 15 2010, 16:22:56) [GCC 4.4.5] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import py_compile >>> py_compile.compile("BackLinksMenu.py") >>> exit() johannjc@clover plugins > sudo apachectl restart
- Also some plugins are packaged as egg files in which case you have to first unzip them - make the changes, compile, and then rezip them. Some egg files are not rezippable, but should have a setup.py script that you can pass to easy_install that should create the egg.
- See ProcessingObjects for an example using
- I couldn't find an easy way to add
[[BackLinksMenu]]
to every wiki page, but did learn about page templates.- Template pages can be created in the Template sub directory. Currently there are three template pages
- PageTemplates/Templates - for any random page
- PageTemplates/TestTemplate - for test problems
- PageTemplates/ProjectPage - for new project pages
- When you create a NewPage you can select any of the existing templates as a starting point.
- All of these templates have the
[[PageOutline]]
and[[BackLinksMenu]]
macros already.
- Template pages can be created in the Template sub directory. Currently there are three template pages
- Reopened ticket #127 after running higher res versions of the field loop advection.
- Added references to IICool to external links page as well as a link to the old astrobear site at the bottom of WikiStart
- Fixed carbuncles appearing in RadiativeInstability test problems #178 and updated the bear2fix.data.img file to show more than just the coarsest grid and modified the problem module to trigger refinement before the first step near the right edge.
- Made various code improvements
- Added CellsPerJeansLength variable which allows users to adjust the jeans-length based refinement. Was previously a parameter = 4. Can be adjusted inside of a users problem module init - not in physics.data
- Duplicated the TestMacro that Matt wrote and had it use a directory on clover instead of the test repository when generating images etc… This enables folks to immediately see the output of the current test on the web before the code gets committed etc… Also added the test results from the last version of the code checked in to the official repository to AstroBearTesting.
- Pondered the Fedderrath particle creation tests
- Met with Eric to discuss GravoTurbulent sim params. Also have a 1024 core reservation that just started this afternoon.
- Running colliding flow sims on kure.
GravoTurbulent Params
Just talked with Eric and we decided to try the following set of runs…
Constant params
- rho=100 particles/cc
- Mach = 5
- cells per jeans length 25 (or higher if computationally feasible)
Variable params
Run | alpha | Gravity | EOS |
A | 1 | SG | II |
B | .2 | SG | II |
C | 1 | SG | ISO |
D | 1 | Static | II |
E | .2 | Static | II |
where
- II is Inoue and Inutsuka
- ISO is isothermal at 41 K
- All runs will start with windowed velocity field of fully developed turbulence.
To achieve a virial parameter of .2, we need to make the density higher or make the box larger… Adjusting the box size is the easiest thing as it will keep the temperature and the velocities the same. Before we were using R = 50 pc and a Mach number of 8.5
We also discussed a steady state gravo-turbulence where mass is taken from small scale collapsed structures and put back in at the large scale with or without additional velocity forcing.
Fully Developed Turbulence
Run | 128^{3} | 256^{3} | 512^{3} | Xmu | M_turb | v_turb | alpha | alpha t_cross/v_turb | t_cross | Temp | SoundSpeed | n | lBox |
ISO | 10 | 2 | 2 | 1.4 | 12.5 | 96.8272888743251 | 2.5 | 0.01333 | 0.516383352062004 | 20 | 4.472135955 | 100 | 50 |
IICool | 10 | 2 | 2 | 1.27 | 12.5 | 66.4996345290655 | 1.1792 | 0.01333 | 0.751883831 | 16.9812533471972 | 5.31997076232524 | 100 | 50 |
Note for the isothermal run, only the mach number of the turbulence is important - everything else is scalable
For the IICool simulation the mach number, the mean density, and the box size all become important considerations…
So now we have turbulence initial conditions that we need to map onto a cloud
The Total Mass of the system in cu is 0.42836E+07 - and since the ambient has a volume of 100^{3-(4/3*pi*20}3^{) and density of 1, this gives a cloud mass of approximately 3.3e6 or a mean density of 51 though it should be closer to 100.. It is likely less because of the softening… If we consider the radius to only be 80% then the mean density would be 99 particles/cc. There is still a net velocity of order 9 cu and a velocity dispersion (mass weighted) of 33. Given that the sound speed is of order 4.5 this is a mach number of 8 }
For the isothermal run - the density can be arbitrarily scaled without effecting the dynamics… For the IICooling run - this is not the case. If should be possible to adjust the window to raise the min density…
Due to an accounting error… the mach number ended up not being the same for the two runs…
Run | 128^{3} | 256^{3} | 512^{3} | Xmu | M_turb | v_turb | alpha | alpha t_cross/v_turb | t_cross | Temp | SoundSpeed | n | lBox |
ISO | 10 | 2 | 2 | 1.4 | 21.5 | 96.8272888743251 | 2.5 | 0.01333 | 0.516383352062004 | 20 | 7.746183109946008 | 100 | 50 |
IICool | 10 | 2 | 2 | 1.27 | 12.5 | 66.4996345290655 | 1.1792 | 0.01333 | 0.751883831 | 16.9812533471972 | 5.31997076232524 | 100 | 50 |
Pondering the Fedderath particle formation requirements
Particle creation requirements
- violates the truelove criterion
- is on the highest level of refinement,
- surrounding flow is converging in each direction,
- has a central gravitational potential minimum,
- is Jeans-unstable,
- is bound, and
- is not within racc of an existing sink particle
- and 2. are the same as Krumholz
- seems reasonable to avoid particle formation in shearing layers
- seems like it would favor particle creation in the center of the grid…
Consider a 1D system with two sink 'planes'. The Greens function is just
and let's consider three points at x=-1,0, & 1. The potential will then beAs you can see only the central point would form a sink particle even though there are 3 equal masses. The paper is not entirely clear about the local gravitational potential and could possibly be referring to a gas potential calculated using only the local mass as opposed to the local gas potential of the entire domain.
Indeed, for check # 5 it calculates the self-energy of the control volume and does calculate a potential using only the local contributions… So they probably are referring to this same local-mass potential with regards to check 4. However, any cell surrounded by a uniform region that violates the true love criterion will have a minimum for the potential at the center of the region because the control volume will be centered on that cell. One could use the periodic approach where the mean density of the region is first subtracted off - in which case the potential would be zero everywhere… and would require the density to be the largest at the center… However, a better check - and one that wouldn't be dependent on the centering of the control volume might just be to make sure that the density is the highest in that cell compared to the surrounding cells in the control volume. This will keep particles from forming in very close proximity. We should update the code with one of these two approaches.
- This requires that the surrounding control volume be as dense as the central cell - otherwise it will violate the true-love criterion without forming a particle… we have…
where
Jeans density criterion is that
which is equivalent toJeans energy criterion is that
or that or in other words, if a region violates the density criterion, we would expect
I've checked in the code that when the density is uniform and just above the jeans criterion that the ratio is
which is consistent with the various assumptions and discretizations… at least when gamma=5/3. When gamma = 1 this changes slightly to closer to .So if there are only a few cells inside of the control volume near the jeans density, it may not satisfy the energy constraints initially… But it seems unlikely that it would artificially fragment without having enough mass present to satisfy this criterion.
After criterion 4, criterion 6 tends to be the most restrictive - as a collapsing region will not form a particle until enough kinetic energy has been converted into thermal energy at the center, and then cooled.
BackLinksMenu macro patched
I patched the BackLinksMenu macro so now pages that link to a given wiki page will be listed in a collapsible section at the top of the page provided the page has [[BackLinksMenu]]
at the top. Also created some new page templates for new pages, test problem pages, and project pages. See PageTemplates/Templates, PageTemplates/TestTemplate and PageTemplates/ProjectPage
Now when you create a new page you can select one of the templates as a starting point. To create new templates just create a page in the PageTemplates directory
Thought about the wiki
So I was thinking it would be useful on the wiki if every page had a sort of directory heading that showed how one would get to the page? For example see ProcessingObjects
Also created a graph showing connectivity between different wiki pages. It's a little too large to see at once…
It looks like there are a couple of useful
Self Gravity sims
I've started running some hires 3D simulations of gravitationally induced turbulence. Here are a couple images of the column density at 17 Myr
as well as a movie
I went with a cloud at a density of 10 particles/cc and a temperature of 160 K and a radius of 20 pc for a total cloud volume of about 8500 M_{Solar}. The cloud material initially has a Jeans Length of 80 pc. I was able to get pressure equilibrium by using the Inuoue and Inutsuka cooling/heating curve which gives an equilibrium temperature (pressure) for a given density plotted below. x-axis is log10 number density and y-axis is log10 pressure in units of (Kelvin/cc) (ignore the horizontal line)
This curve allows for gas at 10 particles/cc and 160 K to be in pressure balance with gas at .23 particles/cc and 6784 K.
It would be possible to raise the density to 100 particles/cc and still be in pressure balance with gas at 1 particle/cc if we wanted to make the cloud more massive.
The cloud is initially kicked with a solenoidal velocity perturbation with a k^{-2} spectra and a total kinetic energy that is only about 1/3 the gravitational binding energy although there is plenty of internal energy to support the cloud
Gravitational Energy | -17.3e46 ergs |
Thermal Energy | 32.6e46 ergs |
Kinetic Energy | 6.14e46 ergs |
Since the Kinetic energy is 1/5 of the thermal energy, the initial perturbations are sub sonic. There is still structure that initially appears, but given the scales, I believe it is driven by thermal instabilities rather than gravitational.
Here are some histograms of density at approximately 8 Myr (red line) and 16 Myr (green line).
and density weighted histograms of the Jeans Length
Eventually the central part of the cloud collapses and forms 4 sink particles at around 17 Myr, while the rest of the cloud is still fairly whispy.
Given that the resolution is 100 pc / 2048 the code should form particles when the jeans length is .2 pc which should occur at densities of 1e5 with a corresponding jeans mass of about 10 solar masses. Given that the cloud contains 8000 M_{Solar} eventually there should form 100s of 10 solar mass cores - provided there are mechanisms for increasing the density to 1e5. Global contraction may be necessary to provide the density enhancement - as the cooling instabilities can only compress the material to 120 particles/cc or so…
In the colliding flows, we are embedding dense small clumps in the colliding streams to provide perturbations. We want these clumps to be in pressure equilibrium with the stream, and for the stream and the clumps to be in a thermally stable region (n < 1 or n > 10). It is also good if the clumps are large enough to be resolved, but be light enough to not plow through the opposing flow. The two stable densities allowed with the minimum density contrast occurs at arround 10 and .23 particles/cc. We decided to have a mean density of 1 particle/cc and a flow radius of 20 pc which gives a Mass flux of 523 m_sun/Myr for about 15 Myr or a total mass about 8000 M_{Solar} which would be comparable to the mass in the spherical cloud.
It would also be interesting comparing this to an isothermal run at 160 K or perhaps much lower… I was also thinking it would be good to always refine the Jeans length by 64 cells or so - until some maximum resolution. Also I've noticed that these runs are rotationally symmetric which is due to a difference in bluegene's compiler handling certain random number generating functions differently. That issue has been fixed.
Storage space
Just a heads up that alfalfa ran out of space - and we only have 3.7 TB of free space left on the three current raid arrays. Here is a summary of current data usage:
user | alfalfa | grass | clover | total |
blin | 1.84 | 0 | 0 | 1.84 |
johannjc | 1.49 | 0.71 | 2.85 | 5.05 |
martinhe | 0.99 | 0.66 | 0 | 1.65 |
shuleli | 0.14 | 0.57 | 0 | 0.71 |
erica | 0.01 | 0.26 | 0.29 | 0.56 |
trac_backup | 0 | 0 | 0 | 0 |
ehansen | 0 | 0 | 0 | 0 |
chaig | 0 | 0.58 | 0 | 0.58 |
yirak | 0 | 0.31 | 0 | 0.31 |
stanny | 0 | 0.24 | 0 | 0.24 |
eschroe3 | 0 | 0.22 | 1.93 | 2.15 |
bshroyer | 0 | 0 | 1.81 | 1.81 |
used | 4.7 | 4.6 | 7.7 | 14.9 |
free | 0.5 | 0.4 | 2.8 | 3.7 |
total | 5.4 | 5.3 | 11 | 21.7 |
Meeting update
- Implemented Bondi Accretion… Updated SinkParticles with description of Bondi accretion routine… Someone should update documentation on SinkParticle
- Debugging Testing script… Making sure it works and catches failures etc…
Procedures for updating scrambler repo
So i checked in what should be a newer version of the code to clover:/data/scrambler_devel that should be tested and (if successful) pushed onto clover:/data/scrambler.
Before testing on clover, someone should run hg update to update the source files to match the newly pushed head. I also was looking through all of the scripts and found…
- buildtests.s - Brandon's old suite of three tests intended to be run on bluehive.
- testbuilds.s - My script that attempts to compile every problem module
- test_scripts/test_Xproc.pbs - pbs submission scripts for Brandon's old suite of tests.
- modules/BuildModuleControl.pl - empty… I believe at one point Brandon or Matt intended to use this as a type of configure script to turn on/off the compilations of various modules
- tests/go.s - script that does post processing on an individual test problem passed as an argument
- postprocess.s - script that cycles through a list of tests and runs go.s (currently only setup to test FieldLoopAdvection)
But i currently don't see Matt's script or any script that would cycle through the problems and actually run them - or update the output on the wiki in astrobear in the scrambler repository
I did find an svn repository cloverdata/repositories/tests/ that has the images and the logs from the last test run that show up if you use matt's wiki plugin
https://clover.pas.rochester.edu/trac/astrobear/wiki/u/noyesma#no1
And there is a copy of the repository at cloverdata/tests in which there is a process_chombos script which I imagine is supposed to update the test logs etc… and then check in the new images to the tests repository so they can be displayed on the wiki… but there should be another script somewhere that actually runs each test problem before calling this processchombos script…
Update
So the buildproblem script can be found in /home/noyesma/bin and Matt's screencast is available at
http://dl.dropbox.com/u/2018623/Test%20Suite%20Tutorial.mov
Meeting update
- Cleaned up sweep scheme - and self gravity (soon to be checked into main repo)
- Implementing Bondi accretion
- Improved protection options. See ProtectionOptions
- Cleaned up cruft in PhysicsData namelist…
!================================================================================ ! AstroBEAR physics.data file. ! Input data describing built-in physics-related parameters of the code. ! Created for problem: Template ! Default values for variables in [brackets] !================================================================================ &PhysicsData !================================================================================ ! Equation of State stuff !================================================================================ gamma = 1.66666666666667 ! Adiabatic index for ideal gas [1.66666666667] Xmu = 1d0 ! mean molecular weight [1d0] iEOS = 0 ! [0]-ideal gas, 4-Isothermal (Equation of state) !================================================================================ ! Field based refinement control !================================================================================ InterpOpts = 0,0,0,0,0,0,0,0,0,0,0 ! [0]-constant, 1-minmod, 2-superbee, 3-vanLeer, 4-mc, 5-Parabolic (not conservative), 6-Linear (not limited). It is recommended to use constant interpolation for any fields in q that correspond to aux fields (iBx, iBy, iBz, etc...) refineVariableFactor = 1d0,1d0,1d0,1d0,1d0,1d0,1d0,1d0,1d0,1d0,1d0 ! weight factors for each field used in triggering refinement. Combined with qtolerance for that level. [1d0] !================================================================================ ! Global source term switches !================================================================================ lSelfGravity = .false. ! Turns on self-gravity if true [.false.] UniformGravity = 0d0 ! Gravitational acceleration in the -y direction [0d0] iCylindrical = 0 ! [0]-No cylindrical geometry, 1-Cylindrical with no angular momentum, 2-Cylindrical with angular momentum !================================================================================ ! Density Protection Options !================================================================================ lRestartOnDensityProtections = .false. ! Do density protections trigger restarts? [.false.] iDensityProtect = 2 ! 0-Set to MinDensity, 1-Set to minimum nearby density, [2]-Average nearby densities iMomentumProtect = 2 ! 0-Conserve momentum, 1-Set to zero, [2]-Average nearby velocities MinDensity = 1d-10 ! Minimum computational density before protection is triggered [1d-10] !================================================================================ ! Pressure Protection Options !================================================================================ lRestartOnPressureProtections = .false. ! Do pressure protections trigger restarts? [.false.] iPressureProtect = 4 ! 0-Set to MinTemp, 1-Set to minimim nearby pressure, [2]-Set to average nearby pressure, 3-Set to minimum nearby temperature, 4-Set to average nearby temperature MinTemp = 1d-10 ! [1d-10] minimum allowed temperature for the system in Kelvin before protection is triggered !================================================================================ ! Description of various scaling parameters ! Define one of each of the following: [nScale/rScale], [TempScale/pScale], and set the other to 0d0. !================================================================================ nScale = 1d0, ! number density scale parameter [particles/cc] rScale = 0d0, ! density scale [g/cc] TempScale = 1d0, ! temperature scale parameter [Kelvin] pScale = 0d0, ! pressure scale [dynes/cm^2] lScale = 1.4959800E+015, ! length scale parameter [cm] (AU=1.49598e13, pc=3.08568025e18, R_sun=6.955e10 !================================================================================ ! MHD related section !================================================================================ lMHD = .false. ! Magnetic Fields present? [.false.] lCheckDivergence = .false. ! Turn on divergence checking [.false.] /
Running visit remotely
So… For folks interested in visualizing data remotely, I recommend installing visit on their local machine and using a visit data server on the host… I also recommend using the same version of visit on your local machine as on the remote host machine.
Here is a sample configuration for alfalfa which you can add under Options→Host Profiles
Hydrostatic Isothermal disk
The goal is to set up a hydrostatic isothermal rotating disk. Need to solve
as well as
for
and where is the gravitational acceleration due to a point mass at the origin
First we can simplify the equations
as well as
Since omega is a free parameter, we can choose a density field that satisfies the first equation and then integrate the second equation to find
as long asLet's try
where
describes the density along the midplane. We can set the height of the disk to be where the density drops to the ambient value - or we can set the background density distribution and replace it with the disk to the height where the pressures balance.Also, if we don't use softening, then [[latex($f_{z}=-\frac{GMz}{r^{3} }
and
We can then solve the second equation for
If we now consider the density along the midplane
Since we expect
and since we have
or
or
Of course this blows up at the origin, so we can either avoid the origin or use some form of softening.
Centrifugal acceleration can only support things from going inward. To keep things from going outward we have to use pressure or gravitational force - Since we want the density to decrease as we go outward, we have to use gravitational force - but this limits how quickly the density can drop. Lets choose
where
. If then we expect the equivalent of a BE sphere (but with a point particle instead of self-gravity).. Solving for along the midplane…
Combining all of this we get a solution for the density field
and rotational velocity goes like
Now for softening
First lets consider plummer softening. We can repeat the above procedure but instead of
we have where is a softening radiusLet's again try
Now the integrand is
which is very similar to before except that andso
And the integral along the midplane is also modified…
Combining these gives the solution for density
and angular velocity
Spline Softening
If the softening function is compact (only modifies gravity inside of a certain radius) then
So we need the integral of the softened function…
Interpolation Verification
After seeing strange behavior - especially wrt pressure I decided to check that the interpolation methods used in prolongation were accurate. It is a little difficult to visualize this but it is possible.
- Start a run in fixed grid.
- Select the frame you want to see interpolated.
- Set the restart_frame to be that frame in global.data
- Set the final frame to be 1 frame later.
- Now the trick is to get the time between those two frames to be very small, so that there is no temporal development. To do this…
- Get the time stamp of the restart frame (restart_time = restart_frame/final_frame*final_time) assuming start_time was originally 0
- Set the start_time to be restart_time-(restart_frame*epsilon)
- Set the final_time to be restart_time+epsilon
- Then increase maxlevel to whatever desired level and set qtolerances to trigger refinement.
- Then restart the simulation
- It should take one root level step and then output a single frame. When viewed it should show the interpolated values.
Here are images from an example where 4 levels of refinement were added to a BE sphere profile.
Constant (0) | MinMod (1) | SuperBee (2) |
VanLeer (3) | Monotone Centered (4) | Parabolic* (5) |
Linear (6) |
The "plateaus" in all but the linear and the parabolic are due to the slope being limited so as not to introduce any new extrema. The Parabolic does not conserve the total quantity and is only used for the gas potential (and time deriv). The 'linear' was just added but could likely have issues near shocks since it is not limited. But all seem to be working 'correctly'. There is a shortcoming in that each limiter is considered in each direction independently. A more multidimensional reconstruction that used a multidimensional limiter would ease the plateaus from persisting along the coordinate directions.
Issues with Self gravity
For a fixed grid run we need phi before and after every advance. So for 4 time steps, we would need to solve for phi 5 times. With no additional elliptic solves we can then use phi before and after to achieve second order time accuracy.
For an AMR patch we need to take 2 time steps so we need 3 solutions for phi. We could use the prolongated parents version of phi, or the previous overlapping grids latest version of phi. In principal they should agree at the coarse boundaries, since the previous grids solution to phi was subject to the same boundary conditions (due to phidot at coarse boundaries exactly agreeing with parent grids). Furthermore the previous grids solution to phi should be more accurate since it accounts for the distribution of gas on the current level. Otherwise level n would end up using values for phigas that were successively prolongated from level 0. So it is important to use the overlapped version of phigas and not the prolongated.
However, phidot should not be copied from old grids since that phidot is from a previous time step. This can cause problems if phidot changes quickly. There will also be a mismatch between the prolongated and overlapped phidot which can cause differences between adjacent cells leading to strong gradients. So it is important to use the prolongated value of phidot and not the overlapped. Additionally since phidot is recalculated before successive prolongations, the same issue with repeated prolongation does not occur.
This is implemented in the code by setting
ProlongateFields = (/…,iPhiGas,iPhiDot,…/)
and by setting EGCopyFields = (/iPhiGas/)
And then following a poisson update ghosting iPhiDot.
To make it 2nd order we need to ghost iPhiGas as well after the poisson update, however this means that we don't need to copy iPhiGas before the 2nd step.
Simpler Self gravity implementation seems to work
So, I gave up on trying to integrate the potential due to particles with that due to the gas to get a better gas potential at coarse fine boundaries. Now, only
and are stored in q. Gas that is accreted during a finest level step will not be reflected in the boundary values to the gas potential at the coarse boundaries, until after the coarser grid updates. Fortunately, since the accretion is incremental this should not pose a significant problem. Additionally if the accretion is steady, then the time derivative of the gas potential should include this effect. I have performed a suite of 2D simulations with various levels of refinement and all seem to handle the collapse nicely.movie |
And the 1 3D TrueLove problem I ran also seems to handle the collapse well. In the absence of accretion or particles, the scheme should be second order in time - although I have yet to verify this.
Also, it was necessary to reduce the softening length to 1 cell, instead of 4 cells, since the particle would form at a cell corner, but the gas would collect in a cell center. Due to the softening the gas cell would exert a greater pull than the particle which would cause the gas to wonder away from the particle.
New ways to be notified of changes to the wiki
I installed a couple of plugins so now…
- Blog authors will be e-mailed when someone comments on their blog posts (although this e-mail can only be sent to their username@pas.rochester.edu)
- Wiki pages now have a 'Watch This' link in the upper right corner. If you select to watch a page, you will receive e-mails whenever the page is modified.
- Under the Preferences menu (right next to logout), their is an announcement tab where you can adjust various settings related to ticket notifications. (I would suggest folks check all of the boxes under "ticket component subscriptions" if they want to stay on top of all tickets)
Resources
Some sample numbers for 3D runs
- 32^{3}+4 = 512^{3} effective with a high filling fraction
- 2.2 GB / frame x 100 frames = .2 TB
- 64 procs @ 1.2 GB of mem / processor = 76 GB of active memory
- 12800 SU's
- 32^{3}+5 = 1024^{3} effective with a smaller filling fraction
- 1.5 GB / frame x 100 frames = .15 TB
- 64 procs @ 1 GB of mem / processor = 64 GB of active memory
- 19200 SU's
- 32^{3}+5 = 1024^{3} effective with a high filling fraction (projected)
- 12 GB / frame x 100 frames = 1.2 TB
- 512 procs @ 1 GB of mem / processor = 512 GB of active memory
- 102400 SU's
Projected needs for research
- 6-10 runs for colliding flows
- 6-10 runs for gravo collapse
- 12-20 runs at 19200 SU's = 230400 - 384000 SU's
- 12-20 runs at .4 TB = 4.8 - 8 TB of data
This is effectively the entire 192 infiniband nodes for 50-80 days or all of bluegene for 1 to 2 weeks (assuming memory issues are not a problem). If we had exclusive access to a quarter of bluegene, we would be able to finish these runs in 1 to 2 months.
Projected needs for groups research
This is of course, just for the research I would like to do in the next few months. When you consider the rest of the group we are looking at maybe 3-4 times that?
As far as storage space is concerned, we currently have 21 TB of storage space of which 5 TB is available. If each grad student/post doc has 5-8 TB of active data that needs to be stored until a paper is published, then we would need at least 30 TB. (That's assuming folks get rid of their old data as soon as the results are published). At some point (when we have the new bluegene machine) we will be generating too much raw data to effectively store and we'll have to figure out how to generate useful output (ie spectra, column density, only coarse grids, etc…)
Working on implementing conservative momentum conservation in 1D with self gravity
Need to:
Remove self gravity source term from source moduleAdd self gravity source terms to reconstructed left and right statesmodify final fluxes to include self gravity source terms.
Here's a plot showing momentum conservation with the new approach
However there were striations
These were due to the approximation it uses to calculate the density… Turns out when the errors in phi are large compared to the density this can cause these striations.
Here's a plot showing the relative error in the derived rho used in the source calculation and the actual rho.
Reducing hypres error tolerance from 1e-4 to 1e-8 improved the situation and lessend the striations
Finally here's a comparison of the old and new methods
Modify gravitational energy source terms.These modified fluxes will need to be stored in the fixup fluxes so they can be properly differenced between levels to ensure strict momentum conservation.this will require a second call to store fixup fluxes. Perhaps a generic routine in data declarations would be better for storing fixup fluxes, and emfs, etc…
extend same modifications to 2D and 3D by adding additional source terms in transverse flux update.
What to do about phi
- To be second order we need the gas potential before and after a hydro step.
- Normally the gas potential after a hydro step would be the same as that before the next hydro step (requiring 1 poisson solve per hydro step) however with accretion, the gas in each cell (and the corresponding gas potential) can change.
- Accretion, however, should not change the total potential (except for differences in softening) - so recalculating the particle potential should allow for modifying the gas potential without another poisson solve.
So a fixed grid algorithm would look like the following:
- Do a hydro step using original gas potential and predicted time centered particle potential.
- calculate new gas potential (poisson solve) and ghost
- momentum and energy flux correction
- Advance particles using back force from gas
- Calculate new sink potential
- Store new total potential
- Perform accretion
- Update new sink potential using new particle masses and difference gas potential keeping total potential fixed
- Repeat
For AMR we need a way to update the gas potential on level l at coarse grid boundaries independent of a level l poisson solve. This can be done using phi and phidot. Then whenever we set the gas potential for the poisson solve we use phi, phidot, and phisinks. So the algorithm looks like the following:
Root Level
- Do a hydro step using original gas potential and predicted time centered particle potential.
- Calculate new sink potential using predicted particle positions
- calculate new gas potential (poisson solve) and ghost
- momentum and energy flux correction using phi_Gas and phi_gas_old (stored in phi)
- Update total potential and time deriv using sink potential and gas potential. phi_new=phi_gas+phi_sinks phi_dot=phi_new-phi+old_phi_sinks
- Prolongate old fields and new total potential time deriv
- After finer level steps, particle positions and masses will be different. So update phisinks as well as phigas keeping phi constant.
Intermediate Level
- Do a hydro step using original gas potential and predicted time centered particle potential.
- Calculate new sink potential using predicted particle positions
- Update gas potential at coarse fine boundaries using phi, phidot, and predicted phi sinks
- calculate new gas potential (poisson solve) and ghost
- momentum and energy flux correction
- Update total potential and time deriv using sink potential and gas potential.
- Prolongate old fields and new total potential time deriv
- After finer level steps, particle positions and masses will be different. So update phisinks as well as phigas keeping phi constant.
- Repeat
Finest Level
- Check for new particles (do this after ghosting)
- Perform accretion
- After accretion, particle positions and masses will be different. So update phisinks as well as phigas keeping phi constant.
- Do a hydro step using gas potential and predicted time centered particle forces.
- Advance particles using back force from gas
- Calculate new sink potential using advanced particle positions
Calculate new sink potential using predicted particle positions- Update gas potential at coarse fine boundaries using phi, phidot, and phi sinks
- calculate new gas potential (poisson solve) and ghost
- momentum and energy flux correction
- Update total potential and time deriv using sink potential and gas potential.
After finer level steps, particle positions and masses will be different. So update phisinks as well as phigas keeping phi constant.- Repeat
Of course the calculation of the predicted particle positions on each level (other than the maxlevel) can be as sophisticated as necessary. One could use the original particle positions and velocities alone, or could advance the level's own version of the particles using the same advance algorithm as well as gas forces etc… Note this is necessary if one wants to thread these advances, since the particle masses may change due to accretion in the middle of a coarse level update. But this threading would still require use of the old time derivative to update the ghost zones instead of the forward time derivative.
Added 1D functionality to code
The code can now run 1D AMR problems and produce output to chombo. Just set nDim = 1. The chombo files will look funny with AMR turned on, but that is just because chombo has to believe they are 2D AMR data sets. Because the data is in fact 1D, it thinks that some data is missing, and leaves these areas white.
You an do line-outs along the y=0 boundary to generate curves.
Here I've plotted the data on levels 0 and 1 (redline) with data on levels 0-4 (blueline).
(Also have not tested self gravity or other source terms, or several of the 'objects', but most will need very minor modification to work).
It is not checked in, but interested folks can pull from alfalfadata/johannjc/scrambler_1Dworks
Meeting Update
- Looked at how Athena handles CTU w/ source terms
https://clover.pas.rochester.edu/trac/astrobear/wiki/CodeExplanation/SourceTerms
- Looked at embedding clumps in flows where both clump and background flow are in pressure equilibrium. (Also updated Clump objects to now support velocities and can appear from off the grid.)
https://clover.pas.rochester.edu/trac/astrobear/attachment/wiki/CollidingFlows/clumpflow10b.gif
- Talked with Rich about mounting cloverdata, grassdata, and alfalfadata, across machines to hopefully avoid folks running sims and having data be dumped to the department disk. Also discovered we have another machine besides alfalfa (botwin) that was intended to be used as a backup? Dave is still waiting to hear back from LSI about raid controller software
Running 32/64 processor jobs on the standard queue...
It seems that running 32/64 processor jobs on the standard queue is no longer practically feasible… I believe this is due to the backfill policy favoring the plethora of single processor jobs. The only 32+ job I've gotten to run lately on the standard queue was after the CRC downtime when backfill would not have been a factor. Steve recommended using the ib queue for larger jobs although this queue restricts us to using 14 nodes plus the 8 afrank nodes instead of the 94 nodes available to the standard queue… I have however, gotten a job to run on the ib queue there with a reasonable wait time although the wall time is limited to 2 days instead of 5 so this means more frequent restarts and more babysitting…
Also talked with Eric and we decided to do 3D runs of non-linear initial density perturbations and virialized velocity perturbations
A little confused about the conversion of gravitational energy to kinetic energy
The last run shows the collapse of a slightly perturbed sphere. I would expect that the total energy is conserved or that KE+U+iE = C… The run is however, isothermal - so any heating due to compression is lost… but does this explain the essentially flat curve that is 2KE+V? https://clover.pas.rochester.edu/trac/astrobear/wiki/GravoTurbulence#VII
So I reran the simulation with a gamma of 1.001 to track this energy. Sure enough the rate of increase in the internal energy is very similar to the kinetic energy during the initial stages of the simulation.
Submitted Scrambler Paper!
Hooray
Meeting Update
- Still working with Dave on getting raid scripts up and running. Clover is optimal - but not sure about grass or alfalfa
- Submitted paper to astro-ph
- Worked on setting up clumplets
- First 3D results from colliding flows
- Perhaps with postprocessing scripts and various image production capabilities it would be possible to have automated output from simulations go to a web directory with a html template that displays pertinent information from the run
- computational scales
- resolution
- plots of density
- histograms of various thingies
- a curve showing the time development of mass, energy, momentum, etc…
We could even have pull down menus for showing various plots quantities - and Kris's handy image scroller ? I think this would save a lot of time in the long run and would provide a convenient way to view data. Of course it would be no substitute for a nice annotated visit plot and having pages on the wiki that discuss the various runs would be good - but folks wouldn't have to generate a lot of content for each run… but it would still all be accessible.
I know it doesn't look like much...
Here's a snapshot of the first frame. There is a parent clump and three child clumps all with radii equal to the jeans length for the density… Additionally the parent clump has a radius equal to the jeans length for the average density of the clump (including it's children)… Soon will extend this to several generations of clumps.
Web-monitoring of simulations
I found a very easy way to generate jpeg images of various processing objects (ie projections of column density) at runtime (or post processing). There is an image format called ppm that is essentially ascii text that can be converted to a jpeg using system calls. Also added code that checks for a postprocess script in the run directory and if present executes it passing in the frame number as a 5 character string. This can be useful for scping chombo files to a longer term storage, or sending jpeg files to a web-accessible folder for monitoring simulations in realtime…
#!/bin/bash mv out/Mass_along_3_$1.jpeg /home/johannjc/public_html/CollidingFlow.jpeg
Also created cron jobs to back up the wiki daily, weekly, and monthly. In the process discovered that the raid_array scripts that e-mailed folks if the arrays were in trouble was lost in the reconfiguring of grass and clover. (And alfalfa is in need of one). Currently working with Dave to get these back up.
Documentation for Processing Objects
I've updated the wiki with documentation on processing objects linked from ScramblerIntro
I've also updated the GravoTurbulence and the CollidingFlows pages with more results.
Added picture of the day to the main wiki page.
Go to Image of the Day for instructions on adding your own png images.
Working on adding processing routines to AstroBEAR
There are several plots I need to generate for the colliding flows problem (and the gravo-turbulence problem) that cannot be generated in visit - 2D distribution functions, column densities, total particle quantities. There is also no easy way in Visit to plot the total particle mass as a function of time… So I've added some processing routines to astroBEAR that can be done during the simulation or after in a post-processing mode. So far I am not generating any images directly - so no new libraries (pgplot, matplotlib, etc… although this might be a nice thing to add eventually), but instead generating the raw data for the plot so the actual plot can still be customized in visit.
A comment about mercurial
Just wanted to encourage people to never copy files from one astrobear repository to another. This essentially corrupts the repository since mercurial will assume the changes to the file were intended. It is difficult to debug corrupted repositories since it is no longer clear what modifications to the source code have been made, and a corrupted repository should never be committed since that will potentially (and quietly) undo bug fixes. If there is a repository with modifications that you would like to replicate in another repository, you should commit those changes to the first repository and then pull that changeset to your own repository. This way mercurial knows where the changes came from and can correctly merge those changes with any that you might make down the road.
Copying data files is not as 'bad' as long as the copied data files are never committed. (Just use hg revert *.data
before you commit). If you do want to update the data files in the repository, then you should not copy them from another directory, but instead commit the changes in the directory where you first modified them - or revert them and replicate the changes manually.
Also, you can commit individual files so if there is a single bug fix in a single file, just run hg commit somefile.f90
and then pull that change to the other repository and update or merge if necessary.
Zooming movies
So the uniform density collapse problem has a self-similar solution where
and
Since
If we plot
as a function of we get the following.Since the curve looks like a circle we try a polynomial expansion
The best fit (plotted above) has coefficients
If we then use this to scale the size of the window, the apparent radius of the cloud remains fixed as it collapses. At a finite time the radius would go to zero, so once the size of the window becomes of order 32 finest level cells, the box stops shrinking and the cloud finally appears to collapse. The dynamical time also decreases as the cloud shrinks, so in the following movie, the frame rate increases 10 fold and the collapse appears to proceed more slowly. This happens before the box stops shrinking (before the cloud appears to collapse).
I also restarted the run but with sink particle formation suppressed. The results are similar…
New refinement criteria for self gravity based on gradient of potential
If the gradients in Phi are much larger then the gradients in pressure, then discretization errors in the potential can lead to large errors in pressure, which can lead to negative temperatures, and perhaps explosions?
The new criteria normalizes the potential by the sound speed squared and then looks at differences between adjacent cells. This is essentially
Also found a bug that caused particles to sometimes not accrete when they should. Explains why I was seeing zero-mass particles being occasionally created. I will check these in as soon as the code passes the suite of tests… In the meantime folks can pull from grassdata/johannjc/scrambler_1110007 if they are running simulations where accretion is expected to happen.
Time Derivatives of Potentials...
So I think I figured out why the particles generated little bubbles when they are initially formed and it has to do with a missing factor of 2 in the force . At first I thought it was the way potentials were stored/updated. and what follows explains how to better do this, but the potential of a line charge is
so the force is proportional to not . And here is the explanation of how potentials are stored/updated…When there are no sink particles, the gas potential
, and it's time derivative are used to interpolate the solution for at coarse fine boundaries. When a sink particle is formed the gas potential should respond to the missing material. The elliptic solver will make this adjustment provided that the boundary value for is correct. However this is problematic at the coarse fine boundaries since the value of at the coarse fine boundary was prolongated from the parent grid pre-accretion. If instead of storing and using at the coarse-fine boundaries, we use then this time derivative of the combined potential should be unchanged under the conversion of gas material to sink particles - however (and this is the tricky part), this is only true if both potentials use the same constant of integration… The elliptic solver normalizes the potential so that the sum is zero. However, the potential of a sink cylinder goes like ln(r) and cannot be conveniently normalized. One could calculate the integral of the potential over the entire box and subtract off the mean value analytically (or numerically), however if a multipole expansion is used instead then this becomes unecessary. A simpler normalization could be used for the potential. where a is chosen so that the integral over the box of is zero. For a circle of radius , the constant . If we compare the potential of a particle at the center of the grid with an overdense cell at the center of the grid we see that there is fairly good agreement if we useHere are two lineouts through the center of the grid with a particle and an overdense cell
And here is the same but with a softened particle potential
And here is the result if the particle is placed at [.125, .125]
And here is the full 2D image of the potentials
Not the periodic nature of the gas potential causes the turnover. Doing this with image particles is not generally feasible.
And here is the result of a simulation with the correct particle force. Click here for a movie.
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…
1^{st} order reconstruction | 2^{nd} order | 3^{rd} 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 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).
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 likeThe 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)
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 .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.
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
- Also updated DevelopmentProjects page to link to those tickets
- 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 ../
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.
Cleaning up portable thread implementation
I have a reservation for a 8 hour interactive session on 32 of the afrank nodes tomorrow during which time I will be doing some final testing of the new dynamic load balancing scheme and making some minor tweaks. Today I am cleaning up the code changes made for the Portable threads and checking the implementation so I will be ready for tomorrow and also so that they can be checked into the main repo. Also working on copying and pasting the paper into the proceedings format for the AstroNUM 2011 conference.
Encouraging results
While the reservation was for 4 nodes with 8 procs per node, the afrank nodes had hyperthreading enabled - which caused mpirun to put 16 processes on a single node or two processors per proc. The idea behind hyper threading is that you overtask nodes so that they can switch back and forth between various jobs. By giving each processor two instances of the code, when one process is stuck waiting for a message, the other process can still work. The processor time per cell update would appear to double, although the actual core time per cell update should remain constant - or drop slightly if the hyperthreading gives an improvement. Going from 8 to 16, the cputime per cell update more then doubles indicating that hyper-threading actually hurts performance. it also makes it more difficult to schedule tasks at the application level, since the processor is switching back and forth between two instances of the code. The take away is that hyperthreading should not be used with AstroBEAR 2.0. Notice that the dynamic load balancing with the portable threads still gives a better performance whether or not hyperthreading is enabled.
Threading Mode | No Threading | Scheduling | Portable Threads | |||
---|---|---|---|---|---|---|
Level Balanced | T | T | F | T | F | |
8 | 3.0360618888 | 3.5354251260 | 3.9989763675 | 2.4490875460 | 2.3749391790 | |
16 (HyperThreading) | 7.1959148653 | 7.1250127389 | 5.3294605455 | |||
32 (HyperThreading?) | 12.1093259489 | 14.4666908237 | 9.2395890449 |
(Not sure if the 32 core run used both nodes or not…)
I was also somewhat surprised to see that global load balancing finally improved performance. Global load balancing is the opposite of having every level balanced. Previously performance was better when each level was balanced via splitting. While this created more smaller grids and therefore increased the overall workload - the improved distribution outweighed the cost - although perhaps because the load balancing scheme was not correctly compensating for imbalances.
Here is an image with the distribution showing the additional fragmentation on the right when Level Balanced is true. The top panels show every patch while the bottom panels show the mpi_id of the patches.
Finished looking at results from 90^{2 cells per proc per level test… Tomorrow I have a reservation and will generate similar data for 16}2 32^{2 64}2 128^{2 and 256}2 if there is time. As the number of cells per proc per level increases, all of the lines will flatten out - but i should have all the necessary data for the paper.
Load balancing algorithm
I made some minor adjustments to the load balancing algorithm to match the description in the paper and the initial results are encouraging… Now the threading / dynamic load balancing is 10% faster on 8 cores… Need to do some scaling runs on 32 / 64 to confirm that this trend persists. Also trying to figure out the right functional form to fit the scaling data.
Weak Scaling
Weak scaling results for various levels of AMR for 128^{2} cells per proc per level
And the weak scaling results for 3 levels of AMR with 64^{2}, 128^{2}, and 512^{2} cells per proc per level
If we let
represent the total number of cells to update, and represent the computational cost per cell and be the additional computational work that is incurred with processors, and let be the wall time thenFurthermore if
then for we can taylor expand which combined with the scaling results givesAlso
However, if we assume that
thenThe problem is we don't have enough data points to determine which functional fit is accurate so it is difficult to project.
So I went back and looked at different possible fit functions and plotted the results…
The best fit is the first one
and if it holds, implies performance drops to 50% at or 1.74 KProcs and to 25% at 5.3 GProcs! Granted this is weak scaling and your resolution at 5.3 GProcs would be 36 MCells^{2} and simulating one crossing time would take 12 yrs instead of 3 yrs… but that's why strong scaling is more important… or being able to weakly scale when the workload per processor is only a few cells.Why would the scaling be proportional to
? If we consider the distribution of grids at each time step and assume that the workload distribution follows a normal distribution with some standard deviation , then the time to complete each step is not given by the mean of the workload but by the max of the workload . The average maximum value is given by the last order statistic where is the number of processors. I couldn't find an expression for the last order statistic for a normal distribution, but monte-carlo experiments give a fitting function . Given the best fit this implies thatExtended Disk simulations to 3D
Did analysis of disks looking at how the gravitational softening length effects the presence of a jet and the loss of angular momentum. Also confirmed that 2.5D produces similar results to 3D axi-symmetric runs. See the results here.
Meeting Update September 6th
Scaling tests on infiniband completed successfully. Just need to parse the data to produce figures for the paper. Also ran a 2D disk simulation to test angmom conservation and posted to the wind disk interaction project page
Working on Scaling runs on bluehive's infiniband nodes and example of load balancing across threads
Currently running two sets of scaling tests. All simulations are 2D
Weak scaling
filling fraction | .25 |
MaxLevel | 0,1,2,3,4 |
Threaded/Every Level Balanced | N/Y, Y/N, Y/Y |
nProcs | 1 2 4 8 16 32 48 64 80 96 104 112 |
Base Res^{1} (per processor) | 64 128 512^{2} |
- ^{1} Domains are kept square with number of cells per processor approximately constant
- ^{2} Only with a MaxLevel of 3
Strong Scaling
filling fraction | .1 .2 .3 .5 .7 1 |
MaxLevel | 4 |
Threaded/Every Level Balanced | N/Y, Y/N, Y/Y |
nProcs | 1 2 4 8 16 32 48 64 80 96 104 112 |
Base Res | 32 64 |
Meeting Update for 8/29/2011
I've implemented GNU portable threads and have done more scaling tests on bluehive (over ethernet). I've been combing through the data making figures for the paper's results section.
Scaling Tests / MultiClump Runs
I put some useful scripts on the wiki for running multiple jobs with various parameters on different numbers of processors.
Both Hydro runs are almost completed on bluehive. The magnetized runs were dying (see comment:58:ticket:121) but I've reduced the default message size and resubmitted the jobs. I could also move these over to bluehive once the scaling tests are completed…
Meeting Update 8/15/2011
Parsing data from scaling tests. Collected 20+ statistics from 2800 runs… Preliminary plots are attached to ScramblerPaper
Truelove Uniform Collapse
Confirmed truelove collapse problem still works correctly
Uniform Collapse page and movie
The particle wanders and there are a few things I want to try to determine the cause…
- Switch from using from previous coarse advance, to the current coarse advance. This prevents level threading but should still give good performance.
- And do a resolution study to see how it affects the outcome.
- And check the code to make sure the accretion forces and the gas forces on the particle have the proper signs.
(Just finished a level 4 run with no significant drifting… Might be a restart issue)
Paper figures
I made a few additional paper figures for the ScramblerPaper. Now just need to integrate them into the descriptions in the text.
.
Documented the ambient object and process.data file
Documentation on the ambient object is now available https://clover.pas.rochester.edu/trac/astrobear/wiki/AmbientObjects
I also added documentation about the optional process.data file https://clover.pas.rochester.edu/trac/astrobear/wiki/ProcessDataExplained
MultiClump Runs - finished - sort of...
Updated the MHD Clumps Page with information about the By_sparse run… There seem to be some issues with the PPM and magnetized cooling rarefactions that produce ripples…
Back on bluehive and bluegene
Apparently my account just had a lock on it - similar to other users following the attack on the crc clusters… but because I was out of town and didn't discover or report the issue until a week later - the dots weren't connected… (It didn't help that my ticket was re-assigned to someone on vacation)
Anyways, I should be able to get the multiclump runs back up and running on bluegene… (That and work on the paper)
Visited UNC to work with Christina on setting up Colliding Flows
I had a good time in Chapel Hill working with Christina and Fabian on the Colliding Flows module. We were able to get 2D AMR Hydro and MHD working before the queuing system stopped working. We did manage to perform a 3D AMR Hydro 64+4 simulation on 128 processors for an effective resolution of 512. The code seemed to perform quite well tracking the NTSI interface and Christina should be posting some nice screenshots soon
Fun with Thermal Instabilities
I created a page dealing with Thermally induced collapse of clouds…
More Fun with Colliding Flows
I added several new plots showing cooling lengths, jeans lengths, destruction times, etc… to the CollidingFlows page.
CollidingFlows project page
I've created a project page for CollidingFlows with low resolution studies to try and find the parameter regime that creates sink particles…
About blog posts
Blog posts should be used by folks to log their work for meetings etc… Rarely should any actual work or documentation be put directly into a blog. Instead try to find a home on the wiki for the information be it in a ticket, project page, discussion page, or in the general wiki documentation on the code and then link to the work from the blog (or just embed the work in the blog using the Image macro) Below are some examples of linking/embedding images from discussion topics, wiki pages, or tickets…
Discussion topics
[attachment:AdvanceTimes_1.jpg:discussion:topic/23 SampleImage] [[Image(discussion:topic/23:AdvanceTimes_1.jpg, 50%)]]
Wiki pages
[attachment:StrongScalingEfficiencies.png:wiki:OptimizingScrambler SampleImage] [[Image(wiki:OptimizingScrambler:StrongScalingEfficiencies.png, 50%)]]
Tickets
[attachment:bug.png:ticket:126 TicketAttachment] [[Image(ticket:126:bug.png, width=50%)]]
Played around with blog plugin
So this new blog plugin is pretty fun… I don't like that one has to create a blog shortname for each entry but I suggest that we all use our login name followed by the date in mmddyy format and then attach letters b, c, d if we make more than one post per day it should (for example this post is johannjc053111b since I made an entry earlier today. We could use more descriptive names, but since these are all in the same blog folder and since we will probably each make several entries each week - it is probably better to use the title to describe what the blog is about and keep the shortnames (URLs) simple. I don't like that we can't currently attach files to the blogs… But if we are working on a ticket, project page, or general documentation we can attach the images there and then link to them from the blog like so….
Discussion attachments
[attachment:AdvanceTimes_1.jpg:discussion:topic/23 SampleImage]
Or within the Image macro as so:
[[Image(discussion:topic/23:AdvanceTimes_1.jpg, 50%)]]
Wiki page attachments
[attachment:StrongScalingEfficiencies.png:wiki:OptimizingScrambler SampleImage]
Or within the Image macro as so:
[[Image(wiki:OptimizingScrambler:StrongScalingEfficiencies.png, 50%)]]
Ticket Attachments
[attachment:bug.png:ticket:126 TicketAttachment]
Or within the image macro as
[[Image(ticket:126:bug.png, width=50%)]]
Worked on Scrambler Paper
I've updated the scrambler paper with several new figures and attached it to the wiki page here ms.pdf
Worked on getting bear2fix on clover to read hdf5 files from bluegene
So I've been able to implement a hack in bear2fix on clover to access the runs from bluegene. This basically involved manipulating the byte ordering of the integers stored within the boxdata that was being read in. After considering the bit representation of 119 and 29 as well as the values it was reading from the chombo file (1996488704 and 486539264) it looked like the order of the bytes was actually reversed - not the order of the bits…
119 | TTTFTTTF | FFFFFFFF | FFFFFFFF | FFFFFFFF |
1996488704 | FFFFFFFF | FFFFFFFF | FFFFFFFF | TTTFTTTF |
29 | TFTTTFFF | FFFFFFFF | FFFFFFFF | FFFFFFFF |
486539264 | FFFFFFFF | FFFFFFFF | FFFFFFFF | TFTTTFFF |
I implemented a function that converts from bluegenes big endian byte order to clover's little endian byte order.
FUNCTION convertfrombluegene(a) INTEGER :: a, convertfrombluegene,i,j convertfrombluegene=0 DO i=0, 3 DO j=0,7 if (btest(a,8*i+j)) convertfrombluegene=ibset(convertfrombluegene, 8*(3-i)+j) END DO END DO END FUNCTION convertfrombluegene
Why is hdf5 having difficulty with these data types? The idea of using hdf5 is to avoid these type of conversion issues… In fact if you run h5dump on the file you can see that the fact that these integers are big endian is stored in the hdf5 file and that h5dump correctly reads there value. However since these integers are stored within a box container data type (which contains 6 integers because of chombo requirements) there values have to be read as a whole. (This is technically only true for attributes and not for datasets but there are box type attributes for which this is still a problem).
Possible solutions:
- It looks like it may be possible to have hdf5 always convert values to little-endian before storing them - however this would cause problems if we ever wanted to do post processing of those data files on a big endian machine.
- One could have bear2fix ask the user if the byte ordering of the box data needs to be reversed but this might be confusing to the naive user.
- Perhaps the best solution would be to see if there is a way bear2fix to determine whether or not these box sub types are little endian or big endian and then check to see if it is running on a big endian or little endian machine.
- (the latter is easy to do by just looking at the first bit of the integer 1 as it will be true for little endian machines)
- The former could be :
- guessed by looking at the last 16 bits to see if there are any 1's… But this would only work for datasets with dimensions less than 2^{16 which is the current case - but won't be for much longer… }
- determined by having bear2fix probe the global.data file to get the global mX which it could then compare with the root level prob_domain attribute - if they didn't match it would know to convert the box data…
- finding an hdf5 command that returns the data type of the box data types component fields… Since visit is able to properly read these data files on clover I assume that this is possible… although I'm not sure if visit needs any of the box data attributes or just the box data datasets.
So - after fixing this and looking at the output chombo file it was apparent that something was still wrong… A little more digging showed that the output of the box data was being byte reversed when creating new fixed grid chombos - when the data type was being recorded in the chombo file as little endian integers… So I had to add an lConvert global logical variable that was only turned on when reading in chombo files - and was turned off when outputting them. After fixing this I was still seeing crazy values for dx and xupper… It appears that the float attributes were also being read in reverse order. I wrote another routine for converting floats from big endian to little endian order:
FUNCTION convertfrombluegenefloat(a) REAL(8) :: a, convertfrombluegenefloat INTEGER(8) :: ai, converti INTEGER :: i,j converti=0 ai=transfer(a,ai) !Cast the 8 byte real as an 8 byte integer since btest and ibset require integers DO i=0, 8 DO j=0,7 if (btest(ai,8*i+j)) converti=ibset(converti, 8*(7-i)+j) END DO END DO convertfrombluegenefloat=transfer(converti,convertfrombluegenefloat) END FUNCTION convertfrombluegenefloat