[firedrake] Calling PETSc within PyOP2 kernels

[Lawrence Mitchell]

Hi Justin,

> On 3 Nov 2015, at 06:59, Justin Chang <jychang48 at gmail.com> wrote:
> 
> Lawrence,
> 
> So I have implemented PETSc into the PYOP2 kernel, but I am sometimes getting strange errors.
> 
> Attached is the working code plus the reactions kernel file.
> 
> Run the code by typing "python 1D_analytical_ADR_ex1.py <seed> <0/1> <0/1>". The last option is for using either Newtons Method (0) or Variational Inequality (1).

I think the following is the problem.

Your kernel does:

void solve_kernel(const double *psiVar, double *cVar, const double *kVar) {
    ...
    for (i = 0; i < numPsi; i++)
        user.psi[i] = psiVar[i];
    for (i = 0; i < numKeq; i++)
        user.keq[i] = kVar[i];
...
}

Where numPsi is 2 (and later on copies out 3 values into cVar.

However, the way data is passed to the kernel is like so:

void wrap_solve_kernel(int start, int end,
                       double *arg0_0, double *arg0_1,
                       double *arg1_0, double *arg1_1,
                       double *arg1_2, double *arg2_0) {
    for ( int n = start; n < end; n++ ) {
        int i = n;
        solve_kernel(arg0_0 + (i * 1), arg1_0 + (i * 1), arg2_0);
    }
}

Which comes from the par_loop invocation:

      op2.par_loop(reactions_kernel,Q.dof_dset,
                   u0.dat(op2.READ),
                   c.dat(op2.WRITE),
                   k1.dat(op2.READ))

So as we can see, the mixed Function u0 (with two components, each scalar-valued) is passed to the calling wrapper in two parts (arg0_0 and arg0_1).  Similarly c (with three parts) is passed in three bits).

However, only the first components are passed on to the kernel itself (which then reads past the end of arrays).

Unfortunately, it looks like the correct code generation doesn't work for this use case.  Instead, you can make things work by explicitly passing the components to the par_loop and accepting them in the kernel.  Something like:

void solve_kernel(const double *psi0, const double *psi1, double *c0, double *c1, double *c2, ...)

{

...

}

And then write the par_loop as:

u0_A, u0_B = u0.split()
c0, c1, c2 = c.split()

op2.par_loop(reactions_kernel, Q.dof_dset,
             u0_A(op2.READ), u0_B(op2.READ),
             c0(op2.WRITE), c1(op2.WRITE), c2(op2.WRITE), ...)


However, this is clearly not very natural, so we should probably fix the code generation for your use case (which I think is the natural thing).

In fact, if you run in a more debug-friendly mode in firedrake, where PyOP2 does additional run-time type checking, your code throws a more friendly error telling you this is unsupported behaviour.  This should probably be switched back to being the default in firedrake (see https://github.com/firedrakeproject/firedrake/issues/626)

Cheers,

Lawrence
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: Message signed with OpenPGP using GPGMail
URL: <http://mailman.ic.ac.uk/pipermail/firedrake/attachments/20151103/7aee4db1/attachment.sig>