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