Use svd from numpy

ukaea · Dec 11, 2024 · 5bf2e20 · 5bf2e20
1 parent 60571f3
commit 5bf2e20
Show file tree

Hide file tree

Showing 3 changed files with 19 additions and 459 deletions.
diff --git a/process/pfcoil.py b/process/pfcoil.py
@@ -8,7 +8,6 @@
 from process.fortran import fwbs_variables as fwbsv
 from process.fortran import constants
 from process.fortran import cs_fatigue_variables as csfv
-from process.fortran import maths_library as ml
 from process.fortran import process_output as op
 from process.fortran import numerics
 from process.fortran import rebco_variables as rcv
@@ -23,6 +22,7 @@
 import logging
 from scipy import optimize
 from scipy.special import ellipk, ellipe
+from scipy.linalg import svd
 
 logger = logging.getLogger(__name__)
 
@@ -892,7 +892,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(
@@ -965,13 +965,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
@@ -980,15 +977,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.
@@ -1189,7 +1185,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,
@@ -1212,7 +1208,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,