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