Distribution being written but distro packed with 0's )-:

[rui-coelho]

Using the repo branch hotfix/48-bbnbi-beam-divergence-and-verification i ran both BBNBI5 and then ASCOT5. The BBNBI5 run as i mentioned before is now free of the beam "bending" but i'm facing something different: dist%5d or dist%rho5d just have zeros inside

1. See below the stdout i get from SLURM (i launched the ASCOT5 run from slum...).

ASCOT5_MAIN
Tag b562fa8-dirty
Branch hotfix/48-bbnbi-beam-divergence-and-verification

Initialized MPI, rank 0, size 1.

Reading and initializing input.

Input file is ascot5_200k_pietro.h5.

Reading options input.
Active QID is 1455651413
Options read and initialized.

Reading magnetic field input.
Active QID is 0194925115

2D magnetic field (B_2DS)
Grid: nR =  301 Rmin = 1.673 Rmax = 4.178
      nz =  301 zmin = -2.209 zmax = 2.200
Psi at magnetic axis (3.121 m, -0.004 m)
-1.377 (evaluated)
-1.377 (given)
Magnetic field on axis:
B_R = -0.000 B_phi = -2.166 B_z = 0.000
Estimated memory usage 11.1 MB
Magnetic field read and initialized.

Reading electric field input.
Active QID is 0106688659

Trivial Cartesian electric field (E_TC)
E_x = 0.000000e+00, E_y = 0.000000e+00, E_z = 0.000000e+00
Estimated memory usage 0.0 MB
Electric field read and initialized.

Reading plasma input.
Active QID is 0212218664

1D plasma profiles (P_1D)
Min rho = -0.00e+00, Max rho = 1.58e+00, Number of rho grid points = 601, Number of ion species = 2
Species Z/A  charge [e]/mass [amu] Density [m^-3] at Min/Max rho    Temperature [eV] at Min/Max rho
   1  /  2     1  /  2.014             7.29e+19/1.00e+10                9.04e+03/1.00e-02
   6  /  6     6  / 12.011             1.42e+18/1.00e+10                9.04e+03/1.00e-02
[electrons]   -1  /  0.001             7.99e+19/1.00e+10                8.27e+03/1.00e-02
Quasi-neutrality is (electron / ion charge density) 7.00
Estimated memory usage 0.0 MB
Plasma data read and initialized.

Reading neutral input.
Active QID is 2392108801

1D neutral density and temperature (N0_1D)
Grid:  nrho =    3   rhomin = 0.000   rhomax = 2.000
 Number of neutral species = 1
Species Z/A   (Maxwellian)
        1/  1 (1)
Estimated memory usage 0.0 MB
Neutral data read and initialized.

Reading wall input.
Active QID is 2061706387

3D wall model (wall_3D)
Number of wall elements 16560
Spanning xmin = -4.316 m, xmax = 4.316 m
         ymin = -4.316 m, ymax = 4.316 m
         zmin = -2.927 m, zmax = 3.069 m
Estimated memory usage 1.1 MB
Wall data read and initialized.

Reading boozer input.
Active QID is 2147410529

Boozer input
psi grid: n =    6 min = 0.000 max = 1.000
thetageo grid: n =   18
thetabzr grid: n =   10
Boozer data read and initialized.

Reading MHD input.
Active QID is 2941662183

MHD (stationary) input
Grid: nrho =    6 rhomin = 0.000 rhomax = 1.000

Modes:
(n,m) = ( 1, 3) Amplitude = 0.1 Frequency =   1 Phase =   0
(n,m) = ( 2, 4) Amplitude =   2 Frequency = 1.5 Phase = 0.785
Estimated memory usage 0.0 MB
MHD data read and initialized.

Reading atomic reaction input.
Active QID is 1923177355

Found data for 0 atomic reactions:
Estimated memory usage 0.0 MB
Atomic reaction data read and initialized.

Reading marker input.
Active QID is 1590347115

Loaded 199761 guiding centers.
Marker data read and initialized.

All input read and initialized.

Initializing marker states.
Estimated memory usage 1.5 MB.
Marker states initialized.

Preparing output.
Note: Output file ascot5_200k_pietro.h5 is already present.

The qid of this run is 1900133304

Inistate written.
Initialized diagnostics, 101.1 MB.
host: All fields initialized. Simulation begins, 48 threads.
host: Simulation complete.
gpu 0.000000 s, host 11920.335599 s
Endstate written.

Combining and writing diagnostics.

Writing diagnostics output.

Writing 5D distribution.

Writing rho 5D distribution.

Diagnostics output written.
Diagnostics written.

Summary of results:
   193721 markers had end condition Thermalization
      738 markers had end condition Aborted
     5302 markers had end condition Min energy

      738 markers were aborted with an error message:
          Input evaluation failed (marker could be outside input data grid)
          at line  404 in B_2DS.c

Done.

2. The data.ls() yields these Results:

   Results:
   run        1900133304 2023-12-01 18:45:02. [active]
   SD_NBI
   bbnbi      0007956606 2023-11-30 00:14:32.
   "BEAMS"

3. To check the distro was filled with zeros i did:

   
   dist=a5.data.SDNBI.getdist('rho5d')
   a5.input_init(bfield=True)
   mom1d = a5.data.active.getdist_moments(dist, "density")
   a5.input_free()
   mom1d.ordinate("density", toravg=True, polavg=True);
   fig = plt.figure(figsize=(4,4));ax1 = fig.add_subplot(1,1,1); mom1d.plot("density", axes=ax1)
   plt.show()

4. The actual marker distribution is just fine as i show in the skin and mileage figures...
![ASCOT_SD_summary_ekin](https://github.com/ascot4fusion/ascot5/assets/77900716/b347b811-18a4-462f-b895-8eb4bcf19593)
![ASCOT_SD_summary_mileage](https://github.com/ascot4fusion/ascot5/assets/77900716/274ed7ff-0bfc-4be3-baff-8fd88a05a046)

[miekkasarki]

Could be also that something's wrong in mom1d = a5.data.active.getdist_moments(dist, "density").

Do you also get zeros with:

dist=a5.data.SDNBI.getdist('5d')
mom2d = a5.data.active.getdist_moments(dist, "density")
mom2d.ordinate("density", toravg=True);
mom2d.plot("density")

and is np.sum(a5.data.SDNBI.getdist("5drho").distribution().ravel()) == 0?

[rui-coelho]

In [252]: np.sum(a5.data.SDNBI.getdist("rho5d").distribution().ravel()) == 0 Out[252]: array(True)

)-:

[miekkasarki]

Is this only for the "large" distributions? I wonder if there is some place that I missed in the fix...

[rui-coelho]

in this case it was a "small/medium" case since i used grids of size 60 except 1 for NPHI. When i was struggling i had NPHI also 60 but decided to go just for 1 since the plasma is in effect axisymmetric in the case of interest. Basically the same settings that Pietro used in ASCOT4 a couple of years back. I currently have a higher resolution case running (750k markers and +100 on grids except 1 on NPHI). Takes roughly 12.6hr to run at the GW.

[rui-coelho]

And the 200k marker grid-60 sized case in the end is a ~360Mb hdf5 file whereas the bbnbi5 output is just ~154Mb. So, data was definitively written and markers seem ok (even if for 738 markers i got 738 markers were aborted with an error message: Input evaluation failed (marker could be outside input data grid)

[rui-coelho]

the high resolution (grid~100 points and 750k markers) also ended up with

In [258]: np.sum(a5.data.SDNBI.getdist("5d").distribution().ravel())
Out[258]: unyt_quantity(0., 'particles*s/(degree*e*kg**2*m**4)')

I will try the tutorials to check if i get the same result....

[rui-coelho]

in the tutorial Dealing with distributions i did manage to get a 5d distribution not packed with 0's. I just tried a new run of bbnbi5 with 30 point grid with 1k markers and ascot5 with much lower mileage and cpu time and got the same list packed with 0's.....so it's clearly not related to the mesh size nor something of the sort...... Here is the options used:

BBNBI5

opt = Opt.get_default()
opt.update({
    "ENABLE_DIST_5D":1,
    "DIST_MIN_R":min_R,    "DIST_MAX_R":max_R,    "DIST_NBIN_R":30,
    "DIST_MIN_PHI":0,    "DIST_MAX_PHI":360,  "DIST_NBIN_PHI":1,
    "DIST_MIN_Z":min_Z,   "DIST_MAX_Z":max_Z,    "DIST_NBIN_Z":30,
    "DIST_MIN_PPA":-max_p, "DIST_MAX_PPA":max_p,  "DIST_NBIN_PPA":30,
    "DIST_MIN_PPE":0,    "DIST_MAX_PPE":max_p, "DIST_NBIN_PPE":30,
    "DIST_MIN_TIME":0,   "DIST_MAX_TIME":1.0, "DIST_NBIN_TIME":1,
    "DIST_MIN_CHARGE":0.5,   "DIST_MAX_CHARGE":1.5,  "DIST_NBIN_CHARGE":1
    })
a5.data.create_input("opt", **opt, desc="BEAMDEPOSITION")

ASCOT5

a5 = Ascot(ascot_filename)
ids = a5.data.active.getstate("ids", endcond="IONIZED")
mrk = a5.data.active.getstate_markers("gc", ids=ids)
a5.data.create_input('gc',**mrk,desc="PNBI4ASCOT",activate=True)
a5.data.create_input('E_TC',desc="DUMMY_EF",activate=True)
a5.data.create_input('Boozer',desc="DUMMY_BOOZER",activate=True)
a5.data.create_input('MHD_STAT',desc="DUMMY_MHD",activate=True)

#-----Update options for ASCOT simulation-------
#      
opt=a5.data.options.active.read()
opt.update({
    "ENABLE_DIST_RHO5D":1,
    "DIST_MIN_RHO":0,    "DIST_MAX_RHO":1,    "DIST_NBIN_RHO":30,
    "DIST_MIN_THETA":0,    "DIST_MAX_THETA":360,    "DIST_NBIN_THETA":20,
    # Simulation mode
    "SIM_MODE":2, "ENABLE_ADAPTIVE":1,
    # Setting max mileage above slowing-down time is a good safeguard to ensure
    # simulation finishes even with faulty inputs. Same with the CPU time limit.
    "ENDCOND_SIMTIMELIM":1, "ENDCOND_MAX_MILEAGE":0.01,
    "ENDCOND_CPUTIMELIM":1, "ENDCOND_MAX_CPUTIME":0.1e3,
    # The energy limit which separates a fast ion from thermal bulk is not well defined,
    # but we usually use Emin = 2 x Tion as the limit. Setting also a fixed minimum energy
    # is advised sine plasma temperature at the pedestal can be low.
    "ENDCOND_ENERGYLIM":1, "ENDCOND_MIN_ENERGY":2.0e3, "ENDCOND_MIN_THERMAL":2.0,
    # Physics
    "ENABLE_ORBIT_FOLLOWING":1, "ENABLE_COULOMB_COLLISIONS":1
    })
a5.data.create_input("opt", **opt, desc="ASCOT_OPTIONS",activate=True)

[rui-coelho]

Continuing my detective work, now i used a prescribed set of markers from the distributions tutorial i.e.

from a5py.ascot5io.marker import Marker
nmrk = 5000
mrk = Marker.generate("gc", n=nmrk, species="alpha")
mrk["energy"][:] = 85.0e3
mrk["pitch"][:]  = 0.99 - 1.98 * np.random.rand(nmrk,)
mrk["r"][:]      = 3.0 + 0.8 * np.random.rand(nmrk,)
a5.data.create_input("gc", **mrk)

merging it with the rest of the dataset instead of my bbnbi5 results. The SIM settings for CPU and MILEAGE were minimal to run for 80sec only and got 4086 markers had end condition Sim time limit and 914 markers had end condition Thermalization.

Still.......

In [6]: np.sum(a5.data.active.getdist("5d").distribution().ravel())
Out[6]: unyt_quantity(0., 'particles*s/(degree*e*kg**2*m**4)')

/&%$#"#$%......there must be something stupidly wrong

[rui-coelho]

Merged now also the plasma1D and the rudimentary 2D wall box...same &%$#"!"#$....

a5.data.create_input("plasma flat", density=1e21,activate=True)
a5.data.create_input("wall_2D",activate=True)

from a5py.ascot5io.marker import Marker
nmrk = 500
mrk = Marker.generate("gc", n=nmrk, species="alpha")
mrk["energy"][:] = 80.0e3
mrk["pitch"][:]  = 0.99 - 1.98 * np.random.rand(nmrk,)
mrk["r"][:]      = 2.5 + 1 * np.random.rand(nmrk,)
a5.data.create_input("gc", **mrk,activate=True)

Only thing left changing is the equilibrium......which should be the least of problems.....

[rui-coelho]

and with your a5.data.create_input("bfield analytical iter circular",activate=True) it works... this is very very suspicious.....the EQDSK i imported is a perfectly valid one as far as the mapping of the plasma1d to (R,Z) is concerned otherwise i would get no decent beam deposition with BBNBI5....nor marker evolution towards thermalisation with ASCOT5. There is likely some limitation/constraint to generate the distribution5d and its variants....

[rui-coelho]

revising my steps since i need to be consistent with the DIST_MIN_CHARGE":0.5, "DIST_MAX_CHARGE":1.5 only if using deuterium and i was using your markers of alphas and this is why i couldn't then get any distribution. Changing the range to 1.5,2.5 i could get non-zero values. This with your mock-up markers.

[rui-coelho]

Ok, tried so many options to see what the heck was going on. This is what triggers the "packed with 0s" and "regular" distributions. If i do NOT include the options for DIST_MIN_CHARGE and DIST_MAX_CHARGE (0.5 and 1.5) in the settings for ascot5 it returns 0's in the distribution. This is all a bit absurd for two reasons:

1. Even if i set DIST_MIN_CHARGE":0.5, "DIST_MAX_CHARGE":1.5 in BBNBI5 options, when accessing these same options in the object it actually returns 0 and 1 respectively. If the type is INT then the user should get clearer instructions on how/what to set when we wish to feature energetic particles of charge Z=1 (and what to do for another Z).
2. After BBNBI5 has ran, my interpretation is that the update method on the options object will update options that were previously set and written in the object of the hdf5 that is read. All options are actually set by the constructor so effectively we only update the options (and not add new ones...). Therefore, it is absurd that i should need to update again with DIST_MIN_CHARGE":0.5, "DIST_MAX_CHARGE":1.5 since these values should have already been set (and remained as such...).

Solution i found: rather than setting (updating) AGAIN the DIST_MIN_CHARGE":0.5, "DIST_MAX_CHARGE":1.5 on the options of ASCOT5 object, i managed to trick the "stupid code" by using instead DIST_MIN_CHARGE":1, "DIST_MAX_CHARGE":2 in the options of BBNBI5. Then these were inherited by ASCOT5 as expected and there was NO NEED to update them on the options of the object.

Request i make: the tutorials all contain the 0.5 and 1.5 values for the MIN and MAX charges. Unless there is a valid reason to do so (very misleading as i showed above), please make sure my solution is actually the correct way to proceed and in this case change the tutorials accordingly.

[miekkasarki]

Thanks for detective work! Using fractional values DIST_MIN_CHARGE":0.5, "DIST_MAX_CHARGE":1.5 should be supported. The reason why they weren't is that those parameters were erronously typecasted as integers before they were written to the HDF5 file. So in the code the interval would be open (1,2) and no physical particle could fit there.

I'll commit the fix and close this issue once it has been merged.

[rui-coelho]

so one should expect in the future release to get the same result either using (0.5 and 1.5) or (1 and 2) ? Right now it does work if i use 1 and 2 as limits since i the lower value is retained. Will this change ?

[miekkasarki]

Yeah (1,2) should work even after the fix but note that this is not really supported. The reason is that in the code, the correct bin is determined like this:

i_q[i] = floor((p_f->charge[i]/CONST_E - dist->min_q) / ((dist->max_q - dist->min_q) / dist->n_q));

So while in principle the interval is semi-open [1,2), the floating point arithmetic could lead to p_f->charge[i]/CONST_E being either X.000...1 or X.999...9. I guess this could be platform dependent or depend on test particle charge. This is why I'd advise to treat (min_q, max_q) as open interval.

[rui-coelho]

fair enough. Once fixed please mind updating the documentation to highlight this "feature" with simple examples e.g. just D, just Alpha and multi species.