Using the athena code for non-ideal MHD simulations

[rahulgaur104]

Hi everyone,

I am a fifth-year graduate student at the University of Maryland, College Park. I work on fusion plasma turbulence. I want to simulate a plasma equilibrium for a rapidly(sound Mach number ~ 5) rotating plasma in cylindrical geometry using the non-ideal(resistive and viscous) MHD equations. Is it possible to solve non-ideal MHD equations using this repo of athena? If so, could you point me towards an example or tell me how to turn on resistivity and viscosity? All the problem generator files I have looked at either work with ideal MHD or use resistivity(resist.cpp) only during the initial setup.

Best regards,
Rahul Gaur

[c-white]

Non-ideal physics is described under Diffusion Processes. In order to turn on constant momentum, thermal, ohmic, and/or ambipolar diffusivities (units of code length squared per code time), just set the variables nu_iso, kappa_iso, eta_ohm, and/or eta_ad in the <problem> section of the input file. There are no configuration changes, nor does the problem generator need to be modified. Note that there are suggestive variables nu_aniso, kappa_aniso, and eta_hall that the code will read and process, but I do not believe anything will be done with them.

Non-constant diffusivities should also be possible, but that requires writing and enrolling certain functions in the problem generator.

I don't know if there are any test problems in the code for non-ideal physics + curvilinear coordinates. These should work together, but I recommend running a simple problem to make sure the code is doing what you think it's doing.

By default, the code will take explicit steps with diffusive terms. These might become much smaller than the sound crossing times, slowing down the calculation. If that is an issue, you should look into turning on super time stepping.

[rahulgaur104]

Thanks so much, Chris! I'll try to run a simple test case and share my input file with you by the end of this week.

[rahulgaur104]

Hi Chris and Kyle(@c-white, @felker),

Here are all the details of the simple test problem that I am trying to solve:

Problem definition

I am trying to solve an axisymmetric, axially symmetric MHD problem in cylindrical coordinates

I have turned on resistivity and viscosity. The system is isothermal. The initial conditions are as follows:

where a=0.15 and b=0.5 as you can see the attached problem file. Assuming radial force balance during the initialization, we must satisfy

I have assumed m = 1. Integrating the radial force balance equation, we obtain the initial density

which has been calculated exactly and is given in the attached pgen file. I have assumed n_0 = 1. The user-defined boundary conditions are as follows:

and periodic boundary conditions in theta. z_0=0.5 is the position of the Z boundary.

Outputs

Upon configuring the custom problem, compiling the code and running with the attached input file I get the following results.

• magnitude of the magnetic field becomes NaN
  

• the azimuthal flow speed vel2 becomes NaN
  

• density goes to 0
  

Related queries

The density should satisfy Neumann BC according to the problem file but it's 0 at r = a, b. Similarly, I want B_z =B_0 everywhere, even at the boundaries. I don't know how to do that. Finally, there could be a problem with my understanding of the normalizations. It would be helpful if you could explain how normalizations work in Athena++.

I would be immensely grateful for any help/suggestion. I have attached the problem generator and input file.
files.zip

[felker]

The hydro and MHD scaling constraints you listed seem correct to me, since we indeed calculate the Alfven speed in dimensionless internal code units as v_a = |B|/sqrt(dens).

I dont have any experience converting the output to SI units, however.

[c-white]

I'm not sure I understand all the notation above, so I'll try to explain normalizations from scratch.

The first thing to remember is that Athena++ is meant for a broader range of astrophysics problems than most other codes (e.g. Gandalf, Enzo, Arepo), so there is no natural physical scale to default to. You might be simulating a cubic millimeter of plasma or a galaxy.

For hydrodynamics, to convert between physical units and code units, you need three independent scales. It's probably easiest to choose length [L], time [t], and density [ρ]. So for example, if [L] = 1 km, then a code simulation running from x = 0 to x = 10 covers a length of 10 km. In any case we have the following dictionary:

Quantity	Code value	Physical value
Position	x	x [L]
Time	t	t [t]
Density	ρ	ρ [ρ]
Velocity	v	v [L] [t]-1
Pressure	p	p [ρ] [L]2 [t]-2
Energy density	u	u [ρ] [L]2 [t]-2
Momentum density	m	m [ρ] [L] [t]-1
Mass	M	M [ρ] [L]3
Diffusivity	η	η [L]2 [t]-1

There is no temperature or entropy, because these quantities do not enter into the equations of fluid dynamics. (If you have thermal diffusivity or radiation, then you need to know that the code treats p/ρ as "temperature." That is, k_B / (μ m_p) = 1 in code units.)

For the magnetic field, Athena++ is based on Lorentz-Heaviside units (rather than CGS, which makes the MHD equations slightly more complicated, or SI, which makes a mess of E&M entirely). For hydro, this distinction doesn't matter, since the above relations are true in all decent unit systems. To convert magnetic fields back to CGS, the dictionary includes:

Quantity	Code value	Physical value in CGS
Magnetic field squared	B2	4 π B2 [ρ] [L]2 [t]-2
Magnetic field	B	(4 π)1/2 B [ρ]1/2 [L] [t]-1

As a check, if [L] = 1 m, [t] = 1 s, and [ρ] = 1 kg/m³, then 1 code unit of B corresponds to about 11.2 gauss.

Wavespeeds don't come into this, except accidentally. Sound waves in Newtonian physics move at c_s = (Γ p / ρ)^1/2, both in LH/CGS/SI and in Athena++, so I suppose if p = 1 and ρ = 1 in code units, and if Γ = 1 (i.e., you are using an isothermal EOS), then c_s = 1 in code units. Similarly, Alfven waves move at v_A = (B^2 / ρ)^1/2 in LH and Athena++ (but not CGS or SI), so for B = 1 and ρ = 1 in code units, v_A = 1 in code units. But the above dictionary is derived based on dimensional analysis with no unnecessary constants; it just so happens that the formulas for some speeds also have no extra constants.

[rahulgaur104]

Hi Chris,

Thanks for such a detailed reply!

Let me convert all the variables from the code output to their dimensional values for a given set of values and you can tell me if I understand it correctly. At a given region in the magnetic mirror, say I have the following values in code units

Simulation run time = 1
Length of the cylinder = 1
Density of the plasma = 1
Pressure = 10
Magnetic field = 1
Sound speed(gamma = 1) = sqrt(10)
Alfven speed = 1

Since, we can choose [L], [t] and [ρ] freely, I choose [L] = 1m, [t] = 1s, [ρ] = 0.001 Kg/m^3. This would mean that the dimensional values are

Simulation run time = 1s
Length of the cylinder = 1m
Density of the plasma = 0.001 Kg/m^3
Pressure = 0.01 Pa
Magnetic field = 0.35 G = 3.5E-5 T
Sound speed(gamma = 1) = sqrt(10) m/s= 3.16 m/s
Afven speed = 3.5E-5/sqrt(4pi rho) = 3.1E-4 m/s

Is this conversion correct?

Thanks!

[c-white]

Everything is right except the Alfven speed. All velocities are converted with the same factor, here [L]/[t] = 1 m/s, so 1 becomes 1 m/s. Alternatively, calculate it in CGS: 0.35 gauss / (4 * pi * 1e-6 g/cc) = 100 cm/s. You can calculate it in SI using 3.5e-5 T and 0.001 kg/m^3, but the formula changes and probably has ε₀ somewhere in it.

[rahulgaur104]

Got it. Thanks, Chris!