Special relativity enabled versus disabled

[kurtsansom]

I implemented a new problem generator targeting solar weather farther than 0.1 AU from the sun.
Working with collaborators they suggested enabled special relativity. When I run the case with it disabled I appear to
get results that appear to make sense, but with special relativity enabled the results make it seem like there is a bug (only the ghost nodes appear correct) and things don't appear to be propagating. The collaborator appears to be getting results that do make sense.

So I am mainly confused on what enabling special relativity does, and what might be some gotchas that I am missing?

Caveat, I am not formally trained in astrophysics, so if someone can break it down "barney style" for me or point me to resources that would be greatly appreciated. I looked at wiki page but I don't know how to fit that into that larger scheme of things.

[c-white]

Simply switching on SR is something many people request of their hydro codes, but, as you've found out, it's not always that straightforward. In particular, there is no practical case where you can turn on SR without making other changes and still get a physically meaningful result.

In Athena++, enabling SR (configuring with -s) has the following primary effects:

• The relationship between primitive and conserved variables changes dramatically. The primitives are rest-frame, rest-mass density ρ; rest-frame pressure; and spatial components of 4-velocity. The conserved density is not ρ, but rather D = γ ρ, where γ is the Lorentz factor. Similarly, the momentum density has a Lorentz factor. The energy density is a very different beast, as it includes (and is sometimes dominated by) the rest-mass energy density. For sufficiently small velocities, the primitives reduce to the Newtonian primitives, modulo units (see next point). Note that in SR, magnetic fields contribute to total momentum, as well as total energy.
• Velocities are always in units of the speed of light. In the Newtonian version, you are free to interpret the physical units for 3 of {length, time, density, pressure, velocity} any way you want, with the others constrained by [velocity] = [length] / [time] and [pressure] = [density] * [velocity]^2. In SR, [velocity] = c, so you only have 2 degrees of freedom in interpreting the inputs and outputs to the code. In all cases you have ([B] / gauss)^2 = 4 π [pressure] / (dyne/cm^2).
• You are locked in to an adiabatic EoS. Isothermal is physically incompatible with SR, and while a general EoS could work, this hasn't (yet) been implemented
• In general, code robustness is reduced. There are simply more things that can go wrong with more complex equations. For example, numerical issues could cause a cell to have a superluminal velocity. Athena++ tries to trap as many issues as possible before they NaN the whole simulation, but with SR you should always check that the results look sensible. There are extra limits gamma_max, sigma_max, and beta_min in the <hydro> section of the input file to help with robustness.
• Accuracy could also be reduced if velocities are always very low. There are no hard limits here, but the wisdom is that for very sub-relativistic problems, using the full power of SR is not necessarily more accurate (you can imagine there are equations that care about γ - 1, which will become the difference of very close numbers).
• The code will probably be slower. The one exception is cases where the Newtonian code was being limited by superluminal speeds (which often happens when a large magnetic field seeps into a diffuse atmosphere near the density floor).
• Timesteps are now fixed (unless using AMR). Athena++ assumes that if you are using SR, something is traveling near the speed of light, so timesteps are set by light crossing times, not sound crossing times. Thus your simulation will never slow to a halt, and you can very accurately estimate runtimes.
• The internal algorithms (Riemann solvers, variable inversion) dramatically change in detail, but this may be of no concern to users.

The first two points above are probably where most issues arise. When setting up initial conditions, I recommend setting the primitive array phydro->w, keeping in mind that (1) velocities should be in units of c, and (2) velocities include the Lorentz factor. Thus if you want the gas to move at 3e9 cm/s in both the x- and y-directions, you would calculate vx = 0.1, vy = 0.1, vz = 0. Then gamma = (1 - vx^2 - vy^2 - vz^2)^(-1/2), so ux = gamma * vx, uy = gamma * uy, and uz = 0. These u's are the primitives. After you set the primitives (and magnetic fields if necessary), call peos->PrimitiveToConserved to set the conserved values phydro->u. The code only needs phydro->u, but this is easier to get by deriving it from phydro->w.

For the outputs, you have the reverse when interpreting velocities. The output values vel1, etc., are in fact ux, etc. Thus you need to compute gamma = (1 + ux^2 + uy^2 + uz^2)^(+1/2), vx = ux / gamma, etc., and note vx will be in units of c.

When setting up a problem and interpreting results, remember that if you interpret 1 code length as [length], then 1 code time is [time] = [length] / c.

[felker]

As an illustration of the second part of @c-white's answer, check out https://github.com/PrincetonUniversity/athena/blob/master/inputs/hydro_sr/athinput.mb_1 and the corresponding problem generator, https://github.com/PrincetonUniversity/athena/blob/master/src/pgen/gr_shock_tube.cpp.

You can see that the input file takes the left and right states of the shock tube as 3-velocities relative to the speed of light. E.g. the left state V_x =0.9c. Then, the problem setup translates that into the proper four-velocity vector:

athena/src/pgen/gr_shock_tube.cpp

Lines 165 to 173 in 5261e49

	// Construct 4-vectors
	Real ut = std::sqrt(1.0 / (1.0 - (SQR(vx)+SQR(vy)+SQR(vz))));
	Real ux = ut * vx;
	Real uy = ut * vy;
	Real uz = ut * vz;
	Real bt = bbx*ux + bby*uy + bbz*uz;
	Real bx = (bbx + bt * ux) / ut;
	Real by = (bby + bt * uy) / ut;
	Real bz = (bbz + bt * uz) / ut;

which are the actual primitive variables used by the code to compute fluxes:

athena/src/pgen/gr_shock_tube.cpp

Line 212 in 5261e49

phydro->w(IVX,k,j,i) = phydro->w1(IM1,k,j,i) = u1;