On Sun, Jun 18, 2023 at 2:49 PM Duan Junming <junming.d...@epfl.ch> wrote:
> *From:* knep...@gmail.com <knep...@gmail.com> > *Sent:* Sunday, June 18, 2023 20:35 > *To:* Duan Junming > *Cc:* petsc-users@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@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/>