DM and DMDA interface

Project.toml

@@ -16,7 +16,8 @@ julia = "1.3"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [targets]
-test = ["ForwardDiff", "UnicodePlots"]
+test = ["ForwardDiff", "UnicodePlots", "Printf"]

examples/ksp/ex50.jl

@@ -0,0 +1,178 @@
+# INCLUDE IN MPI TEST
+
+using PETSc, MPI, Printf
+MPI.Initialized() || MPI.Init()
+PETSc.initialize()
+
+# Set up the RHS for a cell centered scheme for PDE
+#
+#  Δu = -cos(π * x) * cos(π * y)
+#
+# with zeros Neumann boundary conditions
+function rhs!(
+    ksp::PETSc.AbstractKSP{PetscScalar},
+    b_vec::PETSc.AbstractVec{PetscScalar},
+) where {PetscScalar}
+    dm = PETSc.DMDA(ksp)
+    comm = PETSc.getcomm(ksp)
+    corners = PETSc.getcorners(dm)
+    global_size = PETSc.getinfo(dm).global_size[1:2]
+
+    # Grid spacing in each direction
+    h = PetscScalar(1) ./ global_size
+
+    # Build the RHS vector for the cell-centered grid. Since the b_vec is not a
+    # julia vector but a PETSc vector we need to convert it before operating on
+    # it
+    PETSc.map_unsafe_localarray!(b_vec; read = false) do b
+        b = reshape(b, Int64(corners.size[1]), Int64(corners.size[2]))
+        for (iy, y) in enumerate(corners.lower[2]:corners.upper[2])
+            y = (2y - 1) * h[2] / 2
+            for (ix, x) in enumerate(corners.lower[1]:corners.upper[1])
+                x = (2x - 1) * h[1] / 2
+                b[ix, iy] = -cos(π * x) * cos(π * y) * h[1] * h[2]
+            end
+        end
+    end
+
+    # Assemble the PETSc vector
+    PETSc.assemblybegin(b_vec)
+    PETSc.assemblyend(b_vec)
+
+    # Hack to remove constant from the vector since we are using all Neumann
+    # boundary conditions
+    nullspace = PETSc.MatNullSpace{PetscScalar}(comm, PETSc.PETSC_TRUE)
+    PETSc.MatNullSpaceRemove!(nullspace, b_vec)
+    PETSc.destroy(nullspace)
+
+    # zero is the PETSc sucess code
+    return 0
+end
+
+function jacobian!(
+    ksp::PETSc.AbstractKSP{PetscScalar},
+    J::PETSc.AbstractMat{PetscScalar},
+    jac::PETSc.AbstractMat{PetscScalar},
+) where {PetscScalar}
+    dm = PETSc.DMDA(ksp)
+    corners = PETSc.getcorners(dm)
+    PetscInt = eltype(corners.size)
+
+    global_size = PETSc.getinfo(dm).global_size[1:2]
+
+    # Grid spacing in each direction
+    h = PetscScalar(1) ./ global_size
+
+    # XXX: Would Mvectors be better?
+    Sten = PETSc.MatStencil{PetscInt}
+    col = Vector{Sten}(undef, 5)
+    row = Vector{Sten}(undef, 1)
+    val = Vector{PetscScalar}(undef, 5)
+
+    for j in corners.lower[2]:corners.upper[2]
+        for i in corners.lower[1]:corners.upper[1]
+            row[1] = Sten(i = i, j = j)
+            num = 1
+            fill!(val, 0)
+            if i > 1
+                val[num] = -h[2] / h[1]
+                col[num] = Sten(i = i - 1, j = j)
+                num += 1
+            end
+            if i < global_size[1]
+                val[num] = -h[2] / h[1]
+                col[num] = Sten(i = i + 1, j = j)
+                num += 1
+            end
+            if j > 1
+                val[num] = -h[1] / h[2]
+                col[num] = Sten(i = i, j = j - 1)
+                num += 1
+            end
+            if j < global_size[2]
+                val[num] = -h[1] / h[2]
+                col[num] = Sten(i = i, j = j + 1)
+                num += 1
+            end
+            val[num] = -sum(val)
+            col[num] = Sten(i = i, j = j)
+            PETSc.MatSetValuesStencil!(
+                jac,
+                row,
+                col,
+                val,
+                PETSc.INSERT_VALUES;
+                num_cols = num,
+            )
+        end
+    end
+
+    PETSc.assemblybegin(jac)
+    PETSc.assemblyend(jac)
+
+    # Since we have all Neumann there is constant nullspace
+    comm = PETSc.getcomm(ksp)
+    nullspace = PETSc.MatNullSpace{PetscScalar}(comm, PETSc.PETSC_TRUE)
+    PETSc.MatSetNullSpace!(J, nullspace)
+    PETSc.destroy(nullspace)
+    return 0
+end
+
+function main(petsclib; comm = MPI.COMM_WORLD, options...)
+    boundary_type = (PETSc.DM_BOUNDARY_NONE, PETSc.DM_BOUNDARY_NONE)
+    stencil_type = PETSc.DMDA_STENCIL_STAR
+    global_size = (11, 11)
+    procs = (PETSc.PETSC_DECIDE, PETSc.PETSC_DECIDE)
+    dof_per_node = 1
+    stencil_width = 1
+    points_per_proc = (nothing, nothing)
+
+    opts = if MPI.Comm_size(comm) == 1
+        (
+            ksp_monitor = true,
+            ksp_view = true,
+            da_grid_x = 100,
+            da_grid_y = 100,
+            pc_type = "mg",
+            pc_mg_levels = 1,
+            mg_levels_0_pc_type = "ilu",
+            mg_levels_0_pc_factor_levels = 1,
+        )
+    else
+        (
+            da_grid_x = 3,
+            da_grid_y = 3,
+            pc_type = "mg",
+            da_refine = 10,
+            ksp_monitor = nothing,
+            ksp_view = nothing,
+            log_view = nothing,
+        )
+    end
+
+    da = PETSc.DMDACreate2d(
+        petsclib,
+        comm,
+        boundary_type...,
+        stencil_type,
+        global_size...,
+        procs...,
+        dof_per_node,
+        stencil_width,
+        points_per_proc...;
+        opts...,
+    )
+
+    ksp = PETSc.KSP(da; opts...)
+
+    PETSc.KSPSetComputeRHS!(ksp, rhs!)
+    PETSc.KSPSetComputeOperators!(ksp, jacobian!)
+    PETSc.setfromoptions!(ksp)
+    PETSc.solve!(ksp)
+end
+
+for petsclib in PETSc.petsclibs
+    if PETSc.scalartype(petsclib) <: Real
+        main(petsclib)
+    end
+end

src/PETSc.jl

@@ -28,6 +28,8 @@ include("options.jl")
 include("vec.jl")
 include("mat.jl")
 include("matshell.jl")
+include("dm.jl")
+include("dmda.jl")
 include("ksp.jl")
 include("pc.jl")
 include("snes.jl")

src/const.jl

@@ -12,6 +12,7 @@ macro chk(expr)
 end
 
 const PETSC_DEFAULT = -2
+const PETSC_DECIDE = PETSC_DETERMINE = -1
 
 
 @enum PetscBool PETSC_FALSE PETSC_TRUE
@@ -115,4 +116,17 @@ end
     MATOP_SOLVE_ADD=8
     MATOP_SOLVE_TRANSPOSE=9
     MATOP_SOLVE_TRANSPOSE_ADD=10
-end 
+end 
+
+@enum DMBoundaryType begin
+  DM_BOUNDARY_NONE=0
+  DM_BOUNDARY_GHOSTED=1
+  DM_BOUNDARY_MIRROR=2
+  DM_BOUNDARY_PERIODIC=3
+  DM_BOUNDARY_TWIST=4
+end
+
+@enum DMDAStencilType begin
+  DMDA_STENCIL_STAR = 0
+  DMDA_STENCIL_BOX = 1
+end