Follow the convention used in the paper to allow for td and tf
Merge branch 'master' into develop
Implement an extrapolation scheme
Properly compute FLOPS according to the new rhs()
Properly compute RWSZ according to the new rhs()
Fix bug: FLOP_RHS actually depends on n_nu
Turn off extrapolation by default
Accretion flow is sub-Keplerian; also we keep the magnetic fields
Set B=0 to void unrealistic synchrotron radiation
Change variable names for consistency
Turn N_RX into a constant c.n_rx and initialize it to c.n_r - 1
Compute c.n_rx dynamically from the rx runtime parameter
Use a bigger time step cap in large radius
Give up placing theta-grid in constant memory for simplicity
Use shorter names for local variables
Use single character names for local indices
No need to keep the parameters "Rd" and "Rf"

We are refitting everything in the paper.  The two parameters "Rd" and
"Rf" are outdated.  Use their inverses "td" and "tf" instead.
Use a portable way to initialize structs to zero
Set d.I[] = d.tau[] = 0 for s.r < c.r[0], instead of extrapolating
Find i and I so that c.r[i] <= s.r < c.r[I]; compute weight f