astrobear - Blog
https://bluehound2.circ.rochester.edu/astrobear/blog
About blog postsen-USTrac 1.4.1Meeting updateericaMon, 04 Mar 2019 15:11:56 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica03042019
https://bluehound2.circ.rochester.edu/astrobear/blog/erica03042019<p>
Updates are posted on the accretion project page, here:
</p>
<p>
<a class="ext-link" href="https://astrobear.pas.rochester.edu/trac/wiki/u/erica/SteadyStateAccretionProblemPage"><span class="icon"></span>https://astrobear.pas.rochester.edu/trac/wiki/u/erica/SteadyStateAccretionProblemPage</a>
</p>
Localized Bondi dev updateericaThu, 31 Jan 2019 20:18:15 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica01312019
https://bluehound2.circ.rochester.edu/astrobear/blog/erica01312019<p>
Developed a solver for determining which branch of the local Bondi flow kernel edge corresponds to, and integrates the equations inward, and pastes the soln on the grid.
</p>
<p>
See the different possible soln classes in attached PDF.
</p>
<p>
To test the solver have a problem module that specifies an outer boundary condition of the flow (constant values at r≥1), which is internal to the grid.
</p>
<p>
For this boundary, the code figures out which branch to integrate the equations on and pastes the steady state soln within r<1.
</p>
<p>
I let the sim run for 20 frames using this set up (holding the boundary values of the problem fixed for r>1), using the Krumholz accretion module (see attached rho.gif for a movie, and the pdf for the initial conditions).
</p>
<p>
I need to finalize some things with the new sink accretion algorithm to utilize this framework. Things are mostly built with that so shouldnt be much longer. Need to come up with test problems.
</p>
meeting updateericaMon, 15 Oct 2018 14:05:37 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica10152018
https://bluehound2.circ.rochester.edu/astrobear/blog/erica10152018<p>
increased resolution of supernova core collapse problem… am getting funky boundary issues and poor HSE in center core.
</p>
<p>
possible fixes:
</p>
<p>
-integrate outward from center instead of inward from ambient for HSE routine
</p>
meeting notesericaMon, 08 Oct 2018 13:01:44 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica10082018
https://bluehound2.circ.rochester.edu/astrobear/blog/erica10082018<p>
Am trying to interpolate the following 1D radial profiles onto Astrobear's grid in 3D octant symmetry:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%203.17.27%20PM.png" style="padding:0; border:none"><img alt="progenitor profiles" crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%203.17.27%20PM.png" title="progenitor profiles" width="35%" /></a>
</p>
<p>
These profiles represent a ~5 solar mass progenitor for a core collapse supernova. The red vertical line is the radial cut-off I am taking for initialization onto the grid. Beyond this radius am switching the profiles to constant ambient values.
</p>
<p>
These profiles are assumed to be in ~HSE. Thus, I used the HSE module in AstroBEAR for self-gravitating gas to produce the following grid-interpolated profiles:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%206.35.41%20PM.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%206.35.41%20PM.png" width="25%" /></a>
</p>
<p>
Neutrino cooling is being mocked up in the code using the empirical formula:
</p>
<p>
<span class="trac-mathjax" style="display:none">L =[ 2e+18*T(MeV)^6*\rho (g/cm^3) + 1.9e+25* T(MeV)^9/\rho (g/cm^3)] (erg/s)</span>
</p>
<p>
For these initial conditions, this works out to be ~ <span class="trac-mathjax" style="display:none">L=10^{27}</span> erg/s for the peak density. In computational units, this is <span class="trac-mathjax" style="display:none">L\approx 4.5</span>.
</p>
<p>
The initial internal energy of the central zone is,
</p>
<p>
<span class="trac-mathjax" style="display:none">iE=\frac{P}{\rho(\gamma-1)} = 1.467</span>
</p>
<p>
So if the cooling rate were constant over <span class="trac-mathjax" style="display:none">dt=.003</span>, would expect to see the internal
energy drop to <span class="trac-mathjax" style="display:none">iE=1.4535</span>. Here is a plot of the internal energy over dt:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica10082018/Screen%20Shot%202018-10-08%20at%208.00.27%20AM.png" style="padding:0; border:none"><img alt="iE" crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica10082018/Screen%20Shot%202018-10-08%20at%208.00.27%20AM.png" title="iE" width="35%" /></a>
</p>
<p>
In other words, if the cooling time is <span class="trac-mathjax" style="display:none">1/L\approx.25</span>, then within this timestep, <span class="trac-mathjax" style="display:none">L^{-1}/dt\approx 83</span>.
</p>
<p>
The HSE module produces this:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%208.22.43%20PM.png" style="padding:0; border:none"><img alt="HSE" crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%208.22.43%20PM.png" title="HSE" width="35%" /></a>
</p>
<p>
Pretty close, except at origin…
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%209.11.40%20PM.png" style="padding:0; border:none"><img alt="grav vs. press, initial frame" crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%209.11.40%20PM.png" title="grav vs. press, initial frame" width="35%" /></a>
</p>
<p>
After freefall time (or 38 cooling times):
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%209.17.53%20PM.png" style="padding:0; border:none"><img alt="rho decrease in 1 ff time" crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%209.17.53%20PM.png" title="rho decrease in 1 ff time" width="35%" /></a>
</p>
<p>
The mesh shows strange behavior on the reflecting sides
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%209.01.14%20PM.png" style="padding:0; border:none"><img alt="boundaries awry" crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica10082018/Screen%20Shot%202018-10-07%20at%209.01.14%20PM.png" title="boundaries awry" width="35%" /></a>
</p>
<p>
And no sink particles being used, <span class="trac-mathjax" style="display:none">\gamma=1.667</span>
</p>
Bondi update - April 30ericaMon, 30 Apr 2018 14:24:50 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica04302018
https://bluehound2.circ.rochester.edu/astrobear/blog/erica04302018<p>
Two Bondi runs: gam=1.0001, gam=1.4
</p>
<p>
In each case, the Bondi radius (GM/cs_inf<sup>2</sup>) is resolved by the same number of zones (~160).
</p>
<p>
Gam=1 doesn't blow up, but gam=1.4 does. It seems this can be attributed to adiabatic pressure increase in the kernel (due to mass pile up) that sends out a shock wave upstream into the Bondi flow. Does the pressure increase in the kernel support this?
</p>
<p>
Movies are here:
</p>
<p>
gam=1.4 <br />
Density: <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/gam1pt4_rho.gif" title="Attachment 'gam1pt4_rho.gif' in Blog: Bondi update - April 30">gam1pt4_rho.gif</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/gam1pt4_rho.gif" title="Download"></a><br />
Pressure (ram vs. thermal): <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/gam1pt4_press.gif" title="Attachment 'gam1pt4_press.gif' in Blog: Bondi update - April 30">gam1pt4_press.gif</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/gam1pt4_press.gif" title="Download"></a><br />
Speeds (sound vs. vmag): <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/gam1pt4_speeds.gif" title="Attachment 'gam1pt4_speeds.gif' in Blog: Bondi update - April 30">gam1pt4_speeds.gif</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/gam1pt4_speeds.gif" title="Download"></a><br />
</p>
<p>
gam=1 <br />
Density: <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/gam1_rho.gif" title="Attachment 'gam1_rho.gif' in Blog: Bondi update - April 30">gam1_rho.gif</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/gam1_rho.gif" title="Download"></a><br />
Pressure (ram vs. thermal): <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/gam1_press.gif" title="Attachment 'gam1_press.gif' in Blog: Bondi update - April 30">gam1_press.gif</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/gam1_press.gif" title="Download"></a><br />
Speeds (sound vs. vmag): <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/gam1_speeds.gif" title="Attachment 'gam1_speeds.gif' in Blog: Bondi update - April 30">gam1_speeds.gif</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/gam1_speeds.gif" title="Download"></a><br />
</p>
<p>
A comparison of the density and pressure at frame 0 and 1 for the different gamma cases is here:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/rho_comparison.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/rho_comparison.png" width="35%" /></a>
</p>
<p>
And pressure is here:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/press_comparison.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/press_comparison.png" width="35%" /></a>
</p>
<p>
These are the same physical time. I think the dynamical time is set by the gas properties at infinity, in which case these should also be the same dynamical time…
</p>
<p>
Lastly, mdot converges to its theoretical value in the gam=1 case, but not gam=1.4 case:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/mdot_actual_gam1.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/mdot_actual_gam1.png" width="35%" /></a><br />
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/mdot_actual_gam1pt4.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/mdot_actual_gam1pt4.png" width="35%" /></a>
</p>
<p>
The Bondi and accretion algorithm mdots are calculated from different expressions. Classical Bondi mdot is:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/bondi_acc.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/bondi_acc.png" width="50%" /></a>
</p>
<p>
whereas the Krumholz prescription uses:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica04302018/krum_acc.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica04302018/krum_acc.png" width="50%" /></a>
</p>
Page on SN shocked clumpericaMon, 22 Jan 2018 15:59:48 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica01222018
https://bluehound2.circ.rochester.edu/astrobear/blog/erica01222018<p>
<a class="ext-link" href="https://astrobear.pas.rochester.edu/trac/wiki/u/erica/scratch10"><span class="icon"></span>https://astrobear.pas.rochester.edu/trac/wiki/u/erica/scratch10</a>
</p>
<p>
<sup></sup> Work in progress
</p>
Review: Conserving angular momentum over accretionericaThu, 14 Dec 2017 20:40:07 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica12142017
https://bluehound2.circ.rochester.edu/astrobear/blog/erica12142017<p>
Enforcing mass and linear momentum conservation means that angular momentum will generally <em>not</em> be conserved over the accretion step. We can see this by first noting that (by mass conservation) the sum of the mass to be accreted from the kernel (LHS) is equal to the mass of the sink particle after the accretion takes place (RHS; post-accretion quantities are denoted by a prime):
</p>
<div class="trac-mathjax" style="display:none">\Sigma m_i = M'_{sink}$</div><p>
Similarly, momentum conservation says that:
</p>
<div class="trac-mathjax" style="display:none">\Sigma m_i v_{i} = P'_{sink} = M'_{sink} v'_{sink}$</div><p>
Note, this last equation gives the velocity of the center of mass of the accreted material, which must equal the velocity of the sink particle, given momentum conservation.
</p>
<p>
From here, we might be tempted to write:
</p>
<p>
<span class="trac-mathjax" style="display:none">\Sigma m_i r_i \times v_i = L'_{sink}</span>
</p>
<p>
However, we already know that the mass of the sink and velocity of the sink are <span class="trac-mathjax" style="display:none">M'_{sink}</span> and <span class="trac-mathjax" style="display:none">V'_{sink}</span>, respectively. Thus we have to write:
</p>
<div class="trac-mathjax" style="display:none">\Sigma m_i r_i \times v_i = M'_{sink} R' \times v'_{sink}$</div><p>
and ask for what R' is this equation valid. Rearranging this equation (and calling the LHS <span class="trac-mathjax" style="display:none">L_{acc}</span>), we can write:
</p>
<div class="trac-mathjax" style="display:none">\frac{L_{acc}}{M'_{sink}}= R' \times v'_{sink}$</div><p>
In general, <span class="trac-mathjax" style="display:none">L_{acc}</span> and <span class="trac-mathjax" style="display:none">v'_{sink}</span> will not be perpendicular, thus there will be no solution for <span class="trac-mathjax" style="display:none">R'</span>.
</p>
<p>
What if we instead let <span class="trac-mathjax" style="display:none">R'</span> equal the center of mass of the accreted material? Then, the equation for angular momentum conservation becomes:
</p>
<p>
<span class="trac-mathjax" style="display:none">\Sigma m_i r_i \times v_i = \Sigma m_i r_i \times v_{sink}</span>
</p>
<p>
Equality then requires that <span class="trac-mathjax" style="display:none">v_i = v_{sink}</span>, which also will not generally be true.
</p>
<p>
Thus, the sink angular momentum cannot strictly be set equal to the <em>accreted</em> angular momentum, if the angular momentum is to be conserved across the accretion step. Instead, we need an additional vector which can absorb the difference in the angular momentum. This is the reason behind devising a 'spin' angular momentum vector for the sink particle. With the spin vector, the <em>total</em> angular momentum across the accretion step can be conserved.
</p>
<p>
<strong>Using a spin angular momentum vector to conserve angular momentum</strong>
</p>
<p>
As discussed above, the sink angular momentum following accretion (<span class="trac-mathjax" style="display:none">L'_{sink}</span>) is set by the accreted mass and linear momentum from the kernel. It is given by:
</p>
<p>
<span class="trac-mathjax" style="display:none">L'_{sink} = M'_{sink} R' \times v'_{sink} </span>
</p>
<p>
As shown above, this angular momentum vector will generally differ from the total angular momentum accreted from the kernel (<span class="trac-mathjax" style="display:none">L_{acc}</span>):
</p>
<p>
<span class="trac-mathjax" style="display:none">L_{acc} = \Sigma m_i r_i \times v_i \neq L'_{sink} </span>
</p>
<p>
Thus, we can't simply set the particle's spin axis equal to <span class="trac-mathjax" style="display:none">L_{acc}</span>. Instead, we want to look for some function <span class="trac-mathjax" style="display:none">s</span> such that the <em>total</em> angular momentum across accretion is conserved. We call this <span class="trac-mathjax" style="display:none">s</span> vector the particle's spin. The spin can be found by enforcing conservation of the total angular momentum across the accretion step:
</p>
<p>
<span class="trac-mathjax" style="display:none">L'_{sink} + S'_{sink} = L_{sink} + S_{sink} + L_{acc}</span>
</p>
<p>
This equation shows that the updated total angular momentum of the sink = the old total angular momentum + the accreted angular momentum. By rearranging this equation, we can solve for <span class="trac-mathjax" style="display:none">s'</span>, which absorbs any excess angular momentum over <span class="trac-mathjax" style="display:none">L_{acc}</span>:
</p>
<p>
<span class="trac-mathjax" style="display:none">S'_{sink}=S_{sink}+ L_{sink}-L'_{sink} +L_{acc}</span>
</p>
Bondi accretion module modsericaTue, 28 Nov 2017 17:48:46 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica11282017
https://bluehound2.circ.rochester.edu/astrobear/blog/erica11282017<p>
I'll be implementing two changes to the Bondi Accretion module. The first will be a flag which will force sink particles to only accrete the *radial* momentum of gas (see <a href="https://bluehound2.circ.rochester.edu/astrobear/blog/erica11152017">previous blog post</a> for motivation). The second is a check that the accretion volume does not contain any cells that violate the Jeans condition following Bondi accretion. If so, the excess gas will be removed.
</p>
<p>
<strong>Radial momentum accretion fix: description and equations</strong>
</p>
<p>
Here is a diagram defining the relevant variables. Note, in what follows primes (') will denote accreted quantities.
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica11282017/sink_diagram.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica11282017/sink_diagram.png" width="500" /></a>
</p>
<p>
The momentum vector of each cell in the accretion volume is decomposed into its radial and transverse components (relative to the location of the sink particle):
</p>
<div class="trac-mathjax" style="display:none">p=p_r + p_\perp $</div><p>
Following accretion, only the radial component of the momentum will be accreted:
</p>
<div class="trac-mathjax" style="display:none">p'=p_r' + p_\perp $</div><p>
where
</p>
<div class="trac-mathjax" style="display:none">p_r' = f_m \frac {p_i r_i}{r} \hat{r_i} $</div><p>
and
</p>
<div class="trac-mathjax" style="display:none">f_m = \frac {q(1) drho}{q(1)} $</div><p>
(i.e. fm = the fraction of mass to be left in the cell following accretion; the radial momentum will be accreted proportionally to mass).
</p>
<p>
Subtracting the first two equations gives us an equation for the update of the momentum vector:
</p>
<div class="trac-mathjax" style="display:none"> p' = (p_r' - p_r) + p $</div><p>
which in cartesian coordinates is:
</p>
<div class="trac-mathjax" style="display:none"> p' = p_i r_i (\frac{f_m-1}{r}) \hat{r_i} + p_i \hat{e_i} $</div><p>
This last equation gives the accreted momentum in each cell within the accretion volume, for the case of radial momentum accretion.
</p>
'The Angular Momentum Problem' of Sink ParticlesericaWed, 15 Nov 2017 20:56:59 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica11152017
https://bluehound2.circ.rochester.edu/astrobear/blog/erica11152017<p>
So I think I digested the Krumholz accretion paper sufficiently well and will attempt to distill its key points regarding angular momentum accretion onto a sub-grid sink particle, as well as how this may ultimately connect to outflows and jets launched from sink particles.
</p>
<p>
The Krumholz accretion algorithm is fairly sophisticated in that it calculates a modified Bondi accretion rate for rotating gas within the accretion volume, and then only removes gas from each grid point that is bound to the sink particle (<span class="trac-mathjax" style="display:none">drho = Mdot ~dt</span>).
</p>
<p>
Currently, in AstroBEAR, the Bondi accretion routine also removes linear momentum and specific internal energy from each grid point based on the fraction of mass this drho constitutes. In doing so, <em>each component</em> of the linear momentum is reduced by the same fraction. This reduces the magnitude of the momentum vector in the ith cell, but not its orientation. In other words, this results in the loss of angular momentum within that zone, and thus over the kernel:
</p>
<div class="trac-mathjax" style="display:none">L=\Sigma L_i = \Sigma \vec{p}_i \times \vec{r_i}$</div><p>
where <span class="trac-mathjax" style="display:none">L</span> is the total angular momentum over the kernel, and the subscripts denote the given quantity in the ith cell. Note the radial vector points along a line connecting the cell center to the sink particle.
</p>
<p>
This loss in angular momentum is stored in the sink particle, which enables the accretion algorithm to effectively conserve mass, linear momentum, and angular momentum between the grid and the sink particle over an accretion step.
</p>
<p>
Now, onto the angular momentum problem of this routine…
</p>
<p>
Krumholz makes the argument (and after some thought I think I agree with him) that given the sink particle represents an object much smaller than the grid scale (i.e. think single star compared to protostellar core) the angular momentum of the sink particle is negligible to that of the gas on the scale of the accretion volume. For a rough idea, let's assume the accretion volume represents a protostellar core, than its angular momentum is approximately:
</p>
<p>
<span class="trac-mathjax" style="display:none">L_{core}=(100 Solar Masses)(1 pc)(km/s)</span>
</p>
<p>
whereas the angular momentum of the sun is roughly:
</p>
<p>
<span class="trac-mathjax" style="display:none">L_{sun}=(1 Solar Masses)(10^{-8} pc)(km/s)</span>
</p>
<p>
This argument (known as the 'angular momentum problem of star formation,' c.f. this <a class="ext-link" href="http://www.astro.yale.edu/larson/papers/AngMom09.pdf"><span class="icon"></span>review</a>) implies single stars simply cannot acquire the enormous angular momentum of their host cores and clouds: else they would rotate faster than their break up speed.
</p>
<p>
Thus, Krumholz decomposes the linear momentum vectors of each cell into a parallel component and a transverse components (wrt to the line connecting the sink particle and cell center) and <em>only accretes the radial component of the momentum</em>. That is, since the angular momentum of the object the sink particle represents is negligible compared to that of the accretion volume, it should be ignored. By only accreting the radial component of the momentum vector the angular momentum of the gas is preserved (as well as the radial velocity).
</p>
<p>
<strong> Connection to outflows </strong>
</p>
<p>
If one does not accrete angular momentum from the grid, then one cannot inject angular momentum into the grid from sink particles upon the jet launching step. Thus, if we choose to adopt this way of looking at things, our injected protostellar jets will not be rotating on their own (but perhaps will spin up due to the bulk motion of gas within the accretion volume).
</p>
<p>
This contrasts with the other school of thought, namely Federrath's sink particle algorithms (which are currently implemented in the code). This method DOES accrete angular momentum from the grid (as described above for the Bondi accretion routine), and thus under those conditions, injecting a rotating jet back into the grid should (hypothetically) be self-consistent and conservative.
</p>
Bondi Accretion Routine Review ericaWed, 08 Nov 2017 23:38:07 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica11082017
https://bluehound2.circ.rochester.edu/astrobear/blog/erica11082017<p>
The Bondi accretion routine begins by calculating the Bondi radius of the sink particle's host cell using the following equation:
</p>
<p>
RB = GM/(cs<sup>2</sup>+vs<sup>2</sup>)
</p>
<p>
This equation gives the sonic surface of infalling gas onto a moving point particle.
</p>
<p>
Next the routine defines a region around the sink particle (on the finest AMR level) from which gas will be accreted. Within this region, the accuracy of the hydrodynamical solution is lost and so we want to minimize its size. However, we want to accurately-enough resolve inflowing material into this accretion region. To balance these two effects we choose a radius = 4dx. This is called r_acc in the code.
</p>
<p>
The amount of gas removed from the cells within r_acc is given by the accretion rate * the timestep. This accretion rate should be carefully chosen when the sonic radius lies within r_acc. This is because the accretion rate sets the resultant pressure in the kernel, and so if the accretion rate is too high (the resultant pressure too low) the hydrodynamical solution in cells surrounding the kernel could be degraded in the case of subsonic inflow.
</p>
<p>
In the other limit: R_BH > R_acc, then too high of an accretion rate shouldn't matter so much. The only issue I can foresee is if the accretion rate is too low, in which case gas would pile up within the kernel and eventually artificially fragment since new sink particles are not created within the accretion region of already existing particles. For this reason it might make most sense to allow sinks to form within the accretion regions of already existing sinks (if R_BH> R_acc) and then immediately merge (within a time-step).
</p>
<p>
Next the accretion rate is calculated using a semi-analytical expression found by Ruffert & Arnett (1994). See Krumholz 2004, eqn. 11. To calculate rho_inf in that equation, we calculate an average rho within the kernel using assigned weights to each cell within the kernel. These weights depend on the ratio R_BH : r_acc and effectively allow cells with higher inflow speeds to contribute more mass to the accreting sink particle. Correspondingly, these weights are then used to calculate the mass removed by any cell within r_acc during the accretion step.
</p>
<p>
Now this is all well and good if the gas is not rotating within the kernel. However, If it is rotating then the mass should not be accreted spherically symmetrically from all cells within the kernel. To calculate a correction due to rotation, each cell within the kernel is broken up into 8<sup>3</sup> sub-cells each with 1/8<sup>3</sup> the mass, momentum, and energy of the host cell. The minimum distance between the sink particle and these sub-cells is then calculated. If this distance is not less than dx/4, then that sub cell does not contribute to the accreted mass (it has too high an angular momentum to ever fall onto the surface of the sink). Additionally, if any of the cells within the accretion region is unbound, it does not contribute mass to the sink particle. Note, there is a maximum amount of mass any cell can contribute to the sink particle of 25%.
</p>
<p>
Now with regard to momentum changes to cells within the accretion radius, I was surprised to learn that the Krumholz's model for Bondi accretion reduces the radial momentum only (along a line between cell center and sink particle), in order to preserve the *radial velocity and angular momentum of the gas*. This is because the sink particle is envisioned as being so small in comparison to the grid size that its accreted angular momentum is negligible.
</p>
Meeting update for April 18th, 2017ericaTue, 18 Apr 2017 17:54:17 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica04182017
https://bluehound2.circ.rochester.edu/astrobear/blog/erica04182017<p>
<strong>Feedback</strong>
</p>
<p>
Almost ready to start pulling out my feedback development branch and finalize testing/development for that. At this point we are at a crossroads — where we can either
</p>
<p>
a) improve the algorithm to simultaneously inject the target amount of linear and angular momentum into the outflow codes (currently the code has trouble with getting the right amount of angular momentum into the cones),
</p>
<p>
b) ignore angular momentum injection altogether (would limit accurate simulation of outflows on small scales, I suppose),
</p>
<p>
c) leave the code as is, which can introduce some potential wonky behavior of the outflows under certain circumstances (such as a counter rotating accretion disk/star system)
</p>
<p>
<strong>Reorientation Paper</strong>
</p>
<p>
The referee had a few comments for the paper before it could be accepted for publication. Namely, he thought there should be a convergence study in the paper (sims are currently running for this), suggested adding a infinite cooling case to the paper (which we are making a counter argument against, see attached draft of referee responses), and remaking some of the plots so that they have larger text (currently working on this now, should be finished by the end of the day).
</p>
<p>
The referee responses are due no later than May 9th, but would like to get them in by the *end of the week*. Since the report is so short, thought I would definitely run it by Jonathan at the very least before submitting.
</p>
<p>
<strong>Thesis Defense</strong>
</p>
<p>
I still need to finalize the last person on my committee and to check with Segev again about the scheduled defense date of August 4th.
</p>
<p>
In order to defend my thesis by August 4th, I need to have a copy of my finished thesis to the committee by *June 26* (almost two months away). How long do these things typically take to write? If it takes 3 weeks, then that leaves me a little over a month to get the feedback sims up and running for the next paper (before I have to hunker down and write my thesis).
</p>
<p>
<strong>Thermal Instability Presentation</strong>
</p>
<p>
I recently read Koyama & Inutsuka 2004 on the "Field Condition" (convergence condition for the thermal instability), and really enjoyed it. Put together a short presentation on this paper to share with you all (see link here: <a class="ext-link" href="https://docs.google.com/presentation/d/1uaQ5ukZL6sNtdwxfubdH_EHZIQfWY8EESwC1P0QtV9I/pub?start=false&loop=false&delayms=60000"><span class="icon"></span>https://docs.google.com/presentation/d/1uaQ5ukZL6sNtdwxfubdH_EHZIQfWY8EESwC1P0QtV9I/pub?start=false&loop=false&delayms=60000</a>).
</p>
How cooling is called/used in the codeericaFri, 31 Mar 2017 19:44:57 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica03312017
https://bluehound2.circ.rochester.edu/astrobear/blog/erica03312017<p>
The code's main program is scrambler.f90, of course. In it are the initialization sequence calls, like 'physics init'. One of the things physics init does is call cooling init.
</p>
<p>
Cooling Init is in Cooling.f90. It is commented to say "initializes cooling table", but actually doesn't do anything… (It is a stub). Perhaps it is the relic from an older version of the code.
</p>
<p>
Now, the main meat of the cooling functionality goes like this.
</p>
<p>
In source_control.f90, there is a function called "SrcDerivs". This function's main goal is to populate the dqdt array with the code's source functions. Dqdt array contains geometrical source terms (like cylindrical corrections), energy source terms (like cooling), etc. At the end of this function, we have, "SrcDerivs=dqdt" (arrays).
</p>
<p>
If cooling is turned on, this function SrcDerivs calls the function "Cooling" (located in the cooling module: cooling.f90). "Cooling" returns a dqdt for SrcDerivs, specific for the type of cooling that is switched on.
</p>
<p>
For instance, take IICooling. "Cooling" selects the type of cooling that is present with a case statement:
</p>
<p>
If Case(IICool), Call II_Cooling
</p>
<p>
Now, IICooling is located in !IICooling.f90 and it calculates a dqdt(iE) based on a heating/cooling rate which themselves are functions of the II Cooling Parameters (stored in the array !IICoolPars). Once dqdt(iE) is calculated for this cooling choice, it is stored in the iE's element of the array SrcDerivs.
</p>
<p>
SrcDerivs array is used to update the q array in the subroutine, "RKorder2" back in source_control.f90. This is something like: q(i)=SrcDerivs(i)*dt + q(i).
</p>
<p>
The procedure would be similar (i.e. SrcDeriv and q would be updated in a similar manner) if another cooling process was selected instead of IIcooling.
</p>
<p>
Thus, it would seem, the only things necessary for turning cooling on in the code are the following two steps:
</p>
<ol><li>Select the type of cooling in physics.data, e.g. iCooling=3 (for II cool) (as long as iCooling/=0, cooling should be turned on in the code)
</li><li>Define any cooling params to be used by cooling module by adding these to both the problemdata namelist in problem.f90 (e.g. IICoolPars() array for the case of II Cooling), as well as define their corresponding values in problem.data.
</li></ol><p>
You should now be ready to cool!
</p>
coolingNov. 2nd, 2016 Meeting UpdateericaWed, 02 Nov 2016 19:20:23 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica11022016
https://bluehound2.circ.rochester.edu/astrobear/blog/erica11022016<p>
Need letter to YCAA today
</p>
<p>
Questions on Hubble Proposal Draft
</p>
<p>
Paper list on zoom-in: <a class="ext-link" href="http://adsabs.harvard.edu/cgi-bin/nph-abs_connect?library&libname=Zoom-in&libid=56b0f4173b"><span class="icon"></span>http://adsabs.harvard.edu/cgi-bin/nph-abs_connect?library&libname=Zoom-in&libid=56b0f4173b</a>
</p>
<p>
Have an interview next week at Univ. of Hertfordshire. This is with:
Jim Dale, Martin Krause, Martin Hardcastle (as Director of the Centre), and Janet Drew (panel chair, and PI of the grant)
</p>
<p>
Dr Strange this weekend?
</p>
JOBS!ericaWed, 19 Oct 2016 19:56:09 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica10192016
https://bluehound2.circ.rochester.edu/astrobear/blog/erica10192016<p>
Stuff is CRAZY right now — I'm in the thick of job applying.
</p>
<p>
In terms of science, the referee report seems good for the reorientation paper. Will be working on that in the coming weeks.
</p>
meeting updateericaWed, 21 Sep 2016 20:21:43 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica09212016
https://bluehound2.circ.rochester.edu/astrobear/blog/erica09212016<p>
Rotating, supercritical Bonnor Ebert sphere forms a disk, protostar, and outflow.
</p>
<p>
Here is the density right before the outflow leaves the sphere:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09212016/Feedback_rotBE.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09212016/Feedback_rotBE.png" width="35%" /></a>
</p>
<p>
right after:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09212016/Feedback_rotBE2.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09212016/Feedback_rotBE2.png" width="35%" /></a>
</p>
<p>
Looks like a RT instability along jet axis.
</p>
<p>
Radial velocity (ignore the legend label — it's wrong) shows that the kernel is not directed outward, but beyond that, it is… Not sure why yet. Infall is happening in the disk though, which is as expected:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09212016/Vrad10233.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09212016/Vrad10233.png" width="35%" /></a>
</p>
<p>
The temperature shows the outflow is decades larger than the core. This may need to be adjusted:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09212016/temp0000.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09212016/temp0000.png" width="35%" /></a>
</p>
<p>
A movie of density can be found here:
</p>
<p>
<a class="ext-link" href="https://www.youtube.com/watch?v=rpUXAqMHjCs"><span class="icon"></span>https://www.youtube.com/watch?v=rpUXAqMHjCs</a>
</p>
Some results from outflow testing.ericaTue, 13 Sep 2016 17:22:16 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica09132016
https://bluehound2.circ.rochester.edu/astrobear/blog/erica09132016<p>
The following are plots of the outflow kernel, that is, the profiles of injected density and velocity of the two-component sub-grid model:
</p>
<p>
In fixed grid we see a very nice outflow. Here is the density (rho_outflow=rho-rho_amb):
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09132016/fixedgrid.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09132016/fixedgrid.png" width="35%" /></a>
</p>
<p>
But in AMR, the outflow looks like:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09132016/16zones.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09132016/16zones.png" width="35%" /></a>
</p>
<p>
This particle has a buffer around it of finest cells that is r=16 zones=kernel of outflow. We still see edge effects however. I tried making the buffer larger, to r=32 zones, but this doesn't impact the mesh strangely:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09132016/32zones.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09132016/32zones.png" width="35%" /></a>
</p>
<p>
Looking at the outflow velocity (which is along z, the particle's spin axis), we have in fixed grid:
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09132016/fixedgrid_vz.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09132016/fixedgrid_vz.png" width="35%" /></a>
</p>
<p>
And in AMR with 'r_sink=32',
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09132016/32zones_vz.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09132016/32zones_vz.png" width="35%" /></a>
</p>
<p>
I also played with a collapsing BE sphere. Here is a <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica09132016/BEcollapse_rho.gif" title="Attachment 'BEcollapse_rho.gif' in Blog: Some results from outflow testing.">movie</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica09132016/BEcollapse_rho.gif" title="Download"></a> of a marginally stable BE sphere that is perturbed to collapse by a density enhancement. The outflow of the particle is not fixed, and so I think it is aligning to random noise in the accreted spin ,which is why it looks off-center by the end of the movie. We should probably initialized the BE sphere with some rotation.
</p>
Movies of high res outflowericaThu, 11 Aug 2016 18:17:18 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica08112016
https://bluehound2.circ.rochester.edu/astrobear/blog/erica08112016<p>
Attached to this post.
</p>
Tests of the feedback routineericaSat, 30 Jul 2016 17:32:43 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica07302016
https://bluehound2.circ.rochester.edu/astrobear/blog/erica07302016<p>
Using the Bondi routine with the particle's position fixed:
</p>
<table class="wiki">
<tr><td> Mass conserving? </td><td> yes
</td></tr><tr><td> Can it run in AMR? </td><td>
</td></tr><tr><td> Can it run on multiple processors? </td><td>
</td></tr></table>
<p>
Using the Bondi routine with the particle's position fixed:
</p>
<table class="wiki">
<tr><td> Mass conserving? </td><td> yes
</td></tr><tr><td> Momentum conserving? </td><td>
</td></tr><tr><td> Angular momentum conserving? </td><td>
</td></tr></table>
<p>
Using the uniform collapse module:
</p>
Changing box size changes outflow characteristicsericaWed, 27 Jul 2016 19:59:24 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica07272016
https://bluehound2.circ.rochester.edu/astrobear/blog/erica07272016<p>
<strong>Case 1. Box radius= bondi radius</strong>
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/rho_slice_niceoutflow0084.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/rho_slice_niceoutflow0084.png" width="35%" /></a>
</p>
<p>
<a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/rho_slice_niceoutflow.gif" title="Attachment 'rho_slice_niceoutflow.gif' in Blog: Changing box size changes outflow characteristics">movie</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/rho_slice_niceoutflow.gif" title="Download"></a>
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/contour_rho_niceoutflow0084.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/contour_rho_niceoutflow0084.png" width="35%" /></a>
</p>
<p>
<a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/contour_rho_biggerbox.gif" title="Attachment 'contour_rho_biggerbox.gif' in Blog: Changing box size changes outflow characteristics">movie</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/contour_rho_biggerbox.gif" title="Download"></a>
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/spin_niceoutflow0084.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/spin_niceoutflow0084.png" width="35%" /></a>
</p>
<p>
<a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/spin_biggerbox.gif" title="Attachment 'spin_biggerbox.gif' in Blog: Changing box size changes outflow characteristics">movie</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/spin_biggerbox.gif" title="Download"></a>
</p>
<p>
<strong>Case 2. Box radius= 2xbondi radius</strong>
</p>
<p>
Movies:
</p>
<table class="wiki">
<tr><td> slice of rho </td><td> <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/slice_rho_biggerbox.gif" title="Attachment 'slice_rho_biggerbox.gif' in Blog: Changing box size changes outflow characteristics">movie</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/slice_rho_biggerbox.gif" title="Download"></a>
</td></tr><tr><td> contour of rho </td><td> <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/contour_rho_biggerbox.gif" title="Attachment 'contour_rho_biggerbox.gif' in Blog: Changing box size changes outflow characteristics">movie</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/contour_rho_biggerbox.gif" title="Download"></a>
</td></tr><tr><td> spin </td><td> <a class="attachment" href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07272016/spin_biggerbox.gif" title="Attachment 'spin_biggerbox.gif' in Blog: Changing box size changes outflow characteristics">movie</a><a class="trac-rawlink" href="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07272016/spin_biggerbox.gif" title="Download"></a>
</td></tr></table>
Meeting updateericaWed, 20 Jul 2016 19:46:30 GMT
https://bluehound2.circ.rochester.edu/astrobear/blog/erica07202016
https://bluehound2.circ.rochester.edu/astrobear/blog/erica07202016<p>
Fixed a few bugs with the outflow feedback routine and test module. Outflow looks to be qualitatively reasonable. Next have a suite of tests to run to try and break the new routines.
</p>
<p>
For now, am starting with the bondi module which prescribes a bondi flow onto the origin given a central mass. The central sink particle is .4 solar masses. It is accreting via bondi accretion routines and is launching an outflow at .50 efficiency (that is, 50% of the accreted mass is launched in the outflow). The opening angle is 30 degrees.
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07202016/rho0043.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07202016/rho0043.png" width="35%" /></a>
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07202016/rho_contour0043.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07202016/rho_contour0043.png" width="35%" /></a>
</p>
<p>
<a href="https://bluehound2.circ.rochester.edu/astrobear/attachment/blog/erica07202016/vrad_contour0043.png" style="padding:0; border:none"><img crossorigin="anonymous" src="https://bluehound2.circ.rochester.edu/astrobear/raw-attachment/blog/erica07202016/vrad_contour0043.png" width="35%" /></a>
</p>