From 5dd8549c614af9dbb0a4e92215163f3a89e877f0 Mon Sep 17 00:00:00 2001 From: filoferra <66474175+filoferra@users.noreply.github.com> Date: Mon, 28 Sep 2026 18:16:38 +0200 Subject: [PATCH] DMStag: set_values! from prefix views without a copy, length(is), on_lower_side --- CHANGELOG.md | 6 ++- docs/src/man/dmstag.md | 5 +++ docs/src/man/lowlevel_intro.md | 3 +- scripts/renames.jl | 1 + src/dmstag.jl | 73 ++++++++++++++++++++++++++++++++-- src/handles.jl | 14 ++++++- src/public_names.jl | 43 ++++++++++---------- test/dm/dmstag_stencils.jl | 53 +++++++++++++++++++++++- test/dm/dmstag_views.jl | 40 +++++++++++++++++++ 9 files changed, 208 insertions(+), 30 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2bf96981e..d149a3f1e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,12 +29,14 @@ - `LibPETSc.PETSC_UNLIMITED` and `LibPETSc.PETSC_CURRENT`, which PETSc's manual pages ask for. - DMStag locations by axis: `vertex_location(dm)`, `edge_location(dm, a, b)`, `face_location(dm, axis)` and `element_location(dm)`, which mean the same in 2D and 3D, unlike `DMSTAG_DOWN`. - `stencil(dm, loc, I; dof = 0)` builds a `DMStagStencil` from a 1-based element index, allocation free. -- Assembly with DMStag stencils: `set_values!(J, dm, rows, cols, vals, mode)`, `set_values!(v, dm, positions, vals, mode)` and `zero_rows_local!(J, dm, rows, diag)`. They take any `AbstractVector`. +- Assembly with DMStag stencils: `set_values!(J, dm, rows, cols, vals, mode)`, `set_values!(v, dm, positions, vals, mode)` and `zero_rows_local!(J, dm, rows, diag)`. They take any `AbstractVector`, and pass a `Vector` or a prefix view `view(buf, 1:n)` of one without a copy. - `LibPETSc.IS(dm, loc => dof, ...)`, the index set of whole DMStag fields, for `set_fieldsplit_is!`. +- `length(is)` on an `IS`: its global size, as `length` of a `PetscVec`. +- `on_lower_side(loc, axis)`: whether a DMStag location sits on the lower side of its element along `axis`, so on the axes where it has one more index than there are elements. - `local_to_local!(dst, dm, src)` and the in-place `local_to_local!(v, dm)`, to refresh ghost points. - `with_product_coordinates(f, dm)`: read-only access to a DMStag's per-axis coordinates, handed back when `f` returns. - `set_matrix_preallocate_only!(dm, flag)`, so a matrix from `dm` gets its nonzero pattern from the first assembly. -- `with_field_views!(f, dm, vecs...; fields, read, write)`: concretely typed views of DMStag local vectors by `location => dof`, indexed like `stencil`, handed back when `f` returns. About 9 allocations per vector, whatever the grid size. +- `with_field_views!(f, dm, vecs...; fields, read, write)`: concretely typed views of DMStag local vectors by `location => dof`, indexed like `stencil`, handed back when `f` returns. About 9 allocations per vector, whatever the grid size. The docstring shows two DMs checked out by nesting, with one vector handed back mid-scope for a halo exchange. ### Changed diff --git a/docs/src/man/dmstag.md b/docs/src/man/dmstag.md index 37e11d35b..75a671b58 100644 --- a/docs/src/man/dmstag.md +++ b/docs/src/man/dmstag.md @@ -129,6 +129,8 @@ PETSc.element_location(dm) # the element interior `face_location(dm, ndims(dm))` is the last axis in any dimension. An axis outside `1:ndims(dm)` throws an `ArgumentError`. +[`on_lower_side`](@ref)`(loc, axis)` says whether a location sits on the lower side of its element along `axis` (`LEFT`, `DOWN`, `BACK`). Such points have one more index along that axis than there are elements, which is what an "is this column inside the domain" check needs. It takes no DM, so it can't check `axis` against the dimension. + A stencil addresses one unknown: a location, the element index and a component. [`stencil`](@ref) takes the 1-based element index that `corners` and `ghost_corners` use, and PETSc's 0-based component, as [`dof_slot`](@ref) does. Indices are not bounds-checked, so ghost elements work, and the call allocates nothing: ```julia @@ -151,6 +153,8 @@ PETSc.set_values!(b, dm, wall_rows, wall_values) Under `ADD_VALUES`, repeated (row, column) entries in one call are summed. `zero_rows_local!` is collective: a rank that owns no wall passes an empty vector. +A `Vector`, or a prefix view `view(buf, 1:n)` of one, is passed to PETSc without a copy, so one scratch buffer per kind can serve blocks of every size. Any other vector is copied first. + An index set of whole fields, for a field split, takes (location, component) pairs: ```julia @@ -158,6 +162,7 @@ flow = LibPETSc.IS(dm, PETSc.face_location(dm, 1) => 0, PETSc.face_location(dm, 2) => 0, PETSc.element_location(dm) => 0) PETSc.set_fieldsplit_is!(PETSc.pc(ksp), "flow", flow) +length(flow) # the global number of indices ``` ## Setting Coordinates diff --git a/docs/src/man/lowlevel_intro.md b/docs/src/man/lowlevel_intro.md index 71f5c50e9..49d18b0ae 100644 --- a/docs/src/man/lowlevel_intro.md +++ b/docs/src/man/lowlevel_intro.md @@ -163,7 +163,7 @@ LibPETSc.MatDestroy(petsclib, mat) For serial (single-process) objects, the high-level interface handles this automatically via finalizers. -`PETSc.destroy!` works on every handle, including the ones with no high-level layer (`IS`, `AO`, `PF`, `Tao`). Unlike the raw `Destroy` call, it does nothing on a borrowed handle (see `PETSc.owns`), and it is safe to call twice or after a `finalize`/`initialize` cycle. +`PETSc.destroy!` works on every handle, including the ones with no high-level layer (`IS`, `AO`, `PF`, `Tao`). Unlike the raw `Destroy` call, it does nothing on a borrowed handle (see `PETSc.owns`), and it is safe to call twice or after a `finalize`/`initialize` cycle. `length(is)` is an index set's global size. Some functions hand back memory PETSc allocated and say the caller frees it with `PetscFree()`. `PetscFree` is a C macro; `LibPETSc.PetscFree(petsclib, ptr)` does the same from Julia. `DMCreateFieldIS` and `ISColoringGetIS` (with `PETSC_OWN_POINTER`) copy what they return into Julia `Vector`s and free the C arrays themselves. @@ -172,6 +172,7 @@ PETSc.destroy!(::PETSc.LibPETSc.AbstractIS) PETSc.destroy!(::PETSc.LibPETSc.AbstractAO) PETSc.destroy!(::PETSc.LibPETSc.AbstractPF) PETSc.destroy!(::PETSc.LibPETSc.AbstractTao) +Base.length(::PETSc.LibPETSc.AbstractIS) PETSc.LibPETSc.PetscFree ``` diff --git a/scripts/renames.jl b/scripts/renames.jl index bed1505ff..37f16f96e 100644 --- a/scripts/renames.jl +++ b/scripts/renames.jl @@ -489,6 +489,7 @@ const UNCHANGED_PUBLIC = Symbol[ :edge_location, # new in v0.5.2 :face_location, # new in v0.5.2 :element_location, # new in v0.5.2 + :on_lower_side, # new in v0.5.2 :stencil, # new in v0.5.2 :local_to_local!, # new in v0.5.2 :set_matrix_preallocate_only!, # new in v0.5.2 diff --git a/src/dmstag.jl b/src/dmstag.jl index b7b6cc014..f7fb0ce78 100644 --- a/src/dmstag.jl +++ b/src/dmstag.jl @@ -499,6 +499,38 @@ $(doc_external("DMStag/DMStagStencilLocation")) """ element_location(::DMStag) = LibPETSc.DMSTAG_ELEMENT +""" + on_lower_side(loc::LibPETSc.DMStagStencilLocation, axis::Integer) + +Whether `loc` sits on the lower side of its element along `axis`: `true` for +`DMSTAG_LEFT`, `DMSTAG_DOWN_LEFT` and every other location with `LEFT` along axis +1, `DOWN` along axis 2 or `BACK` along axis 3; `false` for `DMSTAG_ELEMENT` and for +the upper and centre sides. Points on the lower side of `axis` have one more index +along it than there are elements, since the upper boundary is stored as the lower +side of element `n + 1`. + +It takes no DM, so it cannot check `axis` against the dimension: any `axis` in +`1:3` is accepted, and a 2D location such as `DMSTAG_LEFT` is not on the lower side +along axis 3. + +```julia +PETSc.on_lower_side(face_location(dm, 1), 1) # true +PETSc.on_lower_side(face_location(dm, 1), 2) # false +``` + +See also [`face_location`](@ref), [`stencil`](@ref). + +# External Links +$(doc_external("DMStag/DMStagStencilLocation")) +""" +function on_lower_side(loc::LibPETSc.DMStagStencilLocation, axis::Integer) + 1 <= axis <= 3 || throw(ArgumentError("axis $axis is not an axis of a DMStag")) + loc == LibPETSc.DMSTAG_NULL_LOCATION && + throw(ArgumentError("DMSTAG_NULL_LOCATION is on no side")) + # the locations count up in base 3, axis 1 fastest: 0 lower, 1 centre, 2 upper + return div(Int(loc) - 1, 3^(axis - 1)) % 3 == 0 +end + function check_axis(::DMStag{PetscLib, N}, axis::Integer) where {PetscLib, N} 1 <= axis <= N || throw(ArgumentError("axis $axis is not an axis of a $(N)D DMStag")) return nothing @@ -555,7 +587,9 @@ end Write the dense block `vals` into `J` at the stencils `rows` × `cols` and return `J`. `vals` is row-major, with `length(rows) * length(cols)` entries. Under `ADD_VALUES`, entries whose row and column repeat are summed. `rows`, `cols` and -`vals` can be any `AbstractVector`; one that is not a `Vector` is copied first. +`vals` can be any `AbstractVector`. A `Vector` or a prefix view of one, +`view(buf, 1:n)`, is passed without a copy, so a scratch buffer can be reused for +blocks of any size; any other vector is copied first. # External Links $(doc_external("DMStag/DMStagMatSetValuesStencil")) @@ -588,7 +622,8 @@ end set_values!(v::AbstractPetscVec, dm::DMStag, positions, vals, mode = INSERT_VALUES) Write `vals` into `v` at the stencils `positions` and return `v`. Both can be any -`AbstractVector` of the same length; one that is not a `Vector` is copied first. +`AbstractVector` of the same length. A `Vector` or a prefix view of one, +`view(buf, 1:n)`, is passed without a copy; any other vector is copied first. # External Links $(doc_external("DMStag/DMStagVecSetValuesStencil")) @@ -662,11 +697,21 @@ function LibPETSc.IS(dm::DMStag{PetscLib}, pairs::AbstractVector{<:Pair}) where return LibPETSc.DMStagCreateISFromStencils(getlib(PetscLib), dm, T(length(stencils)), stencils) end -# The C calls take a `Vector`; any other vector is copied into one +# The C calls take a `Vector` and a count, so a `Vector` passes as it is, and so does +# the parent of a prefix view `view(buf, 1:n)`, the count being `n`. Any other vector +# is copied into one. +const PrefixView{T} = SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true} + stencil_vector(v::Vector{LibPETSc.DMStagStencil}) = v +stencil_vector(v::PrefixView{LibPETSc.DMStagStencil}) = + first(parentindices(v)[1]) == 1 ? parent(v) : collect(v) stencil_vector(v::AbstractVector{LibPETSc.DMStagStencil}) = collect(v) + value_vector(::Type{PetscLib}, v::Vector) where {PetscLib} = eltype(v) === scalartype(PetscLib) ? v : Vector{scalartype(PetscLib)}(v) +value_vector(::Type{PetscLib}, v::PrefixView) where {PetscLib} = + eltype(v) === scalartype(PetscLib) && first(parentindices(v)[1]) == 1 ? parent(v) : + Vector{scalartype(PetscLib)}(v) value_vector(::Type{PetscLib}, v::AbstractVector) where {PetscLib} = Vector{scalartype(PetscLib)}(v) @@ -741,6 +786,28 @@ with_field_views!(dm, x, r; fields = flow, write = (false, true)) do (Vx, Vy, P) end ``` +Vectors of two DMs are checked out by nesting, and a vector is handed back early by +closing its block. Here the stress `τ`, on a second DM `dm_τ` with the same elements, +is computed, handed back for [`local_to_local!`](@ref) to fill its ghosts, then read +by the momentum equation, all inside one checkout of the flow: + +```julia +xx = (element_location(dm_τ) => 0,) +with_field_views!(dm, x, r; fields = flow, write = (false, true)) do (Vx, Vy, P), (Rx, Ry, Rp) + with_field_views!(dm_τ, τ; fields = xx) do (Txx,) + for I in c.lower:c.upper + Txx[I] = 2η * (Vx[I + CartesianIndex(1, 0)] - Vx[I]) / h + end + end + local_to_local!(τ, dm_τ) + with_field_views!(dm_τ, τ; fields = xx, write = false) do (Txx,) + for I in c.lower:c.upper + Rx[I] = (Txx[I] - Txx[I - CartesianIndex(1, 0)] - P[I] + P[I - CartesianIndex(1, 0)]) / h + end + end +end +``` + # External Links $(doc_external("DMStag/DMStagVecGetArray")) $(doc_external("DMStag/DMStagGetLocationSlot")) diff --git a/src/handles.jl b/src/handles.jl index 33bbee3d0..9b3c4a41e 100644 --- a/src/handles.jl +++ b/src/handles.jl @@ -1,7 +1,7 @@ # src/handles.jl # Release for the handles that have no high-level layer of their own: IS, AO, PF # and Tao. They carry `own` like every handle, so `destroy!` behaves as it does -# on a Vec or a KSP. +# on a Vec or a KSP. An IS also has a `length`. import .LibPETSc: AbstractIS, AbstractAO, AbstractPF, AbstractTao @@ -24,6 +24,18 @@ function destroy!(is::AbstractIS{PetscLib}) where {PetscLib} return nothing end +""" + length(is::LibPETSc.AbstractIS) + +The global number of indices in `is`, summed over its ranks, as `length` of a +`PetscVec` is its global length. + +# External Links +$(doc_external("IS/ISGetSize")) +""" +Base.length(is::AbstractIS{PetscLib}) where {PetscLib} = + Int(LibPETSc.ISGetSize(getlib(PetscLib), is)) + """ destroy!(ao::LibPETSc.AbstractAO) diff --git a/src/public_names.jl b/src/public_names.jl index 0cdda16d4..6dfaa657c 100644 --- a/src/public_names.jl +++ b/src/public_names.jl @@ -21,24 +21,25 @@ public @bd_fn, @jacobian_fn, @residual_fn, @simple_fn, AbstractPetscDS, ksp_iterations, l2diff, label, library_info, local_arrays, local_coordinate_array, local_coordinates, local_indices, local_interior_linear_index, local_to_global!, local_to_local!, local_vec, mat_nullspace_create, max_steps, max_time, memtype, - memtype_backend, narrow, ndofs, options_prefix, ownership_range, owns, parse_option, - parse_options, pc, prev_time, project_field!, project_function!, release_local_array, - reset!, reshape_local_array, restore_local_arrays!, save_vtk, scalartype, - set_adapt_type!, set_auxiliary_vec!, set_compute_operators!, set_compute_rhs!, - set_constants!, set_convergence_test!, set_dm!, set_dm_active!, set_equation_type!, - set_error_if_step_fails!, set_exact_final_time!, set_exact_solution!, set_field!, - set_fieldsplit_is!, set_from_options!, set_function!, set_function_domain_error!, - set_ifunction!, set_ijacobian!, set_jacobian!, set_jacobian_nullspace!, - set_jacobian_preconditioner!, set_library!, set_matrix_preallocate_only!, - set_max_snes_failures!, set_max_step_rejections!, set_max_steps!, set_max_time!, - set_monitor!, set_name!, set_nullspace!, set_operators!, set_option!, - set_options_prefix!, set_petscint!, set_post_step!, set_pre_step!, set_problem_type!, - set_residual!, set_rhs_function!, set_rhs_jacobian!, set_shell_apply!, - set_shell_setup!, set_snes_jacobian!, set_snes_local_fem!, set_solution!, - set_step_number!, set_time!, set_timestep!, set_tolerances!, set_type!, - set_uniform_coordinates!, set_user_ctx!, set_values!, setup!, snes, snes_failures, - snes_iterations, solution, solve!, solve_time, star_fd_coloring, stencil, step!, - step_number, step_rejections, tao_usable_after_reinitialize, timestep, tolerances, - type_name, unsafe_local_array, unset_library!, user_ctx, vertex_location, - vtk_merge_tensor, with_field_views!, with_local_array!, with_product_coordinates, - wrap_local_array, zero_rows!, zero_rows_local! + memtype_backend, narrow, ndofs, on_lower_side, options_prefix, ownership_range, owns, + parse_option, parse_options, pc, prev_time, project_field!, project_function!, + release_local_array, reset!, reshape_local_array, restore_local_arrays!, save_vtk, + scalartype, set_adapt_type!, set_auxiliary_vec!, set_compute_operators!, + set_compute_rhs!, set_constants!, set_convergence_test!, set_dm!, set_dm_active!, + set_equation_type!, set_error_if_step_fails!, set_exact_final_time!, + set_exact_solution!, set_field!, set_fieldsplit_is!, set_from_options!, + set_function!, set_function_domain_error!, set_ifunction!, set_ijacobian!, + set_jacobian!, set_jacobian_nullspace!, set_jacobian_preconditioner!, set_library!, + set_matrix_preallocate_only!, set_max_snes_failures!, set_max_step_rejections!, + set_max_steps!, set_max_time!, set_monitor!, set_name!, set_nullspace!, + set_operators!, set_option!, set_options_prefix!, set_petscint!, set_post_step!, + set_pre_step!, set_problem_type!, set_residual!, set_rhs_function!, + set_rhs_jacobian!, set_shell_apply!, set_shell_setup!, set_snes_jacobian!, + set_snes_local_fem!, set_solution!, set_step_number!, set_time!, set_timestep!, + set_tolerances!, set_type!, set_uniform_coordinates!, set_user_ctx!, set_values!, + setup!, snes, snes_failures, snes_iterations, solution, solve!, solve_time, + star_fd_coloring, stencil, step!, step_number, step_rejections, + tao_usable_after_reinitialize, timestep, tolerances, type_name, unsafe_local_array, + unset_library!, user_ctx, vertex_location, vtk_merge_tensor, with_field_views!, + with_local_array!, with_product_coordinates, wrap_local_array, zero_rows!, + zero_rows_local! diff --git a/test/dm/dmstag_stencils.jl b/test/dm/dmstag_stencils.jl index 75e9a9929..f1d5ccfeb 100644 --- a/test/dm/dmstag_stencils.jl +++ b/test/dm/dmstag_stencils.jl @@ -37,6 +37,21 @@ MPI.Initialized() || MPI.Init() @test_throws ArgumentError PETSc.edge_location(dm1, 1, 1) @test_throws ArgumentError PETSc.edge_location(dm3, 2, 2) @test_throws ArgumentError PETSc.edge_location(dm2, 1, 3) + + # on_lower_side against the names: LEFT, DOWN and BACK are the lower sides + lower = ("LEFT", "DOWN", "BACK") + for loc in instances(LibPETSc.DMStagStencilLocation) + loc == LibPETSc.DMSTAG_NULL_LOCATION && continue + name = string(loc) + for axis in 1:3 + @test PETSc.on_lower_side(loc, axis) == occursin(lower[axis], name) + end + end + @test PETSc.on_lower_side(PETSc.face_location(dm2, 2), 2) + @test !any(a -> PETSc.on_lower_side(PETSc.element_location(dm3), a), 1:3) + @test all(a -> PETSc.on_lower_side(PETSc.vertex_location(dm3), a), 1:3) + @test_throws ArgumentError PETSc.on_lower_side(LibPETSc.DMSTAG_LEFT, 4) + @test_throws ArgumentError PETSc.on_lower_side(LibPETSc.DMSTAG_NULL_LOCATION, 1) end @testset "stencil" begin @@ -86,15 +101,49 @@ MPI.Initialized() || MPI.Init() is = LibPETSc.IS(dm, e => 0) @test PETSc.owns(is) - @test LibPETSc.ISGetSize(petsclib, is) == 6 + @test length(is) == 6 + @test length(is) isa Int foreach(PETSc.destroy!, (is, v, A, dm)) # a split of several locations and components is = LibPETSc.IS(dm2, PETSc.face_location(dm2, 1) => 0, PETSc.face_location(dm2, 2) => 0) - @test LibPETSc.ISGetSize(petsclib, is) == 4 * 2 + 3 * 3 + @test length(is) == 4 * 2 + 3 * 3 PETSc.destroy!(is) end + @testset "set_values! from reused buffers" begin + dm = PETSc.DMStag(petsclib, comm, none(2), (3, 2), (0, 0, 1), 1) + e = PETSc.element_location(dm) + at(I...) = PETSc.stencil(dm, e, I) + rows = [at(1, 1), at(2, 1), at(3, 2)] + vals = [1.0, 2.0, 3.0, 4.0, 5.0] + + # a prefix view writes the same block as the Vector it stands for + A = PETSc.PetscMat(dm) + B = PETSc.PetscMat(dm) + PETSc.set_values!(A, dm, view(rows, 1:2), view(rows, 1:2), view(vals, 1:4)) + PETSc.set_values!(B, dm, rows[1:2], rows[1:2], vals[1:4]) + PETSc.assemble!(A) + PETSc.assemble!(B) + @test [A[i, j] for i in 1:2, j in 1:2] == [B[i, j] for i in 1:2, j in 1:2] == [1.0 2.0; 3.0 4.0] + + # and passes the buffer without a copy; a view that is not a prefix is copied + set!(M, n) = PETSc.set_values!(M, dm, view(rows, 1:n), view(rows, 1:n), + view(vals, 1:(n * n)), LibPETSc.ADD_VALUES) + set!(A, 1) + @test @allocated(set!(A, 2)) == 0 + w = PETSc.global_vec(dm) + setv!(x) = PETSc.set_values!(x, dm, view(rows, 1:2), view(vals, 1:2)) + setv!(w) + @test @allocated(setv!(w)) == 0 + PETSc.set_values!(w, dm, view(rows, 2:3), view(vals, 4:5)) + PETSc.assemble!(w) + @test w[1] == 1.0 + @test w[2] == 4.0 + @test w[6] == 5.0 + foreach(PETSc.destroy!, (w, B, A, dm)) + end + foreach(PETSc.destroy!, (dm3, dm2, dm1)) PETSc.finalize(petsclib) end diff --git a/test/dm/dmstag_views.jl b/test/dm/dmstag_views.jl index be02b7a24..a4b0dfa06 100644 --- a/test/dm/dmstag_views.jl +++ b/test/dm/dmstag_views.jl @@ -80,6 +80,46 @@ MPI.Initialized() || MPI.Init() end end + @testset "two DMs, one handed back mid-scope" begin + # the flow DM, and a stress DM with the same elements and one element dof + n = (4, 3) + dm = stag(n, (1, 1, 1)) + dm_τ = stag(n, (0, 0, 1)) + flow = (PETSc.face_location(dm, 1) => 0, PETSc.face_location(dm, 2) => 0, + PETSc.element_location(dm) => 0) + xx = (PETSc.element_location(dm_τ) => 0,) + x = PETSc.local_vec(dm) + r = PETSc.local_vec(dm) + τ = PETSc.local_vec(dm_τ) + c = PETSc.corners(dm) + PETSc.with_field_views!(dm, x; fields = flow[1:1]) do (Vx,) + for I in CartesianIndices(Vx) + Vx[I] = I[1]^2 + end + end + + PETSc.with_field_views!(dm, x, r; fields = flow, write = (false, true)) do (Vx, Vy, P), (Rx, Ry, Rp) + PETSc.with_field_views!(dm_τ, τ; fields = xx) do (Txx,) + for I in c.lower:c.upper + Txx[I] = Vx[I + CartesianIndex(1, 0)] - Vx[I] + end + end + PETSc.local_to_local!(τ, dm_τ) + PETSc.with_field_views!(dm_τ, τ; fields = xx, write = false) do (Txx,) + for I in c.lower:c.upper + Rx[I] = Txx[I] + end + end + end + + # Rx holds Vx[i + 1] - Vx[i] = 2i + 1 on every owned element + got = PETSc.with_field_views!(dm, r; fields = flow[1:1], write = false) do (Rx,) + [Rx[I] for I in c.lower:c.upper] + end + @test got == [2I[1] + 1 for I in c.lower:c.upper] + foreach(PETSc.destroy!, (τ, r, x, dm_τ, dm)) + end + @testset "errors hand the vectors back" begin dm = stag((3, 2), (0, 0, 1)) l = PETSc.local_vec(dm)