Merge pull request #23 from ProjectTorreyPines/dep_update

Updating dependency Interpolations.kl

Project.toml

@@ -9,7 +9,6 @@ Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
 EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17"
 GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043"

src/SD4SOLPS.jl

@@ -115,7 +115,7 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1)
     psin1d = (psi .- g.simag) ./ (g.sibry - g.simag)
     gq.magnetic_axis.b_field_tor = g.bcentr * g.rcentr / g.rmaxis
     gq.q_axis = g.qpsi[1]
-    gq.q_95 = Interpolations.LinearInterpolation(psin1d, g.qpsi)(0.95)
+    gq.q_95 = Interpolations.linear_interpolation(psin1d, g.qpsi)(0.95)
     qmin_idx = argmin(abs.(g.qpsi))
     gq.q_min.value = g.qpsi[qmin_idx]
     if hasproperty(g, :rhovn)
@@ -202,15 +202,15 @@ function core_profile_2d(dd, prof_time_idx, eq_time_idx, quantity)
     # r and z : coordinates of output points where values of p are desired
 
     # psi_at_requested_points =
-    #     Interpolations.LinearInterpolation((r_eq, z_eq), psinrz).(r, z)
-    # rhonpsi = Interpolations.LinearInterpolation(psin_eq_ext, rhon_eq_ext)
+    #     Interpolations.linear_interpolation((r_eq, z_eq), psinrz).(r, z)
+    # rhonpsi = Interpolations.linear_interpolation(psin_eq_ext, rhon_eq_ext)
     # rho_at_requested_points = rhonpsi.(psi_at_requested_points)
-    # itp = Interpolations.LinearInterpolation(rho_prof_ext, p_ext)
+    # itp = Interpolations.linear_interpolation(rho_prof_ext, p_ext)
     # p_at_requested_points = itp.(rho_at_requested_points)
     # return p_at_requested_points
-    rz2psin = Interpolations.LinearInterpolation((r_eq, z_eq), psinrz)
-    psin2rhon = Interpolations.LinearInterpolation(psin_eq_ext, rhon_eq_ext)
-    rhon2prof = Interpolations.LinearInterpolation(rho_prof_ext, p_ext)
+    rz2psin = Interpolations.linear_interpolation((r_eq, z_eq), psinrz)
+    psin2rhon = Interpolations.linear_interpolation(psin_eq_ext, rhon_eq_ext)
+    rhon2prof = Interpolations.linear_interpolation(rho_prof_ext, p_ext)
     return (r, z) -> rhon2prof(psin2rhon(rz2psin(r, z)))
 end
 

src/actuator_model.jl

@@ -192,7 +192,7 @@ function model_gas_valve(
             prepend!(t_ext, minimum(t) - config["delay"])
             flow0_ext = copy(flow0)
             prepend!(flow0_ext, flow0[1])
-            interp = Interpolations.LinearInterpolation(t_ext, flow0_ext)
+            interp = Interpolations.linear_interpolation(t_ext, flow0_ext)
             delayed_flow = interp.(t .- config["delay"])
             return lowpass_filter(t, delayed_flow, config["tau"])
         end

src/repair_eq.jl

@@ -134,8 +134,8 @@ function add_rho_to_equilibrium!(dd::OMAS.dd)
                 # Vacuum field simplification: B = R0 B0 / R
                 Interpolations.deduplicate_knots!(r1)
                 Interpolations.deduplicate_knots!(r2)
-                z1i = Interpolations.LinearInterpolation(r1, z1)
-                z2i = Interpolations.LinearInterpolation(r2, z2)
+                z1i = Interpolations.linear_interpolation(r1, z1)
+                z2i = Interpolations.linear_interpolation(r2, z2)
                 rr = LinRange(rmin, rmax, 101)
                 rc = (rr[1:end-1] + rr[2:end]) / 2.0
                 integral_part_ = [

src/supersize_profile.jl

@@ -5,13 +5,33 @@ Utilities for extrapolating profiles
 # import CalculusWithJulia
 using OMAS: OMAS
 using Interpolations: Interpolations
-using NumericalIntegration: NumericalIntegration
 using GGDUtils: GGDUtils
 
 export extrapolate_core
 export fill_in_extrapolated_core_profile
 export mesh_psi_spacing
 
+"""
+    cumul_integrate(x::AbstractVector, y::AbstractVector)
+
+Computes cumulative integral of y(x) using trapezoidal rule.
+Source code from obsolete NumericalIntegration.jl package.
+https://github.com/dextorious/NumericalIntegration.jl/blob/540b6c4bee089dfef7b9ae46e8ff188382e7c42e/src/NumericalIntegration.jl#L290
+Modified the code to remove use of @inbounds, @fastmath macros that Julia documentation
+recommends to use with caution.
+"""
+function cumul_integrate(x::AbstractVector, y::AbstractVector)
+    init = (x[2] - x[1]) * (y[1] + y[2])
+    n = length(x)
+    retarr = Vector{typeof(init)}(undef, n)
+    retarr[1] = init
+    for i ∈ 2:n
+        retarr[i] = retarr[i-1] + (x[i] - x[i-1]) * (y[i] + y[i-1])
+    end
+
+    return 0.5 * retarr
+end
+
 """
     extrapolate_core(edge_rho, edge_quantity, rho_output)
 
@@ -43,22 +63,22 @@ function extrapolate_core(edge_rho, edge_quantity, rho_output)
     gg = [0, gmid, gped_enforce, gf]
     rr = [0, rmid, rped_enforce, rf]
     # https://www.matecdev.com/posts/julia-interpolation.html
-    itp = Interpolations.LinearInterpolation(
+    itp = Interpolations.linear_interpolation(
         rr,
         gg;
         extrapolation_bc=Interpolations.Line(),
     )
     #itp = Interpolations.extrapolate(itp, Interpolations.Line())
     g = itp(rho_output)
-    q_extend = NumericalIntegration.cumul_integrate(rho_output, g)
+    q_extend = cumul_integrate(rho_output, g)
     q_offset =
-        edge_quantity[1] - Interpolations.LinearInterpolation(rho_output, q_extend)(rf)
+        edge_quantity[1] - Interpolations.linear_interpolation(rho_output, q_extend)(rf)
     q_extend .+= q_offset
 
     output_profile = Array{Float64}(undef, length(rho_output))
     output_profile[rho_output.<rf] = q_extend[rho_output.<rf]
     output_profile[rho_output.>=rf] =
-        Interpolations.LinearInterpolation(
+        Interpolations.linear_interpolation(
             edge_rho,
             edge_quantity,
         ).(rho_output[rho_output.>=rf])
@@ -180,7 +200,7 @@ function fill_in_extrapolated_core_profile!(
         psi1_eq = (eq_time_slice.profiles_1d.psi .- psia) ./ (psib - psia)
         psi2_eq = (eq_prof_2d.psi .- psia) ./ (psib - psia)
         # println(size(psi2_eq), ", ", size(r_eq), ", ", size(z_eq))
-        rzpi = Interpolations.LinearInterpolation((r_eq, z_eq), psi2_eq)
+        rzpi = Interpolations.linear_interpolation((r_eq, z_eq), psi2_eq)
         in_bounds =
             (r .< maximum(r_eq)) .& (r .> minimum(r_eq)) .& (z .> minimum(z_eq)) .&
             (z .< maximum(z_eq))
@@ -195,7 +215,7 @@ function fill_in_extrapolated_core_profile!(
         prepend!(dpsi, [0.0])
         prepend!(drho, [0.0])
         rho_for_quantity[in_bounds] =
-            Interpolations.LinearInterpolation(
+            Interpolations.linear_interpolation(
                 psi1_eq,
                 rho1_eq,
             ).(psi_for_quantity[in_bounds])
@@ -322,7 +342,7 @@ function mesh_psi_spacing(
     psia = eqt.global_quantities.psi_axis
     psib = eqt.global_quantities.psi_boundary
     psin_eq = (psi .- psia) ./ (psib - psia)
-    rzpi = Interpolations.LinearInterpolation((r_eq, z_eq), psin_eq)
+    rzpi = Interpolations.linear_interpolation((r_eq, z_eq), psin_eq)
     println(minimum(r_eq), ", ", maximum(r_eq))
     println(minimum(z_eq), ", ", maximum(z_eq))
 

test/runtests.jl

@@ -369,7 +369,7 @@ if args["geqdsk_to_imas"]
             psin2d = (p2.psi .- gq.psi_axis) ./ (gq.psi_boundary - gq.psi_axis)
             tolerance = 2.0e-3  # It's not always a high res contour so cut some slack
             psin_bry =
-                Interpolations.LinearInterpolation((r_eq, z_eq), psin2d).(r_bry, z_bry)
+                Interpolations.linear_interpolation((r_eq, z_eq), psin2d).(r_bry, z_bry)
             @test maximum(psin_bry) < (1.0 + tolerance)
             @test minimum(psin_bry) > (1.0 - tolerance)