[petsc-users] A flag for FormFunction.
Barry Smith
bsmith at mcs.anl.gov
Thu Jan 14 12:39:46 CST 2010
I would never have that file output from the FormFunction(); it
doesn't have much meaning. because FormFunction() can be called for
many reasons. If you want to, for example, save the solution at the
end of each Newton step then use SNESMonitorSet() and provide a
function that does the output, this will be called after each Newton
step.
On Jan 14, 2010, at 11:39 AM, (Rebecca) Xuefei YUAN wrote:
> 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,
This is not really accurate. If the line search is needed then the
FormFunction() will be called by SNESComputeFunction() several times
during a single Newton step. The input to SNESComputeFunction() is not
always meaningful.
Barry
> 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
>
More information about the petsc-users
mailing list