From bdedc5b799984d63309ab2f7d1f7ae1155eefc47 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 4 Dec 2024 15:52:52 +0000 Subject: [PATCH 01/13] Use scipy elliptic integrals --- process/pfcoil.py | 7 ++-- source/fortran/maths_library.f90 | 55 -------------------------------- 2 files changed, 5 insertions(+), 57 deletions(-) diff --git a/process/pfcoil.py b/process/pfcoil.py index f8ed6928c3..ae16de1778 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -4,6 +4,7 @@ import numba import numpy as np from scipy import optimize +from scipy.special import ellipe, ellipk import process.superconductors as superconductors from process import fortran as ft @@ -1599,13 +1600,15 @@ def axial_stress(self): axial_term_1 = -(constants.rmu0 / 2.0e0) * (ni / (2.0e0 * hl)) ** 2 # term 2 - ekb2_1, ekb2_2 = ml.ellipke(kb2) + ekb2_1 = ellipk(kb2) + ekb2_2 = ellipe(kb2) axial_term_2 = ( 2.0e0 * hl * (math.sqrt(4.0e0 * b**2 + hl**2)) * (ekb2_1 - ekb2_2) ) # term 3 - ek2b2_1, ek2b2_2 = ml.ellipke(k2b2) + ek2b2_1 = ellipk(k2b2) + ek2b2_2 = ellipe(k2b2) axial_term_3 = ( 2.0e0 * hl * (math.sqrt(4.0e0 * b**2 + 4.0e0 * hl**2)) * (ek2b2_1 - ek2b2_2) ) diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 766e6d8b6a..5980ad2ac7 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -100,61 +100,6 @@ function find_y_nonuniform_x(x0,x,y,n) end function find_y_nonuniform_x - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine ellipke(sqk,kk,ek) - - !! Routine that calculates the complete elliptic integral - !! of the first and second kinds - !! author: P J Knight, CCFE, Culham Science Centre - !! sqk : input real : square of the elliptic modulus - !! kk : output real : complete elliptic integral of the first kind - !! ek : output real : complete elliptic integral of the second kind - !! This routine calculates the complete elliptic integral - !! of the first and second kinds. - !!

The method used is that described in the reference, and - !! the code is taken from the Culham maglib library routine FN02A. - !! Approximations for Digital Computers, C. Hastings, - !! Princeton University Press, 1955 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - ! Arguments - - real(dp), intent(in) :: sqk - real(dp), intent(out) :: kk,ek - - ! Local variables - - real(dp) :: a,b,d,e - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - d = 1.0D0 - sqk - e = log(d) - - ! Evaluate series for integral of first kind - - a = (((0.014511962D0*d + 0.037425637D0)*d + 0.035900924D0)*d & - + 0.096663443D0)*d + 1.386294361D0 - b = (((0.004417870D0*d + 0.033283553D0)*d + 0.06880249D0)*d & - + 0.12498594D0)*d + 0.5D0 - - kk = a - b*e - - ! Evaluate series for integral of second kind - - a = (((0.017365065D0*d + 0.047573835D0)*d + 0.06260601D0)*d & - + 0.44325141D0)*d + 1.0D0 - b = (((0.005264496D0*d + 0.040696975D0)*d + 0.09200180D0)*d & - + 0.24998368D0)*d - - ek = a - b*e - - end subroutine ellipke - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function binomial(n,k) result(coefficient) ! This outputs a real approximation to the coefficient From 1ec35d8720eceb300d2715e036b53f8eeb906a90 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 4 Dec 2024 15:59:46 +0000 Subject: [PATCH 02/13] Use scipy comb over binomial --- process/availability.py | 6 +-- source/fortran/maths_library.f90 | 19 ---------- tests/unit/test_maths_library.py | 64 -------------------------------- 3 files changed, 3 insertions(+), 86 deletions(-) delete mode 100644 tests/unit/test_maths_library.py diff --git a/process/availability.py b/process/availability.py index 3690d24afd..5581078438 100644 --- a/process/availability.py +++ b/process/availability.py @@ -1,13 +1,14 @@ import logging import math +from scipy.special import comb as combinations + from process import fortran as ft from process.fortran import constraint_variables as ctv from process.fortran import cost_variables as cv from process.fortran import divertor_variables as dv from process.fortran import fwbs_variables as fwbsv from process.fortran import ife_variables as ifev -from process.fortran import maths_library from process.fortran import physics_variables as pv from process.fortran import process_output as po from process.fortran import tfcoil_variables as tfv @@ -953,10 +954,9 @@ def calc_u_unplanned_vacuum(self, output: bool) -> float: for n in range(cv.redun_vac + 1, total_pumps + 1): # Probability for n failures in the operational period, n > number of redundant pumps - # vac_fail_p.append(maths_library.binomial(total_pumps,n) * (cryo_nfailure_rate**(total_pumps-n)) *(cryo_failure_rate**n)) # calculate sum in formula for downtime - sum_prob = sum_prob + maths_library.binomial(total_pumps, n) * ( + sum_prob = sum_prob + combinations(total_pumps, n) * ( cryo_nfailure_rate ** (total_pumps - n) ) * (cryo_failure_rate**n) * (n - cv.redun_vac) diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 5980ad2ac7..273f1054bd 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -100,25 +100,6 @@ function find_y_nonuniform_x(x0,x,y,n) end function find_y_nonuniform_x - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - function binomial(n,k) result(coefficient) - ! This outputs a real approximation to the coefficient - ! http://en.wikipedia.org/wiki/Binomial_coefficient#Multiplicative_formula - implicit none - integer, intent(in) :: n, k - real(dp) :: coefficient - integer :: numerator, i - if (k == 0) then - coefficient = 1 - else - coefficient = 1.0D0 - do i = 1, k - numerator = n + 1 -i - coefficient = coefficient * real(numerator)/real(i) - end do - end if - end function binomial - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1) diff --git a/tests/unit/test_maths_library.py b/tests/unit/test_maths_library.py deleted file mode 100644 index ceec508589..0000000000 --- a/tests/unit/test_maths_library.py +++ /dev/null @@ -1,64 +0,0 @@ -"""Unit tests for maths_library.f90.""" - -import pytest - -from process import fortran - - -def binomial_param(**kwargs): - """Make parameters for a single binomial() test. - - :return: Parameters, including expected result - :rtype: dict - """ - # Default parameters - defaults = {"n": 1, "k": 1, "expected": 1.0} - - # Merge default dict with any optional keyword arguments to override values - return {**defaults, **kwargs} - - -def binomial_params(): - """Create a list of parameter dicts for the calc_u_planned fixture. - - Case 1: 1, 1 - Case 2: 3, 1 - Case 3: 3, 3 - - :return: List of parameter dicts - :rtype: list - """ - return [ - binomial_param(), - binomial_param(n=3, expected=3.0), - binomial_param(n=3, k=3, expected=1.0), - ] - - -@pytest.fixture(params=binomial_params(), ids=["11", "31", "33"]) -def binomial_fix(request): - """Fixture for the binomial() variables. - - :param request: Request object for accessing parameters - :type request: object - :return: Parameterised arguments and expected return value for binomial() - :rtype: dict - """ - return request.param - - -def test_binomial(binomial_fix): - """Test binomial(). - - :param binomial_fix: Parameterised arguments and expected return value for - binomial() - :type binomial_fix: dict - """ - # Run binomial() with the current fixture, - # then assert the result is the expected one - n = binomial_fix["n"] - k = binomial_fix["k"] - expected = binomial_fix["expected"] - - result = fortran.maths_library.binomial(n, k) - assert result == expected From 956723862d3562adb7b6430ae2a612463369b7c2 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 4 Dec 2024 16:19:22 +0000 Subject: [PATCH 03/13] Interpolate instead of find_y_nonuniform_x --- process/stellarator.py | 17 ++++--- source/fortran/maths_library.f90 | 77 -------------------------------- 2 files changed, 8 insertions(+), 86 deletions(-) diff --git a/process/stellarator.py b/process/stellarator.py index d549b477ef..d8af567401 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -19,7 +19,6 @@ global_variables, heat_transport_variables, impurity_radiation_module, - maths_library, neoclassics_module, numerics, physics_module, @@ -3123,8 +3122,8 @@ def j_max_protect_am2(self, tau_quench, t_detect, fcu, fcond, temp, acs, aturn): 1.32773e14, ] - q_he = maths_library.find_y_nonuniform_x(temp, temp_k, q_he_array_sa2m4, 13) - q_cu = maths_library.find_y_nonuniform_x(temp, temp_k, q_cu_array_sa2m4, 13) + q_he = np.interp(temp, temp_k, q_he_array_sa2m4) + q_cu = np.interp(temp, temp_k, q_cu_array_sa2m4) # This leaves out the contribution from the superconductor fraction for now return (acs / aturn) * np.sqrt( @@ -3348,8 +3347,8 @@ def intersect(self, x1, y1, x2, y2, xin): for _i in range(100): # Find difference in y values at x - y01 = maths_library.find_y_nonuniform_x(x, x1, y1, n1) - y02 = maths_library.find_y_nonuniform_x(x, x2, y2, n2) + y01 = np.interp(x, x1, y1) + y02 = np.interp(x, x2, y2) y = y01 - y02 if abs(y) < epsy: @@ -3357,14 +3356,14 @@ def intersect(self, x1, y1, x2, y2, xin): # Find difference in y values at x+dx - y01 = maths_library.find_y_nonuniform_x(x + dx, x1, y1, n1) - y02 = maths_library.find_y_nonuniform_x(x + dx, x2, y2, n2) + y01 = np.interp(x + dx, x1, y1) + y02 = np.interp(x + dx, x2, y2) yright = y01 - y02 # Find difference in y values at x-dx - y01 = maths_library.find_y_nonuniform_x(x - dx, x1, y1, n1) - y02 = maths_library.find_y_nonuniform_x(x - dx, x2, y2, n2) + y01 = np.interp(x - dx, x1, y1) + y02 = np.interp(x - dx, x2, y2) yleft = y01 - y02 # Adjust x using Newton-Raphson method diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 273f1054bd..ba8b9b5957 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -23,83 +23,6 @@ module maths_library implicit none contains - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function find_y_nonuniform_x(x0,x,y,n) - - !! Routine to find y0 such that y0 = y(x0) given a set of - !! values x(1:n), y(1:n) - !! author: P J Knight, CCFE, Culham Science Centre - !! x0 : input real : x value at which we want to find y - !! x(1:n) : input real array : monotonically de/increasing x values - !! y(1:n) : input real array : y values at each x - !! n : input integer : size of array - !! This routine performs a simple linear interpolation method - !! to find the y value at x = x0. If x0 lies outside the - !! range [x(1),x(n)], the y value at the nearest 'end' of the data - !! is returned. - !! None - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - real(dp) :: find_y_nonuniform_x - - ! Arguments - - integer, intent(in) :: n - real(dp), intent(in) :: x0 - real(dp), dimension(n), intent(in) :: x - real(dp), dimension(n), intent(in) :: y - - ! Local variables - integer :: i,j - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Step through arrays until x crosses the value of interest - - j = 0 - rough_search: do i = 1,n-1 - - if (((x(i)-x0)*(x(i+1)-x0)) <= 0.0D0) then - j = i - exit rough_search - end if - - end do rough_search - - if (j /= 0) then - - ! Simply do a linear interpolation between the two grid points - ! spanning the point of interest - - find_y_nonuniform_x = y(j) + (y(j+1)-y(j))*(x0-x(j))/(x(j+1)-x(j)) - - else ! No points found, so return the 'nearest' y value - - if (x(n) > x(1)) then ! values are monotonically increasing - if (x0 > x(n)) then - find_y_nonuniform_x = y(n) - else - find_y_nonuniform_x = y(1) - end if - - else ! values are monotonically decreasing - if (x0 < x(n)) then - find_y_nonuniform_x = y(n) - else - find_y_nonuniform_x = y(1) - end if - - end if - - end if - - end function find_y_nonuniform_x - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1) From 0f2a66fa240a582ece1a0a75af9b2ff46d66f7b8 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 4 Dec 2024 17:35:49 +0000 Subject: [PATCH 04/13] Use svd from numpy --- process/pfcoil.py | 18 +- source/fortran/maths_library.f90 | 411 --------------------------- tests/integration/test_pfcoil_int.py | 69 ++--- 3 files changed, 31 insertions(+), 467 deletions(-) diff --git a/process/pfcoil.py b/process/pfcoil.py index ae16de1778..cab2f3ecdf 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -4,6 +4,7 @@ import numba import numpy as np from scipy import optimize +from scipy.linalg import svd from scipy.special import ellipe, ellipk import process.superconductors as superconductors @@ -14,7 +15,6 @@ from process.fortran import cs_fatigue_variables as csfv from process.fortran import error_handling as eh from process.fortran import fwbs_variables as fwbsv -from process.fortran import maths_library as ml from process.fortran import pfcoil_module as pf from process.fortran import pfcoil_variables as pfv from process.fortran import physics_variables as pv @@ -896,7 +896,7 @@ def efc( ) # Solve matrix equation - ccls, umat, vmat, sigma, work2 = self.solv(pfv.ngrpmx, ngrp, nrws, gmat, bvec) + ccls = self.solv(pfv.ngrpmx, ngrp, nrws, gmat, bvec) # Calculate the norm of the residual vectors brssq, brnrm, bzssq, bznrm, ssq = rsid( @@ -975,13 +975,10 @@ def solv(self, ngrpmx, ngrp, nrws, gmat, bvec): :rtype: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray] """ - truth = True - eps = 1.0e-10 ccls = np.zeros(ngrpmx) + work2 = np.zeros(ngrpmx) - sigma, umat, vmat, ierr, work2 = ml.svd( - nrws, np.asfortranarray(gmat), truth, truth - ) + umat, sigma, vmat = svd(gmat) for i in range(ngrp): work2[i] = 0.0e0 @@ -990,15 +987,14 @@ def solv(self, ngrpmx, ngrp, nrws, gmat, bvec): # Compute currents for i in range(ngrp): - ccls[i] = 0.0e0 zvec = 0.0e0 for j in range(ngrp): - if sigma[j] > eps: + if sigma[j] > 1.0e-10: zvec = work2[j] / sigma[j] - ccls[i] = ccls[i] + vmat[i, j] * zvec + ccls[i] = ccls[i] + vmat[j, i] * zvec - return ccls, umat, vmat, sigma, work2 + return ccls def ohcalc(self): """Routine to perform calculations for the Central Solenoid. diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index ba8b9b5957..79e0c7ac3d 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -25,417 +25,6 @@ module maths_library contains ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1) - - !! Singular Value Decomposition - !! author: P J Knight, CCFE, Culham Science Centre - !! author: B. S. Garbow, Applied Mathematics Division, Argonne National Laboratory - !! nm : input integer : Max number of rows of arrays a, u, v; >= m,n - !! m : input integer : Actual number of rows of arrays a, u - !! n : input integer : Number of columns of arrays a, u, and the order of v - !! a(nm,n) : input/output real array : On input matrix to be decomposed; - !! on output, either unchanged or overwritten with u or v - !! w(n) : output real array : The n (non-negative) singular values of a - !! (the diagonal elements of s); unordered. If an error exit - !! is made, the singular values should be correct for indices - !! ierr+1,ierr+2,...,n. - !! matu : input logical : Set to .true. if the u matrix in the - !! decomposition is desired, and to .false. otherwise. - !! u(nm,n) : output real array : The matrix u (orthogonal column vectors) - !! of the decomposition if matu has been set to .true., otherwise - !! u is used as a temporary array. u may coincide with a. - !! If an error exit is made, the columns of u corresponding - !! to indices of correct singular values should be correct. - !! matv : input logical : Set to .true. if the v matrix in the - !! decomposition is desired, and to .false. otherwise. - !! v(nm,n) : output real array : The matrix v (orthogonal) of the - !! decomposition if matv has been set to .true., otherwise - !! v is not referenced. v may also coincide with a if u is - !! not needed. If an error exit is made, the columns of v - !! corresponding to indices of correct singular values - !! should be correct. - !! ierr : output integer : zero for normal return, or k if the - !! k-th singular value has not been determined after 30 iterations. - !! rv1(n) : output real array : work array - !! This subroutine is a translation of the algol procedure SVD, - !! Num. Math. 14, 403-420(1970) by Golub and Reinsch, - !! Handbook for Auto. Comp., vol II - Linear Algebra, 134-151(1971). - !!

It determines the singular value decomposition - !! a=usvt of a real m by n rectangular matrix. - !! Householder bidiagonalization and a variant of the QR - !! algorithm are used. - !! None - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - ! Arguments - - integer, intent(in) :: nm, m, n - logical, intent(in) :: matu, matv - real(dp), dimension(nm,n), intent(inout) :: a - real(dp), dimension(nm,n), intent(out) :: u, v - real(dp), dimension(n), intent(out) :: w, rv1 - integer, intent(out) :: ierr - - ! Local variables - - integer :: i,j,k,l,ii,i1,kk,k1,ll,l1,mn,its - real(dp) :: c,f,g,h,s,x,y,z,scale,anorm - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ierr = 0 - - u = a - - ! Householder reduction to bidiagonal form - - g = 0.0D0 - scale = 0.0D0 - anorm = 0.0D0 - - do i = 1, n - - l = i + 1 - rv1(i) = scale * g - g = 0.0D0 - s = 0.0D0 - scale = 0.0D0 - - if (i <= m) then - - do k = i, m - scale = scale + abs(u(k,i)) - end do - - if (scale /= 0.0D0) then - - do k = i, m - u(k,i) = u(k,i) / scale - s = s + u(k,i)**2 - end do - - f = u(i,i) - g = -sign(sqrt(s),f) - h = f * g - s - u(i,i) = f - g - - if (i /= n) then - do j = l, n - s = 0.0D0 - do k = i, m - s = s + u(k,i) * u(k,j) - end do - f = s / h - do k = i, m - u(k,j) = u(k,j) + f * u(k,i) - end do - end do - end if - - do k = i, m - u(k,i) = scale * u(k,i) - end do - - end if - - end if - - w(i) = scale * g - g = 0.0D0 - s = 0.0D0 - scale = 0.0D0 - - if (.not.((i > m) .or. (i == n))) then - - do k = l, n - scale = scale + abs(u(i,k)) - end do - - if (scale /= 0.0D0) then - - do k = l, n - u(i,k) = u(i,k) / scale - s = s + u(i,k)**2 - end do - - f = u(i,l) - g = -sign(sqrt(s),f) - h = f * g - s - u(i,l) = f - g - - do k = l, n - rv1(k) = u(i,k) / h - end do - - if (i /= m) then - do j = l, m - s = 0.0D0 - do k = l, n - s = s + u(j,k) * u(i,k) - end do - do k = l, n - u(j,k) = u(j,k) + s * rv1(k) - end do - end do - end if - - do k = l, n - u(i,k) = scale * u(i,k) - end do - - end if - - end if - - anorm = max(anorm,abs(w(i))+abs(rv1(i))) - - end do ! i - - ! Accumulation of right-hand transformations - - if (matv) then - - ! For i=n step -1 until 1 do - do ii = 1, n - i = n + 1 - ii - if (i /= n) then - - if (g /= 0.0D0) then - do j = l, n - ! Double division avoids possible underflow - v(j,i) = (u(i,j) / u(i,l)) / g - end do - do j = l, n - s = 0.0D0 - do k = l, n - s = s + u(i,k) * v(k,j) - end do - do k = l, n - v(k,j) = v(k,j) + s * v(k,i) - end do - end do - end if - - do j = l, n - v(i,j) = 0.0D0 - v(j,i) = 0.0D0 - end do - - end if - - v(i,i) = 1.0D0 - g = rv1(i) - l = i - end do - - end if - - ! Accumulation of left-hand transformations - - if (matu) then - - ! For i=min(m,n) step -1 until 1 do - mn = n - if (m < n) mn = m - - do ii = 1, mn - i = mn + 1 - ii - l = i + 1 - g = w(i) - if (i /= n) then - do j = l, n - u(i,j) = 0.0D0 - end do - end if - - if (g /= 0.0D0) then - - if (i /= mn) then - do j = l, n - s = 0.0D0 - do k = l, m - s = s + u(k,i) * u(k,j) - end do - f = (s / u(i,i)) / g ! Double division avoids possible underflow - do k = i, m - u(k,j) = u(k,j) + f * u(k,i) - end do - end do - end if - - do j = i, m - u(j,i) = u(j,i) / g - end do - - else - do j = i, m - u(j,i) = 0.0D0 - end do - end if - - u(i,i) = u(i,i) + 1.0D0 - - end do - - end if - - ! Diagonalization of the bidiagonal form - ! For k=n step -1 until 1 do - - do kk = 1, n - k1 = n - kk - k = k1 + 1 - its = 0 - - ! Test for splitting. - ! For l=k step -1 until 1 do - - do - do ll = 1, k - l1 = k - ll - l = l1 + 1 - if ((abs(rv1(l)) + anorm) == anorm) goto 470 - - ! rv1(1) is always zero, so there is no exit - ! through the bottom of the loop - - !+**PJK 23/05/06 Prevent problems from the code getting here with l1=0 - if (l1 == 0) then - write(*,*) 'SVD: Shouldn''t get here...' - goto 470 - end if - - if ((abs(w(l1)) + anorm) == anorm) exit - end do - - ! Cancellation of rv1(l) if l greater than 1 - - c = 0.0D0 - s = 1.0D0 - - do i = l, k - f = s * rv1(i) - rv1(i) = c * rv1(i) - if ((abs(f) + anorm) == anorm) exit - g = w(i) - h = sqrt(f*f+g*g) - w(i) = h - c = g / h - s = -f / h - if (.not. matu) cycle - - do j = 1, m - y = u(j,l1) - z = u(j,i) - u(j,l1) = y * c + z * s - u(j,i) = -y * s + z * c - end do - end do - -470 continue - - ! Test for convergence - - z = w(k) - if (l == k) exit - - ! Shift from bottom 2 by 2 minor - - if (its == 30) then - ! Set error - no convergence to a - ! singular value after 30 iterations - ierr = k - return - end if - - its = its + 1 - x = w(l) - y = w(k1) - g = rv1(k1) - h = rv1(k) - f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.D0 * h * y) - g = sqrt(f*f+1.D0) - f = ((x - z) * (x + z) + h * (y / (f + sign(g,f)) - h)) / x - - ! Next QR transformation - - c = 1.0D0 - s = 1.0D0 - - do i1 = l, k1 - i = i1 + 1 - g = rv1(i) - y = w(i) - h = s * g - g = c * g - z = sqrt(f*f+h*h) - rv1(i1) = z - c = f / z - s = h / z - f = x * c + g * s - g = -x * s + g * c - h = y * s - y = y * c - - if (matv) then - do j = 1, n - x = v(j,i1) - z = v(j,i) - v(j,i1) = x * c + z * s - v(j,i) = -x * s + z * c - end do - end if - - z = sqrt(f*f+h*h) - w(i1) = z - - ! Rotation can be arbitrary if z is zero - - if (z /= 0.0D0) then - c = f / z - s = h / z - end if - - f = c * g + s * y - x = -s * g + c * y - if (.not. matu) cycle - - do j = 1, m - y = u(j,i1) - z = u(j,i) - u(j,i1) = y * c + z * s - u(j,i) = -y * s + z * c - end do - - end do - - rv1(l) = 0.0D0 - rv1(k) = f - w(k) = x - - end do - - ! Convergence - - if (z >= 0.0D0) cycle - - ! w(k) is made non-negative - w(k) = -z - if (.not. matv) cycle - - do j = 1, n - v(j,k) = -v(j,k) - end do - - end do - - end subroutine svd - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine eshellvol(rshell,rmini,rmino,zminor,drin,drout,dz,vin,vout,vtot) !! Routine to calculate the inboard, outboard and total volumes diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py index 7b21838949..aeaf748eb4 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -126,7 +126,7 @@ def test_pfcoil(monkeypatch, pfcoil): monkeypatch.setattr(pv, "vsind", 3.497e2) monkeypatch.setattr(pv, "aspect", 3.1) monkeypatch.setattr(pv, "itart", 0) - monkeypatch.setattr(pv, "beta_poloidal", 6.313e-1) + monkeypatch.setattr(pv, "betap", 6.313e-1) monkeypatch.setattr(tfv, "tftmp", 4.750) monkeypatch.setattr(tfv, "dcond", np.full(9, 9.0e3)) monkeypatch.setattr(tfv, "i_tf_sup", 1) @@ -446,12 +446,13 @@ def test_efc(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): ) assert pytest.approx(ssq) == 4.208729e-4 - assert pytest.approx(ccls[0:4]) == np.array([ - 12846165.42893886, - 16377261.02000236, - 579111.6216917, - 20660782.82356247, - ]) + assert ccls[0:3] == pytest.approx( + np.array([ + 12846165.42893886, + 16377261.02000236, + 579111.6216917, + ]) + ) def test_mtrx(pfcoil: PFCoil): @@ -1642,31 +1643,9 @@ def test_solv(pfcoil: PFCoil): gmat = np.full((3, 3), 2.0, order="F") bvec = np.full(3, 1.0) - ccls, umat, vmat, sigma, work2 = pfcoil.solv(ngrpmx, ngrp, nrws, gmat, bvec) + ccls = pfcoil.solv(ngrpmx, ngrp, nrws, gmat, bvec) - assert_array_almost_equal(ccls, np.array([0.16666667, 0.37079081, -0.03745748])) - assert_array_almost_equal( - umat, - np.array([ - [-0.81649658, -0.57735027, 0.0], - [0.40824829, -0.57735027, -0.70710678], - [0.40824829, -0.57735027, 0.70710678], - ]), - ) - assert_array_almost_equal( - vmat, - np.array([ - [-0.81649658, -0.57735027, 0.0], - [0.40824829, -0.57735027, -0.70710678], - [0.40824829, -0.57735027, 0.70710678], - ]), - ) - assert_array_almost_equal( - sigma, np.array([5.1279005e-16, 6.0000000e00, 0.0000000e00]) - ) - assert_array_almost_equal( - work2, np.array([-2.22044605e-16, -1.73205081e00, 0.00000000e00]) - ) + assert_array_almost_equal(ccls, np.array([-0.069036, 0.488642, 0.080394])) def test_fixb(pfcoil: PFCoil): @@ -2451,13 +2430,13 @@ def test_peakb(monkeypatch: pytest.MonkeyPatch, pfcoil: PFCoil): pfv, "curpfb", np.array([ - 14.742063826112622, - 20.032681634901664, - 0.58040662653667285, - 0.58040662653667285, - 0.42974674788703021, - 0.42974674788703021, - 174.22748790786324, + 0.067422231232391661, + -2.9167273287450968, + -8.1098913365453491, + -8.1098913365453491, + -5.5984385047179153, + -5.5984385047179153, + -186.98751599968148, 0, 0, 0, @@ -2507,13 +2486,13 @@ def test_peakb(monkeypatch: pytest.MonkeyPatch, pfcoil: PFCoil): pfv, "curpfs", np.array([ - 0.067422231232391661, - -2.9167273287450968, - -8.1098913365453491, - -8.1098913365453491, - -5.5984385047179153, - -5.5984385047179153, - -186.98751599968148, + 14.742063826112622, + 20.032681634901664, + 0.58040662653667285, + 0.58040662653667285, + 0.42974674788703021, + 0.42974674788703021, + 174.22748790786324, 0, 0, 0, From 76a0c28e686a56ea0eedd3b69ee0b683db281b3c Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 11 Dec 2024 11:08:44 +0000 Subject: [PATCH 05/13] Convert eshellarea to Python --- process/blanket_library.py | 32 ++++++++++++++++++++-- process/build.py | 10 ++----- source/fortran/maths_library.f90 | 46 -------------------------------- 3 files changed, 32 insertions(+), 56 deletions(-) diff --git a/process/blanket_library.py b/process/blanket_library.py index 953b854c40..9d7c33c472 100644 --- a/process/blanket_library.py +++ b/process/blanket_library.py @@ -270,13 +270,13 @@ def elliptical_component(self, icomponent: int): build_variables.blareaib, build_variables.blareaob, build_variables.blarea, - ) = maths_library.eshellarea(r1, r2, r3, blanket_library.hblnkt) + ) = eshellarea(r1, r2, r3, blanket_library.hblnkt) if icomponent == 1: ( build_variables.shareaib, build_variables.shareaob, build_variables.sharea, - ) = maths_library.eshellarea(r1, r2, r3, blanket_library.hshld) + ) = eshellarea(r1, r2, r3, blanket_library.hshld) # Calculate volumes, assuming 100% coverage if icomponent == 0: @@ -2710,3 +2710,31 @@ def pumppower( po.ovarre(self.outfile, "Mass flow rate in (kg/s) = ", "(mf)", mf, "OP ") return pumppower + + +def eshellarea(rshell, rmini, rmino, zminor): + """Routine to calculate the inboard, outboard and total surface areas + of a toroidal shell comprising two elliptical sections + author: P J Knight, CCFE, Culham Science Centre + rshell : input real : major radius of centre of both ellipses (m) + rmini : input real : horizontal distance from rshell to + inboard elliptical shell (m) + rmino : input real : horizontal distance from rshell to + outboard elliptical shell (m) + zminor : input real : vertical internal half-height of shell (m) + ain : output real : surface area of inboard section (m3) + aout : output real : surface area of outboard section (m3) + atot : output real : total surface area of shell (m3) + This routine calculates the surface area of the inboard and outboard + sections of a toroidal shell defined by two co-centred semi-ellipses. + """ + + # Inboard section + elong = zminor / rmini + ain = 2.0 * np.pi * elong * (np.pi * rshell * rmini - 2.0 * rmini * rmini) + + # Outboard section + elong = zminor / rmino + aout = 2.0 * np.pi * elong * (np.pi * rshell * rmino + 2.0 * rmino * rmino) + + return ain, aout, ain + aout diff --git a/process/build.py b/process/build.py index d33b42d8f8..e1351ab7bd 100644 --- a/process/build.py +++ b/process/build.py @@ -2,6 +2,7 @@ import numpy as np +from process.blanket_library import eshellarea from process.fortran import ( blanket_library, build_variables, @@ -2001,9 +2002,6 @@ def calculate_radial_build(self, output: bool) -> None: + build_variables.dr_fw_plasma_gap_outboard ) - r1 # Calculate surface area, assuming 100% coverage - # maths_library.eshellarea was not working across - # the interface so has been reimplemented here - # as a test ( build_variables.fwareaib, @@ -2038,15 +2036,11 @@ def calculate_radial_build(self, output: bool) -> None: # Calculate surface area, assuming 100% coverage - # maths_library.eshellarea was not working across - # the interface so has been reimplemented here - # as a test - ( build_variables.fwareaib, build_variables.fwareaob, build_variables.fwarea, - ) = maths_library.eshellarea(r1, r2, r3, hfw) + ) = eshellarea(r1, r2, r3, hfw) # Apply area coverage factor diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 79e0c7ac3d..30d6434475 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -224,52 +224,6 @@ end subroutine dshellarea ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine eshellarea(rshell,rmini,rmino,zminor,ain,aout,atot) - - !! Routine to calculate the inboard, outboard and total surface areas - !! of a toroidal shell comprising two elliptical sections - !! author: P J Knight, CCFE, Culham Science Centre - !! rshell : input real : major radius of centre of both ellipses (m) - !! rmini : input real : horizontal distance from rshell to - !! inboard elliptical shell (m) - !! rmino : input real : horizontal distance from rshell to - !! outboard elliptical shell (m) - !! zminor : input real : vertical internal half-height of shell (m) - !! ain : output real : surface area of inboard section (m3) - !! aout : output real : surface area of outboard section (m3) - !! atot : output real : total surface area of shell (m3) - !! This routine calculates the surface area of the inboard and outboard - !! sections of a toroidal shell defined by two co-centred semi-ellipses. - !!

See also eshellvol - !! Internal CCFE note T&M/PKNIGHT/PROCESS/009, P J Knight: - !! Surface Area and Volume Calculations for Toroidal Shells - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use constants, only: pi, twopi - implicit none - - ! Arguments - real(dp), intent(in) :: rshell,rmini,rmino,zminor - real(dp), intent(out) :: ain,aout,atot - - ! Local variables - real(dp) :: elong - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Inboard section - elong = zminor/rmini - ain = twopi * elong * (pi*rshell*rmini - 2.0D0*rmini*rmini) - - ! Outboard section - elong = zminor/rmino - aout = twopi * elong * (pi*rshell*rmino + 2.0D0*rmino*rmino) - - ! Total surface area - atot = ain + aout - - end subroutine eshellarea - ! ------------------------------------------------------------------------ pure function variable_error(variable) real(dp), intent(in) ::variable From 29eeda5946acc9fe802a4b496af1782e6ace0348 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 11 Dec 2024 11:11:49 +0000 Subject: [PATCH 06/13] Convert dshellarea to Python --- process/blanket_library.py | 31 +++++++++++++++++++-- process/build.py | 5 ++-- source/fortran/maths_library.f90 | 48 -------------------------------- 3 files changed, 31 insertions(+), 53 deletions(-) diff --git a/process/blanket_library.py b/process/blanket_library.py index 9d7c33c472..0d9afe8df0 100644 --- a/process/blanket_library.py +++ b/process/blanket_library.py @@ -180,13 +180,13 @@ def dshaped_component(self, icomponent: int): build_variables.blareaib, build_variables.blareaob, build_variables.blarea, - ) = maths_library.dshellarea(r1, r2, blanket_library.hblnkt) + ) = dshellarea(r1, r2, blanket_library.hblnkt) if icomponent == 1: ( build_variables.shareaib, build_variables.shareaob, build_variables.sharea, - ) = maths_library.dshellarea(r1, r2, blanket_library.hshld) + ) = dshellarea(r1, r2, blanket_library.hshld) # Calculate volumes, assuming 100% coverage if icomponent == 0: @@ -2738,3 +2738,30 @@ def eshellarea(rshell, rmini, rmino, zminor): aout = 2.0 * np.pi * elong * (np.pi * rshell * rmino + 2.0 * rmino * rmino) return ain, aout, ain + aout + + +def dshellarea(rmajor, rminor, zminor): + """Routine to calculate the inboard, outboard and total surface areas + of a D-shaped toroidal shell + author: P J Knight, CCFE, Culham Science Centre + rmajor : input real : major radius of inboard straight section (m) + rminor : input real : horizontal width of shell (m) + zminor : input real : vertical half-height of shell (m) + ain : output real : surface area of inboard straight section (m3) + aout : output real : surface area of outboard curved section (m3) + atot : output real : total surface area of shell (m3) + This routine calculates the surface area of the inboard and outboard + sections of a D-shaped toroidal shell defined by the above input + parameters. + The inboard section is assumed to be a cylinder. + The outboard section is defined by a semi-ellipse, centred on the + major radius of the inboard section. + """ + # Area of inboard cylindrical shell + ain = 4.0 * zminor * np.pi * rmajor + + # Area of elliptical outboard section + elong = zminor / rminor + aout = 2.0 * np.pi * elong * (np.pi * rmajor * rminor + 2.0 * rminor * rminor) + + return ain, aout, ain + aout diff --git a/process/build.py b/process/build.py index e1351ab7bd..8578d0fe80 100644 --- a/process/build.py +++ b/process/build.py @@ -2,7 +2,7 @@ import numpy as np -from process.blanket_library import eshellarea +from process.blanket_library import dshellarea, eshellarea from process.fortran import ( blanket_library, build_variables, @@ -12,7 +12,6 @@ divertor_variables, error_handling, fwbs_variables, - maths_library, numerics, pfcoil_variables, physics_variables, @@ -2007,7 +2006,7 @@ def calculate_radial_build(self, output: bool) -> None: build_variables.fwareaib, build_variables.fwareaob, build_variables.fwarea, - ) = maths_library.dshellarea(r1, r2, hfw) + ) = dshellarea(r1, r2, hfw) else: # Cross-section is assumed to be defined by two ellipses # Major radius to centre of inboard and outboard ellipses diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 30d6434475..35d36f61a9 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -174,54 +174,6 @@ subroutine dshellvol(rmajor,rminor,zminor,drin,drout,dz,vin,vout,vtot) end subroutine dshellvol -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine dshellarea(rmajor,rminor,zminor,ain,aout,atot) - - !! Routine to calculate the inboard, outboard and total surface areas - !! of a D-shaped toroidal shell - !! author: P J Knight, CCFE, Culham Science Centre - !! rmajor : input real : major radius of inboard straight section (m) - !! rminor : input real : horizontal width of shell (m) - !! zminor : input real : vertical half-height of shell (m) - !! ain : output real : surface area of inboard straight section (m3) - !! aout : output real : surface area of outboard curved section (m3) - !! atot : output real : total surface area of shell (m3) - !! This routine calculates the surface area of the inboard and outboard - !! sections of a D-shaped toroidal shell defined by the above input - !! parameters. - !! The inboard section is assumed to be a cylinder. - !! The outboard section is defined by a semi-ellipse, centred on the - !! major radius of the inboard section. - !!

See also dshellvol - !! Internal CCFE note T&M/PKNIGHT/PROCESS/009, P J Knight: - !! Surface Area and Volume Calculations for Toroidal Shells - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use constants, only: pi, twopi - implicit none - - ! Arguments - real(dp), intent(in) :: rmajor,rminor,zminor - real(dp), intent(out) :: ain,aout,atot - - ! Local variables - real(dp) :: elong - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Area of inboard cylindrical shell - ain = 4.0D0*zminor*pi*rmajor - - ! Area of elliptical outboard section - elong = zminor/rminor - aout = twopi * elong * (pi*rmajor*rminor + 2.0D0*rminor*rminor) - - ! Total surface area - atot = ain + aout - - end subroutine dshellarea - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ------------------------------------------------------------------------ From 7db19677530c2b24fb0f4499f8b3d1bfc4e0fad1 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 11 Dec 2024 11:17:40 +0000 Subject: [PATCH 07/13] Move integer to string routines to main module --- source/fortran/maths_library.f90 | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 35d36f61a9..9b6dfc9c46 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -189,19 +189,4 @@ pure function variable_error(variable) end function variable_error - ! ------------------------------------------------------------------------ - pure function integer2string(value) - ! Convert an integer value to a 2-digit string with leading zero if required. - integer, intent(in) :: value - character(len=2) integer2string - write (integer2string,'(I2.2)') value - end function integer2string - - pure function integer3string(value) - ! Convert an integer value to a 3-digit string with leading zero if required. - integer, intent(in) :: value - character(len=3) integer3string - write (integer3string,'(I3.3)') value - end function integer3string - end module maths_library From 7f6b1d3736dbf4e167c3002796357863f20b48ec Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 11 Dec 2024 11:23:09 +0000 Subject: [PATCH 08/13] Convert eshellvol to Python --- process/blanket_library.py | 70 +++++++++++++++++++++++++-- source/fortran/maths_library.f90 | 82 -------------------------------- 2 files changed, 67 insertions(+), 85 deletions(-) diff --git a/process/blanket_library.py b/process/blanket_library.py index 0d9afe8df0..96093bbf03 100644 --- a/process/blanket_library.py +++ b/process/blanket_library.py @@ -284,7 +284,7 @@ def elliptical_component(self, icomponent: int): fwbs_variables.volblkti, fwbs_variables.volblkto, fwbs_variables.volblkt, - ) = maths_library.eshellvol( + ) = eshellvol( r1, r2, r3, @@ -298,7 +298,7 @@ def elliptical_component(self, icomponent: int): blanket_library.volshldi, blanket_library.volshldo, fwbs_variables.volshld, - ) = maths_library.eshellvol( + ) = eshellvol( r1, r2, r3, @@ -312,7 +312,7 @@ def elliptical_component(self, icomponent: int): blanket_library.volvvi, blanket_library.volvvo, fwbs_variables.vdewin, - ) = maths_library.eshellvol( + ) = eshellvol( r1, r2, r3, @@ -2765,3 +2765,67 @@ def dshellarea(rmajor, rminor, zminor): aout = 2.0 * np.pi * elong * (np.pi * rmajor * rminor + 2.0 * rminor * rminor) return ain, aout, ain + aout + + +def eshellvol(rshell, rmini, rmino, zminor, drin, drout, dz): + """Routine to calculate the inboard, outboard and total volumes + of a toroidal shell comprising two elliptical sections + author: P J Knight, CCFE, Culham Science Centre + rshell : input real : major radius of centre of both ellipses (m) + rmini : input real : horizontal distance from rshell to outer edge + of inboard elliptical shell (m) + rmino : input real : horizontal distance from rshell to inner edge + of outboard elliptical shell (m) + zminor : input real : vertical internal half-height of shell (m) + drin : input real : horiz. thickness of inboard shell at midplane (m) + drout : input real : horiz. thickness of outboard shell at midplane (m) + dz : input real : vertical thickness of shell at top/bottom (m) + vin : output real : volume of inboard section (m3) + vout : output real : volume of outboard section (m3) + vtot : output real : total volume of shell (m3) + This routine calculates the volume of the inboard and outboard sections + of a toroidal shell defined by two co-centred semi-ellipses. + Each section's internal and external surfaces are in turn defined + by two semi-ellipses. The volumes of each section are calculated as + the difference in those of the volumes of revolution enclosed by their + inner and outer surfaces. + """ + # Inboard section + + # Volume enclosed by outer (higher R) surface of elliptical section + # and the vertical straight line joining its ends + a = rmini + b = zminor + elong = b / a + v1 = 2.0 * np.pi * elong * (0.5 * np.pi * rshell * a**2 - (2.0 / 3.0) * a**3) + + # Volume enclosed by inner (lower R) surface of elliptical section + # and the vertical straight line joining its ends + a = rmini + drin + b = zminor + dz + elong = b / a + v2 = 2.0 * np.pi * elong * (0.5 * np.pi * rshell * a**2 - (2.0 / 3.0) * a**3) + + # Volume of inboard section of shell + vin = v2 - v1 + + # Outboard section + + # Volume enclosed by inner (lower R) surface of elliptical section + # and the vertical straight line joining its ends + a = rmino + b = zminor + elong = b / a + v1 = 2.0 * np.pi * elong * (0.5 * np.pi * rshell * a**2 + (2.0 / 3.0) * a**3) + + # Volume enclosed by outer (higher R) surface of elliptical section + # and the vertical straight line joining its ends + a = rmino + drout + b = zminor + dz + elong = b / a + v2 = 2.0 * np.pi * elong * (0.5 * np.pi * rshell * a**2 + (2.0 / 3.0) * a**3) + + # Volume of outboard section of shell + vout = v2 - v1 + + return vin, vout, vin + vout diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 9b6dfc9c46..266c83f5bd 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -23,88 +23,6 @@ module maths_library implicit none contains - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine eshellvol(rshell,rmini,rmino,zminor,drin,drout,dz,vin,vout,vtot) - - !! Routine to calculate the inboard, outboard and total volumes - !! of a toroidal shell comprising two elliptical sections - !! author: P J Knight, CCFE, Culham Science Centre - !! rshell : input real : major radius of centre of both ellipses (m) - !! rmini : input real : horizontal distance from rshell to outer edge - !! of inboard elliptical shell (m) - !! rmino : input real : horizontal distance from rshell to inner edge - !! of outboard elliptical shell (m) - !! zminor : input real : vertical internal half-height of shell (m) - !! drin : input real : horiz. thickness of inboard shell at midplane (m) - !! drout : input real : horiz. thickness of outboard shell at midplane (m) - !! dz : input real : vertical thickness of shell at top/bottom (m) - !! vin : output real : volume of inboard section (m3) - !! vout : output real : volume of outboard section (m3) - !! vtot : output real : total volume of shell (m3) - !! This routine calculates the volume of the inboard and outboard sections - !! of a toroidal shell defined by two co-centred semi-ellipses. - !! Each section's internal and external surfaces are in turn defined - !! by two semi-ellipses. The volumes of each section are calculated as - !! the difference in those of the volumes of revolution enclosed by their - !! inner and outer surfaces. - !!

See also eshellarea - !! Internal CCFE note T&M/PKNIGHT/PROCESS/009, P J Knight: - !! Surface Area and Volume Calculations for Toroidal Shells - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use constants, only: pi, twopi - implicit none - - ! Arguments - real(kind=dp), intent(in) :: rshell, rmini, rmino, zminor, drin, drout, dz - real(kind=dp), intent(out) :: vin, vout, vtot - - ! Local variables - real(kind=dp) :: a, b, elong, v1, v2 - - ! Global shared variables - ! Input: pi,twopi - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! #TODO - Review both equations containing dz and attempt to separate - ! top and bottom of vacuum vessel thickness - ! See issue #433 for explanation - - ! Inboard section - - ! Volume enclosed by outer (higher R) surface of elliptical section - ! and the vertical straight line joining its ends - a = rmini ; b = zminor ; elong = b/a - v1 = twopi * elong * (0.5D0*pi*rshell*a*a - 2.0D0/3.0D0*a*a*a) - - ! Volume enclosed by inner (lower R) surface of elliptical section - ! and the vertical straight line joining its ends - a = rmini+drin ; b = zminor+dz ; elong = b/a - v2 = twopi * elong * (0.5D0*pi*rshell*a*a - 2.0D0/3.0D0*a*a*a) - - ! Volume of inboard section of shell - vin = v2 - v1 - - ! Outboard section - - ! Volume enclosed by inner (lower R) surface of elliptical section - ! and the vertical straight line joining its ends - a = rmino ; b = zminor ; elong = b/a - v1 = twopi * elong * (0.5D0*pi*rshell*a*a + 2.0D0/3.0D0*a*a*a) - - ! Volume enclosed by outer (higher R) surface of elliptical section - ! and the vertical straight line joining its ends - a = rmino+drout ; b = zminor+dz ; elong = b/a - v2 = twopi * elong * (0.5D0*pi*rshell*a*a + 2.0D0/3.0D0*a*a*a) - - ! Volume of outboard section of shell - vout = v2 - v1 - - ! Total shell volume - vtot = vin + vout - - end subroutine eshellvol ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From afce6c58492dd2cb83a9d37e6e974cedf797b3e3 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 11 Dec 2024 11:33:34 +0000 Subject: [PATCH 09/13] Convert dshellvol to Python --- process/blanket_library.py | 52 ++++++++++++++++++++++-- source/fortran/maths_library.f90 | 68 -------------------------------- 2 files changed, 48 insertions(+), 72 deletions(-) diff --git a/process/blanket_library.py b/process/blanket_library.py index 96093bbf03..31fbfa42ec 100644 --- a/process/blanket_library.py +++ b/process/blanket_library.py @@ -15,7 +15,6 @@ error_handling, fwbs_variables, heat_transport_variables, - maths_library, pfcoil_variables, physics_variables, primary_pumping_variables, @@ -194,7 +193,7 @@ def dshaped_component(self, icomponent: int): fwbs_variables.volblkti, fwbs_variables.volblkto, fwbs_variables.volblkt, - ) = maths_library.dshellvol( + ) = dshellvol( r1, r2, blanket_library.hblnkt, @@ -207,7 +206,7 @@ def dshaped_component(self, icomponent: int): blanket_library.volshldi, blanket_library.volshldo, fwbs_variables.volshld, - ) = maths_library.dshellvol( + ) = dshellvol( r1, r2, blanket_library.hshld, @@ -220,7 +219,7 @@ def dshaped_component(self, icomponent: int): blanket_library.volvvi, blanket_library.volvvo, fwbs_variables.vdewin, - ) = maths_library.dshellvol( + ) = dshellvol( r1, r2, blanket_library.hvv, @@ -2829,3 +2828,48 @@ def eshellvol(rshell, rmini, rmino, zminor, drin, drout, dz): vout = v2 - v1 return vin, vout, vin + vout + + +def dshellvol(rmajor, rminor, zminor, drin, drout, dz): + """Routine to calculate the inboard, outboard and total volumes + of a D-shaped toroidal shell + author: P J Knight, CCFE, Culham Science Centre + rmajor : input real : major radius to outer point of inboard + straight section of shell (m) + rminor : input real : horizontal internal width of shell (m) + zminor : input real : vertical internal half-height of shell (m) + drin : input real : horiz. thickness of inboard shell at midplane (m) + drout : input real : horiz. thickness of outboard shell at midplane (m) + dz : input real : vertical thickness of shell at top/bottom (m) + vin : output real : volume of inboard straight section (m3) + vout : output real : volume of outboard curved section (m3) + vtot : output real : total volume of shell (m3) + This routine calculates the volume of the inboard and outboard sections + of a D-shaped toroidal shell defined by the above input parameters. + The inboard section is assumed to be a cylinder of uniform thickness. + The outboard section's internal and external surfaces are defined + by two semi-ellipses, centred on the outer edge of the inboard section; + its volume is calculated as the difference in those of the volumes of + revolution enclosed by the two surfaces. + """ + # Volume of inboard cylindrical shell + vin = 2.0 * (zminor + dz) * np.pi * (rmajor**2 - (rmajor - drin) ** 2) + + # Volume enclosed by inner surface of elliptical outboard section + # and the vertical straight line joining its ends + a = rminor + b = zminor + elong = b / a + v1 = 2.0 * np.pi * elong * (0.5 * np.pi * rmajor * a**2 + (2.0 / 3.0) * a**3) + + # Volume enclosed by outer surface of elliptical outboard section + # and the vertical straight line joining its ends + a = rminor + drout + b = zminor + dz + elong = b / a + v2 = 2.0 * np.pi * elong * (0.5 * np.pi * rmajor * a**2 + (2.0 / 3.0) * a**3) + + # Volume of elliptical outboard shell + vout = v2 - v1 + + return vin, vout, vin + vout diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 266c83f5bd..a4fccc1e5e 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -26,74 +26,6 @@ module maths_library ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine dshellvol(rmajor,rminor,zminor,drin,drout,dz,vin,vout,vtot) - - !! Routine to calculate the inboard, outboard and total volumes - !! of a D-shaped toroidal shell - !! author: P J Knight, CCFE, Culham Science Centre - !! rmajor : input real : major radius to outer point of inboard - !! straight section of shell (m) - !! rminor : input real : horizontal internal width of shell (m) - !! zminor : input real : vertical internal half-height of shell (m) - !! drin : input real : horiz. thickness of inboard shell at midplane (m) - !! drout : input real : horiz. thickness of outboard shell at midplane (m) - !! dz : input real : vertical thickness of shell at top/bottom (m) - !! vin : output real : volume of inboard straight section (m3) - !! vout : output real : volume of outboard curved section (m3) - !! vtot : output real : total volume of shell (m3) - !! This routine calculates the volume of the inboard and outboard sections - !! of a D-shaped toroidal shell defined by the above input parameters. - !! The inboard section is assumed to be a cylinder of uniform thickness. - !! The outboard section's internal and external surfaces are defined - !! by two semi-ellipses, centred on the outer edge of the inboard section; - !! its volume is calculated as the difference in those of the volumes of - !! revolution enclosed by the two surfaces. - !!

See also dshellarea - !! Internal CCFE note T&M/PKNIGHT/PROCESS/009, P J Knight: - !! Surface Area and Volume Calculations for Toroidal Shells - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use constants, only: pi, twopi - implicit none - - ! Arguments - real(kind=dp), intent(in) :: rmajor, rminor, zminor, drin, drout, dz - real(kind=dp), intent(out) :: vin, vout, vtot - - ! Local variables - real(kind=dp) :: a, b, elong, v1, v2 - - ! Global shared variables - ! Input: pi,twopi - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! #TODO - Review both equations containing dz and attempt to separate - ! top and bottom of vacuum vessel thickness - ! See issue #433 for explanation - - ! Volume of inboard cylindrical shell - vin = 2.0D0*(zminor+dz) * pi*(rmajor**2 - (rmajor-drin)**2) - - ! Volume enclosed by inner surface of elliptical outboard section - ! and the vertical straight line joining its ends - a = rminor ; b = zminor ; elong = b/a - v1 = twopi * elong * (0.5D0*pi*rmajor*a*a + 2.0D0/3.0D0*a*a*a) - - ! Volume enclosed by outer surface of elliptical outboard section - ! and the vertical straight line joining its ends - a = rminor+drout ; b = zminor+dz ; elong = b/a - v2 = twopi * elong * (0.5D0*pi*rmajor*a*a + 2.0D0/3.0D0*a*a*a) - - ! Volume of elliptical outboard shell - vout = v2 - v1 - - ! Total shell volume - vtot = vin + vout - - end subroutine dshellvol - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! ------------------------------------------------------------------------ pure function variable_error(variable) real(dp), intent(in) ::variable From fd84c5abfe1ebc5b46aa340f791ce921be1a204e Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 11 Dec 2024 11:41:10 +0000 Subject: [PATCH 10/13] Remove variable_error checking from maths library --- CMakeLists.txt | 1 - source/fortran/constraint_equations.f90 | 5 +-- source/fortran/iteration_variables.f90 | 9 ++---- source/fortran/maths_library.f90 | 42 ------------------------- source/fortran/output.f90 | 4 --- 5 files changed, 3 insertions(+), 58 deletions(-) delete mode 100644 source/fortran/maths_library.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index c6b7786366..9089a71125 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -88,7 +88,6 @@ LIST(APPEND PROCESS_SRCS ife_variables.f90 heat_transport_variables.f90 buildings_variables.f90 - maths_library.f90 iteration_variables.f90 water_usage_variables.f90 constraint_equations.f90 diff --git a/source/fortran/constraint_equations.f90 b/source/fortran/constraint_equations.f90 index 3ab65fdea9..1218847a0e 100755 --- a/source/fortran/constraint_equations.f90 +++ b/source/fortran/constraint_equations.f90 @@ -47,7 +47,6 @@ subroutine constraint_eqns(m,ieqn,cc,con,err,symbol,units) !! !! 1. use numerics, only: icc - use maths_library, only: variable_error implicit none @@ -307,9 +306,7 @@ subroutine constraint_eqns(m,ieqn,cc,con,err,symbol,units) if (present(units)) units(i) = tmp_units end if - ! Crude method of catching NaN errors - !if ((abs(cc(i)) > 9.99D99).or.(cc(i) /= cc(i))) then - if (variable_error(cc(i))) then + if (isnan(cc(i)).or.(abs(cc(i))>9.99D99)) then ! Add debugging lines as appropriate... select case (icc(i)) diff --git a/source/fortran/iteration_variables.f90 b/source/fortran/iteration_variables.f90 index 8a6ef5058e..d535d577a9 100755 --- a/source/fortran/iteration_variables.f90 +++ b/source/fortran/iteration_variables.f90 @@ -3933,7 +3933,6 @@ subroutine loadxc ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use constants, only: nout - use maths_library, only: variable_error use error_handling, only: idiags, fdiags, report_error use numerics, only: nvar, xcm, ixc, name_xc, lablxc, scafc, scale use physics_variables, only: i_plasma_current @@ -4152,10 +4151,7 @@ subroutine loadxc call report_error(55) end if - ! Crude method of catching NaN errors - - !if ( (abs(xcm(i)) > 9.99D99).or.(xcm(i) /= xcm(i)) ) then - if ( variable_error(xcm(i)) ) then + if ((isnan(xcm(i)).or.(abs(xcm(i))>9.99D99))) then idiags(1) = i ; idiags(2) = ixc(i) ; fdiags(1) = xcm(i) call report_error(56) end if @@ -4190,7 +4186,6 @@ subroutine convxc(xc,nn) use error_handling, only: idiags, fdiags, report_error use numerics, only: ipnvars, scale, ixc, lablxc - use maths_library, only: variable_error #ifndef dp use, intrinsic :: iso_fortran_env, only: dp=>real64 #endif @@ -4406,7 +4401,7 @@ subroutine convxc(xc,nn) ! Crude method of catching NaN errors !if ((abs(xc(i)) > 9.99D99).or.(xc(i) /= xc(i))) then - if(variable_error(xc(i)))then + if (isnan(xc(i)).or.(abs(xc(i))>9.99D99)) then idiags(1) = i ; idiags(2) = ixc(i) ; fdiags(1) = xc(i) call report_error(59) end if diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 deleted file mode 100644 index a4fccc1e5e..0000000000 --- a/source/fortran/maths_library.f90 +++ /dev/null @@ -1,42 +0,0 @@ -module maths_library - - !! Library of mathematical and numerical routines - !! author: P J Knight, CCFE, Culham Science Centre - !! N/A - !! This module contains a large number of routines to enable - !! PROCESS to perform a variety of numerical procedures, including - !! linear algebra, zero finding, integration and minimisation. - !! The module is an amalgamation of the contents of several - !! different pre-existing PROCESS source files, which themselves - !! were derived from a number of different numerical libraries - !! including BLAS and MINPAC. - !! http://en.wikipedia.org/wiki/Gamma_function - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - - !use process_output - - implicit none - -contains - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! ------------------------------------------------------------------------ - pure function variable_error(variable) - real(dp), intent(in) ::variable - logical::variable_error - - if((variable/=variable).or.(variable<-9.99D99).or.(variable>9.99D99))then - variable_error = .TRUE. - else - variable_error = .FALSE. - end if - - end function variable_error - -end module maths_library diff --git a/source/fortran/output.f90 b/source/fortran/output.f90 index 5e5200db7f..c40763bb7a 100755 --- a/source/fortran/output.f90 +++ b/source/fortran/output.f90 @@ -422,7 +422,6 @@ subroutine ovarre(file,descr,varnam,value,output_flag) use numerics, only: name_xc use global_variables, only: icase, vlabel use constants, only: mfile, nout - use maths_library, only: variable_error implicit none ! Arguments @@ -502,7 +501,6 @@ subroutine ovarin(file,descr,varnam,value,output_flag) use numerics, only: name_xc, icc, ioptimz use global_variables, only: xlabel_2, iscan_global use constants, only: mfile, nout - use maths_library, only: variable_error implicit none ! Arguments @@ -619,7 +617,6 @@ subroutine ocosts(file,ccode,descr,value) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use constants, only: pi, mfile, nplot, twopi - use maths_library, only: variable_error implicit none ! Arguments @@ -667,7 +664,6 @@ subroutine obuild(file,descr,thick,total,variable_name) use numerics, only: boundl, boundu use constants, only: electron_charge - use maths_library, only: variable_error implicit none ! Arguments From 6ced5926f57ec2e245f48611f584e725662add3b Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 8 Jan 2025 17:20:01 +0000 Subject: [PATCH 11/13] Reverted some new fixes on main --- process/pfcoil.py | 4 ++-- tests/integration/test_pfcoil_int.py | 30 ++++++++++++++-------------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/process/pfcoil.py b/process/pfcoil.py index cab2f3ecdf..cd49aa9571 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -1195,7 +1195,7 @@ def ohcalc(self): # Allowable coil overall current density at EOF # (superconducting coils only) - (jcritwp, pfv.jcableoh_eof, pfv.jscoh_eof, tmarg1) = self.superconpf( + jcritwp, pfv.jcableoh_eof, pfv.jscoh_eof, tmarg1 = self.superconpf( pfv.bmaxoh, pfv.vfohc, pfv.fcuohsu, @@ -1218,7 +1218,7 @@ def ohcalc(self): # Allowable coil overall current density at BOP - (jcritwp, pfv.jcableoh_bop, pfv.jscoh_bop, tmarg2) = self.superconpf( + jcritwp, pfv.jcableoh_bop, pfv.jscoh_bop, tmarg2 = self.superconpf( pfv.bmaxoh0, pfv.vfohc, pfv.fcuohsu, diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py index aeaf748eb4..400d74be5e 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -126,7 +126,7 @@ def test_pfcoil(monkeypatch, pfcoil): monkeypatch.setattr(pv, "vsind", 3.497e2) monkeypatch.setattr(pv, "aspect", 3.1) monkeypatch.setattr(pv, "itart", 0) - monkeypatch.setattr(pv, "betap", 6.313e-1) + monkeypatch.setattr(pv, "beta_poloidal", 6.313e-1) monkeypatch.setattr(tfv, "tftmp", 4.750) monkeypatch.setattr(tfv, "dcond", np.full(9, 9.0e3)) monkeypatch.setattr(tfv, "i_tf_sup", 1) @@ -2430,13 +2430,13 @@ def test_peakb(monkeypatch: pytest.MonkeyPatch, pfcoil: PFCoil): pfv, "curpfb", np.array([ - 0.067422231232391661, - -2.9167273287450968, - -8.1098913365453491, - -8.1098913365453491, - -5.5984385047179153, - -5.5984385047179153, - -186.98751599968148, + 14.742063826112622, + 20.032681634901664, + 0.58040662653667285, + 0.58040662653667285, + 0.42974674788703021, + 0.42974674788703021, + 174.22748790786324, 0, 0, 0, @@ -2486,13 +2486,13 @@ def test_peakb(monkeypatch: pytest.MonkeyPatch, pfcoil: PFCoil): pfv, "curpfs", np.array([ - 14.742063826112622, - 20.032681634901664, - 0.58040662653667285, - 0.58040662653667285, - 0.42974674788703021, - 0.42974674788703021, - 174.22748790786324, + 0.067422231232391661, + -2.9167273287450968, + -8.1098913365453491, + -8.1098913365453491, + -5.5984385047179153, + -5.5984385047179153, + -186.98751599968148, 0, 0, 0, From d92a8ed62de67f295c08fcb8d9dbde9473f7c23f Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 29 Jan 2025 10:59:16 +0000 Subject: [PATCH 12/13] Fix ruff formatting in PF coil --- process/pfcoil.py | 24 ------------------------ 1 file changed, 24 deletions(-) diff --git a/process/pfcoil.py b/process/pfcoil.py index cd49aa9571..2941dbbb2a 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -274,14 +274,6 @@ def pfcoil(self): bfix, gmat, bvec, - rc, - zc, - cc, - xc, - umat, - vmat, - sigma, - work2, ) # Equilibrium coil currents determined by SVD targeting B @@ -426,14 +418,6 @@ def pfcoil(self): bfix, gmat, bvec, - rc, - zc, - cc, - xc, - umat, - vmat, - sigma, - work2, ) for ccount in range(ngrp0): @@ -796,14 +780,6 @@ def efc( bfix, gmat, bvec, - _rc, - _zc, - _cc, - _xc, - umat, - vmat, - sigma, - work2, ): """Calculates field coil currents. From 972dd210e208663f51cd202f68ceb01536256f62 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 5 Feb 2025 15:37:08 +0000 Subject: [PATCH 13/13] Dont remove 4th ccls in test_efc --- tests/integration/test_pfcoil_int.py | 20 ++------------------ 1 file changed, 2 insertions(+), 18 deletions(-) diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py index 400d74be5e..4ceb207b7d 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -281,7 +281,6 @@ def test_efc(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): :type monkeypatch: MonkeyPatch """ ngrpmx = 10 - nclsmx = 2 nptsmx = 32 nfixmx = 64 lrow1 = 2 * nptsmx + ngrpmx @@ -408,14 +407,6 @@ def test_efc(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): bfix = np.full(lrow1, 0.0) gmat = np.full([lrow1, lcol1], 0.0, order="F") bvec = np.full(lrow1, 0.0) - rc = np.full(nclsmx, 0.0) - zc = np.full(nclsmx, 0.0) - cc = np.full(nclsmx, 0.0) - xc = np.full(nclsmx, 0.0) - umat = np.full([lrow1, lcol1], 0.0, order="F") - vmat = np.full([lrow1, lcol1], 0.0, order="F") - sigma = np.full(ngrpmx, 0.0) - work2 = np.full(ngrpmx, 0.0) ssq, ccls = pfcoil.efc( npts, @@ -435,22 +426,15 @@ def test_efc(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): bfix, gmat, bvec, - rc, - zc, - cc, - xc, - umat, - vmat, - sigma, - work2, ) assert pytest.approx(ssq) == 4.208729e-4 - assert ccls[0:3] == pytest.approx( + assert ccls[0:4] == pytest.approx( np.array([ 12846165.42893886, 16377261.02000236, 579111.6216917, + 0.0, ]) )