Actual source code: posindep.c

  1: /*$Id: posindep.c,v 1.60 2001/09/11 16:34:22 bsmith Exp $*/
  2: /*
  3:        Code for Timestepping with implicit backwards Euler.
  4: */
 5:  #include src/ts/tsimpl.h

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

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

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

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

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

 26: /* ------------------------------------------------------------------------------*/

 28: #undef __FUNCT__  
 30: /*@
 31:     TSPseudoComputeTimeStep - Computes the next timestep for a currently running
 32:     pseudo-timestepping process.

 34:     Collective on TS

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

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

 42:     Level: advanced

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

 48: .keywords: timestep, pseudo, compute

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

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


 65: /* ------------------------------------------------------------------------------*/
 66: #undef __FUNCT__  
 68: /*@C
 69:    TSPseudoDefaultVerifyTimeStep - Default code to verify the quality of the last timestep.

 71:    Collective on TS

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

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

 82:    Level: advanced

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

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

 90: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep()
 91: @*/
 92: int TSPseudoDefaultVerifyTimeStep(TS ts,Vec update,void *dtctx,PetscReal *newdt,int *flag)
 93: {
 95:   *flag = 1;
 96:   return(0);
 97: }


100: #undef __FUNCT__  
102: /*@
103:     TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.

105:     Collective on TS

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

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

115:     Level: advanced

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

121: .keywords: timestep, pseudo, verify 

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

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

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

135:   return(0);
136: }

138: /* --------------------------------------------------------------------------------*/

140: #undef __FUNCT__  
142: static int TSStep_Pseudo(TS ts,int *steps,PetscReal *ptime)
143: {
144:   Vec       sol = ts->vec_sol;
145:   int       ierr,i,max_steps = ts->max_steps,its,ok,lits;
146:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
147:   PetscReal current_time_step;
148: 
150:   *steps = -ts->steps;

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

174:   *steps += ts->steps;
175:   *ptime  = ts->ptime;
176:   return(0);
177: }

179: /*------------------------------------------------------------*/
180: #undef __FUNCT__  
182: static int TSDestroy_Pseudo(TS ts)
183: {
184:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
185:   int       ierr;

188:   if (pseudo->update) {VecDestroy(pseudo->update);}
189:   if (pseudo->func) {VecDestroy(pseudo->func);}
190:   if (pseudo->rhs)  {VecDestroy(pseudo->rhs);}
191:   PetscFree(pseudo);
192:   return(0);
193: }


196: /*------------------------------------------------------------*/

198: /* 
199:     This defines the nonlinear equation that is to be solved with SNES

201:               (U^{n+1} - U^{n})/dt - F(U^{n+1})
202: */
203: #undef __FUNCT__  
205: int TSPseudoFunction(SNES snes,Vec x,Vec y,void *ctx)
206: {
207:   TS     ts = (TS) ctx;
208:   PetscScalar mdt = 1.0/ts->time_step,*unp1,*un,*Funp1;
209:   int    ierr,i,n;

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

226:   return(0);
227: }

229: /*
230:    This constructs the Jacobian needed for SNES 

232:              J = I/dt - J_{F}   where J_{F} is the given Jacobian of F.
233: */
234: #undef __FUNCT__  
236: int TSPseudoJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
237: {
238:   TS            ts = (TS) ctx;
239:   int           ierr;
240:   PetscScalar   mone = -1.0,mdt = 1.0/ts->time_step;

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

246:   /* shift and scale Jacobian */
247:   MatScale(&mone,*AA);
248:   MatShift(&mdt,*AA);
249:   if (*BB != *AA && *str != SAME_PRECONDITIONER) {
250:     MatScale(&mone,*BB);
251:     MatShift(&mdt,*BB);
252:   }

254:   return(0);
255: }


258: #undef __FUNCT__  
260: static int TSSetUp_Pseudo(TS ts)
261: {
262:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
263:   int       ierr;

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

275: #undef __FUNCT__  
277: int TSPseudoDefaultMonitor(TS ts,int step,PetscReal ptime,Vec v,void *ctx)
278: {
279:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
280:   int       ierr;

283:   (*PetscHelpPrintf)(ts->comm,"TS %d dt %g time %g fnorm %gn",step,ts->time_step,ptime,pseudo->fnorm);
284:   return(0);
285: }

287: #undef __FUNCT__  
289: static int TSSetFromOptions_Pseudo(TS ts)
290: {
291:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
292:   int        ierr;
293:   PetscTruth flg;


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

309:   return(0);
310: }

312: #undef __FUNCT__  
314: static int TSView_Pseudo(TS ts,PetscViewer viewer)
315: {
317:   return(0);
318: }

320: /* ----------------------------------------------------------------------------- */
321: #undef __FUNCT__  
323: /*@
324:    TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the 
325:    last timestep.

327:    Collective on TS

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

335:    Level: advanced

337:    Calling sequence of func:
338: .  func (TS ts,Vec update,void *ctx,PetscReal *newdt,int *flag);

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

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

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

351: .seealso: TSPseudoDefaultVerifyTimeStep(), TSPseudoVerifyTimeStep()
352: @*/
353: int TSPseudoSetVerifyTimeStep(TS ts,int (*dt)(TS,Vec,void*,PetscReal*,int*),void* ctx)
354: {
355:   int ierr,(*f)(TS,int (*)(TS,Vec,void*,PetscReal *,int *),void *);


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

367: #undef __FUNCT__  
369: /*@
370:     TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to 
371:     dt when using the TSPseudoDefaultTimeStep() routine.

373:    Collective on TS

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

379:     Options Database Key:
380: $    -ts_pseudo_increment <increment>

382:     Level: advanced

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

386: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
387: @*/
388: int TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
389: {
390:   int ierr,(*f)(TS,PetscReal);


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

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

411:    Collective on TS

413:     Input Parameter:
414: .   ts - the timestep context

416:     Options Database Key:
417: $    -ts_pseudo_increment_dt_from_initial_dt

419:     Level: advanced

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

423: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
424: @*/
425: int TSPseudoIncrementDtFromInitialDt(TS ts)
426: {
427:   int ierr,(*f)(TS);


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


440: #undef __FUNCT__  
442: /*@
443:    TSPseudoSetTimeStep - Sets the user-defined routine to be
444:    called at each pseudo-timestep to update the timestep.

446:    Collective on TS

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

454:    Level: intermediate

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

459: .  newdt - the newly computed timestep
460: .  ctx - [optional] timestep context

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

466: .keywords: timestep, pseudo, set

468: .seealso: TSPseudoDefaultTimeStep(), TSPseudoComputeTimeStep()
469: @*/
470: int TSPseudoSetTimeStep(TS ts,int (*dt)(TS,PetscReal*,void*),void* ctx)
471: {
472:   int ierr,(*f)(TS,int (*)(TS,PetscReal *,void *),void *);


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

484: /* ----------------------------------------------------------------------------- */

486: EXTERN_C_BEGIN
487: #undef __FUNCT__  
489: int TSPseudoSetVerifyTimeStep_Pseudo(TS ts,int (*dt)(TS,Vec,void*,PetscReal*,int*),void* ctx)
490: {
491:   TS_Pseudo *pseudo;

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

501: EXTERN_C_BEGIN
502: #undef __FUNCT__  
504: int TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
505: {
506:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

509:   pseudo->dt_increment = inc;
510:   return(0);
511: }
512: EXTERN_C_END

514: EXTERN_C_BEGIN
515: #undef __FUNCT__  
517: int TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
518: {
519:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

522:   pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
523:   return(0);
524: }
525: EXTERN_C_END

527: EXTERN_C_BEGIN
528: #undef __FUNCT__  
530: int TSPseudoSetTimeStep_Pseudo(TS ts,int (*dt)(TS,PetscReal*,void*),void* ctx)
531: {
532:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

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

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

543: EXTERN_C_BEGIN
544: #undef __FUNCT__  
546: int TSCreate_Pseudo(TS ts)
547: {
548:   TS_Pseudo  *pseudo;
549:   int        ierr;

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:   ierr     = 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: }
593: EXTERN_C_END

595: #undef __FUNCT__  
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: int 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;
620:   int       ierr;

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: }