[Documentation] Improving documentation (#24)

This PR adds the first version of documentation for Chmy.jl

Project.toml

@@ -1,7 +1,7 @@
 name = "Chmy"
 uuid = "33a72cf0-4690-46d7-b987-06506c2248b9"
 authors = ["Ivan Utkin <[email protected]>, Ludovic Raess <[email protected]>, and contributors"]
-version = "0.1.16"
+version = "0.1.17"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

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

@@ -8,14 +8,26 @@ makedocs(
     authors="Ivan Utkin, Ludovic Räss and contributors",
     format = Documenter.HTML(
         prettyurls=get(ENV, "CI", nothing) == "true", # easier local build
-        ansicolor=true
+        ansicolor=true,
+        assets = ["assets/favicon.ico"],
         ),
     modules = [Chmy],
     warnonly = [:missing_docs],
     pages = Any[
         "Home" => "index.md",
-        "Usage" => Any["usage/runtests.md"],
-        "Library" => Any["lib/modules.md"]
+        "Getting Started with Chmy.jl" => "getting_started.md",
+        "Concepts" => Any["concepts/architectures.md",
+                        "concepts/grids.md",
+                        "concepts/fields.md",
+                        "concepts/bc.md",
+                        "concepts/grid_operators.md",
+                        "concepts/kernels.md"
+        ],
+        "Examples" => Any["examples/overview.md"
+        ],
+        "Library" => Any["lib/modules.md"],
+        "Developer documentation" => Any["developer_documentation/running_tests.md",
+                                         "developer_documentation/workers.md"],
     ]
 )
 

docs/src/concepts/architectures.md

@@ -0,0 +1,42 @@
+# Architectures
+
+## Backend Selection & Architecture Initialization
+
+Chmy.jl supports CPUs, 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 objects such as grids and kernel launchers.
+
+## Specifying the device ID and stream priority
+
+On systems with multiple GPUs, passing the keyword argument `device_id` to the `Arch` constructor will select and set the selected device as a current device.
+
+For advanced users, we provide a function `activate!(arch; priority)` 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.
+
+## 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.
+
+To make the `Architecture` object aware of MPI topology, user can pass an MPI communicator object and dimensions of the Cartesian topology to the `Arch` constructor:
+
+```julia
+using MPI
+
+arch = Arch(CPU(), MPI.COMM_WORLD, (0, 0, 0))
+```
+
+Passing zeros as the last argument will automatically spread the dimensions to be as close as possible to each other, see [MPI.jl documentation](https://juliaparallel.org/MPI.jl/stable/reference/topology/#MPI.Dims_create) for details.

docs/src/concepts/bc.md

@@ -0,0 +1,67 @@
+# 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. Additional initial and/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 encounters. In the following, 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_{\boldsymbol{n}} u = g$ on $\partial \Omega$, where $\boldsymbol{n}$ 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 through the boundary |
+| 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 with `bc!()`
+
+In the following, 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 on $\partial \Omega$ | `bc!(arch, grid, field => Dirichlet())` | `bc!(arch, grid, field => Dirichlet(v))` |
+| Neumann on $\partial \Omega$ | `bc!(arch, grid, field => Neumann())` | `bc!(arch, grid, field => Neumann(v))` |
+
+Note that the syntax shown in the table above is a **fused expression** of both _specifying_ and _applying_ the boundary conditions.
+
+!!! warning "$\partial \Omega$ Refers to the Entire Domain Boundary!"
+    By specifying `field` to a single boundary condition, we impose the boundary condition on the entire domain boundary by default. See the section for "Mixed Boundary Conditions" below for specifying different BC on different parts of the domain boundary.
+
+Alternatively, one could also define the boundary conditions beforehand using `batch()` provided the `grid` information as well as the `field` variable. This way the boundary condition to be prescibed is **precomputed**.
+
+```julia
+# pre-compute batch
+bt = batch(grid, field => Neumann()) # specify Neumann BC for the variable `field`
+bc!(arch, grid, bt)                  # apply the boundary condition
+```
+
+In the script [batcher.jl](https://github.com/PTsolvers/Chmy.jl/blob/main/examples/batcher.jl), we provide a MWE using both **fused** and **precomputed** expressions for BC update.
+
+## Specifying BC within a `launch`
+
+When using `launch` to specify the execution of a kernel (more see section [Kernels](./kernels.md)), one can pass the specified boundary condition(s) as an optional parameter using `batch`, provided the grid information of the discretized space. This way we can gain efficiency from making good use of already cached values.
+
+In the 2D diffusion example as introduced in the tutorial ["Getting Started with Chmy.jl"](../getting_started.md), we need to update the temperature field `C` at k-th iteration using the values of heat flux `q` and physical time step size `Δt` from (k-1)-th iteration. When launching the kernel `update_C!` with `launch`, we simultaneously launch the kernel for the BC update using:
+
+```julia
+launch(arch, grid, update_C! => (C, q, Δt, grid); bc=batch(grid, C => Neumann(); exchange=C))
+```
+
+### 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 = Neumann(), y = (Dirichlet(b), Dirichlet(a))))
+```

docs/src/concepts/fields.md

@@ -0,0 +1,138 @@
+# 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 defined 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 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)
+```
+
+Use the function `location` to get the location of the field as a tuple. Vector and tensor fields are currently defined as `NamedTuple`'s (likely to change in the future), so one could query the locations of individual components, e.g. `location(V.x)` or `location(τ.xy)`
+
+!!! 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.
+
+### Initialising `Field`
+
+Chmy.jl provides functionality to set the values of the fields as a function of spatial coordinates:
+
+```julia
+C = Field(backend, grid, Center())
+
+# 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` allowing dispatch for special casess. 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 `OneField` type, providing type information about the contents of the field.
+
+```julia
+# Defines a field with constant value 1.0
+onefield = Chmy.OneField{Float64}()
+```
+
+Notably, these two fields shall equal to each other as expected.
+
+```julia
+julia> field == onefield
+true
+```