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@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake