fix types;add core remake script

.github/workflows/CI.yml

@@ -11,7 +11,7 @@ jobs:
       matrix:
         version:
           # - '1.6'
-          - '1.8'
+          - '1.10'
         os:
           - ubuntu-latest
         arch:

Manifest.toml

@@ -333,10 +333,16 @@ deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.Interpolations]]
-deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
-git-tree-sha1 = "b7bc05649af456efc75d178846f47006c2c4c3c7"
+deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
+git-tree-sha1 = "88a101217d7cb38a7b481ccd50d21876e1d1b0e0"
 uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-version = "0.13.6"
+version = "0.15.1"
+
+    [deps.Interpolations.extensions]
+    InterpolationsUnitfulExt = "Unitful"
+
+    [deps.Interpolations.weakdeps]
+    Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [[deps.InvertedIndices]]
 git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038"

src/create_exact_core.jl

@@ -0,0 +1,155 @@
+## upgrade of Interpolations or JLD2 may necessitate rerunning 
+# to rebuild interpolation in new package environment
+# Make this a build script?
+
+using Interpolations, JLD2, LinearAlgebra, ToeplitzMatrices
+Λ = 0.8249
+
+#Chebyshev methods
+function gpecore_exact(K,L=2,N=100,R = K)
+    #currently r does  nothing!
+    #N = 100
+    #L = 2
+    #Κ = 1
+    #R = 2 # Stretching coordinate.  R ~ kappa seems to be a good ballpark
+
+    # Convert Chebyshev grid to [0,L] from [1 -1];
+    z,Dz = getChebDMatrix(N)
+    blank,D2z = getChebD2Matrix(N)
+    z = vec(reverse(z,dims=2))
+    z = (z)*L./2
+    Dz = -2*Dz./L
+    D2z = 4*D2z./L.^2
+
+    # y is the physical coordinate, range [0 infinity] (i.e y = r)
+    y =  @. R*(1+z)/(1-z)
+
+    #initial guess based on ansatz
+    ψ = @. y/hypot(1/Λ,y)
+    ψ[1] = 0
+    ψ[end] = 1
+    ψ0 = vec(ψ)
+
+    Q = vec(z.^2 .- 2*z .+ 1)
+    Qmat = repeat(Q,1,N)
+    Zmat = repeat(2*(z .-1),1,N)
+    Ymat = repeat(1 ./y,1,N);
+
+    # Second Derivative
+    residuals = -0.5*( (Q.*(D2z*ψ) + 2*(z .-1).*(Dz*ψ) ).*Q/(4*R^2)
+            + (Q/(2*R) ./y).*(Dz*ψ) )+ 0.5*K^2 *ψ ./y.^2 + ψ.^3 .- ψ
+    residuals[1] = 0; residuals[end] =0
+
+    while sum(abs.(residuals).^2) > 1e-12
+        # Second Derivative
+        residuals = -0.5*( (Q.*(D2z*ψ) + 2*(z .-1).*(Dz*ψ) ).*Q/(4*R^2)
+                + (Q/(2*R) ./y).*(Dz*ψ) )+ 0.5*K^2 *ψ ./y.^2 + ψ.^3 .- ψ
+        residuals[1] = 0; residuals[end] =0
+
+        Jacobi = -0.5*( (Qmat.*(D2z) + Zmat.*(Dz) ).*(Qmat/(4*R^2))
+        + (Qmat/(2*R).*Ymat).*(Dz) ) + diagm(0 => 0.5*K^2 ./y.^2)+ diagm(0 => 3*ψ.^2 .- 1)
+
+
+        Jacobi[1,:] = [1 zeros(1,N-1)]
+        Jacobi[N,:] = [zeros(1,N-1) 1]
+
+
+        Δ = vec(-4/7*(Jacobi\residuals))
+
+        ψ = ψ + Δ
+        ψ[1] = 0
+        ψ[end] =1
+    end
+    res = norm(residuals)^2
+    return y,ψ,res
+end
+
+"""
+`x, DM =  chebdif(N,M)`
+
+computes the differentiation matrices `D1, D2, ..., DM` on Chebyshev nodes.
+
+Input:
+
+`N`         Size of differentiation matrix.
+
+`M`         Number of derivatives required (integer).
+
+Note      `0 < M <= N-1`.
+
+Output:
+
+DM        `DM(1:N,1:N,ell)` contains ell-th derivative matrix, ell=1..M.
+
+The code implements two strategies for enhanced
+accuracy suggested by W. Don and S. Solomonoff in
+SIAM J. Sci. Comp. Vol. 6, pp. 1253--1268 (1994).
+The two strategies are (a) the use of trigonometric
+identities to avoid the computation of differences
+x(k)-x(j) and (b) the use of the "flipping trick"
+which is necessary since sin t can be computed to high
+relative precision when t is small whereas sin (pi-t) cannot.
+Note added May 2003:  It may, in fact, be slightly better not to
+implement the strategies (a) and (b).   Please consult the following
+paper for details:   "Spectral Differencing with a Twist", by
+R. Baltensperger and M.R. Trummer, to appear in SIAM J. Sci. Comp.
+
+J.A.C. Weideman, S.C. Reddy 1998.  Help notes modified by
+JACW, May 2003."""
+function chebdif(N, M)
+   @assert 0<M<=N-1
+
+   n1,n2 = (floor(N/2)|> Int,ceil(N/2)|> Int)  # Indices used for flipping trick.
+
+   k = collect(0.0:N-1)' # Compute theta vector.
+   th = k*π/(N-1)
+
+   x = sin.(π*collect(N-1:-2:1-N)'./(2*(N-1))) # Compute Chebyshev points.
+
+   T = repeat(th/2,N,1)'
+   DX = 2*sin.(T'.+T).*sin.(T'.-T) # Trigonometric identity.
+   DX = [DX[1:n1,:]; -reverse(reverse(DX[1:n2,:],dims=2),dims=1)] # Flipping trick.
+   DX[diagind(DX)] .= 1 # Put 1's on the main diagonal of DX.
+
+   C = Matrix(Toeplitz((-1).^k',(-1).^k'))  # C is the matrix with
+   C[1,:] = C[1,:]*2; 
+   C[N,:] = C[N,:]*2  # entries c(k)/c(j)
+   C[:,1] = C[:,1]/2; 
+   C[:,N] = C[:,N]/2;
+
+   Z = 1 ./DX  # Z contains entries 1/(x(k)-x(j))
+   Z[diagind(Z)] .= 0  # with zeros on the diagonal.
+
+   D = Matrix{Float64}(I, N, N)  # D contains diff. matrices.
+
+   DM = zeros(N,N,M)
+
+   for ell = 1:M
+       D = ell*Z.*(C.*repeat(diag(D),1,N) - D) # Off-diagonals
+       D[diagind(D)] = -sum(D',dims=1)  # Correct main diagonal of D
+       DM[:,:,ell] = D  # Store current D in DM
+   end
+
+   return x, DM
+end
+
+function getChebDMatrix(n)
+   z, M = chebdif(n, 1)
+   Dz = M[:,:,1]
+   return z, Dz
+end
+
+function getChebD2Matrix(n)
+   z, M = chebdif(n, 2)
+   D2z = M[:,:,2]
+   return z, D2z
+end
+
+y,ψ,res = gpecore_exact(1)
+ψi = interpolate((y,),ψ,Gridded(Linear()))
+
+## ansatz core 
+scalar_ansatz(x) = sqrt(x^2/(1+x^2))
+ψa = interpolate((y,),scalar_ansatz.(y),Gridded(Linear()))
+
+@save joinpath(@__DIR__,"cores.jld2") ψi ψa

src/creation.jl

@@ -213,22 +213,7 @@ function periodic_dipole!(psi::F,dip::Vector{ScalarVortex{T}}) where {T <: CoreS
     @pack! psi = ψ
 end
 
-
-
-
-# TODO: dispatch on dipole type
-# function periodic_dipole!(psi::F,dip::Dipole) where F <: Field
-#     @unpack ψ,x,y = psi
-#     rp = vortex_array(dip.vp)[1:2]
-#     rn = vortex_array(dip.vn)[1:2]
-#     @. ψ *= abs(dip.vp(x,y')*dip.vn(x,y'))
-#     ψ .*= exp.(im*dipole_phase(x,y,rp...,rn...))
-#     @pack! psi = ψ
-# end
-
-
 #Chebyshev methods
-
 function gpecore_exact(K,L=2,N=100,R = K)
     #currently r does  nothing!
     #N = 100
@@ -367,3 +352,16 @@ function getChebD2Matrix(n)
    D2z = M[:,:,2]
    return z, D2z
 end
+
+
+# TODO: dispatch on dipole type
+# function periodic_dipole!(psi::F,dip::Dipole) where F <: Field
+#     @unpack ψ,x,y = psi
+#     rp = vortex_array(dip.vp)[1:2]
+#     rn = vortex_array(dip.vn)[1:2]
+#     @. ψ *= abs(dip.vp(x,y')*dip.vn(x,y'))
+#     ψ .*= exp.(im*dipole_phase(x,y,rp...,rn...))
+#     @pack! psi = ψ
+# end
+
+

src/types.jl

@@ -46,13 +46,13 @@ mutable struct PointVortex <: Vortex
 end
 
 struct Ansatz <: CoreShape
-    f::Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}
+    f::Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}
     ξ::Float64
     Λ::Float64
 end
 
 struct Exact <: CoreShape
-    f::Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}
+    f::Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}
     ξ::Float64
 end