Automatically wrapping with Clang.jl

.gitignore

@@ -1,10 +1,8 @@
-deps/ComplexDouble/
-deps/RealDouble/
-deps/RealSingle/
-deps/deps.jl
-deps/petsc-*.tar.gz
-*.jl.cov
-*.jl.mem
+*.swp
+*.DS_Store
+*~
 docs/build/
 docs/site/
-Manifest.toml
+Manifest.toml
+*.vscode/
+examples/viz2D_out/

Project.toml

@@ -4,16 +4,21 @@ authors = ["Simon Byrne <[email protected]>"]
 version = "0.1.3"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+MPICH_jll = "7cb0a576-ebde-5e09-9194-50597f1243b4"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PETSc_jll = "8fa3689e-f0b9-5420-9873-adf6ccf46f2d"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 
 [compat]
-MPI = "0.15, 0.16, 0.17, 0.18, 0.19"
-PETSc_jll = "3.15.2"
-julia = "1.3"
+MPI = "0.20"
+MPICH_jll = "4.1.2"
+PETSc_jll = "3.21"
+julia = "1.9"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

README.md

@@ -3,8 +3,22 @@
 [![Build Status](https://github.com/JuliaParallel/PETSc.jl/workflows/CI/badge.svg)](https://github.com/JuliaParallel/PETSc.jl/actions/workflows/ci.yml)
 [![doc dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliaparallel.github.io/PETSc.jl/dev/)
 
+PETSc.jl provides an interface to the [PETSc](https://www.mcs.anl.gov/petsc/) library,
+allowing the combination of Julia features (such as automatic differentiation) with the PETSc's infrastructure, including
+linear, nonlinear, and optimization solvers, timesteppers, domain management (DM), and more, in a distributed-memory (MPI) environment.
+
+This package comprises two main components:
+
+1. An automatically generated, low-level interface for large parts of the PETSc API.
+2. A curated, high-level, more Julianic interface for selected functionality.
+
+The low-level interface covers the entire PETSc API, but may be awkward to work with and likely requires
+previous experience with PETSc to use effectively. The high level interface is designed to be
+more familiar and convenient for Julia users but exposes only a small portion of the functionality
+of the underlying library. This high-level interface should be considered unstable, as its
+implementation involves design decisions which may be revisited - note the low version number
+of this package if relying on these high-level features.
 
-This package provides a low level interface for [PETSc](https://www.mcs.anl.gov/petsc/) and allows combining julia features (such as automatic differentiation) with the PETSc infrastructure and nonlinear solvers.
 
 ## Installation
 

docs/src/index.md

@@ -1,7 +1,51 @@
 # PETSc.jl
 
- [PETSc.jl](https://github.com/JuliaParallel/PETSc.jl) is a Julia wrapper for the Portable, Extensible Toolkit for Scientific Computation [PETSc](https://petsc.org/release/documentation/manual/) package, which allows solving ordinary and partial differential equations in parallel on laptops or massively parallel high-performance systems.
+ [PETSc.jl](https://github.com/JuliaParallel/PETSc.jl) is a Julia wrapper for the Portable, Extensible Toolkit for Scientific Computation [PETSc](https://petsc.org/) package, which allows solving ordinary and partial differential equations in parallel on laptops or massively parallel high-performance systems.
 
  The use of Julia greatly simplifies the code that developers have to write, while allowing to employ Julia features such as automatic differentiation. The Julia wrapper also comes with a pre-built library, which greatly simplifies the process of getting your first code working in parallel, on different operating systems. In many cases, the Julia code is significantly shorter than its C counterpart.
 
- This wrapper mimics the PETSc-functionality as closely as possible, but remains work in progress (meaning that not everything has been translated yet). See the official [user guide](https://petsc.org/release/overview/) if you want to learn more about PETSc in general. For Julia-specific examples, have a look at our [examples](https://github.com/JuliaParallel/PETSc.jl/tree/main/examples) or [tests](https://github.com/JuliaParallel/PETSc.jl/tree/main/test). 
+ This wrapper mimics the PETSc-functionality as closely as possible, but remains work in progress (meaning that not everything has been translated yet). See the official [user guide](https://petsc.org/release/overview/) if you want to learn more about PETSc in general. For Julia-specific examples, have a look at our [examples](https://github.com/JuliaParallel/PETSc.jl/tree/main/examples) or [tests](https://github.com/JuliaParallel/PETSc.jl/tree/main/test).
+
+This package includes a low-level, automatically-generated wrapper layer, upon which a higher-level interface is built.
+
+
+## The High-Level Interface
+
+The high level interface is designed to be familiar and convenient for Julia users, but exposes only a small portion of the functionality
+of the underlying PETSc library.  This interface should be considered *unstable*, as its implementation involves design decisions which are likely to be revisited - note the low version number
+of this package if considering relying on these features.
+
+For example, with this interface, PETSc's [KSP](https://petsc.org/release/docs/manual/ksp) linear solvers (including Krylov methods) can be used in a way similar to solvers from other Julia packages. See the example in [Getting started](@ref) and the API in [KSP](@ref).
+
+## The Low-Level Interface
+
+The low-level interface covers more of the PETSc API, but may be awkward to work with and likely requires
+previous experience with PETSc to use effectively. It is automatically generated with [Clang.jl](https://github.com/JuliaInterop/Clang.jl).
+
+The high-level interface described in [KSP](@ref) creates a [KSP](https://petsc.org/release/docs/manual/ksp) linear solver object via the low-level interface to [`KSPSolve()`](https://petsc.org/release/docs/manualpages/KSP/KSPCreate.html), with the use of constructs which require knowledge of PETSc's nature as a [C library](https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/). Expert users are of course free to directly use the low level interface, as in this simple example which directly calls [`PetscGetVersionNumber()`](https://petsc.org/release/docs/manualpages/PetscGetVersionNumber.html).
+
+```@example
+using MPI
+MPI.Initialized() || MPI.Init()
+using PETSc
+
+petsclib = PETSc.petsclibs[1]
+PetscInt = petsclib.PetscInt
+PetscErrorCode = Cint
+
+PETSc.initialize(petsclib)
+major = Ref{PetscInt}(0)
+minor = Ref{PetscInt}(0)
+subminor = Ref{PetscInt}(0)
+release = Ref{PetscInt}(0)
+error_code = ccall(
+      (:PetscGetVersionNumber, petsclib.petsc_library),
+      PetscErrorCode,
+      (Ptr{PetscInt}, Ptr{PetscInt}, Ptr{PetscInt}, Ptr{PetscInt}),
+      major, minor, subminor, release
+      )
+
+version = (major[], minor[], subminor[])
+println("PETSc $(version[1]).$(version[2]).$(version[3])")
+PETSc.finalize(petsclib)
+```

examples/Liouville_Bratu_Gelfand.jl

@@ -0,0 +1,276 @@
+# INCLUDE IN MPI TEST
+#=
+In this example we solve the [Liouville–Bratu–Gelfand
+equation](https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation):
+```math
+∇² ψ + λ exp(ψ)
+```
+with zero Dirichlet boundary conditions.
+
+To solve the problem we use the standard central, finite difference
+approximation of the Laplacian in `d`-dimensions
+
+This example is motivated by the following PETSc examples:
+- [`snes/ex5`](https://petsc.org/release/src/snes/tutorials/ex5.c.html)
+- [`snes/ex14`](https://petsc.org/release/src/snes/tutorials/ex14.c.html)
+=#
+
+using MPI
+using PETSc
+using OffsetArrays: OffsetArray
+using LinearAlgebra: norm
+
+opts = if !isinteractive()
+    PETSc.parse_options(ARGS)
+else
+    (ksp_monitor = true, ksp_view = true)
+end
+
+# Set our MPI communicator
+comm = MPI.COMM_WORLD
+
+# Set our PETSc Scalar Type
+PetscScalar = Float64
+
+# get the PETSc lib with our chosen `PetscScalar` type
+petsclib = PETSc.getlib(; PetscScalar = PetscScalar)
+
+# Initialize PETSc
+PETSc.initialize(petsclib)
+
+# dimensionality of the problem
+dim = haskey(opts, :dim) ? opts.dim : 3
+
+# Set the total number of grid points in each direction
+Nq = ntuple(_ -> 10, dim)
+
+# Set the boundary conditions on each side
+bcs = ntuple(_ -> PETSc.DM_BOUNDARY_NONE, dim)
+
+# Set parameter
+λ = PetscScalar(6)
+
+# Create the PETSC dmda object
+da = PETSc.DMDA(
+    petsclib,
+    comm,
+    bcs,                     # boundary conditions
+    Nq,                      # Global grid size
+    1,                       # Number of DOF per node
+    1,                       # Stencil width
+    PETSc.DMDA_STENCIL_STAR; # Stencil type
+    opts...,
+)
+
+# Create the PETSC snes object
+snes = PETSc.SNES(petsclib, comm; opts...)
+
+# add the da to the snes
+PETSc.setDM!(snes, da)
+
+# Set up the initial guess
+x = PETSc.DMGlobalVec(da)
+xl = PETSc.DMLocalVec(da)
+PETSc.withlocalarray!(xl; read = false) do l_x
+    corners = PETSc.getcorners(da)
+
+    # Get the global grid dimensions
+    Nq = PETSc.getinfo(da).global_size
+
+    # Figure out the interior points 
+    int_min = min(CartesianIndex(corners.size), CartesianIndex(2, 2, 2))
+    int_max = max(CartesianIndex(corners.size .- 1), CartesianIndex(1, 1, 1))
+    interior = (int_min):(int_max)
+
+    # Allows us to adress the local array with global indexing
+    ox = @view PETSc.reshapelocalarray(l_x, da)[1, :, :, :]
+
+    # Set up the global coordinates in each direction
+    # -1 to 1 when Nq > 1 and 0 otherwise
+    coords = map(
+        Nq ->
+            Nq == 1 ? range(PetscScalar(0), stop = 0, length = Nq) :
+            range(-PetscScalar(1), stop = 1, length = Nq),
+        Nq,
+    )
+
+    scaling = λ / (λ + 1)
+
+    # Loop over all the points on the processor and set the initial condition to
+    # be a hat function
+    for i in ((corners.lower):(corners.upper))
+        ox[i] =
+            scaling * sqrt(
+                min(
+                    (1 - abs(coords[1][i[1]])),
+                    (1 - abs(coords[2][i[2]])),
+                    (1 - abs(coords[3][i[3]])),
+                ),
+            )
+    end
+end
+PETSc.update!(x, xl, PETSc.INSERT_VALUES)    # update local -> global
+
+# Set up the nonlinear function
+r = similar(x)
+PETSc.setfunction!(snes, r) do g_fx, snes, g_x
+    # Get the DMDA associated with the snes
+    da = PETSc.getDMDA(snes)
+
+    # Get a local vector and transfer the data from the global vector into it
+    l_x = PETSc.DMLocalVec(da)
+    PETSc.update!(l_x, g_x, PETSc.INSERT_VALUES)
+
+    ghostcorners = PETSc.getghostcorners(da)
+    corners = PETSc.getcorners(da)
+
+    # Global grid size
+    Nq = PETSc.getinfo(da).global_size
+
+    # grid spacing in each dimension
+    Δx, Δy, Δz = PetscScalar(1) ./ Nq
+
+    # Get local arrays
+    PETSc.withlocalarray!(
+        (g_fx, l_x);
+        read = (false, true),
+        write = (true, false),
+    ) do fx, x
+
+        # reshape the array and allow for global indexing
+        x = @view PETSc.reshapelocalarray(x, da)[1, :, :, :]
+        fx = @view PETSc.reshapelocalarray(fx, da)[1, :, :, :]
+
+        # Store a tuple of stencils in each direction
+        stencils = (
+            (
+                CartesianIndex(-1, 0, 0),
+                CartesianIndex(0, 0, 0),
+                CartesianIndex(1, 0, 0),
+            ),
+            (
+                CartesianIndex(0, -1, 0),
+                CartesianIndex(0, 0, 0),
+                CartesianIndex(0, 1, 0),
+            ),
+            (
+                CartesianIndex(0, 0, -1),
+                CartesianIndex(0, 0, 0),
+                CartesianIndex(0, 0, 1),
+            ),
+        )
+        # Weights for each direction
+        weights = (Δy * Δz / Δx, Δx * Δz / Δy, Δx * Δy / Δz)
+
+        # loop over indices and set the function value
+        for ind in ((corners.lower):(corners.upper))
+            # If on the boundary just set equal to the incoming data
+            # otherwise apply the finite difference operator
+            if any(ntuple(j -> ind[j] == 1 || ind[j] == Nq[j], dim))
+                fx[ind] = x[ind]
+            else
+                # Apply the source
+                u = -Δx * Δy * Δz * λ * exp(x[ind])
+
+                # Apply the finite diffference stencil
+                for (s, w) in zip(stencils[1:dim], weights[1:dim])
+                    u +=
+                        w * (-x[ind + s[1]] + 2 * x[ind + s[2]] - x[ind + s[3]])
+                end
+                fx[ind] = u
+            end
+        end
+    end
+
+    # Clean up the local vector
+    PETSc.destroy(l_x)
+    return 0
+end
+
+J = PETSc.MatAIJ(da)
+PETSc.setjacobian!(snes, J) do J, snes, g_x
+    # Get the DMDA associated with the snes
+    da = PETSc.getDMDA(snes)
+
+    # Get the corners of the points we own
+    corners = PETSc.getcorners(da)
+
+    # Global grid size
+    Nq = PETSc.getinfo(da).global_size
+
+    # grid spacing in each dimension
+    Δx, Δy, Δz = PetscScalar(1) ./ Nq
+
+    # Store a tuple of stencils in each direction
+    stencils = (
+        (
+            CartesianIndex(-1, 0, 0),
+            CartesianIndex(0, 0, 0),
+            CartesianIndex(1, 0, 0),
+        ),
+        (
+            CartesianIndex(0, -1, 0),
+            CartesianIndex(0, 0, 0),
+            CartesianIndex(0, 1, 0),
+        ),
+        (
+            CartesianIndex(0, 0, -1),
+            CartesianIndex(0, 0, 0),
+            CartesianIndex(0, 0, 1),
+        ),
+    )
+    # Weights for each direction
+    weights = (Δy * Δz / Δx, Δx * Δz / Δy, Δx * Δy / Δz)
+
+    # Get a local array of the solution vector
+    PETSc.withlocalarray!(g_x; write = false) do l_x
+        # reshape so we can use multi-D indexing
+        x = @view PETSc.reshapelocalarray(l_x, da)[1, :, :, :]
+
+        # loop over indices and set the function value
+        for ind in ((corners.lower):(corners.upper))
+            # If on the boundary just set equal to the incoming data
+            # otherwise apply the finite difference operator
+            if any(ntuple(j -> ind[j] == 1 || ind[j] == Nq[j], dim))
+                J[ind, ind] = 1
+            else
+                # We accumulate the diagonal and add it at the end
+                Jii = -Δx * Δy * Δz * λ * exp(x[ind]) # Apply the source
+
+                # Apply the finite diffference stencil
+                for (s, w) in zip(stencils[1:dim], weights[1:dim])
+                    Jii += w * 2
+                    J[ind, ind + s[1]] = -w
+                    J[ind, ind + s[3]] = -w
+                end
+                J[ind, ind] = Jii
+            end
+        end
+    end
+
+    # Assemble the Jacobian matrix
+    PETSc.assemble!(J)
+    return 0
+end
+
+if MPI.Comm_rank(comm) == 0
+    println(@elapsed(PETSc.solve!(x, snes)))
+else
+    PETSc.solve!(x, snes)
+end
+g = similar(x)
+snes.f!(g, snes, x)
+nm = norm(g)
+if MPI.Comm_rank(comm) == 0
+    @show nm
+end
+
+# Do some clean up
+PETSc.destroy(J)
+PETSc.destroy(x)
+PETSc.destroy(g)
+PETSc.destroy(r)
+PETSc.destroy(da)
+PETSc.destroy(snes)
+
+PETSc.finalize(petsclib)