Convert superconductors.f90 to Python

CMakeLists.txt

@@ -81,7 +81,6 @@ LIST(APPEND PROCESS_SRCS
     reinke_module.f90
     sctfcoil.f90
     plasma_profiles.f90
-    superconductors.f90
     final_module.f90
     cost_variables.f90
     divertor_variables.f90

process/pfcoil.py

@@ -11,12 +11,12 @@
 from process.fortran import maths_library as ml
 from process.fortran import process_output as op
 from process.fortran import numerics
-from process.fortran import superconductors as sc
 from process.fortran import rebco_variables as rcv
 from process.fortran import constraint_variables as ctv
 from process import maths_library as pml
 from process.utilities.f2py_string_patch import f2py_compatible_to_string
 from process import fortran as ft
+import process.superconductors as superconductors
 import math
 import numpy as np
 import numba
@@ -175,7 +175,6 @@ def pfcoil(self):
 
         # N.B. Problems here if k=ncls(group) is greater than 2.
         for j in range(pfv.ngrp):
-
             if pfv.ipfloc[j] == 1:
                 # PF coil is stacked on top of the Central Solenoid
                 for k in pfv.ncls[j]:
@@ -333,7 +332,6 @@ def pfcoil(self):
                 )
 
             else:
-
                 # Conventional aspect ratio scaling
                 nfxf0 = 0
                 ngrp0 = 0
@@ -554,7 +552,6 @@ def pfcoil(self):
 
         for ii in range(pfv.ngrp):
             for ij in range(pfv.ncls[ii]):
-
                 if pfv.ipfloc[ii] == 1:
                     # PF coil is stacked on top of the Central Solenoid
                     dx = 0.5e0 * bv.ohcth
@@ -588,7 +585,6 @@ def pfcoil(self):
                         pfv.zh[i] = pfv.zpf[i] - dz
 
                 else:
-
                     # Other coils. N.B. Current density RJCONPF[i] is defined in
                     # routine INITIAL for these coils.
                     area = abs(pfv.ric[i] * 1.0e6 / pfv.rjconpf[i])
@@ -624,7 +620,6 @@ def pfcoil(self):
         for ii in range(pfv.ngrp):
             iii = ii
             for ij in range(pfv.ncls[ii]):
-
                 # Peak field
 
                 if ij == 0:
@@ -1023,7 +1018,6 @@ def ohcalc(self):
             sgn = 1.0e0
             pfv.ric[pfv.nohc - 1] = sgn * 1.0e-6 * pfv.cohbop * pfv.areaoh
         else:
-
             sgn = -1.0e0
             pfv.ric[pfv.nohc - 1] = sgn * 1.0e-6 * pfv.coheof * pfv.areaoh
 
@@ -1104,7 +1098,6 @@ def ohcalc(self):
 
         # Stress ==> cross-sectional area of supporting steel to use
         if pfv.ipfres == 0:
-
             # Superconducting coil
 
             # New calculation from M. N. Wilson for hoop stress
@@ -1192,7 +1185,6 @@ def ohcalc(self):
             )
 
         if pfv.ipfres == 0:
-
             # Allowable coil overall current density at EOF
             # (superconducting coils only)
 
@@ -1306,7 +1298,6 @@ def peakb(self, i, ii, it):
         jj = 0
         for iii in range(pfv.ngrp):
             for jjj in range(pfv.ncls[iii]):
-
                 jj = jj + 1
                 # Radius, z-coordinate and current for each coil
                 if iii == ii - 1:
@@ -1419,32 +1410,27 @@ def bfmax(self, rj, a, b, h):
         )
 
         if beta > 3.0:
-
             b1 = constants.rmu0 * rj * (b - a)
             f = (3.0 / beta) ** 2
             bfmax = f * b0 * (1.007 + (alpha - 1.0) * 0.0055) + (1.0 - f) * b1
 
         elif beta > 2.0:
-
             rat = (1.025 - (beta - 2.0) * 0.018) + (alpha - 1.0) * (
                 0.01 - (beta - 2.0) * 0.0045
             )
             bfmax = rat * b0
 
         elif beta > 1.0:
-
             rat = (1.117 - (beta - 1.0) * 0.092) + (alpha - 1.0) * (beta - 1.0) * 0.01
             bfmax = rat * b0
 
         elif beta > 0.75:
-
             rat = (1.30 - 0.732 * (beta - 0.75)) + (alpha - 1.0) * (
                 0.2 * (beta - 0.75) - 0.05
             )
             bfmax = rat * b0
 
         else:
-
             rat = (1.65 - 1.4 * (beta - 0.5)) + (alpha - 1.0) * (
                 0.6 * (beta - 0.5) - 0.20
             )
@@ -1769,7 +1755,6 @@ def induct(self, output):
                 ]
 
         if bv.iohcl != 0:
-
             # Central Solenoid self inductance
             a = pfv.rohc  # mean radius of coil
             b = 2.0e0 * pfv.zh[pfv.nohc - 1]  # length of coil
@@ -2704,7 +2689,6 @@ def waveform(self):
             pfv.waves[nplas - 1, it] = 1.0e0
 
         for ic in range(pfv.nohc):
-
             # Find where the peak current occurs
             # Beginning of pulse, t = tramp
             if (abs(pfv.curpfs[ic]) >= abs(pfv.curpfb[ic])) and (
@@ -2794,7 +2778,7 @@ def deltaj_nbti(temperature):
             :return: difference in current density
             :rtype: float
             """
-            jcrit0, t = sc.jcrit_nbti(temperature, bmax, c0, bc20m, tc0m)
+            jcrit0, _ = superconductors.jcrit_nbti(temperature, bmax, c0, bc20m, tc0m)
             if ml.variable_error(jcrit0):  # superconductors.jcrit_nbti has failed.
                 print(f"superconductors.jcrit_nbti: {bmax=} {temperature=} {jcrit0=}")
 
@@ -2809,7 +2793,7 @@ def deltaj_wst(temperature):
             :return: difference in current density
             :rtype: float
             """
-            jcrit0, b, t = sc.wstsc(temperature, bmax, strain, bc20m, tc0m)
+            jcrit0, _, _ = superconductors.wstsc(temperature, bmax, strain, bc20m, tc0m)
             if ml.variable_error(jcrit0):  # superconductors.wstsc has failed.
                 print(f"deltaj_wst: {bmax=} {temperature=} {jcrit0=}")
 
@@ -2824,7 +2808,9 @@ def deltaj_gl_nbti(temperature):
             :return: difference in current density
             :rtype: float
             """
-            jcrit0, b, t = sc.gl_nbti(temperature, bmax, strain, bc20m, tc0m)
+            jcrit0, _, _ = superconductors.gl_nbti(
+                temperature, bmax, strain, bc20m, tc0m
+            )
             if ml.variable_error(jcrit0):  # GL_Nbti has failed.
                 print(f"deltaj_GL_nbti: {bmax=} {temperature=} {jcrit0=}")
 
@@ -2839,7 +2825,9 @@ def deltaj_gl_rebco(temperature):
             :return: difference in current density
             :rtype: float
             """
-            jcrit0, b, t = sc.gl_rebco(temperature, bmax, strain, bc20m, tc0m)
+            jcrit0, _, _ = superconductors.gl_rebco(
+                temperature, bmax, strain, bc20m, tc0m
+            )
             if ml.variable_error(jcrit0):  # superconductors.GL_REBCO has failed.
                 print(f"deltaj_gl_REBCO: {bmax=} {temperature=} {jcrit0=}")
 
@@ -2854,7 +2842,9 @@ def deltaj_hijc_rebco(temperature):
             :return: difference in current density
             :rtype: float
             """
-            jcrit0, b, t = sc.hijc_rebco(temperature, bmax, strain, bc20m, tc0m)
+            jcrit0, _, _ = superconductors.hijc_rebco(
+                temperature, bmax, strain, bc20m, tc0m
+            )
             if ml.variable_error(jcrit0):  # superconductors.GL_REBCO has failed.
                 print(f"deltaj_hijc_REBCO: {bmax=} {temperature=} {jcrit0=}")
 
@@ -2870,7 +2860,7 @@ def deltaj_hijc_rebco(temperature):
             # jcritsc returned by superconductors.itersc is the critical current density in the
             # superconductor - not the whole strand, which contains copper
 
-            jcritsc, bcrit, tcrit = sc.itersc(thelium, bmax, strain, bc20m, tc0m)
+            jcritsc, _, _ = superconductors.itersc(thelium, bmax, strain, bc20m, tc0m)
             jcritstr = jcritsc * (1.0e0 - fcu)
 
         elif isumat == 2:
@@ -2885,23 +2875,22 @@ def deltaj_hijc_rebco(temperature):
 
             jstrand = jwp / (1.0e0 - fhe)
 
-            jcritstr, tmarg = sc.bi2212(bmax, jstrand, thelium, fhts)
+            jcritstr, tmarg = superconductors.bi2212(bmax, jstrand, thelium, fhts)
             jcritsc = jcritstr / (1.0e0 - fcu)
-            tcrit = thelium + tmarg
 
         elif isumat == 3:
             # NbTi data
             bc20m = 15.0e0
             tc0m = 9.3e0
             c0 = 1.0e10
-            jcritsc, tcrit = sc.jcrit_nbti(thelium, bmax, c0, bc20m, tc0m)
+            jcritsc, _ = superconductors.jcrit_nbti(thelium, bmax, c0, bc20m, tc0m)
             jcritstr = jcritsc * (1.0e0 - fcu)
 
         elif isumat == 4:
             # As (1), but user-defined parameters
             bc20m = bcritsc
             tc0m = tcritsc
-            jcritsc, bcrit, tcrit = sc.itersc(thelium, bmax, strain, bc20m, tc0m)
+            jcritsc, _, _ = superconductors.itersc(thelium, bmax, strain, bc20m, tc0m)
             jcritstr = jcritsc * (1.0e0 - fcu)
 
         elif isumat == 5:
@@ -2912,12 +2901,12 @@ def deltaj_hijc_rebco(temperature):
             # jcritsc returned by superconductors.itersc is the critical current density in the
             # superconductor - not the whole strand, which contains copper
 
-            jcritsc, bcrit, tcrit = sc.wstsc(thelium, bmax, strain, bc20m, tc0m)
+            jcritsc, _, _ = superconductors.wstsc(thelium, bmax, strain, bc20m, tc0m)
             jcritstr = jcritsc * (1.0e0 - fcu)
 
         elif isumat == 6:
             # "REBCO" 2nd generation HTS superconductor in CrCo strand
-            jcritsc, validity = sc.jcrit_rebco(thelium, bmax, 0)
+            jcritsc, _ = superconductors.jcrit_rebco(thelium, bmax)
             jcritstr = jcritsc * (1.0e0 - fcu)
 
             # The CS coil current at EOF
@@ -2931,7 +2920,7 @@ def deltaj_hijc_rebco(temperature):
             # Durham Ginzburg-Landau critical surface model for Nb-Ti
             bc20m = tfv.b_crit_upper_nbti
             tc0m = tfv.t_crit_nbti
-            jcritsc, bcrit, tcrit = sc.gl_nbti(thelium, bmax, strain, bc20m, tc0m)
+            jcritsc, _, _ = superconductors.gl_nbti(thelium, bmax, strain, bc20m, tc0m)
             jcritstr = jcritsc * (1.0e0 - fcu)
 
             # The CS coil current at EOF
@@ -2941,7 +2930,7 @@ def deltaj_hijc_rebco(temperature):
             # Durham Ginzburg-Landau critical surface model for REBCO
             bc20m = 429e0
             tc0m = 185e0
-            jcritsc, bcrit, tcrit = sc.gl_rebco(thelium, bmax, strain, bc20m, tc0m)
+            jcritsc, _, _ = superconductors.gl_rebco(thelium, bmax, strain, bc20m, tc0m)
             # A0 calculated for tape cross section already
             jcritstr = jcritsc * (1.0e0 - fcu)
 
@@ -2954,7 +2943,9 @@ def deltaj_hijc_rebco(temperature):
             # Hazelton experimental data + Zhai conceptual model for REBCO
             bc20m = 138
             tc0m = 92
-            jcritsc, bcrit, tcrit = sc.hijc_rebco(thelium, bmax, strain, bc20m, tc0m)
+            jcritsc, _, _ = superconductors.hijc_rebco(
+                thelium, bmax, strain, bc20m, tc0m
+            )
             # A0 calculated for tape cross section already
             jcritstr = jcritsc * (1.0e0 - fcu)
 
@@ -2995,14 +2986,20 @@ def deltaj_hijc_rebco(temperature):
                 ttestp = ttest + delt
 
                 if isumat in [1, 4]:
-                    jcrit0, b, t = sc.itersc(ttest, bmax, strain, bc20m, tc0m)
+                    jcrit0, _, _ = superconductors.itersc(
+                        ttest, bmax, strain, bc20m, tc0m
+                    )
                     if (abs(jsc - jcrit0) <= jtol) and (
                         abs((jsc - jcrit0) / jsc) <= 0.01
                     ):
                         break
 
-                    jcritm, b, t = sc.itersc(ttestm, bmax, strain, bc20m, tc0m)
-                    jcritp, b, t = sc.itersc(ttestp, bmax, strain, bc20m, tc0m)
+                    jcritm, _, _ = superconductors.itersc(
+                        ttestm, bmax, strain, bc20m, tc0m
+                    )
+                    jcritp, _, _ = superconductors.itersc(
+                        ttestp, bmax, strain, bc20m, tc0m
+                    )
 
                 # Kludge to avoid divide by 0
                 if jcritm == jcritp:
@@ -3021,7 +3018,9 @@ def deltaj_hijc_rebco(temperature):
                 deltaj_nbti, x1, x2, 100e0
             )
             tmarg = current_sharing_t - thelium
-            jcrit0, t = sc.jcrit_nbti(current_sharing_t, bmax, c0, bc20m, tc0m)
+            jcrit0, _ = superconductors.jcrit_nbti(
+                current_sharing_t, bmax, c0, bc20m, tc0m
+            )
             if ml.variable_error(
                 current_sharing_t
             ):  # current sharing secant solver has failed.
@@ -3039,7 +3038,9 @@ def deltaj_hijc_rebco(temperature):
                 deltaj_wst, x1, x2, 100e0
             )
             tmarg = current_sharing_t - thelium
-            jcrit0, b, t = sc.wstsc(current_sharing_t, bmax, strain, bc20m, tc0m)
+            jcrit0, _, _ = superconductors.wstsc(
+                current_sharing_t, bmax, strain, bc20m, tc0m
+            )
             if ml.variable_error(
                 current_sharing_t
             ):  # current sharing secant solver has failed.
@@ -3049,7 +3050,7 @@ def deltaj_hijc_rebco(temperature):
 
         # Temperature margin: An alternative method using secant solver
         elif isumat == 6:
-            current_sharing_t = sc.current_sharing_rebco(bmax, jsc)
+            current_sharing_t = superconductors.current_sharing_rebco(bmax, jsc)
             tmarg = current_sharing_t - thelium
             tfv.temp_margin = tmarg
 
@@ -3063,7 +3064,9 @@ def deltaj_hijc_rebco(temperature):
                 deltaj_gl_nbti, x1, x2, 100e0
             )
             tmarg = current_sharing_t - thelium
-            jcrit0, b, t = sc.gl_nbti(current_sharing_t, bmax, strain, bc20m, tc0m)
+            jcrit0, _, _ = superconductors.gl_nbti(
+                current_sharing_t, bmax, strain, bc20m, tc0m
+            )
             if ml.variable_error(
                 current_sharing_t
             ):  # current sharing secant solver has failed.
@@ -3081,7 +3084,9 @@ def deltaj_hijc_rebco(temperature):
                 deltaj_gl_rebco, x1, x2, 100e0
             )
             tmarg = current_sharing_t - thelium
-            jcrit0, b, t = sc.gl_rebco(current_sharing_t, bmax, strain, bc20m, tc0m)
+            jcrit0, _, _ = superconductors.gl_rebco(
+                current_sharing_t, bmax, strain, bc20m, tc0m
+            )
             if ml.variable_error(
                 current_sharing_t
             ):  # current sharing secant solver has failed.
@@ -3098,7 +3103,9 @@ def deltaj_hijc_rebco(temperature):
                 deltaj_hijc_rebco, x1, x2, 100e0
             )
             tmarg = current_sharing_t - thelium
-            jcrit0, b, t = sc.hijc_rebco(current_sharing_t, bmax, strain, bc20m, tc0m)
+            jcrit0, _, _ = superconductors.hijc_rebco(
+                current_sharing_t, bmax, strain, bc20m, tc0m
+            )
             if ml.variable_error(
                 current_sharing_t
             ):  # current sharing secant solver has failed.