Skip to content

3131 errror in calculation of temperature margin bug #3142

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jonmaddock merged 2 commits into from
May 13, 2024
Merged

3131 errror in calculation of temperature margin bug #3142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
319965b
replaced the Newton-Raphson method in by scipy.optimize.newton
Apr 12, 2024
45e0187
New function Bottura_scaling to simplify supercon routines.
Apr 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 0 additions & 64 deletions process/maths_library.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
276 changes: 33 additions & 243 deletions process/pfcoil.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,15 @@
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.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
import logging
from scipy import optimize

logger = logging.getLogger(__name__)

Expand Down Expand Up @@ -126,7 +127,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,
Expand Down Expand Up @@ -624,6 +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]))

pfv.rjpfalw[i], jstrand, jsc, tmarg = self.superconpf(
bmax,
pfv.vf[i],
Expand Down Expand Up @@ -982,12 +984,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

Expand Down Expand Up @@ -2755,97 +2754,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
Expand Down Expand Up @@ -2970,155 +2878,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

if isumat == 3:
arguments = (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
arguments = (isumat, jsc, bmax, strain, bc20m, tc0m)

another_estimate = 2 * thelium
t_zero_margin, root_result = optimize.newton(
Comment thread
mkovari marked this conversation as resolved.
superconductors.current_density_margin,
thelium,
fprime=None,
args=arguments,
tol=1.0e-06,
maxiter=50,
fprime2=None,
x1=another_estimate,
rtol=1.0e-6,
full_output=True,
disp=False,
)
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
)
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

Expand Down
Loading