Merge branch 'master' of github.com:gridap/GridapPETSc.jl into adding…

…_mumps_petsc_functions

Project.toml

@@ -21,7 +21,7 @@ GridapDistributed = "0.2"
 MPI = "0.14, 0.15, 0.16, 0.17, 0.18, 0.19"
 PETSc_jll = "3.13"
 PartitionedArrays = "0.2.4"
-SparseMatricesCSR = "0.6.1"
+SparseMatricesCSR = "0.6.6"
 julia = "1.3"
 
 [extras]

src/Environment.jl

@@ -35,8 +35,9 @@ end
 
 function with(f;kwargs...)
   Init(;kwargs...)
-  f()
+  out = f()
   Finalize()
+  out
 end
 
 # In an MPI environment context,

src/PETScArrays.jl

@@ -24,7 +24,7 @@ function Finalize(a::PETScVector)
     if a.comm == MPI.COMM_SELF
        @check_error_code PETSC.VecDestroy(a.vec)
     else
-       @check_error_code PETSC.PetscObjectRegisterDestroy(a.vec[].ptr)
+       @check_error_code PETSC.PetscObjectRegisterDestroy(a.vec[])
     end
     a.initialized = false
     @assert Threads.threadid() == 1
@@ -120,7 +120,7 @@ function _copy!(a::Vector,b::Vec)
   ni = length(a)
   ix = collect(PetscInt,0:(ni-1))
   v = convert(Vector{PetscScalar},a)
-  @check_error_code PETSC.VecGetValues(b.ptr,ni,ix,v)
+  @check_error_code PETSC.VecGetValues(b,ni,ix,v)
   if !(v === a)
     a .= v
   end
@@ -135,12 +135,12 @@ function _copy!(a::Vec,b::Vector)
   ni = length(b)
   ix = collect(PetscInt,0:(ni-1))
   v = convert(Vector{PetscScalar},b)
-  @check_error_code PETSC.VecSetValues(a.ptr,ni,ix,v,PETSC.INSERT_VALUES)
+  @check_error_code PETSC.VecSetValues(a,ni,ix,v,PETSC.INSERT_VALUES)
 end
 
 function _get_local_oh_vector(a::Vec)
   v=PETScVector(MPI.COMM_SELF)
-  @check_error_code PETSC.VecGhostGetLocalForm(a.ptr,v.vec)
+  @check_error_code PETSC.VecGhostGetLocalForm(a,v.vec)
   if v.vec[] != C_NULL  # a is a ghosted vector
     v.ownership=a
     Init(v)
@@ -198,7 +198,7 @@ function Finalize(a::PETScMatrix)
     if a.comm == MPI.COMM_SELF
       @check_error_code PETSC.MatDestroy(a.mat)
     else
-      @check_error_code PETSC.PetscObjectRegisterDestroy(a.mat[].ptr)
+      @check_error_code PETSC.PetscObjectRegisterDestroy(a.mat[])
     end
     a.initialized = false
     @assert Threads.threadid() == 1
@@ -308,20 +308,54 @@ function _copy!(petscmat::Mat,mat::Matrix)
    for i=1:size(mat)[1]
      row[1]=PetscInt(i-1)
      vals .= view(mat,i,:)
-     PETSC.MatSetValues(petscmat.ptr,
+     PETSC.MatSetValues(petscmat,
                         PetscInt(1),
                         row,
                         n,
                         cols,
                         vals,
                         PETSC.INSERT_VALUES)
    end
-   @check_error_code PETSC.MatAssemblyBegin(petscmat.ptr, PETSC.MAT_FINAL_ASSEMBLY)
-   @check_error_code PETSC.MatAssemblyEnd(petscmat.ptr, PETSC.MAT_FINAL_ASSEMBLY)
+   @check_error_code PETSC.MatAssemblyBegin(petscmat, PETSC.MAT_FINAL_ASSEMBLY)
+   @check_error_code PETSC.MatAssemblyEnd(petscmat, PETSC.MAT_FINAL_ASSEMBLY)
+end
+
+function _copy!(petscmat::Mat,mat::AbstractSparseMatrix)
+  Tm  = SparseMatrixCSR{0,PetscScalar,PetscInt}
+  csr = convert(Tm,mat)
+  ia  = csr.rowptr
+  ja  = csr.colval
+  a   = csr.nzval
+  m   = csr.m
+  n   = csr.n
+  maxnnz = maximum( ia[i+1]-ia[i] for i=1:m )
+  row    = Vector{PetscInt}(undef,1)
+  cols   = Vector{PetscInt}(undef,maxnnz)
+  for i=1:size(mat,1)
+    row[1]=PetscInt(i-1)
+    current=1
+    for j=ia[i]+1:ia[i+1]
+      col=ja[j]+1
+      cols[current]=PetscInt(col-1)
+      current=current+1
+    end
+    vals = view(a,ia[i]+1:ia[i+1])
+    PETSC.MatSetValues(
+      petscmat,
+      PetscInt(1),
+      row,
+      ia[i+1]-ia[i],
+      cols,
+      vals,
+      PETSC.INSERT_VALUES)
+  end
+  @check_error_code PETSC.MatAssemblyBegin(petscmat, PETSC.MAT_FINAL_ASSEMBLY)
+  @check_error_code PETSC.MatAssemblyEnd(petscmat, PETSC.MAT_FINAL_ASSEMBLY)
 end
 
 
 
+
 function Base.convert(::Type{PETScMatrix},a::PETScMatrix)
   a
 end

test/PoissonTests.jl

@@ -1,3 +1,4 @@
+using SparseMatricesCSR
 using Gridap
 using Gridap.Algebra
 using Gridap.FESpaces
@@ -52,7 +53,8 @@ function main(parts)
 
       a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ
       l(v) = ∫( v*f )dΩ + ∫( v*g )dΓn
-      op = AffineFEOperator(a,l,U,V)
+      assem=SparseMatrixAssembler(SparseMatrixCSR{0,PetscScalar,PetscInt},Vector{Float64},U,V)
+      op = AffineFEOperator(a,l,U,V,assem)
 
       ls = PETScLinearSolver(mykspsetup)
       fels = LinearFESolver(ls)

test/sequential/PETScArraysTests.jl

@@ -8,7 +8,8 @@ using GridapPETSc: PetscScalar, PetscInt
 using LinearAlgebra
 
 options = "-info"
-GridapPETSc.with(args=split(options)) do
+out_1 = "some output"
+out_2 = GridapPETSc.with(args=split(options)) do
 
   n = 10
   v = PETScVector(n)
@@ -90,6 +91,9 @@ GridapPETSc.with(args=split(options)) do
   @test typeof(C*aj) == PETScVector
   @test C*ap == C*aj
 
+  out_1
 end
 
+@test out_1 === out_2
+
 end # module

test/sequential/PLaplacianDriver.jl

@@ -0,0 +1,56 @@
+module PLaplacianDriver
+
+using Test
+using Gridap
+using Gridap.FESpaces
+using GridapPETSc
+using GridapPETSc: PetscScalar, PetscInt
+using SparseArrays
+using SparseMatricesCSR
+
+options = "-snes_type newtonls -snes_rtol 1.0e-14 -snes_atol 0.0 -snes_monitor -ksp_converged_reason -pc_type ilu -ksp_type gmres -snes_converged_reason"
+
+msg = "some output"
+
+out = GridapPETSc.with(args=split(options)) do
+
+  domain = (0,4,0,4)
+  cells = (100,100)
+  model = CartesianDiscreteModel(domain,cells)
+
+  k = 1
+  u((x,y)) = (x+y)^k
+  σ(∇u) =(1.0+∇u⋅∇u)*∇u
+  dσ(∇du,∇u) = 2*∇u⋅∇du + (1.0+∇u⋅∇u)*∇du
+  f(x) = -divergence(y->σ(∇(u,y)),x)
+
+  Ω = Triangulation(model)
+  dΩ = Measure(Ω,2*k)
+  r(u,v) = ∫( ∇(v)⋅(σ∘∇(u)) - v*f )dΩ
+  j(u,du,v) = ∫( ∇(v)⋅(dσ∘(∇(du),∇(u))) )dΩ
+
+  reffe = ReferenceFE(lagrangian,Float64,k)
+  V = TestFESpace(model,reffe,dirichlet_tags="boundary")
+  U = TrialFESpace(u,V)
+
+  nls = PETScNonlinearSolver()
+
+  op = FEOperator(r,j,U,V)
+  uh = solve(nls,op)
+  eh = u - uh
+  @test sqrt(sum(∫( abs2(eh) )dΩ)) < 1.0e-9
+
+  Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
+  Tv = Vector{PetscScalar}
+  assem = SparseMatrixAssembler(Tm,Tv,U,V)
+  op = FEOperator(r,j,U,V,assem)
+  uh = solve(nls,op)
+  eh = u - uh
+  @test sqrt(sum(∫( abs2(eh) )dΩ)) < 1.0e-9
+
+  msg
+end
+
+@test out === msg
+
+end # module

test/sequential/runtests.jl

@@ -18,6 +18,8 @@ using Test
 
 @time @testset "ElasticityDriver" begin include("ElasticityDriver.jl") end
 
+@time @testset "PLaplacianDriver" begin include("PLaplacianDriver.jl") end
+
 end # module