[petsc-users] Matrix-free method in PETSc

[Matthew Knepley]

On Tue, Feb 18, 2020 at 8:20 AM Yuyun Yang <yyang85 at stanford.edu> wrote:

> Thanks for the clarification.
>
> Got one more question: if I have variable coefficients, my stencil will be
> updated at every time step, so will the coefficients in myMatMult. In that
> case, is it necessary to destroy the shell matrix and create it all over
> again, or can I use it as it is, only calling the stencil update function,
> assuming the result will be passed into the matrix operation automatically?
>

You update the information in the context associated with the shell matrix.
No need to destroy it.

  Thanks,

    Matt


> Thanks,
> Yuyun
>
> On 2/18/20, 7:34 AM, "Smith, Barry F." <bsmith at mcs.anl.gov> wrote:
>
>
>
>     > On Feb 17, 2020, at 7:56 AM, Yuyun Yang <yyang85 at stanford.edu>
> wrote:
>     >
>     > Hello,
>     >
>     > I actually have a question about the usage of DMDA since I'm quite
> new to this. I wonder if the DMDA suite of functions can be directly called
> on vectors created from VecCreate?
>
>        Yes, but you have to make sure the ones you create have the same
> sizes and parallel layouts. Generally best to get them from the DMDA or
> VecDuplicate() than the hassle of figuring out sizes.
>
>     > Or the vectors have to be formed by DMDACreateGlobalVector? I'm also
> not sure about what the dof and stencil width arguments do.
>     >
>     > I'm still unsure about the usage of MatCreateShell and
> MatShellSetOperation, since it seems that MyMatMult should still have 3
> inputs just like MatMult (the matrix and two vectors). Since I'm not
> forming the matrix, does that mean the matrix input is meaningless but
> still needs to exist for the sake of this format?
>
>         Well the matrix input is your shell matrix so it likely has
> information you need to do your multiply routine. MatShellGetContext() (No
> you do not want to put your information about the matrix stencil inside
> global variables!)
>
>
>     >
>     > After I create such a shell matrix, can I use it like a regular
> matrix in KSP and utilize preconditioners?
>     >
>     > Thanks!
>     > Yuyun
>     > From: petsc-users <petsc-users-bounces at mcs.anl.gov> on behalf of
> Yuyun Yang <yyang85 at stanford.edu>
>     > Sent: Sunday, February 16, 2020 3:12 AM
>     > To: Smith, Barry F. <bsmith at mcs.anl.gov>
>     > Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
>     > Subject: Re: [petsc-users] Matrix-free method in PETSc
>     >
>     > Thank you, that is very helpful information indeed! I will try it
> and send you my code when it works.
>     >
>     > Best regards,
>     > Yuyun
>     > From: Smith, Barry F. <bsmith at mcs.anl.gov>
>     > Sent: Saturday, February 15, 2020 10:02 PM
>     > To: Yuyun Yang <yyang85 at stanford.edu>
>     > Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
>     > Subject: Re: [petsc-users] Matrix-free method in PETSc
>     >
>     >   Yuyun,
>     >
>     >     If you are speaking about using a finite difference stencil on a
> structured grid where you provide the Jacobian vector products yourself by
> looping over the grid doing the stencil operation we unfortunately do not
> have exactly that kind of example.
>     >
>     >     But it is actually not difficult. I suggest starting with
> src/ts/examples/tests/ex22.c It computes the sparse matrix explicitly with
> FormIJacobian()
>     >
>     >     What you need to do is instead in main() use MatCreateShell()
> and MatShellSetOperation(,MATOP_MULT,(void (*)(void))MyMatMult) then
> provide the routine MyMatMult() to do your stencil based matrix free
> product; note that you can create this new routine by taking the structure
> of IFunction() and reorganizing it to do the Jacobian product instead. You
> will need to get the information about the shell matrix size on each
> process by calling DMDAGetCorners().
>     >
>     >     You will then remove the explicit computation of the Jacobian,
> and also remove the Event stuff since you don't need it.
>     >
>     >      Extending to 2 and 3d is straight forward.
>     >
>     >      Any questions let us know.
>     >
>     >    Barry
>     >
>     >    If you like this would make a great merge request with your code
> to improve our examples.
>     >
>     >
>     > > On Feb 15, 2020, at 9:42 PM, Yuyun Yang <yyang85 at stanford.edu>
> wrote:
>     > >
>     > > Hello team,
>     > >
>     > > I wanted to apply the Krylov subspace method to a matrix-free
> implementation of a stencil, such that the iterative method acts on the
> operation without ever constructing the matrix explicitly (for example,
> when doing backward Euler).
>     > >
>     > > I'm not sure whether there is already an example for that
> somewhere. If so, could you point me to a relevant example?
>     > >
>     > > Thank you!
>     > >
>     > > Best regards,
>     > > Yuyun
>
>
>
>

-- 
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/20200218/ccf391ad/attachment-0001.html>