TORAX - TokaMaker Coupling

[nelsonand]

Hi all! (tagging @jcitrin for visibility)

We (with @john4255, @hansec and @CarlosPazSoldan) have been putting together a workflow to couple the Grad-Shafranov solving capabilities from TokaMaker (part of the OpenFUSIONToolkit) with the transport capabilities of TORAX. In sum, since TokaMaker can handle the magnetic evolution of a pulse, including current and voltage trajectories, this lets a user do a full, time-dependent simulation of an entire discharge.

An example of some of the outputs can be seen here for a baseline ITER case:

For the time being, all of this coupling has been done in a separate repo (https://github.com/john4255/simple-pulse-design/tree/main). Feel free to take a look and let us know if you have any suggestions!

I talked about this briefly at the CPS-FR summer school with @tamaranorman, who suggested that it could be nice to have a call sometime to chat about this and potentially what to do with this tool. There are many potential use cases, including potentially eventual integration with MOSAIC (@djbattagl) and as a stand-alone feature. Before that though, there are a couple things that we still have questions on regarding the TORAX implementation, and we should probably find the best way to host, clean and present the code.

Looking forward to the discussion! I'll send an email too in case that's a better form to communicate.

[jcitrin]

Hi @nelsonand ! Thanks for this, extremely interesting. We will respond to the email and set up a chat for the more detailed questions on implementation.

[nelsonand]

Often in our initial tests here we find hollow density or temperature profiles, resulting from either too high of a pedestal specification or two low of a heating/fueling actuator. Is there any way in TORAX to help deal with this? (I.e. have TORAX throw an error with hollow profiles, or automatically try and adjust the edge boundary condition or heating actuator?)

Eg, density profiles for a sim:

[jcitrin]

Would it be possible to take a look at your configs such that we could reproduce this, or otherwise you generate a minimal reproducible example?

Without seeing the config in detail, this looks like a case where there was some initial profile set, and then perturbed very strongly with the pedestal model setting a discontinuous (in time) internal boundary condition at rho_norm ~ 0.9?

If so it's important to either (or both):

• Set a physically consistent initial condition

• Ensure that pedestal formation is carried out over a physically realistic time window instead of instantaneously (by setting ne_ped, Ti_ped, etc with TimeVaryingArrays). Soon we'll have improved an pedestal_model API to more conveniently set a LH transition time window, but for now it must be done by hand. The issue is that if it's instantaneous, you get these very hollow profiles that then lock in since it often stabilizes the turbulence.

I.e. have TORAX throw an error with hollow profiles

This is not recommended since there are physically realistic cases of hollow profiles, e.g. with inner core impurity radiation, or after LH transitions (albeit less extreme hollowing than here).

or automatically try and adjust the edge boundary condition or heating actuator?)

Setting up a mechanism to automatically adjust the fueling source based on a desired line or volume averaged density is something we are considering. But I believe that these issues you are having may stem from setting a "too instantaneous pedestal" so I'd try to fix that first.

Apologies if my guess of what's happening is way off, but it's hard to tell without more info on the config.

[john4255]

Hi all, Thanks for your help. The full config with profiles and everything is attached as a JSON file. It is fairly long so you can also take a look at the code here: https://github.com/john4255/simple-pulse-design/ The config is initialized in baseconfig.py and model.py (function name "_get_torax_config"). The scenario we are looking at is example5.py. If anything here is confusing or seems wrong let me know, the code is still under development. John

…

On Tue, Oct 14, 2025 at 6:47 PM Jonathan Citrin ***@***.***> wrote: Would it be possible to take a look at your configs such that we could reproduce this, or otherwise you generate a minimal reproducible example? Without seeing the config in detail, this looks like a case where there was some initial profile set, and then perturbed very strongly with the pedestal model setting a discontinuous (in time) internal boundary condition at rho_norm ~ 0.9? If so it's important to either (or both): - Set a physically consistent initial condition - Ensure that pedestal formation is carried out over a physically realistic time window instead of instantaneously (by setting ne_ped, Ti_ped, etc with TimeVaryingArrays). Soon we'll have improved an pedestal_model API to more conveniently set a LH transition time window, but for now it must be done by hand. The issue is that if it's instantaneous, you get these very hollow profiles that then lock in since it often stabilizes the turbulence. I.e. have TORAX throw an error with hollow profiles This is not recommended since there are physically realistic cases of hollow profiles, e.g. with inner core impurity radiation, or after LH transitions (albeit less extreme hollowing than here). or automatically try and adjust the edge boundary condition or heating actuator?) Setting up a mechanism to automatically adjust the fueling source based on a desired line or volume averaged density is something we are considering. But I believe that these issues you are having may stem from setting a "too instantaneous pedestal" so I'd try to fix that first. Apologies if my guess of what's happening is way off, but it's hard to tell without more info on the config. — Reply to this email directly, view it on GitHub <#1587 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BPSKCZ4WCJMOVQAV3UG5DAD3XV4O3AVCNFSM6AAAAACHKIF2EKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTINRYGA4TSOI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[jcitrin]

@john4255 , apologies for dropping the ball here.

What is the status here? Is help with your setup and interpretation still needed? Is there any update to the config + code if so?

[john4255]

Hi, Thanks for reaching out. I am still seeing some of the same issues with hollow profiles near the core region. I have tried interpolating across the LH transition period to smooth the pedestal increase, but I am still getting very low temperature profiles (example shown below). Do you have any suggestions on what I could try to debug this? [image: 174657-ti.png]

…

On Tue, Nov 11, 2025 at 12:41 PM Jonathan Citrin ***@***.***> wrote: @john4255 <https://github.com/john4255> , apologies for dropping the ball here. What is the status here? Is help with your setup and interpretation still needed? Is there any update to the config + code if so? — Reply to this email directly, view it on GitHub <#1587 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BPSKCZ22KSXXXGAHVZTWTW334INVXAVCNFSM6AAAAACHKIF2EKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTIOJTHA4DKNA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[theo-brown]

Had a quick look at the config - you don't appear to have any particle sources, is this intentional?

[john4255]

Hi, Sorry, the config I sent earlier is for a different example (we decided to just use fixed density for that one, for the time being). The example I am testing currently (with evolving density) has gas puffing with S = 2.5e21 and decay length = 0.2. https://github.com/john4255/simple-pulse-design/blob/main/174657.py Let me know if I should send over any additional information. Thanks, John

…

On Thu, Nov 20, 2025 at 6:23 PM Theo Brown ***@***.***> wrote: Had a quick look at the config - you don't appear to have any particle sources, is this intentional? — Reply to this email directly, view it on GitHub <#1587 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BPSKCZZTTBCAAN2SMISTWF335ZENZAVCNFSM6AAAAACHKIF2EKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKMBTGI3TKMA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[theo-brown]

Could you add the JSON dump of the config for this run to the repo?

[john4255]

Hi, I added the JSON config to the repo. https://github.com/john4255/simple-pulse-design/blob/main/174657.json I also attached the temperature and density profile graphs to this email. Let me know if you have any questions. From, John

…

On Fri, Nov 21, 2025 at 2:43 AM Theo Brown ***@***.***> wrote: Could you add the JSON dump of the config for this run to the repo? — Reply to this email directly, view it on GitHub <#1587 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BPSKCZ76KAVASMLEHRIYXWL3527C3AVCNFSM6AAAAACHKIF2EKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKMBTGU2DMNY> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[john4255]

Hi all, I had a question about the way torax handles plasma current. I am trying to calculate the non-inductive current in a shot, and found that integrating j_total over area results in a value that is about ~5% off from the Ip value. Is this expected behavior? Or am I doing something wrong? Thanks, John

…

On Mon, Nov 24, 2025 at 1:00 PM John ***@***.***> wrote: Hi, I added the JSON config to the repo. https://github.com/john4255/simple-pulse-design/blob/main/174657.json I also attached the temperature and density profile graphs to this email. Let me know if you have any questions. From, John On Fri, Nov 21, 2025 at 2:43 AM Theo Brown ***@***.***> wrote: > Could you add the JSON dump of the config for this run to the repo? > > — > Reply to this email directly, view it on GitHub > <#1587 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/BPSKCZ76KAVASMLEHRIYXWL3527C3AVCNFSM6AAAAACHKIF2EKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKMBTGU2DMNY> > . > You are receiving this because you were mentioned.Message ID: > ***@***.*** > com> >

[jcitrin]

Hi John,

5% seems too much. How did you calculate the integral exactly?

The exact way to do it is with a Riemann sum. Since in TORAX we use the finite volume method, the cell values are considered averages over the cell, and a Riemann sum thus makes sense.

Using the test files (e.g. tests/test_data/test_iterhybrid_predictor_corrector.nc) and taking the last time point from the xarray DataTree (dt), let us take:

j_total = dt.profiles.j_total.to_numpy()[-1,:]
spr = dt.profiles.spr.to_numpy()[-1,:]
rho_norm = dt.rho_norm.to_numpy()

Note that all these profiles are on the output "cell + edges" grid, i.e. the cell grid plus the inner and outer faces at rho_norm=0,1 ,

Therefore the Riemann sum is:

Ip = np.sum(j_total[1:-1]*spr[1:-1])*drho_norm

where drho_norm = rho_norm[-2]-rho_norm[-3] , i.e. the cell grid spacing

This provides the same value as dt.scalars.Ip[-1] within floating point precision

Comparing to trapezoid integration, we can do:

Ip = np.trapz(j_total*spr, rho_norm) , and we get a number around 0.5% off the true value.

[john4255]

Hi Jonathan, Thanks for this. I tried running the code you sent, and it seems to be working. For some reason the error was very large at t=0.1 when I was using that as t_initial, but when I changed t_initial to 0.0 the code worked as expected. I am now trying to calculate a non-inductive FF' profile from the torax output, but the code I am using is not working. Could you provide an equation I can use to calculate non-inductive FF' from the torax outputs and/or would it be possible to make this an output field in the torax data tree output? Thanks, John

…

On Thu, Jan 15, 2026 at 2:26 PM Jonathan Citrin ***@***.***> wrote: Hi John, 5% seems too much. How did you calculate the integral exactly? The exact way to do it is with a Riemann sum. Since in TORAX we use the finite volume method, the cell values are considered averages over the cell, and a Riemann sum thus makes sense. Using the test files (e.g. tests/test_data/ test_iterhybrid_predictor_corrector.nc) and taking the last time point from the xarray DataTree (dt), let us take: j_total = dt.profiles.j_total.to_numpy()[-1,:] spr = dt.profiles.spr.to_numpy()[-1,:] rho_norm = dt.rho_norm.to_numpy() Note that all these profiles are on the output "cell + edges" grid, i.e. the cell grid plus the inner and outer faces at rho_norm=0,1 , Therefore the Riemann sum is: Ip = np.sum(j_total[1:-1]*spr[1:-1])*drho_norm where drho_norm = rho_norm[-2]-rho_norm[-3] , i.e. the cell grid spacing This provides the same value as dt.scalars.Ip[-1] within floating point precision Comparing to trapezoid integration, we can do: Ip = np.trapz(j_total*spr, rho_norm) , and we get a number around 0.5% off the true value. — Reply to this email directly, view it on GitHub <#1587 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BPSKCZ4X4NAVAHD3OWHNHDD4G7SV5AVCNFSM6AAAAACHKIF2EKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKNJRGA2DMNI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[jcitrin]

Hi John. Glad to hear it's improving.

We already calculate ffprime in the output, it should be in profiles.ffprime

Do you mean something different by "non-inductive FF'", as in calculate FF' purely with Ohmic current and neglect current drive and bootstrap current sources?

If so, you can calculate it with a variant of calc_FFprime

torax/torax/_src/physics/formulas.py

Line 112 in a2919e0

def calc_FFprime(

what is needed is to subtract the non-inductive contributions from j_total_face . Ensure that it's the toroidal projections. Everything is available for this calculation in the output.

I would first verify that it does what you expect before considering it as a native TORAX output. What is the use-case if specifically a non-inductive version of FF' ?

[john4255]

Hi, We're trying to calculate the portion of the FF' profile that comes from only non-inductive current (so bootstrap + current from ECCD/NBI but not ohmic current). We want to calculate the consumed flux but TokaMaker needs to know how much of the FF' profile is coming from non-inductive current to determine the contribution from the central solenoid. I will take a look at the function you linked to. Thanks, John

…

On Mon, Jan 26, 2026 at 1:54 PM Jonathan Citrin ***@***.***> wrote: Hi John. Glad to hear it's improving. We already calculate ffprime in the output, it should be in profiles.ffprime Do you mean something different by "non-inductive FF'", as in calculate FF' purely with Ohmic current and neglect current drive and bootstrap current sources? If so, you can calculate it with a variant of calc_FFprime https://github.com/google-deepmind/torax/blob/a2919e03a96cfa10ee7bdc7e56438e8d3bcef247/torax/_src/physics/formulas.py#L112 what is needed is to subtract the non-inductive contributions from j_total_face . Ensure that it's the toroidal projections. Everything is available for this calculation in the output. I would first verify that it does what you expect before considering it as a native TORAX output. What is the use-case if specifically a non-inductive version of FF' ? — Reply to this email directly, view it on GitHub <#1587 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BPSKCZ7SMCVEWQYNNADUEK34IZPHBAVCNFSM6AAAAACHKIF2EKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKNRQHE3DQNQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[hansec]

To clarify, we don't need it in the form of F*F'. Is it possible to get the non-inductive <J_phi> from TORAX? Or alternatively the inductive <J_phi>?

We just need to be able to distinguish the inductive from the non-inductive current.

[jcitrin]

Inductive j_phi is already there in there as profiles.j_ohmic (toroidal projection) and profiles.j_parallel_ohmic (parallel projection). The parallel currents (and a few other things) were not documented in the docs. I just added them in #1942

I'm now doing a quick PR to add non-inductive current density (toroidal and parallel projections), as well as total non-inductive current and non-inductive current scalars. These are all useful to have.

[jcitrin]

#1943 for new output profiles of non-inductive current density