Dear Barry,
For example, I am solving a time-dependent problem as
dfdt = F(f(xi,eta),g(xi,eta));
dgdt = G(f(xi,eta),g(xi,eta))
on some moving grid, i.e., at each time step, the grid is not the same
as the one in the previous step.
In FormFunction, there are two parts, the first part is giving the
spatial only residual function F(f,g) and G(f,g); the second part is
to add the time-dependent part into it.
There are some quantity h I would like to save, for example, the
velocity of the moving grid,
h = h(x(xi,eta),y(xi,eta),f(xi,eta))
and I would like to save this h as an output file for tracking purpose.
This h is calculated and stored in the second part, but if the
FormFunction is called by functions other than SNESComputeJacobian,
f,x,y are changed from previous nonliner/linear iterations and h has
been changed, too.
I use the flag "callNumber" to open two different files in which this
h is stored.
if (callNumber==0){
open "useful file" and save current h.
}else{
open "abandon file" and save current h.
}
If FormFunction is called by SNESComputeFunction, it means the f,x,y
are just the content from previous time iteration, on the contrary,
they are the content from previous nonlinear/linear iterations.
The short version of the code for FormFunction is:
#undef __FUNCT__
#define __FUNCT__ "FormFunction"
PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void*dummg)
{
// begin of FormFunctionLocal
// begin of part 1: prepare FormFunctionLocal without temporal terms
GETTING SPATIAL TERMS FOR RESIDUAL FUNCTION
// end of part 1: prepare FormFunctionLocal without temporal terms
// begin of part 2: add temporal terms
{
if (parameters->timeAccuracyOrder == 1){
// getting the velocity of the grids
char v_fileName[127];
PetscViewer v_viewer;
FILE *v_fp;
if(parameters->callNumber==0){
sprintf(v_fileName,
"twgcqt2unffni_grid_velocity_tx%i_ty%i_x%i_y%i_nl%i_nt%i.m", info.mx,
info.my,parameters->mxfield,parameters->myfield,parameters->numberOfLevels,parameters->currentTimeStep);
}else{
sprintf(v_fileName,
"unwanted_twgcqt2unffni_grid_velocity_tx%i_ty%i_x%i_y%i_nl%i_nt%i.m",
info.mx,
info.my,parameters->mxfield,parameters->myfield,parameters->numberOfLevels,parameters->currentTimeStep);
}
ierr = PetscViewerASCIIOpen
(PETSC_COMM_WORLD,v_fileName ,
&v_viewer); CHKERRQ (ierr);
ierr = PetscViewerSetFormat (v_viewer,
PETSC_VIEWER_ASCII_MATLAB); CHKERRQ (ierr);
ierr = PetscViewerASCIIPrintf(v_viewer,"function
[V1,V2,Vr1,Vr2,Vc1,Vc2,ct] =
twgcqt2unffni_grid_velocity_tx%i_ty%i_x%i_y%i_nl%i_nt%i", info.mx,
info.my,parameters->mxfield,parameters->myfield,parameters->numberOfLevels,parameters->currentTimeStep);CHKERRQ(ierr);
ierr =
PetscFOpen(PETSC_COMM_WORLD,v_fileName,"a",&v_fp); CHKERRQ(ierr);
if(v_fp){
PetscViewerASCIIPrintf(v_viewer,"ct =
%1.5f;\n",parameters->currentTime);
for (j = jFirst; j <= jLast; j++){
for (i = iFirst; i <= iLast; i++){
if(j!=pye){
V1 = tempiJ*hx*hy*(x_xi
*0.5*ihy*(field[j+1][i].phi - field[j-1][i].phi) -
x_eta*0.5*ihx*(field[j][i+1].phi - field[j][i-1].phi));
V2 =
-tempiJ*hx*hy*(y_eta*0.5*ihx*(field[j][i+1].phi - field[j][i-1].phi) -
y_xi *0.5*ihy*(field[j+1][i].phi - field[j-1][i].phi));
// begin of adding temporal terms
ADDING TEMPORAL TERMS
// end of adding temporal terms
PetscViewerASCIIPrintf(v_viewer,"V1(%i,%i)=%1.20f;\n",j+1,i+1,V1);
PetscViewerASCIIPrintf(v_viewer,"V2(%i,%i)=%1.20f;\n",j+1,i+1,V2);
}else if (j==pye){
V1 = tempiJ*hx*hy*(x_xi
*0.5*ihy*((2*field[j][i].phi - field[j-1][i].phi) - field[j-1][i].phi)
- x_eta*0.5*ihx*(field[j][i+1].phi - field[j][i-1].phi));
V2 =
-tempiJ*hx*hy*(y_eta*0.5*ihx*(field[j][i+1].phi - field[j][i-1].phi) -
y_xi *0.5*ihy*((2*field[j][i].phi - field[j-1][i].phi) -
field[j-1][i].phi));
PetscViewerASCIIPrintf(v_viewer,"V1(%i,%i)=%1.20f;\n",j+1,i+1,V1);
PetscViewerASCIIPrintf(v_viewer,"V2(%i,%i)=%1.20f;\n",j+1,i+1,V2);
}
}
}
}
ierr = PetscFClose(PETSC_COMM_WORLD,v_fp);CHKERRQ(ierr);
PetscViewerDestroy(v_viewer);
}
}
// end of part 2: add temporal terms
// end of FormFunctionLocal
parameters->callNumber = parameters->callNumber+1;
PetscFunctionReturn(0);
}
And in main loop, the callNumber has been set to zero at the beginning
of each time iteration.
The analytic Jacobian is not provided as
ierr = DMMGSetSNES(dmmg, FormFunction,0);CHKERRQ(ierr);
Thanks very much!
Rebecca
Quoting Barry Smith <bsmith at mcs.anl.gov>:
>
> This isn't really a good solution.
>
> Could you explain what you want to do differently with the function
> and when, then we may have a better solution.
>
> Barry
>
> On Jan 14, 2010, at 10:50 AM, (Rebecca) Xuefei YUAN wrote:
>
>> Dear Barry,
>>
>> Yes, I would like the FormFunction to do different things depending
>> on if it is called by SNESComputeFunction.
>>
>> I found a solution to it, but I am not sure whether it is good/easy
>> enough to track.
>>
>> I put up a integer "callNumber" into my user-defined contents, and
>> set it to be 0 at the beginning of each time iteration.
>>
>> When FormFunction is called, it is first called by
>> SNESComputeFunction, followed by other PETSc functions, so my if
>> condition is
>> if (callNumber==0){
>> do A
>> }else{
>> do B
>> }
>>
>> Before PetscFunctionReturn(0),
>>
>> callNumber + = 1;
>>
>> So far it is good, and it did what I would like to.
>>
>> Is there a better way of doing it?
>>
>> Thanks a lot!
>>
>> Rebecca
>>
>> Quoting Barry Smith <bsmith at mcs.anl.gov>:
>>
>>>
>>> There really isn't a mechanism to know by whom a function is called.
>>>
>>> Do you want this so you know when the function is being called
>>> directly vs when it is being called in computing a Jacobian derivative
>>> with finite differencing? Do you want the function to do something
>>> different in the two cases?
>>>
>>> Barry
>>>
>>> On Jan 14, 2010, at 9:14 AM, (Rebecca) Xuefei YUAN wrote:
>>>
>>>> Dear all,
>>>>
>>>> I would like to setup a flag for the FormFunction like:
>>>>
>>>> if (FormFunction is called by SNESComputeFunction){
>>>> flag = 1;
>>>> }else{
>>>> flag = 0;
>>>> }
>>>>
>>>> Is there any way I can know which function called FormFunction? I
>>>> can tell very clear in debug mode, but I need this information
>>>> in the code.
>>>>
>>>> Thanks very much!
>>>>
>>>> Rebecca
>>>>
>>
>>
>>
>> --
>> (Rebecca) Xuefei YUAN
>> Department of Applied Physics and Applied Mathematics
>> Columbia University
>> Tel:917-399-8032
>> www.columbia.edu/~xy2102
>>
--
(Rebecca) Xuefei YUAN
Department of Applied Physics and Applied Mathematics
Columbia University
Tel:917-399-8032
www.columbia.edu/~xy2102
