From de24e8173047ca0448ccdae6cb4fefdb2a01cf11 Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Wed, 28 Apr 2021 14:06:51 +0300 Subject: [PATCH] Fix issue#18 --- Manifest.toml | 16 ++++++ Project.toml | 4 +- docs/src/man/log.md | 4 -- src/plot/pyplot.jl | 124 +++++++++++++++++++++----------------------- src/plot/utility.jl | 40 +++++++------- 5 files changed, 99 insertions(+), 89 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index f016e827..e07a1f52 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -27,6 +27,10 @@ git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" version = "0.12.8" +[[CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" + [[Conda]] deps = ["JSON", "VersionParsing"] git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421" @@ -43,6 +47,18 @@ version = "1.1.0" deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +[[Dierckx]] +deps = ["Dierckx_jll"] +git-tree-sha1 = "5fefbe52e9a6e55b8f87cb89352d469bd3a3a090" +uuid = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" +version = "0.5.1" + +[[Dierckx_jll]] +deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"] +git-tree-sha1 = "a580560f526f6fc6973e8bad2b036514a4e3b013" +uuid = "cd4c43a9-7502-52ba-aa6d-59fb2a88580b" +version = "0.0.1+0" + [[Downloads]] deps = ["ArgTools", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" diff --git a/Project.toml b/Project.toml index f73c1a36..f82cad66 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "Batsrus" uuid = "e74ebddf-6ac1-4047-a0e5-c32c99e57753" authors = ["Hongyang Zhou "] -version = "0.3.1" +version = "0.3.2" [deps] +Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" @@ -16,6 +17,7 @@ UnitfulRecipes = "42071c24-d89e-48dd-8a24-8a12d9b8861f" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] +Dierckx = "0.5" FortranFiles = "0.6" Glob = "1.3" LightXML = "0.9" diff --git a/docs/src/man/log.md b/docs/src/man/log.md index cd115183..d8a001aa 100644 --- a/docs/src/man/log.md +++ b/docs/src/man/log.md @@ -5,10 +5,6 @@ Because of the embarrassing parallelism nature of postprocessing, it is quite ea For the plotting, streamline tracing and particle tracing, a common problem is the grid and related interpolation process. Now I have [FieldTracer.jl](https://github.com/henry2004y/FieldTracer.jl) and [TestParticle.jl](https://github.com/henry2004y/TestParticle.jl) designed specifically for these tasks. -## Array Storage Ordering - -I have already made a lot of mistakes by mixing the row-major and column-major codes. Explicitly list all the parts that require extra care! - ### VTK AMR Grid Structure `vtkOverlappingAMR` implements a somewhat strict Berger-Collela AMR scheme: diff --git a/src/plot/pyplot.jl b/src/plot/pyplot.jl index 177760a8..07842498 100644 --- a/src/plot/pyplot.jl +++ b/src/plot/pyplot.jl @@ -1,6 +1,7 @@ # Plotting functionalities. using PyPlot +using Dierckx: Spline2D export plotdata, plotlogdata, animatedata, plot, scatter, contour, contourf, plot_surface, tricontourf, plot_trisurf, streamplot, streamslice, cutplot @@ -11,7 +12,7 @@ import PyPlot: plot, scatter, contour, contourf, plot_surface, tricontourf, include("utility.jl") """ - plotlogdata(data, head, func, plotmode="line") + plotlogdata(data, head, func; plotmode="line") Plot information from log file. # Input arguments @@ -142,19 +143,19 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar", if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar", "contlog","contbarlog") - xi, yi, wi = getdata(data, var, plotrange, plotinterval) + Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval) # More robust method needed! if plotmode[ivar] ∈ ["contbar", "contbarlog"] if level == 0 - c = contourf(xi, yi, wi') + c = contourf(Xi, Yi, Wi) else - c = contourf(xi, yi, wi', level) + c = contourf(Xi, Yi, Wi, level) end elseif plotmode[ivar] ∈ ["cont", "contlog"] - c = contour(xi, yi, wi) + c = contour(Xi, Yi, Wi) elseif plotmode[ivar] ∈ ["surfbar", "surfbarlog"] - c = plot_surface(xi, yi, wi) + c = plot_surface(Xi, Yi, Wi) end occursin("bar", plotmode[ivar]) && colorbar() @@ -196,12 +197,13 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar", VarIndex2_ = findindex(data, VarStream[2]) if data.head.gencoord # Generalized coordinates - X, Y = vec(x[:,:,1]), vec(x[:,:,2]) + xrange = @view x[:,:,1] + yrange = @view x[:,:,2] if any(isinf.(plotrange)) - if plotrange[1] == -Inf plotrange[1] = minimum(X) end - if plotrange[2] == Inf plotrange[2] = maximum(X) end - if plotrange[3] == -Inf plotrange[3] = minimum(Y) end - if plotrange[4] == Inf plotrange[4] = maximum(Y) end + if plotrange[1] == -Inf plotrange[1] = minimum(xrange) end + if plotrange[2] == Inf plotrange[2] = maximum(xrange) end + if plotrange[3] == -Inf plotrange[3] = minimum(yrange) end + if plotrange[4] == Inf plotrange[4] = maximum(yrange) end end # Create grid values first. @@ -210,9 +212,9 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar", # The PyCall here can be potentially replaced with Spline2D. # Perform linear interpolation of the data (x,y) on grid(xi,yi) - triang = matplotlib.tri.Triangulation(X,Y) - Xi = [y for _ in xi, y in yi] - Yi = [x for x in xi, _ in yi] + triang = matplotlib.tri.Triangulation(xrange[:], yrange[:]) + + Xi, Yi = meshgrid(xi, yi) W = w[:,1,VarIndex1_] interpolator = matplotlib.tri.LinearTriInterpolator(triang, W) @@ -223,32 +225,31 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar", v2 = interpolator(Xi, Yi) else # Cartesian coordinates - X = x[:,1,1] - Y = x[1,:,2] + xrange = @view x[:,1,1] + yrange = @view x[1,:,2] if all(isinf.(plotrange)) - Xi, Yi = X, Y + Xi, Yi = meshgrid(xrange, yrange) v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]' else - if plotrange[1] == -Inf plotrange[1] = minimum(X) end - if plotrange[2] == Inf plotrange[2] = maximum(X) end - if plotrange[3] == -Inf plotrange[3] = minimum(Y) end - if plotrange[4] == Inf plotrange[4] = maximum(Y) end + if plotrange[1] == -Inf plotrange[1] = minimum(xrange) end + if plotrange[2] == Inf plotrange[2] = maximum(xrange) end + if plotrange[3] == -Inf plotrange[3] = minimum(yrange) end + if plotrange[4] == Inf plotrange[4] = maximum(yrange) end w1, w2 = w[:,:,VarIndex1_], w[:,:,VarIndex2_] xi = range(plotrange[1], stop=plotrange[2], step=plotinterval) yi = range(plotrange[3], stop=plotrange[4], step=plotinterval) - Xi = [i for i in xi, _ in yi] - Yi = [j for _ in xi, j in yi] + Xi, Yi = meshgrid(xi, yi) - spline = Spline2D(X, Y, w1) + spline = Spline2D(xrange, yrange, w1) v1 = spline(Xi[:], Yi[:]) - v1 = reshape(v1, size(Xi))' + v1 = reshape(v1, size(Xi)) - spline = Spline2D(X, Y, w2) + spline = Spline2D(xrange, yrange, w2) v2 = spline(Xi[:], Yi[:]) - v2 = reshape(v2, size(Xi))' + v2 = reshape(v2, size(Xi)) end end @@ -360,12 +361,11 @@ function plotdata(data::Data, func::AbstractString; cut="", plotmode="contbar", #yi = range(cut2[1,1], stop=cut2[end,1], length=size(cut2)[1]) xi = range(cut1[1,1], stop=cut1[1,end], - step=(cut1[1,end]-cut1[1,1])/(size(cut1,2)-1)) + step=(cut1[1,end] - cut1[1,1]) / (size(cut1,2) - 1)) yi = range(cut2[1,1], stop=cut2[end,1], - step=(cut2[end,1]-cut2[1,1])/(size(cut2,1)-1)) + step=(cut2[end,1] - cut2[1,1]) / (size(cut2,1) - 1)) - Xi = [i for _ in yi, i in xi] - Yi = [j for j in yi, _ in xi] + Xi, Yi = meshgrid(xi, yi) s = streamplot(Xi,Yi,v1,v2; color="w", linewidth=1.0, density) end @@ -492,12 +492,11 @@ function streamslice(data::Data, var::AbstractString; end xi = range(cut1[1,1], stop=cut1[end,1], - step = (cut1[end,1]-cut1[1,1])/(size(cut1,1)-1)) + step = (cut1[end,1] - cut1[1,1]) / (size(cut1,1) - 1)) yi = range(cut2[1,1], stop=cut2[1,end], - step = (cut2[1,end]-cut2[1,1])/(size(cut2,2)-1)) + step = (cut2[1,end] - cut2[1,1]) / (size(cut2,2) - 1)) - Xi = [i for _ in yi, i in xi] - Yi = [j for j in yi, _ in xi] + Xi, Yi = meshgrid(xi, yi) s = streamplot(Xi,Yi,v1',v2'; kwargs...) @@ -548,15 +547,13 @@ Wrapper over the contour function in matplotlib. function contour(data::Data, var::AbstractString, levels=0; plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...) - xi, yi, wi = getdata(data, var, plotrange, plotinterval) + Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval) if levels != 0 - c = plt.contour(xi, yi, wi', levels; kwargs...) + c = plt.contour(Xi, Yi, Wi, levels; kwargs...) else - c = plt.contour(xi, yi, wi'; kwargs...) + c = plt.contour(Xi, Yi, Wi; kwargs...) end - - return c::PyCall.PyObject end """ @@ -567,12 +564,12 @@ Wrapper over the contourf function in matplotlib. function contourf(data::Data, var::AbstractString, levels=0; plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...) - xi, yi, wi = getdata(data, var, plotrange, plotinterval) + Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval) if levels != 0 - c = plt.contourf(xi, yi, wi', levels; kwargs...) + c = plt.contourf(Xi, Yi, Wi, levels; kwargs...) else - c = plt.contourf(xi, yi, wi'; kwargs...) + c = plt.contourf(Xi, Yi, Wi; kwargs...) end end @@ -600,7 +597,7 @@ function tricontourf(data::Data, var::AbstractString; W = W[xyIndex] end - c = tricontourf(X, Y, W) + tricontourf(X, Y, W) end """ @@ -628,7 +625,7 @@ function plot_trisurf(data::Data, var::AbstractString; W = W[xyIndex] end - c = plot_trisurf(X, Y, W) + plot_trisurf(X, Y, W) end @@ -640,12 +637,9 @@ Wrapper over the plot_surface function in matplotlib. function plot_surface(data::Data, var::AbstractString; plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs=Dict()) - xi, yi, wi = getdata(data, var, plotrange, plotinterval) + Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval) - Xi = [y for _ in xi, y in yi] - Yi = [x for x in xi, _ in yi] - - c = plot_surface(Xi, Yi, wi; kwargs...) + plot_surface(Xi, Yi, Wi; kwargs...) end """ @@ -679,8 +673,7 @@ function streamplot(data::Data, var::AbstractString; # Is there a triangulation method in Julia? tr = matplotlib.tri.Triangulation(X, Y) - Xi = [y for _ in xi, y in yi] - Yi = [x for x in xi, _ in yi] + Xi, Yi = meshgrid(xi, yi) interpolator = matplotlib.tri.LinearTriInterpolator(tr, w[:,1,VarIndex1_]) v1 = interpolator(Xi, Yi) @@ -689,34 +682,33 @@ function streamplot(data::Data, var::AbstractString; v2 = interpolator(Xi, Yi) else # Cartesian coordinates - X = x[:,1,1] - Y = x[1,:,2] + xrange = @view x[:,1,1] + yrange = @view x[1,:,2] if all(isinf.(plotrange)) - Xi, Yi = X, Y + Xi, Yi = x[:,:,1]', x[:,:,2]' v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]' else - if plotrange[1] == -Inf plotrange[1] = minimum(X) end - if plotrange[2] == Inf plotrange[2] = maximum(X) end - if plotrange[3] == -Inf plotrange[3] = minimum(Y) end - if plotrange[4] == Inf plotrange[4] = maximum(Y) end + if plotrange[1] == -Inf plotrange[1] = minimum(xrange) end + if plotrange[2] == Inf plotrange[2] = maximum(xrange) end + if plotrange[3] == -Inf plotrange[3] = minimum(yrange) end + if plotrange[4] == Inf plotrange[4] = maximum(yrange) end w1, w2 = w[:,:,VarIndex1_], w[:,:,VarIndex2_] xi = range(plotrange[1], stop=plotrange[2], step=plotinterval) yi = range(plotrange[3], stop=plotrange[4], step=plotinterval) - Xi = [i for i in xi, _ in yi] - Yi = [j for _ in xi, j in yi] + Xi, Yi = meshgrid(xi, yi) - spline = Spline2D(X, Y, w1) + spline = Spline2D(xrange, yrange, w1) v1 = spline(Xi[:], Yi[:]) - v1 = reshape(v1, size(Xi))' + v1 = reshape(v1, size(Xi)) - spline = Spline2D(X, Y, w2) + spline = Spline2D(xrange, yrange, w2) v2 = spline(Xi[:], Yi[:]) - v2 = reshape(v2, size(Xi))' + v2 = reshape(v2, size(Xi)) end end - c = streamplot(Xi, Yi, v1, v2; kwargs...) + streamplot(Xi, Yi, v1, v2; kwargs...) end \ No newline at end of file diff --git a/src/plot/utility.jl b/src/plot/utility.jl index 87d751f1..a3e50731 100644 --- a/src/plot/utility.jl +++ b/src/plot/utility.jl @@ -1,6 +1,6 @@ # Utility functions for plotting. -"Prepare data for passing to plotting functions." +"Prepare 2D data arrays for passing to plotting functions." function getdata(data::Data, var::AbstractString, plotrange, plotinterval) x, w = data.x, data.w ndim = data.head.ndim @@ -24,35 +24,32 @@ function getdata(data::Data, var::AbstractString, plotrange, plotinterval) # Perform linear interpolation of the data (x,y) on grid(xi,yi) triang = matplotlib.tri.Triangulation(X,Y) interpolator = matplotlib.tri.LinearTriInterpolator(triang, W) - Xi = [x for x in xi, _ in yi] - Yi = [y for _ in xi, y in yi] - wi = interpolator(Xi, Yi) + Xi, Yi = meshgrid(xi, yi) + Wi = interpolator(Xi, Yi) else # Cartesian coordinates + xrange = @view x[:,1,1] + yrange = @view x[1,:,2] if all(isinf.(plotrange)) - xi = x[:,1,1] - yi = x[1,:,2] - wi = w[:,:,VarIndex_] + Xi, Yi = meshgrid(xrange, yrange) + Wi = w[:,:,VarIndex_]' else - X = x[:,1,1] - Y = x[1,:,2] - if plotrange[1] == -Inf plotrange[1] = minimum(X) end - if plotrange[2] == Inf plotrange[2] = maximum(X) end - if plotrange[3] == -Inf plotrange[3] = minimum(Y) end - if plotrange[4] == Inf plotrange[4] = maximum(Y) end + if plotrange[1] == -Inf plotrange[1] = minimum(xrange) end + if plotrange[2] == Inf plotrange[2] = maximum(xrange) end + if plotrange[3] == -Inf plotrange[3] = minimum(yrange) end + if plotrange[4] == Inf plotrange[4] = maximum(yrange) end W = w[:,:,VarIndex_] xi = range(plotrange[1], stop=plotrange[2], step=plotinterval) yi = range(plotrange[3], stop=plotrange[4], step=plotinterval) - spline = Spline2D(X, Y, W) - Xi = [x for x in xi, _ in yi] - Yi = [y for _ in xi, y in yi] + spline = Spline2D(xrange, yrange, W) + Xi, Yi = meshgrid(xi, yi) wi = spline(Xi[:], Yi[:]) - wi = reshape(wi, size(Xi))' + Wi = reshape(wi, size(Xi)) end end - xi, yi, wi + Xi, Yi, Wi end @@ -62,3 +59,10 @@ function findindex(data::Data, var::AbstractString) isnothing(VarIndex_) && error("$(var) not found in file header variables!") VarIndex_ end + +"Generating consistent 2D arrays for passing to plotting functions." +function meshgrid(x, y) + X = [x for _ in y, x in x] + Y = [y for y in y, _ in x] + X, Y +end \ No newline at end of file