From 4d8cfd4f7f4cbaa4a26e6450bb29f47137306c63 Mon Sep 17 00:00:00 2001 From: tduretz Date: Wed, 25 Sep 2024 15:40:14 +0200 Subject: [PATCH 1/5] add restart example --- scripts/SolCx/pureshear_ALE_restart.jl | 215 +++++++++++++++++++++++++ scripts/pureshear_ALE.jl | 100 +++++++----- 2 files changed, 278 insertions(+), 37 deletions(-) create mode 100644 scripts/SolCx/pureshear_ALE_restart.jl diff --git a/scripts/SolCx/pureshear_ALE_restart.jl b/scripts/SolCx/pureshear_ALE_restart.jl new file mode 100644 index 00000000..85e78cd0 --- /dev/null +++ b/scripts/SolCx/pureshear_ALE_restart.jl @@ -0,0 +1,215 @@ +using Statistics, LinearAlgebra, Printf, Base.Threads, CairoMakie, JLD2 +const year = 365*3600*24 +const USE_GPU = false + +using JustPIC, JustPIC._2D +const backend = JustPIC.CPUBackend + +const ALE = true +last_step = 100 +restart = false + +using ParallelStencil +using ParallelStencil.FiniteDifferences2D +@static if USE_GPU + @init_parallel_stencil(CUDA, Float64, 2) +else + @init_parallel_stencil(Threads, Float64, 2) +end + +@parallel_indices (I...) function InitialFieldsParticles!(phases, px, py, index) + for ip in cellaxes(phases) + # quick escape + @index(index[ip, I...]) == 0 && continue + x = @index px[ip, I...] + y = @index py[ip, I...] + if x 0 + extrema phase ratio @ vertices = $(extrema(sum(phase_ratios.vertex.data, dims=2))) + extrema phase ratio @ centers = $(extrema(sum(phase_ratios.center.data, dims=2))) + ") + + # Time step + Δt = C * min(Δ...) / max(maximum(abs.(V.x)), maximum(abs.(V.y))) + @show Δt + + # Create necessary tuples + grid_vx = (verts.x, cents_ext.y) + grid_vy = (cents_ext.x, verts.y) + Vxc = 0.5*(V.x[1:end-1,2:end-1] .+ V.x[2:end-0,2:end-1]) + Vyc = 0.5*(V.y[2:end-1,1:end-1] .+ V.y[2:end-1,2:end-0]) + Vmag = sqrt.(Vxc.^2 .+ Vyc.^2) + + for it=it0:Nt + + @printf("Step %05d, #particles = %06d\n", it, Npart) + + t += Δt + + # advection!(particles, RungeKutta2(), values(V), (grid_vx, grid_vy), Δt) + # advection_LinP!(particles, RungeKutta2(), values(V), (grid_vx, grid_vy), Δt) + advection_MQS!(particles, RungeKutta2(), values(V), (grid_vx, grid_vy), Δt); + move_particles!(particles, values(verts), particle_args); + inject_particles_phase!(particles, phases, (), (), values(verts)) + phase_ratios_vertex!(phase_ratios, particles, values(verts), phases); + phase_ratios_center!(phase_ratios, particles, values(cents), phases); + Npart = sum(particles.index.data) + + if ALE + xlims[1] += xlims[1]*ε̇bg*Δt + xlims[2] += xlims[2]*ε̇bg*Δt + ylims[1] -= ylims[1]*ε̇bg*Δt + ylims[2] -= ylims[2]*ε̇bg*Δt + L = ( x =(xlims[2]-xlims[1]), y =(ylims[2]-ylims[1]) ) + Δ = (x=L.x/Nc.x, y=L.y/Nc.y ) + cents_ext = ( + x = LinRange(xlims[1]-Δ.x/2, xlims[2]+Δ.x/2, Nc.x+2), + y = LinRange(ylims[1]-Δ.y/2, ylims[2]+Δ.y/2, Nc.y+2), + ) + verts = ( + x = LinRange(xlims[1], xlims[2], Nc.x+1), + y = LinRange(ylims[1], ylims[2], Nc.y+1), + ) + grid_vx = (verts.x, cents_ext.y) + grid_vy = (cents_ext.x, verts.y) + Δt = C * min(Δ...) / max(maximum(abs.(V.x)), maximum(abs.(V.y))) + @parallel SetVelocity(V, verts, ε̇bg) + move_particles!(particles, values(verts), (phases,)) + phase_ratios_vertex!(phase_ratios, particles, values(verts), phases) + end + + if mod(it, Nout) == 0 || it==1 + @show + particle_density = [sum(p) for p in particles.index] + # visualisation + p = particles.coords + ppx, ppy = p + pxv = ppx.data[:] + pyv = ppy.data[:] + clr = phases.data[:] + idxv = particles.index.data[:] + f = Figure() + ax = Axis(f[1, 1], title="Particles", aspect=L.x/L.y) + scatter!(ax, Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]), colormap=:roma, markersize=2) + xlims!(ax, verts.x[1], verts.x[end]) + ylims!(ax, verts.y[1], verts.y[end]) + display(f) + end + + # Save breakpoint + if it==100 + @show file = @sprintf("./Breakpoint%05d.jld2", Nt) + jldsave(file; particles, phases, phase_ratios, particle_args, xlims, ylims, t) + end + end + + return nothing +end + +main() \ No newline at end of file diff --git a/scripts/pureshear_ALE.jl b/scripts/pureshear_ALE.jl index 8a831778..078ce0ef 100644 --- a/scripts/pureshear_ALE.jl +++ b/scripts/pureshear_ALE.jl @@ -1,6 +1,6 @@ using Statistics, LinearAlgebra, Printf, Base.Threads, CairoMakie const year = 365*3600*24 -const USE_GPU = true +const USE_GPU = false using JustPIC, JustPIC._2D const backend = JustPIC.CPUBackend @@ -30,38 +30,52 @@ end return nothing end +@parallel_indices (I...) function SetVelocity(V, verts, ε̇bg) + + if I[1]<=size(V.x,1) && I[2]<=size(V.x,2) + V.x[I...] = verts.x[I[1]]*ε̇bg + end + + if I[1]<=size(V.y,1) && I[2]<=size(V.y,2) + V.y[I...] = -verts.y[I[2]]*ε̇bg + end + +return nothing +end + + + function main() @printf("Running on %d thread(s)\n", nthreads()) - L = (x=1., y=1.) - Nc = (x=41, y=41 ) - Nv = (x=Nc.x+1, y=Nc.y+1 ) - Δ = (x=L.x/Nc.x, y=L.y/Nc.y ) + + # Parameters + L = (x=1., y=1.) + Nc = (x=41, y=41 ) + Nv = (x=Nc.x+1, y=Nc.y+1 ) + Δ = (x=L.x/Nc.x, y=L.y/Nc.y ) Nt = 100 - Nout = 1 + Nout = 10 C = 0.25 + # Model extent verts = (x=LinRange(-L.x/2, L.x/2, Nv.x), y=LinRange(-L.y/2, L.y/2, Nv.y)) cents = (x=LinRange(-Δ.x/2+L.x/2, L.x/2-Δ.x/2, Nc.x), y=LinRange(-Δ.y/2+L.y/2, L.y+Δ.y/2-L.y/2, Nc.y)) cents_ext = (x=LinRange(-Δ.x/2-L.x/2, L.x/2+Δ.x/2, Nc.x+2), y=LinRange(-Δ.y/2-L.y/2, L.y+Δ.y/2+L.y/2, Nc.y+2)) + xlims = [verts.x[1] , verts.x[end]] + ylims = [verts.y[1] , verts.y[end]] - size_x = (Nc.x+1, Nc.y+2) - size_y = (Nc.x+2, Nc.y+1) - + # Arrays + size_x = (Nc.x+1, Nc.y+2) + size_y = (Nc.x+2, Nc.y+1) V = ( x = @zeros(size_x), y = @zeros(size_y), ) # Set velocity field - ε̇bg = -1.0 - for i=1:size(V.x,1), j=1:size(V.x,2) - CUDA.@allowscalar V.x[i,j] = verts.x[i]*ε̇bg - end - - for i=1:size(V.y,1), j=1:size(V.y,2) - CUDA.@allowscalar V.y[i,j] = -verts.y[j]*ε̇bg - end + ε̇bg = 1. + @parallel SetVelocity(V, verts, ε̇bg) # Initialize particles ------------------------------- nxcell, max_xcell, min_xcell = 12, 24, 5 @@ -112,39 +126,51 @@ function main() move_particles!(particles, values(verts), particle_args); inject_particles_phase!(particles, phases, (), (), values(verts)) phase_ratios_vertex!(phase_ratios, particles, values(verts), phases); - phase_ratios_center!(phase_ratios, particles, values(verts), phases); + phase_ratios_center!(phase_ratios, particles, values(cents), phases); if ALE - @show L = (x=1.0+ε̇bg*t, y=1.0-ε̇bg*t) + xlims[1] += xlims[1]*ε̇bg*Δt + xlims[2] += xlims[2]*ε̇bg*Δt + ylims[1] -= ylims[1]*ε̇bg*Δt + ylims[2] -= ylims[2]*ε̇bg*Δt + @show L = ( x =(xlims[2]-xlims[1]), y =(ylims[2]-ylims[1]) ) Δ = (x=L.x/Nc.x, y=L.y/Nc.y ) - verts = (x=LinRange(-L.x/2, L.x/2, Nv.x), y=LinRange(-L.y/2, L.y/2, Nv.y)) - cents_ext = (x=LinRange(-Δ.x/2-L.x/2, L.x/2+Δ.x/2, Nc.x+2), y=LinRange(-Δ.y/2-L.y/2, L.y+Δ.y/2+L.y/2, Nc.y+2)) + cents_ext = ( + x = LinRange(xlims[1]-Δ.x/2, xlims[2]+Δ.x/2, Nc.x+2), + y = LinRange(ylims[1]-Δ.y/2, ylims[2]+Δ.y/2, Nc.y+2), + ) + verts = ( + x = LinRange(xlims[1], xlims[2], Nc.x+1), + y = LinRange(ylims[1], ylims[2], Nc.y+1), + ) grid_vx = (verts.x, cents_ext.y) grid_vy = (cents_ext.x, verts.y) Δt = C * min(Δ...) / max(maximum(abs.(V.x)), maximum(abs.(V.y))) + @parallel SetVelocity(V, verts, ε̇bg) + move_particles!(particles, values(verts), (phases,)) + phase_ratios_vertex!(phase_ratios, particles, values(verts), phases) end - if mod(it,Nout) == 0 || it==1 + if mod(it, Nout) == 0 || it==1 @show Npart = sum(particles.index.data) particle_density = [sum(p) for p in particles.index] @show size(particle_density) - end + # visualisation + p = particles.coords + ppx, ppy = p + pxv = ppx.data[:] + pyv = ppy.data[:] + clr = phases.data[:] + idxv = particles.index.data[:] + f = Figure() + ax = Axis(f[1, 1], title="Particles", aspect=L.x/L.y) + scatter!(ax, Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]), colormap=:roma, markersize=2) + xlims!(ax, verts.x[1], verts.x[end]) + ylims!(ax, verts.y[1], verts.y[end]) + display(f) + end end - # Plots - p = particles.coords - ppx, ppy = p - pxv = ppx.data[:] - pyv = ppy.data[:] - clr = phases.data[:] - idxv = particles.index.data[:] - f = Figure() - ax = Axis(f[1, 1], title="Particles", aspect=L.x/L.y) - scatter!(ax, Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]), colormap=:roma, markersize=2) - xlims!(ax, verts.x[1], verts.x[end]) - ylims!(ax, verts.y[1], verts.y[end]) - display(f) - return nothing end From 28e98ccf5e527beedbd7026df1ccd3a46f9d8ddd Mon Sep 17 00:00:00 2001 From: Thibault Duretz <48455309+tduretz@users.noreply.github.com> Date: Thu, 26 Sep 2024 08:06:00 +0200 Subject: [PATCH 2/5] Update scripts/SolCx/pureshear_ALE_restart.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- scripts/SolCx/pureshear_ALE_restart.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/scripts/SolCx/pureshear_ALE_restart.jl b/scripts/SolCx/pureshear_ALE_restart.jl index 85e78cd0..3cfce5c6 100644 --- a/scripts/SolCx/pureshear_ALE_restart.jl +++ b/scripts/SolCx/pureshear_ALE_restart.jl @@ -49,18 +49,18 @@ function main() # Load Breakpoint data if restart - file = @sprintf("./Breakpoint%05d.jld2", last_step) - data = load(file) + file = @sprintf("./Breakpoint%05d.jld2", last_step) + data = load(file) particles = data["particles"] phases = data["phases"] phase_ratios = data["phase_ratios"] particle_args = data["particle_args"] - xlims = data["xlims"] - ylims = data["ylims"] - t = data["t"] - Nt = last_step + 100 - it0 = last_step + 1 - ε̇bg = -1. + xlims = data["xlims"] + ylims = data["ylims"] + t = data["t"] + Nt = last_step + 100 + it0 = last_step + 1 + ε̇bg = -1. else xlims = [-0.5 , 0.5] ylims = [-0.5 , 0.5] From 21ef2d568c2f951f42fc889a4840856a94f72d7b Mon Sep 17 00:00:00 2001 From: Thibault Duretz <48455309+tduretz@users.noreply.github.com> Date: Thu, 26 Sep 2024 08:06:51 +0200 Subject: [PATCH 3/5] Update scripts/SolCx/pureshear_ALE_restart.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- scripts/SolCx/pureshear_ALE_restart.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/SolCx/pureshear_ALE_restart.jl b/scripts/SolCx/pureshear_ALE_restart.jl index 3cfce5c6..3cec3cc4 100644 --- a/scripts/SolCx/pureshear_ALE_restart.jl +++ b/scripts/SolCx/pureshear_ALE_restart.jl @@ -126,7 +126,7 @@ function main() phase_ratios = JustPIC._2D.PhaseRatios(backend, 2, values(Nc)); phase_ratios_vertex!(phase_ratios, particles, values(verts), phases) - phase_ratios_center!(phase_ratios, particles, values(verts), phases) + phase_ratios_center!(phase_ratios, particles, values(cents), phases) Npart = sum(particles.index.data) println(" From 0d8b7ac3014b03395a83e957245d854f8c1b4a68 Mon Sep 17 00:00:00 2001 From: Thibault Duretz <48455309+tduretz@users.noreply.github.com> Date: Thu, 26 Sep 2024 08:07:27 +0200 Subject: [PATCH 4/5] Update scripts/SolCx/pureshear_ALE_restart.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- scripts/SolCx/pureshear_ALE_restart.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/scripts/SolCx/pureshear_ALE_restart.jl b/scripts/SolCx/pureshear_ALE_restart.jl index 3cec3cc4..f7c234d7 100644 --- a/scripts/SolCx/pureshear_ALE_restart.jl +++ b/scripts/SolCx/pureshear_ALE_restart.jl @@ -129,11 +129,6 @@ function main() phase_ratios_center!(phase_ratios, particles, values(cents), phases) Npart = sum(particles.index.data) - println(" - it => 0 - extrema phase ratio @ vertices = $(extrema(sum(phase_ratios.vertex.data, dims=2))) - extrema phase ratio @ centers = $(extrema(sum(phase_ratios.center.data, dims=2))) - ") # Time step Δt = C * min(Δ...) / max(maximum(abs.(V.x)), maximum(abs.(V.y))) From 22be8cdc63fc25cff3a55c34f9a7a9e6a7d82ec6 Mon Sep 17 00:00:00 2001 From: tduretz Date: Thu, 26 Sep 2024 09:18:04 +0200 Subject: [PATCH 5/5] add all suggestions --- scripts/{SolCx => }/pureshear_ALE_restart.jl | 31 ++++++++++---------- 1 file changed, 15 insertions(+), 16 deletions(-) rename scripts/{SolCx => }/pureshear_ALE_restart.jl (92%) diff --git a/scripts/SolCx/pureshear_ALE_restart.jl b/scripts/pureshear_ALE_restart.jl similarity index 92% rename from scripts/SolCx/pureshear_ALE_restart.jl rename to scripts/pureshear_ALE_restart.jl index f7c234d7..b6c644fc 100644 --- a/scripts/SolCx/pureshear_ALE_restart.jl +++ b/scripts/pureshear_ALE_restart.jl @@ -5,10 +5,6 @@ const USE_GPU = false using JustPIC, JustPIC._2D const backend = JustPIC.CPUBackend -const ALE = true -last_step = 100 -restart = false - using ParallelStencil using ParallelStencil.FiniteDifferences2D @static if USE_GPU @@ -43,13 +39,13 @@ return nothing end -function main() +function main(ALE, restart, last_step) @printf("Running on %d thread(s)\n", nthreads()) - # Load Breakpoint data + # Load checkpoint data if restart - file = @sprintf("./Breakpoint%05d.jld2", last_step) + file = @sprintf("./Checkpoint%05d.jld2", last_step) data = load(file) particles = data["particles"] phases = data["phases"] @@ -135,11 +131,8 @@ function main() @show Δt # Create necessary tuples - grid_vx = (verts.x, cents_ext.y) - grid_vy = (cents_ext.x, verts.y) - Vxc = 0.5*(V.x[1:end-1,2:end-1] .+ V.x[2:end-0,2:end-1]) - Vyc = 0.5*(V.y[2:end-1,1:end-1] .+ V.y[2:end-1,2:end-0]) - Vmag = sqrt.(Vxc.^2 .+ Vyc.^2) + grid_vx = verts.x, cents_ext.y + grid_vy = cents_ext.x, verts.y for it=it0:Nt @@ -190,16 +183,16 @@ function main() clr = phases.data[:] idxv = particles.index.data[:] f = Figure() - ax = Axis(f[1, 1], title="Particles", aspect=L.x/L.y) + ax = Axis(f[1, 1], title="Particles", aspect=L.x/L.y, xlabel="x", ylabel="y") scatter!(ax, Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]), colormap=:roma, markersize=2) xlims!(ax, verts.x[1], verts.x[end]) ylims!(ax, verts.y[1], verts.y[end]) display(f) end - # Save breakpoint + # Save checkpoint if it==100 - @show file = @sprintf("./Breakpoint%05d.jld2", Nt) + @show file = @sprintf("./Checkpoint%05d.jld2", Nt) jldsave(file; particles, phases, phase_ratios, particle_args, xlims, ylims, t) end end @@ -207,4 +200,10 @@ function main() return nothing end -main() \ No newline at end of file +################################### + +ALE = true +last_step = 100 +restart = false + +main(ALE, restart, last_step) \ No newline at end of file