Merge pull request #119 from JuliaGeodynamics/adm/moreDocs

Add more docs

docs/make.jl

@@ -3,13 +3,15 @@ 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",
+        "Interpolations" =>"particles.md",
         "CellArrays" =>"CellArrays.md",
         "Examples" => [
             "field_advection2D.md",

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 

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

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.
+
+<img src="assets/lerp.png" width="250"  />
+
+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))
+```
+
+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 arbitrary 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)
+```

docs/src/particles.md

@@ -1,11 +1,47 @@
-# 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:
 
-## Particles
+- `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 arbitrary field(s) throughout the simulation.
+- `MarkerChain` is a one or two dimensional chain of particles, used to track surfaces / interfaces.
 
-## Passive markers
+### Simulation 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 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
+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 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
+### 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.