Vector assembly

[fredrikekre]

I wonder what the best way to assemble the RHS is. I have managed to hack something together, but it is not so nice so would like a better way as in the following example:

Consider this 2-part system:

# Full system
1 0 0 | 2
0 2 0 | 4
0 0 1 | 2

# Part 1:
1 0 0 | 2  # row 1
0 1 0 | 2  # row 2

# Part 2:
0 1 0 | 2  # row 2
0 0 1 | 2  # row 3

I would like to do something like:

parts = get_part_ids(sequential, 2)
rows = PRange(parts, 3)
cols = PRange(parts, 3)

I, J, V, rhs_I, rhs_V = map_parts(parts) do p
    if p == 1
        i = [1, 2]
        j = [1, 2]
        v = [1.0, 1.0]
        rhs_i = [1, 2]
        rhs_v = [2.0, 2.0]
    elseif p == 2
        i = [2, 3]
        j = [2, 3]
        v = [1.0, 1.0]
        rhs_i = [2, 3]
        rhs_v = [2.0, 2.0]
    end
    return i, j, v, rhs_i, rhs_v
end

add_gids!(rows, I)
assemble!(I, J, V, rows)
assemble!(rhs_I, rhs_V, rows) # <----------------- MethodError

add_gids!(cols, J)
A = PSparseMatrix(I, J, V, rows, cols; ids=:global)

# This works already, but since vectors are not assembled it gives the wrong answer
b = PVector(rhs_I, rhs_V, rows; ids=:global)

which works for the matrix, but not for the RHS since the marked line is a MethodError. Would it make sense to support this "workflow" and add the missing method for synchronizing 2 vectors like this (similar to how I, J, V is synchronized in assemble!(I, J, V, rows))? Or is there a better way that I have not understood yet. Thanks.

[fredrikekre]

Bump, any thoughts on this?

[fverdugo]

Hi @fredrikekre ! Tanks for resurrecting this. I was not aware of this discussion.

[fverdugo]

Ok! Now I see it. The workflow is different since the vector is dense. Try assemble!(b) just after creating b. Note that there is also the asynchronous versions of assemble! They are pretty equivalent to the x_assembly_begin and x_assembly_end from petsc.

[fredrikekre]

Thanks, that makes sense.