From b5cda1566f0a26ada450cdc3108e5ce8c604dc6e Mon Sep 17 00:00:00 2001 From: Ashton Bradley Date: Tue, 30 Apr 2024 12:14:55 +1200 Subject: [PATCH] fix types;add core remake script --- .github/workflows/CI.yml | 2 +- Manifest.toml | 12 ++- src/cores.jld2 | Bin 9568 -> 8734 bytes src/create_exact_core.jl | 155 +++++++++++++++++++++++++++++++++++++++ src/creation.jl | 28 ++++--- src/types.jl | 4 +- 6 files changed, 180 insertions(+), 21 deletions(-) create mode 100644 src/create_exact_core.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 83f9fe4..81eff29 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -11,7 +11,7 @@ jobs: matrix: version: # - '1.6' - - '1.8' + - '1.10' os: - ubuntu-latest arch: diff --git a/Manifest.toml b/Manifest.toml index a70cd80..6ea9045 100644 --- a/Manifest.toml +++ b/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" diff --git a/src/cores.jld2 b/src/cores.jld2 index da92d9577a089559941d2b6acf3b514beb1f12bb..8ce84fa9e4a20a50d047e3eb691584b24abc017d 100644 GIT binary patch delta 1116 zcmZ{jYc$(;7{~wrB!aEBQq{RsXi`00Hj~tPT#h|X%d{Aq%T|}N64OK4QN*Ry|Ejlb zrXe{_OIJl@3`c{onO0kCMOj_Oq)I1JcadM*#-W;}OZT$Ro6qw+=X{>?y!pz#i#&4o zQdkc5>{{nGTtKoiwcdgKSA3d|`szMR=s@Z*($mdYP)Y(frhY5A8OK4`OQGRY_4AZT zqr5K-gpW=3<9b4?Ei}sILjPbPZpN&S+hR_p&u`PYjh8%3{SSBqknqpRdL#@Z8{(TV zjD%x2M_^$?WsornW>bG>yH?Q&YQpBhnBMCZLZ^^ly@t8*N*m}-OlZ1Sq87ROfUk34 z*3F+IfHBp6toC3$@NPJ9yF*eyHuQq!_2@LDv3q|vBAo{xu}DVE$_A=SB!50R4=4wm zsv3&&0Sgl}9(OMUMoL{CFZ??A+R3bIGfIKGl(m9g6an#>>n@(SzXJSgtA$Z9)lhDV zUyt7|hUy{x$P{KBuxN_7j>8fl)RnDC4UN!#0;{{QrwK|)lHcfq%^QGUk#j!@Itj>k z^!E0vdjt~_Dvg-!22D|>y&tU?7>jef)!;v&BhpyZU(g4|la3)1GyTvRf2Z5x-~g=N zF!;qoB7?CBRJrZFK^XryHij;K0-UnQTU{^=-^cXi20VER8DcadIXwcy`@XVN90R`Rlx^(H7}okv#{$*Wo2Q< zEUZ~Yk(NW1z}+2}Gx}5s>)=6jjMCYtBq!?buymYE8o@72>)92G2P{&+K> zLIsj8R(isy3L5y+{6u2}<>%PDzja2CxbqA#5{e*vA;XBzLvU;`X=h#z0`7cq;_YGt zCWqQ0q+$eP1Ba%e1;LV~-{XRZ2#yE*>heJ!f0SoIHrl7;LZ`6_~arJTa@H3Ye9`aPSUBS54! zE^`CH&`DV-#CHN?0-CFpm*&e#UQn$)+^St`{aQ&Z0b6az= vHl{CYRpUPJmq56eG5k*8s~c1k48u3$)X3!YgqWlyYvEM}n^ODt*)#tFw6W^` delta 1128 zcmb8uYfMuI6bJBods`G|K@_UUn+6wF$8(U1TaF1h5TPiK35D?xi3HdNQ(g|}KsTJw zwcQOWn012dD6mO+Zxq?WC=3wFD?A*8mY1{zE(8Wd!LD1guO>TRekUjYoNs@jkHVkt ziFrD?)_HZIxE8fCjdei&ogsU@Q5p{i+=VoVAH%5oL(@}J5_Iqw?=aPX9?kVz!bXYXqcJ=MggX;Od->IZ8fZr}0Hz+6q$}pR^eq8{5BW^WK z!V=&}i<&<^QwCAR)q;$u@1QGp*t+gU1sL4U-9`#Uz&ZTVHbXKh28=7hFDLv6@9*Ce zCR<6MGXRT8XVk-0%5+lihmG)wx2WSeS8+HM)t`I;sf3cEoS8GY8IxCb;{>{~y4)5ifO z)7zg_-w(!qT5FC*4ksnqPAr=taDS@G*a;hkJ1A8&S~>!0AubV%>QQKT)G9xIZVc+H zu@;}kaaj9RU48WAByhLl66lg=(B;f-8-OXeKita?el`u86!lDF-VE>);w@%N6wq=x zeAzHn2{m~|IlBRzSt#rP(wl1dK(-Qk=4}ml%!#YZdNuIQ zomkRe>}5FJlE_s~FGC3TQO|aFg&^L8+FL6efu()|q1aWZXte0P-?^?z+yR9fop}_W#8I75k_bm5jc&<3 zA^3D#XZX%+z@E-xadW0t=ZSA~xBs!Kw((v30*%2au+SG}%35#Gj1~Dc?fO4W{4XoF gGRMfZf?_7VX*-2_<3o04Bu{%h!p@B>yTp(A8?P?ve*gdg diff --git a/src/create_exact_core.jl b/src/create_exact_core.jl new file mode 100644 index 0000000..62004a6 --- /dev/null +++ b/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 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 \ No newline at end of file diff --git a/src/creation.jl b/src/creation.jl index fbdc008..301be66 100644 --- a/src/creation.jl +++ b/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 + + diff --git a/src/types.jl b/src/types.jl index e1a7014..4bba8f7 100644 --- a/src/types.jl +++ b/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