Uneven time interval between outputs.

[mishraabhijeet0]

I am trying to implement a non-uniform magnetic field (instead of the uniform constant field provided with the jet.cpp pgen) around the jet inlet in the ambient medium. Here is the snippet of the implementation.

if (MAGNETIC_FIELDS_ENABLED) {
    Real alpha=10;
    AthenaArray<Real> psi;
    int nx1 = block_size.nx1 + 2*NGHOST;
    int nx2 = block_size.nx2 + 2*NGHOST;
    int nx3 = block_size.nx3 + 2*NGHOST;
    psi.NewAthenaArray(nx3, nx2, nx1);
    
    for(int k=ks; k<=ke+1; ++k) {
      for(int j=js; j<=je+1; ++j) {
        for(int i=is; i<=ie+1; ++i) {
          if((SQR(pcoord->x2f(j)) + SQR(pcoord->x3f(k)))<(r_jet*r_jet + 3) &&
          (SQR(pcoord->x2f(j)) + SQR(pcoord->x3f(k)))>(r_jet*r_jet))
          {
            psi(k,j,i) = (SQR(pcoord->x2f(j))+SQR(pcoord->x3f(j)))*
                          exp(-SQR(pcoord->x2f(j))-SQR(pcoord->x3f(k))-SQR(pcoord->x1f(i)));
          }
          else{
            psi(k,j,i) = 0;
          }
          
        }
      }
    }
    for (int k=ks; k<=ke; ++k) {
      for (int j=js; j<=je; ++j) {
        for (int i=is; i<=ie+1; ++i) {
          
          
            pfield->b.x3f(k,j,i) = -(1/(SQR(pcoord->x2f(j))+SQR(pcoord->x3f(k))))
                                   *(pcoord->x3f(k)*((psi(k,j,i+1)-psi(k,j,i))/pcoord->dx1f(i))+
                                   (alpha*pcoord->x2f(j)*psi(k,j,i)));
          
        }
      }
    }
    for (int k=ks; k<=ke; ++k) {
      for (int j=js; j<=je+1; ++j) {
        for (int i=is; i<=ie; ++i) {
          
          
            pfield->b.x2f(k,j,i) = (1/(SQR(pcoord->x2f(j))+SQR(pcoord->x3f(k))))*
                                   ((-pcoord->x2f(j)*
                                   ((psi(k,j,i+1)-psi(k,j,i))/pcoord->dx1f(i)))+
                                   (alpha*pcoord->x3f(k)*psi(k,j,i)));
          
        }
      }
    }
    for (int k=ks; k<=ke+1; ++k) {
      for (int j=js; j<=je; ++j) {
        for (int i=is; i<=ie; ++i) {
          
          
            pfield->b.x1f(k,j,i) = (1/(SQR(pcoord->x2f(j))+SQR(pcoord->x3f(k))))*
                                   ((pcoord->x3f(k)*
                                   ((psi(k+1,j,i)-psi(k,j,i))/pcoord->dx3f(k)))+
                                   (pcoord->x2f(j)*
                                   ((psi(k,j+1,i)-psi(k,j,i))/pcoord->dx2f(j))));
          
        }
      }
    }
    if (NON_BAROTROPIC_EOS) {
      for (int k=ks; k<=ke; ++k) {
        for (int j=js; j<=je; ++j) {
          for (int i=is; i<=ie; ++i) {
            phydro->u(IEN,k,j,i) += 0.5*(SQR(bx_amb) + SQR(by_amb) + SQR(bz_amb));
          }
        }
      }
    }
  }

I have tested for the divergence condition and it vanishes. Also, my jet inlet is from the inner x-boundary so I have placed the magnetic field around the inlet coiled over a cylindrical surface. However, if I try to run the simulation, with following inputs,

<output2>
file_type  = hdf5      # Binary data dump
variable   = prim     # variables to be output
dt         = 1.0      # time increment between outputs

<time>
cfl_number = 0.4      # The Courant, Friedrichs, & Lewy (CFL) Number
nlim       = -1       # cycle limit
tlim       = 300.0    # time limit
integrator  = vl2       # time integration algorithm
xorder      = 2         # order of spatial reconstruction
ncycle_out  = 1         # interval for stdout summary info

I expect 300 output files but I get only 19 files while the sequence terminates on time limit. Please explain what is the fault in the simulation. Also, the magnetic field array in the output has a shape (3,1,128,128,128 ) in this case. If you could explain the meaning of these dimensions it would be very helpful for me to plot the streamlines of this vector field. Thank you.

[felker]

What is printed to stdout when you run this simulation? How many timesteps does it execute?

(3,1,128,128,128 ) = 3 magnetic field vector components, 1 meshblock, 128 grid cells in x, y, and z dimensions

[felker]

the reason you only get 19 outputs is that your timestep dt > 1.0 in the last 9 cycles, so there is no way to dump ~290 outputs roughly every second of simulation time

[mishraabhijeet0]

I understand the problem. But I am not sure how dt is being updated given that I am using the vl2 integrator. In other simulations that I have created so far, say, for dt=0.05 and tlim=2.0 , the simulation generated exactly 40 output files. But here this is not the case. Can you please elaborate on this?

[c-white]

The dt being printed is not the dt from the input file. The former is the stepsize chosen by the code to advance the simulation one step further in time; the latter is your requested output cadence. The actual time between outputs can only be the maximum of these two values.

You need to look at your output data. The fact that the timestep changed by over 5 orders of magnitude means the state of the simulation (in particular the fast magnetosonic speed) changed by that much. The simulation almost certainly has garbage values, possibly nans, and it looks like it started going wrong early on.

[mishraabhijeet0]

Thank you for the insight. I found out one reason for this to be the incorrect equation of state. The issue was still not solved, so I was trying to implement the jet.cpp file as provided with the athena++ distribution in 3D by changing nx3 in the input file from 1 to some higher value. This did not work as you can see below:

cycle=0 time=0.0000000000000000e+00 dt=3.7499999999996252e+00
cycle=1 time=3.7499999999996252e+00 dt=3.7499999999996252e+00
cycle=2 time=7.4999999999992504e+00 dt=3.7499999999996252e+00
cycle=3 time=1.1249999999998876e+01 dt=3.7499999999996252e+00
cycle=4 time=1.4999999999998501e+01 dt=3.7499999999996252e+00
cycle=5 time=1.8749999999998124e+01 dt=3.7499999999996252e+00
cycle=6 time=2.2499999999997748e+01 dt=3.7499999999996252e+00
cycle=7 time=2.6249999999997371e+01 dt=3.7499999999996252e+00
cycle=8 time=2.9999999999996994e+01 dt=3.7499999999996252e+00
cycle=9 time=3.3749999999996618e+01 dt=3.7499999999996252e+00
cycle=10 time=3.7499999999996241e+01 dt=3.7499999999996252e+00
cycle=11 time=4.1249999999995865e+01 dt=3.7499999999996252e+00
cycle=12 time=4.4999999999995488e+01 dt=3.7499999999996252e+00
cycle=13 time=4.8749999999995111e+01 dt=3.7499999999996252e+00
cycle=14 time=5.2499999999994735e+01 dt=3.7499999999996252e+00
cycle=15 time=5.6249999999994358e+01 dt=3.7499999999996252e+00
cycle=16 time=5.9999999999993982e+01 dt=3.7499999999996252e+00
cycle=17 time=6.3749999999993605e+01 dt=3.7499999999996252e+00
cycle=18 time=6.7499999999993236e+01 dt=3.7499999999996252e+00
cycle=19 time=7.1249999999992866e+01 dt=3.7499999999996252e+00
cycle=20 time=7.4999999999992497e+01 dt=3.7499999999996252e+00
cycle=21 time=7.8749999999992127e+01 dt=3.7499999999996252e+00
cycle=22 time=8.2499999999991758e+01 dt=3.7499999999996252e+00
cycle=23 time=8.6249999999991388e+01 dt=3.7499999999996252e+00
cycle=24 time=8.9999999999991019e+01 dt=3.7499999999996252e+00
cycle=25 time=9.3749999999990649e+01 dt=3.7499999999996252e+00
cycle=26 time=9.7499999999990280e+01 dt=3.7499999999996252e+00
cycle=27 time=1.0124999999998991e+02 dt=3.7499999999996252e+00
cycle=28 time=1.0499999999998954e+02 dt=3.7499999999996252e+00
cycle=29 time=1.0874999999998917e+02 dt=3.7499999999996252e+00
cycle=30 time=1.1249999999998880e+02 dt=3.7499999999996252e+00
cycle=31 time=1.1624999999998843e+02 dt=3.7499999999996252e+00
cycle=32 time=1.1999999999998806e+02 dt=3.7499999999996252e+00
cycle=33 time=1.2374999999998769e+02 dt=3.7499999999996252e+00
cycle=34 time=1.2749999999998732e+02 dt=3.7499999999996252e+00
cycle=35 time=1.3124999999998695e+02 dt=3.7499999999996252e+00
cycle=36 time=1.3499999999998658e+02 dt=3.7499999999996252e+00
cycle=37 time=1.3874999999998622e+02 dt=3.7499999999996252e+00
cycle=38 time=1.4249999999998585e+02 dt=3.7499999999996252e+00
cycle=39 time=1.4624999999998548e+02 dt=3.7499999999996252e+00
cycle=40 time=1.4999999999998511e+02 dt=3.7499999999996252e+00
cycle=41 time=1.5374999999998474e+02 dt=3.7499999999996252e+00
cycle=42 time=1.5749999999998437e+02 dt=3.7499999999996252e+00
cycle=43 time=1.6124999999998400e+02 dt=3.7499999999996252e+00
cycle=44 time=1.6499999999998363e+02 dt=3.7499999999996252e+00
cycle=45 time=1.6874999999998326e+02 dt=3.7499999999996252e+00
cycle=46 time=1.7249999999998289e+02 dt=3.7499999999996252e+00
cycle=47 time=1.7624999999998252e+02 dt=3.7499999999996252e+00
cycle=48 time=1.7999999999998215e+02 dt=3.7499999999996252e+00
cycle=49 time=1.8374999999998178e+02 dt=3.7499999999996252e+00
cycle=50 time=1.8749999999998141e+02 dt=3.7499999999996252e+00
cycle=51 time=1.9124999999998104e+02 dt=3.7499999999996252e+00
cycle=52 time=1.9499999999998067e+02 dt=3.7499999999996252e+00
cycle=53 time=1.9874999999998030e+02 dt=3.7499999999996252e+00
cycle=54 time=2.0249999999997993e+02 dt=3.7499999999996252e+00
cycle=55 time=2.0624999999997956e+02 dt=3.7499999999996252e+00
cycle=56 time=2.0999999999997920e+02 dt=3.7499999999996252e+00
cycle=57 time=2.1374999999997883e+02 dt=3.7499999999996252e+00
cycle=58 time=2.1749999999997846e+02 dt=3.7499999999996252e+00
cycle=59 time=2.2124999999997809e+02 dt=3.7499999999996252e+00
cycle=60 time=2.2499999999997772e+02 dt=3.7499999999996252e+00
cycle=61 time=2.2874999999997735e+02 dt=3.7499999999996252e+00
cycle=62 time=2.3249999999997698e+02 dt=3.7499999999996252e+00
cycle=63 time=2.3624999999997661e+02 dt=3.7499999999996252e+00
cycle=64 time=2.3999999999997624e+02 dt=3.7499999999996252e+00
cycle=65 time=2.4374999999997587e+02 dt=3.7499999999996252e+00
cycle=66 time=2.4749999999997550e+02 dt=3.7499999999996252e+00
cycle=67 time=2.5124999999997513e+02 dt=3.7499999999996252e+00
cycle=68 time=2.5499999999997476e+02 dt=3.7499999999996252e+00
cycle=69 time=2.5874999999997436e+02 dt=3.7499999999996252e+00
cycle=70 time=2.6249999999997397e+02 dt=3.7499999999996252e+00
cycle=71 time=2.6624999999997357e+02 dt=3.7499999999996252e+00
cycle=72 time=2.6999999999997317e+02 dt=3.7499999999996252e+00
cycle=73 time=2.7374999999997277e+02 dt=3.7499999999996252e+00
cycle=74 time=2.7749999999997237e+02 dt=3.7499999999996252e+00
cycle=75 time=2.8124999999997198e+02 dt=3.7499999999996252e+00
cycle=76 time=2.8499999999997158e+02 dt=3.7499999999996252e+00
cycle=77 time=2.8874999999997118e+02 dt=3.7499999999996252e+00
cycle=78 time=2.9249999999997078e+02 dt=3.7499999999996252e+00
cycle=79 time=2.9624999999997038e+02 dt=3.7499999999996252e+00
cycle=80 time=2.9999999999996999e+02 dt=3.0013325158506632e-11
cycle=81 time=3.0000000000000000e+02 dt=6.0026650317013264e-11

Can you tell me what could be the reason?

[c-white]

By changing nx3, you changed from a 2D simulation to a 3D one.

I'm not sure what exactly is supposed to be wrong with the output. It looks like the code was requested to step at exactly time 300. Due to floating point rounding, each timestep is slightly less than 3.75, so the code reaches a time of 3 - 3e-11. The next step must be 3e-11 in order to end at exactly 3. The (hypothetical) step after that would be 6e-11, and then 1.2e-10, etc., doubling each time until it quickly reached 3.75. This is because then the code never takes a step more than twice its last step, to prevent numerical issues. Still, the simulation ended at time 3, so the next step size is largely irrelevant.

In any event, whenever trying to understand the output of the code, you should inspect your data to see what the state of the simulation is. The vis/python/plot_slice.py script is designed for this purpose.