Follow-up from "Resolve "convert stellarator.f90 to Python""

[jonmaddock]

In GitLab by @timothy-nunn on Apr 26, 2023, 13:55

The following discussion from !847 should be addressed:

• @jonmaddock started a discussion: (+3 comments)

  I don' think this is a good idea: the model running differently because it's aware of which constraints are being applied to it. I believe the returned values will be zero in the case where the constraint isn't active, which I find confusing, and this isn't done anywhere else to my knowledge. If it's purely performance reasoning, then I think this should probably be removed.

[jonmaddock]

Has this been dealt with already @timothy-nunn ?

[timothy-nunn]

@jonmaddock no. https://github.com/ukaea/PROCESS/blob/main/process/stellarator.py#L161

[timothy-nunn]

I will say this issue was annoying when we first encountered it. power_at_ignition_point had some side effects that changed the results for a Stellarator that did not enforce icc 91, hence why myself and JL opted to leave this logic alone.

[timothy-nunn]

Removing this condition check changes the stellarator_helias test:

WARNING  test_process_input_files:test_process_input_files.py:133 Variable                     Ref             New        % Change
WARNING  test_process_input_files:test_process_input_files.py:136 ------------------------------------------------------------
WARNING  test_process_input_files:test_process_input_files.py:138 res_eq_con016           2.09e-06       -4.31e-06         -306.65
WARNING  test_process_input_files:test_process_input_files.py:138 eq_con016              -2.09e-09        4.31e-09          306.65
WARNING  test_process_input_files:test_process_input_files.py:138 ineq_con083             1.85e-09       -2.88e-09         -255.51
WARNING  test_process_input_files:test_process_input_files.py:138 res_eq_con002           3.77e-11        1.26e-10          233.37
WARNING  test_process_input_files:test_process_input_files.py:138 eq_con002               1.64e-10        5.45e-10          233.36
WARNING  test_process_input_files:test_process_input_files.py:138 ineq_con024              5.7e-07        1.76e-06          208.25
WARNING  test_process_input_files:test_process_input_files.py:138 ineq_con035              3.3e-08        8.21e-08          148.94

Do we want to remove this condition if it is changing the solution (albeit not by much, since only the constraint tolerances are changing)? Also, (as the comment notes) running this model takes ages and massively increases the runtime of the tests.

[timothy-nunn]

@jjwalkowiak maybe you have a thought on this? For context, we are talking about making this small change which produces the above changes in one of the Stellarator regression tests. Do you think this is worth doing?

[jjwalkowiak]

TBH we usually don't use icc 91. If you don't like the condition checking icc, maybe replace it with flag in the setting, which will turn it off by default? Then it can be stated, that this flag have to be turned on for icc 91. I am a bit worried that this function changes the helias results - it should not do this. In general I don't trust this part of the code - I had no time analyze it in detail, but the general feeling is that it was some kind of workaround, which violated at least few good programming practices.

[timothy-nunn]

Hi @jjwalkowiak

Thanks for your quick reply. If it is something that you don't use and is a bit of hack, is it something we want to keep in PROCESS? If you want me to go down the setting route, I am happy to do this, but I also want to try and root out dubious bits of code that cannot be trusted.

Cheers,
Tim

[jjwalkowiak]

Well, I have mixed feelings. Technically speaking it works (or at least I hope so). Part of the problem is that current implementation is inefficient (it calls physics module few times, save some results and restores some variables to previous values etc.).
I can't say that I really trust it, but at the same time getting rid of this completely doesn't sound as a good solution. I would prefer to fix it, but I don't have time for this right now (or in predictable future).
If you want to be strict, you should probably remove it for now and I can restore it later. But maybe some intermediate solution can be used, to place it as non-default option, with remark that you should know what you do when you use this.

[timothy-nunn]

So having looked at the code a bit more, power_at_ignition_point is essentially trying to calculate:

• physics_variables.p_plasma_loss_mw
• physics_variables.pscalingmw

The first is calculated in:

PROCESS/process/models/physics/physics.py

Lines 1268 to 1301 in 67340b1

	(
	physics_variables.pden_electron_transport_loss_mw,
	physics_variables.pden_ion_transport_loss_mw,
	physics_variables.t_electron_energy_confinement,
	physics_variables.t_energy_confinement,
	physics_variables.t_ion_energy_confinement,
	physics_variables.p_plasma_loss_mw,
	) = self.calculate_confinement_time(
	physics_variables.m_fuel_amu,
	physics_variables.p_alpha_total_mw,
	physics_variables.aspect,
	physics_variables.b_plasma_toroidal_on_axis,
	physics_variables.nd_plasma_ions_total_vol_avg,
	physics_variables.nd_plasma_electrons_vol_avg,
	physics_variables.nd_plasma_electron_line,
	physics_variables.eps,
	physics_variables.hfact,
	physics_variables.i_confinement_time,
	physics_variables.i_plasma_ignited,
	physics_variables.kappa,
	physics_variables.kappa95,
	physics_variables.p_non_alpha_charged_mw,
	current_drive_variables.p_hcd_injected_total_mw,
	physics_variables.plasma_current,
	physics_variables.pden_plasma_core_rad_mw,
	physics_variables.rmajor,
	physics_variables.rminor,
	physics_variables.temp_plasma_electron_density_weighted_kev,
	physics_variables.temp_plasma_ion_density_weighted_kev,
	physics_variables.q95,
	physics_variables.qstar,
	physics_variables.vol_plasma,
	physics_variables.n_charge_plasma_effective_vol_avg,
	)

And the second is calculated with stuff that is also calculated above/wouldn't change.

So we can replace the calling of the physics model with essentially a little bit of copy pasting from physics. We can test this by adding icc=91 before and after.

If this fails, the new data structure should allow the physics model to run in isolation with a 'mock' data structure without side-effects to the main data structure.

[jjwalkowiak]

In general I think that this is the way. The problem is that you need to copy-paste more than a little bit of physics module to make it work. If you look into p_plasma_loss_mw calculation, you will see that it depends on other parameters which have to be recalculated first, like fusion reaction rate and radiation. But if you are familiar with the content of the physics module, maybe it would be easier for you - I would have to first analyze what is calculated where and so on.
BTW I think it would be nice to separate p_plasma_loss_mw into another staticmethod, because now it is a bit hidden inside confinement time.

[grmtrkngtn]

@timothy-nunn could you link to the triaged issue to resolve this?