[petsc-users] A flag for FormFunction.

Thu Jan 14 11:39:35 CST 2010

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