From 94477a8dd78abe2b8d10f64c2270d74fe4357918 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Sun, 7 Jul 2024 15:33:49 +0200 Subject: [PATCH 1/7] particles docs --- docs/src/particles.md | 38 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 36 insertions(+), 2 deletions(-) diff --git a/docs/src/particles.md b/docs/src/particles.md index 85ed8304..49d5bfa5 100644 --- a/docs/src/particles.md +++ b/docs/src/particles.md @@ -4,8 +4,42 @@ Particles stored in `CellArray`s objects from [`CellArrays.jl`](https://github.c # Particle objects -## Particles +There are three type of `AbstractParticles` types: + +- `Particles` is the basic type used to advect and track information on the whole domain of our model. The dimension of the `CellArrays` has to match the dimension of the model. +- `PassiveMarkers` is a set of particles where their initial position is defined by the user. These particles are only advected and are not meant to have any feedback with the simulation; their purpose is to track the history of any arbitray field(s) throughout the simulation. +- `MarkerChain` is a one or two dimensional chain of particles, used to track surfaces / interfaces. +## Particles +```julia +struct Particles{Backend,N,M,I,T1,T2} <: AbstractParticles + coords::NTuple{N,T1} + index::T2 + nxcell::I + max_xcell::I + min_xcell::I + np::I +end +``` +Where `coords` is a tuple containing the coordinates of the particles; `index` is an `BitArray` where `true` if in the correspondant `CellArray` there is an active particle, otherwise false; `nxcell`, `min_xcell`, `max_xcell` are the initial, minimum and maximum number of particles per cell; and, `np` is the initial number of particles. + ## Passive markers +```julia +struct MarkerChain{Backend,N,M,I,T1,T2,TV} <: AbstractParticles + coords::NTuple{N,T1} + index::T2 + cell_vertices::TV + max_xcell::I + min_xcell::I +end +``` +Where `coords` is a tuple containing the coordinates of the particles; `index` is an `BitArray` where `true` if in the correspondant `CellArray` there is an active particle, otherwise false; `cell_vertices` is the lower-left corner (`(x,)` in 2D, `(x,y)` in 3D) of the cell containing those particles;and, `min_xcell`, `max_xcell` are theminimum and maximum number of particles per cell. -## Marker chain \ No newline at end of file +## Marker chain +```julia +struct PassiveMarkers{Backend,T} <: AbstractParticles + coords::T + np::Int64 +end +``` +Where `coords` is a tuple containing the coordinates of the particles; and `np` is the number of passive markers. From a69902efe1391dfd384d924c235335fbdeea9d52 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Sun, 7 Jul 2024 15:34:13 +0200 Subject: [PATCH 2/7] typo --- docs/src/particles.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/particles.md b/docs/src/particles.md index 49d5bfa5..9e4953c4 100644 --- a/docs/src/particles.md +++ b/docs/src/particles.md @@ -21,7 +21,7 @@ struct Particles{Backend,N,M,I,T1,T2} <: AbstractParticles np::I end ``` -Where `coords` is a tuple containing the coordinates of the particles; `index` is an `BitArray` where `true` if in the correspondant `CellArray` there is an active particle, otherwise false; `nxcell`, `min_xcell`, `max_xcell` are the initial, minimum and maximum number of particles per cell; and, `np` is the initial number of particles. +Where `coords` is a tuple containing the coordinates of the particles; `index` is a `BitArray` where `true` if in the correspondant `CellArray` there is an active particle, otherwise false; `nxcell`, `min_xcell`, `max_xcell` are the initial, minimum and maximum number of particles per cell; and, `np` is the initial number of particles. ## Passive markers ```julia From 9182b5b3f1c38e1b91fef9ea6acd8b9b27846f46 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Sun, 7 Jul 2024 15:35:49 +0200 Subject: [PATCH 3/7] polishing --- docs/src/field_advection2D_MPI.md | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/docs/src/field_advection2D_MPI.md b/docs/src/field_advection2D_MPI.md index a7889adf..a1acf06f 100644 --- a/docs/src/field_advection2D_MPI.md +++ b/docs/src/field_advection2D_MPI.md @@ -99,20 +99,24 @@ Now start the simulation ```julia niter = 250 for it in 1:niter - advection!(particles, RungeKutta2(), V, (grid_vx, grid_vy), dt) # advect particles - move_particles!(particles, xvi, particle_args) # move particles in the memory - inject_particles!(particles, (pT, ), xvi) # inject particles if needed - particle2grid!(T, pT, xvi, particles) # interpolate particles to the grid + # advect particles + advection!(particles, RungeKutta2(), V, (grid_vx, grid_vy), dt) + # move particles in the memory + move_particles!(particles, xvi, particle_args) + # inject particles if needed + inject_particles!(particles, (pT, ), xvi) + # interpolate particles to the grid + particle2grid!(T, pT, xvi, particles) end ``` ## Visualization To visualize the results, we need to allocate a global array `T_v` and buffer arrays `T_nohalo` without the overlapping halo (here with `width = 1`) ```julia -nx_v = (size(T, 1)-2)*dims[1] # global size of `T` without halos -ny_v = (size(T, 2)-2)*dims[2] # global size of `T` without halos -T_v = zeros(nx_v, ny_v) # initialize global `T` -T_nohalo = @zeros(size(T).-2) # local `T` without overlapping halo +nx_v = (size(T, 1) - 2) * dims[1] # global size of `T` without halos +ny_v = (size(T, 2) - 2) * dims[2] # global size of `T` without halos +T_v = zeros(nx_v, ny_v) # initialize global `T` +T_nohalo = @zeros(size(T).-2) # local `T` without overlapping halo ``` Visualization with GLMakie.jl From 8456722e932a3e454aae9c62b25f89fc2a3c7e6d Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Sun, 7 Jul 2024 21:01:30 +0200 Subject: [PATCH 4/7] docs --- docs/make.jl | 11 ++++++----- docs/src/CellArrays.md | 5 ++++- docs/src/particles.md | 12 +++++++----- 3 files changed, 17 insertions(+), 11 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 47fd191d..a9a22dcd 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -3,13 +3,14 @@ push!(LOAD_PATH, "../src/") @info "Making documentation..." makedocs(; - sitename="JustPIC.jl", - authors="Albert de Montserrat and contributors", - modules=[JustPIC, JustPIC._2D, JustPIC._3D], - format=Documenter.HTML(; prettyurls=get(ENV, "CI", nothing) == "true"), # easier local build + sitename = "JustPIC.jl", + authors = "Albert de Montserrat and contributors", + modules = [JustPIC, JustPIC._2D, JustPIC._3D], + format = Documenter.HTML(; prettyurls=get(ENV, "CI", nothing) == "true"), # easier local build warnonly = Documenter.except(:footnote), - pages=[ + pages = [ "Home" => "index.md", + "Particles" =>"particles.md", "CellArrays" =>"CellArrays.md", "Examples" => [ "field_advection2D.md", diff --git a/docs/src/CellArrays.md b/docs/src/CellArrays.md index a1abee32..e95fa940 100644 --- a/docs/src/CellArrays.md +++ b/docs/src/CellArrays.md @@ -1,6 +1,7 @@ # Working with CellArrays -## Instantiating a CellArray +## Instantiating a `CellArray` + With the help of `ParallelStencil.jl` we can easily create a `CellArray` object. The `CellArray` object is a container that holds the data of a grid. The data is stored in small nD-arrays, and the grid is divided into cells. Each cell contains a number of elements. The `CellArray` object is used to store the data of the particles in the simulation. ```julia @@ -25,6 +26,8 @@ julia> CA = @fill(x, ni..., celldims = ncells, eltype = Float64) [20.0, 20.0] [20.0, 20.0] ``` +## Indexing a `CellArray` + We can access to the data of one `CellArray` by indexing a given grid cell. This will however instantiate a `StaticArray` object with the data of the cell. ```julia-repl diff --git a/docs/src/particles.md b/docs/src/particles.md index 9e4953c4..c2fd344b 100644 --- a/docs/src/particles.md +++ b/docs/src/particles.md @@ -1,8 +1,10 @@ -# Memory layout +# Particles + +## Memory layout Particles stored in `CellArray`s objects from [`CellArrays.jl`](https://github.com/omlins/CellArrays.jl) and are constantly sorted by their parent cell to avoid a loss of spatial locality with time. -# Particle objects +## Particle objects There are three type of `AbstractParticles` types: @@ -10,7 +12,7 @@ There are three type of `AbstractParticles` types: - `PassiveMarkers` is a set of particles where their initial position is defined by the user. These particles are only advected and are not meant to have any feedback with the simulation; their purpose is to track the history of any arbitray field(s) throughout the simulation. - `MarkerChain` is a one or two dimensional chain of particles, used to track surfaces / interfaces. -## Particles +### Simulation particles ```julia struct Particles{Backend,N,M,I,T1,T2} <: AbstractParticles coords::NTuple{N,T1} @@ -23,7 +25,7 @@ end ``` Where `coords` is a tuple containing the coordinates of the particles; `index` is a `BitArray` where `true` if in the correspondant `CellArray` there is an active particle, otherwise false; `nxcell`, `min_xcell`, `max_xcell` are the initial, minimum and maximum number of particles per cell; and, `np` is the initial number of particles. -## Passive markers +### Passive markers ```julia struct MarkerChain{Backend,N,M,I,T1,T2,TV} <: AbstractParticles coords::NTuple{N,T1} @@ -35,7 +37,7 @@ end ``` Where `coords` is a tuple containing the coordinates of the particles; `index` is an `BitArray` where `true` if in the correspondant `CellArray` there is an active particle, otherwise false; `cell_vertices` is the lower-left corner (`(x,)` in 2D, `(x,y)` in 3D) of the cell containing those particles;and, `min_xcell`, `max_xcell` are theminimum and maximum number of particles per cell. -## Marker chain +### Marker chain ```julia struct PassiveMarkers{Backend,T} <: AbstractParticles coords::T From 8d4a97565fb5c3a967b63ebfa7f75b85d003a91a Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Sun, 7 Jul 2024 22:16:44 +0200 Subject: [PATCH 5/7] docs --- docs/make.jl | 1 + docs/src/interpolations.md | 55 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 56 insertions(+) create mode 100644 docs/src/interpolations.md diff --git a/docs/make.jl b/docs/make.jl index a9a22dcd..12b0e4c9 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,6 +11,7 @@ makedocs(; pages = [ "Home" => "index.md", "Particles" =>"particles.md", + "Interpolations" =>"particles.md", "CellArrays" =>"CellArrays.md", "Examples" => [ "field_advection2D.md", diff --git a/docs/src/interpolations.md b/docs/src/interpolations.md new file mode 100644 index 00000000..1264648c --- /dev/null +++ b/docs/src/interpolations.md @@ -0,0 +1,55 @@ +# Particles interpolations + +For now, interpolations support only equidistant regular rectangular/cubic grids. Please file and issue if you would be interested in supporting other kind of grids. + +## Grid to particle + +Information from the grid is linearly interpolated to the particles. The one dimensional linear interpolation kernel (also referred as _lerp_) is defined as + +$v_{\text{p}} = t v_0 + (1 -t) v_1$ + +where the $t$, $v_0$, and $v_1$ are graphically described below. + + + +Numerically, it is more appropiately implemented as a double [fma](https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation) as it is slightly more accurate than a naive implementation: + +```julia +v_p = fma(t, v1, fma(-t, v0, v0)) +``` + +Bi- and tri-linear interpolation over a rectangular or cubic cells is thus nothing else than a chain of lerp kernels along the different dimensions of the cell. For example, the bilinear interpolation requires two lerps along the left and right sides of the cell, followed by a lerp on the horizontal direction; and trilinear interpolation requires two bilinear kernels and one lerp. + +N-linear interpolation is implemented in a recursive way to exploit compiler optimizations and reduce boilerplate code. Since in practical terms we will do up to tri-linear interpolation the maximum recursion depth is five, so that stack will never overflow. + +We can interpolate an arbitray field $F$ onto the particles with the `grid2particle` function: + +```julia +using JustPIC, JustPIC._2D +# define model domain +nxcell, max_xcell, min_xcell = 24, 30, 12 +nx = ny = 128 +Lx = Ly = 1.0 +xvi = range(0, Lx, length=n), range(0, Ly, length=n) +# field F at the grid +F = [y for x in xv, y in yv] +# instantiate empty `CellArray` +Fp, = init_cell_arrays(particles, Val(1)); +# interpolate F onto Fp +grid2particle!(Fp, xvi, F, particles); +``` + +## Particle to grid + +Information on the particles is currently interpolated to the vertices of the grid cells with an inverse distance weighting interpolant + +$v_{i,j} = \frac{\sum^N_{k=1} \omega_k v_k}{\sum^n_{k=1} \omega_k}$ + +where the weight is $\omega_i = d^{-n}$, with $d$ being the distance between the particle and the node, and $n$ a integer number. + +In shared memory environments (e.g. OpenMP parallelism or GPU programming), this particular interpolation often requires the use of atomic operations due to race conditions. In `JustPIC` we do avoid using atomic operations by parallelising over the grid nodes instead of the particles, i.e. we iterated over the nodes, with an inner iteration over the `CellArrays` neighbouring the $i$-th grid node. + +This inteprolation is done with the `particle2grid!` function: +```julia-repl +julia> particle2grid!(F, Fp, xvi, particles) +``` \ No newline at end of file From 36d069460ca13dd0cbd8c660c90f3df8bb150b3e Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Sun, 7 Jul 2024 22:18:51 +0200 Subject: [PATCH 6/7] typos --- docs/src/interpolations.md | 4 ++-- docs/src/particles.md | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/interpolations.md b/docs/src/interpolations.md index 1264648c..622772c7 100644 --- a/docs/src/interpolations.md +++ b/docs/src/interpolations.md @@ -12,7 +12,7 @@ where the $t$, $v_0$, and $v_1$ are graphically described below. -Numerically, it is more appropiately implemented as a double [fma](https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation) as it is slightly more accurate than a naive implementation: +Numerically, it is more appropriately implemented as a double [fma](https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation) as it is slightly more accurate than a naive implementation: ```julia v_p = fma(t, v1, fma(-t, v0, v0)) @@ -22,7 +22,7 @@ Bi- and tri-linear interpolation over a rectangular or cubic cells is thus nothi N-linear interpolation is implemented in a recursive way to exploit compiler optimizations and reduce boilerplate code. Since in practical terms we will do up to tri-linear interpolation the maximum recursion depth is five, so that stack will never overflow. -We can interpolate an arbitray field $F$ onto the particles with the `grid2particle` function: +We can interpolate an arbitrary field $F$ onto the particles with the `grid2particle` function: ```julia using JustPIC, JustPIC._2D diff --git a/docs/src/particles.md b/docs/src/particles.md index c2fd344b..83ae58f9 100644 --- a/docs/src/particles.md +++ b/docs/src/particles.md @@ -23,7 +23,7 @@ struct Particles{Backend,N,M,I,T1,T2} <: AbstractParticles np::I end ``` -Where `coords` is a tuple containing the coordinates of the particles; `index` is a `BitArray` where `true` if in the correspondant `CellArray` there is an active particle, otherwise false; `nxcell`, `min_xcell`, `max_xcell` are the initial, minimum and maximum number of particles per cell; and, `np` is the initial number of particles. +Where `coords` is a tuple containing the coordinates of the particles; `index` is a `BitArray` where `true` if in the correspondent `CellArray` there is an active particle, otherwise false; `nxcell`, `min_xcell`, `max_xcell` are the initial, minimum and maximum number of particles per cell; and, `np` is the initial number of particles. ### Passive markers ```julia @@ -35,7 +35,7 @@ struct MarkerChain{Backend,N,M,I,T1,T2,TV} <: AbstractParticles min_xcell::I end ``` -Where `coords` is a tuple containing the coordinates of the particles; `index` is an `BitArray` where `true` if in the correspondant `CellArray` there is an active particle, otherwise false; `cell_vertices` is the lower-left corner (`(x,)` in 2D, `(x,y)` in 3D) of the cell containing those particles;and, `min_xcell`, `max_xcell` are theminimum and maximum number of particles per cell. +Where `coords` is a tuple containing the coordinates of the particles; `index` is an `BitArray` where `true` if in the correspondent `CellArray` there is an active particle, otherwise false; `cell_vertices` is the lower-left corner (`(x,)` in 2D, `(x,y)` in 3D) of the cell containing those particles;and, `min_xcell`, `max_xcell` are theminimum and maximum number of particles per cell. ### Marker chain ```julia From 30147238f38ccb3f49149e3521f3ab7fa75dd8d5 Mon Sep 17 00:00:00 2001 From: albert-de-montserrat Date: Mon, 8 Jul 2024 09:21:04 +0200 Subject: [PATCH 7/7] typo --- docs/src/particles.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/particles.md b/docs/src/particles.md index 83ae58f9..840b22a1 100644 --- a/docs/src/particles.md +++ b/docs/src/particles.md @@ -9,7 +9,7 @@ Particles stored in `CellArray`s objects from [`CellArrays.jl`](https://github.c There are three type of `AbstractParticles` types: - `Particles` is the basic type used to advect and track information on the whole domain of our model. The dimension of the `CellArrays` has to match the dimension of the model. -- `PassiveMarkers` is a set of particles where their initial position is defined by the user. These particles are only advected and are not meant to have any feedback with the simulation; their purpose is to track the history of any arbitray field(s) throughout the simulation. +- `PassiveMarkers` is a set of particles where their initial position is defined by the user. These particles are only advected and are not meant to have any feedback with the simulation; their purpose is to track the history of any arbitrary field(s) throughout the simulation. - `MarkerChain` is a one or two dimensional chain of particles, used to track surfaces / interfaces. ### Simulation particles