Alternative energy conservation (or entropy consservation) for pressure positivity in low beta regions

[kyriakosTapinou]

Has anyone implemented an alternative energy conservation for the non-relativistic MHD solver for use when pressure positivity becomes an issue in low beta regions?

Maintaining a physical pressure while using a conservation of total energy can be problematic when the magnetic energy and/or kinetic energy are orders of magnitude greater than the internal energy.

I’m simulating a high-magnetic field region (solar active region with a flux tube emerging into the corona) with low densities. It seems the magnetic energy dominates and pressure suffers. This is a common problem but most approaches simply recommend ditching the total energy conservation for a conservation of internal energy or hydrodynamic energy only. I've also seen people implement an entropy conservation to try resolve the issue.

Do you have any ideas that don't include a different energy conservation or using a larger pfloor?

Thanks in advance

[kyriakosTapinou]

I think in general to do this there must be a rigid problem setting i.e.

• a specific flux solver/s (unless the change is made to all flux solvers)
• a new hydro class e.g. adiabatic_mhd_low_beta

then the obvious places to make the required changes are in

• the specific flux solver to remove the magnetic energy flux contribution
• conserved to primitive functions

I'm sure there are many more details which I'll post here as I go.

[tomidakn]

Sorry for my late response. My colleague, Shinsuke Takasao at Osaka University has implemented so-called dual energy formula, which solves the total energy and internal energy, and uses the internal energy formula as a floor when the total energy formula fails. I think we can merge it to the master at some point.

[kyriakosTapinou]

You're right on time. That's a good feature, props to Shinsuke. When do you think that would be merged in?

Alternatively, is Shinsuke's fork of Athena public so that I could grab a patch from it?

[tomidakn]

I asked Shinsuke to join this discussion.

He said he is able to provide his code but based on an older version now, or you need to wait for a while (he is moving to another institute this spring) if you want it on the latest version of the master.

[kyriakosTapinou]

I'd gladly take Shinsuke's old code and update it myself. Happy to make a pull request if it works well. I think this will be useful for many people - it'll open up more use cases for athena in general.

[apbailey]

this would be potentially useful for me as well, just to add to whatever demand there is for it

[msbc]

I just realized that my fork of Athena where I've integrated passive scalars into the EOS may also be useful here. I've used the pasive scalers to store the last retrieved cell temperature to use as an initial guess for the EOS root find.

[kyriakosTapinou]

Is this to avoid the current time step temperature, which may be wrong/unstable, from being used? I'm not familiar with the general equations of state and the root finding used in them. My personal use case is only working with ideal gas equation of state at the moment. Is the general equation of state still not activated for MHD (I read in the docs)?

[msbc]

The general EOS module is turned on by a configuration flag. I was using an EOS that is a function of temperature and density. Athena doesn't have/keep track of temperature, so a made a fork of the code where the pasive scalers also get passed to the EOS calls and the Riemann problem. I store temperature in the scalars every time the EOS is called, and I use that as an initial guess the next time I have to convert either the primitives or conservatives to density/temperature space. Also the EOS module handles prim cons conversations, regardless of which EOS you use.

[takasao]

Sorry for the late response. I am currently very busy because I am moving to another institute. I have quickly created the public github repository of athena++ with the implementation of dual energy formalism.
https://github.com/takasao/athena_dualen
Please check pgen/blast_dual_en.cpp and inputs/mhd/athinput.blast_dual_en for a demonstration. The PDF file Dual_energy_formalism.pdf generally describes the method. Please note that the current version of this code is based on an older version of athena++ master (Commit: 850eaa7).

[kyriakosTapinou]

Thank you, I'm looking through this now. I'll hopefully get it working in the latest version without any troubles!

[kyriakosTapinou]

@takasao I've implemented your dual-energy formalism in my work and it seems to be function as advertised.

Some regions of high velocity field divergence lead to errors of negative internal energies and eventually crashes when the RK4 or RK3 time integrators are used. The VL2 integrator seems to work well but I need to test more thoroughly. This is for a MHD problem.

@tomidakn I'm not sure if you want my re-implementation of Shinsuke's work or not, let me know. If so, I also have some data-driving boundary stuff (B-driving and E-driving) which may also be useful.