Skip to content

Stochastic PiCLES pull request#42

Open
tomprotin wants to merge 41 commits intomochell:stochastic_PiCLESfrom
tomprotin:main-fork-pull-request
Open

Stochastic PiCLES pull request#42
tomprotin wants to merge 41 commits intomochell:stochastic_PiCLESfrom
tomprotin:main-fork-pull-request

Conversation

@tomprotin
Copy link
Copy Markdown
Collaborator

This works in PiCLES 1.10.0. You can try to run tests/test_case_1.jl (but I advise decreasing the number of particles in this test case, at line 114, you can set it to 150 or 1500 to see if it works).
I don't know how your test cases work, but I changed the behavior of some critical functions like :

  • run()
  • advance() (when the actual stepping is computed for an individual particle, at lines 173 -175 in mapping_2D.jl)
  • time_step!()

…each node when using the "nonparametric" spreading type
…uctures for the validation of the 1/r decay of energy density
@mochell mochell self-assigned this Nov 5, 2025
Copy link
Copy Markdown
Owner

@mochell mochell left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • move SIR_analysis maybe into analysis/SIR_analysis
  • remove the listed files in .gignore from the branch and then remove them from .gitignore
  • adjust Abstract model defition and the timestepping methods (as discussed)
  • time_step!(model::Abstract2DModel, ... ):
    • Look at the version of the function in this branch (my version), I refactors it with time_step!_advance, you could do the same here.
  • core_2D_spread.jl It maybe the easierst to rename it back to core_2D.jl to see the differences better. I have not looked at this file in more detail
  • Mapping_2d.jl the ParticleToNode! routines changed. I deleted the definitions that are not used any more.. I hope it still works..
  • Do you actually use particle_waves_v6.jl ? I think you can just use v5 and and then we can get rid of v6

nPart,
wnds,
cur,
Mstat} <: Abstract2DModel where {Mstat<:Union{Nothing,stat}, PCollection<:Union{Vector,Array}}
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change Abstract2DModel to something like Abstract2DStochasticModel. Define the model Type in Architectures.jl, and then make a timestepper for that model type

Suggested change
Mstat} <: Abstract2DModel where {Mstat<:Union{Nothing,stat}, PCollection<:Union{Vector,Array}}
Mstat} <: Abstract2DStochasticModel where {Mstat<:Union{Nothing,stat}, PCollection<:Union{Vector,Array}}

Comment on lines +45 to +72
function write_particles_to_csv(wave_model, Δt)
sec=string(Int64(floor((wave_model.clock.time+Δt)/60)))
dec=string(Int64(floor(10*((wave_model.clock.time+Δt)/60-floor((wave_model.clock.time+Δt)/60)))))
save_path = wave_model.plot_savepath

nParticles = wave_model.n_particles_launch
@info wave_model.clock.time/60

parts = wave_model.ParticleCollection[(end-nParticles+1):end]

logE = zeros(nParticles)
cx = zeros(nParticles)
cy = zeros(nParticles)
x = zeros(nParticles)
y = zeros(nParticles)
for i in 1:nParticles
logE[i] = parts[i].ODEIntegrator[1]
cx[i] = parts[i].ODEIntegrator[2]
cy[i] = parts[i].ODEIntegrator[3]
x[i] = parts[i].ODEIntegrator[4]
y[i] = parts[i].ODEIntegrator[5]
end

data = DataFrame(id=1:nParticles, logE = logE, cx = cx, cy = cy, x = x, y = y)
data2 = Tables.table(transpose(wave_model.State[:, :, 1]))
CSV.write(save_path*"/data/particles_"*sec*","*dec*".csv", data)
CSV.write(save_path*"/data/mesh_values_"*sec*","*dec*".csv", data2)
end
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should move to Utils/InputOutput.jl

Comment on lines 74 to 75
"""
time_step!(model, Δt; callbacks=nothing)
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

delete this method since you don't do 1D.

Comment on lines +92 to +124
if sim.model.save_particles
save_path = sim.model.plot_savepath
filename = save_path*"/data/simu_infos.csv"

Nx = sim.model.grid.Nx
Ny = sim.model.grid.Ny
xmin = sim.model.grid.xmin
xmax = sim.model.grid.xmax
ymin = sim.model.grid.ymin
ymax = sim.model.grid.ymax
lne_source = sim.model.ODEdefaults.lne
c_x_source = sim.model.ODEdefaults.c̄_x
c_y_source = sim.model.ODEdefaults.c̄_y
x_source = sim.model.ODEdefaults.x
y_source = sim.model.ODEdefaults.y
angular_spread_source = sim.model.ODEdefaults.angular_σ
Δt = sim.Δt
stop_time = sim.stop_time

data = DataFrame(Nx=Nx, Ny=Ny, xmin=xmin,
xmax=xmax, ymin=ymin, ymax=ymax,
lne_source=lne_source, c_x_source=c_x_source,
c_y_source=c_y_source, x_source=x_source,y_source=y_source,
angular_spread_source=angular_spread_source,
Δt=Δt, stop_time=stop_time)
CSV.write(filename, data)

filename2 = save_path*"/data/sigma.csv"

covariance_init = sim.model.proba_covariance_init
data2 = DataFrame(covariance_init, :auto)
CSV.write(filename2, data2)
end
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Write this in a function and put it in the same place were write_particle_to_cvs is. then just call the function here.

Comment on lines +209 to +214
if sim.model.plot_steps
plot_state_and_error_points(sim, gridnotes)
sec=string(Int64(floor((sim.model.clock.time)/60)))
dec=string(Int64(floor(10*(sim.model.clock.time/60-floor((sim.model.clock.time)/60)))))
plt.savefig(joinpath([save_path*"/plots/", "energy_plot_no_spread_"*sec*","*dec*".png"]))
end
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

write this in a function as well.

usually the initilization uses wind constitions to seed the particles.
"""
function init_particles!(model::Abstract2DModel; defaults::PP=nothing, verbose::Bool=false) where {PP<:Union{ParticleDefaults1D,ParticleDefaults2D,Nothing}}
function init_particles!(model::Abstract2DModel; defaults::PP=nothing, verbose::Bool=false) where {PP<:Union{ParticleDefaults,Array{Any,1},Nothing}}
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

change model type here so this is then a new function type and doesn't conflict with the other code

# @info ParticleCollection
model.ParticleCollection = ParticleCollection

if defaults isa ParticleDefaults
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe clean this up and comment your code a bit so one has an idea what you do

Comment on lines +312 to +319
gridnotes = TwoDGridNotes(model.grid)

SeedParticle2D(
model.State, ij,
model.ODEsystem, defaults, model.ODEsettings,
model.grid.stats, model.grid.ProjetionKernel, model.grid.PropagationCorrection,
ij_mesh, ij_wind,
model.ODEsettings.timestep,
model.boundary, model.periodic_boundary)
ParticleCollection = []
SeedParticle_i = SeedParticle_mapper(SeedParticle2D!,
ParticleCollection, model.State,
model.ODEsystem, nothing, model.ODEsettings,
gridnotes, model.winds, model.ODEsettings.timestep,
model.boundary, model.periodic_boundary)
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This may not work anymore since I changed the type of the ParticleCollection.. we have to see ..

@mochell mochell self-requested a review November 5, 2025 19:19
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants