ascot4fusion / ascot5 Goto Github PK

The keyboard controls provided by VTK/Pyvista are very unintuitive. We need something like what FPS games have.

I've already implemented this and I'll proceed with merging. However, there is an issue which anyone is welcome to fix. I couldn't figure out how to toggle off (all) the pre-existing commands meaning I had to use whatever keys where free instead what I'd wanted (wasd + arrow keys).

So the current bindings are:

    p.add_key_event('w', lambda : control_camera('move_forward'))
    p.add_key_event('s', lambda : control_camera('move_backward'))
    p.add_key_event('a', lambda : control_camera('move_left'))
    p.add_key_event('d', lambda : control_camera('move_right'))
    p.add_key_event('n', lambda : control_camera('move_up'))
    p.add_key_event('m', lambda : control_camera('move_down'))
    p.add_key_event('r', lambda : control_camera('rotate_cw'))
    p.add_key_event('y', lambda : control_camera('rotate_ccw'))
    p.add_key_event('t', lambda : control_camera('look_up'))
    p.add_key_event('g', lambda : control_camera('look_down'))
    p.add_key_event('f', lambda : control_camera('look_left'))
    p.add_key_event('h', lambda : control_camera('look_right'))

A strange bug that seems compiler/platform specific. I've removed the first data point in GC transform test because of this.

This test simulates GOs but collects GC data (using the option parameter). When I run this on my desktop the test runs properly and all collected orbit points are valid. When it is run on Github action the test fails because the first collected orbit point contains dummy data (e.g. R=0 which leads to NaNs in postprocessing).

This is a minor bug (no hotfix needed) but I just want to record it here.

As i recall, to have a complete experience using the ASCOT5 environment we are currently forced to compile with MPI=0, This translates to reasonable running times (~3.5h using low resolution 60 point grid in r,z,vpar,vperp and 200k markers) in a single node (48 cores) at the Gateway as long as the beam energies are low i.e. 85keV. However, the same run using 500keV beams is telling me it will require something like ~77h (!!!)

My end conditions are "safe" to get as much as possible all markers ionised (i used the same for both 85 and 500keV beams) i.e.

    "ENDCOND_SIMTIMELIM":1, "ENDCOND_MAX_MILEAGE":0.7,
    "ENDCOND_CPUTIMELIM":1, "ENDCOND_MAX_CPUTIME":1.7e3,
    "ENDCOND_ENERGYLIM":1, "ENDCOND_MIN_ENERGY":2.0e3, "ENDCOND_MIN_THERMAL":1.5,

Related question: is the CPU time used to simulate each marker saved somewhere in the state ? For the 85keV beams the mileage is below 200ms (peak at ~50ms). But i don't know how much CPU it took unless i estimate it from the stdout but that's very imprecise (not enough digits in the fractional hours spent and not very meaningful anyway since don't discriminate end condition reached)...

The original ASCOT5 was written for Xeon PHI with quite gigantic offload kernels. The typical GPU offloading software uses numerous small kernel calls. Let us try to split the kernels.

Allocating large distributions fail and leads to segfault since the array sizes (and indexing) are ints. These should be replaced by size_t instead.

For CI physics tests we should have one for AFSI as well. I think a valid and easy test would be to have uniform DT plasma and then compare the AFSI result to the analytical one e.g. from NRL. The test should verify i) number of reactions in a volume ii) energy distribution and iii) pitch distribution.

For beam-thermal and beam-beam we would just provide a Maxwellian distribution as an input 5D distribution.

Right now BBNBI ignores the halo component and assumes isotropic divergence using only the horizontal value. All relevant parameters are in the input already, so implementing this is just a few lines of code. (I label this as a bug right now.)

For the verification test, we could have a single beamlet and verify from the beam deposition that the divergence and halo were properly implemented (the distribution should be bivariate double Gaussian).

We could also inject the beam into an uniform plasma and using the mean-free-path for the beam ionization we would have analytical estimate for the beam deposition.

Usually one wishes to present the distribution function in velocity space (Vpar,Vperp) and/or (R,Z) space. However, it is clear from the documentation that dist5d comes in (R, phi, z, ppar, pperp)...so in momentum rather than velocity space. Now, if when running ASCOT we use a SINGLE mass species in each one of the injectors used, there is no problem in converting to velocity space. If however i use 8 injectors, some with D and others with T beam ions, i have no way to map to velocity space since the dist5d knows nothing about the markers !
Incidentally, i note that dist5d "raw" is actually 8d, possibly (?,nR,nPhi,nZ,nPpa,nPpe,nTime,nCharge) but the first dimension (typically one) i have no clue what it is, nor i am sure if the nTime and nCharge are in the correct order....

Help.....?!

1. Write a conda recipe that builds ascot5_main (and other binaries), libascot.so and installs a5py along with any dependencies.

• This requires that path to libascot.so is not hard-coded (but I guess we want to retain hard-coded path when using ascot from source or not using it in conda)
• Any non-code files listed in MANIFEST.IN and accessed via setuptools during runtime.
• All libraries should be dynamically linked

2. Test that the resulting tarball can be used to install ASCOT5 in other platforms (we can probably ignore Windows for now).
3. GitHub action to build the tarball (but not upload it except as a temporary artifact).
4. Make a PR to conda-forge to share the recipe and distribute via official channels.

Release 5.5

• Source code made open and is now hosted on GitHub.
• Major overhaul of the Python interface a5py.
• Online documentation and tutorials now available.
• New GUI (usable but unstable).
• BBNBI integrated to the rest of the code.
• AFSI is now available.
• BioSaw integrated to the rest of the code.
• Implemented constant of motion distribution output.

The conversion between different COCOS is not working in the current ascot5 version. This is important when loading eqdsk files with a COCOS that ASCOT cannot use, i.e. any COCOS other than 3.

I implemented a fix by adding a fix in a5py/templates/importdata.py (line 228)

and a new COCOS class and functions in a5py/physlib/cocos.py

def fromCocosNtoCocosM(eqd, cocos_m, phiclockwise=None, weberperrad=None):
    """Transform equilibrium from cocos_n (determined from eqd, see below) to cocos_m.

    Parameters
    ----------
    eqd : dict
        Dictionary from reading the EQDSK file.
    cocosm : int
        Target COCOS.

    Returns
    -------
    eqdout : dict
        Equilibrium data converted to cocos_m.
    """

    cocos_n = assign(eqd["qpsi"][0], eqd["cpasma"], eqd["bcentr"], eqd["simagx"], eqd["sibdry"], phiclockwise, weberperrad)

    transform_dict = transform_cocos(cocos(cocos_n), cocos(cocos_m))

    # Define output
    eqdout = copy.deepcopy(eqd)
    # eqdout["nx"] = eqd["nx"] # For clarity
    # eqdout["ny"] = eqd["ny"] # -||-
    eqdout["rdim"] = eqd["rdim"]*transform_dict["R"]
    eqdout["zdim"] = eqd["zdim"]*transform_dict["Z"]
    eqdout["rcentr"] = eqd["rcentr"]*transform_dict["R"]
    eqdout["rleft"] = eqd["rleft"]*transform_dict["R"]
    eqdout["zmid"] = eqd["zmid"]*transform_dict["Z"]
    eqdout["rmagx"] = eqd["rmagx"]*transform_dict["R"]
    eqdout["zmagx"] = eqd["zmagx"]*transform_dict["Z"]
    eqdout["simagx"] = eqd["simagx"]*transform_dict["PSI"]
    eqdout["sibdry"] = eqd["sibdry"]*transform_dict["PSI"]
    eqdout["bcentr"] = eqd["bcentr"]*transform_dict["B"]
    eqdout["cpasma"] = eqd["cpasma"]*transform_dict["I"]
    eqdout["fpol"] = eqd["fpol"]*transform_dict["F"]
    eqdout["pres"] = eqd["pres"]*transform_dict["PRES"]
    eqdout["ffprime"] = eqd["ffprime"]*transform_dict["FFPRIME"]
    eqdout["pprime"] = eqd["pprime"]*transform_dict["PPRIME"]
    eqdout["psi"] = eqd["psi"]*transform_dict["PSI"]
    eqdout["qpsi"] = eqd["qpsi"]*transform_dict["Q"]
    #eqdout["psiaxis"] = psiaxis # Wrong key? psiaxis is called simagx according to freeqdsk docs
    #eqdout["psiedge"] = psiedge # Wrong key? psiedge is called sibdry according to freeqdsk docs
    #eqdout["nbdry"] = eqd["nbdry"]
    #eqdout["nlim"] = eqd["nlim"]
    eqdout["rbdry"] = eqd["rbdry"]*transform_dict["R"]
    eqdout["zbdry"] = eqd["zbdry"]*transform_dict["Z"]
    eqdout["rlim"] = eqd["rlim"]*transform_dict["R"]
    eqdout["zlim"] = eqd["zlim"]*transform_dict["Z"]

    return eqdout

def transform_cocos(cc_in: COCOS, cc_out: COCOS, sigma_Ip = None, sigma_B0 = None, ld = (1, 1), lB = (1, 1), exp_mu0 = (0, 0)):
    """
    Returns a dictionary of the multiplicative factors to transform COCOS from `cc_in` to `cc_out`

    Parameters
    ----------
    sigma_Ip : Union[Tuple[int, int], None]
        A tuple of the (Input, Output) current sign or nothing
    sigma_B0 : Union{NTuple{2,Int},Nothing}` - A tuple of the (Input, Output) toroidal field sign or nothing
    ld : NTuple{2,<:Real}
        A tuple of the (Input, Output) length scale factor. Default = (1,1)
    lB : NTuple{2,<:Real}
        A tuple of the (Input, Output) magnetic field scale factor. Default = (1,1)
    exp_mu0 : NTuple{2,<:Real}
        A tuple of the (Input, Output) mu0 exponent (0, 1). Default = (0,0)

    Returns
    -------
    transforms : dict
        Transform multiplicative factors to be able to convert from `cc_in to `cc_out`
    """

    ld_eff = ld[1] / ld[0]
    lB_eff = lB[1] / lB[0]
    exp_mu0_eff = exp_mu0[1] - exp_mu0[0]

    sigma_RpZ_eff = cc_in.sigma_RpZ * cc_out.sigma_RpZ

    if sigma_Ip is None:
        sigma_Ip_eff = cc_in.sigma_RpZ * cc_out.sigma_RpZ
    else:
        sigma_Ip_eff = sigma_Ip[0] * sigma_Ip[1]

    if sigma_B0 is None:
        sigma_B0_eff = cc_in.sigma_RpZ * cc_out.sigma_RpZ
    else:
        sigma_B0_eff = sigma_B0[0] * sigma_B0[1]

    sigma_Bp_eff = cc_in.sigma_Bp * cc_out.sigma_Bp
    exp_Bp_eff = cc_out.exp_Bp - cc_in.exp_Bp
    sigma_rhotp_eff = cc_in.sigma_rhotp * cc_out.sigma_rhotp

    mu0 = 4 * 3.14159265358979323846 * 1e-7  # pi is used directly for more precision

    transforms = {}
    transforms["R"] = ld_eff
    transforms["Z"] = ld_eff
    transforms["PRES"] = (lB_eff ** 2) / (mu0 ** exp_mu0_eff)
    transforms["PSI"] = lB_eff * (ld_eff ** 2) * sigma_Ip_eff * sigma_Bp_eff * ((2 * 3.14159265358979323846) ** exp_Bp_eff) * (ld_eff ** 2) * lB_eff
    transforms["TOR"] = lB_eff * (ld_eff ** 2) * sigma_B0_eff
    transforms["PPRIME"] = (lB_eff / ((ld_eff ** 2) * (mu0 ** exp_mu0_eff))) * sigma_Ip_eff * sigma_Bp_eff / ((2 * 3.14159265358979323846) ** exp_Bp_eff)
    transforms["FFPRIME"] = lB_eff * sigma_Ip_eff * sigma_Bp_eff / ((2 * 3.14159265358979323846) ** exp_Bp_eff)
    transforms["B"] = lB_eff * sigma_B0_eff
    transforms["F"] = sigma_B0_eff * ld_eff * lB_eff
    transforms["I"] = sigma_Ip_eff * ld_eff * lB_eff / (mu0 ** exp_mu0_eff)
    transforms["J"] = sigma_Ip_eff * lB_eff / ((mu0 ** exp_mu0_eff) * ld_eff)
    transforms["Q"] = sigma_Ip_eff * sigma_B0_eff * sigma_rhotp_eff

    return transforms

The fix is based on the equations that describe how to transform between any COCOS in O. Sauter et al, CPC, 2013.
I will create a new branch locally and submit a pull request to fix this issue, in line with the developer docs found in https://ascot4fusion.github.io/ascot5/developing.html.

Please let me know if you have any comments!

rho5d distribution from bbnbi5 only populates the first theta bin, even though the bbnbi resulting markers would suggest other theta bins should also be populated.

Does the current bbnbi not support rho5d? If so, it would be nice to be implemented.

For the 3D wall it is possible to provide flags for wall elements e.g. to separate divertor in post-processing. So far these have not been used in the actual simulation but this will change with BMC, where flag indicates the elements of interest.

This issue implements:

• Reading wall flags from HDF5
• Function to return flag of a given element
• Flags for the 2D wall as well

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...
   
   

Now all MPI processes are writing the input sanity checks to the output log which makes it messy and unreadable when there are several processes. However, fixing this issue requires passing the MPI rank to the init_offload routines which would require quite a lot of changes so need to think this through.

I've also noticed that the simulation progress file was not written with VERBOSE=2 at least on the system where I was running. Need to investigate this.

One missing feature from ASCOT4 is the background plasma rotation and its effect on beam ionization, coulomb collisions and beam-thermal fusion.

For the current 1D plasma inputs, the usual way of representing rotation as toroidal and poloidal angular velocities should be sufficient, but we probably want to include the possibility of specifying arbitrary plasma flows with 2D and 3D plasma inputs. The interfaces and implementation should be done with this in mind so that the rotation can be included in all velocity components.

When generating markers for ascot (from a bbnbi results distribution) using a5.markergen.generate as follows:
nmrk = 4000
anum = 1
znum = 1
mass = 1.007unyt.amu
charge = 1.0unyt.e
dist = a5_run_tcv.data.active.getdist("5d")
dist.integrate(charge=np.s_[:], time=np.s_[:])
mrk = a5_run_tcv.markergen.generate(nmrk, mass, charge, anum, znum
plt.scatter(mrk["pitch"], mrk["energy"])

The markers pitch and energy distribution appears to be changed from the bbnbi5 result markers:
mrk = a5_run_tcv.data.active.getstate_markers("gc")
plt.scatter(mrk["pitch"], mrk["energy"])

Atomic reactions are currently divided between three modules:

• suzuki which calculates beam-stopping cross section analytically
• asigma which calculates any cross-sections using ADAS
• atomic which is used in simulation to alter the charge state based on values given by asigma.

Issue is that Suzuki is not used by the atomic module making it impossible to use the Suzuki BMS data and e.g. OpenADAS CX data simultaneously. These should be combined somehow.

Hi ASCOT Team,

I am currently testing ASCOT5 on Ubuntu 22.04.1 LTS within a virtual machine. During my initial test, while running the Live Simulations part of the introduction under the ascotenv environment, I encountered an AttributeError.

AttributeError Traceback (most recent call last)
Cell In[18], line 10
8 # Marker input can be anything, but here we just use the one from ascot.h5
9 mrk = a5.data.marker.active.read()
---> 10 a5.simulation_initmarkers(**mrk)
12 # Options input can also be anything, but here we just use the one from ascot.h5
13 opt = a5.data.options.active.read()

File ~/ascot5/ascot5/a5py/ascotpy/libsimulate.py:231, in LibSimulate.simulation_initmarkers(self, **mrk)
229 p.pitch = mrk["pitch"][i]
230 p.zeta = mrk["zeta"][i]
--> 231 p.mass = mrk["mass"][i] * unyt.atomic_mass_unit.value
232 p.charge = mrk["charge"][i] * unyt.elementary_charge.value
233 p.anum = mrk["anum"][i]

AttributeError: 'Unit' object has no attribute 'value'

Why did this happen? Does the unyt package have a wrongly defined method for unyt.atomic_mass_unit.value? I hope you can shed light on the reason.

Thanks,
Saic

Impurity transport can be modelled with the atomic module once the following issues have been implemented:

• In atomic.c, the atomic_rates() function should deal separately cases with q=0, q=znum and 0 < q < znum.
• atomic_rates() function should calculate effective rates collectively for all reactions that increase the charge state and for all reactions that decrease the charge state.
• Proper treatment of partially ionized plasma species in atomic_rates().
• asigma.c and asigma_loc.c should work as is or with minor changes.
• Neutral module should include anum and znum of the species (propably not needed but can be done simultaneously).

Nobody is currently working on this, so this is just for the record.

Add new option to a5manage to show contents of the file:

a5manage -s ascot.h5

Equivalent to a5.data.ls().

a5manage -s ascot.h5 bfield

Equivalent to a5.data.bfield.ls().

Let's have something to guide people when we open the repo, but before we have set the 5.5 branch properly.

@jvarje what issue/branch is this?

Updating this list as new features are merged to develop.

Release notes

• Made the code Conda compatible:
  • Included environment.yaml for creating an environment where ASCOT5 can be built and run.
  • Updated installation instructions to make use of Conda.
  • All dependencies are found in conda-forge.
• Several new import templates:
  • wall data from EQDSK and VTK
  • plasma profiles from text files with option to extrapolate
  • perturbation from MARS-F
• a5manage can now display the contents of the file.
• Camera controls for the interactivate 3D wall (load) plot.
• New plots and improved some others.

The ALD force for synchrotron radiation for RE studies.

It would be nice to implement a new flag logscale, to show the distributions in logarithmic scale.

Implement support for interfacing ASCOT5 with Stellerator-Tools to generate magnetic field components needed for ASCOT5 simulations on an extended (beyond the LCFS) cylindrical grid.

Using VMEC + MAKEGRID + BMW this can be achieved, and would be nice to have interfaced tools that do all this directly from a5py.

5.5 release is mostly focused on Python side of things. New features are:

• 1D/2D distributions evaluated in post-processing from 5D distributions.
• Constant-Of-Motion (COM) distribution collected in run-time.
• Ascot and Ascotpy combined as a single entity "Ascot".
  • Seamless interface to libascot.so.
  • Can run simulations from Python.
  • Many pre/post-processing scripts and tools are now methods of this object.
  • Inputs can be generated via templates.
  • Far from complete, but sets the framework for further development.
• Document generation from source and tutorials.
• New GUI.

This update will break many existing scripts that users have, but here's how to fix most of them.

# Replace this line
from a5py.ascot5io.ascot5 import Ascot
# with this
from a5py import Ascot

# Now the "data" attribute in new Ascot object works in a very similar way as the old Ascot object
a5 = Ascot("ascot.h5")
a5.data.ls()

# You don't need to import Ascotpy anymore so remove these lines
from a5py.ascotpy import Ascotpy
# and use the Ascot object instead:
import unyt
a5.input_init(bfield=True)

# Notice that now multiple quantities can be evaluated simultaneously
bphi, rho = a5.input_eval(6.2*unyt.m, 0*unyt.deg, 0*unyt.m, 0*unyt.s, "bphi", "rho")
a5.input_free()

# Previously output data was accessed like this:
endr = run.endstate.get("r")
orbr = run.orbit.get("r")

# Now output is accessed from the run node (use mode="prt" for full-orbit quantities)
endr, endz = a5.data.active.getstate("r", "z", mode="gc", state="end")
orbr, orbz = a5.data.active.getorbit("r", "z", ids=1)

# Some output requires magnetic field or other inputs, so initialize them first if you get an error
a5.input_init(run=True, bfield=True)
inipsi = a5.data.active.getstate("psi")
a5.input_free()

# Run and Ascot objects contain more useful functions for evaluation and plotting which you can use to simplify your scripts!

# Inputs were previously created like this
from a5py.ascot5io.mrk_gc import write_hdf5
write_hdf5("ascot.h5", **markerdata, desc="Markers")
# Now this is done via Ascot object
a5.data.create_input("gc", **markerdata, desc="Markers", activate=True)

# Premade inputs and importing data are implemented as templates. E.g. here is how to construct and write boozer data
a5.input_init(bfield=True)
a5.data.create_input("boozer_tokamak", desc="Boozer data", activate=True)
a5.input_free()

Setting logscale=True in method hist2D currently does nothing.

The sanity checks (i.e. the stuff printed in simulation log, see below for example) for plasma 1DS and 1Dt are plotting non-sensical values and could potentially result in segfault. This bug is not present in plasma_1D input and even in 1DS and 1Dt does not affect the simulation results in any way. The bug is due to someone mindlessly copy-pasting the sanity check from the 1D case so all the pointers etc. used to print the data are off.

1D plasma profiles (P_1DS)
Min rho = 0.00e+00, Max rho = 1.00e+02, Number of rho grid points = 3, Number of ion species = 1
Species Z/A  charge [e]/mass [amu] Density [m^-3] at Min/Max rho    Temperature [eV] at Min/Max rho
   1  /  1     1  /  1.000             4.61e+01/4.61e+01                -2.22e+20/-2.22e+20       
[electrons]   -1  /  0.001             4.61e+01/4.61e+01                -2.22e+20/-2.22e+20       
Quasi-neutrality is (electron / ion charge density) 1.00
Estimated memory usage 0.0 MB

It would be nice if Python would use an exponential fit to extrapolate plasma profiles.

TRANSP outputs the plasma and field data that it uses internally during a given time window, in a "plasma state file", provided that the namelist setting FI_OUTTIM is set. The markers generated by TRANSP's NBI module, NUBEAM, are output in the marker birth file which is part of standard output.

It should be possible to use this to generate completely self-consistent input for ASCOT5 from TRANSP.
(currently trying to fork and add a feature branch to include existing scripts which do this)

The BBNBI output now consists of ionized markers, shinethrough, and the NBI ion birth rate as a 5D distribution. Do we also need neutral density, i.e., should we collect the distribution data while the beam neutrals are in flight?

The goal of this work is to implement a time-dependent magnetic (and electric) field background in ASCOT5.
The main issue here is how I am going to interpolate a 3D magnetic field along the time variable, B(R,phi,z,t). My idea is to extend the cubic spline interpolation methods from PSPLINE to a 4 parameters cubic spline. This is something straightforward but should be done with care, because this hasn't been constructed in the PSPLINE library. Both "compact" and "explicit" splines can be extended to 4 parameters. A comparison between the two methods can be of interest.

  The "explicit" method is claimed to be faster but the memory consumption increases a lot in 4D (N = 256N_RN_phiN_zN_t)

    
  The "compact" method is claimed to be somewhat slower but it doesn't consume so much memory (N = 16N_RN_phiN_zN_t).

Does anybody know what differences I can expect in terms of accuracy?
I am also still missing good control cases against which I can benchmark this interpolation method. With this, I mean problems (ideally, with an analytical solution) to check the errors that I am making with this approach.
I am thinking that one way could be to test this new module against the MHD module in ASCOT4. However, such test would be unfair for my module because in the MHD module the perturbations come from analytical equations and doesn't really interpolate anything. Therefore, I consider this a "test case", rather than a "control case".

I have implemented the possibility to load wall data from a G-EQDSK file (magnetic equilibrium files with file endings .eqdsk, or sometimes .geqdsk). G-EQDSK files (almost always) contain (R,z) points that define the first wall. You can therefore choose to load this data as a wall_2D or a wall_3D (by rotating the 2D data around toroidally). This can be specified via the return3Dwall keyword argument. The data is loaded into a wall_3D by default.

You use this new feature as follows:
**define your ASCOT object as a5**
a5.data.create_input("wall_geqdsk",fn="/path/to/your/G-EQDSK/file.eqdsk") # will be loaded as a wall_3D

Or, if you would rather have a wall_2D
**define your ASCOT object as a5**
a5.data.create_input("wall_geqdsk",fn="/path/to/your/G-EQDSK/file.eqdsk", return3Dwall=False) # will be loaded as a wall_2D

P.S. COCOS is taken into account. As a sidenote, when loading .eqdsk files (in general, not just walls), you have to specify explicitly if the .eqdsk file has a phi-direction clockwise (viewed from above). You do this via the phiclockwise keyword argument. You don't have to specify explicitly if the .eqdsk file has a phi-direction counter-clockwise (viewed from above)(phiclockwise=False is assumed by default, as it is true for most magnetic equilibria). Continuing, you also have to specify explicitly if the .eqdsk file has a poloidal flux that has NOT been divided by 2pi. You do this by setting the weberperrad keyword argument to False (True by default, since most magnetic equilibria have a poloidal flux that HAS been divided by 2pi). Therefore, the above a5.data.create_input() calls would look like this in these cases:

Phi-direction in the clockwise direction (viewed from above):
**define your ASCOT object as a5**
a5.data.create_input("wall_geqdsk",fn="/path/to/your/G-EQDSK/file.eqdsk", phiclockwise=True)

Poloidal flux NOT already divided by 2*pi:
**define your ASCOT object as a5**
a5.data.create_input("wall_geqdsk",fn="/path/to/your/G-EQDSK/file.eqdsk", weberperrad=False)

Phi-direction in the clockwise direction (viewed from above) AND poloidal flux NOT already divided by 2pi:
**define your ASCOT object as a5**
a5.data.create_input("wall_geqdsk",fn="/path/to/your/G-EQDSK/file.eqdsk", phiclockwise=True, weberperrad=False)

Right now Python repository is just located in folder a5py, but we should merge these two repositories properly. I propose following file structure:

ascot5
|-> src/ (C-repo)
|-> a5py/ (__ init __.py for a5py)
|-> docs/
|-> bin/
|-> scripts/ (wild west for Ascot5 related scripts; TBD)
|-> setup.py
|-> Makefile etc.

Weights should actually be in 1/s, not in 1.

Hi ASCOT Team,

I am testing ASCOT5 on a MacBook Pro with Apple Silicon running Sonoma 14.0, for this I had to use GNU compilers V12, instead of V10 as suggested in the documentation. GNU compilers V10 won't install well using MacPorts.
The issue popping up while compiling is that gcc won't find immintrin.h library in my system, this library is used in ascot5/simulate/simulate_.c code files. If I comment the line for including this library in all ascot5/simulate/simulate_.c files, the error goes away and I am able to build ascot5_main and libascot, even run a test simulation.
The question is: is immintrin.h used in some way so I shouldn't be dropping that dependency? or is OK to don't use it?
Thanks.
Leo

I think there is something fishy with the dict_5d units and tools to integrate/plot.

d = a5.data.active.getdist("5d")
In [44]: d.distribution().units
Out[44]: particles*s/(degree*e*kg**2*m**4)

1. I may understand the degree since it stands for the phi coordinate (let alone jacobian issues this poses when then integrating over (R,Z) since i suspect integrating over Phi for you is just multiply by 2Pi if the distribution is axisymmetric...).
2. But i don't understand what the electron charge e is doing there. At most, since you are using the particle momentum (and not the velocity) for the velocity phase space, i could expect an atomic number but the effective mass should be already included in the value of the distribution so even this is a bit weird.....
3. When we integrate the d.distribution() over the PHI, PPAR, PPER, CHARGE or over only PHI, PPAR, PPER we get the same result which doesn't make sense unless the charge is actually the charge number (Z=1 in my case) but again i don't get it why it is part of the units of the distribution in the first place...
4. The units are very much different from what ASCOT4 uses since when i grab a HDF5 output file from a run i was given in the past by Pietro/Matteo i see that the units are s^2/m^5 which makes much more sense to me. They are the units of the 6d distribution BUT already multiplied by the Vperp units resulting from the integration to go from 3d velocity to Vpar and Vperp. But the units in ASCOT5, using the same logic, if i would use P instead of V should be, unless i'm wrong, s^2/(Kg^2*m^5).
5. Using the very same post processing script to analyse the HDF5 (using h5py and not the toolset you developed), incidentally, the distribution function over (R,Z) and over (Vpar,Vperp) comes at values which are of the order (electron charge / proton mass) higher that the values i obtained from the old ASCOT4 data.

So....enlighten me if i'm clearly missing something.... )-:

While GUI works, it's really easy to make it crash and there are plots that don't show anything and buttons that don't do anything.

I'll make an issue here for fixing these things, and for others to report any new problems or provide feature requests. I'll also make a branch that I will merge from time to time since this GUI is more like a hobby project which I don't expect to be finished soon.

As the title says, this feature is missing. My plan is to implement this in post-processing as the endstate store the marker position(on the wall) and the velocity vector. The angle of incidence is defined with respect to the surface normal, but here we can decide that the normal always points in the same direction as where the marker came. (The other option would be to fix the order in which triangle vertices must be given in input).

BBNBI currently uses Suzuki and only Suzuki. To make it use ADAS data #29 must be resolved first. Then all that is needed are few changes to the BBNBI simulation loop.

Since ADAS data is really needed only for He beams or plasma, we should also decide how BBNBI operates. Does it stop at the ionization or trace markers until they are fully ionized?

Having psi0 different than what is evaluated with splines is a constant headache. Here's what we need to do:

1. Read EQDSK and convert it to B_2DS dictionary
2. Initialize splines using the dictionary
3. Find psi0 with gradien descent method (in libascot.so)
4. Add little tiny bit of padding and replace psi0 in the dictionary
5. Write the data to HDF5

This of course should be an optional argument in the import_geqdsk template.

With many markers ( more than 1e5 ) it is useful to compress the datasets.

Feature to plot the 3D wall loads by adding the option of only markers, to illustrate the statistics in each triangle.
Feature to plot only interested 3D walls.
Feature to plot heat loads from only specific markers.
Bug fix: wall.number_of_elements was never set, and therefore wall.getNumberOfElements() always returns None.

I understand the concept of the "stupid code" being straightforward in it's output but, again, when trying to profit from the completeness of the HDF5 file output i find myself empty handed regarding some fundamental equilibrium data (not in the EQDSK) that is useful downstream: the toroidal magnetic field on axis ! This is useful since it is a common used normaliser e.g. omega_c for a given EP in case the Pphi comes in eV (i know...i know...not SI...). Of course i could get the R,Z of the axis (available in bfield group) and interpolate bphi (also available in bfield) on the magnetic axis.......but surely the result would depend on MY interpolator and not correspond on what internally the stupid code might actually have used.......

comments ?

ASCOT5 should be an IMAS actor or so.

The ITER The Integrated Modelling & Analysis Suite (IMAS)

This issue is for merging (and finishing) the BMC code to the main.