[petsc-users] Advice on small block matrix vector multiplication

[Matthew Knepley]

On Sun, Jun 18, 2023 at 2:49 PM Duan Junming <junming.duan at epfl.ch> wrote:

> *From:* knepley at gmail.com <knepley at gmail.com>
> *Sent:* Sunday, June 18, 2023 20:35
> *To:* Duan Junming
> *Cc:* petsc-users at mcs.anl.gov
> *Subject:* Re: [petsc-users] Advice on small block matrix vector
> multiplication
>
>
> On Sun, Jun 18, 2023 at 2:13 PM Duan Junming via petsc-users <
> petsc-users at mcs.anl.gov> wrote:
>
>> Dear all,
>>
>
>> I am using DMPlex to manage the unknowns, two fields, one for pressure,
>> and one for velocities with two/three components, defined in each cell.
>> They're represented by polynomials, with N (10~50) dofs for each component.
>>
> I have an operator which can be written in a matrix form (N-by-N, dense),
>> to be applied on the pressure field or each component of the velocities in
>> each cell (the same for each cell and also for each component).
>>
> I was wondering which matrix should be defined to implement the block
>> matrix-vector multiplication, here block means the pressure or the
>> component of the velocities. Maybe a sequential block mat? Could you
>> recommend any example?
>>
> Or I just implement this matrix-vector multiplication by hand?
>>
>
> Dear Matt,
>
> Thank you for your quick reply!
>
>
> 1) It sounds like you have a collocated discretization, meaning p,u,v,w
> are all at the same spots. Is this true?
>
> You're right. They're collocated at the same position.
>
>
> 2) You have a dense operator, like FFT, that can act on each component
>
> Right, a dense operator applied on each component.
>
>
> 3) I think you should make a vector with blocksize d+1 and extract the
> components with
>
>   https://petsc.org/main/manualpages/Vec/VecStrideGather/
>
> then act on them, then restore with
>
>   https://petsc.org/main/manualpages/Vec/VecStrideScatter/
>
> You can use the *All() versions to do all the components at once.
>
> Does this function work with the global/local vector generated from DMPlex?
>

It should. By default, Plex orders all unknowns on a mesh point by field.


> Now the vector is like: p_1, p_2, ..., p_N, u_1, v_1, w_1, ..., u_N, v_N,
> w_N.
>

When you extract a component, it will have one field in order, as above.

  Thanks,

    Matt


>   Thanks,
>
>      Matt
>
>
> Thanks!
>> Junming
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230618/573ad706/attachment.html>