From 319965bc6348b8aac3c9495a56d1dc6c67e68ba0 Mon Sep 17 00:00:00 2001 From: Michael Kovari Date: Fri, 12 Apr 2024 11:28:37 +0100 Subject: [PATCH 1/2] replaced the Newton-Raphson method in by scipy.optimize.newton --- process/maths_library.py | 104 ++++++------- process/pfcoil.py | 306 ++++++++----------------------------- process/sctfcoil.py | 142 ++++------------- process/superconductors.py | 147 +++++++++++++----- 4 files changed, 252 insertions(+), 447 deletions(-) diff --git a/process/maths_library.py b/process/maths_library.py index 34d7520dac..27e57fc703 100644 --- a/process/maths_library.py +++ b/process/maths_library.py @@ -4,66 +4,66 @@ Fortran; this is a temporary measure during the Python conversion process. """ -import numpy as np +# import numpy as np -def secant_solve(f, x1, x2, opt_tol=None): - """Solve an equation f(x)=0. +# def secant_solve(f, x1, x2, opt_tol=None): +# """Solve an equation f(x)=0. - Only requires one function evaluation per iteration (plus two to start with) - Does not require any derivatives. - https://en.wikipedia.org/wiki/Secant_method - Requires two initial values, x0 and x1, which should ideally be chosen to lie close to the root. +# Only requires one function evaluation per iteration (plus two to start with) +# Does not require any derivatives. +# https://en.wikipedia.org/wiki/Secant_method +# Requires two initial values, x0 and x1, which should ideally be chosen to lie close to the root. - TODO This is a duplicate of maths_library.f90:secant_solve(). The reason for - this is because Python functions can't be passed to Fortran subroutines via - the f2py interface, hence a Python function can't be passed into a Fortran - secant_solve(). Implementing a Python version (and temporarily leaving the - Fortran one) allows Python and Fortran functions to both use secant_solve() - as the Python conversion progresses. +# TODO This is a duplicate of maths_library.f90:secant_solve(). The reason for +# this is because Python functions can't be passed to Fortran subroutines via +# the f2py interface, hence a Python function can't be passed into a Fortran +# secant_solve(). Implementing a Python version (and temporarily leaving the +# Fortran one) allows Python and Fortran functions to both use secant_solve() +# as the Python conversion progresses. - :param f: function to solve - :type f: funnction, check - :param x1: first initial value of x - :type x1: float - :param x2: second initial value of x - :type x2: float - :param opt_tol: optional tolerance, defaults to None - :type opt_tol: float, optional - :return: the solution, error boolean and the residual - :rtype: tuple[float, bool, float] - """ - error = False - tol = 0.001 - if opt_tol: - tol = opt_tol +# :param f: function to solve +# :type f: funnction, check +# :param x1: first initial value of x +# :type x1: float +# :param x2: second initial value of x +# :type x2: float +# :param opt_tol: optional tolerance, defaults to None +# :type opt_tol: float, optional +# :return: the solution, error boolean and the residual +# :rtype: tuple[float, bool, float] +# """ +# error = False +# tol = 0.001 +# if opt_tol: +# tol = opt_tol - x = np.zeros(20) - x[0] = x1 - x[1] = x2 +# x = np.zeros(20) +# x[0] = x1 +# x[1] = x2 - # Calculate the first two values before the loop - fximinus1 = f(x[1]) - fximinus2 = f(x[0]) - for i in range(2, 20): - # Kludge to avoid zero in denominator - if fximinus1 == fximinus2: - fximinus1 = fximinus1 + (fximinus1 * 1e-6) +# # Calculate the first two values before the loop +# fximinus1 = f(x[1]) +# fximinus2 = f(x[0]) +# for i in range(2, 20): +# # Kludge to avoid zero in denominator +# if fximinus1 == fximinus2: +# fximinus1 = fximinus1 + (fximinus1 * 1e-6) - x[i] = x[i - 1] - fximinus1 * ((x[i - 1] - x[i - 2]) / (fximinus1 - fximinus2)) - # Store values for the *next iteration - fximinus2 = fximinus1 - fximinus1 = f(x[i]) - residual = fximinus1 +# x[i] = x[i - 1] - fximinus1 * ((x[i - 1] - x[i - 2]) / (fximinus1 - fximinus2)) +# # Store values for the *next iteration +# fximinus2 = fximinus1 +# fximinus1 = f(x[i]) +# residual = fximinus1 - # test for convergence - if abs(residual) < tol: - solution = x[i] - return solution, error, residual +# # test for convergence +# if abs(residual) < tol: +# solution = x[i] +# return solution, error, residual - # Convergence not achieved. Return the best value and a flag. - error = True - solution = x[i] - print(f"Secant solver not converged. {solution=} {residual=}") +# # Convergence not achieved. Return the best value and a flag. +# error = True +# solution = x[i] +# print(f"Secant solver not converged. {solution=} {residual=}") - return solution, error, residual +# return solution, error, residual diff --git a/process/pfcoil.py b/process/pfcoil.py index 63451837cb..fe136e7f11 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -13,7 +13,8 @@ from process.fortran import numerics 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 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 @@ -21,6 +22,10 @@ import numpy as np import numba import logging +from scipy import optimize + +# from process.sctfcoil import current_density_margin + logger = logging.getLogger(__name__) @@ -126,7 +131,7 @@ def pfcoil(self): # Set up call to MHD scaling routine for coil currents. # First break up Central Solenoid solenoid into 'filaments' - # Central Solenoid radius + # Central Solenoid mean radius pfv.rohc = bv.bore + 0.5e0 * bv.ohcth # nfxf is the total no of filaments into which the Central Solenoid is split, @@ -624,6 +629,7 @@ def pfcoil(self): # Allowable current density (for superconducting coils) for each coil, index i if pfv.ipfres == 0: bmax = max(abs(pfv.bpf[i]), abs(pf.bpf2[i])) + print("bmax at line 632 =", bmax) pfv.rjpfalw[i], jstrand, jsc, tmarg = self.superconpf( bmax, pfv.vf[i], @@ -982,12 +988,9 @@ def solv(self, ngrpmx, ngrp, nrws, gmat, bvec): return ccls, umat, vmat, sigma, work2 def ohcalc(self): - """Routine to perform calculations for the Central Solenoid solenoid. + """Routine to perform calculations for the Central Solenoid. author: P J Knight, CCFE, Culham Science Centre - This subroutine performs the calculations for the - Central Solenoid solenoid coil. - AEA FUS 251: A User's Guide to the PROCESS Systems Code """ hohc = bv.hmax * pfv.ohhghf @@ -1004,6 +1007,19 @@ def ohcalc(self): # Radius of inner edge pfv.ra[pfv.nohc - 1] = pfv.rb[pfv.nohc - 1] - bv.ohcth + print( + "bv.bore = ", + bv.bore, + "pfv.rohc = ", + pfv.rohc, + " pfv.ra[pfv.nohc - 1] = ", + pfv.ra[pfv.nohc - 1], + pfv.rohc, + "pfv.rb[pfv.nohc - 1] = ", + pfv.rb[pfv.nohc - 1], + "bv.ohcth = ", + bv.ohcth, + ) # Total cross-sectional area pfv.areaoh = 2.0e0 * hohc * bv.ohcth @@ -1063,6 +1079,16 @@ def ohcalc(self): pfv.rb[pfv.nohc - 1], hohc, ) + print( + "pfv.coheof", + pfv.coheof, + "pfv.ra[pfv.nohc - 1]", + pfv.ra[pfv.nohc - 1], + "pfv.rb[pfv.nohc - 1]", + pfv.rb[pfv.nohc - 1], + "hohc", + hohc, + ) # Peak field due to other PF coils plus plasma timepoint = 5 @@ -1182,7 +1208,7 @@ def ohcalc(self): if pfv.ipfres == 0: # Allowable coil overall current density at EOF # (superconducting coils only) - + print("bmaxoh at line 1192 =", pfv.bmaxoh) (jcritwp, pfv.jstrandoh_eof, pfv.jscoh_eof, tmarg1,) = self.superconpf( pfv.bmaxoh, pfv.vfohc, @@ -1199,6 +1225,7 @@ def ohcalc(self): pfv.rjohc = jcritwp * pfv.awpoh / pfv.areaoh # Allowable coil overall current density at BOP + print("bmaxoh0 at line 1208 =", pfv.bmaxoh0) (jcritwp, pfv.jstrandoh_bop, pfv.jscoh_bop, tmarg2,) = self.superconpf( pfv.bmaxoh0, @@ -2755,97 +2782,6 @@ def superconpf( :rtype: tuple[float, float, float, float] """ - """TODO maths_library.secant_solve() requires a function of one variable, - e.g. f(x). However, this function can require other variables as arguments - e.g. constants. Access to these variables (e.g. bmax, bc20m, tc0m) has - previously been provided through nested functions with implicit access - to the parent scope of superconpf() when the functions are passed to - secant_solve(). Now in Python, it might be better to explcitly pass - these variables as optional argument(s) to secant_solve() and remove - this nested function requirement. - """ - - def deltaj_nbti(temperature): - """Critical current density and current density difference in NbTi. - - :param temperature: temperature - :type temperature: float - :return: difference in current density - :rtype: float - """ - 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=}") - - deltaj_nbti = jcrit0 - jsc - return deltaj_nbti - - def deltaj_wst(temperature): - """Critical current density and current density difference for WST Nb3Sn. - - :param temperature: temperature - :type temperature: float - :return: difference in current density - :rtype: float - """ - jcrit0, _, _ = superconductors.wstsc(temperature, bmax, strain, bc20m, tc0m) - if ml.variable_error(jcrit0): # superconductors.wstsc has failed. - print(f"deltaj_wst: {bmax=} {temperature=} {jcrit0=}") - - deltaj_wst = jcrit0 - jsc - return deltaj_wst - - def deltaj_gl_nbti(temperature): - """Critical current density and current density difference in GL NbTi. - - :param temperature: temperature - :type temperature: float - :return: difference in current density - :rtype: float - """ - 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=}") - - deltaj_gl_nbti = jcrit0 - jsc - return deltaj_gl_nbti - - def deltaj_gl_rebco(temperature): - """Critical current density and current density difference in GL REBCO. - - :param temperature: temperature - :type temperature: float - :return: difference in current density - :rtype: float - """ - 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=}") - - deltaj_gl_rebco = jcrit0 - jsc - return deltaj_gl_rebco - - def deltaj_hijc_rebco(temperature): - """Critical current density and current density difference in high current density REBCO. - - :param temperature: temperature - :type temperature: float - :return: difference in current density - :rtype: float - """ - 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=}") - - deltaj_hijc_rebco = jcrit0 - jsc - return deltaj_hijc_rebco - # Find critical current density in superconducting strand, jcritstr if isumat == 1: # ITER Nb3Sn critical surface parameterization @@ -2970,155 +2906,37 @@ def deltaj_hijc_rebco(temperature): jsc = jstrand / (1.0e0 - fcu) # Temperature margin (already calculated in superconductors.bi2212 for isumat=2) - if (isumat == 1) or (isumat == 4): - # Newton-Raphson method; start at requested minimum temperature margin - ttest = thelium + tfv.tmargmin_cs - delt = 0.01e0 - jtol = 1.0e4 - - # Actual current density in superconductor, which should be equal to jcrit(thelium+tmarg) - # when we have found the desired value of tmarg - for lap in range(100): - if ttest <= 0.0: - eh.idiags[0] = lap - eh.fdiags[0] = ttest - eh.report_error(158) - break - - ttestm = ttest - delt - ttestp = ttest + delt - - if isumat in [1, 4]: - jcrit0, _, _ = superconductors.itersc( - ttest, bmax, strain, bc20m, tc0m - ) - if (abs(jsc - jcrit0) <= jtol) and ( - abs((jsc - jcrit0) / jsc) <= 0.01 - ): - break - - 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: - jcritp = jcritp + (jcritp * 1e-6) - ttest = ttest - 2.0e0 * delt * (jcrit0 - jsc) / (jcritp - jcritm) + if ( + (isumat == 1) + or (isumat == 3) + or (isumat == 4) + or (isumat == 5) + or (isumat == 7) + or (isumat == 8) + or (isumat == 9) + ): # Find temperature at which current density margin = 0 + # Temperature range ("bracketing interval") (K) + a = 4 + b = tc0m # tc0m = critical temperature (K) for superconductor at zero field and strain + + if isumat == 3: + args = (isumat, jsc, bmax, strain, bc20m, tc0m, c0) else: - eh.idiags[0] = lap - eh.fdiags[0] = ttest - eh.report_error(158) - - tmarg = ttest - thelium - - # MDK 13/7/18 Use secant solver for NbTi. - elif isumat == 3: - x1 = 4e0 # Initial values of temperature - x2 = 6e0 - # Solve for deltaj_nbti = 0 - current_sharing_t, error, residual = pml.secant_solve( - deltaj_nbti, x1, x2, 100e0 - ) - tmarg = current_sharing_t - thelium - jcrit0, _ = superconductors.jcrit_nbti( - current_sharing_t, bmax, c0, bc20m, tc0m - ) - if ml.variable_error( - current_sharing_t - ): # current sharing secant solver has failed. - print( - f"NbTi: {current_sharing_t=} {tmarg=} {jsc=} {jcrit0=} {residual=}" - ) - - # MDK 13/7/18 Use secant solver for WST. - elif isumat == 5: - # Current sharing temperature for WST Nb3Sn - x1 = 4e0 # Initial values of temperature - x2 = 6e0 - # Solve for deltaj_wst = 0 - current_sharing_t, error, residual = pml.secant_solve( - deltaj_wst, x1, x2, 100e0 - ) - tmarg = current_sharing_t - thelium - jcrit0, _, _ = superconductors.wstsc( - current_sharing_t, bmax, strain, bc20m, tc0m - ) - if ml.variable_error( - current_sharing_t - ): # current sharing secant solver has failed. - print( - f"WST: {current_sharing_t=} {tmarg=} {jsc=} {jcrit0=} {residual=}" - ) - - # Temperature margin: An alternative method using secant solver - elif isumat == 6: - current_sharing_t = superconductors.current_sharing_rebco(bmax, jsc) - tmarg = current_sharing_t - thelium - tfv.temp_margin = tmarg - - # SCM 16/03/20 Use secant solver for GL_nbti. - elif isumat == 7: - # Current sharing temperature for Durham Ginzburg-Landau Nb-Ti - x1 = 4.0e0 # Initial values of temperature - x2 = 6.0e0 - # Solve for deltaj_GL_nbti = 0 - current_sharing_t, error, residual = pml.secant_solve( - deltaj_gl_nbti, x1, x2, 100e0 - ) - tmarg = current_sharing_t - thelium - jcrit0, _, _ = superconductors.gl_nbti( - current_sharing_t, bmax, strain, bc20m, tc0m - ) - if ml.variable_error( - current_sharing_t - ): # current sharing secant solver has failed. - print( - f"Gl_nbti: {current_sharing_t=} {tmarg=} {jsc=} {jcrit0=} {residual=}" - ) - - # SCM 10/08/20 Use secant solver for superconductors.GL_REBCO. - elif isumat == 8: - # Current sharing temperature for Durham Ginzburg-Landau REBCO - x1 = 4.0e0 # Initial values of temperature - x2 = 6.0e0 - # Solve for deltaj_superconductors.GL_REBCO = 0 - current_sharing_t, error, residual = pml.secant_solve( - deltaj_gl_rebco, x1, x2, 100e0 - ) - tmarg = current_sharing_t - thelium - jcrit0, _, _ = superconductors.gl_rebco( - current_sharing_t, bmax, strain, bc20m, tc0m + args = (isumat, jsc, bmax, strain, bc20m, tc0m) + + t_zero_margin, root_result = optimize.brentq( + superconductors.current_density_margin, + a, + b, + args, + xtol=1e-4, + rtol=1e-4, + maxiter=100, + full_output=True, + disp=True, ) - if ml.variable_error( - current_sharing_t - ): # current sharing secant solver has failed. - print( - f"Gl_REBCO: {current_sharing_t=} {tmarg=} {jsc=} {jcrit0=} {residual=}" - ) - - elif isumat == 9: - # Current sharing temperature for Hazelton REBCO - x1 = 19.0e0 # Initial values of temperature - x2 = 21.0e0 - # Solve for deltaj_superconductors.HIJC_REBCO = 0 - current_sharing_t, error, residual = pml.secant_solve( - deltaj_hijc_rebco, x1, x2, 100e0 - ) - tmarg = current_sharing_t - thelium - jcrit0, _, _ = superconductors.hijc_rebco( - current_sharing_t, bmax, strain, bc20m, tc0m - ) - if ml.variable_error( - current_sharing_t - ): # current sharing secant solver has failed. - print( - f"HIJC_REBCO: {current_sharing_t=} {tmarg=} {jsc=} {jcrit0=} {residual=}" - ) + tmarg = t_zero_margin - thelium return jcritwp, jcritstr, jcritsc, tmarg diff --git a/process/sctfcoil.py b/process/sctfcoil.py index 521291ab3e..92a31da9df 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -21,6 +21,7 @@ import process.superconductors as superconductors from process.utilities.f2py_string_patch import f2py_compatible_to_string +from scipy import optimize logger = logging.getLogger(__name__) @@ -186,7 +187,7 @@ def supercon_croco(self, aturn, bmax, iop, thelium, output: bool): # when we have found the desired value of tmarg jsc = iooic * jcritsc - # Temperature margin using secant solver + # Temperature margin current_sharing_t = superconductors.current_sharing_rebco(bmax, jsc) tmarg = current_sharing_t - thelium tfcoil_variables.temp_margin = ( @@ -588,6 +589,8 @@ def supercon( 4 = ITER Nb3Sn, user-defined parameters 5 = WST Nb3Sn parameterisation 7 = Durham Ginzburg-Landau Nb-Ti parameterisation + 8 = Durham Ginzburg-Landau critical surface model for REBCO + 9 = Hazelton experimental data + Zhai conceptual model for REBCO fhts : input real : Adjustment factor (<= 1) to account for strain, radiation damage, fatigue or AC losses tdmptf : input real : Dump time (sec) @@ -789,6 +792,7 @@ def supercon( error_handling.report_error(266) # Temperature margin (already calculated in superconductors.bi2212 for isumat=2) + if ( (isumat == 1) or (isumat == 3) @@ -797,112 +801,32 @@ def supercon( or (isumat == 7) or (isumat == 8) or (isumat == 9) - ): - # Newton-Raphson method; start approx at requested minimum temperature margin - ttest = thelium + tfcoil_variables.tmargmin_tf + 0.001e0 - delt = 0.01e0 - jtol = 1.0e4 - - for lap in range(100): - if ttest <= delt: - error_handling.idiags[0] = lap - error_handling.fdiags[0] = ttest - error_handling.report_error(157) - break - - # Calculate derivative numerically - ttestm = ttest - delt - ttestp = ttest + delt - - # Issue #483 to be on the safe side, check the fractional as well as the absolute error - if isumat in (1, 4): - jcrit0, _, _ = superconductors.itersc( - ttest, bmax, strain, bc20m, tc0m - ) - if (abs(jsc - jcrit0) <= jtol) and ( - abs((jsc - jcrit0) / jsc) <= 0.01 - ): - break - jcritm, _, _ = superconductors.itersc( - ttestm, bmax, strain, bc20m, tc0m - ) - jcritp, _, _ = superconductors.itersc( - ttestp, bmax, strain, bc20m, tc0m - ) - elif isumat == 3: - jcrit0, _ = superconductors.jcrit_nbti(ttest, bmax, c0, bc20m, tc0m) - if (abs(jsc - jcrit0) <= jtol) and ( - abs((jsc - jcrit0) / jsc) <= 0.01 - ): - break - jcritm, _ = super.jcrit_nbti(ttestm, bmax, c0, bc20m, tc0m) - jcritp, _ = superconductors.jcrit_nbti( - ttestp, bmax, c0, bc20m, tc0m - ) - elif isumat == 5: - jcrit0, _, _ = superconductors.wstsc( - ttest, bmax, strain, bc20m, tc0m - ) - if (abs(jsc - jcrit0) <= jtol) and ( - abs((jsc - jcrit0) / jsc) <= 0.01 - ): - break - jcritm, _, _ = superconductors.wstsc( - ttestm, bmax, strain, bc20m, tc0m - ) - jcritp, _, _ = superconductors.wstsc( - ttestp, bmax, strain, bc20m, tc0m - ) - elif isumat == 7: - jcrit0, _, _ = superconductors.gl_nbti( - ttest, bmax, strain, bc20m, tc0m - ) - if (abs(jsc - jcrit0) <= jtol) and ( - abs((jsc - jcrit0) / jsc) <= 0.01 - ): - break - jcritm, _, _ = superconductors.gl_nbti( - ttestm, bmax, strain, bc20m, tc0m - ) - jcritp, _, _ = superconductors.gl_nbti( - ttestp, bmax, strain, bc20m, tc0m - ) - elif isumat == 8: - jcrit0, _, _ = superconductors.gl_rebco( - ttest, bmax, strain, bc20m, tc0m - ) - if (abs(jsc - jcrit0) <= jtol) and ( - abs((jsc - jcrit0) / jsc) <= 0.01 - ): - break - jcritm, _, _ = superconductors.gl_rebco( - ttestm, bmax, strain, bc20m, tc0m - ) - jcritp, _, _ = superconductors.gl_rebco( - ttestp, bmax, strain, bc20m, tc0m - ) - elif isumat == 9: - jcrit0, _, _ = superconductors.hijc_rebco( - ttest, bmax, strain, bc20m, tc0m - ) - if (abs(jsc - jcrit0) <= jtol) and ( - abs((jsc - jcrit0) / jsc) <= 0.01 - ): - break - jcritm, _, _ = superconductors.hijc_rebco( - ttestm, bmax, strain, bc20m, tc0m - ) - jcritp, _, _ = superconductors.hijc_rebco( - ttestp, bmax, strain, bc20m, tc0m - ) + ): # Find temperature at which current density margin = 0 + # Temperature range ("bracketing interval") (K) + # Setting the upper limit of the interval slightly lower than tc0m ensures + # the reduced temperature < 1 after strain is taken into accoun + a = 4 + b = ( + 0.94 * tc0m + ) # tc0m = critical temperature (K) for superconductor at zero field and strain - ttest = ttest - 2.0e0 * delt * (jcrit0 - jsc) / (jcritp - jcritm) + if isumat == 3: + args = (isumat, jsc, bmax, strain, bc20m, tc0m, c0) else: - error_handling.idiags[0] = lap - error_handling.fdiags[0] = ttest - error_handling.report_error(157) - - tmarg = ttest - thelium + args = (isumat, jsc, bmax, strain, bc20m, tc0m) + + t_zero_margin, root_result = optimize.brentq( + superconductors.current_density_margin, + a, + b, + args, + xtol=1e-4, + rtol=1e-4, + maxiter=100, + full_output=True, + disp=True, + ) + tmarg = t_zero_margin - thelium tfcoil_variables.temp_margin = tmarg # Find the current density limited by the protection limit @@ -914,17 +838,13 @@ def supercon( ) if output: # Output -------------------------- - if ttest <= 0.0e0: + if tmarg <= 0.0e0: logger.warning( """Negative TFC temperature margin - ttest: {ttest} + tmarg: {tmarg} bmax: {bmax} jcrit0: {jcrit0} jsc: {jsc} - ttestp: {ttestp} - ttestm: {ttestm} - jcritp: {jcritp} - jcritm: {jcritm} """ ) diff --git a/process/superconductors.py b/process/superconductors.py index 1e0efdf9da..2be52f2d1d 100644 --- a/process/superconductors.py +++ b/process/superconductors.py @@ -2,7 +2,8 @@ import numpy as np from process.fortran import error_handling as eh, rebco_variables -from process.maths_library import secant_solve + +from scipy import optimize logger = logging.getLogger(__name__) @@ -77,7 +78,29 @@ def deltaj_rebco(temperature): jcritical, _ = jcrit_rebco(temperature, bfield) return jcritical - j - current_sharing_t, _, _ = secant_solve(deltaj_rebco, 4, 20, 1e7) + # Replace secant_solve by scipy.optimize.brentq. + # Find temperature at which current density margin = 0 + # Temperature range ("bracketing interval") (K) + # For brentq to work, the current density margin must be + # negative at the upper temperature of the interval. + + a = 4 # Bottom of temperature interval + b = 90 # Top of temperature interval + + # No additional arguments are required for deltaj_rebco since it only has one argument. + current_sharing_t, root_result = optimize.brentq( + deltaj_rebco, + a, + b, + xtol=1e-4, + rtol=1e-4, + maxiter=100, + full_output=True, + disp=True, + ) + + print(root_result) + return current_sharing_t @@ -126,45 +149,43 @@ def itersc(thelium, bmax, strain, bc20max, tc0max): strfun = strfun * ca1 - ca2 * strain strfun = 1.0 + (1 / (1.0 - ca1 * eps0a)) * strfun - # $B^*_{C2} (0,\epsilon)$ + # Critical field at zero temperature and current, corrected for strain bc20eps = bc20max * strfun - # $T^*_C (0,\epsilon)$ + # Critical temperature at zero field and current, corrected for strain tc0eps = tc0max * strfun ** (1 / 3) - # Reduced temperature, restricted to be < 1 - # Should remain < 1 for thelium < 0.94*tc0max (i.e. 15 kelvin for i_tf_sc_mat=1) - if thelium / tc0eps >= 1.0: eh.fdiags[0] = thelium eh.fdiags[1] = tc0eps eh.report_error(159) - t = min(thelium / tc0eps, 0.9999) - - # Reduced magnetic field at zero temperature - # Should remain < 1 for bmax < 0.83*bc20max (i.e. 27 tesla for i_tf_sc_mat=1) + # Reduced temperature at zero field, corrected for strain + # I think t > 0 unless something has gone badly wrong. + t = thelium / tc0eps if bmax / bc20eps >= 1.0: eh.fdiags[0] = bmax eh.fdiags[1] = bc20eps eh.report_error(160) - bzero = min(bmax / bc20eps, 0.9999) + # Reduced field at zero temperature, taking account of strain + bzero = bmax / bc20eps - # Critical temperature (K) - tcrit = tc0eps * (1.0 - bzero) ** (1 / 1.52) # bzero must be < 1 to avoid NaNs + # Critical temperature at given strain and field + # tcrit is not used in the calculation of jcrit. + if bzero < 1.0: + tcrit = tc0eps * (1.0 - bzero) ** (1 / 1.52) + else: + tcrit = -tc0eps * abs(1.0 - bzero) ** ( + 1 / 1.52 + ) # Prevents NaNs. Sets tcrit negative - # Critical field (T) + # Critical field (T) at given strain and temperature bcrit = bc20eps * (1.0 - t**1.52) - # Reduced magnetic field, restricted to be < 1 - if bmax / bcrit >= 1.0: - eh.fdiags[0] = bmax - eh.fdiags[1] = bcrit - eh.report_error(161) - - bred = min(bmax / bcrit, 0.9999) + # Reduced field at given strain and temperature + bred = bmax / bcrit # Critical current density in superconductor (A/m2) # ITER parameterization is for the current in a single strand, @@ -173,10 +194,14 @@ def itersc(thelium, bmax, strain, bc20max, tc0max): scalefac = np.pi * (0.5 * diter) ** 2 * (1.0 - cuiter) jc1 = (csc / bmax) * strfun - jc2 = (1.0 - t**1.52) * (1.0 - t**2) # t must be < 1 to avoid NaNs - jc3 = bred**p * (1.0 - bred) ** q # bred must be < 1 to avoid NaNs + jc2 = (1.0 - t**1.52) * (1.0 - t**2) - jcrit = jc1 * jc2 * jc3 / scalefac + if (bred > 0) and (bred < 1): # Normal case + jc3 = bred**p * (1.0 - bred) ** q + jcrit = jc1 * jc2 * jc3 / scalefac + else: + jc3 = abs(bred) ** p * abs(1.0 - bred) ** q # Prevents NaNs + jcrit = -abs(jc1 * jc2 * jc3 / scalefac) # Sets jcrit negative return jcrit, bcrit, tcrit @@ -192,32 +217,26 @@ def jcrit_nbti(temperature, bmax, c0, bc20max, tc0max): tc0max : input real : Critical temperature (K) at zero field and strain jcrit : output real : Critical current density in superconductor (A/m2) tcrit : output real : Critical temperature (K) - This routine calculates the critical current density and - temperature in superconducting TF coils using NbTi - as the superconductor. + This routine calculates the critical current density and temperature in + superconducting TF coils using NbTi superconductor. """ bratio = bmax / bc20max - if bmax < bc20max: - # Critical temperature (K) + if bratio < 1: + # Critical temperature (K) at field tcrit = tc0max * (1.0 - bratio) ** 0.59 else: - # Allow bmax > bc20max but set error flag - # Fudge to give real (negative) value if bratio < 1 + # Allow bmax > bc20max + # Fudge to give real (negative) value if bratio > 1 tcrit = tc0max * (1.0 - bratio) - # Allow tbar to be negative but set error flag + # If temperature > tcrit then the critical surface has been exceeded and tbar is negative tbar = 1.0 - temperature / tcrit # Critical current density (A/m2) jcrit = c0 * (1.0 - bratio) * tbar - # if ((temperature > tcrit).or.(bmax > bc20max))then - # write(*,*)'jcrit_nbti: out of range: ', 'bmax =', bmax, ' bc20max =', bc20max, & - # ' temperature =',temperature, ' tcrit =',tcrit - # end if - return jcrit, tcrit @@ -445,13 +464,16 @@ def hijc_rebco(thelium, bmax, strain, bc20max, t_c0): described in: M. J. Wolf, N. Bagrets, W. H. Fietz, C. Lange and K. Weiss, "Critical Current Densities of 482 A/mm2 in HTS CrossConductors at 4.2 K and 12 T," in IEEE Transactions on Applied Superconductivity, vol. 28, no. 4, pp. 1-4, - June 2018, Art no. 4802404, doi: 10.1109/TASC.2018.2815767. And on the experimental + June 2018, Art no. 4802404, doi: 10.1109/TASC.2018.2815767. + + And on the experimental data presented here: "2G HTS Wire Development at SuperPower", Drew W. Hazelton, February 16, 2017 https://indico.cern.ch/event/588810/contributions/2473740/ The high Ic parameterization is a result of modifications based on Ic values observed in: "Conceptual design of HTS magnets for fusion nuclear science facility", Yuhu Zhai, Danko van der Laan, Patrick Connolly, Charles Kessel, 2021, https://doi.org/10.1016/j.fusengdes.2021.112611 + The parameter A is transformed into a function A(T) based on a Newton polynomial fit considering A(4.2 K) = 2.2e8, A(20 K) = 2.3e8 and A(65 K) = 3.5e8. These values were selected manually. A good fit to the pubished data can be seen in the 4-10 T range @@ -487,8 +509,21 @@ def hijc_rebco(thelium, bmax, strain, bc20max, t_c0): A_t = A_0 + (u * thelium**2) + (v * thelium) # Critical current density (A/m2) + # In the original formula bcrit must be > bmax to prevent NaNs. + # However, negative jcrit is permissible (I think). + # So when bcrit < bmax, I reverse the sign of the bracket, + # giving a negative but real value of jcrit. + + if bcrit > bmax: + jcrit = ( + (A_t / bmax) * bcrit**b * (bmax / bcrit) ** p * (1 - bmax / bcrit) ** q + ) + else: + jcrit = ( + (A_t / bmax) * bcrit**b * (bmax / bcrit) ** p * (bmax / bcrit - 1) ** q + ) - jcrit = (A_t / bmax) * bcrit**b * (bmax / bcrit) ** p * (1 - bmax / bcrit) ** q + # print("thelium = ", thelium, " bcrit = ", bcrit, " bmax = ", bmax, " 1 - bmax / bcrit = ", 1 - bmax / bcrit) # Jc times HTS area: default area is width 4mm times HTS layer thickness 1 um, # divided by the tape area to provide engineering Jc per tape,! @@ -710,3 +745,35 @@ def croco(jcritsc, conductor_area, croco_od, croco_thick): conductor_rebco_fraction, conductor_critical_current, ) + + +def current_density_margin(ttest, isumat, jsc, bmax, strain, bc20m, tc0m, c0=None): + """Current density margin is the difference between the operating current density and + the critical current density of a superconductor, at a given temperature and field. + It is zero at the current-sharing temperature. + ttest : input real : Temperature + isumat : input real : Switch for superconductor material + (This has different global names depending on which coil is referred to.) + jsc : input real : actual current density + bmax : input real : magnetic field (T) + strain : input real : superconductor strain + bc20m, tc0m : input real : superconductor parameters + """ + + # 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, _, _ = wstsc(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(ttest, bmax, strain, bc20m, tc0m) + + return jcrit - jsc From 45e01879cac2c85b0a053bc5d2caf47a0cc1387d Mon Sep 17 00:00:00 2001 From: Michael Kovari Date: Tue, 16 Apr 2024 15:31:03 +0100 Subject: [PATCH 2/2] New function Bottura_scaling to simplify supercon routines. --- process/maths_library.py | 64 -------- process/pfcoil.py | 58 ++----- process/sctfcoil.py | 30 ++-- process/superconductors.py | 228 ++++++++++----------------- tests/integration/test_pfcoil_int.py | 26 +-- tests/unit/test_sctfcoil.py | 12 +- tests/unit/test_superconductors.py | 6 +- 7 files changed, 133 insertions(+), 291 deletions(-) diff --git a/process/maths_library.py b/process/maths_library.py index 27e57fc703..d21952ac38 100644 --- a/process/maths_library.py +++ b/process/maths_library.py @@ -3,67 +3,3 @@ It is possible that some functions here are duplicated in both Python and Fortran; this is a temporary measure during the Python conversion process. """ - -# import numpy as np - - -# def secant_solve(f, x1, x2, opt_tol=None): -# """Solve an equation f(x)=0. - -# Only requires one function evaluation per iteration (plus two to start with) -# Does not require any derivatives. -# https://en.wikipedia.org/wiki/Secant_method -# Requires two initial values, x0 and x1, which should ideally be chosen to lie close to the root. - -# TODO This is a duplicate of maths_library.f90:secant_solve(). The reason for -# this is because Python functions can't be passed to Fortran subroutines via -# the f2py interface, hence a Python function can't be passed into a Fortran -# secant_solve(). Implementing a Python version (and temporarily leaving the -# Fortran one) allows Python and Fortran functions to both use secant_solve() -# as the Python conversion progresses. - -# :param f: function to solve -# :type f: funnction, check -# :param x1: first initial value of x -# :type x1: float -# :param x2: second initial value of x -# :type x2: float -# :param opt_tol: optional tolerance, defaults to None -# :type opt_tol: float, optional -# :return: the solution, error boolean and the residual -# :rtype: tuple[float, bool, float] -# """ -# error = False -# tol = 0.001 -# if opt_tol: -# tol = opt_tol - -# x = np.zeros(20) -# x[0] = x1 -# x[1] = x2 - -# # Calculate the first two values before the loop -# fximinus1 = f(x[1]) -# fximinus2 = f(x[0]) -# for i in range(2, 20): -# # Kludge to avoid zero in denominator -# if fximinus1 == fximinus2: -# fximinus1 = fximinus1 + (fximinus1 * 1e-6) - -# x[i] = x[i - 1] - fximinus1 * ((x[i - 1] - x[i - 2]) / (fximinus1 - fximinus2)) -# # Store values for the *next iteration -# fximinus2 = fximinus1 -# fximinus1 = f(x[i]) -# residual = fximinus1 - -# # test for convergence -# if abs(residual) < tol: -# solution = x[i] -# return solution, error, residual - -# # Convergence not achieved. Return the best value and a flag. -# error = True -# solution = x[i] -# print(f"Secant solver not converged. {solution=} {residual=}") - -# return solution, error, residual diff --git a/process/pfcoil.py b/process/pfcoil.py index fe136e7f11..6c7b0761bc 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -14,7 +14,6 @@ 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 @@ -24,9 +23,6 @@ import logging from scipy import optimize -# from process.sctfcoil import current_density_margin - - logger = logging.getLogger(__name__) RMU0 = ft.constants.rmu0 @@ -629,7 +625,7 @@ def pfcoil(self): # Allowable current density (for superconducting coils) for each coil, index i if pfv.ipfres == 0: bmax = max(abs(pfv.bpf[i]), abs(pf.bpf2[i])) - print("bmax at line 632 =", bmax) + pfv.rjpfalw[i], jstrand, jsc, tmarg = self.superconpf( bmax, pfv.vf[i], @@ -1007,19 +1003,6 @@ def ohcalc(self): # Radius of inner edge pfv.ra[pfv.nohc - 1] = pfv.rb[pfv.nohc - 1] - bv.ohcth - print( - "bv.bore = ", - bv.bore, - "pfv.rohc = ", - pfv.rohc, - " pfv.ra[pfv.nohc - 1] = ", - pfv.ra[pfv.nohc - 1], - pfv.rohc, - "pfv.rb[pfv.nohc - 1] = ", - pfv.rb[pfv.nohc - 1], - "bv.ohcth = ", - bv.ohcth, - ) # Total cross-sectional area pfv.areaoh = 2.0e0 * hohc * bv.ohcth @@ -1079,16 +1062,6 @@ def ohcalc(self): pfv.rb[pfv.nohc - 1], hohc, ) - print( - "pfv.coheof", - pfv.coheof, - "pfv.ra[pfv.nohc - 1]", - pfv.ra[pfv.nohc - 1], - "pfv.rb[pfv.nohc - 1]", - pfv.rb[pfv.nohc - 1], - "hohc", - hohc, - ) # Peak field due to other PF coils plus plasma timepoint = 5 @@ -1208,7 +1181,7 @@ def ohcalc(self): if pfv.ipfres == 0: # Allowable coil overall current density at EOF # (superconducting coils only) - print("bmaxoh at line 1192 =", pfv.bmaxoh) + (jcritwp, pfv.jstrandoh_eof, pfv.jscoh_eof, tmarg1,) = self.superconpf( pfv.bmaxoh, pfv.vfohc, @@ -1225,7 +1198,6 @@ def ohcalc(self): pfv.rjohc = jcritwp * pfv.awpoh / pfv.areaoh # Allowable coil overall current density at BOP - print("bmaxoh0 at line 1208 =", pfv.bmaxoh0) (jcritwp, pfv.jstrandoh_bop, pfv.jscoh_bop, tmarg2,) = self.superconpf( pfv.bmaxoh0, @@ -2916,25 +2888,25 @@ def superconpf( or (isumat == 8) or (isumat == 9) ): # Find temperature at which current density margin = 0 - # Temperature range ("bracketing interval") (K) - a = 4 - b = tc0m # tc0m = critical temperature (K) for superconductor at zero field and strain if isumat == 3: - args = (isumat, jsc, bmax, strain, bc20m, tc0m, c0) + arguments = (isumat, jsc, bmax, strain, bc20m, tc0m, c0) else: - args = (isumat, jsc, bmax, strain, bc20m, tc0m) + arguments = (isumat, jsc, bmax, strain, bc20m, tc0m) - t_zero_margin, root_result = optimize.brentq( + another_estimate = 2 * thelium + t_zero_margin, root_result = optimize.newton( superconductors.current_density_margin, - a, - b, - args, - xtol=1e-4, - rtol=1e-4, - maxiter=100, + thelium, + fprime=None, + args=arguments, + tol=1.0e-06, + maxiter=50, + fprime2=None, + x1=another_estimate, + rtol=1.0e-6, full_output=True, - disp=True, + disp=False, ) tmarg = t_zero_margin - thelium diff --git a/process/sctfcoil.py b/process/sctfcoil.py index 92a31da9df..d29f8f054b 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -802,30 +802,28 @@ def supercon( or (isumat == 8) or (isumat == 9) ): # Find temperature at which current density margin = 0 - # Temperature range ("bracketing interval") (K) - # Setting the upper limit of the interval slightly lower than tc0m ensures - # the reduced temperature < 1 after strain is taken into accoun - a = 4 - b = ( - 0.94 * tc0m - ) # tc0m = critical temperature (K) for superconductor at zero field and strain if isumat == 3: - args = (isumat, jsc, bmax, strain, bc20m, tc0m, c0) + arguments = (isumat, jsc, bmax, strain, bc20m, tc0m, c0) else: - args = (isumat, jsc, bmax, strain, bc20m, tc0m) + arguments = (isumat, jsc, bmax, strain, bc20m, tc0m) - t_zero_margin, root_result = optimize.brentq( + another_estimate = 2 * thelium + t_zero_margin, root_result = optimize.newton( superconductors.current_density_margin, - a, - b, - args, - xtol=1e-4, - rtol=1e-4, - maxiter=100, + thelium, + fprime=None, + args=arguments, + # args=(isumat, jsc, bmax, strain, bc20m, tc0m,), + tol=1.0e-06, + maxiter=50, + fprime2=None, + x1=another_estimate, + rtol=1.0e-6, full_output=True, disp=True, ) + # print(root_result) # Diagnostic for newton method tmarg = t_zero_margin - thelium tfcoil_variables.temp_margin = tmarg diff --git a/process/superconductors.py b/process/superconductors.py index 2be52f2d1d..68db6bb073 100644 --- a/process/superconductors.py +++ b/process/superconductors.py @@ -78,37 +78,29 @@ def deltaj_rebco(temperature): jcritical, _ = jcrit_rebco(temperature, bfield) return jcritical - j - # Replace secant_solve by scipy.optimize.brentq. - # Find temperature at which current density margin = 0 - # Temperature range ("bracketing interval") (K) - # For brentq to work, the current density margin must be - # negative at the upper temperature of the interval. - - a = 4 # Bottom of temperature interval - b = 90 # Top of temperature interval - # No additional arguments are required for deltaj_rebco since it only has one argument. - current_sharing_t, root_result = optimize.brentq( + + estimate = 10.0 + another_estimate = 20.0 + current_sharing_t, root_result = optimize.newton( deltaj_rebco, - a, - b, - xtol=1e-4, - rtol=1e-4, - maxiter=100, + estimate, + tol=1e-6, + rtol=1e-6, + maxiter=50, + x1=another_estimate, full_output=True, disp=True, ) - print(root_result) - return current_sharing_t -def itersc(thelium, bmax, strain, bc20max, tc0max): +def itersc(temperature, bmax, strain, bc20max, tc0max): """Implementation of ITER Nb3Sn critical surface implementation author: R Kemp, CCFE, Culham Science Centre author: P J Knight, CCFE, Culham Science Centre - thelium : input real : Coolant/SC temperature (K) + temperature : input real : superconductor temperature (K) bmax : input real : Magnetic field at conductor (T) strain : input real : Strain in superconductor bc20max : input real : Upper critical field (T) for superconductor @@ -117,19 +109,23 @@ def itersc(thelium, bmax, strain, bc20max, tc0max): jcrit : output real : Critical current density in superconductor (A/m2) bcrit : output real : Critical field (T) tcrit : output real : Critical temperature (K) - This routine calculates the critical current density and + Critical current density and temperature in the superconducting TF coils using the ITER Nb3Sn critical surface model. $J_C(B,T,\\epsilon)$ Parameterization for ITER Nb3Sn production, L. Bottura, CERN-ITER Collaboration Report, Version 2, April 2nd 2008 - (distributed by Arnaud Devred, ITER, 10th April 2008) + (distributed by Arnaud Devred, ITER, 10th April 2008), now published as + Jc(B,T,epsilon) Parameterization for the ITER Nb3Sn Production, + Luca Bottura and Bernardo Bordini, + IEEE TRANSACTIONS ON APPLIED SUPERCONDUCTIVITY, VOL. 19, NO. 3, JUNE 2009. + ITER Nb3Sn critical surface parameterization (2MMF7J) (2008), https://user.iter.org/?uid=2MMF7J&action=get_document ITER DDD 11-7: Magnets - conductors (2NBKXY) (2009), https://user.iter.org/?uid=2NBKXY&action=get_document """ - csc = 19922.0e6 # scaling constant C + csc = 19922.0 # scaling constant C [AT/mm2] p = 0.63 # low field exponent p q = 2.1 # high field exponent q ca1 = 44.48 # strain fitting constant C_{a1} @@ -138,70 +134,17 @@ def itersc(thelium, bmax, strain, bc20max, tc0max): diter = 0.82 # ITER strand diameter (mm) cuiter = 0.5 # ITER strand copper fraction - # $\epsilon_{sh}$ - epssh = (ca2 * eps0a) / (np.sqrt(ca1**2 - ca2**2)) - - # Strain function $s(\epsilon)$ - # 0.83 < s < 1.0, for -0.005 < strain < 0.005 - strfun = np.sqrt(epssh**2 + eps0a**2) - np.sqrt( - (strain - epssh) ** 2 + eps0a**2 + jscaling, bcrit, tcrit = Bottura_scaling( + csc, p, q, ca1, ca2, eps0a, temperature, bmax, strain, bc20max, tc0max ) - strfun = strfun * ca1 - ca2 * strain - strfun = 1.0 + (1 / (1.0 - ca1 * eps0a)) * strfun - - # Critical field at zero temperature and current, corrected for strain - bc20eps = bc20max * strfun - - # Critical temperature at zero field and current, corrected for strain - tc0eps = tc0max * strfun ** (1 / 3) - - if thelium / tc0eps >= 1.0: - eh.fdiags[0] = thelium - eh.fdiags[1] = tc0eps - eh.report_error(159) - - # Reduced temperature at zero field, corrected for strain - # I think t > 0 unless something has gone badly wrong. - t = thelium / tc0eps - - if bmax / bc20eps >= 1.0: - eh.fdiags[0] = bmax - eh.fdiags[1] = bc20eps - eh.report_error(160) - - # Reduced field at zero temperature, taking account of strain - bzero = bmax / bc20eps - - # Critical temperature at given strain and field - # tcrit is not used in the calculation of jcrit. - if bzero < 1.0: - tcrit = tc0eps * (1.0 - bzero) ** (1 / 1.52) - else: - tcrit = -tc0eps * abs(1.0 - bzero) ** ( - 1 / 1.52 - ) # Prevents NaNs. Sets tcrit negative - - # Critical field (T) at given strain and temperature - bcrit = bc20eps * (1.0 - t**1.52) - - # Reduced field at given strain and temperature - bred = bmax / bcrit # Critical current density in superconductor (A/m2) - # ITER parameterization is for the current in a single strand, + # ITER parameters are for the current in a single strand, # not per unit area, so scalefac converts to current density + # Convert from mm2 to m2. + scalefac = np.pi * (0.5 * diter) ** 2 * (1.0 - cuiter) / 1.0e6 - scalefac = np.pi * (0.5 * diter) ** 2 * (1.0 - cuiter) - - jc1 = (csc / bmax) * strfun - jc2 = (1.0 - t**1.52) * (1.0 - t**2) - - if (bred > 0) and (bred < 1): # Normal case - jc3 = bred**p * (1.0 - bred) ** q - jcrit = jc1 * jc2 * jc3 / scalefac - else: - jc3 = abs(bred) ** p * abs(1.0 - bred) ** q # Prevents NaNs - jcrit = -abs(jc1 * jc2 * jc3 / scalefac) # Sets jcrit negative + jcrit = jscaling / scalefac return jcrit, bcrit, tcrit @@ -570,102 +513,93 @@ def wstsc(temperature, bmax, strain, bc20max, tc0max): # epsilon_{0,a} eps0a = 0.00312 - # $\epsilon_{sh}$ + jscaling, bcrit, tcrit = Bottura_scaling( + csc, p, q, ca1, ca2, eps0a, temperature, bmax, strain, bc20max, tc0max + ) + + # scale from mm2 to m2 + scalefac = 1.0e6 + jcrit = jscaling * scalefac + return jcrit, bcrit, tcrit + + +def Bottura_scaling( + csc, p, q, ca1, ca2, eps0a, temperature, bmax, strain, bc20max, tc0max +): + """ + This implements the scaling from + Jc(B,T,epsilon) Parameterization for the ITER Nb3Sn Production, + Luca Bottura and Bernardo Bordini, + IEEE TRANSACTIONS ON APPLIED SUPERCONDUCTIVITY, VOL. 19, NO. 3, JUNE 2009. + + The parameters and scale factors vary from one wire type to another. + """ + epssh = (ca2 * eps0a) / (np.sqrt(ca1**2 - ca2**2)) - # Strain function $s(\epsilon)$ - # 0.83 < s < 1.0, for -0.005 < strain < 0.005 + # Strain function + # 0.83 < s < 1.0, for -0.005 < strain < 0.005 strfun = np.sqrt(epssh**2 + eps0a**2) - np.sqrt( (strain - epssh) ** 2 + eps0a**2 ) strfun = strfun * ca1 - ca2 * strain - strfun = 1.0 + (1.0 / (1.0 - ca1 * eps0a)) * strfun - if strfun < 0: - logger.warning(f"Strain function < 0: {strfun} for {strain = }") + strfun = 1.0 + (1 / (1.0 - ca1 * eps0a)) * strfun - # $B^*_{C2} (0,\epsilon)$ + # Critical field at zero temperature and current, corrected for strain bc20eps = bc20max * strfun - # $T^*_C (0,\epsilon)$ - tc0eps = tc0max * strfun ** (1.0 / 3.0) - - # Reduced temperature - # Should remain < 1 for temperature < 0.94*tc0max (i.e. 15 kelvin for i_tf_sc_mat=1) + # Critical temperature at zero field and current, corrected for strain + tc0eps = tc0max * strfun ** (1 / 3) if temperature / tc0eps >= 1.0: eh.fdiags[0] = temperature eh.fdiags[1] = tc0eps eh.report_error(159) - # t = min(thelium/tc0eps, 0.9999D0) + # Reduced temperature at zero field, corrected for strain + # t > 1 is permitted, indicating the temperature is above the critical value at zero field. t = temperature / tc0eps - # Reduced magnetic field at zero temperature - # Should remain < 1 for bmax < 0.83*bc20max (i.e. 27 tesla for i_tf_sc_mat=1) - if bmax / bc20eps >= 1.0: eh.fdiags[0] = bmax eh.fdiags[1] = bc20eps eh.report_error(160) - # bzero = min(bmax/bc20eps, 0.9999D0) + # Reduced field at zero temperature, taking account of strain bzero = bmax / bc20eps - if bzero < 1.0: - # Critical temperature (K) - tcrit = tc0eps * (1.0 - bzero) ** (1.0 / 1.52) - else: - # Allow bzero > 1, fudge to give real (negative) value of tcrit - # This generates a real (negative) and continuous (but not differentiable) - # function of bzero. - tcrit = tc0eps - - # Critical field (T). Negative if normalised temperature t>1 - if t > 0.0: - bcrit = bc20eps * (1.0 - t**1.52) - else: - # Allow t<0, fudge to give real value of bcrit - bcrit = bc20eps * (1.0 - t) - - # Reduced magnetic field, restricted to be < 1 - if bmax / bcrit >= 1.0: - eh.fdiags[0] = bmax - eh.fdiags[1] = bcrit - eh.report_error(161) - - # bred = min(bmax/bcrit, 0.9999D0) - bred = bmax / bcrit + # Critical temperature at given strain and field + # tcrit is not used in the calculation of jcrit. + if bzero < 1.0: # Normal case, field is within critical surface + tcrit = tc0eps * (1.0 - bzero) ** (1 / 1.52) + else: # Abnormal case, field is too high. + tcrit = -tc0eps * abs(1.0 - bzero) ** ( + 1 / 1.52 + ) # Prevents NaNs. Sets tcrit negative - if (bred > 0.0) and (bred < 1.0): - jc3 = bred**p * (1.0 - bred) ** q # bred must be < 1 to avoid NaNs - else: - # Allow bred > 1 or <0, fudge to give real (negative) value of jc3 - # This generates a real (negative) and continuous (but not differentiable) - # function of bred. - jc3 = bred * (1.0 - bred) - if not np.isfinite(jc3): - raise RuntimeError("jc3 jcrit is NaN.") + # Critical field (T) at given strain and temperature + bcrit = bc20eps * (1.0 - t**1.52) - # Critical current density in superconductor (A/m2) jc1 = (csc / bmax) * strfun - if t > 0.0: - jc2 = (1.0 - t**1.52) * (1.0 - t**2) - else: - # Allow t<0, fudge to give real value of jc2 - # This generates a real and continuous (but not differentiable) function of t. - jc2 = (1.0 - t) * (1.0 - t**2) - - # jc3 = bred**p * (1.0D0-bred)**q ! bred must be < 1 to avoid NaNs + # Check if we are inside the critical surface + if (t > 0) and (t < 1) and (bmax > 0) and (bmax < bcrit) and (bcrit > 0): + # Reduced field at given strain and temperature + bred = bmax / bcrit - # scale from mm2 to m2 - scalefac = 1.0e6 - - jcrit = jc1 * jc2 * jc3 * scalefac - if not np.isfinite(jcrit): - raise RuntimeError("jcrit is NaN.") + jc2 = (1.0 - t**1.52) * (1.0 - t**2) + jc3 = bred**p * (1.0 - bred) ** q + jscaling = jc1 * jc2 * jc3 - return jcrit, bcrit, tcrit + else: + # Outside the critical surface. + # We construct a simple function that is always negative and + # becomes more negative as field and temperature increase. + jc2 = t + jc3 = bmax / max(bcrit, 1.0e-8) + jscaling = -abs(jc1 * jc2 * jc3) + + return jscaling, bcrit, tcrit def croco(jcritsc, conductor_area, croco_od, croco_thick): diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py index 8c462355f6..5242628df8 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -17,8 +17,10 @@ from process.fortran import pfcoil_variables as pfv from process.fortran import physics_variables as pv from process.fortran import error_handling as eh + from process.fortran import fwbs_variables as fwbsv from process.fortran import tfcoil_variables as tfv + from process.fortran import times_variables as tv from process.fortran import constants from process.pfcoil import PFCoil, mtrx, fixb @@ -63,7 +65,7 @@ def test_pfcoil(monkeypatch, pfcoil): monkeypatch.setattr(pfv, "fcohbop", 1.0) monkeypatch.setattr(pfv, "rjconpf", np.full(22, 1.1e7)) monkeypatch.setattr(pfv, "ngrp", 4) - monkeypatch.setattr(pfv, "rohc", 0.0) + monkeypatch.setattr(pfv, "rohc", 3.0) monkeypatch.setattr(pfv, "ncls", np.array([1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0])) monkeypatch.setattr(pfv, "zpf", np.full(22, 0.0)) monkeypatch.setattr(pfv, "cptdin", np.full(22, 4.22e4)) @@ -194,7 +196,7 @@ def test_ohcalc(monkeypatch, reinitialise_error_module, pfcoil): monkeypatch.setattr(pfv, "bmaxoh", 1.4e1) monkeypatch.setattr(pfv, "i_cs_stress", 0) monkeypatch.setattr(pfv, "coheof", 1.693e7) - monkeypatch.setattr(pfv, "rohc", 0.0) + monkeypatch.setattr(pfv, "rohc", 3.0) monkeypatch.setattr(pfv, "vfohc", 3.0e-1) monkeypatch.setattr(pfv, "jstrandoh_bop", 1.069e8) monkeypatch.setattr(pfv, "fcuohsu", 7.000e-1) @@ -265,11 +267,11 @@ def test_ohcalc(monkeypatch, reinitialise_error_module, pfcoil): pfcoil.ohcalc() - assert pytest.approx(pfv.bpf[4]) == 9.299805e2 - assert pytest.approx(pfv.rjohc) == -7.728453e9 + assert pytest.approx(pfv.bpf[4]) == 13.073958753751993 + assert pytest.approx(pfv.rjohc) == 54101481.7685945 -def test_efc(pfcoil, monkeypatch): +def test_efc(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): """Test efc subroutine. efc() requires specific arguments in order to work; these were discovered @@ -459,7 +461,7 @@ def test_efc(pfcoil, monkeypatch): ) -def test_mtrx(pfcoil): +def test_mtrx(pfcoil: PFCoil): """Test mtrx subroutine. mtrx() requires specific arguments in order to work; these were discovered @@ -1640,7 +1642,7 @@ def test_mtrx(pfcoil): assert_array_almost_equal(bvec, bvec_exp) -def test_solv(pfcoil): +def test_solv(pfcoil: PFCoil): """Test solv() with simple arguments. Running baseline_2019 results in 2D array args with 740 elements: unfeasible @@ -1686,7 +1688,7 @@ def test_solv(pfcoil): ) -def test_fixb(pfcoil): +def test_fixb(pfcoil: PFCoil): """Test fixb subroutine. fixb() requires specific arguments in order to work; these were discovered @@ -2026,7 +2028,7 @@ def test_fixb(pfcoil): assert_array_almost_equal(bfix, bfix_exp) -def test_peakb(monkeypatch, pfcoil): +def test_peakb(monkeypatch: pytest.MonkeyPatch, pfcoil: PFCoil): """Test peakb subroutine. peakb() requires specific arguments in order to work; these were discovered @@ -2626,7 +2628,7 @@ def test_peakb(monkeypatch, pfcoil): assert pytest.approx(bzo) == bzo_exp -def test_superconpf(monkeypatch, pfcoil): +def test_superconpf(monkeypatch: pytest.MonkeyPatch, pfcoil: PFCoil): """Test superconpf subroutine. superconpf() requires specific arguments in order to work; these were @@ -2672,7 +2674,7 @@ def test_superconpf(monkeypatch, pfcoil): assert pytest.approx(tmarg) == tmarg_exp -def test_axial_stress(pfcoil, monkeypatch): +def test_axial_stress(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): """Test axial_stress subroutine. axial_stress() requires specific mocked vars in order to work; these were @@ -2815,7 +2817,7 @@ def test_axial_stress(pfcoil, monkeypatch): assert pytest.approx(axial_force) == axial_force_exp -def test_induct(pfcoil, monkeypatch): +def test_induct(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): """Test induct subroutine. induct() requires specific mocked vars in order to work; these were diff --git a/tests/unit/test_sctfcoil.py b/tests/unit/test_sctfcoil.py index 1a1c8c5233..3dc2a64fb5 100644 --- a/tests/unit/test_sctfcoil.py +++ b/tests/unit/test_sctfcoil.py @@ -210,11 +210,11 @@ class SuperconParam(NamedTuple): tmax=150, bcritsc=24, tcritsc=16, - expected_temp_margin=2.3431632224075836, + expected_temp_margin=2.34312129, expected_jwdgpro=17475706.393616617, expected_jwdgcrt=41107234.360397324, expected_vd=9988.2637896807955, - expected_tmarg=2.3431632224075836, + expected_tmarg=2.34312129, ), SuperconParam( tmargmin_tf=1.5, @@ -250,11 +250,11 @@ class SuperconParam(NamedTuple): tmax=150, bcritsc=24, tcritsc=16, - expected_temp_margin=2.3431632224075836, + expected_temp_margin=2.34312129, expected_jwdgpro=17475706.393616617, expected_jwdgcrt=41107234.360397324, expected_vd=10001.287165953383, - expected_tmarg=2.3431632224075836, + expected_tmarg=2.34312129, ), SuperconParam( tmargmin_tf=1.5, @@ -290,11 +290,11 @@ class SuperconParam(NamedTuple): tmax=150, bcritsc=24, tcritsc=16, - expected_temp_margin=2.3431632224075836, + expected_temp_margin=2.34312129, expected_jwdgpro=17475706.393616617, expected_jwdgcrt=41107234.360397324, expected_vd=10001.287165953383, - expected_tmarg=2.3431632224075836, + expected_tmarg=2.34312129, ), ), ) diff --git a/tests/unit/test_superconductors.py b/tests/unit/test_superconductors.py index 78309774a7..5f35508f1e 100644 --- a/tests/unit/test_superconductors.py +++ b/tests/unit/test_superconductors.py @@ -61,7 +61,7 @@ def test_itersc(iterscparam): """ jcrit, bcrit, tcrit = superconductors.itersc( - thelium=iterscparam.thelium, + temperature=iterscparam.thelium, bmax=iterscparam.bmax, strain=iterscparam.strain, bc20max=iterscparam.bc20max, @@ -148,8 +148,8 @@ def test_jcrit_rebco(): def test_current_sharing_rebco(): - assert superconductors.current_sharing_rebco(7.0, 2e7) == pytest.approx( - 75.76286550648135 + assert superconductors.current_sharing_rebco(11.0, 2e7) == pytest.approx( + 71.28702697627514 )