1 /*
2 Code for Timestepping with implicit backwards Euler.
3 */
4 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
5
6 typedef struct {
7 Vec update; /* work vector where new solution is formed */
8 Vec func; /* work vector where F(t[i],u[i]) is stored */
9 Vec xdot; /* work vector for time derivative of state */
10
11 /* information used for Pseudo-timestepping */
12
13 PetscErrorCode (*dt)(TS,PetscReal*,void*); /* compute next timestep, and related context */
14 void *dtctx;
15 PetscErrorCode (*verify)(TS,Vec,void*,PetscReal*,PetscBool*); /* verify previous timestep and related context */
16 void *verifyctx;
17
18 PetscReal fnorm_initial,fnorm; /* original and current norm of F(u) */
19 PetscReal fnorm_previous;
20
21 PetscReal dt_initial; /* initial time-step */
22 PetscReal dt_increment; /* scaling that dt is incremented each time-step */
23 PetscReal dt_max; /* maximum time step */
24 PetscBool increment_dt_from_initial_dt;
25 PetscReal fatol,frtol;
26 } TS_Pseudo;
27
28 /* ------------------------------------------------------------------------------*/
29
30 /*@C
31 TSPseudoComputeTimeStep - Computes the next timestep for a currently running
32 pseudo-timestepping process.
33
34 Collective on TS
35
36 Input Parameter:
37 . ts - timestep context
38
39 Output Parameter:
40 . dt - newly computed timestep
41
42 Level: developer
43
44 Notes:
45 The routine to be called here to compute the timestep should be
46 set by calling TSPseudoSetTimeStep().
47
48 .seealso: TSPseudoTimeStepDefault(), TSPseudoSetTimeStep()
49 @*/
TSPseudoComputeTimeStep(TS ts,PetscReal * dt)50 PetscErrorCode TSPseudoComputeTimeStep(TS ts,PetscReal *dt)
51 {
52 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
53 PetscErrorCode ierr;
54
55 PetscFunctionBegin;
56 ierr = PetscLogEventBegin(TS_PseudoComputeTimeStep,ts,0,0,0);CHKERRQ(ierr);
57 ierr = (*pseudo->dt)(ts,dt,pseudo->dtctx);CHKERRQ(ierr);
58 ierr = PetscLogEventEnd(TS_PseudoComputeTimeStep,ts,0,0,0);CHKERRQ(ierr);
59 PetscFunctionReturn(0);
60 }
61
62
63 /* ------------------------------------------------------------------------------*/
64 /*@C
65 TSPseudoVerifyTimeStepDefault - Default code to verify the quality of the last timestep.
66
67 Collective on TS
68
69 Input Parameters:
70 + ts - the timestep context
71 . dtctx - unused timestep context
72 - update - latest solution vector
73
74 Output Parameters:
75 + newdt - the timestep to use for the next step
76 - flag - flag indicating whether the last time step was acceptable
77
78 Level: advanced
79
80 Note:
81 This routine always returns a flag of 1, indicating an acceptable
82 timestep.
83
84 .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep()
85 @*/
TSPseudoVerifyTimeStepDefault(TS ts,Vec update,void * dtctx,PetscReal * newdt,PetscBool * flag)86 PetscErrorCode TSPseudoVerifyTimeStepDefault(TS ts,Vec update,void *dtctx,PetscReal *newdt,PetscBool *flag)
87 {
88 PetscFunctionBegin;
89 *flag = PETSC_TRUE;
90 PetscFunctionReturn(0);
91 }
92
93
94 /*@
95 TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.
96
97 Collective on TS
98
99 Input Parameters:
100 + ts - timestep context
101 - update - latest solution vector
102
103 Output Parameters:
104 + dt - newly computed timestep (if it had to shrink)
105 - flag - indicates if current timestep was ok
106
107 Level: advanced
108
109 Notes:
110 The routine to be called here to compute the timestep should be
111 set by calling TSPseudoSetVerifyTimeStep().
112
113 .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStepDefault()
114 @*/
TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal * dt,PetscBool * flag)115 PetscErrorCode TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscBool *flag)
116 {
117 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
118 PetscErrorCode ierr;
119
120 PetscFunctionBegin;
121 *flag = PETSC_TRUE;
122 if (pseudo->verify) {
123 ierr = (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);CHKERRQ(ierr);
124 }
125 PetscFunctionReturn(0);
126 }
127
128 /* --------------------------------------------------------------------------------*/
129
TSStep_Pseudo(TS ts)130 static PetscErrorCode TSStep_Pseudo(TS ts)
131 {
132 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
133 PetscInt nits,lits,reject;
134 PetscBool stepok;
135 PetscReal next_time_step = ts->time_step;
136 PetscErrorCode ierr;
137
138 PetscFunctionBegin;
139 if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
140 ierr = VecCopy(ts->vec_sol,pseudo->update);CHKERRQ(ierr);
141 ierr = TSPseudoComputeTimeStep(ts,&next_time_step);CHKERRQ(ierr);
142 for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
143 ts->time_step = next_time_step;
144 ierr = TSPreStage(ts,ts->ptime+ts->time_step);CHKERRQ(ierr);
145 ierr = SNESSolve(ts->snes,NULL,pseudo->update);CHKERRQ(ierr);
146 ierr = SNESGetIterationNumber(ts->snes,&nits);CHKERRQ(ierr);
147 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr);
148 ts->snes_its += nits; ts->ksp_its += lits;
149 ierr = TSPostStage(ts,ts->ptime+ts->time_step,0,&(pseudo->update));CHKERRQ(ierr);
150 ierr = TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,pseudo->update,&stepok);CHKERRQ(ierr);
151 if (!stepok) {next_time_step = ts->time_step; continue;}
152 pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
153 ierr = TSPseudoVerifyTimeStep(ts,pseudo->update,&next_time_step,&stepok);CHKERRQ(ierr);
154 if (stepok) break;
155 }
156 if (reject >= ts->max_reject) {
157 ts->reason = TS_DIVERGED_STEP_REJECTED;
158 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);CHKERRQ(ierr);
159 PetscFunctionReturn(0);
160 }
161
162 ierr = VecCopy(pseudo->update,ts->vec_sol);CHKERRQ(ierr);
163 ts->ptime += ts->time_step;
164 ts->time_step = next_time_step;
165
166 if (pseudo->fnorm < 0) {
167 ierr = VecZeroEntries(pseudo->xdot);CHKERRQ(ierr);
168 ierr = TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);CHKERRQ(ierr);
169 ierr = VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);CHKERRQ(ierr);
170 }
171 if (pseudo->fnorm < pseudo->fatol) {
172 ts->reason = TS_CONVERGED_PSEUDO_FATOL;
173 ierr = PetscInfo3(ts,"Step=%D, converged since fnorm %g < fatol %g\n",ts->steps,pseudo->fnorm,pseudo->frtol);CHKERRQ(ierr);
174 PetscFunctionReturn(0);
175 }
176 if (pseudo->fnorm/pseudo->fnorm_initial < pseudo->frtol) {
177 ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
178 ierr = PetscInfo4(ts,"Step=%D, converged since fnorm %g / fnorm_initial %g < frtol %g\n",ts->steps,pseudo->fnorm,pseudo->fnorm_initial,pseudo->fatol);CHKERRQ(ierr);
179 PetscFunctionReturn(0);
180 }
181 PetscFunctionReturn(0);
182 }
183
184 /*------------------------------------------------------------*/
TSReset_Pseudo(TS ts)185 static PetscErrorCode TSReset_Pseudo(TS ts)
186 {
187 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
188 PetscErrorCode ierr;
189
190 PetscFunctionBegin;
191 ierr = VecDestroy(&pseudo->update);CHKERRQ(ierr);
192 ierr = VecDestroy(&pseudo->func);CHKERRQ(ierr);
193 ierr = VecDestroy(&pseudo->xdot);CHKERRQ(ierr);
194 PetscFunctionReturn(0);
195 }
196
TSDestroy_Pseudo(TS ts)197 static PetscErrorCode TSDestroy_Pseudo(TS ts)
198 {
199 PetscErrorCode ierr;
200
201 PetscFunctionBegin;
202 ierr = TSReset_Pseudo(ts);CHKERRQ(ierr);
203 ierr = PetscFree(ts->data);CHKERRQ(ierr);
204 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",NULL);CHKERRQ(ierr);
205 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",NULL);CHKERRQ(ierr);
206 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",NULL);CHKERRQ(ierr);
207 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",NULL);CHKERRQ(ierr);
208 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",NULL);CHKERRQ(ierr);
209 PetscFunctionReturn(0);
210 }
211
212 /*------------------------------------------------------------*/
213
214 /*
215 Compute Xdot = (X^{n+1}-X^n)/dt) = 0
216 */
TSPseudoGetXdot(TS ts,Vec X,Vec * Xdot)217 static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
218 {
219 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
220 const PetscScalar mdt = 1.0/ts->time_step,*xnp1,*xn;
221 PetscScalar *xdot;
222 PetscErrorCode ierr;
223 PetscInt i,n;
224
225 PetscFunctionBegin;
226 *Xdot = NULL;
227 ierr = VecGetArrayRead(ts->vec_sol,&xn);CHKERRQ(ierr);
228 ierr = VecGetArrayRead(X,&xnp1);CHKERRQ(ierr);
229 ierr = VecGetArray(pseudo->xdot,&xdot);CHKERRQ(ierr);
230 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr);
231 for (i=0; i<n; i++) xdot[i] = mdt*(xnp1[i] - xn[i]);
232 ierr = VecRestoreArrayRead(ts->vec_sol,&xn);CHKERRQ(ierr);
233 ierr = VecRestoreArrayRead(X,&xnp1);CHKERRQ(ierr);
234 ierr = VecRestoreArray(pseudo->xdot,&xdot);CHKERRQ(ierr);
235 *Xdot = pseudo->xdot;
236 PetscFunctionReturn(0);
237 }
238
239 /*
240 The transient residual is
241
242 F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
243
244 or for ODE,
245
246 (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
247
248 This is the function that must be evaluated for transient simulation and for
249 finite difference Jacobians. On the first Newton step, this algorithm uses
250 a guess of U^{n+1} = U^n in which case the transient term vanishes and the
251 residual is actually the steady state residual. Pseudotransient
252 continuation as described in the literature is a linearly implicit
253 algorithm, it only takes this one Newton step with the steady state
254 residual, and then advances to the next time step.
255 */
SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)256 static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)
257 {
258 Vec Xdot;
259 PetscErrorCode ierr;
260
261 PetscFunctionBegin;
262 ierr = TSPseudoGetXdot(ts,X,&Xdot);CHKERRQ(ierr);
263 ierr = TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);CHKERRQ(ierr);
264 PetscFunctionReturn(0);
265 }
266
267 /*
268 This constructs the Jacobian needed for SNES. For DAE, this is
269
270 dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
271
272 and for ODE:
273
274 J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
275 */
SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)276 static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)
277 {
278 Vec Xdot;
279 PetscErrorCode ierr;
280
281 PetscFunctionBegin;
282 ierr = TSPseudoGetXdot(ts,X,&Xdot);CHKERRQ(ierr);
283 ierr = TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,PETSC_FALSE);CHKERRQ(ierr);
284 PetscFunctionReturn(0);
285 }
286
287
TSSetUp_Pseudo(TS ts)288 static PetscErrorCode TSSetUp_Pseudo(TS ts)
289 {
290 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
291 PetscErrorCode ierr;
292
293 PetscFunctionBegin;
294 ierr = VecDuplicate(ts->vec_sol,&pseudo->update);CHKERRQ(ierr);
295 ierr = VecDuplicate(ts->vec_sol,&pseudo->func);CHKERRQ(ierr);
296 ierr = VecDuplicate(ts->vec_sol,&pseudo->xdot);CHKERRQ(ierr);
297 PetscFunctionReturn(0);
298 }
299 /*------------------------------------------------------------*/
300
TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void * dummy)301 static PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
302 {
303 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
304 PetscErrorCode ierr;
305 PetscViewer viewer = (PetscViewer) dummy;
306
307 PetscFunctionBegin;
308 if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
309 ierr = VecZeroEntries(pseudo->xdot);CHKERRQ(ierr);
310 ierr = TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);CHKERRQ(ierr);
311 ierr = VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);CHKERRQ(ierr);
312 }
313 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr);
314 ierr = PetscViewerASCIIPrintf(viewer,"TS %D dt %g time %g fnorm %g\n",step,(double)ts->time_step,(double)ptime,(double)pseudo->fnorm);CHKERRQ(ierr);
315 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr);
316 PetscFunctionReturn(0);
317 }
318
TSSetFromOptions_Pseudo(PetscOptionItems * PetscOptionsObject,TS ts)319 static PetscErrorCode TSSetFromOptions_Pseudo(PetscOptionItems *PetscOptionsObject,TS ts)
320 {
321 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
322 PetscErrorCode ierr;
323 PetscBool flg = PETSC_FALSE;
324 PetscViewer viewer;
325
326 PetscFunctionBegin;
327 ierr = PetscOptionsHead(PetscOptionsObject,"Pseudo-timestepping options");CHKERRQ(ierr);
328 ierr = PetscOptionsBool("-ts_monitor_pseudo","Monitor convergence","",flg,&flg,NULL);CHKERRQ(ierr);
329 if (flg) {
330 ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),"stdout",&viewer);CHKERRQ(ierr);
331 ierr = TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr);
332 }
333 flg = pseudo->increment_dt_from_initial_dt;
334 ierr = PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",flg,&flg,NULL);CHKERRQ(ierr);
335 pseudo->increment_dt_from_initial_dt = flg;
336 ierr = PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,NULL);CHKERRQ(ierr);
337 ierr = PetscOptionsReal("-ts_pseudo_max_dt","Maximum value for dt","TSPseudoSetMaxTimeStep",pseudo->dt_max,&pseudo->dt_max,NULL);CHKERRQ(ierr);
338 ierr = PetscOptionsReal("-ts_pseudo_fatol","Tolerance for norm of function","",pseudo->fatol,&pseudo->fatol,NULL);CHKERRQ(ierr);
339 ierr = PetscOptionsReal("-ts_pseudo_frtol","Relative tolerance for norm of function","",pseudo->frtol,&pseudo->frtol,NULL);CHKERRQ(ierr);
340 ierr = PetscOptionsTail();CHKERRQ(ierr);
341 PetscFunctionReturn(0);
342 }
343
TSView_Pseudo(TS ts,PetscViewer viewer)344 static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
345 {
346 PetscErrorCode ierr;
347 PetscBool isascii;
348
349 PetscFunctionBegin;
350 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
351 if (isascii) {
352 TS_Pseudo *pseudo = (TS_Pseudo*) ts->data;
353 ierr = PetscViewerASCIIPrintf(viewer," frtol - relative tolerance in function value %g\n",(double)pseudo->frtol);CHKERRQ(ierr);
354 ierr = PetscViewerASCIIPrintf(viewer," fatol - absolute tolerance in function value %g\n",(double)pseudo->fatol);CHKERRQ(ierr);
355 ierr = PetscViewerASCIIPrintf(viewer," dt_initial - initial timestep %g\n",(double)pseudo->dt_initial);CHKERRQ(ierr);
356 ierr = PetscViewerASCIIPrintf(viewer," dt_increment - increase in timestep on successful step %g\n",(double)pseudo->dt_increment);CHKERRQ(ierr);
357 ierr = PetscViewerASCIIPrintf(viewer," dt_max - maximum time %g\n",(double)pseudo->dt_max);CHKERRQ(ierr);
358 }
359 PetscFunctionReturn(0);
360 }
361
362 /* ----------------------------------------------------------------------------- */
363 /*@C
364 TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
365 last timestep.
366
367 Logically Collective on TS
368
369 Input Parameters:
370 + ts - timestep context
371 . dt - user-defined function to verify timestep
372 - ctx - [optional] user-defined context for private data
373 for the timestep verification routine (may be NULL)
374
375 Level: advanced
376
377 Calling sequence of func:
378 $ func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscBool *flag);
379
380 + update - latest solution vector
381 . ctx - [optional] timestep context
382 . newdt - the timestep to use for the next step
383 - flag - flag indicating whether the last time step was acceptable
384
385 Notes:
386 The routine set here will be called by TSPseudoVerifyTimeStep()
387 during the timestepping process.
388
389 .seealso: TSPseudoVerifyTimeStepDefault(), TSPseudoVerifyTimeStep()
390 @*/
TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (* dt)(TS,Vec,void *,PetscReal *,PetscBool *),void * ctx)391 PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool*),void *ctx)
392 {
393 PetscErrorCode ierr;
394
395 PetscFunctionBegin;
396 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
397 ierr = PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal*,PetscBool*),void*),(ts,dt,ctx));CHKERRQ(ierr);
398 PetscFunctionReturn(0);
399 }
400
401 /*@
402 TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
403 dt when using the TSPseudoTimeStepDefault() routine.
404
405 Logically Collective on TS
406
407 Input Parameters:
408 + ts - the timestep context
409 - inc - the scaling factor >= 1.0
410
411 Options Database Key:
412 . -ts_pseudo_increment <increment>
413
414 Level: advanced
415
416 .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
417 @*/
TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)418 PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
419 {
420 PetscErrorCode ierr;
421
422 PetscFunctionBegin;
423 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
424 PetscValidLogicalCollectiveReal(ts,inc,2);
425 ierr = PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));CHKERRQ(ierr);
426 PetscFunctionReturn(0);
427 }
428
429 /*@
430 TSPseudoSetMaxTimeStep - Sets the maximum time step
431 when using the TSPseudoTimeStepDefault() routine.
432
433 Logically Collective on TS
434
435 Input Parameters:
436 + ts - the timestep context
437 - maxdt - the maximum time step, use a non-positive value to deactivate
438
439 Options Database Key:
440 . -ts_pseudo_max_dt <increment>
441
442 Level: advanced
443
444 .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
445 @*/
TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)446 PetscErrorCode TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
447 {
448 PetscErrorCode ierr;
449
450 PetscFunctionBegin;
451 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
452 PetscValidLogicalCollectiveReal(ts,maxdt,2);
453 ierr = PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));CHKERRQ(ierr);
454 PetscFunctionReturn(0);
455 }
456
457 /*@
458 TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
459 is computed via the formula
460 $ dt = initial_dt*initial_fnorm/current_fnorm
461 rather than the default update,
462 $ dt = current_dt*previous_fnorm/current_fnorm.
463
464 Logically Collective on TS
465
466 Input Parameter:
467 . ts - the timestep context
468
469 Options Database Key:
470 . -ts_pseudo_increment_dt_from_initial_dt
471
472 Level: advanced
473
474 .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
475 @*/
TSPseudoIncrementDtFromInitialDt(TS ts)476 PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
477 {
478 PetscErrorCode ierr;
479
480 PetscFunctionBegin;
481 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
482 ierr = PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));CHKERRQ(ierr);
483 PetscFunctionReturn(0);
484 }
485
486
487 /*@C
488 TSPseudoSetTimeStep - Sets the user-defined routine to be
489 called at each pseudo-timestep to update the timestep.
490
491 Logically Collective on TS
492
493 Input Parameters:
494 + ts - timestep context
495 . dt - function to compute timestep
496 - ctx - [optional] user-defined context for private data
497 required by the function (may be NULL)
498
499 Level: intermediate
500
501 Calling sequence of func:
502 $ func (TS ts,PetscReal *newdt,void *ctx);
503
504 + newdt - the newly computed timestep
505 - ctx - [optional] timestep context
506
507 Notes:
508 The routine set here will be called by TSPseudoComputeTimeStep()
509 during the timestepping process.
510 If not set then TSPseudoTimeStepDefault() is automatically used
511
512 .seealso: TSPseudoTimeStepDefault(), TSPseudoComputeTimeStep()
513 @*/
TSPseudoSetTimeStep(TS ts,PetscErrorCode (* dt)(TS,PetscReal *,void *),void * ctx)514 PetscErrorCode TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void *ctx)
515 {
516 PetscErrorCode ierr;
517
518 PetscFunctionBegin;
519 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
520 ierr = PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal*,void*),void*),(ts,dt,ctx));CHKERRQ(ierr);
521 PetscFunctionReturn(0);
522 }
523
524 /* ----------------------------------------------------------------------------- */
525
526 typedef PetscErrorCode (*FCN1)(TS,Vec,void*,PetscReal*,PetscBool*); /* force argument to next function to not be extern C*/
TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void * ctx)527 static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void *ctx)
528 {
529 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
530
531 PetscFunctionBegin;
532 pseudo->verify = dt;
533 pseudo->verifyctx = ctx;
534 PetscFunctionReturn(0);
535 }
536
TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)537 static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
538 {
539 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
540
541 PetscFunctionBegin;
542 pseudo->dt_increment = inc;
543 PetscFunctionReturn(0);
544 }
545
TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)546 static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
547 {
548 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
549
550 PetscFunctionBegin;
551 pseudo->dt_max = maxdt;
552 PetscFunctionReturn(0);
553 }
554
TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)555 static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
556 {
557 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
558
559 PetscFunctionBegin;
560 pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
561 PetscFunctionReturn(0);
562 }
563
564 typedef PetscErrorCode (*FCN2)(TS,PetscReal*,void*); /* force argument to next function to not be extern C*/
TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void * ctx)565 static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void *ctx)
566 {
567 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
568
569 PetscFunctionBegin;
570 pseudo->dt = dt;
571 pseudo->dtctx = ctx;
572 PetscFunctionReturn(0);
573 }
574
575 /* ----------------------------------------------------------------------------- */
576 /*MC
577 TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping
578
579 This method solves equations of the form
580
581 $ F(X,Xdot) = 0
582
583 for steady state using the iteration
584
585 $ [G'] S = -F(X,0)
586 $ X += S
587
588 where
589
590 $ G(Y) = F(Y,(Y-X)/dt)
591
592 This is linearly-implicit Euler with the residual always evaluated "at steady
593 state". See note below.
594
595 Options database keys:
596 + -ts_pseudo_increment <real> - ratio of increase dt
597 . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
598 . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
599 - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol
600
601 Level: beginner
602
603 References:
604 + 1. - Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential Algebraic Equations, 2003.
605 - 2. - C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.
606
607 Notes:
608 The residual computed by this method includes the transient term (Xdot is computed instead of
609 always being zero), but since the prediction from the last step is always the solution from the
610 last step, on the first Newton iteration we have
611
612 $ Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
613
614 Therefore, the linear system solved by the first Newton iteration is equivalent to the one
615 described above and in the papers. If the user chooses to perform multiple Newton iterations, the
616 algorithm is no longer the one described in the referenced papers.
617
618 .seealso: TSCreate(), TS, TSSetType()
619
620 M*/
TSCreate_Pseudo(TS ts)621 PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
622 {
623 TS_Pseudo *pseudo;
624 PetscErrorCode ierr;
625 SNES snes;
626 SNESType stype;
627
628 PetscFunctionBegin;
629 ts->ops->reset = TSReset_Pseudo;
630 ts->ops->destroy = TSDestroy_Pseudo;
631 ts->ops->view = TSView_Pseudo;
632 ts->ops->setup = TSSetUp_Pseudo;
633 ts->ops->step = TSStep_Pseudo;
634 ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
635 ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
636 ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
637 ts->default_adapt_type = TSADAPTNONE;
638 ts->usessnes = PETSC_TRUE;
639
640 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
641 ierr = SNESGetType(snes,&stype);CHKERRQ(ierr);
642 if (!stype) {ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr);}
643
644 ierr = PetscNewLog(ts,&pseudo);CHKERRQ(ierr);
645 ts->data = (void*)pseudo;
646
647 pseudo->dt = TSPseudoTimeStepDefault;
648 pseudo->dtctx = NULL;
649 pseudo->dt_increment = 1.1;
650 pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
651 pseudo->fnorm = -1;
652 pseudo->fnorm_initial = -1;
653 pseudo->fnorm_previous = -1;
654 #if defined(PETSC_USE_REAL_SINGLE)
655 pseudo->fatol = 1.e-25;
656 pseudo->frtol = 1.e-5;
657 #else
658 pseudo->fatol = 1.e-50;
659 pseudo->frtol = 1.e-12;
660 #endif
661 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",TSPseudoSetVerifyTimeStep_Pseudo);CHKERRQ(ierr);
662 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",TSPseudoSetTimeStepIncrement_Pseudo);CHKERRQ(ierr);
663 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",TSPseudoSetMaxTimeStep_Pseudo);CHKERRQ(ierr);
664 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",TSPseudoIncrementDtFromInitialDt_Pseudo);CHKERRQ(ierr);
665 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",TSPseudoSetTimeStep_Pseudo);CHKERRQ(ierr);
666 PetscFunctionReturn(0);
667 }
668
669 /*@C
670 TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
671 Use with TSPseudoSetTimeStep().
672
673 Collective on TS
674
675 Input Parameters:
676 + ts - the timestep context
677 - dtctx - unused timestep context
678
679 Output Parameter:
680 . newdt - the timestep to use for the next step
681
682 Level: advanced
683
684 .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
685 @*/
TSPseudoTimeStepDefault(TS ts,PetscReal * newdt,void * dtctx)686 PetscErrorCode TSPseudoTimeStepDefault(TS ts,PetscReal *newdt,void *dtctx)
687 {
688 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
689 PetscReal inc = pseudo->dt_increment;
690 PetscErrorCode ierr;
691
692 PetscFunctionBegin;
693 ierr = VecZeroEntries(pseudo->xdot);CHKERRQ(ierr);
694 ierr = TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);CHKERRQ(ierr);
695 ierr = VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);CHKERRQ(ierr);
696 if (pseudo->fnorm_initial < 0) {
697 /* first time through so compute initial function norm */
698 pseudo->fnorm_initial = pseudo->fnorm;
699 pseudo->fnorm_previous = pseudo->fnorm;
700 }
701 if (pseudo->fnorm == 0.0) *newdt = 1.e12*inc*ts->time_step;
702 else if (pseudo->increment_dt_from_initial_dt) *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
703 else *newdt = inc*ts->time_step*pseudo->fnorm_previous/pseudo->fnorm;
704 if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
705 pseudo->fnorm_previous = pseudo->fnorm;
706 PetscFunctionReturn(0);
707 }
708