Actual source code: tsadapt.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
4: static PetscFunctionList TSAdaptList;
5: static PetscBool TSAdaptPackageInitialized;
6: static PetscBool TSAdaptRegisterAllCalled;
7: static PetscClassId TSADAPT_CLASSID;
9: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt);
10: PETSC_EXTERN PetscErrorCode TSAdaptCreate_None(TSAdapt);
11: PETSC_EXTERN PetscErrorCode TSAdaptCreate_CFL(TSAdapt);
15: /*@C
16: TSAdaptRegister - adds a TSAdapt implementation
18: Not Collective
20: Input Parameters:
21: + name_scheme - name of user-defined adaptivity scheme
22: - routine_create - routine to create method context
24: Notes:
25: TSAdaptRegister() may be called multiple times to add several user-defined families.
27: Sample usage:
28: .vb
29: TSAdaptRegister("my_scheme",MySchemeCreate);
30: .ve
32: Then, your scheme can be chosen with the procedural interface via
33: $ TSAdaptSetType(ts,"my_scheme")
34: or at runtime via the option
35: $ -ts_adapt_type my_scheme
37: Level: advanced
39: .keywords: TSAdapt, register
41: .seealso: TSAdaptRegisterAll()
42: @*/
43: PetscErrorCode TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
44: {
48: PetscFunctionListAdd(&TSAdaptList,sname,function);
49: return(0);
50: }
54: /*@C
55: TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt
57: Not Collective
59: Level: advanced
61: .keywords: TSAdapt, register, all
63: .seealso: TSAdaptRegisterDestroy()
64: @*/
65: PetscErrorCode TSAdaptRegisterAll(void)
66: {
70: TSAdaptRegister(TSADAPTBASIC,TSAdaptCreate_Basic);
71: TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);
72: TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);
73: return(0);
74: }
78: /*@C
79: TSFinalizePackage - This function destroys everything in the TS package. It is
80: called from PetscFinalize().
82: Level: developer
84: .keywords: Petsc, destroy, package
85: .seealso: PetscFinalize()
86: @*/
87: PetscErrorCode TSAdaptFinalizePackage(void)
88: {
92: PetscFunctionListDestroy(&TSAdaptList);
93: TSAdaptPackageInitialized = PETSC_FALSE;
94: TSAdaptRegisterAllCalled = PETSC_FALSE;
95: return(0);
96: }
100: /*@C
101: TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
102: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
103: TSCreate_GL() when using static libraries.
105: Level: developer
107: .keywords: TSAdapt, initialize, package
108: .seealso: PetscInitialize()
109: @*/
110: PetscErrorCode TSAdaptInitializePackage(void)
111: {
115: if (TSAdaptPackageInitialized) return(0);
116: TSAdaptPackageInitialized = PETSC_TRUE;
117: PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
118: TSAdaptRegisterAll();
119: PetscRegisterFinalize(TSAdaptFinalizePackage);
120: return(0);
121: }
125: PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
126: {
127: PetscBool match;
128: PetscErrorCode (*checkstage)(TSAdapt,TS,PetscBool*);
129: PetscErrorCode ierr,(*r)(TSAdapt);
133: PetscObjectTypeCompare((PetscObject)adapt,type,&match);
134: if (match) return(0);
135: PetscFunctionListFind(TSAdaptList,type,&r);
136: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
137: if (adapt->ops->destroy) {(*adapt->ops->destroy)(adapt);}
138: /* Reinitialize function pointers in TSAdaptOps structure */
139: checkstage = adapt->ops->checkstage;
140: PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));
141: adapt->ops->checkstage = checkstage;
142: (*r)(adapt);
143: PetscObjectChangeTypeName((PetscObject)adapt,type);
144: return(0);
145: }
149: PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
150: {
155: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
156: return(0);
157: }
161: /*@C
162: TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView().
164: Collective on PetscViewer
166: Input Parameters:
167: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
168: some related function before a call to TSAdaptLoad().
169: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
170: HDF5 file viewer, obtained from PetscViewerHDF5Open()
172: Level: intermediate
174: Notes:
175: The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.
177: Notes for advanced users:
178: Most users should not need to know the details of the binary storage
179: format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
180: But for anyone who's interested, the standard binary matrix storage
181: format is
182: .vb
183: has not yet been determined
184: .ve
186: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
187: @*/
188: PetscErrorCode TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
189: {
191: PetscBool isbinary;
192: char type[256];
197: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
198: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
200: PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
201: TSAdaptSetType(adapt, type);
202: if (adapt->ops->load) {
203: (*adapt->ops->load)(adapt,viewer);
204: }
205: return(0);
206: }
210: PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer)
211: {
213: PetscBool iascii,isbinary;
217: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);}
220: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
221: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
222: if (iascii) {
223: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
224: PetscViewerASCIIPrintf(viewer," number of candidates %D\n",adapt->candidates.n);
225: if (adapt->ops->view) {
226: PetscViewerASCIIPushTab(viewer);
227: (*adapt->ops->view)(adapt,viewer);
228: PetscViewerASCIIPopTab(viewer);
229: }
230: } else if (isbinary) {
231: char type[256];
233: /* need to save FILE_CLASS_ID for adapt class */
234: PetscStrncpy(type,((PetscObject)adapt)->type_name,256);
235: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
236: } else if (adapt->ops->view) {
237: (*adapt->ops->view)(adapt,viewer);
238: }
239: return(0);
240: }
244: PetscErrorCode TSAdaptDestroy(TSAdapt *adapt)
245: {
249: if (!*adapt) return(0);
251: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return(0);}
252: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
253: PetscViewerDestroy(&(*adapt)->monitor);
254: PetscHeaderDestroy(adapt);
255: return(0);
256: }
260: /*@
261: TSAdaptSetMonitor - Monitor the choices made by the adaptive controller
263: Collective on TSAdapt
265: Input Arguments:
266: + adapt - adaptive controller context
267: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
269: Level: intermediate
271: .seealso: TSAdaptChoose()
272: @*/
273: PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
274: {
280: if (flg) {
281: if (!adapt->monitor) {PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);}
282: } else {
283: PetscViewerDestroy(&adapt->monitor);
284: }
285: return(0);
286: }
290: /*@C
291: TSAdaptSetCheckStage - set a callback to check convergence for a stage
293: Logically collective on TSAdapt
295: Input Arguments:
296: + adapt - adaptive controller context
297: - func - stage check function
299: Arguments of func:
300: $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)
302: + adapt - adaptive controller context
303: . ts - time stepping context
304: - accept - pending choice of whether to accept, can be modified by this routine
306: Level: advanced
308: .seealso: TSAdaptChoose()
309: @*/
310: PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscBool*))
311: {
315: adapt->ops->checkstage = func;
316: return(0);
317: }
321: /*@
322: TSAdaptSetStepLimits - Set minimum and maximum step sizes to be considered by the controller
324: Logically Collective
326: Input Arguments:
327: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
328: . hmin - minimum time step
329: - hmax - maximum time step
331: Options Database Keys:
332: + -ts_adapt_dt_min - minimum time step
333: - -ts_adapt_dt_max - maximum time step
335: Level: intermediate
337: .seealso: TSAdapt
338: @*/
339: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
340: {
344: if (hmin != PETSC_DECIDE) adapt->dt_min = hmin;
345: if (hmax != PETSC_DECIDE) adapt->dt_max = hmax;
346: return(0);
347: }
351: /*@
352: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
354: Collective on TSAdapt
356: Input Parameter:
357: . adapt - the TSAdapt context
359: Options Database Keys:
360: . -ts_adapt_type <type> - basic
362: Level: advanced
364: Notes:
365: This function is automatically called by TSSetFromOptions()
367: .keywords: TS, TSGetAdapt(), TSAdaptSetType()
369: .seealso: TSGetType()
370: @*/
371: PetscErrorCode TSAdaptSetFromOptions(TSAdapt adapt)
372: {
374: char type[256] = TSADAPTBASIC;
375: PetscBool set,flg;
379: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
380: * function can only be called from inside TSSetFromOptions_GL() */
381: PetscOptionsHead("TS Adaptivity options");
382: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,
383: ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
384: if (flg || !((PetscObject)adapt)->type_name) {
385: TSAdaptSetType(adapt,type);
386: }
387: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",adapt->dt_min,&adapt->dt_min,NULL);
388: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",adapt->dt_max,&adapt->dt_max,NULL);
389: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);
390: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
391: if (set) {TSAdaptSetMonitor(adapt,flg);}
392: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(adapt);}
393: PetscOptionsTail();
394: return(0);
395: }
399: /*@
400: TSAdaptCandidatesClear - clear any previously set candidate schemes
402: Logically Collective
404: Input Argument:
405: . adapt - adaptive controller
407: Level: developer
409: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
410: @*/
411: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
412: {
417: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
418: return(0);
419: }
423: /*@C
424: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
426: Logically Collective
428: Input Arguments:
429: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
430: . name - name of the candidate scheme to add
431: . order - order of the candidate scheme
432: . stageorder - stage order of the candidate scheme
433: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
434: . cost - relative measure of the amount of work required for the candidate scheme
435: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
437: Note:
438: This routine is not available in Fortran.
440: Level: developer
442: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
443: @*/
444: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
445: {
446: PetscInt c;
450: if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
451: if (inuse) {
452: if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
453: adapt->candidates.inuse_set = PETSC_TRUE;
454: }
455: /* first slot if this is the current scheme, otherwise the next available slot */
456: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
458: adapt->candidates.name[c] = name;
459: adapt->candidates.order[c] = order;
460: adapt->candidates.stageorder[c] = stageorder;
461: adapt->candidates.ccfl[c] = ccfl;
462: adapt->candidates.cost[c] = cost;
463: adapt->candidates.n++;
464: return(0);
465: }
469: /*@C
470: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
472: Not Collective
474: Input Arguments:
475: . adapt - time step adaptivity context
477: Output Arguments:
478: + n - number of candidate schemes, always at least 1
479: . order - the order of each candidate scheme
480: . stageorder - the stage order of each candidate scheme
481: . ccfl - the CFL coefficient of each scheme
482: - cost - the relative cost of each scheme
484: Level: developer
486: Note:
487: The current scheme is always returned in the first slot
489: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
490: @*/
491: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
492: {
495: if (n) *n = adapt->candidates.n;
496: if (order) *order = adapt->candidates.order;
497: if (stageorder) *stageorder = adapt->candidates.stageorder;
498: if (ccfl) *ccfl = adapt->candidates.ccfl;
499: if (cost) *cost = adapt->candidates.cost;
500: return(0);
501: }
505: /*@C
506: TSAdaptChoose - choose which method and step size to use for the next step
508: Logically Collective
510: Input Arguments:
511: + adapt - adaptive contoller
512: - h - current step size
514: Output Arguments:
515: + next_sc - scheme to use for the next step
516: . next_h - step size to use for the next step
517: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
519: Note:
520: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
521: being retried after an initial rejection.
523: Level: developer
525: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
526: @*/
527: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
528: {
530: PetscReal wlte = -1.0;
538: if (adapt->candidates.n < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"%D candidates have been registered",adapt->candidates.n);
539: if (!adapt->candidates.inuse_set) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n);
540: (*adapt->ops->choose)(adapt,ts,h,next_sc,next_h,accept,&wlte);
541: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
542: /* Reduce time step if it overshoots max time */
543: PetscReal max_time = ts->max_time;
544: PetscReal next_dt = 0.0;
545: if (ts->ptime + ts->time_step + *next_h >= max_time) {
546: next_dt = max_time - (ts->ptime + ts->time_step);
547: if (next_dt > PETSC_SMALL) *next_h = next_dt;
548: else ts->reason = TS_CONVERGED_TIME;
549: }
550: }
551: if (*next_sc < 0 || adapt->candidates.n <= *next_sc) SETERRQ2(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Chosen scheme %D not in valid range 0..%D",*next_sc,adapt->candidates.n-1);
552: if (!(*next_h > 0.)) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
554: if (adapt->monitor) {
555: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
556: if (wlte < 0) {
557: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt '%s': step %3D %s t=%-11g+%10.3e family='%s' scheme=%D:'%s' dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,((PetscObject)ts)->type_name,*next_sc,adapt->candidates.name[*next_sc],(double)*next_h);
558: } else {
559: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt '%s': step %3D %s t=%-11g+%10.3e wlte=%5.3g family='%s' scheme=%D:'%s' dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)wlte,((PetscObject)ts)->type_name,*next_sc,adapt->candidates.name[*next_sc],(double)*next_h);
560: }
561: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
562: }
563: return(0);
564: }
568: /*@
569: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails)
571: Collective
573: Input Arguments:
574: + adapt - adaptive controller context
575: - ts - time stepper
577: Output Arguments:
578: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
580: Level: developer
582: .seealso:
583: @*/
584: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscBool *accept)
585: {
586: PetscErrorCode ierr;
587: SNES snes;
588: SNESConvergedReason snesreason;
594: *accept = PETSC_TRUE;
595: TSGetSNES(ts,&snes);
596: SNESGetConvergedReason(snes,&snesreason);
597: if (snesreason < 0) {
598: PetscReal dt,new_dt;
599: *accept = PETSC_FALSE;
600: TSGetTimeStep(ts,&dt);
601: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
602: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
603: PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
604: if (adapt->monitor) {
605: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
606: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt '%s': step %3D stage rejected t=%-11g+%10.3e, nonlinear solve failures %D greater than current TS allowed\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,dt,ts->num_snes_failures);
607: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
608: }
609: } else {
610: new_dt = dt*adapt->scale_solve_failed;
611: TSSetTimeStep(ts,new_dt);
612: if (adapt->monitor) {
613: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
614: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt '%s': step %3D stage rejected t=%-11g+%10.3e retrying with dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,(double)dt,(double)new_dt);
615: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
616: }
617: }
618: }
619: if (adapt->ops->checkstage) {(*adapt->ops->checkstage)(adapt,ts,accept);}
620: return(0);
621: }
627: /*@
628: TSAdaptCreate - create an adaptive controller context for time stepping
630: Collective on MPI_Comm
632: Input Parameter:
633: . comm - The communicator
635: Output Parameter:
636: . adapt - new TSAdapt object
638: Level: developer
640: Notes:
641: TSAdapt creation is handled by TS, so users should not need to call this function.
643: .keywords: TSAdapt, create
644: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
645: @*/
646: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
647: {
649: TSAdapt adapt;
653: *inadapt = NULL;
654: TSAdaptInitializePackage();
656: PetscHeaderCreate(adapt,_p_TSAdapt,struct _TSAdaptOps,TSADAPT_CLASSID,"TSAdapt","General Linear adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
658: adapt->dt_min = 1e-20;
659: adapt->dt_max = 1e50;
660: adapt->scale_solve_failed = 0.25;
662: *inadapt = adapt;
663: return(0);
664: }