diff --git a/Project.toml b/Project.toml index a28354c2..86cda67b 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" PETSc_jll = "8fa3689e-f0b9-5420-9873-adf6ccf46f2d" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" [compat] MPI = "0.20" @@ -19,12 +20,12 @@ SparseArrays = "1.10" julia = "1.10" [extras] +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [targets] test = ["ForwardDiff", "UnicodePlots", "Test", "Plots", "SparseDiffTools", "Printf"] diff --git a/src/const.jl b/src/const.jl index a59f2f95..0b2f7549 100644 --- a/src/const.jl +++ b/src/const.jl @@ -56,7 +56,38 @@ const PetscLogDouble = Cdouble @enum InsertMode NOT_SET_VALUES INSERT_VALUES ADD_VALUES MAX_VALUES MIN_VALUES INSERT_ALL_VALUES ADD_ALL_VALUES INSERT_BC_VALUES ADD_BC_VALUES - +@enum MatOption begin + MAT_OPTION_MIN = -3 + MAT_UNUSED_NONZERO_LOCATION_ERR = -2 + MAT_ROW_ORIENTED = -1 + MAT_SYMMETRIC = 1 + MAT_STRUCTURALLY_SYMMETRIC = 2 + MAT_FORCE_DIAGONAL_ENTRIES = 3 + MAT_IGNORE_OFF_PROC_ENTRIES = 4 + MAT_USE_HASH_TABLE = 5 + MAT_KEEP_NONZERO_PATTERN = 6 + MAT_IGNORE_ZERO_ENTRIES = 7 + MAT_USE_INODES = 8 + MAT_HERMITIAN = 9 + MAT_SYMMETRY_ETERNAL = 10 + MAT_NEW_NONZERO_LOCATION_ERR = 11 + MAT_IGNORE_LOWER_TRIANGULAR = 12 + MAT_ERROR_LOWER_TRIANGULAR = 13 + MAT_GETROW_UPPERTRIANGULAR = 14 + MAT_SPD = 15 + MAT_NO_OFF_PROC_ZERO_ROWS = 16 + MAT_NO_OFF_PROC_ENTRIES = 17 + MAT_NEW_NONZERO_LOCATIONS = 18 + MAT_NEW_NONZERO_ALLOCATION_ERR = 19 + MAT_SUBSET_OFF_PROC_ENTRIES = 20 + MAT_SUBMAT_SINGLEIS = 21 + MAT_STRUCTURE_ONLY = 22 + MAT_SORTED_FULL = 23 + MAT_FORM_EXPLICIT_TRANSPOSE = 24 + MAT_STRUCTURAL_SYMMETRY_ETERNAL = 25 + MAT_SPD_ETERNAL = 26 + MAT_OPTION_MAX = 27 +end @enum NormType begin NORM_1 = 0 diff --git a/src/ksp.jl b/src/ksp.jl index aa41b5b6..e13100e0 100644 --- a/src/ksp.jl +++ b/src/ksp.jl @@ -216,7 +216,7 @@ struct Fn_KSPComputeOperators{T} end return r_rnorm[] end - function solve!(x::AbstractVec{$PetscScalar}, ksp::KSP{$PetscScalar}, b::AbstractVec{$PetscScalar}) + function solve!(x::AbstractVec{$PetscScalar}, ksp::KSP{$PetscScalar, LT}, b::AbstractVec{$PetscScalar}) where LT with(ksp.opts) do @chk ccall((:KSPSolve, $libpetsc), PetscErrorCode, (CKSP, CVec, CVec), ksp, b, x) @@ -241,19 +241,29 @@ struct Fn_KSPComputeOperators{T} end return x end + function solve!(x::AbstractVec{$PetscScalar}, tksp::LinearAlgebra.AdjointFactorization{T,K}, b::AbstractVec{$PetscScalar}) where {T,K <: KSP{$PetscScalar}} + ksp = parent(tksp) + with(ksp.opts) do + @chk ccall((:KSPSolveTranspose, $libpetsc), PetscErrorCode, + (CKSP, CVec, CVec), ksp, b, x) + end + return x + end + end # no generic Adjoint solve defined, but for Real we can use Adjoint solve!(x::AbstractVec{T}, aksp::Adjoint{T,K}, b::AbstractVec{T}) where {K <: KSP{T}} where {T<:Real} = solve!(x, transpose(parent(aksp)), b) -const KSPAT{T, LT} = Union{KSP{T, LT}, Transpose{T, KSP{T, LT}}, Adjoint{T, KSP{T, LT}}} +const KSPAT{T, LT} = Union{KSP{T, LT}, Transpose{T, KSP{T, LT}}, Adjoint{T, KSP{T, LT}}, LinearAlgebra.AdjointFactorization{T, KSP{T, LT}}} LinearAlgebra.ldiv!(x::AbstractVec{T}, ksp::KSPAT{T, LT}, b::AbstractVec{T}) where {T, LT} = solve!(x, ksp, b) function LinearAlgebra.ldiv!(x::AbstractVector{T}, ksp::KSPAT{T, LT}, b::AbstractVector{T}) where {T, LT} parent(solve!(AbstractVec(x), ksp, AbstractVec(b))) end Base.:\(ksp::KSPAT{T, LT}, b::AbstractVector{T}) where {T, LT} = ldiv!(similar(b), ksp, b) +#Base.:\(kspt::LinearAlgebra.AdjointFactorization{T,KSPT{T, LT}}, b::AbstractVector{T}) where {T, LT} = ldiv!(similar(b), kspt._A', b) """ diff --git a/src/mat.jl b/src/mat.jl index 65471fa3..76f5e394 100644 --- a/src/mat.jl +++ b/src/mat.jl @@ -239,6 +239,26 @@ function MatSetValuesStencil! end return nothing end + function MatSetOption(mat::AbstractMat{$PetscScalar}, + op::MatOption, + flg::Bool) + + fl = PETSC_TRUE + if !flg + fl = PETSC_FALSE + end + @chk ccall((:MatSetOption, $libpetsc), PetscErrorCode, + (CMat, + MatOption, + PetscBool), + mat, + op, + fl + ) + return nothing + end + + @@ -283,6 +303,7 @@ function MatSetValuesStencil! end return val[] end + function assemblybegin(M::AbstractMat{$PetscScalar}, t::MatAssemblyType=MAT_FINAL_ASSEMBLY) @chk ccall((:MatAssemblyBegin, $libpetsc), PetscErrorCode, (CMat, MatAssemblyType), M, t) return nothing @@ -291,6 +312,7 @@ function MatSetValuesStencil! end @chk ccall((:MatAssemblyEnd, $libpetsc), PetscErrorCode, (CMat, MatAssemblyType), M, t) return nothing end + #function view(mat::AbstractMat{$PetscScalar}, viewer::AbstractViewer{$PetscLib}=ViewerStdout($petsclib, mat.comm)) # if assembled(mat) # @chk ccall((:MatView, $libpetsc), PetscErrorCode, diff --git a/test/old_test.jl b/test/old_test.jl index be465004..2454504b 100644 --- a/test/old_test.jl +++ b/test/old_test.jl @@ -39,7 +39,7 @@ PETSc.initialize() PETSc.destroy(pc_extra) # set new pc, check ref counts are modified by ksp - pc_new = PETSc.PC{Float64}(MPI.COMM_SELF) + pc_new = PETSc.PC{Float64}(MPI.COMM_SELF); @test PETSc.nrefs(pc_new) == 1 PETSc.setpc!(ksp, pc_new) @test PETSc.nrefs(pc_new) == 2 @@ -57,10 +57,8 @@ PETSc.initialize() w = M'*x @test w ≈ S'*x - # FIXME - # y = ksp' \ w - # @test S'*y ≈ w rtol=1e-8 - + y = ksp' \ w + @test S'*y ≈ w rtol=1e-8 f!(y,x) = y .= 2 .*x diff --git a/test/options.jl b/test/options.jl index 79c6b064..8c0d6724 100644 --- a/test/options.jl +++ b/test/options.jl @@ -36,24 +36,22 @@ using PETSc opts["nothing_opt"] = nothing @test "" == opts["nothing_opt"] @test "1" == opts["new_opt"] - - # FIXME - #= + # Check that viewer is working _stdout = stdout (rd, wr) = redirect_stdout() @show opts @test readline(rd) == "opts = #PETSc Option Table entries:" - @test readline(rd) == "-da_grid_x 100" - @test readline(rd) == "-da_grid_y 100" - @test readline(rd) == "-ksp_monitor" - @test readline(rd) == "-ksp_view" - @test readline(rd) == "-mg_levels_0_pc_type ilu" - @test readline(rd) == "-new_opt 1" - @test readline(rd) == "-nothing_opt" - @test readline(rd) == "-pc_mg_levels 1" - @test readline(rd) == "-pc_type mg" + @test readline(rd) == "-da_grid_x 100 # (source: code)" + @test readline(rd) == "-da_grid_y 100 # (source: code)" + @test readline(rd) == "-ksp_monitor # (source: code)" + @test readline(rd) == "-ksp_view # (source: code)" + @test readline(rd) == "-mg_levels_0_pc_type ilu # (source: code)" + @test readline(rd) == "-new_opt 1 # (source: code)" + @test readline(rd) == "-nothing_opt # (source: code)" + @test readline(rd) == "-pc_mg_levels 1 # (source: code)" + @test readline(rd) == "-pc_type mg # (source: code)" @test readline(rd) == "#End of PETSc Option Table entries" @test readline(rd) == "" @@ -66,22 +64,22 @@ using PETSc show(stdout, "text/plain", glo_opts) end @test readline(rd) == "#PETSc Option Table entries:" - @test readline(rd) == "-da_grid_x 100" - @test readline(rd) == "-da_grid_y 100" - @test readline(rd) == "-ksp_monitor" - @test readline(rd) == "-ksp_view" - @test readline(rd) == "-mg_levels_0_pc_type ilu" - @test readline(rd) == "-new_opt 1" - @test readline(rd) == "-nothing_opt" - @test readline(rd) == "-pc_mg_levels 1" - @test readline(rd) == "-pc_type mg" + @test readline(rd) == "-da_grid_x 100 # (source: code)" + @test readline(rd) == "-da_grid_y 100 # (source: code)" + @test readline(rd) == "-ksp_monitor # (source: code)" + @test readline(rd) == "-ksp_view # (source: code)" + @test readline(rd) == "-mg_levels_0_pc_type ilu # (source: code)" + @test readline(rd) == "-new_opt 1 # (source: code)" + @test readline(rd) == "-nothing_opt # (source: code)" + @test readline(rd) == "-pc_mg_levels 1 # (source: code)" + @test readline(rd) == "-pc_type mg # (source: code)" @test readline(rd) == "#End of PETSc Option Table entries" show(stdout, "text/plain", glo_opts) @test readline(rd) == "#No PETSc Option Table entries" redirect_stdout(_stdout) - =# + PETSc.finalize(petsclib) end diff --git a/test/test_dmstag.jl b/test/test_dmstag.jl index cff0b766..fb8ba339 100644 --- a/test/test_dmstag.jl +++ b/test/test_dmstag.jl @@ -312,7 +312,6 @@ end # FIXME @testset "DMStag create matrixes" begin - #= comm = MPI.COMM_WORLD mpirank = MPI.Comm_rank(comm) mpisize = MPI.Comm_size(comm) @@ -320,99 +319,103 @@ end PETSc.initialize(petsclib) PetscScalar = PETSc.scalartype(petsclib) PetscInt = PETSc.inttype(petsclib) - - dm_1D = PETSc.DMStagCreate1d(petsclib,comm,PETSc.DM_BOUNDARY_NONE,200,2,2; stag_grid_x=10); - PETSc.setuniformcoordinatesproduct!(dm_1D, (0,), (10,)) - - A = PETSc.creatematrix(dm_1D); # - @test size(A) == (42,42) - PETSc.assembled(A) - - # set some values using normal indices: - A[1,1] = 1.0 - A[1,10] = 1.0 - - pos1 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_LEFT,1,0,0,1) - pos2 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_RIGHT,4,0,0,0) - pos = [pos1, pos2] - val1 = PetscScalar.([2222.2, 3.2]); - PETSc.DMStagMatSetValuesStencil(dm_1D, A, pos1, pos1, 11.1, PETSc.INSERT_VALUES) - PETSc.DMStagMatSetValuesStencil(dm_1D, A, 1, [pos2], 2, pos, val1, PETSc.INSERT_VALUES) - - PETSc.assemble(A) - @test A[1,10] == 1.0 - - # Reads a value from the matrix, using the stencil structure - @test PETSc.DMStagMatGetValuesStencil(dm_1D, A, pos1, pos1)== PetscScalar(11.1) - @test PETSc.DMStagMatGetValuesStencil(dm_1D, A, 1, [pos2], 2, pos)==val1 - - PETSc.destroy(dm_1D); - - dofCenter = 1; - dofEdge = 1; - dofVertex = 0 - stencilWidth = 1; - dm_2D = PETSc.DMStagCreate2d(petsclib,comm, - PETSc.DM_BOUNDARY_GHOSTED, - PETSc.DM_BOUNDARY_GHOSTED, - 10,11, - PETSc.PETSC_DECIDE,PETSc.PETSC_DECIDE, - dofVertex,dofEdge,dofCenter, - PETSc.DMSTAG_STENCIL_BOX,stencilWidth) - - vec_test_2D_global = PETSc.createglobalvector(dm_2D) - vec_test_2D_local = PETSc.createlocalvector(dm_2D) - - corners = PETSc.getcorners(dm_2D) - ghost_corners = PETSc.getghostcorners(dm_2D) - - - for ix=corners.lower[1]:corners.upper[1] - for iy=corners.lower[2]:corners.upper[2] - local dof - # DOF at the center point - dof = 0; - posA = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_DOWN,ix,iy,0,dof) - value = PetscScalar(ix+10); - PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posA, value, PETSc.INSERT_VALUES) - - dof = 0; - posB = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_LEFT,ix,iy,0,dof) - value = PetscScalar(33); - PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posB, value, PETSc.INSERT_VALUES) - - dof = 0; - posC = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_ELEMENT,ix,iy,0,dof) - value = PetscScalar(44); - PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posC, value, PETSc.INSERT_VALUES) - - + if PetscScalar == Float64 || PetscScalar == Float32 + + dm_1D = PETSc.DMStagCreate1d(petsclib,comm,PETSc.DM_BOUNDARY_NONE,200,2,2; stag_grid_x=10); + PETSc.setuniformcoordinatesproduct!(dm_1D, (0,), (10,)) + + A = PETSc.creatematrix(dm_1D); # + PETSc.MatSetOption(A, PETSc.MAT_NEW_NONZERO_ALLOCATION_ERR, false) + @test size(A) == (42,42) + PETSc.assembled(A) + + # set some values using normal indices: + A[1,1] = 1.0 + A[1,10] = 1.0 + + pos1 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_LEFT,1,0,0,1) + pos2 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_RIGHT,4,0,0,0) + pos = [pos1, pos2] + val1 = PetscScalar.([2222.2, 3.2]); + PETSc.DMStagMatSetValuesStencil(dm_1D, A, pos1, pos1, 11.1, PETSc.INSERT_VALUES) + PETSc.DMStagMatSetValuesStencil(dm_1D, A, 1, [pos2], 2, pos, val1, PETSc.INSERT_VALUES) + + PETSc.assemble(A) + @test A[1,10] == 1.0 + + # Reads a value from the matrix, using the stencil structure + @test PETSc.DMStagMatGetValuesStencil(dm_1D, A, pos1, pos1)== PetscScalar(11.1) + @test PETSc.DMStagMatGetValuesStencil(dm_1D, A, 1, [pos2], 2, pos)==val1 + + PETSc.destroy(dm_1D); + + dofCenter = 1; + dofEdge = 1; + dofVertex = 0 + stencilWidth = 1; + dm_2D = PETSc.DMStagCreate2d(petsclib,comm, + PETSc.DM_BOUNDARY_GHOSTED, + PETSc.DM_BOUNDARY_GHOSTED, + 10,11, + PETSc.PETSC_DECIDE,PETSc.PETSC_DECIDE, + dofVertex,dofEdge,dofCenter, + PETSc.DMSTAG_STENCIL_BOX,stencilWidth) + + vec_test_2D_global = PETSc.createglobalvector(dm_2D) + vec_test_2D_local = PETSc.createlocalvector(dm_2D) + + corners = PETSc.getcorners(dm_2D) + ghost_corners = PETSc.getghostcorners(dm_2D) + + + for ix=corners.lower[1]:corners.upper[1] + for iy=corners.lower[2]:corners.upper[2] + local dof + # DOF at the center point + dof = 0; + posA = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_DOWN,ix,iy,0,dof) + value = PetscScalar(ix+10); + PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posA, value, PETSc.INSERT_VALUES) + + dof = 0; + posB = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_LEFT,ix,iy,0,dof) + value = PetscScalar(33); + PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posB, value, PETSc.INSERT_VALUES) + + dof = 0; + posC = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_ELEMENT,ix,iy,0,dof) + value = PetscScalar(44); + PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posC, value, PETSc.INSERT_VALUES) + + + end end - end - PETSc.assemble(vec_test_2D_global) # assemble global vector + PETSc.assemble(vec_test_2D_global) # assemble global vector - # Add the global values to the local values - PETSc.update!(vec_test_2D_local, vec_test_2D_global,PETSc.INSERT_VALUES) + # Add the global values to the local values + PETSc.update!(vec_test_2D_local, vec_test_2D_global,PETSc.INSERT_VALUES) - # retrieve value back from the local array and check that it agrees with global one - dof = 0; - pos = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_DOWN,2,2,0,dof) - @test PETSc.DMStagVecGetValuesStencil(dm_2D, vec_test_2D_local, pos) == 12.0 + # retrieve value back from the local array and check that it agrees with global one + dof = 0; + pos = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_DOWN,2,2,0,dof) + @test PETSc.DMStagVecGetValuesStencil(dm_2D, vec_test_2D_local, pos) == 12.0 - # Extract an array that holds all DOF's - X2D_dofs = PETSc.DMStagVecGetArray(dm_2D,vec_test_2D_local) # extract arrays with all DOF (mostly for visualizing) - @test X2D_dofs[4,4,1] ≈ PetscScalar(12.0) - @test X2D_dofs[4,4,2] ≈ PetscScalar(33.0) - @test X2D_dofs[4,4,3] ≈ PetscScalar(44.0) + # Extract an array that holds all DOF's + X2D_dofs = PETSc.DMStagVecGetArray(dm_2D,vec_test_2D_local) # extract arrays with all DOF (mostly for visualizing) + @test X2D_dofs[4,4,1] ≈ PetscScalar(12.0) + @test X2D_dofs[4,4,2] ≈ PetscScalar(33.0) + @test X2D_dofs[4,4,3] ≈ PetscScalar(44.0) - # Extract an array of a specific DOF (here a face velocity @ the left) - Xarray = PETSc.DMStagGetGhostArrayLocationSlot(dm_2D,vec_test_2D_local, PETSc.DMSTAG_LEFT, 0) - @test sum(X2D_dofs[:,:,2]-Xarray)==0 # check if the local array is identical to the full array + # Extract an array of a specific DOF (here a face velocity @ the left) + Xarray = PETSc.DMStagGetGhostArrayLocationSlot(dm_2D,vec_test_2D_local, PETSc.DMSTAG_LEFT, 0) + @test sum(X2D_dofs[:,:,2]-Xarray)==0 # check if the local array is identical to the full array - Xarray .= 111. # Set a value @ a specific location - @test vec_test_2D_local[2] ≈ PetscScalar(111) # verify that this is changed in the PETSc Vec + Xarray .= 111. # Set a value @ a specific location + @test vec_test_2D_local[2] ≈ PetscScalar(111) # verify that this is changed in the PETSc Vec + + end PETSc.finalize(petsclib) end end @@ -434,165 +437,166 @@ end PETSc.initialize(petsclib) PetscScalar = PETSc.scalartype(petsclib) PetscInt = PETSc.inttype(petsclib) + if PetscScalar == Float64 || PetscScalar == Float32 + # Define a struct that holds data we need in the local SNES routines below + mutable struct Data_1{PetscScalar,PetscInt} + dm + x_l + f_l + end - # Define a struct that holds data we need in the local SNES routines below - mutable struct Data_1{PetscScalar,PetscInt} - dm - x_l - f_l - end - - user_ctx = Data_1{PetscScalar,PetscInt}(nothing, nothing, nothing); # holds data we need in the local - - function FormRes!(ptr_fx_g, ptr_x_g, user_ctx) + user_ctx = Data_1{PetscScalar,PetscInt}(nothing, nothing, nothing); # holds data we need in the local - # Note that in PETSc, ptr_x_g and ptr_fx_g are pointers to global vectors. - # Copy global to local vectors that are stored in user_ctx - PETSc.update!(user_ctx.x_l, ptr_x_g, PETSc.INSERT_VALUES) - PETSc.update!(user_ctx.f_l, ptr_fx_g, PETSc.INSERT_VALUES) + function FormRes!(ptr_fx_g, ptr_x_g, user_ctx) - # Retrieve arrays from the local vectors - ArrayLocal_x = PETSc.DMStagVecGetArrayRead(user_ctx.dm, user_ctx.x_l); # array with all local x-data - ArrayLocal_f = PETSc.DMStagVecGetArray(user_ctx.dm, user_ctx.f_l); # array with all local residual + # Note that in PETSc, ptr_x_g and ptr_fx_g are pointers to global vectors. + # Copy global to local vectors that are stored in user_ctx + PETSc.update!(user_ctx.x_l, ptr_x_g, PETSc.INSERT_VALUES) + PETSc.update!(user_ctx.f_l, ptr_fx_g, PETSc.INSERT_VALUES) - # Compute local residual - ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx) + # Retrieve arrays from the local vectors + ArrayLocal_x = PETSc.DMStagVecGetArrayRead(user_ctx.dm, user_ctx.x_l); # array with all local x-data + ArrayLocal_f = PETSc.DMStagVecGetArray(user_ctx.dm, user_ctx.f_l); # array with all local residual - # Finalize local arrays - Base.finalize(ArrayLocal_x) - Base.finalize(ArrayLocal_f) + # Compute local residual + ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx) - # Copy local into global residual vector - PETSc.update!(ptr_fx_g, user_ctx.f_l, PETSc.INSERT_VALUES) + # Finalize local arrays + Base.finalize(ArrayLocal_x) + Base.finalize(ArrayLocal_f) - end + # Copy local into global residual vector + PETSc.update!(ptr_fx_g, user_ctx.f_l, PETSc.INSERT_VALUES) - function ComputeLocalResidual(dm, ArrayLocal_x, ArrayLocal_f, user_ctx) - # Compute the local residual. The vectors include ghost points + end - T = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_LEFT, 0); - fT = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_LEFT, 0); + function ComputeLocalResidual(dm, ArrayLocal_x, ArrayLocal_f, user_ctx) + # Compute the local residual. The vectors include ghost points - P = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_ELEMENT, 0); - fP = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_ELEMENT, 0); + T = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_LEFT, 0); + fT = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_LEFT, 0); - # compute the FD stencil - indices = PETSc.DMStagGetIndices(dm); # indices of (center/element) points, not including ghost values. - gc = PETSc.getghostcorners(user_ctx.dm); # start and end of loop including ghost points - c = PETSc.getcorners(user_ctx.dm); # start and end of loop including ghost points + P = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_ELEMENT, 0); + fP = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_ELEMENT, 0); - nT = length(T); # array length - dx = 1.0/(c.size[1]-1); - xp = (gc.lower[1]:gc.upper[1]).*dx; # coordinates including ghost points (to define source term) - F = 6.0.*xp .+ (xp .+1.e-12).^6.0; # define source term function + # compute the FD stencil + indices = PETSc.DMStagGetIndices(dm); # indices of (center/element) points, not including ghost values. + gc = PETSc.getghostcorners(user_ctx.dm); # start and end of loop including ghost points + c = PETSc.getcorners(user_ctx.dm); # start and end of loop including ghost points - # Nonlinear equation @ nodal points - ind = indices.vertex.x; # There is one more "vertex" point - i = ind[2:end-1] - fT[ind[1]] = T[ind[1] ]-0.5; # left BC - fT[ind[end]] = T[ind[end]]-2.0; # right BC - fT[i] = (T[i .+ 1] - 2*T[i] + T[i .- 1])/dx^2 + T[i].*T[i] - F[i] # NL diffusion with source term + nT = length(T); # array length + dx = 1.0/(c.size[1]-1); + xp = (gc.lower[1]:gc.upper[1]).*dx; # coordinates including ghost points (to define source term) + F = 6.0.*xp .+ (xp .+1.e-12).^6.0; # define source term function - # second, non-coupled, equation @ center points - ind = indices.center.x; # There is one more "vertex" point - i = ind[2:end-1]; - fP[ind[1]] = P[ind[1]]-30.; # left BC - fP[ind[end]] = P[ind[end]]-20.; # right BC - fP[i] = (P[i .+ 1] - 2*P[i] + P[i .- 1])/dx^2 # steady state diffusion + # Nonlinear equation @ nodal points + ind = indices.vertex.x; # There is one more "vertex" point + i = ind[2:end-1] + fT[ind[1]] = T[ind[1] ]-0.5; # left BC + fT[ind[end]] = T[ind[end]]-2.0; # right BC + fT[i] = (T[i .+ 1] - 2*T[i] + T[i .- 1])/dx^2 + T[i].*T[i] - F[i] # NL diffusion with source term - end + # second, non-coupled, equation @ center points + ind = indices.center.x; # There is one more "vertex" point + i = ind[2:end-1]; + fP[ind[1]] = P[ind[1]]-30.; # left BC + fP[ind[end]] = P[ind[end]]-20.; # right BC + fP[i] = (P[i .+ 1] - 2*P[i] + P[i .- 1])/dx^2 # steady state diffusion - function ForwardDiff_res(x, user_ctx) - - f = zero(x) # vector of zeros, of same type as x (local vector) + end - ArrayLocal_x = PETSc.DMStagVecGetArray(user_ctx.dm, x); # array with all local x-data - ArrayLocal_f = PETSc.DMStagVecGetArray(user_ctx.dm, f); # array with all local residual + function ForwardDiff_res(x, user_ctx) - ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx); - # As the residual vector f is linked with ArrayLocal_f, we don't need to pass ArrayLocal_f back + f = zero(x) # vector of zeros, of same type as x (local vector) - return f; - end + ArrayLocal_x = PETSc.DMStagVecGetArray(user_ctx.dm, x); # array with all local x-data + ArrayLocal_f = PETSc.DMStagVecGetArray(user_ctx.dm, f); # array with all local residual - function FormJacobian!(ptr_x_g, J, P, user_ctx) + ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx); + # As the residual vector f is linked with ArrayLocal_f, we don't need to pass ArrayLocal_f back - # This requires several steps: - # - # 1) Extract local vector from global solution (x) vector - # 2) Compute local jacobian from the residual routine (note that - # this routine requires julia vectors as input) - - # Extract the local vector - PETSc.update!(user_ctx.x_l, ptr_x_g, PETSc.INSERT_VALUES) - x = PETSc.unsafe_localarray(PetscScalar, user_ctx.x_l.ptr; write=false, read=true) - f_Residual = (x -> ForwardDiff_res(x, user_ctx)); # pass additional arguments into the routine - J_julia = ForwardDiff.jacobian(f_Residual,x); - - # Note: since x is the LOCAL vector, J_julia also ends up having the same size. - ind = PETSc.LocalInGlobalIndices(user_ctx.dm); - if PETSc.assembled(J) == false - J = PETSc.MatSeqAIJ(sparse(J_julia[ind,ind])); - else - J .= sparse(J_julia[ind,ind]); + return f; end - Base.finalize(x_g) - return sparse(J_julia[ind,ind]), ind - end - - # Main part + function FormJacobian!(ptr_x_g, J, P, user_ctx) + + # This requires several steps: + # + # 1) Extract local vector from global solution (x) vector + # 2) Compute local jacobian from the residual routine (note that + # this routine requires julia vectors as input) + + # Extract the local vector + PETSc.update!(user_ctx.x_l, ptr_x_g, PETSc.INSERT_VALUES) + x = PETSc.unsafe_localarray(PetscScalar, user_ctx.x_l.ptr; write=false, read=true) + f_Residual = (x -> ForwardDiff_res(x, user_ctx)); # pass additional arguments into the routine + J_julia = ForwardDiff.jacobian(f_Residual,x); + + # Note: since x is the LOCAL vector, J_julia also ends up having the same size. + ind = PETSc.LocalInGlobalIndices(user_ctx.dm); + if PETSc.assembled(J) == false + J = PETSc.MatSeqAIJ(sparse(J_julia[ind,ind])); + else + J .= sparse(J_julia[ind,ind]); + end + + Base.finalize(x_g) + return sparse(J_julia[ind,ind]), ind + end - # Construct a 1D test case for a coupled P-T diffusion solver, with 1 DOF @ the center & 1 DOF @ faces - nx = 21; - user_ctx.dm = PETSc.DMStagCreate1d(petsclib,comm, - PETSc.DM_BOUNDARY_GHOSTED, - nx, - 1, # DOF @ vertex - 1, # DOF @ center - PETSc.DMSTAG_STENCIL_BOX, - 1); # Stencil width + # Main part + # Construct a 1D test case for a coupled P-T diffusion solver, with 1 DOF @ the center & 1 DOF @ faces + nx = 21; + user_ctx.dm = PETSc.DMStagCreate1d(petsclib,comm, + PETSc.DM_BOUNDARY_GHOSTED, + nx, + 1, # DOF @ vertex + 1, # DOF @ center + PETSc.DMSTAG_STENCIL_BOX, + 1); # Stencil width - x_g = PETSc.createglobalvector(user_ctx.dm) - f_g = PETSc.createglobalvector(user_ctx.dm) - user_ctx.x_l = PETSc.createlocalvector(user_ctx.dm) - user_ctx.f_l = PETSc.createlocalvector(user_ctx.dm) + x_g = PETSc.createglobalvector(user_ctx.dm) + f_g = PETSc.createglobalvector(user_ctx.dm) + user_ctx.x_l = PETSc.createlocalvector(user_ctx.dm) + user_ctx.f_l = PETSc.createlocalvector(user_ctx.dm) + PJ = PETSc.creatematrix(user_ctx.dm); # extract (global) matrix from DMStag + PETSc.MatSetOption(PJ, PETSc.MAT_NEW_NONZERO_ALLOCATION_ERR, false) - PJ = PETSc.creatematrix(user_ctx.dm); # extract (global) matrix from DMStag - J_julia, ind = FormJacobian!(x_g, PJ, PJ, user_ctx) - PJ = PETSc.MatSeqAIJ(J_julia) # assemble non-zero structure + J_julia, ind = FormJacobian!(x_g, PJ, PJ, user_ctx) + PJ = PETSc.MatSeqAIJ(J_julia) # assemble non-zero structure + PETSc.MatSetOption(PJ, PETSc.MAT_NEW_NONZERO_ALLOCATION_ERR, false) - S = PETSc.SNES{PetscScalar}(petsclib, comm; - snes_rtol=1e-12, - snes_view=false, - snes_monitor=true, - ksp_view=true, - # pc_type="none", - snes_monitor_true_residual=false, - snes_converged_reason=false); - S.user_ctx = user_ctx; + S = PETSc.SNES{PetscScalar}(petsclib, comm; + snes_rtol=1e-12, + snes_view=false, + snes_monitor=true, + ksp_view=true, + # pc_type="none", + snes_monitor_true_residual=false, + snes_converged_reason=false); + S.user_ctx = user_ctx; - PETSc.setfunction!(S, FormRes!, f_g) - PETSc.setjacobian!(S, FormJacobian!, PJ, PJ) + PETSc.setfunction!(S, FormRes!, f_g) + PETSc.setjacobian!(S, FormJacobian!, PJ, PJ) - # Solve - PETSc.solve!(x_g, S); - # @show x_g + # Solve + PETSc.solve!(x_g, S); - # check - @test x_g[4] ≈ 29.5 - @test x_g[11] ≈ 0.63797 rtol=1e-4 + # check + @test x_g[4] ≈ 29.5 + @test x_g[11] ≈ 0.63797 rtol=1e-4 + end PETSc.finalize(petsclib) - end end + @testset "DMStag: 2D SNES AD" begin comm = MPI.COMM_WORLD @@ -607,164 +611,167 @@ end PETSc.initialize(petsclib) PetscScalar = PETSc.scalartype(petsclib) PetscInt = PETSc.inttype(petsclib) + if PetscScalar == Float64 || PetscScalar == Float32 - mutable struct Data_2D{PetscScalar,PetscInt} - dm - x_l - f_l - end - user_ctx = Data_2D{PetscScalar,PetscInt}(nothing, nothing, nothing); # holds data we need in the local + mutable struct Data_2D{PetscScalar,PetscInt} + dm + x_l + f_l + end + user_ctx = Data_2D{PetscScalar,PetscInt}(nothing, nothing, nothing); # holds data we need in the local - function FormRes!(ptr_fx_g, ptr_x_g, user_ctx) - # Note that in PETSc, ptr_x_g and ptr_fx_g are pointers to global vectors. + function FormRes!(ptr_fx_g, ptr_x_g, user_ctx) + # Note that in PETSc, ptr_x_g and ptr_fx_g are pointers to global vectors. - # Copy global to local vectors - PETSc.update!(user_ctx.x_l, ptr_x_g, PETSc.INSERT_VALUES) - PETSc.update!(user_ctx.f_l, ptr_fx_g, PETSc.INSERT_VALUES) + # Copy global to local vectors + PETSc.update!(user_ctx.x_l, ptr_x_g, PETSc.INSERT_VALUES) + PETSc.update!(user_ctx.f_l, ptr_fx_g, PETSc.INSERT_VALUES) - # Retrieve arrays from the local vectors - ArrayLocal_x = PETSc.DMStagVecGetArrayRead(user_ctx.dm, user_ctx.x_l); # array with all local x-data - ArrayLocal_f = PETSc.DMStagVecGetArray(user_ctx.dm, user_ctx.f_l); # array with all local residual + # Retrieve arrays from the local vectors + ArrayLocal_x = PETSc.DMStagVecGetArrayRead(user_ctx.dm, user_ctx.x_l); # array with all local x-data + ArrayLocal_f = PETSc.DMStagVecGetArray(user_ctx.dm, user_ctx.f_l); # array with all local residual - # Compute local residual - ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx) + # Compute local residual + ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx) - # Finalize local arrays - Base.finalize(ArrayLocal_x) - Base.finalize(ArrayLocal_f) + # Finalize local arrays + Base.finalize(ArrayLocal_x) + Base.finalize(ArrayLocal_f) - # Copy local into global residual vector - PETSc.update!(ptr_fx_g, user_ctx.f_l, PETSc.INSERT_VALUES) + # Copy local into global residual vector + PETSc.update!(ptr_fx_g, user_ctx.f_l, PETSc.INSERT_VALUES) - end - - function ForwardDiff_res(x, user_ctx) - f = zero(x) # vector of zeros, of same type as x (local vector) + end - ArrayLocal_x = PETSc.DMStagVecGetArray(user_ctx.dm, x); # array with all ocal x-data - ArrayLocal_f = PETSc.DMStagVecGetArray(user_ctx.dm, f); # array with all ocal residual + function ForwardDiff_res(x, user_ctx) + f = zero(x) # vector of zeros, of same type as x (local vector) - ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx); + ArrayLocal_x = PETSc.DMStagVecGetArray(user_ctx.dm, x); # array with all ocal x-data + ArrayLocal_f = PETSc.DMStagVecGetArray(user_ctx.dm, f); # array with all ocal residual - # As the residual vector f is linked with ArrayLocal_f, we don't need to pass ArrayLocal_f back + ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx); - return f; - end + # As the residual vector f is linked with ArrayLocal_f, we don't need to pass ArrayLocal_f back - function ComputeLocalResidual(dm, ArrayLocal_x, ArrayLocal_f, user_ctx) - # Compute the local residual. The vectors include ghost points + return f; + end - # Important! Make sure you retrieve the values from the correct locations. In this example we have a - T = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_ELEMENT, 0); - fT = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_ELEMENT, 0); + function ComputeLocalResidual(dm, ArrayLocal_x, ArrayLocal_f, user_ctx) + # Compute the local residual. The vectors include ghost points - # compute the FD stencil - indices = PETSc.DMStagGetIndices(dm); # indices of (center/element) points, not including ghost values. + # Important! Make sure you retrieve the values from the correct locations. In this example we have a + T = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_ELEMENT, 0); + fT = PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_ELEMENT, 0); - sz = size(user_ctx.dm); # array length - dx = 1.0/(sz[1]-1); - dz = 1.0/(sz[2]-1); + # compute the FD stencil + indices = PETSc.DMStagGetIndices(dm); # indices of (center/element) points, not including ghost values. - # set ghost points for BC'S - bnd = PETSc.DMStagGetBoundaryTypes(user_ctx.dm) - if bnd[1] == PETSc.DM_BOUNDARY_GHOSTED - T[1,:] = T[2,:]; # zero flux; dT/dx=0 - T[end,:] = T[end-1,:]; # zero flux - end + sz = size(user_ctx.dm); # array length + dx = 1.0/(sz[1]-1); + dz = 1.0/(sz[2]-1); - # Diffusion @ center points - indx = indices.center.x; - indz = indices.center.y; + # set ghost points for BC'S + bnd = PETSc.DMStagGetBoundaryTypes(user_ctx.dm) + if bnd[1] == PETSc.DM_BOUNDARY_GHOSTED + T[1,:] = T[2,:]; # zero flux; dT/dx=0 + T[end,:] = T[end-1,:]; # zero flux + end - ix = indx[1:end] # use ghost points in x (required GHOSTED x-boundary) - iz = indz[2:end-1] # center points + # Diffusion @ center points + indx = indices.center.x; + indz = indices.center.y; - # upper and lower BC (including corners) - fT[indx,indz[1]] = T[indx,indz[1]] .- 0.5; # bottom BC - fT[indx,indz[end]] = T[indx,indz[end]] .- 2.0; # top BC + ix = indx[1:end] # use ghost points in x (required GHOSTED x-boundary) + iz = indz[2:end-1] # center points + # upper and lower BC (including corners) + fT[indx,indz[1]] = T[indx,indz[1]] .- 0.5; # bottom BC + fT[indx,indz[end]] = T[indx,indz[end]] .- 2.0; # top BC - fT[ix,iz] = (T[ix .+ 1,iz] - 2*T[ix,iz] + T[ix .- 1,iz])/dx^2 + - (T[ix,iz .+ 1] - 2*T[ix,iz] + T[ix,iz .- 1])/dz^2 - end + fT[ix,iz] = (T[ix .+ 1,iz] - 2*T[ix,iz] + T[ix .- 1,iz])/dx^2 + + (T[ix,iz .+ 1] - 2*T[ix,iz] + T[ix,iz .- 1])/dz^2 - function FormJacobian!(ptr_x_g, J, P, user_ctx) - # This requires several steps: - # - # 1) Extract local vector from global solution (x) vector - # 2) Compute local jacobian from the residual routine (note that - # this routine requires julia vectors as input) - - # Extract the local vector - PETSc.update!(user_ctx.x_l, ptr_x_g, PETSc.INSERT_VALUES) - x = PETSc.unsafe_localarray(PetscScalar, user_ctx.x_l.ptr; write=false, read=true) - - f_Residual = (x -> ForwardDiff_res(x, user_ctx)); # pass additional rguments into the routine - J_julia = ForwardDiff.jacobian(f_Residual,x); - - # Note: since x is the LOCAL vector, J_julia also ends up having the same size. - ind = PETSc.LocalInGlobalIndices(user_ctx.dm); - if PETSc.assembled(J) == false - J = PETSc.MatSeqAIJ(sparse(J_julia[ind,ind])); - else - J .= sparse(J_julia[ind,ind]); end - return sparse(J_julia[ind,ind]), ind, sparse(J_julia) - end + function FormJacobian!(ptr_x_g, J, P, user_ctx) + # This requires several steps: + # + # 1) Extract local vector from global solution (x) vector + # 2) Compute local jacobian from the residual routine (note that + # this routine requires julia vectors as input) + + # Extract the local vector + PETSc.update!(user_ctx.x_l, ptr_x_g, PETSc.INSERT_VALUES) + x = PETSc.unsafe_localarray(PetscScalar, user_ctx.x_l.ptr; write=false, read=true) + + f_Residual = (x -> ForwardDiff_res(x, user_ctx)); # pass additional rguments into the routine + J_julia = ForwardDiff.jacobian(f_Residual,x); + + # Note: since x is the LOCAL vector, J_julia also ends up having the same size. + ind = PETSc.LocalInGlobalIndices(user_ctx.dm); + if PETSc.assembled(J) == false + J = PETSc.MatSeqAIJ(sparse(J_julia[ind,ind])); + else + J .= sparse(J_julia[ind,ind]); + end + + return sparse(J_julia[ind,ind]), ind, sparse(J_julia) + end - # Main routine starts here ---- - - dofVertex = 0 - dofEdge = 0 - dofCenter = 1 - nx,nz = 6,25 - user_ctx.dm = PETSc.DMStagCreate2d(petsclib,comm, - PETSc.DM_BOUNDARY_GHOSTED, - PETSc.DM_BOUNDARY_NONE, - nx, - nz, - 1, - 1, - dofVertex, - dofEdge, - dofCenter, - PETSc.DMSTAG_STENCIL_BOX, - 1) - PJ = PETSc.creatematrix(user_ctx.dm); # extract global) matrix from DMStag - x_g = PETSc.createglobalvector(user_ctx.dm) - f_g = PETSc.createglobalvector(user_ctx.dm) - user_ctx.x_l = PETSc.createlocalvector(user_ctx.dm) - user_ctx.f_l = PETSc.createlocalvector(user_ctx.dm) - - S = PETSc.SNES{PetscScalar}(petsclib, comm; - snes_rtol=1e-12, - snes_monitor=true, - pc_type="none", - snes_monitor_true_residual=true, - snes_converged_reason=false); - S.user_ctx = user_ctx; - - J_julia, ind, J_full = FormJacobian!(x_g, PJ, PJ, user_ctx) - PJ = PETSc.MatSeqAIJ(J_julia) # assemble non-zero structure - - PETSc.setfunction!(S, FormRes!, f_g) - PETSc.setjacobian!(S, FormJacobian!, PJ, PJ) - - # Solve 2D system - sol = PETSc.solve!(x_g, S); - - PETSc.update!(user_ctx.x_l,sol, PETSc.INSERT_VALUES); # copy global solution -> local vector - T2d = PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.x_l, PETSc.DMSTAG_ELEMENT, 0); - - @test T2d[5,5] ≈ 0.75 rtol=1e-3 - # - # ----------------- + # Main routine starts here ---- + + dofVertex = 0 + dofEdge = 0 + dofCenter = 1 + nx,nz = 6,25 + user_ctx.dm = PETSc.DMStagCreate2d(petsclib,comm, + PETSc.DM_BOUNDARY_GHOSTED, + PETSc.DM_BOUNDARY_NONE, + nx, + nz, + 1, + 1, + dofVertex, + dofEdge, + dofCenter, + PETSc.DMSTAG_STENCIL_BOX, + 1) + PJ = PETSc.creatematrix(user_ctx.dm); # extract global) matrix from DMStag + PETSc.MatSetOption(PJ, PETSc.MAT_NEW_NONZERO_ALLOCATION_ERR, false) + + x_g = PETSc.createglobalvector(user_ctx.dm) + f_g = PETSc.createglobalvector(user_ctx.dm) + user_ctx.x_l = PETSc.createlocalvector(user_ctx.dm) + user_ctx.f_l = PETSc.createlocalvector(user_ctx.dm) + + S = PETSc.SNES{PetscScalar}(petsclib, comm; + snes_rtol=1e-12, + snes_monitor=true, + pc_type="none", + snes_monitor_true_residual=true, + snes_converged_reason=false); + S.user_ctx = user_ctx; + + J_julia, ind, J_full = FormJacobian!(x_g, PJ, PJ, user_ctx) + PJ = PETSc.MatSeqAIJ(J_julia) # assemble non-zero structure + PETSc.MatSetOption(PJ, PETSc.MAT_NEW_NONZERO_ALLOCATION_ERR, false) + + PETSc.setfunction!(S, FormRes!, f_g) + PETSc.setjacobian!(S, FormJacobian!, PJ, PJ) + + # Solve 2D system + sol = PETSc.solve!(x_g, S); + + PETSc.update!(user_ctx.x_l,sol, PETSc.INSERT_VALUES); # copy global solution -> local vector + T2d = PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.x_l, PETSc.DMSTAG_ELEMENT, 0); + + @test T2d[5,5] ≈ 0.75 rtol=1e-3 + # + # ----------------- + end PETSc.finalize(petsclib) - end - =# end