[Documentation] Improving documentation

README.md

@@ -7,6 +7,8 @@
 
 Hey, that's Chmy.jl 🚀
 
+Chmy.jl provides a backend-agnostic toolkit for finite difference computations on multi-dimensional computational staggered grids. Chmy.jl features task-based distributed memory parallelisation capabilities.
+
 ## Documentation
 Checkout the [documentation](https://PTsolvers.github.io/Chmy.jl/dev) for API reference and usage.
 

docs/make.jl

@@ -14,6 +14,21 @@ makedocs(
     warnonly = [:missing_docs],
     pages = Any[
         "Home" => "index.md",
+        "Concepts" => Any["concepts/overview.md",
+                        "concepts/architectures.md",
+                        "concepts/grids.md",
+                        "concepts/fields.md",
+                        "concepts/bc.md",
+                        "concepts/grid_operators.md",
+                        "concepts/workers.md", 
+                        "concepts/kernels.md"
+        ],
+        "Examples" => Any["examples/overview.md",
+                          "examples/diffusion_2d.md",
+                          "examples/diffusion_2d_mpi.md",
+                          "examples/diffusion_2d_perf.md",
+                          "examples/batcher.md"
+        ],
         "Usage" => Any["usage/runtests.md"],
         "Library" => Any["lib/modules.md"]
     ]

docs/src/concepts/architectures.md

@@ -0,0 +1,43 @@
+# Architectures
+
+The abstract type `Architecture` can be concretized to be either an architecture of a concrete type `SingleDeviceArchitecture` or a `DistributedArchitecture` for MPI usage.
+
+## Single Device Architecture
+
+An object being the type of `SingleDeviceArchitecture` meaning that it contains the handle to a single CPU or GPU device with a respetive backend, where a device can be specified by an integer ID.
+
+### Backend Selection & Architecture Initialization
+
+We currently support CPU architectures, as well as CUDA and ROC backends for Nvidia and AMD GPUs through a thin wrapper around the [`KernelAbstractions.jl`](https://github.com/JuliaGPU/KernelAbstractions.jl) for users to select desirable backends.
+
+```julia
+# Default with CPU
+arch = Arch(CPU())
+```
+
+```julia
+# using CUDA
+arch = Arch(CUDABackend())
+```
+
+```julia
+# using AMDGPU
+arch = Arch(ROCBackend())
+```
+
+At the beginning of program, one may specify the backend and initialize the architecture they desire to use. The initialized `arch` variable will be required explicitly at creation of some entiries such as grids and kernel launchers.
+
+### Specifying Stream Priority
+
+For advanced users, we provide a convenient wrapper `activate!(arch::SingleDeviceArchitecture; priority=:normal)` for specifying the stream priority owned by the task one is executing. The stream priority will be set to `:normal` by default, where `:low` and `:high` are also possible options given that the target backend has priority control over streams implemented.
+
+For an example implementation, see [`AMDGPU.priority!`](https://amdgpu.juliagpu.org/stable/streams/#AMDGPU.priority!).
+
+
+## Distributed Architecture
+
+Our distributed architecture builds upon the abstraction of having GPU clusters that build on the same GPU architecture. Note that in general, GPU clusters may be equipped with hardware from different vendors, incorporating different types of GPUs to exploit their unique capabilities for specific tasks.
+
+With this abstraction of having homogeneous clusters, an object of type `DistributedArchitecture{ChildArch,Topo}` can have an uniform `child_arch` property representing all child architectures that comprise the distributed architecture. We also provide conveninent getters `get_backend` and `get_device` for obtaining the underlying backend and device of child architectures.
+
+Additionally, a distributed architecture also has a `topo` property that specifies that underlying MPI topology used. We currently support the creation of MPI communicator with Cartesian topology information attached.

docs/src/concepts/bc.md

@@ -0,0 +1,46 @@
+# Boundary Conditions
+
+Using [Chmy.jl](https://github.com/PTsolvers/Chmy.jl), we aim to study partial differential equations (PDEs) arising from physical or engineering problems. Generally, PDEs possess infinitely many solutions. Therefore, to obtain a unique solution, additional initial or boundary conditions are necessary for the model problem to be well-posed, ensuring the existence and uniqueness of a stable solution.
+
+We provide a small overview for boundary conditions that one often encounter. Followingly, we consider the unknown function $u : \Omega \mapsto \mathbb{R}$ defined on some  bounded computational domain $\Omega \subset \mathbb{R}^d$ in a $d$-dimensional space. With the domain boundary denoted by $\partial \Omega$, we have some function $g : \partial \Omega \mapsto \mathbb{R}$ prescribed on the boundary.
+
+
+| Type    | Form | Example |
+|:------------|:------------|:---------|
+| Dirichlet | $u = g$ on $\partial \Omega$ | In fluid dynamics, the no-slip condition for viscous fluids states that at a solid boundary the fluid has zero velocity relative to the boundary. |
+| Neumann | $\partial_\nu u = g$ on $\partial \Omega$, where $\nu$ is the outer normal vector to $\Omega$ | It specifies the values in which the derivative of a solution is applied within the boundary of the domain. An application in thermodynamics is a prescribed heat flux from a surface, which serves as boundary condition |
+| Robin  |  $u + \alpha \partial_\nu u = g$ on $\partial \Omega$, where $\alpha \in \mathbb{R}$.  | Also called impedance boundary conditions from their application in electromagnetic problems |
+
+
+
+
+## Applying Boundary Conditions
+
+Followingly, we describe the syntax in [Chmy.jl](https://github.com/PTsolvers/Chmy.jl) for launching kernels that impose boundary conditions on some `field` that is well-defined on a `grid` with backend specified through `arch`. For Dirichlet and Neumann boundary conditions, they are referred to as homogeneous if $g = 0$, otherwise they are non-homogeneous if $g = v$ holds, for some $v\in \mathbb{R}$.
+
+|     | Homogeneous | Non-homogeneous |
+|:------------|:------------|:------------|
+| Dirichlet | `bc!(arch, grid, field => Dirichlet())` | `bc!(arch, grid, field => Dirichlet(v))` |
+| Neumann | `bc!(arch, grid, field => Neumann())` | `bc!(arch, grid, field => Neumann(v))` |
+
+### Mixed Boundary Conditions
+
+In the code example above, by specifying boundary conditions using syntax such as `field => Neumann()`, we essentially launch a kernel that impose the Neumann boundary condition on the entire domain boundary $\partial \Omega$. More often, one may be interested in prescribing different boundary conditions on different parts of $\partial \Omega$.
+
+
+The following figure showcases a 2D square domain $\Omega$ with different boundary conditions applied on each side:
+
+- The top boundary (red) is a Dirichlet boundary condition where $u = a$.
+- The bottom boundary (blue) is also a Dirichlet boundary condition where $u = b$.
+- The left and right boundaries (green) are Neumann boundary conditions where $\frac{\partial u}{\partial y} = 0$.
+
+
+```@raw html
+<img src="../assets/mixed_bc_example.png" width="60%"/>
+```
+
+To launch a kernel that satisfies these boundary conditions in Chmy.jl, you can use the following code:
+
+```julia
+bc!(arch, grid, field => (x = (Dirichlet(a), Dirichlet(b)), y = Neumann()))
+```

docs/src/concepts/fields.md

@@ -0,0 +1,141 @@
+# Fields
+
+With a given grid that allows us to define each point uniquely in a high-dimensional space, we abstract the data values to be defined on the grid under the concept `AbstractField`. Following is the type tree of the abstract field and its derived data types.
+
+```@raw html
+<img src="../assets/field_type_tree.svg" width="70%"/>
+```
+
+## Defining a multi-dimensional `Field`
+
+Consider the following example, where we predefined a variable `grid` of type `Chmy.UniformGrid`, similar as in the previous section [Grids](./grids.md). We can now define physical properties on the grid.
+
+When defining a scalar field `Field` on the grid, we need to specify the arrangement of the field values. These values can either be stored at the cell centers of each control volume `Center()` or on the cell vertices/faces `Vertex()`.
+
+```julia
+# Define geometry, architecture..., a 2D grid
+grid = UniformGrid(arch; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=(nx, ny))
+
+# Define pressure as a scalar field
+Pr   = Field(backend, grid, Center())
+```
+
+With the methods `VectorField` and `TensorField`, we can conveniently construct 2-dimensional and 3-dimensional fields, with predefined locations for each field dimension on a staggered grid.
+
+```julia
+# Define velocity as a vector field on the 2D grid
+V = VectorField(backend, grid)
+
+# Define stress as a tensor field on the 2D grid
+τ = TensorField(backend, grid)
+```
+
+!!! tip "Acquiring Locations on the Grid Cell"
+    One could use a convenient getter for obtaining locations of variable on the staggered-grid. Such as `Chmy.location(Pr)` for scalar-valued pressure field and `Chmy.location(τ.xx)` for a tensor field.
+
+### Setup Initial Conditions on `Field`
+
+We could also modify field values in order to set up initial conditions on the values. We first have an empty field variable `C` defined, by using `interior(C)` we can inspect its values being empty equal to zero.
+
+```julia
+C = Field(backend, grid, Center())
+```
+
+```julia
+# Set initial values of the field randomly
+set!(C, grid, (_, _) -> rand())
+
+# Set initial values to 2D Gaussian
+set!(C, grid, (x, y) -> exp(-x^2 - y^2))
+```
+
+```@raw html
+<img src="../assets/field_set_ic_random.png" width="50%"/>
+```
+
+```@raw html
+<img src="../assets/field_set_ic_gaussian.png" width="50%"/>
+```
+
+
+
+A slightly more complex usage involves passing extra parameters to be used for initial conditions setup.
+
+```julia
+# Define a temperature field with values on cell centers
+T = Field(backend, grid, Center())
+
+# Function for setting up the initial conditions on T
+init_incl(x, y, x0, y0, r, in, out) = ifelse((x - x0)^2 + (y - y0)^2 < r^2, in, out)
+
+# Set up the initial conditions with parameters specified
+set!(T, grid, init_incl; parameters=(x0=0.0, y0=0.0, r=0.1lx, in=T0, out=Ta))
+```
+
+## Defining a parameterized `FunctionField`
+
+A field could also be represented in a parameterized way, having a function that associates a single number to every point in the space.
+
+An object of the concrete type `FunctionField` can be initialized with its constructor. The constructor takes in 
+
+
+1. A function `func`
+2. A `grid`
+3. A location tuple `loc` for specifying the distribution of variables
+
+Optionally, one can also use the boolean variable `discrete` to indicate if the function field is typed `Discrete` or `Continuous`. Any additional parameters to be used in the function `func` can be passed to the optional parameter `parameters`.
+
+### Example: Creation of a parameterized function field
+Followingly, we create a `gravity` variable that is two-dimensional and comprises of two parameterized `FunctionField` objects on a predefined uniform grid `grid`.
+
+**1. Define Functions that Parameterize the Field**
+
+In this step, we specify how the gravity field should be parameterized in x-direction and y-direction, with `η` as the additional parameter used in the parameterization.
+
+```julia
+# forcing terms
+ρgx(x, y, η) = -0.5 * η * π^2 * sin(0.5π * x) * cos(0.5π * y)
+ρgy(x, y, η) = 0.5 * η * π^2 * cos(0.5π * x) * sin(0.5π * y)
+```
+
+**2. Define Locations for Variable Positioning**
+
+We specify the location on the fully-staggered grid as introduced in the _Location on a Grid Cell_ section of the concept [Grids](./grids.md).
+
+```julia
+vx_node = (Vertex(), Center())
+vy_node = (Center(), Vertex())
+```
+
+**3. Define the 2D Gravity Field**
+
+By specifying the locations on which the parameterized field should be calculated, as well as concretizing the value `η = η0` by passing it as the optional parameter `parameters` to the constructor, we can define the 2D gravity field:
+
+```julia
+η0 = 1.0
+gravity = ( x=FunctionField(ρgx, grid, vx_node; parameters=η0),
+            y=FunctionField(ρgy, grid, vy_node; parameters=η0))
+```
+
+## Defining Constant Fields
+
+For completeness, we also provide an abstract type `ConstantField`, which comprises of a generic `ValueField` type, and two special types `ZeroField`, `OneField` for conveninent usage. With such a construct, we can easily define value fields properties and other parameters using constant values in a straightforward and readable manner. Moreover, explicit information about the grid on which the field should be defined can be abbreviated. For example:
+
+```julia
+# Defines a field with constant values 1.0
+field = Chmy.ValueField(1.0)
+```
+
+Alternatively, we could also use the conveninent and lightweight contructor.
+
+```julia
+# Defines a field with constant value 1.0 using the convenient constructor
+onefield = Chmy.OneField{Float64}()
+```
+
+Notably, these two fields shall equal to each other as expected.
+
+```julia
+julia> field == onefield
+true
+```

docs/src/concepts/grid_operators.md

@@ -0,0 +1,110 @@
+# Grid Operators
+
+The gist of the finite difference relies on replacing derivatives by difference quotients anchored on structured grids. We currently support various finite difference operators for fields defined in Cartesian coordinates. The table below summarizes the most common usage of grid operators, with the grid `g::StructuredGrid` and index `I = @index(Global, Cartesian)` defined and `P = Field(backend, grid, location)` is some field defined on the grid `g`.
+
+| Mathematical Formulation | Code |
+|:-------|:------------|
+| $\frac{\partial}{\partial x} P$ | ` ∂x(P, g, I)` |
+| $\frac{\partial}{\partial y} P$ | ` ∂y(P, g, I)` |
+| $\frac{\partial}{\partial z} P$ | ` ∂z(P, g, I)` |
+| $\nabla P$ | ` divg(P, g, I)` |
+
+## Computing the Divergence of a Vector Field
+
+To illustrate the usage of grid operators, we compute the divergence of an vector field $V$ using the `divg` function. We first allocate memory for required fields.
+
+```julia
+V  = VectorField(backend, grid)
+∇V = Field(backend, grid, Center())
+# use set! to set up the initial vector field...
+```
+
+The kernel that computes the divergence needs to have the grid information passed as for other finite difference operators.
+
+```julia
+@kernel inbounds = true function update_∇!(V, ∇V, g::StructuredGrid, O)
+    I = @index(Global, NTuple)
+    I = I + O
+    ∇V[I...] = divg(V, g, I...)
+end
+```
+
+The kernel can then be launched when required as we detailed in section [Kernels](./kernels.md).
+
+```julia
+launch(arch, grid, update_∇! => (V, ∇V, grid))
+```
+
+
+
+## Masking
+
+Masking is particularly important when performing finite differences on GPUs, as it allows for efficient and accurate computations by selectively applying operations only where needed, allowing more flexible control over the grid operators and improving performance. Thus, by providing masked grid operators, we enable more flexible control over the domain on which the grid operators should be applied for advanced users.
+
+In the following example, we first define a mask `ω` on the 2D `StructuredGrid`. Then we specify to **not** mask the center area of all Vx, Vy nodes (accessible through `ω.vc`, `ω.cv`) on the staggered grid.
+
+```julia
+# define the mask
+ω = FieldMask2D(arch, grid) # with backend and grid geometry defined
+
+# define the initial inclusion
+r = 2.0
+init_inclusion = (x,y) -> ifelse(x^2 + y^2 < r^2, 1.0, 0.0)
+
+# mask all other entries other than the initial inclusion 
+set!(ω.vc, grid, init_inclusion)
+set!(ω.cv, grid, init_inclusion)
+```
+
+We can then pass the mask to other grid operators when applying them within the kernel. When computing masked derivatives, a mask being the subtype of `AbstractMask` is premultiplied at the corresponding grid location for each operand:
+
+```julia
+@kernel function update_strain_rate!(ε̇, V, ω::AbstractMask, g::StructuredGrid, O)
+    I = @index(Global, NTuple)
+    I = I + O
+    # with masks ω
+    ε̇.xx[I...] = ∂x(V.x, ω, g, I...)
+    ε̇.yy[I...] = ∂y(V.y, ω, g, I...)
+    ε̇.xy[I...] = 0.5 * (∂y(V.x, ω, g, I...) + ∂x(V.y, ω, g, I...))
+end
+```
+
+The kernel can be launched as follows, with some launcher defined using `launch = Launcher(arch, grid)`:
+
+```julia
+# define fields
+ε̇ = TensorField(backend, grid)
+V = VectorField(backend, grid)
+
+# launch kernel
+launch(arch, grid, update_strain_rate! => (ε̇, V, ω, grid))
+```
+
+
+## Interpolation
+
+Interpolating physical parameters such as permeability and density between various types of nodes is frequently necessary on a staggered grid. Chmy.jl provides conveninent functions for performing interpolations of field values between different types of nodes.
+
+In the following example, we specify to use the linear interpolation rule `lerp` when interpolating nodal values of the density field `ρ`, defined on pressure nodes (with location `(Center(), Center())`) to `ρvx` and `ρvy`, defined on Vx and Vy nodes respectively.
+
+```julia
+# define density ρ on pressure nodes
+ρ   = Field(backend, grid, Center())
+ρ0  = 3.0; set!(ρ, ρ0)
+
+# allocate memory for density on Vx, Vy nodes
+ρvx = Field(backend, grid, (Vertex(), Center()))
+ρvy = Field(backend, grid, (Center(), Vertex()))
+```
+
+The kernel `interpolate_ρ!` performs the actual interpolation of nodal values and requires the grid information passed by `g`.
+
+```julia
+@kernel function interpolate_ρ!(ρ, ρvx, ρvy, g::StructuredGrid, O)
+    I = @index(Global,  NTuple)
+    I = I + O
+    # Interpolate from pressure nodes to Vx, Vy nodes
+    lerp(ρvx, location(ρ), g, I...)
+    lerp(ρvy, location(ρ), g, I...)
+end
+```