Is there a simple way to do an axisymmetric 2D (i.e. RZ) simulation?

[i-z0mbie-wxm]

Hello, I'm Jonny Dadras. I'm a computational chemical physicist by training with a love for PMI.
Previously I was at Zap and I am currently contracting with a new fusion startup. I am not an expert in MHD codes.

I am trying to run a simple toy model of a plasma column in RZ. I made a cpp module that should generate and initialize the problem. From reviewing the code I believe the standard r, phi, z ordering is followed. I have seen several examples of problems in the R Phi plane, I believe (this is consistent with a paper by Skinner).

If I try to define a mesh in an athinput.par such that:
nx1=100, nx2=1, nx3=100; etc
It errors out; stating that 2D problems in x1-x3 plane not supported

Does anyone have any advice on how to proceed?

Thanks for any help

[tomidakn]

We are sorry but we do not support 2D RZ coordinates with nx2=1. This is not fundamental but just due to some historical reason.

If you want to do so, I suggest you to use either 2D spherical polar coordinates (r,theta), or a thin slice of 3D cylindrical coordinates (i.e., nx2=NGHOST). Alternatively, you may implement it by yourself - just swap x2 and x3 components in the Cylindrical coordinate class (note but in this case the coordinate becomes left-handed; you can simply flip the sign of the phi direction as well to keep it right-handed).

[i-z0mbie-wxm]

Thank you so much for your help. I have another naive question.
I believe I have defined a column of plasma in this wedged cylindrical geometry. But, there does not seem to be any dynamics. Even with initializing some velocity in z & theta, when I plot the density everything seems static. Do I need to use EnrollUserExplicitSource to update the system at each timestep? I thought initializing and then letting Athena solve the MHD equations would be sufficient.
Thanks again for your time and any advice.

[tomidakn]

If you just want MHD calculation, you do not need to add any user-defined source term by yourself. However, I cannot tell what is your problem without seeing your problem generator and input file.

[i-z0mbie-wxm]

Hello, sorry for the delay. I would be grateful for some help with setting up this toy model of a pinch in approximate RZ geometry

pinch_cyl.cpp

athinput.pinch.txt

Here is the c++ module I've started, and the input file--changed to *.txt to share here. I guess I'm having two issues I haven't been able to resolve. One is that it seems like my plasma column is virtually static; even adding in random seeds to the density, initializing some velocity, and using the ExplicitSource function. Another issue is that my dt seems bigger than I expect; the program calculates something like ~10 us; hence why I have certain values in the input file. Any suggestions/improvements would be greatly appreciated; thanks in advance for your time.

[i-z0mbie-wxm]

PS. I tried to setup the problem in code units, hence my values for rho_code and pres_code; I wonder if this was a mistake?

[i-z0mbie-wxm]

PPS. I switched to SI units and now the plasma seems to evolve and dt makes more sense. Though I still don't observe typical instabilities. I must have misunderstood the code units...
I look forward to any further assistance and advice.

[i-z0mbie-wxm]

Hello, hope you're doing well. I was wondering if you had any time to go over my module?
I added some random seeds to encourage sausage and RT instabilities; but, I still do not see clear indications of these.
Again, thank you for your time and any advice.

[tomidakn]

I am sorry for my sparse response. I was completely overwhelmed with my teaching duties.

I do not have time to run the calculation on my side, so let me point out some possibilities.

1. The coordinate system is not an input file option but an configuration option. So make sure you configure the code with --coord=cylindrical and rebuild the whole code. The large time step may be caused by this. Also, I suggest you to use x2min = -x2max to improve the symmetry.

2. Using physical (SI / cgs) units may be causing the issue, because there are some implicit floors within the code. I recommend you to set units so that everything (the size of the domain, typical density, pressure, flow/sound/Alfven speeds, etc.) remains near O(1). Also, please note that Athena++ does not use Gauss units but "computational units" for the magnetic flux density B such that the magnetic pressure is B^2/2, not B^2/8\pi (in other words, the magnetic permeability \mu=1). Please refer to the method paper.

3. What are the plasma beta and Mach number you want to calculate? Our scheme is not robust when the flow is strongly dominated by magnetic fields and/or is extremely supersonic.

4. ApplyCurrentPulseInternal() is problematic. I guess you want to change B\theta depending on time using this source function. This is actually very tricky in Athena++ as we use constrained transport for MHD. If you change the face-centered field naively, it may violate the divergence constraint, and the code is very sensitive to divB errors. So we officially do not allow this (that's why our source function interface does not have access to FaceField). If you dare to do this, technically you can, but you have to be very careful and you need to know how the code works internally. In your case, you want to change all B\theta. Technically this does not violate the divB constraint if done correctly. While not recommended, you can directly update pfield->b as you did. However, as this is x2 face, the loop limits should be (is - ie, js - je+1, ks - ke), not (is - ie+1, js - je, ks - ke). And, you need to update the total energy cons(IEN,k,j,i) as you did in the initial condition.

5. You do not need nx2 = 16, just nx2 = 4 (or even 2) should be enough. Also, you do not need to split the domain in the x2 direction, so set <meshblock> nx2 equal to <mesh> nx2.

Please show me some figures; does the magnetic field evolve as you expected?

[i-z0mbie-wxm]

Thanks for getting back to me. I hope your classes are going well. I did some professoring before a long time ago; I understand how challenging that can be.

Yes, I made sure to configure with cylindrical. Here's the head of my config.log:

  Your Athena++ distribution has now been configured with the following options:
  Problem generator:            pinch_cyl
  Coordinate system:            cylindrical
  Equation of state:            adiabatic
  Riemann solver:               hlld
  Magnetic fields:              ON
  Number of scalars:            0
  Number of chemical species:   0
  Special relativity:           OFF
  General relativity:           OFF

Sorry, for clarity, I had just set rho_code = 1 and pres_code = 0.016 by hand. Switching to SI units allowed me to see dynamics. Thanks for pointing again that u0 = 1. I knew that; but I hadn't enforced it after my calculation of B. I have now.

Thanks for noting the loop limits and the grid numerics. I have updated the energy now, by adding these lines to the ApplyCurrentPulse function:

  // correction tomidakn
  // Update the total En (IEN); assuming eint & ke updated via mhd consistently
  idx = 0;
  for (int k = pmb->ks; k <= pmb->ke; ++k) {
    for (int j = pmb->js; j <= pmb->je; ++j) {
      for (int i = pmb->is; i <= pmb->ie; ++i, ++idx) {
        Real B2_new = 0.5 *(SQR(pfield->bcc(IB1, k, j, i)) +
                       SQR(pfield->bcc(IB2, k, j, i)) +
                       SQR(pfield->bcc(IB3, k, j, i)));
        // Update total energy; correction tomidakn
        pmb->phydro->u(IEN, k, j, i) += mu0 * (B2_new - B2_old[idx]);
      }
    }
  }
...  

where I previously compute B2_old. Since B is now in the code units, I make sure to convert back by multiplying by mu0. But, I know this is an energy density; does this mean I should multiply by the volume of the cell--which I don't know here or in ProblemGenerator? I think I have done most everything consistently now; it looks like I'm observing the dynamics I expect.

The simulations are very long now as I started getting dt ~ 1e-14... When something takes a long time it means you're on the right track, right? :)

[i-z0mbie-wxm]

Hello, hope you are doing well.
The dt kept collapsing. I realized with my grid size, dx should be ~390 um, following the CFL condition, this would indicated some velocity O(10^3) c, which is unphysical.
I believe the energy update should already be in the compute units of Athena; so when mu0 = 1, I got simulations that completed and observed some reasonable dynamics. Though even with seeding some sausage-like instabilities, I don't see any clear pronounced indications of such phenomena.
But, when you have some time to maybe go over my code to see if there are an glaring errors, I may have missed, I would be grateful. Thanks again for your time and any advice.

[tomidakn]

In the code unit, you do not need mu0. Please refer to other code examples provided in the src/pgen directory. If there is any doubt, please read the source code as everything is public. For example, read src/eos/adiabatic_mhd.cpp to see how physical quantities are related in the code units.

You should analyze your data what is causing the short time step. Usually such a short time step means something is broken, like magnetic fields exploded, density or pressure gets small or even negative (and floored to very small positive value, but the floor failed to save the situation), or something else. If you want my input, please provide more detailed information with some figures, not just your impression. Again, grid-based codes are not very robust when plasma beta is very low, and you might be exceeding the limit. In this case, you have to change the method.

I am happy to help, but please understand that I do not have time to check every user's code neither run your code on my side. It is already midnight here in Japan and I am still working at my office. I'm very sorry but please help yourself as much as possible.