[petsc-users] Howto evaluate function on grid

[Barry Smith]

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;
> }