diff --git a/Project.toml b/Project.toml index 37b9421..e8344b9 100644 --- a/Project.toml +++ b/Project.toml @@ -31,7 +31,6 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" -SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/test/Project.toml b/test/Project.toml index 5a66e8a..129c2a5 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] +IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/test/manual/T01_test_PIC_1D.jl b/test/manual/T01_test_PIC_1D.jl deleted file mode 100644 index c7504ad..0000000 --- a/test/manual/T01_test_PIC_1D.jl +++ /dev/null @@ -1,187 +0,0 @@ -""" -Manual diagnostic: 1D particle-in-cell test scenario. - -Interactive script for visual checks of PIC mapping behavior in 1D. -""" - -using Statistics -using Plots - -using PiCLES.ParticleMesh: OneDGrid, OneDGridNotes -import PiCLES.ParticleInCell - -import PiCLES - -function check_sum(charges_1d) - Float64(sum(charges_1d)) #/grid1d.nx -end - -using SharedArrays - -# %% define some functions -plot_path_base = "plots/tests/T01_PIC_1D/with_merge_rule/" -#mkdir(plot_path_base) - -function PIC_loop(grid1d, charges_1d, xp; N= 100, cg= 0.2, verbose=false) - x_collect = [] - p_collect = [] - push!(x_collect, xp) - push!(p_collect, charges_1d) - - #xp = grid1d.x .+ cg0 #.* (1 .+ rand()) - i_charges = charges_1d - - State = SharedMatrix{Float64}(grid1d.Nx, 1) - State .= 0 - for ti in 1:1:N - - index_positions, weights = ParticleInCell.compute_weights_and_index(grid1d, xp) - #index_positions = wrap_indexpositons(index_positions, grid1d) - - ParticleInCell.push_to_grid!(State, i_charges, index_positions, weights, grid1d.Nx, true) - - if verbose - @show ip_charges_sum - Float64(sum(State)) - end - - # update charges - i_charges = dropdims(Array{Float64}(State), dims=2)#Array{Float64}(State) #State - - # propagate - xp = grid1dnotes.x .+ cg .* i_charges #.* (0.9 .+ 0.05 *(1 .+ rand(grid1d.Nx))) - #xp = grid1d.x .+ cg0 .* (0.9 .+ 0.05 *(1 .+ rand(grid1d.Nx))) - #xp = grid1d.x .+ 0.1 * (.-0.5 .+ rand(grid1d.Nx)) - #xp = xp .+ 0.1 * (.-0.2 .+ rand(grid1d.Nx)) - - push!(x_collect, xp) - push!(p_collect, copy(State)) - State .= 0 - end - return x_collect, p_collect -end - - -function animate(x_collect, p_collect, grid1d; path, name) - - anim = @animate for i in 1:1:length(p_collect) - - scatter(x_collect[i], p_collect[i], label="new") - if i > 1 - scatter!(x_collect[i-1], p_collect[i-1], color="gray", alpha=0.4, label="old") - end - title!("sum = $( round(sum(p_collect[i]) )) ") - - ylims!(-0.5, 1.9) - xlims!(grid1d.xmin, grid1d.xmax,) - end - gif(anim, joinpath(path, name * ".gif"), fps=10) - #round(sum(p_collect[1]) ) -end - -function convert_store_to_tuple(p_collect, grid1d) - store_data = cat(p_collect..., dims=2) - store_data = permutedims(store_data, (2, 1)) - x = OneDGridNotes(grid1d).x - steps = collect(range(1, size(store_data)[1], step=1)) - return (data=store_data, x=x, steps=steps) -end - -# %% -n_particles = 101 -eta_min, eta_max =0, 20 -grid1d = OneDGrid(eta_min, eta_max, n_particles) -grid1dnotes = OneDGridNotes(grid1d) - -#grid1d = OneGridPars(grid) -#xp = rand(n_particles) * eta_max - -# make State vector: -State = SharedMatrix{Float64}(grid1d.Nx, 1) -#State = grid1dnotes.x * 0 - -# initial charge position -xp = grid1dnotes.x .+ grid1dnotes.dx *1.5 #+ rand(n_particles)/5 - - -charges_1d = xp * 0 .- 0.2 -charges_1d[40:Int(ceil(grid1d.Nx*2/3))] .=1 - -ip_charges_sum = check_sum(charges_1d) -@info "charge sum $(check_sum(charges_1d))" - -# manual step -index_positions, weights = ParticleInCell.compute_weights_and_index(grid1d, xp) - -#index_positions = wrap_indexpositons(index_positions, grid1d) # is part of push_charged_to_gridpoints! -ParticleInCell.push_to_grid!(State, charges_1d, index_positions, weights, grid1d.Nx) -@info "charge sum $(check_sum(State))" - - -plot() -scatter(xp, charges_1d, label ="advected", alpha= 1, s =3 ) -scatter!(grid1dnotes.x, State, color="black", alpha=0.5, s=1, label="new") - - - -# %% -xp = grid1dnotes.x .+ grid1dnotes.dx * 1.5 #+ rand(n_particles)/5 -charges_1d = xp * 0 .- 0.2 -charges_1d[40:Int(ceil(grid1d.Nx * 2 / 3))] .= 1 - -x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; N=50, cg = 0.2, verbose=true); - -data = convert_store_to_tuple(p_collect, grid1d) -plot() -heatmap(data.x, data.steps, data.data) -#save figure -savefig(joinpath(plot_path_base, "T01_PIC_1D_forward.png")) - -animate(x_collect, p_collect, grid1d, path = plot_path_base, name = "T01_PIC_1D_forward") - - - -# %% -xp = grid1dnotes.x .+ grid1dnotes.dx * 1.5 #+ rand(n_particles)/5 -charges_1d = xp * 0 .- 0.2 -charges_1d[40:Int(ceil(grid1d.Nx * 2 / 3))] .= 1 - -x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; cg = -0.2, N=50, verbose=true); - -data = convert_store_to_tuple(p_collect, grid1d) -plot() -heatmap(data.x, data.steps, data.data) -savefig(joinpath(plot_path_base, "T01_PIC_1D_backward.png")) - -animate(x_collect, p_collect, grid1d, path = plot_path_base, name = "T01_PIC_1D_backward") - - -# %% -xp = grid1dnotes.x .+ grid1dnotes.dx * 1.5 #+ rand(n_particles)/5 -#charges_1d = rand(n_particles)* 0 .+ 1 -charges_1d = sin.(xp) *0.2 .+0.2 - -x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; cg =-0.3, N= 50, verbose=true); -animate(x_collect, p_collect, grid1d, path = plot_path_base, name = "T01_PIC_1D_backward_sin") - -data = convert_store_to_tuple(p_collect, grid1d) -plot() -heatmap(data.x, data.steps, data.data) -savefig(joinpath(plot_path_base, "T01_PIC_1D_backward_sin.png")) -animate(x_collect, p_collect, grid1d, path=plot_path_base, name="T01_PIC_1D_backward_sin") - - - -# %% diverging sin -xp = grid1dnotes.x .+ grid1dnotes.dx * 1.5 #+ rand(n_particles)/5 -#charges_1d = rand(n_particles)* 0 .+ 1 -charges_1d = sin.(xp) *0.4 .+0.2 - -x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; cg =-0.3, N= 30, verbose=true); -animate(x_collect, p_collect, grid1d, path=plot_path_base , name="T01_PIC_1D_backward_sin_div") - -data = convert_store_to_tuple(p_collect, grid1d) -plot() -heatmap(data.x, data.steps, data.data) -savefig(joinpath(plot_path_base, "T01_PIC_1D_backward_sin_div.png")) - -animate(x_collect, p_collect, grid1d, path=plot_path_base, name="T01_PIC_1D_backward_sin_div") \ No newline at end of file diff --git a/test/manual/T03_PIC_propagation_1d.jl b/test/manual/T03_PIC_propagation_1d.jl new file mode 100644 index 0000000..771e91f --- /dev/null +++ b/test/manual/T03_PIC_propagation_1d.jl @@ -0,0 +1,285 @@ +""" +Manual diagnostic: 1D particle-in-cell test scenario. + +Interactive script for visual checks of PIC mapping behavior in 1D. +""" + +using Statistics +using Plots + +using PiCLES.Grids.CartesianGrid: CartesianGridMesh1D, gridnotes_1d +import PiCLES.ParticleInCell + +import PiCLES + +function check_sum(charges_1d) + Float64(sum(charges_1d)) #/grid1d.nx +end + +using SharedArrays + +# %% define some functions +plot_path_base = "plots/tests/T01_PIC_1D/with_merge_rule/" +#mkdir(plot_path_base) + +function grid1d_axes(grid1d) + stats = grid1d.stats + return (Nx=stats.Nx.N, xmin=stats.xmin, xmax=stats.xmax, dx=stats.dx, x=gridnotes_1d(grid1d)) +end + +function compute_weights_and_index_1d(grid1d, xp::Vector{Float64}) + axes = grid1d_axes(grid1d) + index_list, weight_list = [], [] + for xi in xp + xp_normed = (xi - axes.xmin) / axes.dx + idx, wtx = ParticleInCell.get_absolute_i_and_w(xp_normed) + push!(index_list, idx) + push!(weight_list, wtx) + end + return index_list, weight_list +end + +function PIC_loop(grid1d, charges_1d, xp; N= 100, cg= 0.2, verbose=false) + axes = grid1d_axes(grid1d) + + x_collect = [] + p_collect = [] + push!(x_collect, xp) + push!(p_collect, charges_1d) + + #xp = grid1d.x .+ cg0 #.* (1 .+ rand()) + i_charges = charges_1d + + State = SharedMatrix{Float64}(axes.Nx, 1) + State .= 0 + for ti in 1:1:N + + index_positions, weights = compute_weights_and_index_1d(grid1d, xp) + #index_positions = wrap_indexpositons(index_positions, grid1d) + + ParticleInCell.push_to_grid!(State, i_charges, index_positions, weights, axes.Nx, true) + + if verbose + @show ip_charges_sum - Float64(sum(State)) + end + + # update charges + i_charges = dropdims(Array{Float64}(State), dims=2)#Array{Float64}(State) #State + + # propagate + xp = axes.x .+ cg .* i_charges #.* (0.9 .+ 0.05 *(1 .+ rand(axes.Nx))) + #xp = grid1d.x .+ cg0 .* (0.9 .+ 0.05 *(1 .+ rand(grid1d.Nx))) + #xp = grid1d.x .+ 0.1 * (.-0.5 .+ rand(grid1d.Nx)) + #xp = xp .+ 0.1 * (.-0.2 .+ rand(grid1d.Nx)) + + push!(x_collect, xp) + push!(p_collect, copy(State)) + State .= 0 + end + return x_collect, p_collect +end + + +function animate(x_collect, p_collect, grid1d; path, name) + + anim = @animate for i in 1:1:length(p_collect) + + scatter(x_collect[i], p_collect[i], label="new") + if i > 1 + scatter!(x_collect[i-1], p_collect[i-1], color="gray", alpha=0.4, label="old") + end + title!("sum = $( round(sum(p_collect[i]) )) ") + + ylims!(-0.5, 1.9) + xlims!(grid1d.stats.xmin, grid1d.stats.xmax,) + end + gif(anim, joinpath(path, name * ".gif"), fps=10) + #round(sum(p_collect[1]) ) +end + +function convert_store_to_tuple(p_collect, grid1d) + store_data = cat(p_collect..., dims=2) + store_data = permutedims(store_data, (2, 1)) + x = gridnotes_1d(grid1d) + steps = collect(range(1, size(store_data)[1], step=1)) + return (data=store_data, x=x, steps=steps) +end + +# %% +n_particles = 101 +eta_min, eta_max =0, 20 +grid1d = CartesianGridMesh1D(eta_min, eta_max, n_particles; Ny=3, periodic_boundary=true) +grid1dnotes = gridnotes_1d(grid1d) +grid1daxes = grid1d_axes(grid1d) + +#grid1d = OneGridPars(grid) +#xp = rand(n_particles) * eta_max + +# make State vector: +State = SharedMatrix{Float64}(grid1daxes.Nx, 1) +#State = grid1dnotes.x * 0 + +# initial charge position +xp = grid1dnotes .+ grid1daxes.dx *1.5 #+ rand(n_particles)/5 + + +charges_1d = xp * 0 .- 0.2 +charges_1d[40:Int(ceil(grid1daxes.Nx*2/3))] .=1 + +ip_charges_sum = check_sum(charges_1d) +@info "charge sum $(check_sum(charges_1d))" + +# manual step +index_positions, weights = compute_weights_and_index_1d(grid1d, xp) + +#index_positions = wrap_indexpositons(index_positions, grid1d) # is part of push_charged_to_gridpoints! +ParticleInCell.push_to_grid!(State, charges_1d, index_positions, weights, grid1daxes.Nx) +@info "charge sum $(check_sum(State))" + + +plot() +scatter(xp, charges_1d, label ="advected", alpha= 1, s =3 ) +scatter!(grid1dnotes, State, color="black", alpha=0.5, s=1, label="new") + + +# %% Box moving to the right, zero base +xp = grid1dnotes .+ grid1daxes.dx * 1.5 #+ rand(n_particles)/5 +charges_1d = xp * 0 +charges_1d[40:Int(ceil(grid1daxes.Nx * 2 / 3))] .= 1 + +case_title = "Box right, zero base" +case_name = "box_r0" + +x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; N=50, cg=0.2, verbose=true); + +data = convert_store_to_tuple(p_collect, grid1d) +plot() +heatmap(data.x, data.steps, data.data) +title!(case_title) +#save figure +savefig(joinpath(plot_path_base, case_name * ".png")) + +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + + + +# %% Box moving to the left, non-zero base, i.e. convergence and divergence +xp = grid1dnotes .+ grid1daxes.dx * 1.5 #+ rand(n_particles)/5 +charges_1d = xp * 0 .- 0.5 +charges_1d[40:Int(ceil(grid1daxes.Nx * 2 / 3))] .= 1 + +case_title = "Box right, non-zero base" +case_name = "box_rb" + +x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; N=100, cg = 0.2, verbose=true); + +data = convert_store_to_tuple(p_collect, grid1d) +plot() +heatmap(data.x, data.steps, data.data) +title!(case_title) +#save figure +savefig(joinpath(plot_path_base, case_name * ".png")) + +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + + +# %% box to the left, zero base +xp = grid1dnotes .+ grid1daxes.dx * 1.5 #+ rand(n_particles)/5 +charges_1d = xp * 0 +charges_1d[40:Int(ceil(grid1daxes.Nx * 2 / 3))] .= 1 + +case_title = "Box left, zero base" +case_name = "box_l0" + +x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; cg=-0.2, N=100, verbose=true); + +data = convert_store_to_tuple(p_collect, grid1d) +plot() +heatmap(data.x, data.steps, data.data) +title!(case_title) +savefig(joinpath(plot_path_base, case_name * ".png")) + +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + + + +# %% box to the left, non-zero base, i.e. convergence and divergence +xp = grid1dnotes .+ grid1daxes.dx * 1.5 #+ rand(n_particles)/5 +charges_1d = xp * 0 .- 0.5 +charges_1d[40:Int(ceil(grid1daxes.Nx * 2 / 3))] .= 1 + +case_title = "Box left, non-zero base" +case_name = "box_lb" + +x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; cg = -0.2, N=100, verbose=true); + +data = convert_store_to_tuple(p_collect, grid1d) +plot() +heatmap(data.x, data.steps, data.data) +title!(case_title) +savefig(joinpath(plot_path_base, case_name * ".png")) + +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + + + +# %% sin wave moving to the left, Speed dependent on charge, i.e. non-linear advection + +xp = grid1dnotes .+ grid1daxes.dx * 0.5 # #+ rand(n_particles)/5 + +charges_1d = sin.(3 * 2 * pi * xp ./ grid1daxes.xmax) * 0.2 .+ 0.2 + +case_title = "Sin left, nonlinear advection" +case_name = "sin_l" + +x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; cg=-0.3, N=50, verbose=true); +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + +data = convert_store_to_tuple(p_collect, grid1d) +plot() +heatmap(data.x, data.steps, data.data) +title!(case_title) +savefig(joinpath(plot_path_base, case_name * ".png")) +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + + + +# %% sin wave moving to the left, Speed dependent on charge, i.e. non-linear advection + +xp = grid1dnotes .+ grid1daxes.dx * 0.5 # #+ rand(n_particles)/5 + +charges_1d = sin.(3* 2 *pi * xp ./ grid1daxes.xmax) * 0.2 .+ 0.2 + +case_title = "Sin left, nonlinear advection" +case_name = "sin_l" + +x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; cg =-0.3, N= 50, verbose=true); +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + +data = convert_store_to_tuple(p_collect, grid1d) +plot() +heatmap(data.x, data.steps, data.data) +title!(case_title) +savefig(joinpath(plot_path_base, case_name * ".png")) +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + + + +# %% diverging sin, moving to left and right until cancelation +xp = grid1dnotes .+ grid1daxes.dx * 1.5 #+ rand(n_particles)/5 +#charges_1d = rand(n_particles)* 0 .+ 1 +charges_1d = sin.(3* 2 *pi * xp ./ grid1daxes.xmax) *0.4 .+0.0 + +case_title = "Diverging sin, canceling" +case_name = "sin_div" + +x_collect, p_collect = PIC_loop(grid1d, charges_1d, xp; cg =-0.3, N= 30, verbose=true); +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) + +data = convert_store_to_tuple(p_collect, grid1d) +plot() +heatmap(data.x, data.steps, data.data) +title!(case_title) +savefig(joinpath(plot_path_base, case_name * ".png")) + +animate(x_collect, p_collect, grid1d, path=plot_path_base, name=case_name) \ No newline at end of file diff --git a/test/manual/T03_PIC_propagation_2d_blob.jl b/test/manual/T03_PIC_propagation_2d_blob.jl index d49fd83..6acab08 100644 --- a/test/manual/T03_PIC_propagation_2d_blob.jl +++ b/test/manual/T03_PIC_propagation_2d_blob.jl @@ -4,16 +4,14 @@ using Pkg Pkg.activate("PiCLES/") """ -Manual diagnostic: 2D PIC propagation in a blob scenario. - -Interactive script for exploratory stepping and visual inspection. +Manual diagnostic: 2D PIC propagation from a seeded blob. This script runs a compact case matrix (upper-right, right-only, bottom-only, lower-left propagation; each with dissipation on/off) on a periodic Cartesian grid using simulation-level forcing, advances each case with repeated remeshing, and produces a 3x2 diagnostic figure (binary on/boundary masks, state heatmaps, and total-energy time series) with descriptive metadata in the suptitle, saving one final plot per case to `save_path`. """ import Plots as plt using Setfield, IfElse -using PiCLES.ParticleSystems: particle_waves_v6 as PW +using PiCLES.ParticleSystems import PiCLES: FetchRelations, ParticleTools using PiCLES.Operators.core_2D: ParticleDefaults, InitParticleInstance, GetGroupVelocity @@ -47,6 +45,7 @@ using StructArrays using BenchmarkTools using PiCLES.Grids +using PiCLES using PiCLES.Operators.TimeSteppers: time_step! #using OrdinaryDiffEq @@ -61,40 +60,26 @@ pwd() # % Parameters U10, V10 = 0.00, 00.0 DT = 20minutes -ODEpars, Const_ID, Const_Scg = PW.ODEParameters(r_g=0.85) +ODEpars, Const_ID, Const_Scg = ODEParameters(r_g=0.85) u(x, y, t) = IfElse.ifelse.(x .< 250e3, U10, 0.00) + y * 0.0 + t * 0.0 v(x, y, t) = IfElse.ifelse.(x .< 250e3, V10, 0.00) + y * 0.0 + t * 0.0 .* cos(t * 5 / (1 * 60 * 60 * 2π)) # u(x, y, t) = U10 + x * 0.0 + y * 0.0 + t * 0.0 # v(x, y, t) = V10 + x * 0.0 + y * 0.0 + t * 0.0 -winds = (u=u, v=v) +forcing = PiCLES.custom_structures.ForcingCollection(u_wind=u, v_wind=v) -grid = TwoDCartesianGridMesh(400e3, 41, 200e3, 21; periodic_boundary=(false, true)) -heatmap(grid.data.x[:, 1], grid.data.y[1, :], transpose(grid.data.mask)) +grid = TwoDCartesianGridMesh(400e3, 41, 300e3, 31; periodic_boundary=(true, true)) # %% -particle_system = PW.particle_equations(u, v, γ=Const_ID.γ, q=Const_ID.q, - propagation=true, - input=false, - dissipation=true, - peak_shift=false, - direction=false, - ); - -default_ODE_parameters = (r_g=ODEpars.r_g, C_α=Const_Scg.C_alpha, - C_φ=Const_ID.c_β, C_e=Const_ID.C_e, g=9.81);#, M=M); - # define setting and standard initial conditions WindSeamin = FetchRelations.MinimalWindsea(U10, V10, DT); lne_local = round( log(WindSeamin["E"]) , digits=4) -cg_u_local = round( WindSeamin["cg_bar_x"] , digits=4) -cg_v_local = round(WindSeamin["cg_bar_y"], digits=4) ParticleMin = FetchRelations.MinimalParticle(2, 0, DT) # get_initial_windsea(2, 2, DT,particle_state=true ) -ODE_settings = PW.ODESettings( +ODE_settings = ODESettings( Parameters=ODEpars, # define mininum energy threshold log_energy_minimum=lne_local,#log(FetchRelations.Eⱼ(0.1, DT)), @@ -111,93 +96,123 @@ ODE_settings = PW.ODESettings( # %% -function plot_particle_collection(wave_model) +function plot_state(wave_model, energy_time, energy_series; case_title="") particles = wave_model.ParticleCollection - p = plot(layout=(3, 2), size=(1200, 1100)) - heatmap!(p, transpose(particles.on), subplot=1, title="on | iter=" * string(wave_model.clock.iteration) * " | total energy = " * string(round(sum(wave_model.State[:, :, 1]), digits=4))) - heatmap!(p, transpose(particles.boundary), subplot=2, title="boundary") + total_energy = sum(wave_model.State[:, :, 1]) + meta_title = "iter=$(wave_model.clock.iteration) | total energy=$(round(total_energy, digits=4))" + suptitle = isempty(case_title) ? meta_title : "$(case_title) | $(meta_title)" + ymax = isempty(energy_series) ? 1.0 : maximum(energy_series) + ymax_plot = max(1e-8, 1.1 * ymax) + ymin_plot = isempty(energy_series) ? 0.0 : 0.95 * minimum(energy_series) + + p = plot(layout=(3, 2), size=(1000, 1200), plot_title=suptitle) + heatmap!(p, transpose(particles.on), subplot=1, + title="on (binary: 1=active, 0=inactive)", + clims=(0, 1), c=:grays, colorbar=true) + heatmap!(p, transpose(particles.boundary), subplot=2, + title="boundary (binary: 1=boundary, 0=interior)", + clims=(0, 1), c=:grays, colorbar=true) heatmap!(p, transpose(wave_model.State[:, :, 1]), subplot=3, title="State: Energy", clims=(0, NaN)) heatmap!(p, transpose(wave_model.State[:, :, 2]), subplot=4, title="State: x momentum ", clims=(0, NaN)) + plot!(p, energy_time, energy_series, subplot=5, + title="Total energy over time", + xlabel="time [h]", ylabel="sum(E)", + ylims=(ymin_plot, ymax_plot), + label="total E", color=:black) heatmap!(p, transpose(wave_model.State[:, :, 3]), subplot=6, title="State: y momentum ") - # title = plot!(title="Plot title", grid=false, showaxis=false, bottom_margin=-50Plots.px) - plot!(p, aspect_ratio=:equal) + # Keep map-style panels square; let the time-series panel use a free aspect ratio. + for s in (1, 2, 3, 4, 6) + plot!(p, subplot=s, aspect_ratio=:equal) + end + plot!(p, subplot=5, aspect_ratio=:none) display(p) + return p end # %% -default_windsea = FetchRelations.get_initial_windsea(2, 2, DT,particle_state=true ) -# default_particle = ParticleDefaults(lne_local, cg_u_local, cg_v_local, 0.0, 0.0) -wave_model = WaveGrowthModels2D.WaveGrowth2D(; - grid = grid, - winds = winds, - ODEsys = particle_system, - ODEsets = ODE_settings, # ODE_settings - #ODEinit_type=ParticleDefaults2D(default_windsea[1]*0.1, default_windsea[2]*0.1, default_windsea[3]*0.1, 0.0, 0.0), - ODEinit_type=ParticleDefaults(ParticleMin), - #ParticleDefaults2D(log(2), 0.0, 0.0, 0.0, 0.0), #"wind_sea", # default_ODE_valuves - periodic_boundary= true, - boundary_type= "same", - #minimal_particle=FetchRelations.MinimalParticle(U10, V10, DT), - movie = true) - -#wave_model.ODEsettings.wind_min_squared = 0.0 -#wave_model.minimal_state = 2 * wave_model.minimal_state - using PiCLES.Operators.mapping_2D: reset_PI_u!, ParticleToNode! -# ### build Simulation -wave_simulation = Simulation(wave_model, Δt=20minutes, stop_time=4hours)#1hours) -initialize_simulation!(wave_simulation) -# plot_particle_collection(wave_model) - - -for PI in wave_simulation.model.ParticleCollection[5:15, 5:15] - #reset_PI_u!(PI, ui= [ log( (5/4)^2) , 5.0 , 10.0, 0.0, 0.0]) - reset_PI_u!(PI, ui= FetchRelations.get_initial_windsea(10.0, 12.0, 4hour, particle_state=true) ) - ParticleToNode!(PI, wave_simulation.model.State, wave_simulation.model.grid, wave_simulation.model.periodic_boundary) +function seed_blob!(wave_simulation; cgx, cgy, patch_i=5:15, patch_j=5:15) + for PI in wave_simulation.model.ParticleCollection[patch_i, patch_j] + reset_PI_u!(PI, ui=FetchRelations.get_initial_windsea(cgx, cgy, 2hour, particle_state=true)) + ParticleToNode!(PI, wave_simulation.model.State, wave_simulation.model.grid, wave_simulation.model.periodic_boundary) + end end -plot_particle_collection(wave_model) - - -# %% Advance only -wave_simulation.model.ParticleCollection[8, 8].ODEIntegrator.u -wave_simulation.model.ParticleCollection[5, 2].ODEIntegrator.u - -FailedCollection = Vector{AbstractMarkedParticleInstance}([]) -TimeSteppers.time_step!_advance(wave_simulation.model, wave_simulation.Δt, FailedCollection) - -wave_simulation.model.ParticleCollection[20, 12].ODEIntegrator.u -wave_simulation.model.ParticleCollection[8, 8].ODEIntegrator.u -wave_simulation.model.ParticleCollection[5, 2].ODEIntegrator.u - - -TimeSteppers.time_step!_remesh(wave_simulation.model, wave_simulation.Δt) - -wave_simulation.model.ParticleCollection[20, 12].ODEIntegrator.u -wave_simulation.model.ParticleCollection[8, 8].ODEIntegrator.u -wave_simulation.model.ParticleCollection[5, 2].ODEIntegrator.u - -wave_simulation.model.State[5, 2, :] -wave_simulation.model.State[:,:,1] - -#TimeSteppers.time_step!(wave_simulation.model, wave_simulation.Δt) - - -# %% timstepper test -plot_particle_collection(wave_model) - -for i in 1:1:100 - TimeSteppers.time_step!(wave_simulation.model, wave_simulation.Δt) - - if i%8 == 0 - plot_particle_collection(wave_simulation.model) - sleep(0.02) +function run_case(; case_title, case_slug, cgx, cgy, dissipation) + particle_system = particle_equations( + γ=Const_ID.γ, + q=Const_ID.q, + propagation=true, + input=false, + dissipation=dissipation, + peak_shift=false, + direction=false, + ) + + wave_model = WaveGrowthModels2D.WaveGrowth2D( + grid=grid, + ODEsys=particle_system, + ODEsets=ODE_settings, + ODEinit_type=ParticleDefaults(ParticleMin), + periodic_boundary=true, + boundary_type="same", + movie=true, + ) + + wave_simulation = Simulation( + wave_model, + forcing=forcing, + Δt=20minutes, + stop_time=16hours, + ) + initialize_simulation!(wave_simulation) + seed_blob!(wave_simulation; cgx=cgx, cgy=cgy) + + energy_time = [wave_simulation.model.clock.time / 3600] + energy_series = [sum(wave_simulation.model.State[:, :, 1])] + + for i in 1:600 + TimeSteppers.time_step!(wave_simulation.model, wave_simulation.Δt; forcing=wave_simulation.forcing) + + push!(energy_time, wave_simulation.model.clock.time / 3600) + push!(energy_series, sum(wave_simulation.model.State[:, :, 1])) + + if i % 8 == 0 + plot_state(wave_simulation.model, energy_time[2:end], energy_series[2:end]; case_title=case_title) + sleep(0.02) + end + + wave_simulation.model.State[:, :, :] .= 0.0 end - wave_simulation.model.State[:, :, :] .= 0.0 + p_final = plot_state(wave_simulation.model, energy_time[2:end], energy_series[2:end]; case_title=case_title) + savefig(p_final, joinpath(save_path, "T03_blob_" * case_slug * ".png")) +end + +# %% Case matrix: 4 propagation directions x dissipation on/off +cases = [ + (title="Upper-right propagation | dissipation ON", slug="ur_diss_on", cgx=10.0, cgy=12.0, dissipation=true), + (title="Upper-right propagation | dissipation OFF", slug="ur_diss_off", cgx=10.0, cgy=12.0, dissipation=false), + (title="Right-only propagation | dissipation ON", slug="r_diss_on", cgx=12.0, cgy=0.0, dissipation=true), + (title="Right-only propagation | dissipation OFF", slug="r_diss_off", cgx=12.0, cgy=0.0, dissipation=false), + (title="Bottom-only propagation | dissipation ON", slug="b_diss_on", cgx=0.0, cgy=-12.0, dissipation=true), + (title="Bottom-only propagation | dissipation OFF", slug="b_diss_off", cgx=0.0, cgy=-12.0, dissipation=false), + (title="Lower-left propagation | dissipation ON", slug="ll_diss_on", cgx=-10.0, cgy=-12.0, dissipation=true), + (title="Lower-left propagation | dissipation OFF", slug="ll_diss_off", cgx=-10.0, cgy=-12.0, dissipation=false), +] + +for case in cases + @info "Running case: $(case.title)" + run_case( + case_title=case.title, + case_slug=case.slug, + cgx=case.cgx, + cgy=case.cgy, + dissipation=case.dissipation, + ) end diff --git a/test/manual/T03_PIC_propagation_2d_land.jl b/test/manual/T03_PIC_propagation_2d_land.jl index 3843d70..3e39cfd 100644 --- a/test/manual/T03_PIC_propagation_2d_land.jl +++ b/test/manual/T03_PIC_propagation_2d_land.jl @@ -12,7 +12,7 @@ Interactive script for boundary and mask behavior inspection. import Plots as plt using Setfield, IfElse -using PiCLES.ParticleSystems: particle_waves_v6 as PW +using PiCLES.ParticleSystems import PiCLES: FetchRelations, ParticleTools using PiCLES.Operators.core_2D: ParticleDefaults, InitParticleInstance, GetGroupVelocity @@ -49,6 +49,7 @@ using StructArrays using BenchmarkTools using PiCLES.Grids +using PiCLES using PiCLES.Operators.TimeSteppers: time_step! #using OrdinaryDiffEq @@ -63,7 +64,7 @@ pwd() U10, V10 = 15.0, -10.0 DT = 20minutes -ODEpars, Const_ID, Const_Scg = PW.ODEParameters(r_g=0.85) +ODEpars, Const_ID, Const_Scg = ODEParameters(r_g=0.85) u(x, y, t) = IfElse.ifelse.(x .< 250e3, U10, 0.00) + y * 0.0 + t * 0.0 @@ -71,7 +72,7 @@ v(x, y, t) = IfElse.ifelse.(x .< 250e3, V10, 0.00) + y * 0.0 + t * 0.0 .* cos(t # u(x, y, t) = U10 + x * 0.0 + y * 0.0 + t * 0.0 # v(x, y, t) = V10 + x * 0.0 + y * 0.0 + t * 0.0 -winds = (u=u, v=v) +forcing = PiCLES.custom_structures.ForcingCollection(u_wind=u, v_wind=v) grid = TwoDCartesianGridMesh(400e3, 41, 200e3, 21, periodic_boundary=(false, true)) heatmap(grid.data.x[:, 1], grid.data.y[1, :], transpose(grid.data.mask)) @@ -89,7 +90,7 @@ heatmap(grid.data.x[:,1], grid.data.y[1,:], transpose(grid.data.mask)) # heatmap(transpose(v.(grid.data.x, grid.data.y, 0))) # %% -particle_system = PW.particle_equations(u, v, γ=Const_ID.γ, q=Const_ID.q, +particle_system = particle_equations(γ=Const_ID.γ, q=Const_ID.q, propagation=true, input=false, dissipation=false, @@ -106,7 +107,7 @@ lne_local = log(WindSeamin["E"]) cg_u_local = WindSeamin["cg_bar_x"] cg_v_local = WindSeamin["cg_bar_y"] -ODE_settings = PW.ODESettings( +ODE_settings = ODESettings( Parameters=ODEpars, # define mininum energy threshold log_energy_minimum=lne_local,#log(FetchRelations.Eⱼ(0.1, DT)), @@ -137,7 +138,6 @@ end default_windsea = FetchRelations.get_initial_windsea(U10, V10, DT,particle_state=true ) default_particle = ParticleDefaults(lne_local, cg_u_local, cg_v_local, 0.0, 0.0) wave_model = WaveGrowthModels2D.WaveGrowth2D(; grid=grid, - winds=winds, ODEsys=particle_system, ODEsets=ODE_settings, # ODE_settings ODEinit_type=ParticleDefaults2D(default_windsea[1], default_windsea[2], default_windsea[3], 0.0, 0.0), @@ -151,7 +151,10 @@ wave_model = WaveGrowthModels2D.WaveGrowth2D(; grid=grid, #wave_model.minimal_state = 2 * wave_model.minimal_state # ### build Simulation -wave_simulation = Simulation(wave_model, Δt=10minutes, stop_time=4hours)#1hours) +wave_simulation = Simulation(wave_model, + forcing=forcing, + Δt=10minutes, + stop_time=4hours) #1hours) initialize_simulation!(wave_simulation) plot_particle_collection(wave_model) @@ -160,8 +163,14 @@ PI = wave_simulation.model.ParticleCollection[2,1] step!(PI.ODEIntegrator, DT, true) model = wave_simulation.model +FailedCollection = Vector{AbstractMarkedParticleInstance}([]) +forcing_xy = TimeSteppers.time_slice_forcing(wave_simulation.forcing, model.grid, model.clock.time) +wind_i = ( + forcing_xy.u(PI.position_xy[1], PI.position_xy[2]), + forcing_xy.v(PI.position_xy[1], PI.position_xy[2]), +) mapping_2D.advance!(PI, model.State, FailedCollection, - model.grid, model.winds, wave_simulation.Δt, + model.grid, wind_i, wave_simulation.Δt, model.ODEsettings.log_energy_maximum, model.ODEsettings.wind_min_squared, model.periodic_boundary, @@ -171,10 +180,15 @@ mapping_2D.advance!(PI, model.State, FailedCollection, # %% Adance only test model = wave_simulation.model FailedCollection = Vector{AbstractMarkedParticleInstance}([]) +forcing_xy = TimeSteppers.time_slice_forcing(wave_simulation.forcing, model.grid, model.clock.time) for a_particle in model.ParticleCollection[findall(model.grid.data.mask .== 1)] @info a_particle.position_ij + wind_i = ( + forcing_xy.u(a_particle.position_xy[1], a_particle.position_xy[2]), + forcing_xy.v(a_particle.position_xy[1], a_particle.position_xy[2]), + ) mapping_2D.advance!(a_particle, model.State, FailedCollection, - model.grid, model.winds, wave_simulation.Δt, + model.grid, wind_i, wave_simulation.Δt, model.ODEsettings.log_energy_maximum, model.ODEsettings.wind_min_squared, model.periodic_boundary, @@ -193,7 +207,7 @@ end plot_particle_collection(wave_model) for i in 1:1:180 - TimeSteppers.time_step!(wave_simulation.model, wave_simulation.Δt) + TimeSteppers.time_step!(wave_simulation.model, wave_simulation.Δt; forcing=wave_simulation.forcing) if i%8 == 0 plot_particle_collection(wave_simulation.model) diff --git a/test/runtests.jl b/test/runtests.jl index beb0f29..9d4161c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,8 @@ using Test include("unit/test_single_particle_2d_local_winds_smoke.jl") include("unit/test_single_particle_1d_alias_smoke.jl") include("unit/test_grids_homogeneous_forcing_smoke.jl") + include("unit/test_pic_1d_charge_conservation.jl") + include("unit/test_pic_2d_propagation_energy_conservation.jl") # include("unit/test_1d_wave_growth_quantitative_periodic.jl") include("unit/test_1d_wave_growth_quantitative_single_particle.jl") end diff --git a/test/unit/test_pic_1d_charge_conservation.jl b/test/unit/test_pic_1d_charge_conservation.jl new file mode 100644 index 0000000..e035418 --- /dev/null +++ b/test/unit/test_pic_1d_charge_conservation.jl @@ -0,0 +1,91 @@ +""" +Quantitative 1D PIC conservation checks. + +This test mirrors the non-divergent manual cases from +test/manual/T01_test_PIC_1D.jl and verifies total charge remains within 2% +of the initial value over an integration that is 4x longer than in the manual +diagnostic script. +""" + +using Test +using SharedArrays + +using PiCLES.Grids.CartesianGrid: CartesianGridMesh1D, gridnotes_1d +import PiCLES.ParticleInCell + +function grid1d_axes(grid1d) + stats = grid1d.stats + return (Nx=stats.Nx.N, xmin=stats.xmin, xmax=stats.xmax, dx=stats.dx, x=gridnotes_1d(grid1d)) +end + +function compute_weights_and_index_1d(grid1d, xp::Vector{Float64}) + axes = grid1d_axes(grid1d) + index_list, weight_list = [], [] + for xi in xp + xp_normed = (xi - axes.xmin) / axes.dx + idx, wtx = ParticleInCell.get_absolute_i_and_w(xp_normed) + push!(index_list, idx) + push!(weight_list, wtx) + end + return index_list, weight_list +end + +function integrate_charge_history(grid1d, charges_1d, xp; N::Int, cg::Float64) + axes = grid1d_axes(grid1d) + q = copy(charges_1d) + x = copy(xp) + + state = SharedMatrix{Float64}(axes.Nx, 1) + charge_hist = Vector{Float64}(undef, N + 1) + charge_hist[1] = Float64(sum(q)) + + for ti in 1:N + state .= 0 + index_positions, weights = compute_weights_and_index_1d(grid1d, x) + ParticleInCell.push_to_grid!(state, q, index_positions, weights, axes.Nx, true) + + q = dropdims(Array{Float64}(state), dims=2) + x = axes.x .+ cg .* q + + charge_hist[ti + 1] = Float64(sum(q)) + end + + return charge_hist +end + +function init_box_case(axes; base::Float64) + xp = axes.x .+ axes.dx * 1.5 + charges_1d = fill(base, axes.Nx) + charges_1d[40:Int(ceil(axes.Nx * 2 / 3))] .= 1.0 + return xp, charges_1d +end + +function init_sin_left_case(axes) + xp = axes.x .+ axes.dx * 0.5 + charges_1d = sin.(3 * 2 * pi * xp ./ axes.xmax) * 0.2 .+ 0.2 + return xp, charges_1d +end + +function max_relative_charge_drift(charge_hist::Vector{Float64}) + q0 = charge_hist[1] + return maximum(abs.(charge_hist .- q0)) / abs(q0) +end + +@testset "PIC 1D total charge conservation (non-divergent cases)" begin + grid1d = CartesianGridMesh1D(0.0, 20.0, 101; Ny=3, periodic_boundary=true) + axes = grid1d_axes(grid1d) + + # N is 4x the manual script values: 50 -> 200, 100 -> 400 + cases = [ + (name="box right, zero base", N=200, cg=0.2, init=axes -> init_box_case(axes; base=0.0)), + (name="box left, zero base", N=200, cg=-0.2, init=axes -> init_box_case(axes; base=0.0)), + (name="sin left, nonlinear advection", N=200, cg=-0.3, init=init_sin_left_case), + ] + + for case in cases + xp, charges_1d = case.init(axes) + charge_hist = integrate_charge_history(grid1d, charges_1d, xp; N=case.N, cg=case.cg) + drift = max_relative_charge_drift(charge_hist) + @test drift <= 0.02 + end +end diff --git a/test/unit/test_pic_2d_propagation_energy_conservation.jl b/test/unit/test_pic_2d_propagation_energy_conservation.jl new file mode 100644 index 0000000..7cc0b28 --- /dev/null +++ b/test/unit/test_pic_2d_propagation_energy_conservation.jl @@ -0,0 +1,109 @@ +""" +2D blob energy-window test (dissipation off). + +Builds a compact version of the manual 2D blob setup and checks that, after 100 +steps, the final total domain energy remains within +/-2% of the energy at step 3. +""" + +using Test + +using PiCLES +using PiCLES.ParticleSystems +import PiCLES: FetchRelations +using PiCLES.Grids.CartesianGrid: TwoDCartesianGridMesh +using PiCLES.Models.WaveGrowthModels2D +using PiCLES.Simulations +using PiCLES.Operators: TimeSteppers +using PiCLES.Operators.core_2D: ParticleDefaults +using PiCLES.Operators.mapping_2D: reset_PI_u!, ParticleToNode! +using Oceananigans.Units + +function seed_blob!(sim; cgx, cgy, patch_i=5:15, patch_j=5:15) + for PI in sim.model.ParticleCollection[patch_i, patch_j] + reset_PI_u!(PI, ui=FetchRelations.get_initial_windsea(cgx, cgy, 2hour, particle_state=true)) + ParticleToNode!(PI, sim.model.State, sim.model.grid, sim.model.periodic_boundary) + end +end + +function energy_series_for_case(; cgx, cgy, steps=100) + U10, V10 = 0.0, 0.0 + DT = Float64(20minutes) + + ODEpars, Const_ID, _ = ODEParameters(r_g=0.85) + u(x, y, t) = ifelse.(x .< 250e3, U10, 0.00) + y * 0.0 + t * 0.0 + v(x, y, t) = ifelse.(x .< 250e3, V10, 0.00) + y * 0.0 + t * 0.0 + forcing = PiCLES.custom_structures.ForcingCollection(u_wind=u, v_wind=v) + + grid = TwoDCartesianGridMesh(400e3, 41, 300e3, 31; periodic_boundary=(true, true)) + + particle_system = particle_equations( + γ=Const_ID.γ, + q=Const_ID.q, + propagation=true, + input=false, + dissipation=false, + peak_shift=false, + direction=false, + ) + + WindSeamin = FetchRelations.MinimalWindsea(U10, V10, DT) + lne_local = log(WindSeamin["E"]) + ODE_settings = ODESettings( + Parameters=ODEpars, + log_energy_minimum=lne_local, + log_energy_maximum=log(27), + saving_step=6000, + timestep=DT, + total_time=12days, + dt=1e-3, + dtmin=1e-4, + dtmax=DT, + ) + + particle_min = FetchRelations.MinimalParticle(2, 0, DT) + + model = WaveGrowthModels2D.WaveGrowth2D( + grid=grid, + ODEsys=particle_system, + ODEsets=ODE_settings, + ODEinit_type=ParticleDefaults(particle_min), + periodic_boundary=true, + boundary_type="same", + movie=false, + ) + + sim = Simulation(model; forcing=forcing, Δt=DT, stop_time=16hours) + initialize_simulation!(sim) + seed_blob!(sim; cgx=cgx, cgy=cgy) + + energies = Float64[] + for _ in 1:steps + TimeSteppers.time_step!(sim.model, sim.Δt; forcing=sim.forcing) + push!(energies, sum(sim.model.State[:, :, 1])) + sim.model.State[:, :, :] .= 0.0 + end + return energies +end + +@testset "2D blob total energy window (dissipation off)" begin + cases = [ + (name="upper-right", cgx=10.0, cgy=12.0), + (name="right-only", cgx=12.0, cgy=0.0), + (name="bottom-only", cgx=0.0, cgy=-12.0), + (name="lower-left", cgx=-10.0, cgy=-12.0), + ] + + for case in cases + e = energy_series_for_case(; cgx=case.cgx, cgy=case.cgy, steps=100) + @test length(e) == 100 + + e_ref = e[3] + e_final = e[end] + + @test isfinite(e_ref) && isfinite(e_final) + @test e_ref != 0.0 + + rel_change = abs(e_final - e_ref) / abs(e_ref) + @test rel_change <= 0.02 + end +end