test

Investigated role of density (similar to magnetic field)

Fiducial

Spitzer

spitzer (.25 density)

spitzer (.2 density)

spitzer (.15 density)

spitzer (.1 density)

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

\eta = 1.15 \times 10^{-14} Z \log \Lambda(ne, T, Z, \mu) \left(\frac{k_BT}{eV}\right)^{-3/2} \mbox{s}

And then the magnetic diffusivity (with units of cm²/s) as

D=\frac{c^2}{4 \pi} \eta

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 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 ( \beta << 1 or M >> 1 ), 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 ( M >> 1 and M_A >> 1 ) the computational time is independent of the flow speed (only a function of resolution N ). 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 \frac{1}{\min \left [M, M_A \right]} (ignoring the extra factor of N^D due to changes in the number of zones with resolution)

This can be combined as

T_c \propto \frac{1}{\min \left [1, M, M_A \right ]}

This makes simulating the RT instability M << 1 relatively computational expensive. Alfven waves can also restrict the time stepping when

M_A = \beta^{1/2} M < 1.

So simulations with modest M but very small \beta 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 \frac{\Delta x^2}{\eta} and the cell crossing time \frac{\Delta x}{U}. The cell diffusion time will be smaller than the cell crossing time when

\frac{\Delta x^2}{\eta} << \frac{\Delta x}{U} \rightarrow R_M < N

So the computational time to simulate a crossing time will go as

T_c = \frac{\max \left[1, \frac{N}{R_M} \right]}{\min \left [1, M, M_A \right]}

So - you pay a penalty with explicit time stepping when R_M < N 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 R_M) - so in the explicit case we have

T_c = \frac{ \max \left[1, \frac{N}{R_M} \right]}{\min \left [1, M, M_A \right]} \times \max \left[1, R_M \right] = \frac{\max \left[R_M, N, \frac{N}{R_M} \right]}{\min \left [1, M, M_A \right]}

Making the resistive solve implicit would reduce this to

T_c = \frac{\max \left[1, R_M \right]}{\min \left [1, M, M_A \right]}

Added some documentation on MHD in AstroBEAR (including Hall terms)

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

\Delta T = \lambda \Delta F

where the climate sensitivity \lambda \approx .8 K/(W/m^2)

and we can approximate the logarithmic response of radiative forcing to changes in CO2

\Delta F = 5.35 \log \left ( \frac{P}{P_0} \right)

Combining these we get \Delta T = 4.3 K \log \left ( \frac{P}{P_0} \right)

or

T = T_0 + 4.3 K \log \left (\frac{P}{P_0} \right)

or

P = P_0 e^\frac{T-T_0}{4.3 K}

and \frac{dT}{dP} = \frac{4.3 K}{P} = \frac{4.3 K}{P_0} e^{\frac{-\left(T-T_0\right)}{4.3 K}}

So the climate sensitivity is not constant - but falls off by a factor of e every 4.3K change in temperature.

Here are lineouts from 4 of the runs with Bz field only. Note M=1.5 and \beta = .2 Also - the time reported should be scaled by 12.3 - so the final time is 426 ns

Each picture shows lineouts at three locations shown here

Each plot shows

• ram + thermal pressure \rho v^2 + P
• magnetic pressure B^2/2
• total pressure \rho v^2 + P + B^2/2
• density
• magnetic field
• velocity (in x)
• Temperature

And the values are scaled as follows

Click for movies

	No Cooling	Cooling
No Diffusion
Diffusion

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.

Diffusion	Cooling	No Cooling
No Diffusion
Tiny Diffusion
Diffusion
Diffusion + No diffusion in wire
Diffusion + 10% By	Cooling	No Cooling
Diffusion + Tiny Wire Diffusion + 10% By
Diffusion + No Wire Diffusion +10% By

testing

test

abc

To solve the ionization equations with scattering, you would want to use a diffusion equation. Something like

\frac{\partial I}{\partial t} - \nabla \cdot \frac{c \lambda}{\sigma n_H} \nabla I=-\sigma n_H c I

\frac{\partial n_H}{\partial t} = -\sigma n_H c I + \alpha n_en_H

\frac{\partial e}{\partial t} = \Lambda \sigma n_H c I

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).

​Reynolds et al

Ran suite of simulations of anisotropic Parker type winds with magnetic fields

Strong field at 45 degrees

​https://astrobear.pas.rochester.edu/trac/wiki/u/johannjc/injection

Model with modest magnetic field

The color shows the neutral density in particles/cc. Also shown are the inner, outer, and tau=1 surface (dashed lines) given by John's python script along with the calculated tau=1 from the code (solid black line). Mach surface is shown in red

movie​

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

​https://astrobear.pas.rochester.edu/trac/wiki/u/johannjc/scratchpad22

Here is an image from visit showing the non-radial velocities (in red) along with pressure forces (in black) projected in the x-y plane.

And here is a lineout of density showing that the scale height is only a few cells across with 7 zones per planet radius

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.

I updated the developer guide with the various graphs showing the stencil components and their dependencies… See SweepScheme

• Galactic code not constrained by memory… Running 10⁵ particles

• Data from new Hot Jupiter Run with solar ionizing flux values..

• ​considering trajectories in lab frame
• Wire turbulence MHD run is at 157 and running well…
• Working on paper…

movie​

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

Equations 18-22 and surrounding text discuss various opacities…

​http://iopscience.iop.org/article/10.1086/520791/pdf;jsessionid=47EB6BA531954D009DF98114B22F802F.c3.iopscience.cld.iop.org

cross section for hydrogen hydrogen collisions can be estimated as 4 \pi R^2 where R is approximately the bohr radius… This gives

\sigma = 3.5 \times 10^{-16} cm^2

The MFP is then l = \frac{1}{n \sigma}

The material leaving the planet is at 1.93 \times 10^{9} cm^{-3} so the MFP is then 1e-7 AU << dx

​https://astrobear.pas.rochester.edu/trac/wiki/u/johannjc/hotjupiter

Blackbody luminosity

We have previously calculated the blackbody luminosity from the Intensity as

L_\nu=\displaystyle \int_S \int_0^{2\pi} \int_0^{\pi/2} I(T, \nu) \cos \theta \sin \theta d\theta d\phi dA

where

I(T,\nu)= \frac{2 h \nu^3}{c^2} \frac{1}{e^{\frac{h \nu}{k T}}-1}

which gives

L_\nu=\frac{8 \pi^2 h R^2}{c^2} \frac{\nu^3}{e^{\frac{\nu}{\nu_S}}-1}

where

\nu_S = \frac{k T}{h}

and L has units of ergs/s/Hz

We 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²/A

After multiplying by 4*pi*47.1 pc ² to get Luminosity per angstrom, we need to convert to Luminostiy per Hz

To go from L_\lambda to L_\nu

we need to multiply by

L_\nu=L_\lambda \frac{d \lambda}{d \nu} = L_\lambda \frac{\lambda}{\nu}

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 -

\sigma_H(v) = \sigma_L \left ( \nu/ \nu_L \right )^{-3}

where

\sigma_L = 6.3 \times 10^{-18} \mbox{cm}^2

and

\nu_L=\frac{13.6 \mbox{eV} }{h} = 3.29e15 Hz

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)

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

The prescription of ​Bourrier et al is to calculate (at each pixel) the average optical depth \tau_i for range of velocity bins with spacing \Delta v (or frequency bins with size \Delta \nu, so that

\sum \tau_i \Delta \nu = \tau

So we have

\tau_i \Delta \nu = \tau_i \frac{\Delta v}{\lambda_0} = \displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \tau (v) dv

and

\tau (v) = \sigma_{v0} \int N(v') \phi ( v-v' ) dv'

Now to get the line profile in terms of velocities, we can start with the line profile in terms of frequencies:

\phi ( \nu - \nu' ) = \frac{\Gamma}{4 \pi^2 (\nu - \nu')^2 + \left ( \Gamma / 2 \right )^2}

and realizing that this is a distribution, we can use

\phi ( v - v') dv' = \phi ( \nu - \nu') d \nu

where

\nu = \nu_0 + \nu_0 \frac{v}{c}

so

\nu - \nu' = \frac{v - v'}{\lambda_0}

and

d \nu = \frac{dv}{\lambda_0}

to get

\phi ( v-v' ) =\frac{\Gamma \lambda_0}{4 \pi^2 (v - v')^2 + \left ( \Gamma \lambda_0 / 2 \right )^2}

Putting this all together we have…

\tau_i = \frac{\lambda_0}{\Delta v}\displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \int N(v') \frac{\Gamma \lambda_0}{4 \pi^2 (v - v')^2 + \left ( \Gamma \lambda_0 / 2 \right )^2} dv' dv

Now if we bin the column densities into velocity bins, we have

N_j = \displaystyle \int_{v_j-\frac{1}{2}\Delta v}^{v_j + \frac{1}{2}\Delta v} N (v') dv'

and

\int N(v') \frac{\Gamma \lambda_0}{4 \pi^2 (v - v')^2 + \left ( \Gamma \lambda_0 / 2 \right )^2} dv' = \sum N_j \frac{\Gamma \lambda_0}{4 \pi^2 (v - v_j)^2 + \left ( \Gamma \lambda_0 / 2 \right )^2}

and we have

\tau_i = \frac{\lambda_0}{\Delta v} \sigma_{v0} \sum N_j \displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \frac{\Gamma \lambda_0}{4 \pi^2 (v - v_j)^2 + \left ( \Gamma \lambda_0 / 2 \right )^2} dv

or

\tau_i = \frac{\lambda_0}{\Delta v} \sigma_{v0} \sum N_j C_{ij}

where

C_{ij}=\displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \frac{\Gamma \lambda_0}{4 \pi^2 (v - v_j)^2 + \left ( \Gamma \lambda_0 / 2 \right )^2} dv

Now \Delta v is of order 20 km/s and \Gamma \lambda_0 is 76 m/s… so if i == j, then the integral may as well go to infinity and we have

C_{ii} = 1

and for i \ne j (wave arms in air for a bit)

C_{ij} = \frac{\Gamma \lambda_0}{4 \pi^2 (v_i-v_j)^2} \Delta v

so we have

\tau_i = \frac{\sigma_{v0} N_i \lambda_0}{\Delta v} + \displaystyle \sum_{i \ne j} N_j \sigma_{v0} \frac{\Gamma \lambda_0^2}{4 \pi^2 (v_i-v_j)^2}

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 \kappa_i where

\kappa_i(n,v) = \delta_{ij}\frac{\sigma_{v0} n X \lambda_0}{\Delta v} + \left ( 1 - \delta_{ij} \right ) n \sigma_{v0} \frac{\Gamma \lambda_0^2}{4 \pi^2 (v_i-v_j)^2}

where j corresponds to the bin that the LOS cell velocity v lies within and X is the ionization fraction computing using either the ​Saha equation - or balancing ​radiative recombination with ​direct ionization… They give very different answers for the ionization fraction at the planet wind temp.

Remaining questions

1. How to chose the ionization fraction? - Need to balance photo ionization with radiative recombination
2. We should use thermal broadening (which dominates over natural broadening)..
3. Binning particle velocities before convolution seems numerically expedient - but binning errors are orders larger than convolution corrections…?
4. 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…

\tau_i = \frac{\lambda_0}{\Delta v}\displaystyle \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \int N(v') \phi(v-v') dv' dv

Now without thermal broadening, we have

N(v') = \sum n_j \delta(v'-v_j) \Delta z

however, with thermal broadening, we have

N(v') = \sum n_j \frac{1}{\sqrt{\pi} v_{th}} e^{-(v'-v_j)^2/v_{th}^2} \Delta z

where

v_{th}=\sqrt{\frac{2k_B T}{m}}

which gives us

\tau_i = \frac{\lambda_0}{\Delta v}\displaystyle \sum n_j \Delta z \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \int \frac{1}{\sqrt{\pi} v_{th}} e^{-(v'-v_j)^2/v_{th}^2} \phi(v-v') dv' dv

Now the integral over v' 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…

\phi(v-v')=\delta(v-v')

which gives

\tau_i = \frac{\lambda_0}{\Delta v}\displaystyle \sum n_j \Delta z \int_{v_i-\frac{1}{2}\Delta v}^{v_i + \frac{1}{2}\Delta v} \sigma_{v0} \frac{1}{\sqrt{\pi} v_{th}} e^{-(v-v_j)^2/v_{th}^2} dv

or upon integrating the exponential, we get

\tau_i = \frac{\sigma_{v0} \lambda_0}{2\Delta v}\displaystyle \sum n_j \Delta z \left [ \mbox{erf}\left ( \frac{v_j-v_i+\frac{1}{2}\Delta v)}{v_{th}} \right ) - \mbox{erf} \left ( \frac{v_j - v_i - \frac{1}{2}\Delta v)}{v_{th}} \right ) \right ]

So for each frequency bin, we can create a LOS integration of a corresponding integrated opacity

\kappa_i=\frac{n \sigma_{v0} \lambda_0}{2\Delta v} \left [ \mbox{erf}\left ( \frac{v_j-v_i+\frac{1}{2}\Delta v)}{v_{th}} \right ) - \mbox{erf} \left ( \frac{v_j - v_i - \frac{1}{2}\Delta v)}{v_{th}} \right ) \right ]

where n is the number density of neutral Hydrogen

Photoionization

Also need to update code to include photoionization rates as well as direct ionization and radiative recombination…

Photoionization rate can be calculated from

A=n_{HI} \displaystyle \int_{\nu_L}^\infty \frac{L_\nu \sigma_H(\nu)}{4 \pi r^2 h \nu} d\nu

We can approximate

\sigma_H(v) = \sigma_L \left ( \nu/ \nu_L \right )^{-3}

where

\sigma_L = 6.3 \times 10^{-18} \mbox{cm}^2

and

\nu_L=\frac{13.6 \mbox{eV} }{h}

and we can calculate the luminosity from Planck's law integrated over the area and solid angle…

L_\nu=\displaystyle \int_S \int_0^{2\pi} \int_0^{\pi/2} I(T, \nu) \cos \theta \sin \theta d\theta d\phi dA

L_\nu=4 \pi R^2 \left (\pi \right ) \frac{2 h \nu^3}{c^2} \frac{1}{e^{\frac{h \nu}{k T}}-1}

L_\nu=4 \pi R^2 \left (\pi \right ) \frac{2 h \nu^3}{c^2} \frac{1}{e^{\frac{\nu}{\nu_S}}-1}

L_\nu=\frac{8 \pi^2 h R^2}{c^2} \frac{\nu^3}{e^{\frac{\nu}{\nu_S}}-1}

where

\nu_S = \frac{k T}{h}

and

\lambda_S = \frac{c}{\nu_S}

Putting this all together, we get

A=2 \pi n_{HI} \nu_S \frac{\nu_L^3}{\nu_S^3} \frac{R^2}{r^2} \frac{\sigma_L}{\lambda_S^2 } \displaystyle \int_{\frac{v_L}{v_S}}^\infty \frac{1}{x} \frac{1}{e^x-1} dx

Now assuming a temperature of 5778 K, this gives

A = \frac{\alpha n_{HI}}{r^2}

where

\alpha = 3.6774 \times 10^{17} \mbox{cm}^2/\mbox{s}

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³

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

​Bzowski et al 2008

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 1.7165405 \times 10^{19} cm^2/s / r^2

currently I am using 3.5808d19

Consider the 2D hydro linearization of the wave equation in each dimension

\partial_t W_i = A_{ij} \partial_x W_j

where A_{ij} = \frac{\partial F_i}{\partial W_j}

The matrix A has eigenvalues u-a, u, u, u+a corresponding to left propagating sound wave, contact, shear discontinuity, and right moving sound wave. And we can construct the matrix of left eigenvalues L (where each row is a left eigenvector)

and the matrix of right eigenvectors R (where each column is a right eigenvector)

and the diagonal matrix of eigenvalues \Lambda where each diagonal term is an eigenvalue

The eigenvectors are orthogonal

\mathbf{LR} = \mathbf{I}

and we also have

\mathbf{L{^-1}R^{-1}} = \mathbf{I}

\mathbf{RL} = \mathbf{I}

\mathbf{AR=R \Lambda}

\mathbf{LA = \Lambda L}

which combining gives

\mathbf{L A R = \Lambda}

\mathbf{A = L^{-1} \Lambda R^{-1}}

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

\partial_t \mathbf{W}=\mathbf{A}\partial_x\mathbf{W}

\partial_t \mathbf{W}=\mathbf{RLA}\partial_x\mathbf{W}

\partial_t \mathbf{W}=\mathbf{R\Lambda L}\partial_x\mathbf{W}

\partial_t \mathbf{W}=\sum \lambda_i \mathbf{r_i} \mathbf{l_i^T} \partial_x \mathbf{W}

Now if we renormalize the left eigenvectors

and the right eigenvectors

then

\partial_t \mathbf{W}=\sum \alpha_i \lambda_i \mathbf{r_i}

where

\alpha_i=\mathbf{l_i^T} \partial_x \mathbf{W}

has units of \frac{1}{dx}

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

Restart of HD2794b with additional refinement

MHD run consistently stalls while outputting frame 109. Output file appears to be 3x larger than previous frames. Must be a bug in the io routines.

​movie

And WireTurbulence run is at frame 106/200

• 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…
  • 
  • 
  • 

Working on implementing Anisotropic Thermal Diffusion

The 4 columns correspond to Isotropic diffusion, Anisotropic but with \chi_\perp = \chi_\parallel, 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, \lambda = 0, 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​

see u/johannjc/scratchpad6

movie​

Movie showing streamlines of the velocity field (magenta), magnetic field lines, the mach surface (black) and log density (pseudocolor) as the simulation relaxes after injecting a dipole field.

movie​

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

• Level 0 resolution is 1600x160x160
• Simulations started at level 2
• Backed off to level 1 and restarted MHD from beginning
• Hydro run continued from frame 62 but at level 1.

Level	Hydro	MHD
2	0-84	0-39
1	62-200	0-42

Hydro run at frame 200

​

movie​

—-

MHD run at frame 38

movie​

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

\frac{R_\Omega}{a} = \frac{v_{esc}}{a\Omega} = \sqrt{\frac{2GM_p}{r_p a^2}\frac{a^3}{G\left(M_p+M_s \right) }} = \sqrt{\frac{2aq}{r_p \left(q+1 \right)}}

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

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)

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…

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

v_{esc}=\sqrt{\frac{2G M_p}{R_p}}

and

\Omega=\sqrt{\frac{G \left ( M_s+M_p \right ) }{a^3}}

We can calculate an coriolis radius R_C where

2 \Omega v_{esc} = \frac{v_{esc}^2}{R_C}

which gives

R_C=\frac{v_{esc}}{2 \Omega}=\sqrt{\frac{M_p a^3}{2 \left ( M_s + M_p \right ) R_p}} \approx \sqrt{\frac{R_H^3}{R_p}}

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
  1. v_poloidal parallel to B_poloidal
  2. rho and P held constant
  3. v set to zero ( v_phi is corotating)
  4. 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

• Code now supports parallel IO
• Cleaned up unused variables

Consider the magnetic field to be similar to an atmosphere where we balance magnetic pressure gradients with magnetic tension forces? Than we have

\frac{d B^2}{2 d r} = -\frac{B^2}{r}

which gives us

B^2=B_0^2\frac{r^2_{o}}{r^2}

We can then integrate this to get the net tension force

F_T=h\displaystyle \int_{r_i}^{r_o}{\frac{B(r)^2}{r} 2\pi r dr} = \pi h B_0^2 r_{o}^2 \left ( \frac{1}{r_i} - \frac{1}{r_o} \right )

and dividing by the inner area gives

P_T=B_0^2 \frac{r_o \Delta r}{r_i^2} = P_{ram}

which we can combine with flux freezing

2 \pi B_0 r_o \Delta r = B_0 \pi \left ( r_o^2 - R^2 \right )

and

\Delta r = r_o-r_i

Which can be reduced to

2 r_o \Delta r = \left ( r_o^2 - R^2 \right )

B_0^2 \frac{r_o \Delta r}{\left ( r_o-\Delta r \right )^2} = P_{ram}

And then can be combined to

\frac{2 P_{ram}}{B_0^2}\left ( r_o-\Delta r \right )^2 = \left ( r_o^2 - R^2 \right )

\frac{\Delta r}{r_o} = \frac{1}{2} \left ( 1 - \frac{R}{r_o}^2 \right )

\frac{\beta_{ram}}{4} \left ( R^2 + r_o^2 \right ) = \left ( r_o^2 - R^2 \right )

which eventually gives (if I haven't made any mistakes…)

r_0 = R \sqrt{\frac{4+\beta_{ram}}{4-\beta_{ram}}}

although clearly I have since \beta_{ram} > 4

I will start with a definition of ram pressure for a cylindrical flow

P_{ram}(r)=\rho_0 (r_o/r)^2 c^2

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

B_0 S_0 = B_1 S_1

where S refers to the area the flux threads.

This gives us

B_0 \pi r_0 \delta r_0 = B_1 \pi r_1 \delta r_1

If we now use the tension "pressure" going as

P_t = \frac{B(r)^2}{r}

then we have ram pressure = magnetic tension leading to

B_0^2 r_0^2 r^{-3} = \rho_0 (r_o/r)^2 c^2

or the balance radius r_b

r_b^3 = (\frac{B_0^2 r_0}{\rho_0 c^2})^{1/2}

Planetary magnetic fields are dipole fields of the form

\frac{r_0^3}{r^5} \left ( 3\mathbf{r}\left ( \mathbf{B_0} \cdot \mathbf{r} \right) - \mathbf{B_0} r^2 \right )

or

B_{\theta}=-B_0 \cos{\theta}\frac{r_0^3}{r^3}

and

B_r=2 B_0 \sin{\theta}\frac{r_0^3}{r^3}

which (if we assume that the dipole is oriented along z) gives a field of

2 B_0 \hat{z}

at the pole and a field of

-B_0 \hat{z}

at the equator.

This also corresponds to a vector potential

\mathbf{A}=r_0^3 \frac{\mathbf{B_0} \times \mathbf{r}}{r^3}

Now if we assume that

\mathbf{B_0} = B_0\hat{z}

Then our vector potential is

\mathbf{A}=r_0^3 \frac{B_0 \sin \theta}{r^2} \hat{\phi}

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 16 B_0 \hat{z} at the center would require a potential of the form

\mathbf{A}=\frac{16 B_0 r \sin{\theta}}{2} \hat{\phi}

which agrees with the potential outside at r=r_0/2

So we don't need any smoothing… just a piecewise definition for the potential that switches at r=r_0/2

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

A_\phi^*=r_0^3\frac{B_0 \sin{\theta}}{r^2} \left (1 -\frac{r}{r_s} \right ) ^2

For r > r_s the potential is just 0.

This gives

B_{\theta}^*=\frac{1}{r}\frac{\partial}{\partial r}\left ( r A_\phi \right ) = -B_0 \sin{\theta}\frac{r_0^3}{r^3} \left ( 1 - \frac{r^2}{r_s^2} \right )

and

B_r^*=\frac{1}{r \sin{\theta}}\frac{\partial}{\partial \theta}\left ( \sin{\theta} A_\phi \right ) = 2B_0 \cos{\theta}\frac{r_0^3}{r^3} \left ( 1 - \frac{r}{r_s} \right )^2

for r < r_s and \mathbf{B}=\mathbf{0} for r>r_s

	Smooth	UnSmoothed
log
linear

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

\psi=v/c_s

when in fact

\psi=v^2/c_s^2

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);

Here is a snapshot of frame 51 of 200 of the hydro run and frame 34/200 of the MHD run.

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…

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

M'=M+\sum{\Delta m_i}

M'V'=MV+\sum{\Delta m_i v_i}

M'R' = M'R + \sum{\Delta m_i \left ( r_i - R \right )}=MR + \sum{\Delta m_i r_i }

But our treatment of angular momentum has glossed over some subtleties

S'=S+\sum{\Delta m_i \left ( r_i -R \right ) \times v_i}

A better treatment is to conserve total angular momentum.

J=L+S

where

L=MR \times V

and S is the spin.

J'=J+\sum{\Delta j_i} = J +\sum{\Delta m_i r_i \times v_i}

S'=J'-L'

S'=S + MR \times V - M'R' \times V' + \sum{\Delta m_i r_i \times v_i}

S' = S+MR \times V - \left( MR + \sum{\Delta m_i r_i} \right ) \times \frac{MV + \sum{\Delta m_i v_i}}{M+\sum{m_i}}+\sum{\Delta m_i r_i \times v_i}

which if we approximate to first order in \frac{\sum{\Delta m_i}}{M}

S' \approx S+MR \times V - \left( MR + \sum{\Delta m_i r_i} \right ) \times \left ( V + \frac{\sum{\Delta m_i v_i}}{M} \right ) \left (1 -\frac{\sum{\Delta m_i}}{M} \right )+\sum{\Delta m_i r_i \times v_i}

S' \approx S+MR \times V - MR \times V - \sum{\Delta m_i r_i } \times V - R \times \sum{\Delta m_i v_i} + R \times \sum{\Delta m_i} V+\sum{\Delta m_i r_i \times v_i}

S' \approx S +\sum{\Delta m_i \left (r_i - R \right ) \times \left (v_i - V \right )}

which implies our method is only accurate if V << v_i and \sum{\Delta m_i} << M

2D slice movie​

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.

2D version with stepped on boundaries

2D version with outflow only wind at top (semi-extrapolating)

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.

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³w² = 5726 N_r³ = (18 N_r)³

So if we want 100 zones in radius, we would be effectively simulating an 1800³ 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)

Worked on calculating convective instability and opacities. Found description of stellar convection ​here

Plot showing \frac{d \log T}{d \log P} as well as 1 - \frac{1}{\gamma}

Plot 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.

• Upper right panel is 2.5D with 2048² effective
• Lower right panel is 2.5 D with 256² effective
• left half is slice from 3D with 128² effective

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.

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

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

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

And yes - they all blow up

These are not all at the same time.

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

Modified temperature profile to avoid enormous density contrasts

Made a movie of the template module

movie​

Running wire turbulence, but shocked flow is not yet turbulent…

• Eddie's images
• Working on Parallel HDF5 routines - and cleaning up io/chombo stuff.
• Photo Ablated Clump post

Ran a simple photo ablated clump model in 2D - AMR - parallel

Click here for a movie​

​movie

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

\frac{||x-x'||}{||x||}= \mbox{cond}(A) \frac{||b-b'||}{||b||}

which would imply that \frac{||x-x'||}{||x||}= \mbox{cond}(A) \times \mbox{tol}

In general we have \mbox{cond}(A) = ||A|| \; ||A^{-1}||

and \frac{1}{\mbox{cond}(A)} 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

\frac{||x-x'||}{||x||} \approx 2 \mbox{cond}(A)\epsilon_M

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⁴, 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…

\max{\frac{|x_i-x'_i|}{|x_i|}}

New time stepping that is based on actual time scale for energy to change.

Movie​

test

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?

• 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.

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

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 profile using density profile and external pressure

Mathematically, we are solving the equation of hydrostatic equilibrium

\nabla P(r) = -\rho(r) \nabla \phi(r)

subject to the constraint that P(R)=P_0.

Since \phi(r) = -\frac{G M(r)}{r} we can rewrite the equation for HSE as

and we can express the enclosed mass as

Now we can discretize equation 2 as

where \Delta r_j = r_{j+1}-r_j

and equation 1 as

but a taylor series expansion shows that the error at each step in the integration goes as

\frac{\Delta r_j}{r_j} + \frac{\Delta \rho_j}{\rho_j} + ...

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

• ​Paper describing microphysics
• ​Paper describing implementation and analytic approximation

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

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

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

movie​

Making the momentum be conserved requires recasting the gravitational force as a total differential

\dot p=-\rho \nabla \phi = -\nabla \cdot F so we need to find F such that

\nabla \cdot F = \rho \nabla \phi

Since \nabla^2 \phi = 4\pi G (\rho-\bar{\rho})

we can substitute for \rho and we have

\nabla \cdot F = (\frac{\nabla^2\phi}{4 \pi G}+\bar{\rho}) \nabla \phi which is equivalent (in 1D) to \nabla \cdot \left [\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G} + \bar{\rho} \phi \right ]

where we can identify the equivalent momentum flux tensor as F=\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G} + \bar{\rho} \phi

• In more than 1D we have

or we can write it as \mathbf{T}=\frac{1}{4\pi G}\left [ \nabla \phi \nabla \phi - \frac{1}{2} \left ( \nabla \phi \right ) \cdot \left ( \nabla \phi \right ) \mathbf{I} \right ]

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

• 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…

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.

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 v_1,v_2,r_1,r_2.

The center of mass residing at the origin implies that m_1r_1=m_2r_2

and we have r_1+r_2 = a which can be solved for

r_1=a\frac{m_2}{M}

r_2=a\frac{m_1}{M}

We also must have zero net momentum so

m_1v_1=m_2v_2

which implies that

\frac{v_1}{r_1} = \frac{v_2}{r_2}

Which is consistent with them having the same angular velocity at aphelion.

Now h=\frac{L}{\mu} = \left( \frac{r_1 v_1}{m_2} + \frac{r_2 v_2}{m_1} \right ) M = a \left ( v_1 + v_2 \right )

And we know that twice the entire area of the ellipse is swept out in one orbital period

A=\pi a b = \pi a^2(1-e^2)

T=2 \pi \sqrt{\frac{a^3}{GM}}

we have

h=a \left ( v_1+v_2 \right ) = \frac{2A}{T}=\sqrt{a(1-e^2)GM}

or

v_1+v_2=\sqrt{\frac{GM(1-e^2)}{a}}

which then gives us

v_1=\frac{\sqrt{\frac{GM(1-e^2)}{a}}}{1+\frac{m_1}{m_2}}

v_2=\frac{\sqrt{\frac{GM(1-e^2)}{a}}}{1+\frac{m_2}{m_1}}

which together with

r_1=a\frac{m_2}{M}

r_2=a\frac{m_1}{M}

define

\mathbf{x_1}=[r_1,0,0] \mathbf{v_1}=[0,v_1,0]

and

\mathbf{x_2}=[-r_2,0,0] \mathbf{v_2}=[0,-v_2,0]

Just because it is pretty

ideal gas equations

p=nkT

e=\frac{p}{\gamma-1}

and definitions of sound speed and enthalpy

c_s = \left ( \frac{\partial P}{\partial \rho} \right )_s = \sqrt{\gamma P}{\rho}

h=e+p

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

• 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)

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

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.

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)

We now have excellent agreement with Krumholz results after modifying boundary conditions to match Su Olson, and correcting the ratio of Gas to Radiation energy density…

So Marshak boundary conditions are mixed

\left . \left ( E-\frac{2}{3\kappa}\nabla{E} \right ) \right |_0 = \frac{4}{c} F

One could instead do

E=\frac{4}{c}F

or

\frac{c}{3\kappa} \nabla E = F

so that

\frac{\partial E}{\partial t} = \nabla \cdot \left ( \frac{c}{3\kappa} \nabla E \right ) = \nabla \cdot F

In any even, if we discretize this mixed equation we get (E_g is ghost zone, E_i is internal zone at left edge)

E_g + E_i - \frac{4}{3\kappa \Delta x} \left ( E_i - E_g \right ) = \frac{8}{c} F

which we can solve for E_g as

E_g = \frac{\left ( 1 - \frac{3\kappa \Delta x}{4} \right ) E_i + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}}

and we have E_g - E_i = \frac{\left ( -\frac{6\kappa \Delta x}{4} \right ) E_i + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}}

If we then plug this into the time evolution equation ignoring portions that are unchanged we get:

E^{n+1}_i-E^{n}_i = \alpha \left ( E_g^*-E_i^* \right )

where \alpha = \frac{c \Delta t}{3 \kappa \Delta x^2} we get

E^{n+1}_i-E^{n}_i = \alpha \left ( \frac{-\frac{6\kappa \Delta x}{4} E_i^* + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}} \right )

where E^* = \psi E^{n+1} + \bar{\psi} E^n is the time averaged energy and depends on the time stepping (crank-nicholson or forward euler)

and finally we arrive at

\left ( 1+\frac{2 \psi \alpha}{1+\frac{4}{3\kappa\Delta x}} \right ) E^{n+1}_i = \left ( 1-\frac{2 \bar{\psi} \alpha}{1+\frac{4}{3\kappa\Delta x}} \right ) E^{n}_i+ \frac{\frac{8\alpha}{c} F}{1+\frac{4}{3\kappa \Delta x}} This is very similar to what we normally get ie

\left ( 1+\psi \alpha \right ) E^{n+1}_i - \psi \alpha E^{n+1}_g = \left ( 1-\bar{\psi} \alpha \right ) E^{n}_i + \bar{\psi} \alpha E^{n}_g

and noting that \alpha = \frac{c \Delta t}{3 \kappa \Delta x^2} for the limiter \lambda=1/3

or flux is just

F \frac{\Delta t}{\Delta x} \rightarrow F \frac{\Delta t}{\Delta x} \frac{2}{1+\frac{3\kappa\Delta x}{4}}

Note that E^{n+1}_i = E^{n}_i for

E = \frac{4F}{c}

and that if 3 \kappa \Delta x = 4 we revert to what we would expect if there was no Radiation energy in the ghost zone.

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

• 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.

So we start with the linearized differential equations

which we Fourier transform in physical space and solve the third equation for \phi 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

\lambda^2 + k^2 c_s^2 -4 \pi G \rho_0 = 0

and we have

\lambda = \pm \sqrt{4 \pi G \rho_0-k^2 c_s^2}

with eigen vectors

However given the presence of self-gravity we see that \lambda can be real as well as imaginary.

So for stable waves we have \lambda = \pm i\omega where \omega is real and positive. And there are two solutions… (a left and right going gravity-modified sound wave. Gravo-accoustic modes)

And if \lambda=\pm\omega where \omega 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.

I ran a 3D 16^{3 field loop advection test on bluestreak one 1 node using {16,32,64} total threads between {1,2,4,8,16,32*,64*} mpi tasks}

The Implicit radiation transport solver is working in parallel for fixed grid.

movie​

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.

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 > 

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 a=b

note that the old latex macro still renders images on the server

[[latex($a=b$)]]

gives a=b

so the only problem will occur if people have used

[[latex($a=b$)]]

a=b

(or removed their dollar signs when we switched from tracmath to tracmathjax)

• Working on Blue Gene Optimization
• Merged development branch with radiative transfer code and implementing explicit radiation source terms

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)

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 field loop advection run, the single node threaded version was 20% faster than the single node mpi version

• Compared astrobear with loop and sweep direction reordering
  • 10% speedup on 2D 16²+4
  • 20% speedup on 3D 64³+0
• Next going to look at ways of using shared memory efficiently.

• 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...)

And here is the pdf​

• 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'
  

• 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.

• MHD Mixed Results
• Radiation Transfer coded, ready for testing.

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 B^x_{i+1/2,j,k} = B^x_{i+1/2,j,k} \left ( 1 - 4\alpha \right ) + \alpha \left ( B^x_{i+1/2,j+1,k} + B^x_{i+1/2,j-1,k} +B^x_{i+1/2,j,k+1}+B^x_{i+1/2,j,k-1} \right ) + \alpha \left ( B^y_{i+1,j+1/2,k}-B^y_{i,j+1/2,k} - B^y_{i+1,j-1/2,k}+B^y_{i,j-1/2,k} \right ) + \alpha \left (B^z_{i+1,j,k+1/2}-B^z_{i,j,k+1/2} - B^z_{i+1,j,k-1/2}+B^z_{i,j,k-1/2} \right )

Now if we are only worried about Bx, and smoothing in y and z we can simplify this to B^x_{i+1/2,j,k} = B^x_{i+1/2,j,k} \left ( 1 - 4\alpha \right ) + \alpha \left ( B^x_{i+1/2,j+1,k} + B^x_{i+1/2,j-1,k} +B^x_{i+1/2,j,k+1}+B^x_{i+1/2,j,k-1} \right )

q_{i,j,k}=q_{i,j,k} \times \left (1 - 4\alpha \right ) + \alpha \left (q_{i,j-1,k}+q_{i,j+1,k} + q_{i,j,k-1}+q_{i,j,k+1} \right )

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.}

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 A = \mbox{"There are 3 jelly beans"} and B = \mbox{"a jelly bean drawn at random is numbered 1"}

Now the probability of drawing the number 1 jelly bean IF there are 3 jelly beans is P(B|A) = 1/3 and the odds of there being 3 jelly beans (as opposed to 4) are P(A)=1/2 but what is P(B) ? Well, we have to add all the ways we could draw a '1'. P(B) = P(B|A)P(A) + P(B|C)P(C) where C=\mbox{"There are 4 jelly beans"}. So then we have

P(A|B) = \frac{P(B|A) P(A)}{P(B|A)P(A) + P(B|C)P(C)} = \frac{1/3 \times 1/2}{1/3 \times 1/2 + 1/4 \times 1/2} = 4/7 .

Now A was drawn from the discrete set of events [A,C], but we could imagine a continuous set of events (\x$)]]. Then Bayes theorem reduces to P(x|B) = \frac{P(B|x) P(x)}{\int{P(B|x') P(x') dx'}}

So 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 P(x) = 1/L . Now, P(B|x) = 1/x for x > 1 and 0 otherwise. So we have

P(x|B)=\frac{1/x \times 1/L}{\int\limits_{1}^{L} 1/x' dx' \times 1/L} = \frac{1/x}{\log(L)-\log(1)}

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 < x > \sim \frac{L}{\log(L)} and the 50th percentile will go as x_{50} \sim L^{1/2}

The above analysis is for a distribution in which there are equal probabilities for being found at any point along the line. P(B|x)=1/x , and an equal likelihood for the line being any length less than L P(x)=1/L , however we can also consider a system in which the probability of being at any given point grows exponentially with distance (or time). P(B|x) \propto exp(x) while still keeping P(x)=1/L

For 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 \exp^{Rt}

In this case we have P(B|x)=\frac{e^{1}}{e^x} and P(x)=\frac{1}{L} . So we have P(x|B)=\frac{e^{1-x} \times 1/L}{\int \limits_1^L e^{1-x'} dx' \times 1/L} = \frac{e^{-x}}{e^{-1}-e^{-L}}

This gives an expected value of < x > = 2 \frac{e^{-1}-(L+1)/2 e^{-L}}{e^{-1}-e^{-L}} and an x_{50}=1.517. In this case the expected value is fairly insensitive to L 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?

• Magnetized Colliding Flow run

• Information on numeric scales
• Numerics of Flux Limited Diffusion

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.

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₁₀₀ ~ 20%. The odds that there are either 1 or 2 jelly beans would then be (1/1 + ½) / H₁₀₀ = H₂/H₁₀₀ ~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') …

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²)
    • 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

Finished the documentation on FluxLimitedDiffusion in the Developer's Guide

New ​movie of Magnetized Colliding Flows

• 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

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 \gamma=\frac{\sum{\gamma_i \rho_i}}{\sum{\rho_i}}

but it turns out that the value that should be weighted is the degrees of freedom

\bar{f}=\frac{\sum{f_i n_i}}{\sum{n_i}}

There is a relation between the degrees of freedom and the adiabatic index

\frac{f}{2}=\frac{1}{\gamma-1}

so an effective gamma can be derived by taking

\bar{\gamma}=\frac{\sum{n_i}}{\sum{\frac{n_i}{\gamma_i-1}}}+1

• 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

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….

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²)

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⁴ due to helium excitation, the majority of the cooling comes from metals excitation by electrons (especially Oxygen at 10⁵). 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:

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 PH_{cool}=n_H^2 PH(ne, \chi, T) while the DM metal curve is only a function of T DM_{cool}=n_H n_e DM(T)

So I divided Pat's table by \chi to get

PH'(n_e, \xi, T)=PH(n_e, \chi, T)/\chi

so that n_H n_e PH'=n_H^2 PH 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 n_e=\chi n_{Hi} 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:

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)

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

• 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

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 P(\rho,T) = \frac{\rho RT}{1-\rho b}-a\rho^2 However to solve the euler equations we need a caloric equation of state

So we need the internal energy density (which can't be derived from the pressure function alone) e=\frac{f}{2} \rho RT-a\rho^2 where f is the degrees of freedom.

This allows us to write our caloric EOS e(p,\rho)=\frac{f(p+a\rho^2)(1-\rho b)}{2}-a\rho^2

as well as solve for the pressure given the internal energy density p(e,\rho)=\frac{2(e+a\rho^2)}{f(1-\rho b)}-a\rho^2

so we can calculate the specific enthalpy

h(p,\rho)=(e+p)/\rho=\frac{f(p+a\rho^2)(1-\rho b)}{2\rho}-a\rho+\frac{p}{\rho}

and then use the specific enthalpy to calculate the sound speed using various thermodynamic identities…

c^2=-\frac{\left (\frac{\partial h}{\partial \rho}\right )_p}{\left ( \frac{\partial h}{\partial p} \right) _\rho - \frac{1}{\rho}}

c^2=-\frac{f\left [ 2a\rho^2(1-\rho b)-b\rho(p+a\rho^2) -(p+a\rho^2)(1-\rho b) \right ] - 2\left(p+a\rho^2 \right)}{f\rho(1-\rho b)}

c^2=\frac{f+2}{f(1-\rho b)}\left(\frac{p}{\rho}+a\rho\right) - 2 a \rho

which when a = 0 reverts to

c^2=\frac{f+2}{f(1-\rho b)}\frac{p}{\rho}=\gamma\frac{p}{\rho}

and when b = 0 reverts to

c^2=\frac{f+2}{f}\frac{p}{\rho} - \frac{f-2}{f}a \rho

which gives damping sound waves when

p < \frac{f-2}{f+2}a\rho^2

and when a=b=0 reverts to

c^2=\frac{f+2}{f}\frac{p}{\rho}=\gamma \frac{p}{\rho}

and finally we need to invert the equation for the specific enthalpy

p(h,\rho)=\frac{2\rho h + 2a\rho^2-f(a\rho^2)(1-\rho b)}{(2+f(1-\rho b))}

So we ran the basic shock tube tests comparing ideal gas and VanDerWaals gas.

Note that as \rho \rightarrow \frac{1}{b} then pressure should begin to dramatically increase with density due to limited volume.

And if \rho \rightarrow \frac{RT}{a} then attractive Van der Waals forces begin to significantly modify the pressure.

The parameters in cgs were

Note since \rho_L=10 \mbox{ so } \rho_L b = .3258 \approx 1

And \rho_L << \frac{RT}{a} 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).

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

\gamma*_s=\gamma_s + 2 \left( 1 - \frac{\hat{\gamma}}{\hat{\Gamma}} \right ) \left ( \hat{\gamma} - 1 \right ) \frac{P* - P_s}{P*+P_s}

where

\hat{\gamma}=\frac{1}{2} \left ( \gamma_L + \gamma_R \right ) \hat{\Gamma}=\frac{1}{2} \left ( \Gamma_L + \Gamma_R \right )

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

(S_R-S_L)U^* - (S_R U_R - S_L U_L) + (F_R - F_L) = 0

And then integrating around the right half we have

S_R U^* - (S_R U_R) + (F_R - F^*) = 0

Or the left half

-S_L U^* - (-S_L U_L) + (F^* - F_L) = 0

which we can solve F*

F^* = \frac{S_R F_L - S_L F_R + S_RS_L(U_R-U_L)}{S_R-S_L}

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

S_M=\frac{\left(\rho u\right)*}{\rho*}

and then the jump conditions can be solved for p* from either the right or left as:

p^*=p_L+\rho_L(S_L-u_L)(S_M-u_L) p^*=p_R+\rho_R(S_R-u_R)(S_M-u_R)

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.

• 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

Working on the MHD documentation and decided to create a nifty diagram

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… \rho=\rho_\infty e^{-\phi/\left( P_\infty / \rho_\infty \right )} T=P_\infty/\rho_\infty
• 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 P_\infty to give pressure matching at the disk edge we end up with a transcendental equation
• P(R_{disk})=P_\infty e^{-\phi(R)/\left( P_\infty/\rho_\infty \right )}
• If instead we fix T or P_\infty / \rho_\infty and solve for \rho_\infty we get an easier equation to solve.
• P(R_{disk})=\rho_{\infty}(P_\infty/\rho_\infty) e^{-\phi(R)/\left( P_\infty/\rho_\infty \right )}

• \rho_{\infty} = P(R_{disk})(\rho_\infty/P_\infty) e^{\phi(R)/\left( P_\infty/\rho_\infty \right )}

movie​

movie​

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

• 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 subroutine
  • Simulation 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 possible
  • Dimensions 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 code
  • Also 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 scales
  • Initializing a Grid says that each infodef has a %dx which is not true
  • Initializing a grid has an incorrect expression for the number of ghost zones
  • Flagging cells for refinement implies that the errflag array extends into ghost zones – which it doesn't
  • AdditionalPhysics should probably be a separate page
  • AdditionalPhysics 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 guide
• Meeting pages are depracated
• Collaborators & Projects, Publications, and Image Gallery probably don't belong under 'NewUsers'

• 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.

• 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 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³+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
  

• 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.

• Looked at BlueStreak performance
• Made new figures for ​paper
• Working on Christina's run #267
• Preparing talk for Berkeley on 12/7…

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.

• 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

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

• 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

• worked on Martin's paper
• Hydrostatic Star Problem
• Applied for Miller fellowship
• Working on Spectra for colliding flow papers

export PYTHONPATH=/data/trac/astrobear/plugins:$PYTHONPATH; easy_install --install-dir /data/trac/astrobear/plugins tracmypageplugin/0.11/

Nora Rose Carroll-Nellenback

Born September 21, 2012

weighing 11 lb 7 oz

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…

CurrentTests(johannjc, 250px)

(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…

Worked with Shule and Ivan on code implementation and prepped for journal club. Left for wedding on Friday

• Completed test modules and documentation for:
  • Filling Fraction
  • Marquee
  • Isotropic Turbulence
  • Gravitational Cascade
  • Molecular Cloud Formation
• Updated Matt's tests macro plugin to show links to test descriptions (see AstroBearTesting)

1. Link and build your problem module.
2. Modifying the data files to follow the format of those in Template
3. 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)
4. make a ref directory parallel to the out directory and copy over chombo00001.hdf
5. make a logs directory
6. make an images directory
7. Run bear2fix and select operations to make a nice image that can be used to verify the solution is 'correct'
8. run ../../bear2fix/ps2png.pl out/pgout00001.ps
9. mv out/pgout00001.png images/ref.png
10. Now rerun the same code with a different number of processors and pipe the output to logs/testlog
11. Rerun bear2fix in batch mode to produce another image for this run
12. run ../../bear2fix/ps2png.pl out/pgout00001.ps
13. mv out/pgout00001.png images/sim.png
14. mv bear2fix.data bear2fix.data.img
15. Run bear2fix and select operation 10
16. mv bear2fix.data bear2fix.data.test
17. mv out/chombo00001.hdf logs/sim.hdf
18. 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
19. add bear2fix.data.img and bear2fix.data.test to the mercurial repository (hg add bear2fix.data.img bear2fix.data.test)
20. 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!
21. Commit the new file additions and the modification of buildproblem to the repo.
22. Now you will want to use svn to check out a new version of the test repo. svn co ​file:///cloverdata/repositories/tests ~/tests
23. Create a directory for your problem module in the test repo parallel to Bondi, etc…
24. Navigate into your newly created directory and move the logs, images, and ref directory from your code directory
25. Now add your test problem directory to the svn repo using svn add MolecularCloudFormation (or whatever yours is called)
26. Now we can commit those additions to the svn repo. svn ci -m 'added my problem module'
27. 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.
28. 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
29. Now go back to your code directory and run buildproblem -np 8 -t to start the problem testing
30. If this passes then you are all set to check in your mercurial repo! If not then seek help
31. Don't forget to add your problem module to ​https://clover.pas.rochester.edu/trac/astrobear/wiki/TestSuite
32. And then create a page for your module using the TestTemplate template

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…

• Worked on static AMR improvements. #242
• Found a bug with self gravity and sink particles.

movie​

• 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

• 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

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.

First every distance is softened to avoid trouble at origin

Distance from origin: r=>\sqrt(r^2+\epsilon^2)

Distance from rotation axis: s=>\sqrt(r_s^2+\epsilon^2)

Disk radius: r_D=>\sqrt(r_D^2+\epsilon^2)

Then angular velocity is set to be keplerian and balance force along the 's' direction

\Omega=\sqrt{\frac{G M}{s^3}}

Then the density within the disk is set to be constant along the midplane \rho_0 with a decaying exponential to give vertical hydrostatic equilibrium.

\rho=\rho_0 \exp{ \left [ h \left ( \frac{1}{r}-\frac{1}{s} \right ) \right ]}

where the scale height is given by h=\frac{GM}{c_s^2}

At the outer rim of the disk when s > r_D Omega goes to zero and the profile changes slightly to

\rho=\rho_0 \exp{ \left [ h \left ( \frac{1}{r}-\frac{1}{r_D} \right ) \right ]}

which is just the profile for an envelope about a point mass in hydrostatic equilibrium ie \rho \propto \exp(1/r)

Then for each point \rho => \max(\rho, \rho_B) where \rho_B 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.

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.

• Found bug in reflecting self-gravity #228
• working on scaling - #202
• implementing new refinement criteria/implementation #225

​http://youtu.be/grhEYDd0HnA

(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

​http://youtu.be/Q4l5B32glbI

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 \phi = -1/r) and still satisfy poisson's equation \nabla^2 \phi = \rho > 0 since \nabla^2 \phi \ne \frac{\partial^2 \phi}{\partial^2 r} but rather \nabla^2 \phi = \frac{1}{r^2}\frac{\partial}{\partial r}r^2\frac{\partial \phi}{\partial r}

So the potential can go like r² 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² and switches to just x¹ 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…

\frac{\partial{|g|}}{\partial (-r)} = -\frac{\partial GM_E(r)/r^2}{\partial r} = 2GM_E(r)/r^3 - \frac{G}{r^2}\frac{\partial M_E(r)}{\partial r}

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, M_E(r)=4/3\pi r^3 \rho so we have

\frac{\partial{|g|}}{\partial (-r)} = \frac{8\pi G}{3} - 4\pi G = -\frac{4 \pi G}{3}

So 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…

M_E(r) = C r^2

and we have

\frac{\partial{|g|}}{\partial (-r)} = \frac{2 G C}{r} - \frac{2 G C}{r} = 0

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 2GM/r^3 - 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…

• 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…

​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!!!

If we modify simple BH accretion about a fixed object by applying an acceleration a 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 \frac{1}{2}at^2. If the acceleration is not perpendicular to the flow but at some angle \alpha 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₀ - we have the secondary accelerating towards the primary at r_s \Omega^2 and the projection relative to the incoming flow velocity as just r_s \Omega^2 \frac{v_s}{\sqrt{v_s^2+v_w^2}}. Given that acceleration of captured material occurs over a time scale \frac{R_{B}}{\sqrt{v_s^2+v_w^2}} we get that the offset should go as \frac{1}{2} \frac{v_s^3R_B^2}{r_s (v_s^2+v_w^2)^{3/2}} 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.

• 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)² level 2 steps… So for the entire simulation there will be X^D+1 (1+F R^D+1 + (F R^D+1)² + … + (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	X^D \propto N	fixed #cells per processor
Strong Scaling	X^D \propto 1	fixed #cells
Hybrid Scaling	X^{D+1} \propto N	fixed wall-time
Hybrid-AMR Scaling	X^{D+1} \left ( \frac{1-\epsilon^L}{1-\epsilon} \right ) \propto N	fixed wall-time (\epsilon=f(2^{D+1}) where f is the filling fraction and L 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⁺¹⁵3)=13% efficiency - without taking into account any latency or communication bottlenecks…. This is just do to having small grids of size 10³ - combined with extended ghost zones… If the grids are 20³ 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…

• Working on scaling tests #202
• Started gravoturb sims on Kraken #168
• Colliding flow restarts on Kraken #203
• Closed tickets #191 (threading), #195 (gfortran), #196 (cameras), & #199 (visit lower bounds)
• CollidingFlows movies

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³ mesh that I am remapping to a the middle half of a level 3 region of a simulation with a 1024³ 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

See blog:johannjc03062012b

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²=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² resolution - be only 6 pc across.. We could lower the density, or raise the temperature …