Replace maths library with Python libraries (#3418)

* Use scipy elliptic integrals * Use scipy comb over binomial * Interpolate instead of find_y_nonuniform_x * Use svd from numpy * Convert eshellarea to Python * Convert dshellarea to Python * Move integer to string routines to main module * Convert eshellvol to Python * Convert dshellvol to Python * Remove variable_error checking from maths library * Reverted some new fixes on main * Fix ruff formatting in PF coil * Dont remove 4th ccls in test_efc
ukaea · Feb 12, 2025 · 6537260 · 6537260
1 parent 6bac4f3
commit 6537260
Show file tree

Hide file tree

Showing 12 changed files with 215 additions and 1,062 deletions.
diff --git 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/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/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,
@@ -180,21 +179,21 @@ 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:
             (
                 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,
@@ -270,21 +269,21 @@ 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:
             (
                 fwbs_variables.volblkti,
                 fwbs_variables.volblkto,
                 fwbs_variables.volblkt,
-            ) = maths_library.eshellvol(
+            ) = eshellvol(
                 r1,
                 r2,
                 r3,
@@ -298,7 +297,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 +311,7 @@ def elliptical_component(self, icomponent: int):
                 blanket_library.volvvi,
                 blanket_library.volvvo,
                 fwbs_variables.vdewin,
-            ) = maths_library.eshellvol(
+            ) = eshellvol(
                 r1,
                 r2,
                 r3,
@@ -2710,3 +2709,167 @@ 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
+
+
+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
+
+
+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
+
+
+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/process/build.py b/process/build.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 
+from process.blanket_library import dshellarea, eshellarea
 from process.fortran import (
     blanket_library,
     build_variables,
@@ -11,7 +12,6 @@
     divertor_variables,
     error_handling,
     fwbs_variables,
-    maths_library,
     numerics,
     pfcoil_variables,
     physics_variables,
@@ -2001,15 +2001,12 @@ 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,
                 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
@@ -2038,15 +2035,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