[firedrake] Calling PETSc within PyOP2 kernels

[Justin Chang]

Got it, thanks!

Would it also be wise to also reuse the Vecs r,x,lb,ub? Or is the
performance difference going to be negligible for these?

Thanks,
Justin

On Thu, Nov 12, 2015 at 4:29 AM, Lawrence Mitchell <
lawrence.mitchell at imperial.ac.uk> wrote:

> On 12/11/15 10:38, Justin Chang wrote:
> > Okay, I have no problems after reinstalling. I must have
> > unintentionally done something taboo to my previous MacOSX firedrake
> > installation.
>
> OK, good.
>
> I note two things about your code.
>
> 1. If VI is on, you never call VecDestroy on ub.
>
> 2. If you don't want to create and destroy a SNES and Mat for every
> dof, you can do the following, and noting that the FormJacobian and
> FormFunction calls are the same for every kernel (it's just the data
> you feed them that changes):
>
>
> import ctypes
>
> class CtxPointer(ctypes.Structure):
>     _fields_ = [("snes", ctypes.c_voidp),
>                 ("J", ctypes.c_voidp),
>                 ("setup", ctypes.c_int)]
>
> reaction_snes = PETSc.SNES().create(comm=PETSc.COMM_SELF)
> reaction_snes.setOptionsPrefix("reactions_")
> reaction_snes.setFromOptions()
> reaction_mat = PETSc.Mat().create(comm=PETSc.COMM_SELF)
> reaction_mat.setOptionsPrefix("reactions_")
> reaction_mat.setType(PETSc.Mat.Type.SEQDENSE)
> reaction_mat.setSizes((3, 3))
> reaction_mat.setUp()
> reaction_mat.setFromOptions()
>
> AppCtx = CtxPointer()
> AppCtx.snes = reaction_snes.handle
> AppCtx.J = reaction_mat.handle
> AppCtx.setup = 0
>
> ctx = op2.Global(1, data=ctypes.addressof(AppCtx), dtype=np.uintp)
>
>
> Add an extra argument to the reactions par_loop:
>       op2.par_loop(reactions3_kernel,Q.dof_dset,u0_A.dat(op2.READ),u0_
> modify your solve_kernel:B.dat(op2.READ),
>              c_A.dat(op2.WRITE),c_B.dat(op2.WRITE),
>              c_C.dat(op2.WRITE),k1.dat(op2.READ),
>              ctx(op2.READ))
>
> Abusing op2.Globals for this purpose is a bit icky.  We should
> probably add the ability to pass in an AppCtx-like object, but I don't
> know exactly how this should interface with the python-level stuff.
>
> Now you need to modify your solve kernel:
>
> typedef struct {
>     SNES snes;
>     Mat J;
>     int setup;
> } AppCtx;
>
> void solve_kernel(const double *psiA, const double *psiB, double *cA,
> double *cB,
>     double *cC, const double *k1, const uintptr_t *ctx_) {
>   AppCtx           *ctx = (AppCtx *)(ctx_[0]);
>   SNES              snes = ctx->snes;
>   Mat               J = ctx->J;
>   Vec               x,r,lb,ub;
>   const int         numC = 3;
>   const int         VI = %(VI)d;
>   const PetscScalar *xx;
>   PetscInt          i;
>   PetscErrorCode    ierr;
>   /* Initialize */
>   problem_type = %(problem)d;
>   VecCreateSeq(PETSC_COMM_SELF, numC, &x);
>   VecDuplicate(x,&r);
>
>   /* Initialize values */
>   VecSet(x,1);
>   psiValues[0] = psiA[0];
>   psiValues[1] = psiB[0];
>   keq[0] = k1[0];
>
>   /* Solver */
>   if (!ctx->setup) {
>     SNESSetFunction(snes,r,FormResidual,NULL);
>     SNESSetJacobian(snes,J,J,FormJacobian,NULL);
>     ctx->setup = 1;
>   }
>   if (VI) {
>     VecDuplicate(x,&lb);
>     VecDuplicate(x,&ub);
>     VecSet(lb,0);
>     VecSet(ub,2);
>     SNESSetType(snes,SNESVINEWTONRSLS);
>     SNESVISetVariableBounds(snes,lb,ub);
>   }
>   SNESSolve(snes,NULL,x);
>
>   /* Store solution */
>   VecGetArrayRead(x,&xx);
>   cA[0] = xx[0];
>   cB[0] = xx[1];
>   cC[0] = xx[2];
>   VecRestoreArrayRead(x,&xx);
>
>   /* Destroy */
>   VecDestroy(&x);
>   VecDestroy(&r);
>   if (VI) {
>     VecDestroy(&lb);
>     VecDestroy(&ub);
>   }
> }
>
> """
>
> Now you reuse the same SNES and Mat every call, this way you don't
> build and destroy thousands of SNESes.  And you can set up the solver
> in the usual way using the -reactions_ options prefix.
>
> Cheers,
>
> Lawrence
>
>
>
>
> _______________________________________________
> firedrake mailing list
> firedrake at imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
>
-------------- next part --------------
HTML attachment scrubbed and removed