Newton solve crashing when calculating `t_zero_margin`

[timothy-nunn]

I get the following traceback

 **************************************************************************************************************
/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/physics.py:7255: RuntimeWarning: divide by zero encountered in divide
  * (neped / n_greenwald) ** -0.174
/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/superconductors.py:754: RuntimeWarning: invalid value encountered in scalar power
  b_critical = b_c20_eps * (1.0 - f_temp_conductor_critical_no_field**1.52)
Traceback (most recent call last):
  File "/home/ubuntu/PROCESS-convergence/.venv/bin/process", line 33, in <module>
    sys.exit(load_entry_point('process==3.1.0', 'console_scripts', 'process')())
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/main.py", line 740, in main
    Process(args)
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/main.py", line 136, in __init__
    self.run_mode()
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/main.py", line 236, in run_mode
    self.run.run()
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/main.py", line 393, in run
    self.run_scan(self.solver)
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/main.py", line 487, in run_scan
    self.scan = Scan(self.models, solver)
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/scan.py", line 157, in __init__
    self.run_scan()
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/scan.py", line 171, in run_scan
    ifail = self.doopt()
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/scan.py", line 194, in doopt
    ifail = self.optimiser.run()
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/optimiser.py", line 54, in run
    ifail = self.solver.solve()
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/solver.py", line 211, in solve
    x, _, _, res = solve(
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/pyvmcon/vmcon.py", line 146, in solve
    result = problem(x)
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/solver.py", line 133, in __call__
    objf, conf = self._evaluator.fcnvmc1(n, self.total_constraints, x, 0)
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/evaluators.py", line 55, in fcnvmc1
    objf, conf = self.caller.call_models(xv, m)
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/caller.py", line 75, in call_models
    self._call_models_once(xc)
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/caller.py", line 270, in _call_models_once
    self.models.sctfcoil.run(output=False)
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/superconducting_tf_coil.py", line 306, in run
    ) = self.supercon(
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/process/superconducting_tf_coil.py", line 1111, in supercon
    t_zero_margin, root_result = optimize.newton(
  File "/home/ubuntu/PROCESS-convergence/.venv/lib/python3.10/site-packages/scipy/optimize/_zeros_py.py", line 391, in newton
    raise RuntimeError(msg)
RuntimeError: Failed to converge after 50 iterations, value is nan.

when running fails_supercon.IN.DAT.txt

Anecdotally, a lot of people experience this crash, especially when doing Monte Carlo on the large tokamak.

We should try and find the cause of this crash and confirm there are no bugs at play here. We should also try and stop this from crashing the code, is there an estimation or alternative method we can fall back on?

[timothy-nunn]

For this input file, the current_density_margin looks as follows on the range t=[0,5]

(Other than at 0) the function monotonically decreases as t increases. Because t~=0 is negative (t=0 behaves weirdly) there is no root for positive t. Therefore, the Newton algorithm attempts to go into negative t to find the point where current_density_margin = 0 however it encounters a NaN for any negative t due to a root.

For reference, normally the curve looks like

[mkovari]

I am a bit depressed that I didn't do something to catch this error when I put in the optimize.newton call.

Can anyone clarify - why does the error message say Failed to converge after 50 iterations, value is nan.? Does this mean that that optimize.newton can iterate even when the function being called returns nan?

[mkovari]

I think the best thing is to add one of our kludges so the relevant functions return a real value of current density even for negative temperature. They are as follows:

    # Critical current density jcrit
    if isumat == 1:
        jcrit, _, _ = itersc(ttest, bmax, strain, bc20m, tc0m)
    elif isumat == 3:
        jcrit, _ = jcrit_nbti(ttest, bmax, c0, bc20m, tc0m)
    if isumat == 4:
        jcrit, _, _ = itersc(ttest, bmax, strain, bc20m, tc0m)
    elif isumat == 5:
        jcrit, _, _ = western_superconducting_nb3sn(ttest, bmax, strain, bc20m, tc0m)
    elif isumat == 7:
        jcrit, _, _ = gl_nbti(ttest, bmax, strain, bc20m, tc0m)
    elif isumat == 8:
        jcrit, _, _ = gl_rebco(ttest, bmax, strain, bc20m, tc0m)
    elif isumat == 9:
        jcrit, _, _ = hijc_rebco(

I am pretty sure this would work, since negative values of t_zero_margin and temp_margin should both be fine. Kludges are described in #3578 or here.

[timothy-nunn]

I am a bit depressed that I didn't do something to catch this error when I put in the optimize.newton call.

Can anyone clarify - why does the error message say Failed to converge after 50 iterations, value is nan.? Does this mean that that optimize.newton can iterate even when the function being called returns nan?

I think that it can continue to run while the return is nan but it obviously has no information to decide what t to try next and so probably just randomly moves around in t until hitting the maximum iterations

[mkovari]

I think it is used by constraint 36, but somewhere its name has changed:

tmp_cc =  1.0D0 - ftmargtf * tmargtf/tmargmin_tf

[timothy-nunn]

Check tmargtf and temp_margin are the same. Can they be consolidated

[mkovari]

Hi Tim, what's the situation with this?

[timothy-nunn]

@mkovari its in my backlog of things to do but no time frame yet