[petsc-users] Howto evaluate function on grid

[Klaus Zimmermann]

On 02/08/2011 06:00 PM, Barry Smith wrote:
> On Feb 8, 2011, at 10:52 AM, Klaus Zimmermann wrote:
>
>> 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))
>>
>     How big are S?
Depending on parameter from 100 to 1000. Also C is dense.
>> 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?
>    NO
Ok.
>>>> 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.
>    You can keep them separate. You can have as many vector inputs and outputs you want (it is only the SNES solvers that need exactly one input and one output). Sometimes storing several vectors "interlaced" gives
> better performance since it uses the cache's better but that is only an optimization.
>
>    If you separate the "iterater" part of the code from the "function" part then you can have a different iterator for structured and unstructured grid but reuse the "function" part.
So is there some general iterator code I could use? With regards to the 
vector layout: After the evaluation I want to calculate quantities like 
PointwiseMult(VecConjugate(u1),u2). I thought that for this it would be 
advantageous to have the output vectors layed out in the same way. Do 
you think the interleaved layout works as well?
>>>> 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;
>>>> }