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…
To convert from to (Gayley 1995),
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.
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.
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.
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.
When using and , 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 to 0.1 from the typical 0.01, I had automatically calculated a new dinflo of .
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 /
When you impose a line strength cutoff of something like , 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 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)
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)
Mmax = Mmax * fdfac * (1.+deltams)
dstar = Mmax / (4 * np.pi * Rstar * Rstar * a)
dinflo = dfac * dstar
In the previous post I calculated wind statistics like and 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.