diff --git a/src/Diagnostic.jl b/src/Diagnostic.jl index 33a89b6..e0e67c8 100644 --- a/src/Diagnostic.jl +++ b/src/Diagnostic.jl @@ -221,105 +221,3 @@ function plt(probe::PSxProbe) heatmap(y, x, array, yflip = true, legend=false, title="Nvx", titlefont=font(10), c=:bwr, clim=(-lim, lim)) end - - - -mutable struct Diagnostic - nt::Int64 - nv::Int64 - nx::Int64 - ts::Vector{Float64} - rhos::Array{Float64, 2} - particle_rhos::Array{Float64, 3} - Us::Array{Float64, 2} - Es::Array{Float64, 2} - vs::Array{Float64, 2} - PSs::Array{Float64, 3} - xmax::Float64 - vmax::Float64 - tmax::Float64 - vbins::Vector{Float64} - N::Int64 -end - -function Diagnostic(nx::Int64, nv::Int64, nt::Int64, nspecies::Int64, xmax::Float64, vmax::Float64, tmax::Float64) - Diagnostic(nt, - nv, - nx, - zeros(nt), - zeros(nt, nx), - zeros(nt, nspecies, nx), - zeros(nt, nx), - zeros(nt, nx), - zeros(nt, nv), - zeros(nt, nx, nv), - xmax, - vmax, - tmax, - LinRange(-vmax-xmax/(nv-1),vmax+xmax/(nv-1),nv+1), - 0 - ) -end - -function histogram_v(particles::ParticleEnsemble, diag::Diagnostic, dir::Int64=1) - vs = @view diag.vs[diag.N,:] - for p in particles.coords - v = p.v[dir] - if abs(v) >= diag.vmax - continue - end - vs[Int64(floor(diag.nv * (v / (2 * diag.vmax) + 0.5))) + 1] += particles.q - end -end - -function histogram_xv(particles::ParticleEnsemble, pic::PIC, diag::Diagnostic) - PSs = @view diag.PSs[diag.N,:,:] - for p in particles.coords - v = p.v[1] - x = p.r[1] - if abs(v) >= diag.vmax - continue - end - PSs[Int64(floor(x*(pic.nx-1)/pic.xmax+0.5)+1), - Int64(floor(diag.nv * (v / (2 * diag.vmax) + 0.5))) + 1] += particles.q - end -end - -function sample(diag::Diagnostic, t::Float64, pic::PIC) - diag.N += 1 - - for p in pic.particles -# diag.vs[diag.N] += fit(Histogram, p.vx, diag.vbins).weights*p.q - histogram_v(p, diag) -# diag.PSs[diag.N] += fit(Histogram, (p.x, p.vx), (pic.xbins, diag.vbins)).weights*p.q - histogram_xv(p, pic, diag) - - end - diag.Us[diag.N,:] .= pic.U - diag.Es[diag.N,:] .= [e[1] for e in pic.E] - diag.rhos[diag.N,:] .= pic.rho - diag.ts[diag.N] = t - -end - -function plt(diag::Diagnostic, pic::PIC) - - function plt_record(x, y, array, title) - lim = maximum(abs.(array)) - return heatmap(x, y, array, yflip = true, legend=false, - title=title, titlefont=font(10), c=:bwr, clim=(-lim, lim)) - end - - p1 = plt_record(pic.x, diag.ts[1:diag.N], diag.rhos[1:diag.N,:], "Charge density") - p2 = plt_record(pic.x, diag.ts[1:diag.N], diag.Us[1:diag.N,:], "Electrical potential") - p3 = plt_record(pic.x, diag.ts[1:diag.N], diag.Es[1:diag.N,:], "Elextrical field") - vx = LinRange(-diag.vmax,diag.vmax,diag.nv) - p4 = plt_record(vx, diag.ts[1:diag.N], diag.vs[1:diag.N,:], "Speed") - - return plot(p1, p2, p3, p4, size = (900, 250), layout=(1,4)) -end - -function plt_PS(diag::Diagnostic, lims=1e-18) - x = dropdims(sum(diag.PSs, dims=1), dims=1)' - return heatmap(x, size=(500, 500), legend=false, c=:bwr, clim=(-lims, lims), title="Phase space") -end \ No newline at end of file diff --git a/src/SimplePIC.jl b/src/SimplePIC.jl index 5643067..6dcc885 100644 --- a/src/SimplePIC.jl +++ b/src/SimplePIC.jl @@ -21,7 +21,7 @@ export init_time, init_monoenergetic, init_thermal export advance, advance!, energy, scatter export epsilon_0, e -export PIC, Diagnostic, RhoProbe, EnergyProbe, NxProbe, NvxProbe, EProbe, PSxProbe, UProbe, sample! +export PIC, RhoProbe, EnergyProbe, NxProbe, NvxProbe, EProbe, PSxProbe, UProbe, sample! export sample, poisson_solve, interpolate, advance, init_leapfrog, solve_init, solve_init_fft, particle_bc export random_maxwell_v, random_maxwell_vcomponent, random_maxwell_vflux @@ -86,12 +86,6 @@ mutable struct PIC{dim, ParticleType, GeometryType, BCType, vdim} #where ParticleType <: AbstractParticle where FieldType <: AbstractField where BCType <: BoundaryCondition particles::Vector{ParticleEnsemble{ParticleType}} interactions::Vector{Interactions{ParticleEnsemble{ParticleType}}} - nx::Int64 - xmax::Float64 - dx::Float64 - x::Vector{Float64} - xbins::Vector{Float64} - dV::Float64 geo::GeometryType BC::BCType rho::Array{Float64, dim} @@ -105,12 +99,6 @@ function PIC(particles::Vector{ParticleEnsemble{ParticleType}}, interactions::Ve geo::AbstractGeometry, BC::BCType, epsilon_0::Float64, vdim::Int64) where ParticleType where BCType <: BoundaryCondition PIC{1, ParticleType, Cartesian1D, BCType, vdim}(particles, interactions, - geo.nx, - geo.xmax, - geo.dx, - LinRange(0, geo.xmax, geo.nx), - [i-0.5 for i in 0:geo.nx]*geo.xmax/(geo.nx-1), - geo.dV, geo, BC, zeros(geo.nx), @@ -239,7 +227,7 @@ end function sample(pic::PIC) pic.rho .= pic.rhobg for p in pic.particles - pic.rho .+= sample_linear(p, pic.nx, pic.dx).*(p.q/pic.dV) + pic.rho .+= sample_linear(p, pic.geo.nx, pic.geo.dx).*(p.q/pic.geo.dV) end if typeof(pic.BC) <: BCPeriodic1D @@ -265,20 +253,20 @@ end function interpolate(pic::PIC) for p in pic.particles - interpolate_linear(p, pic.dx, pic.E) + interpolate_linear(p, pic.geo.dx, pic.E) end end function particle_bc(pic::PIC) for p in pic.particles - particle_bc(p, pic.xmax, pic.BC) + particle_bc(p, pic.geo.xmax, pic.BC) end end function advance(pic::PIC, dt::Float64, B=SVector(0., 0., 0.)) for (p, inter) in zip(pic.particles, pic.interactions) advance(p, inter, dt, B) - particle_bc(p, pic.xmax, pic.BC) + particle_bc(p, pic.geo.xmax, pic.BC) end end @@ -288,7 +276,7 @@ function init_leapfrog(pic::PIC, dt::Float64) interpolate(pic) for p in pic.particles advance_v(p, -dt/2) - particle_bc(p, pic.xmax, pic.BC) + particle_bc(p, pic.geo.xmax, pic.BC) end end