Re: [firedrake] Calling PETSc within PyOP2 kernels

12 Nov 2015

      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

Re: [firedrake] Calling PETSc within PyOP2 kernels

Lawrence Mitchell