[petsc-users] Howto evaluate function on grid

[Klaus Zimmermann]

Hi Barry,

thanks for your response. I guess my code excerpt wasn't very good, nor 
was my description. My apologies.

The evaluation in more details goes like this:
Depending on the coordinates x and y (and only on them) I calculate 4 
vectors S1,...,S4.
I have a constant matrix C in the appctx. The three quantities I am 
interested in are then:

1) u1 = (x+y)*VecTDot(S1, MatMult(C, S2))
2) u2 = u1/(x+y) + (x+y)*VecTDot(S3, MatMult(C, S2))
3) u3 = u1/(x+y) + (x+y)*VecTDot(S1, MatMult(C, S4))

With this (hopefully better) description let me answer to your points 
below individually.
On 02/08/2011 05:29 PM, Barry Smith wrote:
> 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.
This is why I hoped to reuse code for the unstructured version: As long 
as I can call a method with coordinates and context I am fine.
I don't really need any solving. Do I still need different FormFunctions?

>> 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()
I agree. This is mostly because I didn't understand the concepts so well 
at the time I wrote this code and one of the reasons why I would like to 
refactor.
In my case there should in principle be three output vectors. All the 
facilities I have seen in petsc only deal with a single output vector. 
Is this correct?
Of course there is an obvious mapping, but I would prefer to keep the 
vectors apart because that way it is easier to deal with the parallel 
layout.

>> 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.
I am not really using any solving at the moment. Please let me know if 
you need more detail.

Thanks again,
Klaus

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