Improve the content of the new documentation page

docs/make.jl

@@ -29,7 +29,7 @@ makedocs(
            "Warm-start" => "warm-start.md",
            "Matrix-free operators" => "matrix_free.md",
            "Callbacks" => "callbacks.md",
-           "Advanced features" => "internal.md",
+           "Custom Workspaces" => "custom_workspace.md",
            "Performance tips" => "tips.md",
            "Tutorials" => ["CG" => "examples/cg.md",
                            "CAR" => "examples/car.md",

docs/src/internal.md → docs/src/custom_workspaces.md

@@ -1,30 +1,28 @@
-# Custom vector type for the 2D Poisson equation with halo regions
+## Custom workspaces for the Poisson equation with halo regions
 
-## Introduction
+### Introduction
 
-The 2D Poisson equation is a fundamental partial differential equation (PDE) widely used in physics and mathematics to model various phenomena, including temperature distribution and incompressible fluid flow.
-It can be expressed as:
+The Poisson equation is a fundamental partial differential equation (PDE) commonly used in physics and mathematics to model phenomena such as temperature distribution and incompressible fluid flow.
+In a 2D Cartesian domain, it can be expressed as:
 
 ```math
 \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x, y)
 ```
 
-where:
-- $u(x, y)$ is the unknown function we seek.
-- $f(x, y)$ is a known function (the "source term") representing the distribution of sources or sinks in the domain.
+Here, $u(x, y)$ represents the potential function we seek, while $f(x, y)$ is a known function (the "source term") that indicates the distribution of sources or sinks within the domain.
+This equation is typically solved over a rectangular region with boundary conditions specified at the edges.
 
-This equation is typically solved over a rectangular domain with boundary conditions specified at the edges.
+This overview will discuss the numerical methods used to solve the Poisson equation, emphasizing the application of a Laplacian operator on a specialized data structure that incorporates halo regions. 
 
-## Finite difference discretization
+### Finite difference discretization
 
-To numerically solve the Poisson equation, we employ the finite difference method to approximate the second derivatives on a grid.
-This approach involves dividing the domain into a grid of points and using differences between neighboring values to approximate derivatives.
+To numerically solve the Poisson equation, we use the finite difference method to approximate the second derivatives on a grid.
+This approach involves dividing the domain into a grid of points, where each point represents an approximation of the solution $u$ at its corresponding coordinates.
 
-Assuming a square domain $[0, L] \times [0, L]$ discretized with $(N_x+2, N_y+2)$ points along each dimension (with grid spacing $h_x = \frac{L}{N_x+1}$ and $h_y = \frac{L}{N_y+1}$, let $u_{i,j}$ denote the approximation of $u(x_i, y_j)$ at the grid point $(x_i, y_j) = (ih, jh)$.
+Assuming a square domain $[0, L] \times [0, L]$, we discretize it with $(N_x + 2, N_y + 2)$ points in each dimension, leading to grid spacings defined as $h_x = \frac{L}{N_x + 1}$ and $h_y = \frac{L}{N_y + 1}$.
+Let $u_{i,j}$ denote the approximation of $u(x_i, y_j)$ at the grid point $(x_i, y_j) = (ih, jh)$.
 
-### Discretized Laplacian
-
-The 2D Laplacian, $\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}$, can be approximated at each of the $N^2$ interior grid point $(i, j)$ using the central difference formula:
+The 2D Laplacian can be approximated at each of the $N^2$ interior grid points $(i, j)$ using the central difference formula:
 
 ```math
 \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2}
@@ -34,45 +32,46 @@ The 2D Laplacian, $\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\part
 \frac{\partial^2 u}{\partial y^2} \approx \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2}
 ```
 
-Combining these yields the discrete form of the Poisson equation:
+Combining these approximations yields the discrete form of the Poisson equation, which can be expressed as:
 
 ```math
 \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2} = f_{i,j}
 ```
 
-This represents a system of linear equations for the unknowns $u_{i,j}$ at each interior grid point.
+This results in a system of linear equations for the unknowns $u_{i,j}$ at each interior grid point.
 
 ### Boundary conditions
 
-To complete the system, we need boundary conditions for $u$ along the domain edges.
-Common options include:
-- **Dirichlet boundary conditions**: specifying the value of $u$ on the boundary.
-- **Neumann boundary conditions**: specifying the derivative (flux) of $u$ normal to the boundary.
+To complete the system, we need to define boundary conditions for $u$ along the edges of the domain. Common options include:
+
+- **Dirichlet boundary conditions**: These specify the value of $u$ directly on the boundary.
+- **Neumann boundary conditions**: These specify the derivative (or flux) of $u$ normal to the boundary.
 
-## Implementing halo regions with `MyVector`
+### Implementing halo regions with `MyVector`
 
 In practical applications, particularly in parallel computing, it is common to introduce **halo regions** (or ghost cells) around the grid.
-These additional layers store boundary values from neighboring subdomains, allowing each subdomain to compute stencils near its boundaries independently without immediate communication with neighboring domains.
-Halo regions simplify boundary condition management in distributed or multi-threaded environments.
+These additional layers store boundary values from neighboring subdomains, allowing each subdomain to compute stencils near its boundaries independently, without the need for immediate communication with neighboring domains. Halo regions thus simplify boundary condition management in distributed or multi-threaded environments.
 
-In Krylov.jl, while internal storage for each Krylov method typically expects an `AbstractVector`, specific applications can benefit from structured data layouts.
-This is where a specialized vector type called **`MyVector`** comes into play. 
+In the context of Krylov.jl, while internal storage for each Krylov method typically expects an `AbstractVector`, specific applications can benefit from structured data layouts.
+This is where a specialized vector type called **`MyVector`** becomes advantageous.
 
-Using `MyVector` with halo regions enables the implementation of finite difference stencils without boundary condition checks, enhancing both readability and performance.
-The **`OffsetArray`** type from the [OffsetArrays.jl](https://github.com/JuliaArrays/OffsetArrays.jl) package supports custom indexing, making it ideal for grids with halo regions.
-By wrapping an `OffsetArray` within an `MyVector`, we can access elements using custom offsets that align with the grid's physical layout.
-This configuration allows for ``if-less'' stencils, avoiding direct boundary condition checks within the core loop, resulting in cleaner and potentially faster code.
+By using `MyVector` with halo regions, we can implement finite difference stencils without the overhead of boundary condition checks.
+This not only enhances code readability but also improves performance.
+The type **`OffsetArray`** from the package [OffsetArrays.jl](https://github.com/JuliaArrays/OffsetArrays.jl) supports custom indexing, making it ideal for grids that include halo regions.
+By wrapping an `OffsetArray` within a `MyVector`, we can access elements using custom offsets that align with the physical layout of the grid.
+This configuration enables **"if-less"** stencils, effectively avoiding direct boundary condition checks within the core loop, resulting in cleaner and potentially faster code.
 
-The design of `MyVector` can be easily adapted for 1D, 2D, or 3D problems with minimal changes, providing flexibility in handling various grid configurations.
+Moreover, the design of `MyVector` is flexible and can be easily adapted for 1D, 2D, or 3D problems with minimal changes, providing versatility in handling various grid configurations.
 
-## Definition of the `MyVector`
+### Definition of the `MyVector`
 
 The `MyVector` type is a specialized vector designed for efficient handling of grid-based computations, particularly in the context of finite difference methods with halo regions.
 It is parameterized by:
+
 - **`FC`**: The element type of the vector.
-- **`D`**: The type of the data array, which utilizes `OffsetArray` to enable custom indexing.
+- **`D`**: The type of the data array, which utilizes `OffsetArray` to facilitate custom indexing.
 
-Here is the definition of the `MyVector`:
+Below is the definition of `MyVector`:
 
 ```julia
 using OffsetArrays
@@ -86,11 +85,10 @@ struct MyVector{FC, D} <: AbstractVector{FC}
     end
 end
 
-# Constructor
 function MyVector{FC,D}(::UndefInitializer, l::Int64) where {FC,D}
     m = n = sqrt(l) |> Int
-    data = zeros(FC, m+2, n+2)
-    v = OffsetMatrix(data, 0:m+1, 0:n+1)
+    data = zeros(FC, m + 2, n + 2)
+    v = OffsetMatrix(data, 0:m + 1, 0:n + 1)
     return MyVector(v)
 end
 
@@ -107,21 +105,18 @@ end
 
 function Base.getindex(v::MyVector, idx)
     m, n = size(v.data)
-    row = div(idx-1, n-2) + 1
-    col = mod(idx-1, n-2) + 1
+    row = div(idx - 1, n - 2) + 1
+    col = mod(idx - 1, n - 2) + 1
     return v.data[row, col]
 end
 ```
 
-The functions `size` and `getindex` are defined to enable display in the REPL.
-
-## Using `MyVector` for the 2D Poisson equation
-
-By utilizing `MyVector`, we can implement the finite difference stencil for the Laplacian operator efficiently, eliminating the need for conditional checks for boundary elements.
+The `size` and `getindex` functions are defined to facilitate display in the REPL, allowing for easy interaction with `MyVector`, although they are not strictly necessary for the functionality of Krylov.jl.
 
 ### Stencil implementation
 
-Assuming `data` is initialized as an `OffsetArray` with appropriate halo regions, we can define a matrix-free Laplacian operator and apply a typical 5-point stencil operation as follows:
+To efficiently apply the discrete Laplacian operator using a matrix-free approach and a typical 5-point stencil operation, we leverage `OffsetArray` for handling halo regions.
+This allows seamless access to boundary values, enabling clear and performant computation of the Laplacian without the need for direct boundary condition checks within the core stencil loop.
 
 ```julia
 # Define a matrix-free Laplacian operator
@@ -153,23 +148,10 @@ function LinearAlgebra.mul!(y::MyVector{Float64}, A::LaplacianOperator, u::MyVec
 end
 ```
 
-### Benefits of Using `MyVector`
-
-Utilizing `MyVector` offers several significant benefits for solving the 2D Poisson equation:
+### Required methods for Krylov.jl compatibility
 
-- **Simplified Code**: The custom indexing capabilities of `OffsetArray` allow for the incorporation of boundary data from halo regions. This eliminates the need for boundary checks within the core loop, resulting in clearer and more maintainable code.
-
-- **Performance**: By removing boundary checks, we reduce branching in the code, which enhances computational efficiency, particularly for large grids.
-This leads to faster execution times and better overall performance.
-
-- **Flexibility**: `MyVector` can be easily extended to accommodate more complex stencils or additional dimensions (e.g., 3D grids) by simply adjusting the offsets.
-This adaptability makes it a powerful tool for various numerical applications.
-
-By leveraging these advantages, we can efficiently solve the 2D Poisson equation while maintaining a clear and concise code structure.
-
-## Required methods for Krylov.jl compatibility
-
-To integrate `MyVector` with Krylov.jl, the following operations must be defined for compatibility:
+To use `MyVector` within Krylov.jl, we must define several essential vector operations such as dot products, norms, scalar multiplication, and element-wise updates.
+Below, we present the required methods that ensure `MyVector` is compatible with Krylov.jl:
 
 ```julia
 using Krylov
@@ -259,9 +241,9 @@ function Krylov.kfill!(x::MyVector{T}, val::T) where T <: FloatOrComplex
 end
 ```
 
-These methods enable Krylov.jl to use custom vector types, allowing for operations like dot products, norms, scalar multiplication, and element-wise updates, which are essential for Krylov solvers.
+These methods enable Krylov.jl to utilize custom vector types, enhancing the flexibility and performance of Krylov solvers.
 
-## Complete example
+### Solve the 2D poisson equation
 
 ```julia
 using Krylov, LinearAlgebra, OffsetArrays
@@ -299,10 +281,15 @@ u_star = [sin(π * i * Δx) * sin(π * j * Δy) for i=1:Nx, j=1:Ny]
 norm(u_sol.data[1:Nx, 1:Ny] - u_star, Inf)
 ```
 
-## Conclusion
+### Conclusion
+
+The implementation of a 2D Poisson equation solver using `MyVector` offers several notable advantages that enhance both code clarity and efficiency.
+With the custom indexing capabilities of `OffsetArray`, we can seamlessly incorporate boundary data from halo regions, simplifying the code by eliminating boundary checks within the core loop.
+This results in a more readable and maintainable code structure.
 
-This overview illustrates the effective implementation of a 2D Poisson equation solver utilizing a specialized vector type storage with `MyVector` to accommodate halo regions.
-By leveraging custom indexing, we significantly enhance both code readability and performance, enabling a flexible framework that can be adapted to a variety of applications.
+Moreover, by removing boundary checks, we reduce branching in the code, which improves computational efficiency, especially for large grids.
+This approach leads to faster execution times while preserving optimal performance.
+The flexibility of `MyVector` also makes it easy to extend for more complex stencils or additional dimensions, such as 3D grids.
 
 !!! info
     The package [Oceanigans.jl](https://github.com/CliMA/Oceananigans.jl) utilizes a similar approach with its type `Field` to solve linear systems involving millions of variables efficiently with Krylov.jl.