On Feb 8, 2011, at 10:16 AM, Klaus Zimmermann wrote:
> Dear all,
>
> I want to evaluate a function on a grid. Right now I have some custom code
> (see the end of this message).
> I am now thinking of rewriting this using more PETSC facilities because in
> the future we want to extend the
> program to higher dimensions and also unstructured grids.
There is little similarity and little chance of code reuse between using a
structured grid and an unstructured grid at the level of your function
evaluations. So if you truly want to have a structured grid version and
unstructured you should just have different FormFunctions for them.
> As I understand it now there are several
> candidates in PETSC for doing this:
>
> 1) PF
> 2) DALocalFunction
Using the "local" function approach with DMMGSetSNESLocal() is just a way of
hiding the DAVecGetArray() calls from the function code and is good when you
have simple FormFunctions that do no rely on other vectors beside the usual
input and output vectors. In your example below I do not understand why you
have the input and output vectors inside the appctx, usually they are the input
and output arguments to the formfunction set with SNESSetFunction() or
DMMGSetSNES()
> 3) FIAT (?)
>
> Could you please advise on what should be used?
> An additional problem is that several distinct functions should be evaluated
> at the same time due to the
> due to the reuse of intermediate results.
Are these functions associated with different SNES solvers or is the
composition of these "several distinct functions" that defines the equations
you are solving with PETSc? If the later then you just need to write either a
single function that computes everything or have some smart use of inline to
get good performance even with different functions.
Barry
>
> Any help is appreciated!
> Thanks in advance,
> Klaus
>
> ------------------8<------------------8<------------------8<------------------8<------------------8<--------------
>
> PetscErrorCode EvaluatePsiAndGradPsi(AppCtx *user) {
> PetscErrorCode ierr;
> PetscInt i, j, cxs, cxm, cys, cym;
> PetscScalar **lpsi;
> PetscScalar **lGradPsi_x, **lGradPsi_y;
> Vec gc;
> DACoor2d **coors;
> DA cda;
> ierr = DAGetCoordinateDA(user->zGrid, &cda);CHKERRQ(ierr);
> ierr = DAGetCoordinates(user->zGrid, &gc);CHKERRQ(ierr);
> ierr = DAVecGetArray(cda, gc, &coors);CHKERRQ(ierr);
> ierr = DAGetCorners(cda, &cxs, &cys, PETSC_NULL, &cxm, &cym,
> PETSC_NULL);CHKERRQ(ierr);
>
> ierr = DAVecGetArray(user->zGrid, user->psi, &lpsi);CHKERRQ(ierr);
> ierr = DAVecGetArray(user->zGrid, user->gradPsi_x,
> &lGradPsi_x);CHKERRQ(ierr);
> ierr = DAVecGetArray(user->zGrid, user->gradPsi_y,
> &lGradPsi_y);CHKERRQ(ierr);
> for(i=cys; i<cys+cym; i++) {
> for(j=cxs; j<cxs+cxm; j++) {
> PetscReal x = PetscRealPart(coors[i][j].x-coors[i][j].y),
> y = PetscRealPart(coors[i][j].y);
> if(x>=0) {
> ierr = EvaluatePsiAndGradPsi(user, x, y,
> &(lpsi[i][j]),
> &(lGradPsi_x[i][j]),
> &(lGradPsi_y[i][j]));CHKERRQ(ierr);
> }
> }
> }
> ierr = DAVecRestoreArray(user->zGrid, user->gradPsi_x,
> &lGradPsi_x);CHKERRQ(ierr);
> ierr = DAVecRestoreArray(user->zGrid, user->gradPsi_y,
> &lGradPsi_y);CHKERRQ(ierr);
> ierr = DAVecRestoreArray(user->zGrid, user->psi, &lpsi);CHKERRQ(ierr);
> ierr = DAVecRestoreArray(cda, gc, &coors);CHKERRQ(ierr);
> ierr = VecDestroy(gc);CHKERRQ(ierr);
> ierr = DADestroy(cda);CHKERRQ(ierr);
> return 0;
> }
