diff --git a/.gitmodules b/.gitmodules index f84a08b..43b429d 100644 --- a/.gitmodules +++ b/.gitmodules @@ -10,3 +10,6 @@ [submodule "submodules/gmsh"] path = submodules/gmsh url = https://github.com/live-clones/gmsh.git +[submodule "submodules/SchwarzChristoffel.jl"] + path = submodules/SchwarzChristoffel.jl + url = https://github.com/jdeldre/SchwarzChristoffel.jl.git diff --git a/.vscode/settings.json b/.vscode/settings.json index ea71ab1..07723d3 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -47,6 +47,11 @@ "git.fetchOnPull": true, "git.useCommitInputAsStashMessage": true, "cSpell.words": [ - "eltype" + "eltype", + "myblue", + "normp", + "println", + "xlims", + "ylims" ] } \ No newline at end of file diff --git a/Manifest.toml b/Manifest.toml index 4248c4e..7e9715a 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -450,6 +450,12 @@ git-tree-sha1 = "3ebb967819b284dc1e3c0422229b58a40a255649" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" version = "0.26.3" +[[DocumenterTools]] +deps = ["Base64", "DocStringExtensions", "Documenter", "FileWatching", "LibGit2", "Sass"] +git-tree-sha1 = "fd313c343451ff8f88ff711662ea13f9629b9a58" +uuid = "35a29f4d-8980-5a13-9543-d66fff28ecb8" +version = "0.1.10" + [[DoubleFloats]] deps = ["LinearAlgebra", "Polynomials", "Printf", "Quadmath", "Random", "Requires", "SpecialFunctions"] git-tree-sha1 = "2e43ebdac158a9ccb42aed91f0be658d6bb5f81c" @@ -1348,9 +1354,9 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" [[ProgressMeter]] deps = ["Distributed", "Printf"] -git-tree-sha1 = "6e9c89cba09f6ef134b00e10625590746ba1e036" +git-tree-sha1 = "d85d8f0339a9937afac93e152c76f4745b386202" uuid = "92933f4c-e287-5a05-a399-4b506db050ca" -version = "1.5.0" +version = "1.6.0" [[Qt5Base_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"] @@ -1523,6 +1529,12 @@ git-tree-sha1 = "36ebc5622c82eb9324005cc75e7e2cc51181d181" uuid = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" version = "0.0.1" +[[Sass]] +deps = ["libsass_jll"] +git-tree-sha1 = "aa841c3738cec78b5dbccd56dda332710f35f6a5" +uuid = "322a6be2-4ae8-5d68-aaf1-3e960788d1d9" +version = "0.2.0" + [[SciMLBase]] deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"] git-tree-sha1 = "5d91d6098eaffa03eb371f3f3699f01983ba174e" @@ -2107,6 +2119,12 @@ git-tree-sha1 = "6abbc424248097d69c0c87ba50fcb0753f93e0ee" uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" version = "1.6.37+6" +[[libsass_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "663428b7ebaf60c60ee147f0f9466430e9959ad6" +uuid = "47bcb7c8-5119-555a-9eeb-0afcc36cd728" +version = "3.5.5+0" + [[libsodium_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "848ab3d00fe39d6fbc2a8641048f8f272af1c51e" diff --git a/Project-1/example/demo.jl b/Project-1/example/demo.jl index 39d7b40..a45b9da 100644 --- a/Project-1/example/demo.jl +++ b/Project-1/example/demo.jl @@ -15,11 +15,8 @@ using Pkg Pkg.activate(normpath(joinpath(@__DIR__, "../.."))) # using Revise -using BenchmarkTools -using Plots -using OffsetArrays - -myblue = RGB(0.368417, 0.506779, 0.709798) # <==> 94, 129, 181 (the default style in Wolfram Mathematica) +using BenchmarkTools, ProgressMeter +using Plots, OffsetArrays, LinearAlgebra """ yt0012(ξ) @@ -58,7 +55,7 @@ Base.@pure function yt0012(ξ::T where {T<:Real}) end) end -const Mx, My = 100+1, 20+1 +const Mx, My = 100+1, 100+1 include(normpath(joinpath(@__DIR__, "../../src/misc-util.jl"))) using .__CFD2021__misc_util__: tuplejoin @@ -71,6 +68,8 @@ using .__CFD2021__misc_util__: tuplejoin ``` """ +include(normpath(joinpath(@__DIR__, "../src/visualization.jl"))) + include(normpath(joinpath(@__DIR__, "../src/transfinite-interpolate.jl"))) import .transfinite_interpolate: transfinite_interpolate_2d!, transfinite_interpolate_2d @@ -121,15 +120,20 @@ y_hi = (5sinpi.(LinRange(1, 0, Mx)), [LinRange(0, 0, My)...]); ``` """ -p = plot(; xlims=(-5, 5), ylims=(-5, 5), aspect_ratio=:equal) -# p = plot(xlims = (-.6, .6), ylims = (0, .5), aspect_ratio = :equal) interpolated_x = Array{Float64}(undef, (length(v) for v in x_lo)...); interpolated_y = Array{Float64}(undef, (length(v) for v in y_lo)...); -@btime transfinite_interpolate_2d!(interpolated_x, x_lo, x_hi) -@btime transfinite_interpolate_2d!(interpolated_y, y_lo, y_hi) -for u in axes(interpolated_x)[begin] plot!(p,interpolated_x[u,:], interpolated_y[u,:], label = nothing, color = myblue); end -for u in axes(interpolated_x)[ end ] plot!(p,interpolated_x[:,u], interpolated_y[:,u], label = nothing, color = myblue); end -# for u in axes(interpolated_x)[begin] plot!(p,interpolated_x[u,:],-interpolated_y[u,:], label = nothing, color = myblue); end -# for u in axes(interpolated_x)[ end ] plot!(p,interpolated_x[:,u],-interpolated_y[:,u], label = nothing, color = myblue); end -plot!(p, yt0012, 0, 1; label=nothing, color=:red) -display(p) # savefig(normpath(joinpath(@__DIR__, "img/tf-hO.svg"))) +#= @btime =# transfinite_interpolate_2d!(interpolated_x, x_lo, x_hi) +#= @btime =# transfinite_interpolate_2d!(interpolated_y, y_lo, y_hi) +p = plot(; xlims=(-5,5), ylims=(-5,5), aspect_ratio=:equal) +# p = grid_plot!(p, interpolated_x, interpolated_y; title="tf", mirror=true, xlims=(-5, 5), ylims=(-5, 5), aspect_ratio=:equal) +# savefig(normpath(joinpath(@__DIR__, "img/tf-hO.png"))) + +include(normpath(joinpath(@__DIR__, "../src/elliptic-grid-gen.jl"))) +import .elliptic_grid_gen: inverted_poisson_2d_jacobi_step , inverted_poisson_2d_jacobi_step! , + inverted_poisson_2d_jacobi_iterate, inverted_poisson_2d_jacobi_iterate! +P = Q = ((x, y)->(norm([x, y]-[-3., 3.])) < 1. ? 20. : 0. ) +inverted_poisson_2d_jacobi_iterate!(interpolated_x, interpolated_y, 1e-2, nothing, P, Q, 2.) +p = plot(; xlims=(-5,5), ylims=(-5,5), aspect_ratio=:equal) +p = grid_plot!(p, interpolated_x, interpolated_y; title="el", mirror=true, xlims=(-5, 5), ylims=(-5, 5), aspect_ratio=:equal) +display(p) +# savefig(normpath(joinpath(@__DIR__, "img/el-hO.png"))) diff --git a/Project-1/src/elliptic-grid-gen.jl b/Project-1/src/elliptic-grid-gen.jl index c66c0c1..c34a3b5 100644 --- a/Project-1/src/elliptic-grid-gen.jl +++ b/Project-1/src/elliptic-grid-gen.jl @@ -20,83 +20,182 @@ module elliptic_grid_gen -using LinearAlgebra +using LinearAlgebra, ProgressMeter """ -Work in progress. + inverted_poisson_2d_jacobi_step(xs, ys[, P, Q]) + inverted_poisson_2d_jacobi_step!(output_xs, output_ys, xs, ys[, P, Q]) + +Take a Jacobi step for an array-pair represented 2d grid. + +TODO: generalize to higher dimensions. + +# Arguments +- `xs::Matrix{<: Number}`: array of `x`-coordinates. +- `ys::Matrix{<: Number}`: array of `y`-coordinates. +- `P::Function`: (optional) `ξ` related source term. +- `Q::Function`: (optional) `η` related source term. + +# Output +- `output_xs::Matrix{<: Number}`: array of `x`-coordinates, at the next step. +- `output_ys::Matrix{<: Number}`: array of `y`-coordinates, at the next step. + +# Examples +```jldoctest +``` """ -function inverted_poisson_2d_jacobi_step(xs::Matrix{T}, ys::Matrix{T}) where {T <: Number} +function inverted_poisson_2d_jacobi_step( + xs::Matrix{T} , ys::Matrix{T} , + P::S = nothing, Q::S = nothing, + boundaray_ortho_factor::Union{Nothing, Float64}=nothing) where {T <: Number, S <: Union{Nothing, Function}} + (axes(xs) == axes(ys)) || throw(DimensionMismatch("Dimensions of x ($(axes(xs))) and y ($(axes(ys))) does not match.")) output_xs = similar(xs) - output_ys = similar(xs) - inverted_poisson_2d_jacobi_step!(output_xs, output_ys, xs, ys) + output_ys = similar(ys) + inverted_poisson_2d_jacobi_step!(output_xs, output_ys, xs, ys, P, Q, boundaray_ortho_factor) return output_xs, output_ys end -function inverted_poisson_2d_jacobi_step!(output_xs::Matrix{T}, output_ys::Matrix{T}, xs::Matrix{T}, ys::Matrix{T}) where {T <: Number} +function inverted_poisson_2d_jacobi_step!( + output_xs::Matrix{T}, output_ys::Matrix{T}, + xs::Matrix{T}, ys::Matrix{T}, + P::S = nothing , Q::S = nothing , + boundaray_ortho_factor::Union{Nothing, Float64}=nothing) where {T <: Number, S <: Union{Nothing, Function}} + (axes(output_xs) == axes(output_ys) == axes(xs) == axes(ys)) || throw(DimensionMismatch("Dimensions of input ($((axes(xs), axes(ys))) and output ($((axes(output_xs), axes(output_ys)))) are not compatible")) M = [l-1 for l in size(xs)] - - begin # neighboring x's - x_l(u::CartesianIndex) = u[1] == axes(xs)[1][begin] ? zero(eltype(xs)) : xs[u-CartesianIndex(1,0)] - x_r(u::CartesianIndex) = u[1] == axes(xs)[1][ end ] ? zero(eltype(xs)) : xs[u+CartesianIndex(1,0)] - x_d(u::CartesianIndex) = u[2] == axes(xs)[2][begin] ? zero(eltype(xs)) : xs[u-CartesianIndex(0,1)] - x_u(u::CartesianIndex) = u[2] == axes(xs)[2][ end ] ? zero(eltype(xs)) : xs[u+CartesianIndex(0,1)] - xld(u::CartesianIndex) = u[1] == axes(xs)[1][begin] || u[2] == axes(xs)[2][begin] ? zero(eltype(xs)) : xs[u+CartesianIndex(-1,-1)] - xrd(u::CartesianIndex) = u[1] == axes(xs)[1][ end ] || u[2] == axes(xs)[2][begin] ? zero(eltype(xs)) : xs[u+CartesianIndex( 1,-1)] - xlu(u::CartesianIndex) = u[2] == axes(xs)[1][begin] || u[2] == axes(xs)[2][ end ] ? zero(eltype(xs)) : xs[u+CartesianIndex(-1, 1)] - xru(u::CartesianIndex) = u[2] == axes(xs)[1][ end ] || u[2] == axes(xs)[2][ end ] ? zero(eltype(xs)) : xs[u+CartesianIndex( 1, 1)] + "neighboring x's"; begin + @inline x_l(u::CartesianIndex) = u[1] == axes(xs)[1][begin] ? zero(eltype(xs)) : xs[u+CartesianIndex(-1, 0)] + @inline x_r(u::CartesianIndex) = u[1] == axes(xs)[1][ end ] ? zero(eltype(xs)) : xs[u+CartesianIndex( 1, 0)] + @inline x_d(u::CartesianIndex) = u[2] == axes(xs)[2][begin] ? zero(eltype(xs)) : xs[u+CartesianIndex( 0,-1)] + @inline x_u(u::CartesianIndex) = u[2] == axes(xs)[2][ end ] ? zero(eltype(xs)) : xs[u+CartesianIndex( 0, 1)] + @inline xld(u::CartesianIndex) = u[1] == axes(xs)[1][begin] || u[2] == axes(xs)[2][begin] ? zero(eltype(xs)) : xs[u+CartesianIndex(-1,-1)] + @inline xrd(u::CartesianIndex) = u[1] == axes(xs)[1][ end ] || u[2] == axes(xs)[2][begin] ? zero(eltype(xs)) : xs[u+CartesianIndex( 1,-1)] + @inline xlu(u::CartesianIndex) = u[1] == axes(xs)[1][begin] || u[2] == axes(xs)[2][ end ] ? zero(eltype(xs)) : xs[u+CartesianIndex(-1, 1)] + @inline xru(u::CartesianIndex) = u[1] == axes(xs)[1][ end ] || u[2] == axes(xs)[2][ end ] ? zero(eltype(xs)) : xs[u+CartesianIndex( 1, 1)] + end + "neighboring y's"; begin + @inline y_l(u::CartesianIndex) = u[1] == axes(ys)[1][begin] ? zero(eltype(ys)) : ys[u+CartesianIndex(-1, 0)] + @inline y_r(u::CartesianIndex) = u[1] == axes(ys)[1][ end ] ? zero(eltype(ys)) : ys[u+CartesianIndex( 1, 0)] + @inline y_d(u::CartesianIndex) = u[2] == axes(ys)[2][begin] ? zero(eltype(ys)) : ys[u+CartesianIndex( 0,-1)] + @inline y_u(u::CartesianIndex) = u[2] == axes(ys)[2][ end ] ? zero(eltype(ys)) : ys[u+CartesianIndex( 0, 1)] + @inline yld(u::CartesianIndex) = u[1] == axes(ys)[1][begin] || u[2] == axes(ys)[2][begin] ? zero(eltype(ys)) : ys[u+CartesianIndex(-1,-1)] + @inline yrd(u::CartesianIndex) = u[1] == axes(ys)[1][ end ] || u[2] == axes(ys)[2][begin] ? zero(eltype(ys)) : ys[u+CartesianIndex( 1,-1)] + @inline ylu(u::CartesianIndex) = u[1] == axes(ys)[1][begin] || u[2] == axes(ys)[2][ end ] ? zero(eltype(ys)) : ys[u+CartesianIndex(-1, 1)] + @inline yru(u::CartesianIndex) = u[1] == axes(ys)[1][ end ] || u[2] == axes(ys)[2][ end ] ? zero(eltype(ys)) : ys[u+CartesianIndex( 1, 1)] end - begin # neiboring y's - y_l(u::CartesianIndex) = u[1] == axes(ys)[1][begin] ? zero(eltype(ys)) : ys[u-CartesianIndex(1,0)] - y_r(u::CartesianIndex) = u[1] == axes(ys)[1][ end ] ? zero(eltype(ys)) : ys[u+CartesianIndex(1,0)] - y_d(u::CartesianIndex) = u[2] == axes(ys)[2][begin] ? zero(eltype(ys)) : ys[u-CartesianIndex(0,1)] - y_u(u::CartesianIndex) = u[2] == axes(ys)[2][ end ] ? zero(eltype(ys)) : ys[u+CartesianIndex(0,1)] - yld(u::CartesianIndex) = u[1] == axes(ys)[1][begin] || u[2] == axes(ys)[2][begin] ? zero(eltype(ys)) : ys[u+CartesianIndex(-1,-1)] - yrd(u::CartesianIndex) = u[1] == axes(ys)[1][ end ] || u[2] == axes(ys)[2][begin] ? zero(eltype(ys)) : ys[u+CartesianIndex( 1,-1)] - ylu(u::CartesianIndex) = u[2] == axes(ys)[1][begin] || u[2] == axes(ys)[2][ end ] ? zero(eltype(ys)) : ys[u+CartesianIndex(-1, 1)] - yru(u::CartesianIndex) = u[2] == axes(ys)[1][ end ] || u[2] == axes(ys)[2][ end ] ? zero(eltype(ys)) : ys[u+CartesianIndex( 1, 1)] + "1-order diffs"; begin + @inline ∂xξ(u::CartesianIndex) = M[1]/2 * (x_r(u) - x_l(u)) + @inline ∂xη(u::CartesianIndex) = M[2]/2 * (x_u(u) - x_d(u)) + @inline ∂yξ(u::CartesianIndex) = M[1]/2 * (y_r(u) - y_l(u)) + @inline ∂yη(u::CartesianIndex) = M[2]/2 * (y_u(u) - y_d(u)) + @warn "Currently, boundary orthogonalization is not implemented; boundaray_ortho_factor is silently ignored." maxlog=1 #= + doBO = (typeof(boundaray_ortho_factor) <: Nothing) + doBO || (bof = boundaray_ortho_factor) + @inline ∂xξ(u::CartesianIndex) = doBO && u[1] in axes(xs)[1][[begin:begin+1, end-1 : end ]] ? + bof*(-∂xη(u))/sqrt(∂xη(u)^2+∂yη(u)^2) : ∂xξ_b(u) + @inline ∂xη(u::CartesianIndex) = doBO && u[2] in axes(xs)[2][[begin:begin+1, end-1 : end ]] ? + bof*(-∂xξ(u))/sqrt(∂xξ(u)^2+∂yξ(u)^2) : ∂xη_b(u) + @inline ∂yξ(u::CartesianIndex) = doBO && u[1] in axes(ys)[1][[begin:begin+1, end-1 : end ]] ? + bof*(-∂yη(u))/sqrt(∂xη(u)^2+∂yη(u)^2) : ∂yξ_b(u) + @inline ∂yη(u::CartesianIndex) = doBO && u[2] in axes(ys)[2][[begin:begin+1, end-1 : end ]] ? + bof*(-∂yξ(u))/sqrt(∂xξ(u)^2+∂yξ(u)^2) : ∂yη_b(u) + =# end - begin # 1-order diffs - ∂xξ(u::CartesianIndex) = M[1]/2 * (x_r(u) - x_l(u)) - ∂xη(u::CartesianIndex) = M[2]/2 * (x_u(u) - x_d(u)) - ∂yξ(u::CartesianIndex) = M[1]/2 * (y_r(u) - y_l(u)) - ∂yη(u::CartesianIndex) = M[2]/2 * (y_u(u) - y_d(u)) + "2-order diffs"; begin + @inline α( u::CartesianIndex) = ∂xη(u)^2 + ∂yη(u)^2 + @inline β( u::CartesianIndex) = ∂xξ(u)*∂xη(u) + ∂yξ(u)*∂yη(u) + @inline γ( u::CartesianIndex) = ∂xξ(u)^2 + ∂yξ(u)^2 + @inline bh( u::CartesianIndex) = M[1]^2 * α(u) + @inline bv( u::CartesianIndex) = M[2]^2 * γ(u) + @inline bp( u::CartesianIndex) = (bh(u) + bv(u)) * 2 + @inline cx( u::CartesianIndex) = -M[1]*M[2]/2 * β(u) * (xld(u) - xrd(u) - xlu(u) + xru(u)) + @inline cy( u::CartesianIndex) = -M[1]*M[2]/2 * β(u) * (yld(u) - yrd(u) - ylu(u) + yru(u)) end - begin # 2-order diffs - α(u::CartesianIndex) = ∂xη(u)^2 + ∂yη(u)^2 - β(u::CartesianIndex) = ∂xξ(u)*∂xη(u) + ∂yξ(u)*∂yη(u) - γ(u::CartesianIndex) = ∂xξ(u)^2 + ∂yξ(u)^2 - bh(u::CartesianIndex) = M[1]^2 * α(u) - bv(u::CartesianIndex) = M[2]^2 * γ(u) - bp(u::CartesianIndex) = 2 * (bh(u) + bv(u)) - cx(u::CartesianIndex) = - M[1]*M[2]/2 * β(u) * (xld(u) - xrd(u) - xlu(u) + xru(u)) - cy(u::CartesianIndex) = - M[1]*M[2]/2 * β(u) * (yld(u) - yrd(u) - ylu(u) + yru(u)) + "source related terms"; begin + @inline sx( u::CartesianIndex) = (typeof(P) <: Nothing ? 0 : α(u)*P(xs[u], ys[u])*∂xξ(u)) + + (typeof(Q) <: Nothing ? 0 : γ(u)*Q(xs[u], ys[u])*∂xη(u)) + @inline sy( u::CartesianIndex) = (typeof(P) <: Nothing ? 0 : α(u)*P(xs[u], ys[u])*∂yξ(u)) + + (typeof(Q) <: Nothing ? 0 : γ(u)*Q(xs[u], ys[u])*∂yη(u)) end - for u in CartesianIndices(xs) - output_xs[u] = (bh(u)*(x_l(u) + x_r(u)) + bv(u)*(x_d(u) + x_u(u)) + cx(u)) / bp(u) - output_ys[u] = (bh(u)*(y_l(u) + y_r(u)) + bv(u)*(y_d(u) + y_u(u)) + cy(u)) / bp(u) + "actual computation" + for u in CartesianIndices(xs)[[begin, end], :] + output_xs[u] = xs[u] + output_ys[u] = ys[u]; end + for u in CartesianIndices(xs)[:, [begin, end]] + output_xs[u] = xs[u] + output_ys[u] = ys[u] + end + for u in CartesianIndices(xs)[[begin+1, end-1], begin+1:end-1] + output_xs[u] = (bh(u)*(x_l(u) + x_r(u)) + bv(u)*(x_d(u) + x_u(u)) + cx(u) + sx(u)) / bp(u) + output_ys[u] = (bh(u)*(y_l(u) + y_r(u)) + bv(u)*(y_d(u) + y_u(u)) + cy(u) + sy(u)) / bp(u); end + for u in CartesianIndices(xs)[begin+1:end-1, [begin+1, end-1]] + output_xs[u] = (bh(u)*(x_l(u) + x_r(u)) + bv(u)*(x_d(u) + x_u(u)) + cx(u) + sx(u)) / bp(u) + output_ys[u] = (bh(u)*(y_l(u) + y_r(u)) + bv(u)*(y_d(u) + y_u(u)) + cy(u) + sy(u)) / bp(u); + end + Threads.@threads for u in CartesianIndices(xs)[begin+2:end-2, begin+2:end-2] # this is good enough for a fixed dirichlet boundary + output_xs[u] = (bh(u)*(x_l(u) + x_r(u)) + bv(u)*(x_d(u) + x_u(u)) + cx(u) + sx(u)) / bp(u) + output_ys[u] = (bh(u)*(y_l(u) + y_r(u)) + bv(u)*(y_d(u) + y_u(u)) + cy(u) + sy(u)) / bp(u) end return output_xs, output_ys end """ -Work in progress. + inverted_poisson_2d_jacobi_iterate(xs, ys, ε[, P, Q[, boundaray_ortho_factor]]) + inverted_poisson_2d_jacobi_iterate!(xs, ys, ε[, P, Q[, boundaray_ortho_factor[, cache_xs, cache_ys[, normp=2]]]]) + +Full Jacobi relaxation of an array-pair represented 2d grid. + +TODO: generalize to higher dimensions. + +# Arguments +- `xs::Matrix{<: Number}`: array of `x`-coordinates. +- `ys::Matrix{<: Number}`: array of `y`-coordinates. +- `ε::Float64=1e-10`: (optional) tolerance of residue. +- `P::Function`: (optional) `ξ` related source term. +- `Q::Function`: (optional) `η` related source term. +- `boundaray_ortho_factor::Float64`: (optional) controls grid density at boundaries. +- `cache_xs::Matrix{<: Number}`: (optional) buffer array of `x`-coordinates. +- `cache_ys::Matrix{<: Number}`: (optional) buffer array of `y`-coordinates. +- `normp::Real=2`: use `Lₚ` norm for computing residue. + +# Output +- `output_xs::Matrix{<: Number}`: relaxed array of `x`-coordinates. +- `output_ys::Matrix{<: Number}`: relaxed array of `y`-coordinates. + +# Examples +```jldoctest +``` """ -function inverted_poisson_2d_jacobi_iterate(xs::Matrix{T}, ys::Matrix{T}, ε=1e-10) +function inverted_poisson_2d_jacobi_iterate( + xs::Matrix{T} , ys::Matrix{T} , ε::Float64=1e-10, maxiter::Union{Nothing, Int}=nothing, + P::S = nothing, Q::S = nothing, boundaray_ortho_factor::Union{Nothing, Float64}=nothing, + normp::Real = 2) where {T <: Number, S <: Union{Nothing, Function}} + cache_xs = xs cache_ys = ys - inverted_poisson_2d_jacobi_iterate!(cache_xs, cache_ys) + inverted_poisson_2d_jacobi_iterate!(cache_xs, cache_ys, ε, maxiter, P, Q, boundaray_ortho_factor, normp) return cache_xs, cache_ys end -function inverted_poisson_2d_jacobi_iterate!(xs::Matrix{T}, ys::Matrix{T}, ε=1e-10) - cache_xs = similar(xs) - cache_ys = similar(xs) - residue() = max(norm(cache_xs - xs), norm(cache_ys - ys)) +function inverted_poisson_2d_jacobi_iterate!( + xs::Matrix{T} , ys::Matrix{T} , ε::Float64=1e-10, maxiter::Union{Nothing, Int}=nothing, + P::S = nothing, Q::S = nothing, boundaray_ortho_factor::Union{Nothing, Float64}=nothing, + cache_xs::R = nothing, cache_ys::R = nothing, + normp::Real = 2) where {T <: Number, S <: Union{Nothing, Function}, R <: Union{Nothing, Matrix{T}}} + + !(typeof(cache_xs) <: Nothing) || (cache_xs = similar(xs)) + !(typeof(cache_ys) <: Nothing) || (cache_ys = similar(xs)) + @inline residue() = max(norm(cache_xs - xs, normp), norm(cache_ys - ys, normp)) + @inline steppair()=(inverted_poisson_2d_jacobi_step!(cache_xs, cache_ys, xs, ys, P, Q, boundaray_ortho_factor); + inverted_poisson_2d_jacobi_step!(xs, ys, cache_xs, cache_ys, P, Q, boundaray_ortho_factor)) + steppair() + prog = ProgressThresh(ε, desc="Elliptic relaxation:", showspeed=true) + iter = 0; generate_showvalues(iter) = () -> [(:iter,iter)] while residue() > ε - inverted_poisson_2d_jacobi_step!(cache_xs, cache_ys, xs, ys) - inverted_poisson_2d_jacobi_step!(xs, ys, cache_xs, cache_ys) + (typeof(maxiter) <: Nothing) || iter < maxiter || break + ProgressMeter.update!(prog, residue(); showvalues = generate_showvalues(iter)) + steppair(); iter+=1 end return xs, ys end @@ -109,6 +208,6 @@ function dirichlet_boundary_conform end """ To be done. """ -function nuemann_boundary_conform end +function neumann_boundary_conform end end # module elliptic_grid_gen diff --git a/Project-1/src/transfinite-interpolate.jl b/Project-1/src/transfinite-interpolate.jl index 779a7ad..6173e27 100644 --- a/Project-1/src/transfinite-interpolate.jl +++ b/Project-1/src/transfinite-interpolate.jl @@ -90,7 +90,7 @@ TODO: generalize to higher dimensions. # Arguments - `v_lo::Tuple{Vector{<: Real}, Vector{<: Real}}`: a pair of vectors. - `v_hi::Tuple{Vector{<: Real}, Vector{<: Real}}`: a pair of vectors. -- `v_cr::SMatrix{2, 2, <: Real} = nothing`: (optional) values to take at the corners; +- `v_cr::SMatrix{2, 2, <: Real}`: (optional) values to take at the corners; use average of edge values if not provided. Still, in practical use you should make sure the ends of the edges meet. diff --git a/Project-1/src/visualization.jl b/Project-1/src/visualization.jl new file mode 100644 index 0000000..1c981ca --- /dev/null +++ b/Project-1/src/visualization.jl @@ -0,0 +1,86 @@ +# Copyright 2021 Gravifer +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +using Plots + +""" + yt0012(ξ) + +Compute ``yₜ`` (i.e., the half-thickness relative to the chord) +of NACA0012 airfoils from the normalized chord position ``ξ=x/c``. + +Additionally defined upon ``[1, 1.005]`` for a 2-order smooth closure. + +See also: https://en.wikipedia.org/wiki/NACA_airfoil + +TODO: generalize for 00xx airfoils; furthermore, for all 4-digit airfoils. + +# Examples +```jldoctest +julia> yt0012(0.13) +0.051193437032115506 +``` +""" +Base.@pure function yt0012(ξ::T where {T<:Real}) + ε = 0.0050544107511712 + return 5 * 0.12 * ( + if 0 <= ξ <= 1 + 0.2969 * sqrt(ξ) - + 0.1260 * ξ - + 0.3516 * ξ^2 + + 0.2843 * ξ^3 - + 0.1015 * ξ^4 + elseif 1 < ξ <= 1+ε # 2-order smooth joining of the tail edge. + 0.03423616658680458 * sqrt(1+ε - ξ) + + 0.12491972188290168 * (1+ε - ξ) + + 11.579398976610927 * (1+ε - ξ)^2 + + 12.21202728162769 * (1+ε - ξ)^3 + else + 0 + end) +end + +myblue = RGB(0.368417, 0.506779, 0.709798) # <==> 94, 129, 181 (the default style in Wolfram Mathematica) + +function grid_plot( + xs::Matrix{T} , ys::Matrix{T}; + xlims=(-5, 5), ylims=(-5, 5), aspect_ratio=:equal, color=myblue, title=nothing, mirror::Bool=false) where {T <: Number} + + (axes(xs) == axes(ys)) || + throw(DimensionMismatch("Dimensions of x ($(axes(xs))) and y ($(axes(ys))) does not match.")) + p = plot(; xlims=xlims, ylims=ylims, aspect_ratio=aspect_ratio) + p = grid_display!(p, xs, ys; xlims=xlims, ylims=ylims, aspect_ratio=aspect_ratio, color=color, title=title, mirror=mirror) + return p +end + +function grid_plot!(plot_ref, + xs::Matrix{T} , ys::Matrix{T}; + xlims=(-5, 5), ylims=(-5, 5), aspect_ratio=:equal, color=myblue, title=nothing, mirror::Bool=false) where {T <: Number} + + (axes(xs) == axes(ys)) || + throw(DimensionMismatch("Dimensions of x ($(axes(xs))) and y ($(axes(ys))) does not match.")) + p = (typeof(plot_ref) <: Nothing) ? plot(; xlims=xlims, ylims=ylims, aspect_ratio=aspect_ratio) : plot_ref + "normal"; begin + for u in axes(xs)[begin] plot!(p, interpolated_x[u,:], interpolated_y[u,:], label=nothing, color=color); end + for u in axes(ys)[ end ] plot!(p, interpolated_x[:,u], interpolated_y[:,u], label=nothing, color=color); end + end + if mirror + for u in axes(xs)[begin] plot!(p, interpolated_x[u,:],-interpolated_y[u,:], label=nothing, color=color); end + for u in axes(ys)[ end ] plot!(p, interpolated_x[:,u],-interpolated_y[:,u], label=nothing, color=color); end + end + plot!(p, yt0012 , 0, 1; label=nothing, color=:black) + plot!(p, (ξ)-> -yt0012(ξ), 0, 1; label=nothing, color=:black) + title!(p, title) + return p +end diff --git a/Project.toml b/Project.toml index 66425dd..60a243a 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" @@ -40,6 +41,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" PlotThemes = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a" Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Query = "1a8c2f83-1ff3-5112-b086-8aa67b057ba1" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ReferenceFrameRotations = "74f56ac7-18b3-5285-802d-d4bd4f104033" diff --git a/deps/build.jl b/deps/build.jl index acd4cf9..dc3695b 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -12,4 +12,4 @@ # See the License for the specific language governing permissions and # limitations under the License. using Dates -println(now(), "\tEmpty build of the project") \ No newline at end of file +println(now(), "\tEmpty build of the project") diff --git "a/docs/references/\346\244\255\345\234\206\345\236\213\345\201\217\345\276\256\345\210\206\346\226\271\347\250\213\345\212\240\345\257\206\347\275\221\346\240\274\347\232\204\347\224\237\346\210\220\345\217\212\345\256\236\347\216\260.pdf" "b/docs/references/\346\244\255\345\234\206\345\236\213\345\201\217\345\276\256\345\210\206\346\226\271\347\250\213\345\212\240\345\257\206\347\275\221\346\240\274\347\232\204\347\224\237\346\210\220\345\217\212\345\256\236\347\216\260.pdf" new file mode 100644 index 0000000..0ca0372 Binary files /dev/null and "b/docs/references/\346\244\255\345\234\206\345\236\213\345\201\217\345\276\256\345\210\206\346\226\271\347\250\213\345\212\240\345\257\206\347\275\221\346\240\274\347\232\204\347\224\237\346\210\220\345\217\212\345\256\236\347\216\260.pdf" differ diff --git a/submodules/submodules/SchwarzChristoffel.jl b/submodules/SchwarzChristoffel.jl similarity index 100% rename from submodules/submodules/SchwarzChristoffel.jl rename to submodules/SchwarzChristoffel.jl