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