diff --git a/CMakeLists.txt b/CMakeLists.txt index c6b778636..9089a7112 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/process/availability.py b/process/availability.py index 3690d24af..558107843 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/process/blanket_library.py b/process/blanket_library.py index 953b854c4..31fbfa42e 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, @@ -180,13 +179,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: @@ -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, @@ -270,13 +269,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: @@ -284,7 +283,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 +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 index d33b42d8f..8578d0fe8 100644 --- 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 diff --git a/process/pfcoil.py b/process/pfcoil.py index f8ed6928c..2941dbbb2 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -4,6 +4,8 @@ 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 from process import fortran as ft @@ -13,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 @@ -273,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 @@ -425,14 +418,6 @@ def pfcoil(self): bfix, gmat, bvec, - rc, - zc, - cc, - xc, - umat, - vmat, - sigma, - work2, ) for ccount in range(ngrp0): @@ -795,14 +780,6 @@ def efc( bfix, gmat, bvec, - _rc, - _zc, - _cc, - _xc, - umat, - vmat, - sigma, - work2, ): """Calculates field coil currents. @@ -895,7 +872,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( @@ -974,13 +951,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 @@ -989,15 +963,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. @@ -1198,7 +1171,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, @@ -1221,7 +1194,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, @@ -1599,13 +1572,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/process/stellarator.py b/process/stellarator.py index d549b477e..d8af56740 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/constraint_equations.f90 b/source/fortran/constraint_equations.f90 index 3ab65fdea..1218847a0 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 8a6ef5058..d535d577a 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 766e6d8b6..000000000 --- a/source/fortran/maths_library.f90 +++ /dev/null @@ -1,863 +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 - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - 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 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 - ! 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) - - !! 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 - !! 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 - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - 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 - -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - 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 - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - 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 - 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 - - ! ------------------------------------------------------------------------ - 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 diff --git a/source/fortran/output.f90 b/source/fortran/output.f90 index 5e5200db7..c40763bb7 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 diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py index 7b2183894..4ceb207b7 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,23 +426,17 @@ 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 pytest.approx(ccls[0:4]) == np.array([ - 12846165.42893886, - 16377261.02000236, - 579111.6216917, - 20660782.82356247, - ]) + assert ccls[0:4] == pytest.approx( + np.array([ + 12846165.42893886, + 16377261.02000236, + 579111.6216917, + 0.0, + ]) + ) def test_mtrx(pfcoil: PFCoil): @@ -1642,31 +1627,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): diff --git a/tests/unit/test_maths_library.py b/tests/unit/test_maths_library.py deleted file mode 100644 index ceec50858..000000000 --- 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