Burn time consistency loop

[mkovari]

Discussed in #3048

In Power balance output and tburn consistency #17 I added an extra extra loop a long time ago in fcnvm1 (evaluators.py) to ensure that tburn has reached a converged value. This will run all the physics and engineering functions up to 10 times, so it could be contributing to slow run times.

        # Evaluate machine parameters at xv
        self.caller.call_models(xv)

        # Convergence loop to ensure burn time consistency
        if sv.istell == 0:
            for _ in range(10):
                if abs((tv.tburn - tv.tburn0) / max(tv.tburn, 0.01)) <= 0.001:
                    break

                self.caller.call_models(xv)
                if gv.verbose == 1:
                    print("Internal tburn consistency check: ", tv.tburn, tv.tburn0)

See comments below.
Conclusion: no change is required now.
However, if and when PROCESS is changed so that all the function evaluations are repeated several times, this loop should be removed.
@jonmaddock @timothy-nunn

[mkovari]

It seems that the problem arises because the calculation of loop voltage vburn is duplicated and the whole thing is oddly circular.
Pseudocode:

        physics.physics()
        def physics()
              tburn0 = tburn
              tpulse = ... + tburn + ...
              # Calculate Volt-second requirements
              ..., vsstt = vscalc(...tburn...)
             def vscalc():
                 # Volt-second requirements
                 # Loop voltage during flat-top
                 vburn = plascur * rplas * facoh * csawth
                 vsbrn = vburn * (t_fusion_ramp + tburn)
                 vsstt = vsstt + vsbrn

        pfcoil.run()
             pfcoil.py
            def run():
            # Volt-second capability of PF coil set
            self.vsec()
           def vsec():
                   # Calculation of volt-second capability of PF system.
                   # vsbn = total flux swing available for burn (Wb)
                   pfv.vsbn = .....
        
        pulse.run()
        pulse.py
        run()
           burn()
               # Routine to calculate the burn time for a pulsed reactor
               if pulse_variables.lpulse != 1:
                     return
               vsmax = -pfcoil_variables.vsbn
              # Loop voltage during flat-top 
              vburn = plascur * rplas * facoh * csawth
              tburn = vsmax / vburn - t_fusion_ramp

[mkovari]

In principle I would like to eliminate the switch lpulse, since any tokamak burn is limited in duration, and can therefore be considered pulsed. However, the STEP input file uses the "steady-state" option (lpulse = 0), so it's not wise to fix something that's working and needed.

Name	Type	Datatype	Default Value	Description
lpulse	Output	integer	-	Switch for reactor model: =0 continuous operation =1 pulsed operation
tburn	Input	real	1000.0	burn time (s) (calculated if lpulse=1)

[mkovari]

It turns out that the consistency loop above is not used often.

In the large-tokamak scenario, the total number of function evaluations is ft.numerics.ncalls = 1364, while the loop is entered only 13 times. The vast majority of calls to call_models are in the gradient function evaluator fcnvmc2, since a gradient is evaluated for every iteration variable, requiring two calls each.

[mkovari]

I have tried replacing the line that ensures tburn is positive (times_variables.tburn = max(0.0e0, tb)) by a smooth function that is always positive:

Original code:

times_variables.tburn = max(0.0e0, tb)

Smooth function:

        # if calculated burn time is less than 10 seconds use a smooth positive function.
        if(tb < 10):
            times_variables.tburn = 100 / ( 20 - tb )        
        else:
            times_variables.tburn = tb

Allow negative tburn:

times_variables.tburn = tb

(I have used error tolerance for VMCON epsvmc = 1e-4 to save time.)

I tried this out in a version of large-tokamak in which the thickness of the CS was gradually reduced:

nsweep = 63 * ohcth : CS thickness (m) 
isweep = 5
sweep = 0.42, 0.4, 0.39, 0.38, 037

All three options give exactly the same values for nviter and the cumulative number of function calls ncalls.

I have also tried removing the consistency loop just for fun. The scan above crashes on the first scan point.

Changing the scan:

sweep = 0.55, 0.4, 0.39, 0.38, 037

gives:
No loop:

• nviter = 29, 11, 5, 5
• ncalls = 5017, 6866, 7659, 8452, crash.

With loop:

• nviter = 22, 10, 5, 5, crash
• ncalls = 3817, 5499, 6293, 7087, crash

The loop seems to help a little in this one example.