Modifying MHD equations

[zhouhz1992]

Hi All,

I'm new to the code. I'm wondering in which source files can I modify the MHD equations? More specifically I wish to add dynamo terms in GRMHD simulations.

Best regards,
Hongzhe

[c-white]

The short answer is that the code does not have any easy way to modify the magnetic field, by design. Athena is extremely careful to evolve the field correctly, without ever adding a magnetic divergence, and tampering with this can easily break the code.

If you truly want to add a source term for the induction equation, the place to do it is src/field/calculate_corner_e.cpp. The function in that file has two steps. First, it calculates a cell-centered electric field from the cell-centered magnetic field and other quantities. The precise formulas depend on Newtonian/SR/GR. Second, it combines the cell-centered electric fields with the face-centered electric fields returned by the Riemann solver to create edge-centered electric fields in a numerically stable way. This step looks the same with all physics. Afterward, src/field/ct.cpp uses the edge-centered electric fields to update the face-centered magnetic fields, again in a way that is the same with or without relativity.

You can modify the edge-centered electric fields however you want after the code sets them, and the resulting magnetic field evolution will be divergence-free. An alpha dynamo prescription would look like adding to $E^x$ something that scaled with local $yz$-velocity shear, as well as with nearby $B^x$, and had a cutoff (such as doing nothing when plasma $\beta$ is low).

All that said, there is another danger with relativity. Simply maintaining a divergence-free field is not sufficient to guarantee the code will work. Many of the non-ideal MHD terms one sees in Newtonian MHD are in fact parabolic, but relativity demands, both for physical and for numerical reasons, that the equations be strictly hyperbolic. If you want to put a subgrid dynamo into GRMHD, you have to make sure your prescription is truly hyperbolic. I'm not sure I've ever seen such a thing, and in general it is an open question how to correctly add non-ideal terms to relativistic fluid dynamics.

[zhouhz1992]

Hi Chris,
Many thanks for this detailed answer! This is very helpful.