Impurity radiation corrections by jonmaddock · Pull Request #4536 · ukaea/PROCESS

jonmaddock · 2026-08-16T09:12:14Z

Demonstrates the effect of various options for impurity radiation and allows their plotting: this is not intended to be merged in its current form, but as an investigation. Primarily, this demonstrates the importance of separating the reduced core radiation used for the confinement time calculation and the total radiation considered in the plasma power balance, which I currently believe is incorrect.

I hope the following figures highlight a few things and justify the actual changes I want to make. They are all made with the DEMO LAR solution point and varying tungsten fractions. Deliberate horizontal "jitter" is applied to the scatter plots to ensure overlapping points are discernable.

Impurity radiation profiles

The current impurity power densities are incorrect because they are multiplied prematurely by the normalised minor radius $\rho$ in preparation for the power integral. This has the effect of the power density dropping to 0 for $\rho = 0$, which is incorrect. The core radiation reduction is also evident outside of the confinement time scaling calculation, where it is only intended to be used. This is corrected. The edge integral is now performed explicitly in the 4th subfigure, which is used later subsequently.

Radiation power

When the core radiation reduction is removed, the core radiation obviously increases, but the edge radiation decreases, due to the "edge = total - core" calculation that is currently used. The edge is optionally now integrated explicitly. The total radiation is the same in all cases.

Separate radiation reduction for confinement time scaling only

Two plasma core radiation power densities are now calculated: one for the confinement time scaling only, the other for the rest of the code. They only differ when the `tau E rad reduction` option is selected (brown), which is the only difference between the plots. This separates the radiation calculated for the rest of the models from the reduced value required by the scaling.

Effect on plasma power balance (PPB)

The effect of these various options on the absolute plasma power balance constraint residual in MW (i.e. MW m^-3, then multiplied by plasma volume ) is shown. The "original" point being slightly above zero is probably due to using a slightly out-of-date solution vector relative to the current Process version. One real difference occurs from using the edge radiation in the PPB (purple) which should always be done in my view: the plasma is not in thermal equilibrium if the edge radiation is ignored. The second, larger, difference is due to the use of the reduced core radiation in the confinement time calculation only, rather than globally (i.e. in the PPB too). This is how this "correction" was originally intended to be used I believe: as a way of effectively changing the confinement time scaling for highly-radiative regimes. It should not affect the actual calculated radiation elsewhere in the models. Around 30 MW for low tungsten values is fairly significant.

I would like to always use the un-reduced total radiation (core and edge) in the plasma power balance calculation, and the reduced core radiation in the confinement time calculation only. This would mean that i_rad_loss and f_p_plasma_core_rad_reduction would only apply to the the confinement time calculation. I'm keen to hear your thoughts.

… tauE calculation Leave plasma power balance impurity radiation un-reduced.

jonmaddock · 2026-08-16T10:24:55Z

This appears to independently resolve #3968 (comment).

mkovari · 2026-08-17T07:26:12Z

this is not intended to be merged in its current form, but as an investigation.

Can you clarify this?

je-cook · 2026-08-17T07:56:47Z

this is not intended to be merged in its current form, but as an investigation.

Can you clarify this?

I think he means it should probably be split into a few separate PRs after the initial review discussion. Also I guess the switches wont be needed if we agree on the solution

mkovari · 2026-08-17T08:08:24Z

This appears to independently resolve #3968 (comment).

This link is to the comment on volume-integration. Was that intended?

mkovari · 2026-08-17T10:08:55Z

I think he means it should probably be split into a few separate PRs after the initial review discussion. Also I guess the switches wont be needed if we agree on the solution

So do we actually need to approve the PR?

je-cook · 2026-08-17T10:09:57Z

no but a review is needed of the content

grmtrkngtn

Impurity radiation profiles

I agree with the proposed changes to the impurity radiation profiles. Moving the (\rho) weighting from the stored power-density profile into the volume integral is a clear-cut fix.

Core radiation reduction

I also agree with separating the physical core radiation from the reduced value used in the confinement-time calculation. Retaining the option to apply the reduction to the scaling law only seems like a reasonable compromise but the physical power balance must use the full radiation.

Outstanding

The corrected confinement calculation still relies on two input assumptions: radius_plasma_core_norm and
f_p_plasma_core_rad_reduction.

I'm not sure we still still need the latter, why would we ever not subtract the full radiation (why would we use the factor?)

EUROfusion colleagues have also suggested that the effective core boundary could be calculated from the radial power flow rather than prescribed as an input. My understanding is that METIS calculates the cumulative net power, heating minus local radiation, as a function of radius and uses its maximum to determine the effective transport-driving power. Should we investigate this here, or would that be better treated as a separate follow-up to this accounting PR?

jonmaddock · 2026-08-24T08:56:13Z

@mkovari these changes solve #3968 (the issue, not just the comment within the issue) as it separates the power density calculation and its integral. This draft PR is intended to show how the above plots were generated to better understand the problem due to the confusion of multiple issues acting together and the overall effect of the changes I intend to make, which will be a separate PR. These changes are described in the final paragraph of the PR description.

@grmtrkngtn

The corrected confinement calculation still relies on two input assumptions: radius_plasma_core_norm and
f_p_plasma_core_rad_reduction.
I'm not sure we still still need the latter, why would we ever not subtract the full radiation (why would we use the factor?)

That is how the confinement time scaling law has effectively been modified (according to Lux et. al.) to model the highly-radiative regime: of the core region (radius_plasma_core_norm) only use some fraction of the radiation (f_p_plasma_core_rad_reduction). Changing this would require performing Lux et. al.'s analysis again to achieve a different way of modifying the scaling to fit highly-radiative data, which is outside the scope of this task. I'm just interested in using the correct amount of radiation in the plasma power balance calculation here.

Any calculation of the core boundary is an additional feature.

I'd be grateful for a review from @mkovari , @chris-ashe and @ajpearcey .

mkovari

I can't say I have read it all but I am keen to see what this looks like when merged into main.

ajpearcey · 2026-08-24T09:50:38Z

I would be interested to see how the changed PPB changes the final solution vector for the DEMO-LAR. How big is this shift in the Psep value? Also was the edge radiation also previously excluded in the FW heat load calculations?

chris-ashe · 2026-08-24T10:47:35Z

    """total core radiation power per volume (MW/m3)"""

+    pden_plasma_core_rad_tauE_mw: float = 0.0
+    """reduced total core radiation power for tauE calculation (MW/m3)"""


Suggested change

"""reduced total core radiation power for tauE calculation (MW/m3)"""

"""reduced total core radiation power for τₑ calculation [MW/m³]"""

chris-ashe · 2026-08-24T10:47:48Z


-        crfcdr = (fefcdr * self.data.costs.discount_rate) / (fefcdr - 1.0e0)
+        # crfcdr = (fefcdr * self.data.costs.discount_rate) / (fefcdr - 1.0e0)
+        crfcdr = 1.0


Not sure why this is now just set to 1

chris-ashe · 2026-08-24T10:49:21Z

+            # but reduced radiation here in confinement time calculation
            if model == ConfinementRadiationLossModel.FULL_RADIATION:
-                p_plasma_loss_mw -= self.data.physics.pden_plasma_rad_mw * vol_plasma
+                if impurity_radiation.rad_reduction_for_tauE_only:


Should this not be in the ConfinementRadiationLossModel.CORE_ONLY condition below?

chris-ashe · 2026-08-24T10:57:17Z

+        rho = 1.0
+        if rho_fix:
+            rho = self.plasma_profile.neprofile.profile_x
        self.pden_impurity_rad_total_mw = 2.0e-6 * integrate.simpson(


I believe it is better here to rename these variables to have _vol_avg in them as thats the returned value from the integration when multiplying by rho at each point. These values are then multiplied by vol_plasma to get the total radiation power. This should make it more clear as to what is happening

chris-ashe · 2026-08-24T10:59:05Z

                    )
+
+                    # Now idempotent, return
+                    # Pass model caller and opt params for stability constraint evaluation


Not sure about this change, probs best if someone like @timothy-nunn checks it

This PR does not include objective function 20 so these changes are unnecessary at this point.

Also, call_models_and_write_output returns the objective value and normalised residuals, but I cannot see them being used anywhere in this PR.

I suggest all of these changes in caller can be left for a later PR.

chris-ashe · 2026-08-24T11:00:10Z

-    trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd
+    trying to optimise it.
+    """
+    SOLUTION = (-1, "Solution mode (no optimisation)")


Did we accidentaly remove solution mode when setting up the enum @jonmaddock @timothy-nunn ?

Solution mode was removed a long time ago, it was replaced by evaluation mode. This is because solution mode did not produce consistent solutions and so was essentially useless (except for testing, see #4044)

chris-ashe · 2026-08-24T11:01:31Z

+            f_p_plasma_core_rad_reduction = 1.0
+        else:
+            # Reduce radiation in PPB as well
+            f_p_plasma_core_rad_reduction = f_p_plasma_core_rad_reduction_tauE


I assume this will be removed once we put the PR up for the PPB fix?

chris-ashe · 2026-08-24T11:03:34Z

+            )
+
+            # Total radiation profile (core + edge)
+            pden_impurity_rad_total = (


If you multiply the core and edge pden by the total plasma volume and then add them together does it equal pden_impurity_rad_total * vol_plasma?

chris-ashe · 2026-08-24T11:06:08Z

            x=self.plasma_profile.neprofile.profile_x,
            dx=self.plasma_profile.neprofile.profile_dx,
        )
        self.pden_impurity_core_rad_total_mw = 2.0e-6 * integrate.simpson(


For the core value I think we need to be explicit in stating that this is till done as a full plasma volume integral with the mask of 0's where the core would be. This has always caused confusion as to if the core or edge value is only done as a function of its own volume or that of the total plasma

chris-ashe · 2026-08-24T11:07:53Z

-        elif self.data.numerics.i_process_run_mode == PROCESSRunMode.EVALUATION:
-            # No optimisation:
-            # solve equality (consistency) constraints only using fsolve (HYBRD)
+        # ioptimz == -1: solution


Again will leave this one for @timothy-nunn

jonmaddock added 6 commits August 16, 2026 10:29

Add edge radiation profile calculations to ImpurityRadiation class

ab9c152

Add options to fix rho bug and include edge radiation for impurities

89a325f

Reinstate evaluation mode

cbecb1d

Silence silly costs bug

8f8926b

Store constraint values on each iteration for wider access

8282bff

Perform separate radiation reduction for impurity-radiation-corrected…

a9b3913

… tauE calculation Leave plasma power balance impurity radiation un-reduced.

jonmaddock force-pushed the impurity-radiation-corrections branch from 53719e4 to a9b3913 Compare August 16, 2026 09:39

jonmaddock requested review from ajpearcey, chris-ashe, grmtrkngtn and mkovari August 16, 2026 10:23

grmtrkngtn reviewed Aug 19, 2026

View reviewed changes

mkovari approved these changes Aug 24, 2026

View reviewed changes

chris-ashe requested changes Aug 24, 2026

View reviewed changes

	"""reduced total core radiation power for tauE calculation (MW/m3)"""
	"""reduced total core radiation power for τₑ calculation [MW/m³]"""

Uh oh!

Conversation

jonmaddock commented Aug 16, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Impurity radiation profiles

Radiation power

Separate radiation reduction for confinement time scaling only

Effect on plasma power balance (PPB)

Uh oh!

jonmaddock commented Aug 16, 2026

Uh oh!

mkovari commented Aug 17, 2026

Uh oh!

je-cook commented Aug 17, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mkovari commented Aug 17, 2026

Uh oh!

mkovari commented Aug 17, 2026

Uh oh!

je-cook commented Aug 17, 2026

Uh oh!

grmtrkngtn left a comment

Choose a reason for hiding this comment

Impurity radiation profiles

Core radiation reduction

Outstanding

Uh oh!

jonmaddock commented Aug 24, 2026

Uh oh!

mkovari left a comment

Choose a reason for hiding this comment

Uh oh!

ajpearcey commented Aug 24, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

7 participants

jonmaddock commented Aug 16, 2026 •

edited

Loading

je-cook commented Aug 17, 2026 •

edited

Loading