[petsc-users] human-readable examples and PETSc idioms

[Matthew Knepley]

On Thu, Jul 9, 2015 at 11:02 PM, Marco Zocca <zocca.marco at gmail.com> wrote:

> Hello,
>    I'm looking for some non-optimized PETSc code; namely, I struggle a
> bit with generalizing the provided examples. On one hand inlined
> elementary operations help keep track of the FLOPs but make for
> hard-to read code.
>
> For instance, KSP example 2, (
>
> http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex2.c.html
> ), isn't there a better abstraction to defining the finite-difference
> grid Laplacian than the following?
>
> for (Ii=Istart; Ii<Iend; Ii++) {
>       v = -1.0; i = Ii/n; j = Ii - i*n;
>       if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
>       if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
>       if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
>       if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
>       v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
>     }
>
> I see this approach as very error prone, yet this is one of the
> simplest examples.
>

You can use MatSetValuesStencil() as is done here:

http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex50.c.html

on line 162.


> Another problematic trend I see is hardcoding constants and
> assumptions on the functions, e.g. in the following (taken from
>
> http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex2.c.html
> , removed comments for space) we are interested in forming a matrix A
> from the values contained in a vector x; i.e. some 'map' operation
> that sweeps x and prepares rows of A at a time according to some fixed
> stencil, with separate treatment of boundary values .
>
> Isn't there an idiom to abstract out this functionality? It's one of
> the most common operations in numerical PDE codes:
>

Actually, this is only common is toy example. Almost no real PDEs have a
fixed stencil
because you have material coefficients, or nonlinearities, or constraints,
etc.

   Matt


> VecGetArrayRead(x,&xx);
>
> VecGetSize(x,&n);
> d    = (PetscReal)(n - 1); d = d*d;
>
> for (i=1; i<n-1; i++) {
>   j[0] = i - 1; j[1] = i; j[2] = i + 1;
>   A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
>   MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
> }
>
> i = 0;   A[0] = 1.0;
> MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
>
> i = n-1; A[0] = 1.0;
>
> MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
>
> VecRestoreArrayRead(x,&xx);
>
>
> I think it would be extremely helpful if the developers of the
> examples spent a little time in adding theory annotations to the code,
> and I would be very thankful if anybody could point me to relevant
> tutorials and literature.
>
> Thank you in advance and kind regards,
> Marco Zocca
>



-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150710/4912c508/attachment.html>