f90 radiative cooling implementation probably works

Taking a random snapshot of a simulation (lots of structure), I ran the f90 and f77 radiative cooling functions and measured the change in temperature that each one applied to each grid cell. Here, black is f77 and red is f90.

Here’s the cooling, plotted together.


Here’s the difference in cooling, plotted.


The differences are on the order of 5%. I’m not sure this is significant, or where it comes from…

Posted in exploration, figure | Leave a comment

work log on radiative cooling bug

9:51 AM. Radiative cooling implementation has no obvious bugs. I’ll test both functions on the same beta-law wind and see if the output is the same. However, first I’ll recompile and test — maybe I made some mistake incorporating it into the F90 code.
10:07 AM. It appears some code for the Abbot delta has been included in what I received from David. I’ve removed some elements of it, which will help isolate the issue when I’m debugging.
10:19 AM. Spent ten minutes cleaning out french press, and got a snack from the kitchen.
10:25 AM. Code has passed the most basic test — I put a write(*,*) line in the radiate.f90 file and compiled. The radiate function is being called, which is good. There are two possibilities remaining — either the Abbot delta code I removed was causing the issues, or the radiate.f90 function is broken. Hopefully it’s the former. I’m running a test simulation to see if this is indeed the case.
10:31 AM. Simulation time/real time ratio is about 100, so this particular simulation should take about an hour to run. I realized that I failed to back up my emails with Dylan Kee on Swatmail, so my communications with him on the issue of more advanced radiative cooling are forever lost (I lost access to Swatmail about a week ago). I hope I backed up some of the files he send me…
10:41 AM. Briefly distracted by internet. Opened up Authorea and saw “no modifications since 9 days ago”… BibleThump. I’ll fill out Section 3.1, Line Strength, which should discuss Qbar and our CAK line exponent.
11:14 AM. Wrote a subsection in the paper on the relation between Qbar and Z. This could probably merit some more discussion, but I think I hit the important points (Qbar ~ Z, finite Qcut required). I’m going to take a quick break.
11:27 AM. Simulation finished, with the same strange behavior as before. This means we’ll have to dig deeper, and the issue wasn’t the Abbot delta code.
11:37 AM. Ate some quiche, filled up mug.
11:40 AM. The first step is simple — I need to make the f77 and f90 code have the same initial condition. Then I need a conversion between Qbar and CAK kappa. Then I can compare the new and old implementations of radiative cooling.
11:42 AM. I wonder if Jon already has a script to convert vh1.init files between F77 and F90? Ah well… reinventing the wheel vs. sending bothersome emails.
12:06 PM. Initial condition file formats are identical, just need to copy-paste header. Nice! Now to look up the conversion between kappa and Qbar.
12:12 PM. Now that I’ve thought about it for a while, I don’t need to match up Qbar and CAK k. I just need to look at the frame 0 radiative cooling.
12:31 PM. Lunchtime!
3:16 PM. Took a long nap, now back at it. Step 1: make f90 code display radiative cooling information.
3:46 PM. f90 code now displays pressure grid prior to cooling and after, only for the first hydro step. Now to prepare a nice initial condition file.
3:50 PM. The vh1.fin from the f77 run I did earlier this week is a perfect test case. I’ve saved it and now I’ll use this as the vh1.init for both f77 and f90 codes. Next step is to make the f77 code do the same as the f90 does — print out the pressure variables before and after cooling. Also important: I should do my testing with radiate called BEFORE sweep, just because I don’t want sweep to change anything about the initial condition.
4:23 PM. Took a break and wrote a short python script to convert these logs into HTML tables for blog posts. The formatting looks pretty nice! I’ll cook dinner in a few minutes and I’m out of coffee, so I might call it a day.
Posted in bugfix | Leave a comment

To convert from \bar{Q} to k (Gayley 1995),

    \[k=\frac{1}{1-\alpha } \left(\frac{T_{eff}}{5.5 \times 10^{12}}\right)^{\alpha /2} \bar{Q}^{1-\alpha } .\]

Posted on by Zack (Zequn) Li | Leave a comment

weird stuff with radiative cooling

The current F90 implementation of radiative cooling appears to be buggy. Here’s the old behavior,

and here’s the new implementation,

I’m sure some careful debugging will fix this issue.

Posted in bugfix | Leave a comment

mean f_cl in radial ranges — prelim

Preliminary work on the clumping factor in certain radial ranges reveals that the clumping factor is really noisy for runs that last 1e6 seconds. The final result should maybe have 10 times the simulation time. I’m getting my wisdom teeth removed at some point this summer, so maybe that’s when I can run the simulations for a really long time.

Here are two radial ranges (1.15 – 1.5 in blue and 1.5 – 2 stellar radii in green) with mean clumping factor plotted vs. Qbar. It should probably be monotonic, suggesting our noise is pretty large here.


Looking at the simulation videos, it seems like the culprit is “sputtering” going at the base, causing long-time-scale variations.


Posted in figure | Leave a comment

writing is pretty difficult

I’m having a hard time writing, and despite our interesting result, the words are not coming easily. It seems to me that writing is a big part of a scientist’s job, so I really want to get better at this!

I’ve yet again spent the morning with full intent of producing more material, but with no real text being produced. Like any skill, I’m sure I can get better through some deliberate practice. I’ll outline my plans and hopefully a few lessons for myself to remember, in this post.

Outlining. I spent probably a week procrastinating at this stage. I should have eschewed the computer for just pen and paper, since drawing figures and doing other visual note taking probably would have helped me organize my thoughts. Instead, I spent hours going through papers on ADS to try and get some context for writing, which really didn’t help.

Discretizing. After setting up small chunks of the paper with mini-outlines, it’s now easier to know where to start (the hardest part of doing anything I think) whenever I sit down to write.

Introductions. I pooped out probably the worst introduction in the history of scientific writing last week, and I still am grappling with what to do about it. I still don’t really know how to address this.

Pending Results. Even though my simulations are still running, I can still write about the preliminaries and make amendments when time goes on. What’s the best way to approach this? Not sure.

Posted in problems | Leave a comment

I have a PC at home which is substantially faster than my laptop, so I’ve spent the day configuring it so I can log into it through ssh and run simulations.

Posted on by Zack (Zequn) Li | Leave a comment

Qbar = 2000, Qmax = 0.1

When using \overline{Q} = 2000 and Q_{max} = 0.1, we get some really weird behavior. This is reminiscent of the delta factor explorations from a year ago.

Here, I’m using the following Indat file. An important thing to note is that due to the increased Q_{max} to 0.1 from the typical 0.01, I had automatically calculated a new dinflo of 1.7 \times 10^{-11} \text{ g/cm}^3.

 rstrt   = 'no'
 regrid = 0
 prefix  = 'ldi'
 imax    = 2000
 jmax    = 1
 kmax    = 1
 xmin = 1.00
 xmax = 10.0
 ymin = 0.
 ymax = 1.0
 zmin = 0.
 zmax = 1.0
 dxfac = 1.002
 courant = 0.25
 dtmax = 2.3 
 t_stat = 5.e5 
 lum = 8.0e5
 mass = 50.0
 radius = 20.0
 tempw = 40000.
 isotherm = 1
 Qbar   = 2000.0
 Qmax  = 0.1
 ifrc = 4
 amp = 0.0
 phase = 0.0
 period = 4.e3 
 ncycend = 900000000
 ndump   = 900000000
 nprin   = 900000000
 nmovie  = 900000000
 endtime = 5.01e5
 tprin   = 1.e2
fvel = 0.2
dinflo  = 1.72437934367e-11
rcl = 1.3
alpha = 0.65
 tmovie  = 5.0e+05 /
Posted in problems | Leave a comment

Qmax dinflo and mass loss

When you impose a line strength cutoff of something like Q_{max} / \bar{Q} = 0.01, the mass loss rate should be affected. I wrote up a short Python function to compute the mass loss rate, and thus the required density lower boundary condition (“dinflo”).

You would call this function with arguments like “computeDinflo(2000,0.01)”.

import math
import numpy as np 
from scipy.optimize import fsolve

def computeDinflo(Qbar, Qmax):
    dfac = 5#
    k = 1.380658e-16 # (* eV / K *)
    T = 40000# (* K *)
    Rsun = 6.955e10# (* radius of sun *)
    Rstar = 20* Rsun# (* radius of \[Zeta] pup *)
    Lsun = 3.839e33# (* luminosity of sun *)
    Lstar = 8e5 * Lsun# (* luminosity of \[Zeta] pup *)
    Msun = 1.989e33# (* mass of sun *)
    Mstar = 50 * Msun# (* mass of \[Zeta] pup *)
    mH =  1.6737236e-24 # (* mass of H atom *)
    mu = 1.3 # (* mean mass in units of Subscript[m, H] *)
    # Qbar = 2000 # (* bounciness of wind *)
    Qmax =  Qmax * Qbar#
    c = 2.998e10# (* speed of light *)
    csq = c * c
    G =  6.67e-8# (* grav constant *)
    kappae = 0.34# (* electron opacity *)
    alpha = 0.65#(* CAK linelist power law *)

    gammae = kappae * Lstar / (4 * np.pi * G * Mstar * c)
    vesc = np.sqrt(  2. * G * Mstar * (1-gammae) / Rstar )
    Mcak = (Lstar / csq) * (alpha/(1-alpha)) * ( (Qbar * gammae) / (1-gammae) )**((1-alpha)/alpha)
    fdfac = 1. / (1+alpha)**(1./alpha)
    a = np.sqrt(k * T / (mu * mH))
    deltams = 4. *( ( np.sqrt(1.-alpha)) / alpha) * a / vesc

    Mfdfscak = Mcak * fdfac * (1 + deltams)

    def func(m):
        return (1 + Qmax*(m*gammae/( 1 - gammae)*(1 - alpha)/alpha))**(1 - alpha) - (1 + (Qmax/Qbar)**( 1 - alpha) *(m*gammae/( 1 - gammae)*(1 - alpha)/alpha)/(gammae/(1 - gammae)))

    Mmax = (Lstar / csq) * fsolve(func, Mcak * csq / Lstar)[0]
    Mmax = Mmax * fdfac * (1.+deltams)
    dstar = Mmax / (4 * np.pi * Rstar * Rstar * a)
    dinflo = dfac * dstar

    return dinflo
Posted in simulation | Leave a comment

using the wind statistics implemented by Jon

In the previous post I calculated wind statistics like \langle \rho^2 \rangle and \langle v \rangle in Python, outside of the simulation code. This meant I was taking snapshots of information every 500 seconds, which gave some scatter. Jon informed me that his implementation already outputs some of this information, so I’m now using these averages that use information from every timestep.


Posted in figure | Leave a comment