Actual source code: posindep.c

  1: #define PETSCTS_DLL

  3: /*
  4:        Code for Timestepping with implicit backwards Euler.
  5: */
 6:  #include include/private/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  TSPseudoComputeTimeStep(TS ts,PetscReal *dt)
 54: {
 55:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

 60:   (*pseudo->dt)(ts,dt,pseudo->dtctx);
 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  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  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);
230:   return(0);
231: }

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

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

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

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


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

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

274: PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
275: {
276:   TS_Pseudo               *pseudo = (TS_Pseudo*)ts->data;
277:   PetscErrorCode          ierr;
278:   PetscViewerASCIIMonitor viewer = (PetscViewerASCIIMonitor)ctx;

281:   if (!ctx) {
282:     PetscViewerASCIIMonitorCreate(ts->comm,"stdout",0,&viewer);
283:   }
284:   PetscViewerASCIIMonitorPrintf(viewer,"TS %D dt %G time %G fnorm %G\n",step,ts->time_step,ptime,pseudo->fnorm);
285:   if (!ctx) {
286:     PetscViewerASCIIMonitorDestroy(viewer);
287:   }
288:   return(0);
289: }

293: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts)
294: {
295:   TS_Pseudo               *pseudo = (TS_Pseudo*)ts->data;
296:   PetscErrorCode          ierr;
297:   PetscTruth              flg;
298:   PetscViewerASCIIMonitor viewer;

301:   PetscOptionsHead("Pseudo-timestepping options");
302:     PetscOptionsName("-ts_monitor","Monitor convergence","TSPseudoMonitorDefault",&flg);
303:     if (flg) {
304:       PetscViewerASCIIMonitorCreate(ts->comm,"stdout",0,&viewer);
305:       TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void*))PetscViewerASCIIMonitorDestroy);
306:     }
307:     PetscOptionsName("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",&flg);
308:     if (flg) {
309:       TSPseudoIncrementDtFromInitialDt(ts);
310:     }
311:     PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,0);
312:   PetscOptionsTail();
313:   return(0);
314: }

318: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
319: {
321:   return(0);
322: }

324: /* ----------------------------------------------------------------------------- */
327: /*@C
328:    TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the 
329:    last timestep.

331:    Collective on TS

333:    Input Parameters:
334: +  ts - timestep context
335: .  dt - user-defined function to verify timestep
336: -  ctx - [optional] user-defined context for private data
337:          for the timestep verification routine (may be PETSC_NULL)

339:    Level: advanced

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

344: .  update - latest solution vector
345: .  ctx - [optional] timestep context
346: .  newdt - the timestep to use for the next step
347: .  flag - flag indicating whether the last time step was acceptable

349:    Notes:
350:    The routine set here will be called by TSPseudoVerifyTimeStep()
351:    during the timestepping process.

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

355: .seealso: TSPseudoDefaultVerifyTimeStep(), TSPseudoVerifyTimeStep()
356: @*/
357: PetscErrorCode  TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscTruth*),void* ctx)
358: {
359:   PetscErrorCode ierr,(*f)(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal *,PetscTruth *),void *);


364:   PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",(void (**)(void))&f);
365:   if (f) {
366:     (*f)(ts,dt,ctx);
367:   }
368:   return(0);
369: }

373: /*@
374:     TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to 
375:     dt when using the TSPseudoDefaultTimeStep() routine.

377:    Collective on TS

379:     Input Parameters:
380: +   ts - the timestep context
381: -   inc - the scaling factor >= 1.0

383:     Options Database Key:
384: $    -ts_pseudo_increment <increment>

386:     Level: advanced

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

390: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
391: @*/
392: PetscErrorCode  TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
393: {
394:   PetscErrorCode ierr,(*f)(TS,PetscReal);


399:   PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",(void (**)(void))&f);
400:   if (f) {
401:     (*f)(ts,inc);
402:   }
403:   return(0);
404: }

408: /*@
409:     TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
410:     is computed via the formula
411: $         dt = initial_dt*initial_fnorm/current_fnorm 
412:       rather than the default update,
413: $         dt = current_dt*previous_fnorm/current_fnorm.

415:    Collective on TS

417:     Input Parameter:
418: .   ts - the timestep context

420:     Options Database Key:
421: $    -ts_pseudo_increment_dt_from_initial_dt

423:     Level: advanced

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

427: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
428: @*/
429: PetscErrorCode  TSPseudoIncrementDtFromInitialDt(TS ts)
430: {
431:   PetscErrorCode ierr,(*f)(TS);


436:   PetscObjectQueryFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",(void (**)(void))&f);
437:   if (f) {
438:     (*f)(ts);
439:   }
440:   return(0);
441: }


446: /*@C
447:    TSPseudoSetTimeStep - Sets the user-defined routine to be
448:    called at each pseudo-timestep to update the timestep.

450:    Collective on TS

452:    Input Parameters:
453: +  ts - timestep context
454: .  dt - function to compute timestep
455: -  ctx - [optional] user-defined context for private data
456:          required by the function (may be PETSC_NULL)

458:    Level: intermediate

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

463: .  newdt - the newly computed timestep
464: .  ctx - [optional] timestep context

466:    Notes:
467:    The routine set here will be called by TSPseudoComputeTimeStep()
468:    during the timestepping process.

470: .keywords: timestep, pseudo, set

472: .seealso: TSPseudoDefaultTimeStep(), TSPseudoComputeTimeStep()
473: @*/
474: PetscErrorCode  TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void* ctx)
475: {
476:   PetscErrorCode ierr,(*f)(TS,PetscErrorCode (*)(TS,PetscReal *,void *),void *);


481:   PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",(void (**)(void))&f);
482:   if (f) {
483:     (*f)(ts,dt,ctx);
484:   }
485:   return(0);
486: }

488: /* ----------------------------------------------------------------------------- */

494: PetscErrorCode  TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void* ctx)
495: {
496:   TS_Pseudo *pseudo;

499:   pseudo              = (TS_Pseudo*)ts->data;
500:   pseudo->verify      = dt;
501:   pseudo->verifyctx   = ctx;
502:   return(0);
503: }

509: PetscErrorCode  TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
510: {
511:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

514:   pseudo->dt_increment = inc;
515:   return(0);
516: }

522: PetscErrorCode  TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
523: {
524:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

527:   pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
528:   return(0);
529: }

536: PetscErrorCode  TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void* ctx)
537: {
538:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

541:   pseudo->dt      = dt;
542:   pseudo->dtctx   = ctx;
543:   return(0);
544: }

547: /* ----------------------------------------------------------------------------- */

552: PetscErrorCode  TSCreate_Pseudo(TS ts)
553: {
554:   TS_Pseudo      *pseudo;

558:   ts->ops->destroy         = TSDestroy_Pseudo;
559:   ts->ops->view            = TSView_Pseudo;

561:   if (ts->problem_type == TS_LINEAR) {
562:     SETERRQ(PETSC_ERR_ARG_WRONG,"Only for nonlinear problems");
563:   }
564:   if (!ts->Arhs) {
565:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set Jacobian");
566:   }

568:   ts->ops->setup           = TSSetUp_Pseudo;
569:   ts->ops->step            = TSStep_Pseudo;
570:   ts->ops->setfromoptions  = TSSetFromOptions_Pseudo;

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

575:   PetscNew(TS_Pseudo,&pseudo);
576:   PetscLogObjectMemory(ts,sizeof(TS_Pseudo));

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

581:   pseudo->dt_increment                 = 1.1;
582:   pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
583:   pseudo->dt                           = TSPseudoDefaultTimeStep;

585:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",
586:                     "TSPseudoSetVerifyTimeStep_Pseudo",
587:                      TSPseudoSetVerifyTimeStep_Pseudo);
588:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",
589:                     "TSPseudoSetTimeStepIncrement_Pseudo",
590:                      TSPseudoSetTimeStepIncrement_Pseudo);
591:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",
592:                     "TSPseudoIncrementDtFromInitialDt_Pseudo",
593:                      TSPseudoIncrementDtFromInitialDt_Pseudo);
594:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStep_C",
595:                     "TSPseudoSetTimeStep_Pseudo",
596:                      TSPseudoSetTimeStep_Pseudo);
597:   return(0);
598: }

603: /*@C
604:    TSPseudoDefaultTimeStep - Default code to compute pseudo-timestepping.
605:    Use with TSPseudoSetTimeStep().

607:    Collective on TS

609:    Input Parameters:
610: .  ts - the timestep context
611: .  dtctx - unused timestep context

613:    Output Parameter:
614: .  newdt - the timestep to use for the next step

616:    Level: advanced

618: .keywords: timestep, pseudo, default

620: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
621: @*/
622: PetscErrorCode  TSPseudoDefaultTimeStep(TS ts,PetscReal* newdt,void* dtctx)
623: {
624:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
625:   PetscReal      inc = pseudo->dt_increment,fnorm_previous = pseudo->fnorm_previous;

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