A first demo with fancy names
This example does more or less the same things as the first example, but using non regular API function names.
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 difference scheme.
Note that the way we achieve things in the document can be highly improved and the purpose of this example is only demonstrate some method calls to give an overview.
To run this example, execute : mpirun -n your_favorite_positive_integer julia example2.jl
Import package
using MPI
using PetscWrap
Initialize PETSc. Command line arguments passed to Julia are parsed by PETSc. Alternatively, you can also provide "command line arguments by defining them in a string, for instance PetscInitialize("-ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always")
or by providing each argument in separate strings : PetscInitialize(["-ksp_monitor_short", "-ksp_gmres_cgs_refinement_type", "refine_always")
PetscInitialize()
comm = MPI.COMM_WORLD
Number of mesh points and mesh step
n = 11
Δx = 1.0 / (n - 1)
Create a matrix of size (n,n)
and a vector of size (n)
. The autosetup
option triggers a call to setFromOptions
and setUp
A = create_matrix(comm; nrows_glo = n, ncols_glo = n, autosetup = true)
b = create_vector(comm; nrows_glo = n, autosetup = true)
Let's build the right hand side vector. We first get the range of rows of b
handled by the local processor. The rstart, rend = get_range(*)
methods differ a little bit from PETSc API since the two integers it returns are the effective Julia range of rows handled by the local processor. If n
is the total number of rows, then rstart ∈ [1,n]
and rend
is the last row handled by the local processor.
b_start, b_end = get_range(b)
Now let's build the right hand side vector. Their are various ways to do this, this is just one.
n_loc = length(b) ## Note that n_loc = b_end - b_start + 1...
b[b_start:b_end] = 2 * ones(n_loc)
And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices.
A_start, A_end = get_range(A)
for i = A_start:A_end
A[i, (i - 1):i] = [-1.0 1.0] / Δx
end
Set boundary condition (only the proc handling index 1
is acting)
(b_start == 1) && (b[1] = 0.0)
Assemble matrice and vector
assemble!(A)
assemble!(b)
Set up the linear solver
ksp = create_ksp(A; autosetup = true)
Solve the system
x = solve(ksp, b)
Print the solution (here x is still a PetscVec
)
@show x
Convert PetscVec
to Julia Array
(and play with it!)
array = vec2array(x)
Free memory (optional, objects are garbage collected otherwise)
destroy!(A, b, x, ksp)
Note that it's also possible to build a matrix using the COO format as in SparseArrays
:
M = create_matrix(comm; nrows_glo = 3, ncols_glo = 3, autosetup = true)
M_start, M_end = get_range(M)
I = [1, 1, 1, 2, 3]
J = [1, 3, 1, 3, 2]
V = [1, 2, 3, 4, 5]
k = findall(x -> M_start <= x <= M_end, I) # just a trick to allow this example to run in parallel
set_values!(M, I[k], J[k], V[k], ADD_VALUES)
assemble!(M)
@show M
This is very convenient in sequential since you can fill the three vectors I, J, V in your code and decide only at the last moment if you'd like to use SparseArrays
or PetscMat
.
destroy!(M)
This page was generated using Literate.jl.