Actual source code: posindep.c

      1: #define PETSCTS_DLL

      3: /*
      4:        Code for Timestepping with implicit backwards Euler.
      5: */
     6:  #include src/ts/tsimpl.h

      8: typedef struct {
      9:   Vec  update;      /* work vector where new solution is formed */
     10:   Vec  func;        /* work vector where F(t[i],u[i]) is stored */
     11:   Vec  rhs;         /* work vector for RHS; vec_sol/dt */

     13:   /* information used for Pseudo-timestepping */

     15:   PetscErrorCode (*dt)(TS,PetscReal*,void*);              /* compute next timestep, and related context */
     16:   void           *dtctx;
     17:   PetscErrorCode (*verify)(TS,Vec,void*,PetscReal*,PetscTruth*); /* verify previous timestep and related context */
     18:   void           *verifyctx;

     20:   PetscReal  initial_fnorm,fnorm;                  /* original and current norm of F(u) */
     21:   PetscReal  fnorm_previous;

     23:   PetscReal  dt_increment;                  /* scaling that dt is incremented each time-step */
     24:   PetscTruth increment_dt_from_initial_dt;
     25: } TS_Pseudo;

     27: /* ------------------------------------------------------------------------------*/

     31: /*@
     32:     TSPseudoComputeTimeStep - Computes the next timestep for a currently running
     33:     pseudo-timestepping process.

     35:     Collective on TS

     37:     Input Parameter:
     38: .   ts - timestep context

     40:     Output Parameter:
     41: .   dt - newly computed timestep

     43:     Level: advanced

     45:     Notes:
     46:     The routine to be called here to compute the timestep should be
     47:     set by calling TSPseudoSetTimeStep().

     49: .keywords: timestep, pseudo, compute

     51: .seealso: TSPseudoDefaultTimeStep(), TSPseudoSetTimeStep()
     52: @*/
     53: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoComputeTimeStep(TS ts,PetscReal *dt)
     54: {
     55:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

     59:   PetscLogEventBegin(TS_PseudoComputeTimeStep,ts,0,0,0);
     60:   (*pseudo->dt)(ts,dt,pseudo->dtctx);
     61:   PetscLogEventEnd(TS_PseudoComputeTimeStep,ts,0,0,0);
     62:   return(0);
     63: }

     66: /* ------------------------------------------------------------------------------*/
     69: /*@C
     70:    TSPseudoDefaultVerifyTimeStep - Default code to verify the quality of the last timestep.

     72:    Collective on TS

     74:    Input Parameters:
     75: +  ts - the timestep context
     76: .  dtctx - unused timestep context
     77: -  update - latest solution vector

     79:    Output Parameters:
     80: +  newdt - the timestep to use for the next step
     81: -  flag - flag indicating whether the last time step was acceptable

     83:    Level: advanced

     85:    Note:
     86:    This routine always returns a flag of 1, indicating an acceptable
     87:    timestep.

     89: .keywords: timestep, pseudo, default, verify

     91: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep()
     92: @*/
     93: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoDefaultVerifyTimeStep(TS ts,Vec update,void *dtctx,PetscReal *newdt,PetscTruth *flag)
     94: {
     96:   *flag = PETSC_TRUE;
     97:   return(0);
     98: }

    103: /*@
    104:     TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.

    106:     Collective on TS

    108:     Input Parameters:
    109: +   ts - timestep context
    110: -   update - latest solution vector

    112:     Output Parameters:
    113: +   dt - newly computed timestep (if it had to shrink)
    114: -   flag - indicates if current timestep was ok

    116:     Level: advanced

    118:     Notes:
    119:     The routine to be called here to compute the timestep should be
    120:     set by calling TSPseudoSetVerifyTimeStep().

    122: .keywords: timestep, pseudo, verify

    124: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoDefaultVerifyTimeStep()
    125: @*/
    126: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscTruth *flag)
    127: {
    128:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

    132:   if (!pseudo->verify) {*flag = PETSC_TRUE; return(0);}

    134:   (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);

    136:   return(0);
    137: }

    139: /* --------------------------------------------------------------------------------*/

    143: static PetscErrorCode TSStep_Pseudo(TS ts,PetscInt *steps,PetscReal *ptime)
    144: {
    145:   Vec            sol = ts->vec_sol;
    147:   PetscInt       i,max_steps = ts->max_steps,its,lits;
    148:   PetscTruth     ok;
    149:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
    150:   PetscReal      current_time_step;
    151:
    153:   *steps = -ts->steps;

    155:   VecCopy(sol,pseudo->update);
    156:   for (i=0; i<max_steps && ts->ptime < ts->max_time; i++) {
    157:     TSPseudoComputeTimeStep(ts,&ts->time_step);
    158:     TSMonitor(ts,ts->steps,ts->ptime,sol);
    159:     current_time_step = ts->time_step;
    160:     while (PETSC_TRUE) {
    161:       ts->ptime  += current_time_step;
    162:       SNESSolve(ts->snes,PETSC_NULL,pseudo->update);
    163:       SNESGetNumberLinearIterations(ts->snes,&lits);
    164:       SNESGetIterationNumber(ts->snes,&its);
    165:       ts->nonlinear_its += its; ts->linear_its += lits;
    166:       TSPseudoVerifyTimeStep(ts,pseudo->update,&ts->time_step,&ok);
    167:       if (ok) break;
    168:       ts->ptime        -= current_time_step;
    169:       current_time_step = ts->time_step;
    170:     }
    171:     VecCopy(pseudo->update,sol);
    172:     ts->steps++;
    173:   }
    174:   TSComputeRHSFunction(ts,ts->ptime,ts->vec_sol,pseudo->func);
    175:   VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
    176:   TSMonitor(ts,ts->steps,ts->ptime,sol);

    178:   *steps += ts->steps;
    179:   *ptime  = ts->ptime;
    180:   return(0);
    181: }

    183: /*------------------------------------------------------------*/
    186: static PetscErrorCode TSDestroy_Pseudo(TS ts)
    187: {
    188:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

    192:   if (pseudo->update) {VecDestroy(pseudo->update);}
    193:   if (pseudo->func) {VecDestroy(pseudo->func);}
    194:   if (pseudo->rhs)  {VecDestroy(pseudo->rhs);}
    195:   PetscFree(pseudo);
    196:   return(0);
    197: }

    200: /*------------------------------------------------------------*/

    202: /*
    203:     This defines the nonlinear equation that is to be solved with SNES

    205:               (U^{n+1} - U^{n})/dt - F(U^{n+1})
    206: */
    209: PetscErrorCode TSPseudoFunction(SNES snes,Vec x,Vec y,void *ctx)
    210: {
    211:   TS     ts = (TS) ctx;
    212:   PetscScalar mdt = 1.0/ts->time_step,*unp1,*un,*Funp1;
    214:   PetscInt i,n;

    217:   /* apply user provided function */
    218:   TSComputeRHSFunction(ts,ts->ptime,x,y);
    219:   /* compute (u^{n+1) - u^{n})/dt - F(u^{n+1}) */
    220:   VecGetArray(ts->vec_sol,&un);
    221:   VecGetArray(x,&unp1);
    222:   VecGetArray(y,&Funp1);
    223:   VecGetLocalSize(x,&n);
    224:   for (i=0; i<n; i++) {
    225:     Funp1[i] = mdt*(unp1[i] - un[i]) - Funp1[i];
    226:   }
    227:   VecRestoreArray(ts->vec_sol,&un);
    228:   VecRestoreArray(x,&unp1);
    229:   VecRestoreArray(y,&Funp1);

    231:   return(0);
    232: }

    234: /*
    235:    This constructs the Jacobian needed for SNES

    237:              J = I/dt - J_{F}   where J_{F} is the given Jacobian of F.
    238: */
    241: PetscErrorCode TSPseudoJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
    242: {
    243:   TS             ts = (TS) ctx;

    247:   /* construct users Jacobian */
    248:   TSComputeRHSJacobian(ts,ts->ptime,x,AA,BB,str);

    250:   /* shift and scale Jacobian */
    251:   TSScaleShiftMatrices(ts,*AA,*BB,*str);
    252:   return(0);
    253: }

    258: static PetscErrorCode TSSetUp_Pseudo(TS ts)
    259: {
    260:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

    264:   /* SNESSetFromOptions(ts->snes); */
    265:   VecDuplicate(ts->vec_sol,&pseudo->update);
    266:   VecDuplicate(ts->vec_sol,&pseudo->func);
    267:   SNESSetFunction(ts->snes,pseudo->func,TSPseudoFunction,ts);
    268:   SNESSetJacobian(ts->snes,ts->A,ts->B,TSPseudoJacobian,ts);
    269:   return(0);
    270: }
    271: /*------------------------------------------------------------*/

    275: PetscErrorCode TSPseudoDefaultMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
    276: {
    277:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

    281:   (*PetscHelpPrintf)(ts->comm,"TS %D dt %g time %g fnorm %g\n",step,ts->time_step,ptime,pseudo->fnorm);
    282:   return(0);
    283: }

    287: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts)
    288: {
    289:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
    291:   PetscTruth flg;

    295:   PetscOptionsHead("Pseudo-timestepping options");
    296:     PetscOptionsName("-ts_monitor","Monitor convergence","TSPseudoDefaultMonitor",&flg);
    297:     if (flg) {
    298:       TSSetMonitor(ts,TSPseudoDefaultMonitor,PETSC_NULL,PETSC_NULL);
    299:     }
    300:     PetscOptionsName("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",&flg);
    301:     if (flg) {
    302:       TSPseudoIncrementDtFromInitialDt(ts);
    303:     }
    304:     PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,0);
    305:   PetscOptionsTail();

    307:   return(0);
    308: }

    312: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
    313: {
    315:   return(0);
    316: }

    318: /* ----------------------------------------------------------------------------- */
    321: /*@C
    322:    TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
    323:    last timestep.

    325:    Collective on TS

    327:    Input Parameters:
    328: +  ts - timestep context
    329: .  dt - user-defined function to verify timestep
    330: -  ctx - [optional] user-defined context for private data
    331:          for the timestep verification routine (may be PETSC_NULL)

    333:    Level: advanced

    335:    Calling sequence of func:
    336: .  func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscTruth *flag);

    338: .  update - latest solution vector
    339: .  ctx - [optional] timestep context
    340: .  newdt - the timestep to use for the next step
    341: .  flag - flag indicating whether the last time step was acceptable

    343:    Notes:
    344:    The routine set here will be called by TSPseudoVerifyTimeStep()
    345:    during the timestepping process.

    347: .keywords: timestep, pseudo, set, verify

    349: .seealso: TSPseudoDefaultVerifyTimeStep(), TSPseudoVerifyTimeStep()
    350: @*/
    351: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscTruth*),void* ctx)
    352: {
    353:   PetscErrorCode ierr,(*f)(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal *,PetscTruth *),void *);

    358:   PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",(void (**)(void))&f);
    359:   if (f) {
    360:     (*f)(ts,dt,ctx);
    361:   }
    362:   return(0);
    363: }

    367: /*@
    368:     TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
    369:     dt when using the TSPseudoDefaultTimeStep() routine.

    371:    Collective on TS

    373:     Input Parameters:
    374: +   ts - the timestep context
    375: -   inc - the scaling factor >= 1.0

    377:     Options Database Key:
    378: $    -ts_pseudo_increment <increment>

    380:     Level: advanced

    382: .keywords: timestep, pseudo, set, increment

    384: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
    385: @*/
    386: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
    387: {
    388:   PetscErrorCode ierr,(*f)(TS,PetscReal);

    393:   PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",(void (**)(void))&f);
    394:   if (f) {
    395:     (*f)(ts,inc);
    396:   }
    397:   return(0);
    398: }

    402: /*@
    403:     TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
    404:     is computed via the formula
    405: $         dt = initial_dt*initial_fnorm/current_fnorm
    406:       rather than the default update,
    407: $         dt = current_dt*previous_fnorm/current_fnorm.

    409:    Collective on TS

    411:     Input Parameter:
    412: .   ts - the timestep context

    414:     Options Database Key:
    415: $    -ts_pseudo_increment_dt_from_initial_dt

    417:     Level: advanced

    419: .keywords: timestep, pseudo, set, increment

    421: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
    422: @*/
    423: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoIncrementDtFromInitialDt(TS ts)
    424: {
    425:   PetscErrorCode ierr,(*f)(TS);

    430:   PetscObjectQueryFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",(void (**)(void))&f);
    431:   if (f) {
    432:     (*f)(ts);
    433:   }
    434:   return(0);
    435: }

    440: /*@C
    441:    TSPseudoSetTimeStep - Sets the user-defined routine to be
    442:    called at each pseudo-timestep to update the timestep.

    444:    Collective on TS

    446:    Input Parameters:
    447: +  ts - timestep context
    448: .  dt - function to compute timestep
    449: -  ctx - [optional] user-defined context for private data
    450:          required by the function (may be PETSC_NULL)

    452:    Level: intermediate

    454:    Calling sequence of func:
    455: .  func (TS ts,PetscReal *newdt,void *ctx);

    457: .  newdt - the newly computed timestep
    458: .  ctx - [optional] timestep context

    460:    Notes:
    461:    The routine set here will be called by TSPseudoComputeTimeStep()
    462:    during the timestepping process.

    464: .keywords: timestep, pseudo, set

    466: .seealso: TSPseudoDefaultTimeStep(), TSPseudoComputeTimeStep()
    467: @*/
    468: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void* ctx)
    469: {
    470:   PetscErrorCode ierr,(*f)(TS,PetscErrorCode (*)(TS,PetscReal *,void *),void *);

    475:   PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",(void (**)(void))&f);
    476:   if (f) {
    477:     (*f)(ts,dt,ctx);
    478:   }
    479:   return(0);
    480: }

    482: /* ----------------------------------------------------------------------------- */

    488: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void* ctx)
    489: {
    490:   TS_Pseudo *pseudo;

    493:   pseudo              = (TS_Pseudo*)ts->data;
    494:   pseudo->verify      = dt;
    495:   pseudo->verifyctx   = ctx;
    496:   return(0);
    497: }

    503: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
    504: {
    505:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

    508:   pseudo->dt_increment = inc;
    509:   return(0);
    510: }

    516: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
    517: {
    518:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

    521:   pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
    522:   return(0);
    523: }

    530: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void* ctx)
    531: {
    532:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

    535:   pseudo->dt      = dt;
    536:   pseudo->dtctx   = ctx;
    537:   return(0);
    538: }

    541: /* ----------------------------------------------------------------------------- */

    546: PetscErrorCode PETSCTS_DLLEXPORT TSCreate_Pseudo(TS ts)
    547: {
    548:   TS_Pseudo  *pseudo;

    552:   ts->ops->destroy         = TSDestroy_Pseudo;
    553:   ts->ops->view            = TSView_Pseudo;

    555:   if (ts->problem_type == TS_LINEAR) {
    556:     SETERRQ(PETSC_ERR_ARG_WRONG,"Only for nonlinear problems");
    557:   }
    558:   if (!ts->A) {
    559:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set Jacobian");
    560:   }

    562:   ts->ops->setup           = TSSetUp_Pseudo;
    563:   ts->ops->step            = TSStep_Pseudo;
    564:   ts->ops->setfromoptions  = TSSetFromOptions_Pseudo;

    566:   /* create the required nonlinear solver context */
    567:   SNESCreate(ts->comm,&ts->snes);

    569:   PetscNew(TS_Pseudo,&pseudo);
    570:   PetscLogObjectMemory(ts,sizeof(TS_Pseudo));

    572:   PetscMemzero(pseudo,sizeof(TS_Pseudo));
    573:   ts->data = (void*)pseudo;

    575:   pseudo->dt_increment                 = 1.1;
    576:   pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
    577:   pseudo->dt                           = TSPseudoDefaultTimeStep;

    579:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",
    580:                     "TSPseudoSetVerifyTimeStep_Pseudo",
    581:                      TSPseudoSetVerifyTimeStep_Pseudo);
    582:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",
    583:                     "TSPseudoSetTimeStepIncrement_Pseudo",
    584:                      TSPseudoSetTimeStepIncrement_Pseudo);
    585:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",
    586:                     "TSPseudoIncrementDtFromInitialDt_Pseudo",
    587:                      TSPseudoIncrementDtFromInitialDt_Pseudo);
    588:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStep_C",
    589:                     "TSPseudoSetTimeStep_Pseudo",
    590:                      TSPseudoSetTimeStep_Pseudo);
    591:   return(0);
    592: }

    597: /*@C
    598:    TSPseudoDefaultTimeStep - Default code to compute pseudo-timestepping.
    599:    Use with TSPseudoSetTimeStep().

    601:    Collective on TS

    603:    Input Parameters:
    604: .  ts - the timestep context
    605: .  dtctx - unused timestep context

    607:    Output Parameter:
    608: .  newdt - the timestep to use for the next step

    610:    Level: advanced

    612: .keywords: timestep, pseudo, default

    614: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
    615: @*/
    616: PetscErrorCode PETSCTS_DLLEXPORT TSPseudoDefaultTimeStep(TS ts,PetscReal* newdt,void* dtctx)
    617: {
    618:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
    619:   PetscReal inc = pseudo->dt_increment,fnorm_previous = pseudo->fnorm_previous;

    623:   TSComputeRHSFunction(ts,ts->ptime,ts->vec_sol,pseudo->func);
    624:   VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
    625:   if (pseudo->initial_fnorm == 0.0) {
    626:     /* first time through so compute initial function norm */
    627:     pseudo->initial_fnorm = pseudo->fnorm;
    628:     fnorm_previous        = pseudo->fnorm;
    629:   }
    630:   if (pseudo->fnorm == 0.0) {
    631:     *newdt = 1.e12*inc*ts->time_step;
    632:   } else if (pseudo->increment_dt_from_initial_dt) {
    633:     *newdt = inc*ts->initial_time_step*pseudo->initial_fnorm/pseudo->fnorm;
    634:   } else {
    635:     *newdt = inc*ts->time_step*fnorm_previous/pseudo->fnorm;
    636:   }
    637:   pseudo->fnorm_previous = pseudo->fnorm;
    638:   return(0);
    639: }