FUSE Introductory Tutorial

Download this tutorial from the FuseExamples repository

Basic concepts

To make sense of this tutorial, you'll need to know the following organization concepts of FUSE:

  1. 📂 Data storage: All data is stored in the dd structure, which follows the ITER IMAS ontology.
  2. 🧠 Actors: The core components of FUSE simulations are physics and engineering actors.
  3. 🕹️ Control: Actor functionality is governed by act parameters.
  4. 🚀 Initialization: The data structure can be initialized from 0D ini parameters.
  5. 🔧 Use cases: FUSE includes templates for various machines (e.g., FPP, ITER, ARC).
  6. 🔄 Workflows: Self-contained studies and optimizations are conducted via workflows, typically involving multiple FUSE simulations.
  7. 🌍 Interoperability: FUSE interfaces with existing modeling tools like OMFIT/OMAS and the IMAS ecosystem.

A diagram illustrating these concepts is provided below: image.png

Let's get started!


NOTE: Julia is a Just In Time (JIT) programming language. The first time something is executed it will take longer because of the compilation process. Subsequent calls the the same code will be blazingly fast.


Import the necessary packages

using Plots # for plotting
using FUSE # this will also import IMAS in the current namespace

Starting from a use-case

FUSE comes with some predefined use-cases, some of which are used for regression testing.

Note that some use cases are for non-nuclear experiments and certain Actors like Blankets or BalanceOfPlant will not perform any actions.

Here's the list of supported use-cases. These can be customized and you will also be able to build your own.

methods(FUSE.case_parameters)
# 26 methods for generic function case_parameters from FUSE:

Get initial parameters (ini) and actions (act) for a given use-case, let's use KDEMO for example

ini, act = FUSE.case_parameters(:KDEMO);

The ini data structure contains 0D parameters that will be used to bootstrap the dd with plausible data.

The ini parameters can be modified.

ini.equilibrium.B0 = 7.8
ini.equilibrium.R0 = 6.5;

The act data structure contains parameters that define how the actors (ie the models) will behave.

The act parameters can also be modified.

act.ActorCoreTransport.model = :FluxMatcher;

ini and act can now be used to initialize the data dictionary (dd) using the 0D parameters.

NOTE: init() does not return a self-consistent solution, just a plausible starting point to initialize our simulations!

dd = IMAS.dd() # an empty dd
FUSE.init(dd, ini, act);
actors: Equilibrium
actors:  TEQUILA
actors: CXbuild
actors: HCD
actors:  SimpleEC
actors:  SimpleIC
actors:  NeutralFueling
actors: Current
actors:  QED
actors: PassiveStructures

Let's see what we got

plot(dd.build)

plot(dd.equilibrium)

plot(dd.core_profiles)

plot(dd.core_sources)
Example block output

We can @checkin and @checkout variables with an associated tag.

This is handy to save and restore (checkpoint) our progress without having to always start from scratch (we'll use this later).

@checkin :init dd ini act

Running Actors

Let's now run a series of actors and play around with plotting to get a sense of what each individual actor does.

Here's how we can restore things back to after the initialization stage (in case we did anything else in between)

@checkout :init dd ini act

Actors in FUSE can be executed by passing two arguments to them: dd and act.

Let's start by positioning the PF coils, so that we stand a chance to reproduce the desired plasma shape. This will be important to ensure the stability of the ActorStationaryPlasma that we are going to run next.

FUSE.ActorPFdesign(dd, act; do_plot=true); # instead of setting `act.ActorPFdesign.do_plot=true` we can just pass `do_plot=true` as argument without chaning `act`
actors: PFdesign

The ActorStationaryPlasma iterates between plasma transport, pedestal, equilibrium and sources to return a self-consistent plasma solution

peq = plot(dd.equilibrium; label="before")
pcp = plot(dd.core_profiles; color=:gray, label="before")
#act.ActorFluxMatcher.verbose = true
act.ActorFluxMatcher.algorithm = :anderson
#act.ActorFluxMatcher.step_size = 0.1
FUSE.ActorStationaryPlasma(dd, act);
actors: StationaryPlasma
actors:  --------------- 1/5
actors:  HCD
actors:   SimpleEC
actors:   SimpleIC
actors:   NeutralFueling
actors:  Current
actors:   QED
actors:  Pedestal
actors:   EPED
actors:  CoreTransport
actors:   FluxMatcher
actors:    Pedestal
actors:     EPED
actors:    FluxCalculator
actors:     TGLF
actors:     Neoclassical
actors:  Current
actors:   QED
actors:  Equilibrium
actors:   TEQUILA
actors:  --------------- 1/5 @ 674.61%
actors:  HCD
actors:   SimpleEC
actors:   SimpleIC
actors:   NeutralFueling
actors:  Current
actors:   QED
actors:  Pedestal
actors:   EPED
actors:  CoreTransport
actors:   FluxMatcher
actors:    Pedestal
actors:     EPED
actors:    FluxCalculator
actors:     TGLF
actors:     Neoclassical
actors:  Current
actors:   QED
actors:  Equilibrium
actors:   TEQUILA
actors:  --------------- 2/5 @ 645.88%
actors:  HCD
actors:   SimpleEC
actors:   SimpleIC
actors:   NeutralFueling
actors:  Current
actors:   QED
actors:  Pedestal
actors:   EPED
actors:  CoreTransport
actors:   FluxMatcher
actors:    Pedestal
actors:     EPED
actors:    FluxCalculator
actors:     TGLF
actors:     Neoclassical
actors:  Current
actors:   QED
actors:  Equilibrium
actors:   TEQUILA
actors:  --------------- 3/5 @ 290.58%
actors:  HCD
actors:   SimpleEC
actors:   SimpleIC
actors:   NeutralFueling
actors:  Current
actors:   QED
actors:  Pedestal
actors:   EPED
actors:  CoreTransport
actors:   FluxMatcher
actors:    Pedestal
actors:     EPED
actors:    FluxCalculator
actors:     TGLF
actors:     Neoclassical
actors:  Current
actors:   QED
actors:  Equilibrium
actors:   TEQUILA
actors:  --------------- 4/5 @ 169.42%
actors:  HCD
actors:   SimpleEC
actors:   SimpleIC
actors:   NeutralFueling
actors:  Current
actors:   QED
actors:  Pedestal
actors:   EPED
actors:  CoreTransport
actors:   FluxMatcher
actors:    Pedestal
actors:     EPED
actors:    FluxCalculator
actors:     TGLF
actors:     Neoclassical
actors:  Current
actors:   QED
actors:  Equilibrium
actors:   TEQUILA
actors:  --------------- 5/5 @ 104.03%
┌ Warning: Max number of iterations (5) has been reached with convergence error of (1)[0.337, 0.323, 0.145, 0.085, 0.052](5) compared to threshold of 0.05
└ @ FUSE ~/work/FUSE.jl/FUSE.jl/src/actors/compound/stationary_plasma_actor.jl:209

we can compare equilibrium before and after the self-consistency loop

plot!(peq, dd.equilibrium; label="after")
Example block output

we can compare core_profiles before and after the self-consistency loop

plot!(pcp, dd.core_profiles; label="after")
Example block output

here are the sources

plot(dd.core_sources)
Example block output

and the flux-matched transport

plot(dd.core_transport)

HFS sizing actor changes the thickness of the OH and TF layers on the high field side to satisfy current and stresses constraints

plot(dd.build)
FUSE.ActorHFSsizing(dd, act);
plot!(dd.build; cx=false)
Example block output

The stresses on the center stack are stored in the solid_mechanics IDS

plot(dd.solid_mechanics.center_stack.stress)
Example block output

LFS sizing actors change location of the outer TF leg to meet ripple requirements

plot(dd.build)
FUSE.ActorLFSsizing(dd, act);
plot!(dd.build; cx=false)
Example block output

A custom show() method is defined to print the summary of dd.build.layer

dd.build.layer
23×10 DataFrame
 Row │ group   details                            type      ΔR          R_start   R_end     material      area        volume     shape
     │ String  String                             String    Float64     Float64   Float64   String        Float64     Float64    String
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ in                                                   1.24133      0.0       1.24133  steel          20.0726      98.9942  rectangle
   2 │ in                                         oh        0.343278     1.24133   1.58461  nb3sn           6.9632      80.5952  rectangle
   3 │ hfs                                        tf        1.71333      1.58461   3.29794  nb3sn_kdemo    42.9246     889.465   convex hull
   4 │ hfs     gap tf vacuum vessel                         0.0          3.29794   3.29794  vacuum          7.46559    490.871   double ellipse
   5 │ hfs     vacuum  outer                      vessel    0.0982993    3.29794   3.39624  steel           3.06115    124.717   negative offset
   6 │ hfs     gap water                                    0.147449     3.39624   3.54369  water           4.47786    182.552   negative offset
   7 │ hfs     vacuum  inner                      vessel    0.0982993    3.54369   3.64199  steel           2.90933    118.686   negative offset
   8 │ hfs     gap high temp shield vacuum vess…            0.00982993   3.64199   3.65182  vacuum          0.287593    11.7359  negative offset
   9 │ hfs     high temp                          shield    0.196599     3.65182   3.84842  steel           5.62432    229.652   negative offset
  10 │ hfs                                        blanket   0.373538     3.84842   4.22196  lithium_lead   20.8557     951.46    negative offset
  11 │ hfs     first                              wall      0.0196599    4.22196   4.24162  tungsten        1.02635     34.1644  offset
  12 │ lhfs                                       plasma    4.51651      4.24162   8.75812  plasma         32.1335    1252.4
  13 │ lfs     first                              wall      0.0196599    8.75812   8.77778  tungsten        1.02635     34.1644  offset
  14 │ lfs                                        blanket   1.17959      8.77778   9.95738  lithium_lead   20.8557     951.46    negative offset
  15 │ lfs     high temp                          shield    0.196599     9.95738  10.154    steel           5.62432    229.652   negative offset
  16 │ lfs     gap high temp shield vacuum vess…            0.0442026   10.154    10.1982   vacuum          0.287593    11.7359  negative offset
  17 │ lfs     vacuum  inner                      vessel    0.0982993   10.1982   10.2965   steel           2.90933    118.686   negative offset
  18 │ lfs     gap water                                    0.147449    10.2965   10.4439   water           4.47786    182.552   negative offset
  19 │ lfs     vacuum  outer                      vessel    0.0982993   10.4439   10.5422   steel           3.06115    124.717   negative offset
  20 │ lfs     gap tf vacuum vessel                         0.775582    10.5422   11.3178   vacuum          7.46559    490.871   double ellipse
  21 │ lfs                                        tf        1.71333     11.3178   13.0311   nb3sn_kdemo    42.9246     889.465   convex hull
  22 │ out                                                  1.96599      0.0      14.9971   vacuum        164.187     8377.39
  23 │ out                                        cryostat  0.0982993    0.0      15.0954   steel           4.86738    309.732   silo

ActorHFSsizing and ActorLFSsizing only change the layer's thicknesses, so we then need to trigger a build of the 2D cross-sections after them:

FUSE.ActorCXbuild(dd, act);
plot(dd.build)
Example block output

Generate passive structures information (for now the vacuum vessel)

FUSE.ActorPassiveStructures(dd, act)
plot(dd.pf_passive)
Example block output

We can now give the PF coils their final position given the new build

actor = FUSE.ActorPFdesign(dd, act);
plot(actor) # some actors define their own plot
Example block output

With information about both pfactive and pfpassive we can now evaluate vertical stability

FUSE.ActorVerticalStability(dd, act)
IMAS.freeze(dd.mhd_linear)
mhd_linear [DETACHED]
├─ time[0] [s]
└─ time_slice
   └─ 1
      ├─ time ➡ 0 [s]
      └─ toroidal_mode
         ├─ 1
         │  ├─ n_tor0
         │  ├─ perturbation_type
         │  │  ├─ description"Vertical stability margin > 0.15 for stability"
         │  │  └─ name"m_s"
         │  └─ stability_metric ➡ 0.170536
         └─ 2
            ├─ n_tor0
            ├─ perturbation_type
            │  ├─ description"Normalized vertical growth rate < 10 for stability"
            │  └─ name"γτ"
            └─ stability_metric ➡ 7.06988

The ActorNeutronics calculates the heat flux on the first wall

FUSE.ActorNeutronics(dd, act);
p = plot(; layout=2, size=(900, 350))
plot!(p, dd.neutronics.time_slice[].wall_loading, subplot=1)
plot!(p, FUSE.define_neutrons(dd, 100000)[1], dd.equilibrium.time_slice[]; subplot=1, colorbar_entry=false)
plot!(p, dd.neutronics.time_slice[].wall_loading; cx=false, subplot=2, ylabel="")
Example block output

The ActorBlanket will change the thickess of the first wall, breeder, shield, and Li6 enrichment to achieve target TBR

FUSE.ActorBlanket(dd, act);
print_tree(IMAS.freeze(dd.blanket); maxdepth=5)
actors: Blanket
blanket [DETACHED]
├─ module
│  └─ 1
│     ├─ layer
│     │  ├─ 1
│     │  │  ├─ material ➡ "tungsten"
│     │  │  ├─ midplane_thickness ➡ 0.0199587 [m]
│     │  │  └─ name ➡ "lfs first wall"
│     │  ├─ 2
│     │  │  ├─ material ➡ "lithium-lead: Li6/7=90.000%"
│     │  │  ├─ midplane_thickness ➡ 1.29001 [m]
│     │  │  └─ name ➡ "lfs blanket"
│     │  └─ 3
│     │     ├─ material ➡ "steel"
│     │     ├─ midplane_thickness ➡ 0.0858821 [m]
│     │     └─ name ➡ "lfs high temp shield"
│     ├─ name ➡ "blanket"
│     └─ time_slice
│        └─ 1
│           ├─ peak_escape_flux ➡ 181061 [W/m^2]
│           ├─ peak_wall_flux ➡ 1.04024e+06 [W/m^2]
│           ├─ power_incident_neutrons ➡ 9.95942e+06 [W]
│           ├─ power_incident_radiated ➡ 0 [W]
│           ├─ power_thermal_extracted ➡ 1.19513e+07 [W]
│           ├─ power_thermal_neutrons ➡ 1.19513e+07 [W]
│           ├─ power_thermal_radiated ➡ 0 [W]
│           ├─ time ➡ 0 [s]
│           └─ tritium_breeding_ratio ➡ 1.68669
├─ time ➡ [0] [s]
└─ tritium_breeding_ratio ➡ [0.0681484]

The ActorDivertors actor calculates the divertors heat flux

FUSE.ActorDivertors(dd, act);
print_tree(IMAS.freeze(dd.divertors); maxdepth=4)
actors: Divertors
divertors [DETACHED]
├─ divertor
│  └─ 1
│     ├─ power_black_body
│     │  ├─ data ➡ [0] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_conducted
│     │  ├─ data ➡ [1.24193e+08] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_convected
│     │  ├─ data ➡ [0] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_currents
│     │  ├─ data ➡ [0] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_incident
│     │  ├─ data ➡ [3.53054e+07] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_neutrals
│     │  ├─ data ➡ [0] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_radiated
│     │  ├─ data ➡ [0] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_recombination_neutrals
│     │  ├─ data ➡ [0] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_recombination_plasma
│     │  ├─ data ➡ [0] [W]
│     │  └─ time ➡ [0] [s]
│     ├─ power_thermal_extracted
│     │  ├─ data ➡ [3.53054e+07] [W]
│     │  └─ time ➡ [0] [s]
│     └─ target
│        ├─ 1
│        │  ⋮
│        │
│        └─ 2
│           ⋮
│
└─ time ➡ [0] [s]

The ActorBalanceOfPlant calculates the optimal cooling flow rates for the heat sources (breeder, divertor, and wall) and get an efficiency for the electricity conversion cycle

FUSE.ActorBalanceOfPlant(dd, act);
IMAS.freeze(dd.balance_of_plant)
balance_of_plant [DETACHED]
├─ Q_plant[0.133298]
├─ power_electric_net[-1.25672e+08] [W]
├─ power_electric_plant_operation
│  ├─ system
│  │  ├─ 1
│  │  │  ├─ index1
│  │  │  ├─ name"HCD"
│  │  │  ├─ power[1e+08] [W]
│  │  │  └─ subsystem
│  │  │     ├─ 1
│  │  │     │  ├─ index1
│  │  │     │  ├─ name"nbi"
│  │  │     │  └─ power[0] [W]
│  │  │     ├─ 2
│  │  │     │  ├─ index2
│  │  │     │  ├─ name"ec_launchers"
│  │  │     │  └─ power[5e+07] [W]
│  │  │     ├─ 3
│  │  │     │  ├─ index3
│  │  │     │  ├─ name"ic_antennas"
│  │  │     │  └─ power[5e+07] [W]
│  │  │     └─ 4
│  │  │        ├─ index4
│  │  │        ├─ name"lh_antennas"
│  │  │        └─ power[0] [W]
│  │  ├─ 2
│  │  │  ├─ index3
│  │  │  ├─ name"cryostat"
│  │  │  └─ power[3e+07] [W]
│  │  ├─ 3
│  │  │  ├─ index4
│  │  │  ├─ name"tritium_handling"
│  │  │  └─ power[1.5e+07] [W]
│  │  └─ 4
│  │     ├─ index6
│  │     ├─ name"pf_active"
│  │     └─ power[0] [W]
│  └─ total_power[1.45e+08] [W]
├─ power_plant
│  ├─ heat_load
│  │  ├─ breeder[1.19513e+07] [W]
│  │  ├─ divertor[3.53054e+07] [W]
│  │  └─ wall[3.78921e+07] [W]
│  ├─ power_cycle_type"rankine"
│  ├─ power_electric_generated[1.93282e+07] [W]
│  └─ total_heat_supplied[8.51488e+07] [W]
├─ thermal_efficiency_plant[0.226993]
└─ time[0] [s]

ActorCosting will break down the capital and operational costs

FUSE.ActorCosting(dd, act)
plot(dd.costing)
Example block output

Let's checkpoint our results

@checkin :manual dd ini act

Whole facility design

Here we restore the :init checkpoint that we had previously stored. Resetting any changes to dd, ini, and act that we did in the meantime.

@checkout :init dd ini act

Actors can call other actors, creating workflows. For example, the ActorWholeFacility can be used to to get a self-consistent stationary whole facility design.

FUSE.ActorWholeFacility(dd, act);
actors: WholeFacility
actors:  PFdesign
actors:  StationaryPlasma
actors:   --------------- 1/5
actors:   HCD
actors:    SimpleEC
actors:    SimpleIC
actors:    NeutralFueling
actors:   Current
actors:    QED
actors:   Pedestal
actors:    EPED
actors:   CoreTransport
actors:    FluxMatcher
actors:     Pedestal
actors:      EPED
actors:     FluxCalculator
actors:      TGLF
actors:      Neoclassical
actors:   Current
actors:    QED
actors:   Equilibrium
actors:    TEQUILA
actors:   --------------- 1/5 @ 670.38%
actors:   HCD
actors:    SimpleEC
actors:    SimpleIC
actors:    NeutralFueling
actors:   Current
actors:    QED
actors:   Pedestal
actors:    EPED
actors:   CoreTransport
actors:    FluxMatcher
actors:     Pedestal
actors:      EPED
actors:     FluxCalculator
actors:      TGLF
actors:      Neoclassical
actors:   Current
actors:    QED
actors:   Equilibrium
actors:    TEQUILA
actors:   --------------- 2/5 @ 627.55%
actors:   HCD
actors:    SimpleEC
actors:    SimpleIC
actors:    NeutralFueling
actors:   Current
actors:    QED
actors:   Pedestal
actors:    EPED
actors:   CoreTransport
actors:    FluxMatcher
actors:     Pedestal
actors:      EPED
actors:     FluxCalculator
actors:      TGLF
actors:      Neoclassical
actors:   Current
actors:    QED
actors:   Equilibrium
actors:    TEQUILA
actors:   --------------- 3/5 @ 309.99%
actors:   HCD
actors:    SimpleEC
actors:    SimpleIC
actors:    NeutralFueling
actors:   Current
actors:    QED
actors:   Pedestal
actors:    EPED
actors:   CoreTransport
actors:    FluxMatcher
actors:     Pedestal
actors:      EPED
actors:     FluxCalculator
actors:      TGLF
actors:      Neoclassical
actors:   Current
actors:    QED
actors:   Equilibrium
actors:    TEQUILA
actors:   --------------- 4/5 @ 224.43%
actors:   HCD
actors:    SimpleEC
actors:    SimpleIC
actors:    NeutralFueling
actors:   Current
actors:    QED
actors:   Pedestal
actors:    EPED
actors:   CoreTransport
actors:    FluxMatcher
actors:     Pedestal
actors:      EPED
actors:     FluxCalculator
actors:      TGLF
actors:      Neoclassical
actors:   Current
actors:    QED
actors:   Equilibrium
actors:    TEQUILA
actors:   --------------- 5/5 @ 145.64%
┌ Warning: Max number of iterations (5) has been reached with convergence error of (1)[0.335, 0.314, 0.155, 0.112, 0.073](5) compared to threshold of 0.05
└ @ FUSE ~/work/FUSE.jl/FUSE.jl/src/actors/compound/stationary_plasma_actor.jl:209
actors:  HFSsizing
actors:   FluxSwing
actors:   Stresses
actors:  LFSsizing
actors:  CXbuild
actors:  PFdesign
actors:  Equilibrium
actors:   TEQUILA
actors:  CXbuild
actors:  Neutronics
actors:  Blanket
actors:  CXbuild
actors:  PassiveStructures
actors:  Divertors
actors:  PlasmaLimits
actors:   VerticalStability
actors:   TroyonBetaNN
actors:  BalanceOfPlant
actors:   ThermalPlant
actors:   PowerNeeds
actors:  Costing
actors:   CostingARIES

Let's check what we got at a glance with the FUSE.digest(dd) function:

FUSE.digest(dd)
GEOMETRY                                EQUILIBRIUM                             TEMPERATURES
────────────────────────────────────    ────────────────────────────────────    ────────────────────────────────────
R0 → 6.5 [m]                            B0 → 7.8 [T]                            Te0 → 18.6 [keV]
a → 2.01 [m]                            ip → 12.7 [MA]                          Ti0 → 18.1 [keV]
1/ϵ → 3.24                              q95 → 6.12                              <Te> → 8.68 [keV]
κ → 2                                   <Bpol> → 0.802 [T]                      <Ti> → 7.88 [keV]
δ → 0.588                               βpol_MHD → 0.773                        Te0/<Te> → 2.14
ζ → -0.0136                             βtor_MHD → 0.00844                      Ti0/<Ti> → 2.3
Volume → 931 [m³]                       βn_MHD → 1.03
Surface → 759 [m²]

DENSITIES                               PRESSURES                               TRANSPORT
────────────────────────────────────    ────────────────────────────────────    ────────────────────────────────────
ne0 → 9.04e+19 [m⁻³]                    P0 → 0.503 [MPa]                        τe → 2.37 [s]
ne_ped → 6.74e+19 [m⁻³]                 <P> → 0.205 [MPa]                       τe_exp → 1.84 [s]
ne_line → 8.26e+19 [m⁻³]                P0/<P> → 2.46                           H98y2 → 0.864
<ne> → 7.56e+19 [m⁻³]                   βn → 1.04                               H98y2_exp → 0.798
ne0/<ne> → 1.2                          βn_th → 1.04                            Hds03 → 0.643
fGW → 0.821                                                                     Hds03_exp → 0.573
zeff_ped → 2                                                                    τα_thermalization → 0.892 [s]
<zeff> → 2                                                                      τα_slowing_down → 0.99 [s]
impurities → DT Ne20 He4

SOURCES                                 EXHAUST                                 CURRENTS
────────────────────────────────────    ────────────────────────────────────    ────────────────────────────────────
Pec → 50 [MW]                           Psol → 120 [MW]                         ip_bs_aux_ohm → 13 [MA]
rho0_ec → 0.56 [MW]                     PLH → 134 [MW]                          ip_ni → 6.54 [MA]
Pnbi → NaN [MW]                         Bpol_omp → 1.11 [T]                     ip_bs → 3.06 [MA]
Enbi1 → NaN [MeV]                       λq → 0.978 [mm]                         ip_aux → 3.48 [MA]
Pic → 50 [MW]                           qpol → 2.3e+03 [MW/m²]                  ip_ohm → 6.44 [MA]
Plh → NaN [MW]                          qpar → 1.26e+04 [MW/m²]                 ejima → 0.4
Paux_tot → 100 [MW]                     P/R0 → 18.5 [MW/m]                      flattop → 0.7 [Hours]
Pα → 54.8 [MW]                          PB/R0 → 144 [MW T/m]
Pohm → 0.523 [MW]                       PBp/R0 → 14.8 [MW T/m]
Pheat → 155 [MW]                        PBϵ/R0q95 → 7.27 [MW T/m]
Prad_tot → -35.1 [MW]                   neutrons_peak → 0.362 [MW/m²]

BOP                                     BUILD                                   COSTING
────────────────────────────────────    ────────────────────────────────────    ────────────────────────────────────
Pfusion → 274 [MW]                      PF_material → nb3sn                     capital_cost → 6.75 [$B]
Qfusion → 2.74                          TF_material → nb3sn_kdemo               levelized_CoE → Inf [$/kWh]
thermal_cycle_type → rankine            OH_material → nb3sn                     TF_of_total → 17.4 [%]
thermal_efficiency_plant → 22.7 [%]     TF_max_b → 15.4 [T]                     BOP_of_total → 1.78 [%]
thermal_efficiency_cycle → NaN [%]      OH_max_b → 16 [T]                       blanket_of_total → 20.8 [%]
power_electric_generated → 18.2 [MW]    TF_j_margin → 6.68                      cryostat_of_total → 3.01 [%]
Pelectric_net → -127 [MW]               OH_j_margin → 1.4
Qplant → 0.126                          TF_stress_margin → 3.02
TBR → 0.0732                            OH_stress_margin → 1.2

@ time = 0.0 [s]
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
GKS: could not find font middle.ttf
​
​
​

Like before we can checkpoint results for later use

@checkin :awf dd ini act

Getting into the weeds

Saving and loading data

tutorial_temp_dir = tempdir()
filename = joinpath(tutorial_temp_dir, "$(ini.general.casename).json")
"/tmp/K-DEMO.json"

When saving data to be shared outside of FUSE, one can set freeze=true so that all expressions in the dd are evaluated and saved to file.

IMAS.imas2json(dd, filename; freeze=false, strict=false);

Load from JSON

dd1 = IMAS.json2imas(filename);

Exploring the data dictionary

  • FUSE stores data following the IMAS data schema.
  • The root of the data structure is dd, which stands for "Data Dictionary".
  • More details are available in the documentation.

Display part of the equilibrium data in dd

dd.equilibrium.time_slice[1].boundary
boundary
├─ elongation ➡ 1.999
├─ elongation_lowerFunction
├─ elongation_upperFunction
├─ geometric_axis
│  ├─ r ➡ 6.49987 [m]
│  └─ z ➡ 0.0124918 [m]
├─ minor_radius ➡ 2.00734 [m]
├─ outline
│  ├─ r243-element Vector{Float64} [m]
│  │      min:4.51   avg:6.25   max:8.49
│  └─ z243-element Vector{Float64} [m]min:-4.17   avg:-0.0578   max:3.92
├─ ovality ➡ 0.00572603
├─ psi ➡ -0.786039 [Wb]
├─ squareness ➡ -0.0135544
├─ squareness_lower_innerFunction
├─ squareness_lower_outerFunction
├─ squareness_upper_innerFunction
├─ squareness_upper_outerFunction
├─ strike_point
│  ├─ 1
│  │  ├─ last_closed_flux_surface_gap ➡ 1.19218e-08 [m]
│  │  ├─ r ➡ 4.22092 [m]
│  │  └─ z ➡ -4.50703 [m]
│  └─ 2
│     ├─ last_closed_flux_surface_gap ➡ 1.19218e-08 [m]
│     ├─ r ➡ 5.53843 [m]
│     └─ z ➡ -5.16481 [m]
├─ tilt ➡ -0.00162958
├─ triangularity ➡ 0.587596
├─ triangularity_lowerFunction
├─ triangularity_upperFunction
├─ twist ➡ 0.00165424
└─ x_point
   ├─ 1
   │  ├─ r ➡ 5.0906 [m]
   │  └─ z ➡ -4.15032 [m]
   └─ 2
      ├─ r ➡ 4.86063 [m]
      └─ z ➡ 4.37322 [m]

this can be done up to a certain depth with print_tree

print_tree(dd.equilibrium.time_slice[1].boundary; maxdepth=1)
boundary
├─ elongation ➡ 1.999
├─ elongation_lower ➡ Function
├─ elongation_upper ➡ Function
├─ geometric_axis
│  ⋮
│
├─ minor_radius ➡ 2.00734 [m]
├─ outline
│  ⋮
│
├─ ovality ➡ 0.00572603
├─ psi ➡ -0.786039 [Wb]
├─ squareness ➡ -0.0135544
├─ squareness_lower_inner ➡ Function
├─ squareness_lower_outer ➡ Function
├─ squareness_upper_inner ➡ Function
├─ squareness_upper_outer ➡ Function
├─ strike_point
│  ⋮
│
├─ tilt ➡ -0.00162958
├─ triangularity ➡ 0.587596
├─ triangularity_lower ➡ Function
├─ triangularity_upper ➡ Function
├─ twist ➡ 0.00165424
└─ x_point
   ⋮

Plotting data from dd

FUSE uses Plots.jl recipes for visualizing data from dd.

This allows different plots to be shown when calling plot() on different items in the data structure.

Learn more about Plots.jl here

For example plotting the equilibrium...

plot(dd.equilibrium)
Example block output

...or the core profiles

plot(dd.core_profiles)
Example block output

Whant to know what arguments can be passed? use help_plot() function

help_plot(dd.equilibrium; core_profiles_overlay=true, levels_in=21, levels_out=5, show_secondary_separatrix=true, coordinate=:rho_tor_norm)
Example block output

These plots can be composed by calling plot!() instead of plot()

plot(dd.equilibrium; color=:gray, cx=true)
plot!(dd.build.layer)
plot!(dd.pf_active)
plot!(dd.pf_passive)
plot!(dd.pulse_schedule.position_control; color=:red)
Example block output

Plotting an array...

plot(dd.core_profiles.profiles_1d[1].pressure_thermal)
Example block output

...is different from plotting a field from the IDS (which plots the quantity against its coordinate and with units)

plot(dd.core_profiles.profiles_1d[1], :pressure_thermal)
Example block output

Customizing plot attributes:

plot(dd.core_profiles.profiles_1d[1], :pressure_thermal; label="", linewidth=2, color=:red, labelfontsize=25)
Example block output

Use findall(ids, r"...") to search for certain fields. In Julia, string starting with r are regular expressions.

findall(dd, r"\.psi")
[1] dd.core_profiles.profiles_1d[1].grid.psi_norm [101-element Vector{Float64}] (min:0, avg:0.465, max:1)
[2] dd.core_profiles.profiles_1d[1].grid.psi [Wb] [101-element Vector{Float64}] (min:-69.2, avg:-37.4, max:-0.887)
[3] dd.equilibrium.time_slice[1].boundary.psi [Wb] [Float64] (all:-0.786)
[4] dd.equilibrium.time_slice[1].boundary_separatrix.psi [Wb] [Float64] (all:-0.786)
[5] dd.equilibrium.time_slice[1].global_quantities.psi_boundary [Wb] [Float64] (all:-0.786)
[6] dd.equilibrium.time_slice[1].global_quantities.psi_axis [Wb] [Float64] (all:-73.9)
[7] dd.equilibrium.time_slice[1].profiles_1d.psi_norm [129-element Vector{Float64}] (min:0, avg:0.5, max:1)
[8] dd.equilibrium.time_slice[1].profiles_1d.psi [Wb] [129-element Vector{Float64}] (min:-73.9, avg:-37.3, max:-0.786)
[9] dd.equilibrium.time_slice[1].profiles_2d[1].psi [Wb] [31×13 Matrix{Float64}] (min:-0.038, avg:0.634, max:6.65)
[10] dd.equilibrium.time_slice[1].profiles_2d[2].psi [Wb] [66×129 Matrix{Float64}] (min:-73.9, avg:17.8, max:183)

findall(ids, r"...") can be combined with plot() to plot multiple fields

plot(findall(dd, r"\.psi"))
Example block output

Working with time series

The IMAS data structure supports time-dependent data, and IMAS.jl provides ways to handle time data efficiently.

Each dd has a global_time attribute, which is used throughout FUSE and IMAS to indicate the time at which things should be operate.

dd.global_time
0.0

For the sake of demonstrating handling of time, let's add a new time_slice to the equilibrium.

NOTE: in addition to the usual resize!(ids, n::Int), time dependent arrays of structures can be resized with:

  • resize!(ids) which will add a time-slice at the current global time (this is what you want to use in most cases)
  • resize!(ids, time0) which will add a time-slice at time0 seconds

resize the time dependent array of structure

resize!(dd.equilibrium.time_slice, 1.0);

let's just populate it with the data from the previous time slice

dd.equilibrium.time_slice[2] = deepcopy(dd.equilibrium.time_slice[1]);
dd.equilibrium.time_slice[2].time = 1.0
1.0

Here we see that equilibrium has mulitiple time_slices

dd.equilibrium.time
2-element Vector{Float64}:
 0.0
 1.0

We can access time-dependent arrays of structures via integer index...

eqt = dd.equilibrium.time_slice[2]
eqt.time
1.0

...or at a given time, by passing the time as a floating point number (in seconds)

eqt = dd.equilibrium.time_slice[1.0]
eqt.time
1.0

NOTE: If we ask a time that is not exactly in the arrays of structures, we'll get the closest (causal!) time-slice

eqt = dd.equilibrium.time_slice[0.9]
eqt.time

eqt = dd.equilibrium.time_slice[1.1]
eqt.time
1.0

... or at the current dd.global_time by leaving the square brackets empty []

NOTE: using [] is what you want to use in most situations that involve time-dependent arrays of structures!

dd.global_time = 0.0
eqt = dd.equilibrium.time_slice[]
eqt.time

dd.global_time = 1.0
eqt = dd.equilibrium.time_slice[]
eqt.time
1.0

What we described above was for time-dependent arrays of structures.

The other place where time comes in, is when dealing with time-dependent arrays of data.

In this case, we can use the @ddtime macro to manipulate these time-dependent arrays at dd.global_time.

NOTE: Also in this case, @ddtime will operate on the closest (causal!) time point

dd.equilibrium.vacuum_toroidal_field.b0

dd.global_time = 1.0
@ddtime(dd.equilibrium.vacuum_toroidal_field.b0 = 10.0)
dd.equilibrium.vacuum_toroidal_field.b0

dd.global_time = 0.0
@ddtime(dd.equilibrium.vacuum_toroidal_field.b0)

dd.global_time = 1.0
@ddtime(dd.equilibrium.vacuum_toroidal_field.b0)
10.0

Expressions in dd

Some fields in the data dictionary are expressions (ie. Functions). For example dd.core_profiles.profiles_1d[].pressure is dynamically calculated as the product of thermal densities and temperature with addition of fast ions contributions

dd.global_time = 0.0
print_tree(dd.core_profiles.profiles_1d[]; maxdepth=1)
1
├─ conductivity_parallel ➡ Function [ohm^-1.m^-1]
├─ electrons
│  ⋮
│
├─ grid
│  ⋮
│
├─ ion
│  ⋮
│
├─ j_bootstrap ➡ 101-element Vector{Float64} [A/m^2]
│                min:374   avg:1.12e+05   max:5.64e+05
├─ j_non_inductive ➡ 101-element Vector{Float64} [A/m^2]
│                    min:2.4e+03   avg:4.86e+05   max:1.31e+06
├─ j_ohmic ➡ 101-element Vector{Float64} [A/m^2]
│            min:984   avg:3.86e+05   max:6.64e+05
├─ j_tor ➡ 101-element Vector{Float64} [A/m^2]
│          min:8.21e+03   avg:8.66e+05   max:1.98e+06
├─ j_total ➡ 101-element Vector{Float64} [A/m^2]
│            min:9.14e+03   avg:8.73e+05   max:1.96e+06
├─ neutral
│  ⋮
│
├─ pressure ➡ Function [Pa]
├─ pressure_ion_total ➡ Function [Pa]
├─ pressure_parallel ➡ Function [Pa]
├─ pressure_perpendicular ➡ Function [Pa]
├─ pressure_thermal ➡ Function [Pa]
├─ rotation_frequency_tor_sonic ➡ 101-element Vector{Float64} [s^-1]
│                                 all:0
├─ t_i_average ➡ Function [eV]
├─ time ➡ 0 [s]
└─ zeff ➡ 101-element Vector{Float64}
          all:2

accessing a dynamic expression, automatically evaluates it

dd.core_profiles.profiles_1d[].conductivity_parallel
101-element Vector{Float64}:
      1.5971093395933754e9
      1.5545123478521702e9
      1.4972681669431643e9
      1.4377509348652933e9
      1.3774330741251698e9
      1.3167097988600204e9
      1.2672655360694363e9
      1.224180166032008e9
      1.184929712410916e9
      1.1506355756390693e9
      ⋮
      9.489287543462877e7
      7.883895371486433e7
      6.2277285342068225e7
      4.6530348321967095e7
      3.26521948055096e7
      2.1242908687460877e7
      1.2338318542953225e7
      5.568499860385064e6
 661466.7670929075

In addition to evaluating expressions by accessing them, expressions in the tree can be evaluated using IMAS.freeze(ids)

NOTE: IMAS.freeze(ids, field::Symbol) works on a single field and IMAS.refreeze!(ids, field) forces re-evaluation of an expression. Also, IMAS.empty!(ids, field::Symbol) can be used to revert a frozen field back into an expression.

print_tree(IMAS.freeze(dd.core_profiles.profiles_1d[1]); maxdepth=1)
profiles_1d [DETACHED]
├─ conductivity_parallel ➡ 101-element Vector{Float64} [ohm^-1.m^-1]
│                          min:6.61e+05   avg:4.71e+08   max:1.6e+09
├─ electrons
│  ⋮
│
├─ grid
│  ⋮
│
├─ ion
│  ⋮
│
├─ j_bootstrap ➡ 101-element Vector{Float64} [A/m^2]
│                min:374   avg:1.12e+05   max:5.64e+05
├─ j_non_inductive ➡ 101-element Vector{Float64} [A/m^2]
│                    min:2.4e+03   avg:4.86e+05   max:1.31e+06
├─ j_ohmic ➡ 101-element Vector{Float64} [A/m^2]
│            min:984   avg:3.86e+05   max:6.64e+05
├─ j_tor ➡ 101-element Vector{Float64} [A/m^2]
│          min:8.21e+03   avg:8.66e+05   max:1.98e+06
├─ j_total ➡ 101-element Vector{Float64} [A/m^2]
│            min:9.14e+03   avg:8.73e+05   max:1.96e+06
├─ neutral
│  ⋮
│
├─ pressure ➡ 101-element Vector{Float64} [Pa]
│             min:432   avg:2.71e+05   max:5.03e+05
├─ pressure_ion_total ➡ 101-element Vector{Float64} [Pa]
│                       min:206   avg:1.22e+05   max:2.34e+05
├─ pressure_parallel ➡ 101-element Vector{Float64} [Pa]
│                      min:144   avg:9.02e+04   max:1.68e+05
├─ pressure_perpendicular ➡ 101-element Vector{Float64} [Pa]
│                           min:144   avg:9.02e+04   max:1.68e+05
├─ pressure_thermal ➡ 101-element Vector{Float64} [Pa]
│                     min:432   avg:2.71e+05   max:5.03e+05
├─ rotation_frequency_tor_sonic ➡ 101-element Vector{Float64} [s^-1]
│                                 all:0
├─ t_i_average ➡ 101-element Vector{Float64} [eV]
│                min:85.3   avg:1.01e+04   max:1.81e+04
├─ time ➡ 0 [s]
└─ zeff ➡ 101-element Vector{Float64}
          all:2

Comparing two IDSs

We can introduce a change in the dd1 and spot it with the diff function

dd1.equilibrium.time_slice[1].time = -100.0
IMAS.diff(dd.equilibrium, dd1.equilibrium)
Dict{String, String} with 3 entries:
  "time"                     => "length:  2 --  1"
  "vacuum_toroidal_field.b0" => "length:  2 --  1"
  "time_slice"               => "length:  2 --  1"

Summary

Snapshot of dd in 0D quantities (evaluated at dd.global_time).

Extract + plots saved to PDF (printed to screen if filename is omitted). NOTE: For PDF creation to work, one may need to install of DejaVu Sans Mono font.

filename = joinpath(tutorial_temp_dir, "$(ini.general.casename).pdf")
display(filename)
FUSE.digest(dd)#, filename)
GEOMETRY                                EQUILIBRIUM                             TEMPERATURES
────────────────────────────────────    ────────────────────────────────────    ────────────────────────────────────
R0 → 6.5 [m]                            B0 → 7.8 [T]                            Te0 → 18.6 [keV]
a → 2.01 [m]                            ip → 12.7 [MA]                          Ti0 → 18.1 [keV]
1/ϵ → 3.24                              q95 → 6.12                              <Te> → 8.68 [keV]
κ → 2                                   <Bpol> → 0.802 [T]                      <Ti> → 7.88 [keV]
δ → 0.588                               βpol_MHD → 0.773                        Te0/<Te> → 2.14
ζ → -0.0136                             βtor_MHD → 0.00844                      Ti0/<Ti> → 2.3
Volume → 931 [m³]                       βn_MHD → 1.03
Surface → 759 [m²]

DENSITIES                               PRESSURES                               TRANSPORT
────────────────────────────────────    ────────────────────────────────────    ────────────────────────────────────
ne0 → 9.04e+19 [m⁻³]                    P0 → 0.503 [MPa]                        τe → 2.37 [s]
ne_ped → 6.74e+19 [m⁻³]                 <P> → 0.205 [MPa]                       τe_exp → 1.84 [s]
ne_line → 8.26e+19 [m⁻³]                P0/<P> → 2.46                           H98y2 → 0.864
<ne> → 7.56e+19 [m⁻³]                   βn → 1.04                               H98y2_exp → 0.798
ne0/<ne> → 1.2                          βn_th → 1.04                            Hds03 → 0.643
fGW → 0.821                                                                     Hds03_exp → 0.573
zeff_ped → 2                                                                    τα_thermalization → 0.892 [s]
<zeff> → 2                                                                      τα_slowing_down → 0.99 [s]
impurities → DT Ne20 He4

SOURCES                                 EXHAUST                                 CURRENTS
────────────────────────────────────    ────────────────────────────────────    ────────────────────────────────────
Pec → 50 [MW]                           Psol → 120 [MW]                         ip_bs_aux_ohm → 13 [MA]
rho0_ec → 0.56 [MW]                     PLH → 134 [MW]                          ip_ni → 6.54 [MA]
Pnbi → NaN [MW]                         Bpol_omp → 1.11 [T]                     ip_bs → 3.06 [MA]
Enbi1 → NaN [MeV]                       λq → 0.978 [mm]                         ip_aux → 3.48 [MA]
Pic → 50 [MW]                           qpol → 2.3e+03 [MW/m²]                  ip_ohm → 6.44 [MA]
Plh → NaN [MW]                          qpar → 1.26e+04 [MW/m²]                 ejima → 0.4
Paux_tot → 100 [MW]                     P/R0 → 18.5 [MW/m]                      flattop → 0.7 [Hours]
Pα → 54.8 [MW]                          PB/R0 → 144 [MW T/m]
Pohm → 0.523 [MW]                       PBp/R0 → 14.8 [MW T/m]
Pheat → 155 [MW]                        PBϵ/R0q95 → 7.27 [MW T/m]
Prad_tot → -35.1 [MW]                   neutrons_peak → 0.362 [MW/m²]

BOP                                     BUILD                                   COSTING
────────────────────────────────────    ────────────────────────────────────    ────────────────────────────────────
Pfusion → 274 [MW]                      PF_material → nb3sn                     capital_cost → 6.75 [$B]
Qfusion → 2.74                          TF_material → nb3sn_kdemo               levelized_CoE → Inf [$/kWh]
thermal_cycle_type → rankine            OH_material → nb3sn                     TF_of_total → 17.4 [%]
thermal_efficiency_plant → 22.7 [%]     TF_max_b → 15.4 [T]                     BOP_of_total → 1.78 [%]
thermal_efficiency_cycle → NaN [%]      OH_max_b → 16 [T]                       blanket_of_total → 20.8 [%]
power_electric_generated → 18.2 [MW]    TF_j_margin → 6.68                      cryostat_of_total → 3.01 [%]
Pelectric_net → -127 [MW]               OH_j_margin → 1.4
Qplant → 0.126                          TF_stress_margin → 3.02
TBR → 0.0732                            OH_stress_margin → 1.2

@ time = 0.0 [s]
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​