Skip to content

Commit

Permalink
Merge pull request #8 from bmxam/dev
Browse files Browse the repository at this point in the history
Redesign API + use finalizers
  • Loading branch information
bmxam authored Apr 23, 2023
2 parents 52f39fd + b224ba7 commit 6adb596
Show file tree
Hide file tree
Showing 30 changed files with 2,604 additions and 1,469 deletions.
25 changes: 25 additions & 0 deletions .JuliaFormatter.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
# https://github.com/domluna/JuliaFormatter.jl

# Annotates fields in a type definitions with ::Any if no type annotation is provided
annotate_untyped_fields_with_any = false

# Transforms a short function definition to a long function definition if the short
# function definition exceeds the maximum margin.
short_to_long_function_def = true

# Transforms a long function definition to a short function definition if the short
# function definition does not exceed the maximum margin.
long_to_short_function_def = true

# Format code docstrings with the same options used for the code source.
format_docstrings = true

# If true, Markdown files are also formatted. Julia code blocks will be formatted
# in addition to the Markdown being normalized.
format_markdown = true

# If true, whitespace is added for binary operations in indices
whitespace_ops_in_indices = true

# If true, superfluous newlines will be removed
remove_extra_newlines = true
6 changes: 3 additions & 3 deletions .github/workflows/Documentation.yml
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ on:
push:
branches:
- master
tags: '*'
tags: "*"
pull_request:
workflow_dispatch:

Expand All @@ -15,12 +15,12 @@ jobs:
- uses: actions/checkout@v2
- uses: julia-actions/setup-julia@latest
with:
version: '1.5'
version: "1.8"
- name: Install dependencies
run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
- name: Build and deploy
env:
GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
DOC_DEPLOYMENT: true # Workaround to avoid finding PETSc location
run: julia --project=docs/ docs/make.jl
run: julia --project=docs/ docs/make.jl
7 changes: 4 additions & 3 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
name = "PetscWrap"
uuid = "5be22e1c-01b5-4697-96eb-ef9ccdc854b8"
authors = ["bmxam <[email protected]>"]
version = "0.1.5"
version = "0.2.0"

[deps]
Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[compat]
MPI = "0.16, 0.17, 0.18, 0.19"
julia = "1.5"
MPI = "0.20"
julia = "1.5, 1.6, 1.7, 1.8"
99 changes: 58 additions & 41 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,35 +6,53 @@
PetscWrap.jl is a parallel Julia wrapper for the (awesome) [PETSc](https://www.mcs.anl.gov/petsc/) library. It can be considered as a fork from the [GridapPetsc.jl](https://github.com/gridap/GridapPETSc.jl) and [Petsc.jl](https://github.com/JuliaParallel/PETSc.jl) projects : these two projects have extensively inspired this project, and some code has even been directly copied.

The main differences with the two aformentionned projects are:

- parallel support : you can solve linear systems on multiple core with `mpirun -n 4 julia foo.jl`;
- no dependance to other Julia packages except `MPI.jl`;
- possibility to switch from one PETSc "arch" to another;
- less PETSc API functions wrappers for now.

Note that the primary objective of this project is to enable the wrapper of the SLEPc library through the [SlepcWrap.jl](https://github.com/bmxam/SlepcWrap.jl) project.

This project is only a wrapper to PETSc functions, the purpose is not to deliver a julia `Array` (it maybe be one day the purpose of a package `PetscArrays.jl`).

## How to install it

You must have installed the PETSc library on your computer and set the two following environment variables : `PETSC_DIR` and `PETSC_ARCH`.

At run time, PetscWrap.jl looks for the `libpetsc.so` using these environment variables and "load" the library.

To install the package, use the Julia package manager:

```Julia
pkg> add PetscWrap
```

## Contribute

Any contribution(s) and/or remark(s) are welcome! If you need a function that is not wrapped yet but you don't think you are capable of contributing, post an issue with a minimum working example.

Conventions to be applied in future versions ("fancy" stuff is not concerned):

- all PETSc types should have the exact same name in Julia;
- all PETSc functions should have the exact same name in julia, but without the type as a prefix, and with a lower case for the first letter. `VecSetValues` becomes `setValues`. This rule is not applied when the name conflicts with a name from `Base` (for instance `VecView` becomes `vecView` and not `view`);
- all PETSc functions must have the same number of arguments and, if possible the same names in julia, except for out-of-place arguments.
- functions arguments must all be typed. Additional functions, without type or with fewer args, can be defined if the original version is present.

## PETSc compat.
This version of PetscWrap.jl has been tested with petsc-3.13. Complex numbers are supported.

This version of PetscWrap.jl has been tested with petsc-3.19. Complex numbers are supported.

## How to use it
PETSc methods wrappers share the same name as their C equivalent : for instance `MatCreate` or `MatSetValue`. Furthermore, an optional "higher level" API, referred to as "fancy", is exposed : for instance `create_matrix` or `A[i,j] = v`). Note that this second way of manipulating PETSc will evolve according the package's author needs while the first one will try to follow PETSc official API.

PETSc methods wrappers share the almost same name as their C equivalent (with the type) : for instance `MatSetValues` becomes `setValues`. Furthermore, an optional "higher level" API, referred to as "fancy", is exposed : for instance `create_matrix` or `A[i,j] = v`. Note that this second way of manipulating PETSc will evolve according the package's author needs while the first one will try to follow PETSc official API.

You will find examples of use by building the documentation: `julia PetscWrap.jl/docs/make.jl`. Here is one of the examples:

### A first demo

This example serves as a test since this project doesn't have a "testing" procedure yet. In this example,
the equation ``u'(x) = 2`` with ``u(0) = 0`` is solved on the domain ``[0,1]`` using a backward finite
the equation `u'(x) = 2` with `u(0) = 0` is solved on the domain `[0,1]` using a backward finite
difference scheme.

In this example, PETSc classic method names are used. For more fancy names, check the fancy version.
Expand Down Expand Up @@ -63,36 +81,36 @@ Number of mesh points and mesh step

```julia
n = 11
Δx = 1. / (n - 1)
Δx = 1.0 / (n - 1)
```

Create a matrix and a vector
Create a matrix and a vector (you can specify the MPI communicator if you want)

```julia
A = MatCreate()
b = VecCreate()
A = create(Mat)
b = create(Vec)
```

Set the size of the different objects, leaving PETSC to decide how to distribute. Note that we should
set the number of preallocated non-zeros to increase performance.

```julia
MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)
VecSetSizes(b, PETSC_DECIDE, n)
setSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)
setSizes(b, PETSC_DECIDE, n)
```

We can then use command-line options to set our matrix/vectors.

```julia
MatSetFromOptions(A)
VecSetFromOptions(b)
setFromOptions(A)
setFromOptions(b)
```

Finish the set up

```julia
MatSetUp(A)
VecSetUp(b)
setUp(A)
setUp(b)
```

Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor.
Expand All @@ -101,88 +119,87 @@ As in PETSc, the `rstart, rend = *GetOwnershipRange` methods indicate the first
the fancy version of this example for a more Julia-like convention.

```julia
b_start, b_end = VecGetOwnershipRange(b)
b_start, b_end = getOwnershipRange(b)
```

Now let's build the right hand side vector. Their are various ways to do this, this is just one.

```julia
n_loc = VecGetLocalSize(b) # Note that n_loc = b_end - b_start...
VecSetValues(b, collect(b_start:b_end-1), 2 * ones(n_loc))
n_loc = getLocalSize(b) # Note that n_loc = b_end - b_start...
setValues(b, collect(b_start:(b_end - 1)), 2 * ones(n_loc))
```

And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices.

```julia
A_start, A_end = MatGetOwnershipRange(A)
for i in A_start:A_end-1
MatSetValues(A, [i], [i-1, i], [-1. 1.] / Δx, INSERT_VALUES) # MatSetValues(A, I, J, V, INSERT_VALUES)
A_start, A_end = getOwnershipRange(A)
for i = A_start:(A_end - 1)
setValues(A, [i], [i - 1, i], [-1.0 1.0] / Δx, INSERT_VALUES) # setValues(A, I, J, V, INSERT_VALUES)
end
```

Set boundary condition (only the proc handling index `0` is acting)

```julia
(b_start == 0) && VecSetValue(b, 0, 0.)
(b_start == 0) && setValue(b, 0, 0.0)
```

Assemble matrices

```julia
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)
VecAssemblyBegin(b)
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)
VecAssemblyEnd(b)
assemblyBegin(A, MAT_FINAL_ASSEMBLY)
assemblyBegin(b)
assemblyEnd(A, MAT_FINAL_ASSEMBLY)
assemblyEnd(b)
```

At this point, you can inspect `A` or `b` using a viewer (stdout by default), simply call

```julia
MatView(A)
VecView(b)
matView(A)
vecView(b)
```

Set up the linear solver

```julia
ksp = KSPCreate()
KSPSetOperators(ksp, A, A)
KSPSetFromOptions(ksp)
KSPSetUp(ksp)
ksp = create(KSP)
setOperators(ksp, A, A)
setFromOptions(ksp)
setUp(ksp)
```

Solve the system. We first allocate the solution using the `VecDuplicate` method.

```julia
x = VecDuplicate(b)
KSPSolve(ksp, b, x)
x = duplicate(b)
solve(ksp, b, x)
```

Print the solution

```julia
VecView(x)
vecView(x)
```

Access the solution (this part is under development), getting a Julia array; and then restore it

```julia
array, ref = VecGetArray(x) # do something with array
array, ref = getArray(x) # do something with array
@show array
VecRestoreArray(x, ref)
restoreArray(x, ref)
```

Free memory
Free memory. Note that this call is faculative since, by default,
the julia GC will trigger a call to Petsc `destroy` to each object

```julia
MatDestroy(A)
VecDestroy(b)
VecDestroy(x)
destroy.((ksp, A, b, x))
```

Finalize Petsc
Finalize Petsc. This call is faculative : it will be triggered automatically at the end of the script.

```julia
PetscFinalize()

```
```
Loading

0 comments on commit 6adb596

Please sign in to comment.