Actual source code: ts.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
3: #include <petscdmshell.h>
4: #include <petscdmda.h>
5: #include <petscviewer.h>
6: #include <petscdraw.h>
8: /* Logging support */
9: PetscClassId TS_CLASSID, DMTS_CLASSID;
10: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;
12: const char *const TSExactFinalTimeOptions[] = {"STEPOVER","INTERPOLATE","MATCHSTEP","TSExactFinalTimeOption","TS_EXACTFINALTIME_",0};
16: /*
17: TSSetTypeFromOptions - Sets the type of ts from user options.
19: Collective on TS
21: Input Parameter:
22: . ts - The ts
24: Level: intermediate
26: .keywords: TS, set, options, database, type
27: .seealso: TSSetFromOptions(), TSSetType()
28: */
29: static PetscErrorCode TSSetTypeFromOptions(TS ts)
30: {
31: PetscBool opt;
32: const char *defaultType;
33: char typeName[256];
37: if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
38: else defaultType = TSEULER;
40: if (!TSRegisterAllCalled) {TSRegisterAll();}
41: PetscOptionsList("-ts_type", "TS method"," TSSetType", TSList, defaultType, typeName, 256, &opt);
42: if (opt) {
43: TSSetType(ts, typeName);
44: } else {
45: TSSetType(ts, defaultType);
46: }
47: return(0);
48: }
50: struct _n_TSMonitorDrawCtx {
51: PetscViewer viewer;
52: PetscDrawAxis axis;
53: Vec initialsolution;
54: PetscBool showinitial;
55: PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
56: PetscBool showtimestepandtime;
57: int color;
58: };
62: /*@
63: TSSetFromOptions - Sets various TS parameters from user options.
65: Collective on TS
67: Input Parameter:
68: . ts - the TS context obtained from TSCreate()
70: Options Database Keys:
71: + -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSGL, TSSSP
72: . -ts_max_steps maxsteps - maximum number of time-steps to take
73: . -ts_final_time time - maximum time to compute to
74: . -ts_dt dt - initial time step
75: . -ts_monitor - print information at each timestep
76: . -ts_monitor_lg_timestep - Monitor timestep size graphically
77: . -ts_monitor_lg_solution - Monitor solution graphically
78: . -ts_monitor_lg_error - Monitor error graphically
79: . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
80: . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
81: . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
82: . -ts_monitor_draw_solution - Monitor solution graphically
83: . -ts_monitor_draw_solution_phase - Monitor solution graphically with phase diagram
84: . -ts_monitor_draw_error - Monitor error graphically
85: . -ts_monitor_solution_binary <filename> - Save each solution to a binary file
86: - -ts_monitor_solution_vtk <filename.vts> - Save each time step to a binary file, use filename-%%03D.vts
88: Developer Note: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
90: Level: beginner
92: .keywords: TS, timestep, set, options, database
94: .seealso: TSGetType()
95: @*/
96: PetscErrorCode TSSetFromOptions(TS ts)
97: {
98: PetscBool opt,flg;
99: PetscErrorCode ierr;
100: PetscViewer monviewer;
101: char monfilename[PETSC_MAX_PATH_LEN];
102: SNES snes;
103: TSAdapt adapt;
104: PetscReal time_step;
105: TSExactFinalTimeOption eftopt;
106: char dir[16];
110: PetscObjectOptionsBegin((PetscObject)ts);
111: /* Handle TS type options */
112: TSSetTypeFromOptions(ts);
114: /* Handle generic TS options */
115: PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetDuration",ts->max_steps,&ts->max_steps,NULL);
116: PetscOptionsReal("-ts_final_time","Time to run to","TSSetDuration",ts->max_time,&ts->max_time,NULL);
117: PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,NULL);
118: PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&time_step,&flg);
119: if (flg) {
120: TSSetTimeStep(ts,time_step);
121: }
122: PetscOptionsEnum("-ts_exact_final_time","Option for handling of final time step","TSSetExactFinalTime",TSExactFinalTimeOptions,(PetscEnum)ts->exact_final_time,(PetscEnum*)&eftopt,&flg);
123: if (flg) {TSSetExactFinalTime(ts,eftopt);}
124: PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,NULL);
125: PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,NULL);
126: PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",ts->errorifstepfailed,&ts->errorifstepfailed,NULL);
127: PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,NULL);
128: PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,NULL);
130: /* Monitor options */
131: PetscOptionsString("-ts_monitor","Monitor timestep size","TSMonitorDefault","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
132: if (flg) {
133: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),monfilename,&monviewer);
134: TSMonitorSet(ts,TSMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
135: }
136: PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
137: if (flg) {PetscPythonMonitorSet((PetscObject)ts,monfilename);}
139: PetscOptionsName("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",&opt);
140: if (opt) {
141: TSMonitorLGCtx ctx;
142: PetscInt howoften = 1;
144: PetscOptionsInt("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);
145: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
146: TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
147: }
148: PetscOptionsName("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",&opt);
149: if (opt) {
150: TSMonitorLGCtx ctx;
151: PetscInt howoften = 1;
153: PetscOptionsInt("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",howoften,&howoften,NULL);
154: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);
155: TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
156: }
157: PetscOptionsName("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",&opt);
158: if (opt) {
159: TSMonitorLGCtx ctx;
160: PetscInt howoften = 1;
162: PetscOptionsInt("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",howoften,&howoften,NULL);
163: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);
164: TSMonitorSet(ts,TSMonitorLGError,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
165: }
166: PetscOptionsName("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",&opt);
167: if (opt) {
168: TSMonitorLGCtx ctx;
169: PetscInt howoften = 1;
171: PetscOptionsInt("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",howoften,&howoften,NULL);
172: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
173: TSMonitorSet(ts,TSMonitorLGSNESIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
174: }
175: PetscOptionsName("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",&opt);
176: if (opt) {
177: TSMonitorLGCtx ctx;
178: PetscInt howoften = 1;
180: PetscOptionsInt("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",howoften,&howoften,NULL);
181: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
182: TSMonitorSet(ts,TSMonitorLGKSPIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
183: }
184: PetscOptionsName("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",&opt);
185: if (opt) {
186: TSMonitorSPEigCtx ctx;
187: PetscInt howoften = 1;
189: PetscOptionsInt("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",howoften,&howoften,NULL);
190: TSMonitorSPEigCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);
191: TSMonitorSet(ts,TSMonitorSPEig,ctx,(PetscErrorCode (*)(void**))TSMonitorSPEigCtxDestroy);
192: }
193: opt = PETSC_FALSE;
194: PetscOptionsName("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",&opt);
195: if (opt) {
196: TSMonitorDrawCtx ctx;
197: PetscInt howoften = 1;
199: PetscOptionsInt("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",howoften,&howoften,NULL);
200: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);
201: TSMonitorSet(ts,TSMonitorDrawSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
202: }
203: opt = PETSC_FALSE;
204: PetscOptionsName("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",&opt);
205: if (opt) {
206: TSMonitorDrawCtx ctx;
207: PetscReal bounds[4];
208: PetscInt n = 4;
209: PetscDraw draw;
211: PetscOptionsRealArray("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",bounds,&n,NULL);
212: if (n != 4) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Must provide bounding box of phase field");
213: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,1,&ctx);
214: PetscViewerDrawGetDraw(ctx->viewer,0,&draw);
215: PetscDrawClear(draw);
216: PetscDrawAxisCreate(draw,&ctx->axis);
217: PetscDrawAxisSetLimits(ctx->axis,bounds[0],bounds[2],bounds[1],bounds[3]);
218: PetscDrawAxisSetLabels(ctx->axis,"Phase Diagram","Variable 1","Variable 2");
219: PetscDrawAxisDraw(ctx->axis);
220: /* PetscDrawSetCoordinates(draw,bounds[0],bounds[1],bounds[2],bounds[3]); */
221: TSMonitorSet(ts,TSMonitorDrawSolutionPhase,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
222: }
223: opt = PETSC_FALSE;
224: PetscOptionsName("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",&opt);
225: if (opt) {
226: TSMonitorDrawCtx ctx;
227: PetscInt howoften = 1;
229: PetscOptionsInt("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",howoften,&howoften,NULL);
230: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);
231: TSMonitorSet(ts,TSMonitorDrawError,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
232: }
233: opt = PETSC_FALSE;
234: PetscOptionsString("-ts_monitor_solution_binary","Save each solution to a binary file","TSMonitorSolutionBinary",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
235: if (flg) {
236: PetscViewer ctx;
237: if (monfilename[0]) {
238: PetscViewerBinaryOpen(PetscObjectComm((PetscObject)ts),monfilename,FILE_MODE_WRITE,&ctx);
239: TSMonitorSet(ts,TSMonitorSolutionBinary,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);
240: } else {
241: ctx = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)ts));
242: TSMonitorSet(ts,TSMonitorSolutionBinary,ctx,(PetscErrorCode (*)(void**))NULL);
243: }
244: }
245: opt = PETSC_FALSE;
246: PetscOptionsString("-ts_monitor_solution_vtk","Save each time step to a binary file, use filename-%%03D.vts","TSMonitorSolutionVTK",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
247: if (flg) {
248: const char *ptr,*ptr2;
249: char *filetemplate;
250: if (!monfilename[0]) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
251: /* Do some cursory validation of the input. */
252: PetscStrstr(monfilename,"%",(char**)&ptr);
253: if (!ptr) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
254: for (ptr++; ptr && *ptr; ptr++) {
255: PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
256: if (!ptr2 && (*ptr < '0' || '9' < *ptr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03D.vts");
257: if (ptr2) break;
258: }
259: PetscStrallocpy(monfilename,&filetemplate);
260: TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);
261: }
263: PetscOptionsString("-ts_monitor_dmda_ray","Display a ray of the solution","None","y=0",dir,16,&flg);
264: if (flg) {
265: TSMonitorDMDARayCtx *rayctx;
266: int ray = 0;
267: DMDADirection ddir;
268: DM da;
269: PetscMPIInt rank;
271: if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
272: if (dir[0] == 'x') ddir = DMDA_X;
273: else if (dir[0] == 'y') ddir = DMDA_Y;
274: else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
275: sscanf(dir+2,"%d",&ray);
277: PetscInfo2(((PetscObject)ts),"Displaying DMDA ray %c = %D\n",dir[0],ray);
278: PetscNew(TSMonitorDMDARayCtx,&rayctx);
279: TSGetDM(ts,&da);
280: DMDAGetRay(da,ddir,ray,&rayctx->ray,&rayctx->scatter);
281: MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);
282: if (!rank) {
283: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,600,300,&rayctx->viewer);
284: }
285: TSMonitorSet(ts,TSMonitorDMDARay,rayctx,TSMonitorDMDARayDestroy);
286: }
288: TSGetAdapt(ts,&adapt);
289: TSAdaptSetFromOptions(adapt);
291: TSGetSNES(ts,&snes);
292: if (ts->problem_type == TS_LINEAR) {SNESSetType(snes,SNESKSPONLY);}
294: /* Handle specific TS options */
295: if (ts->ops->setfromoptions) {
296: (*ts->ops->setfromoptions)(ts);
297: }
299: /* process any options handlers added with PetscObjectAddOptionsHandler() */
300: PetscObjectProcessOptionsHandlers((PetscObject)ts);
301: PetscOptionsEnd();
302: return(0);
303: }
308: /*@
309: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
310: set with TSSetRHSJacobian().
312: Collective on TS and Vec
314: Input Parameters:
315: + ts - the TS context
316: . t - current timestep
317: - U - input vector
319: Output Parameters:
320: + A - Jacobian matrix
321: . B - optional preconditioning matrix
322: - flag - flag indicating matrix structure
324: Notes:
325: Most users should not need to explicitly call this routine, as it
326: is used internally within the nonlinear solvers.
328: See KSPSetOperators() for important information about setting the
329: flag parameter.
331: Level: developer
333: .keywords: SNES, compute, Jacobian, matrix
335: .seealso: TSSetRHSJacobian(), KSPSetOperators()
336: @*/
337: PetscErrorCode TSComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat *A,Mat *B,MatStructure *flg)
338: {
340: PetscInt Ustate;
341: DM dm;
342: DMTS tsdm;
343: TSRHSJacobian rhsjacobianfunc;
344: void *ctx;
345: TSIJacobian ijacobianfunc;
351: TSGetDM(ts,&dm);
352: DMGetDMTS(dm,&tsdm);
353: DMTSGetRHSJacobian(dm,&rhsjacobianfunc,&ctx);
354: DMTSGetIJacobian(dm,&ijacobianfunc,NULL);
355: PetscObjectStateQuery((PetscObject)U,&Ustate);
356: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.X == U && ts->rhsjacobian.Xstate == Ustate))) {
357: *flg = ts->rhsjacobian.mstructure;
358: return(0);
359: }
361: if (!rhsjacobianfunc && !ijacobianfunc) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
363: if (ts->rhsjacobian.reuse) {
364: MatShift(*A,-ts->rhsjacobian.shift);
365: MatScale(*A,1./ts->rhsjacobian.scale);
366: if (*A != *B) {
367: MatShift(*B,-ts->rhsjacobian.shift);
368: MatScale(*B,1./ts->rhsjacobian.scale);
369: }
370: ts->rhsjacobian.shift = 0;
371: ts->rhsjacobian.scale = 1.;
372: }
374: if (rhsjacobianfunc) {
375: PetscLogEventBegin(TS_JacobianEval,ts,U,*A,*B);
376: *flg = DIFFERENT_NONZERO_PATTERN;
377: PetscStackPush("TS user Jacobian function");
378: (*rhsjacobianfunc)(ts,t,U,A,B,flg,ctx);
379: PetscStackPop;
380: PetscLogEventEnd(TS_JacobianEval,ts,U,*A,*B);
381: /* make sure user returned a correct Jacobian and preconditioner */
384: } else {
385: MatZeroEntries(*A);
386: if (*A != *B) {MatZeroEntries(*B);}
387: *flg = SAME_NONZERO_PATTERN;
388: }
389: ts->rhsjacobian.time = t;
390: ts->rhsjacobian.X = U;
391: PetscObjectStateQuery((PetscObject)U,&ts->rhsjacobian.Xstate);
392: ts->rhsjacobian.mstructure = *flg;
393: return(0);
394: }
398: /*@
399: TSComputeRHSFunction - Evaluates the right-hand-side function.
401: Collective on TS and Vec
403: Input Parameters:
404: + ts - the TS context
405: . t - current time
406: - U - state vector
408: Output Parameter:
409: . y - right hand side
411: Note:
412: Most users should not need to explicitly call this routine, as it
413: is used internally within the nonlinear solvers.
415: Level: developer
417: .keywords: TS, compute
419: .seealso: TSSetRHSFunction(), TSComputeIFunction()
420: @*/
421: PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y)
422: {
424: TSRHSFunction rhsfunction;
425: TSIFunction ifunction;
426: void *ctx;
427: DM dm;
433: TSGetDM(ts,&dm);
434: DMTSGetRHSFunction(dm,&rhsfunction,&ctx);
435: DMTSGetIFunction(dm,&ifunction,NULL);
437: if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
439: PetscLogEventBegin(TS_FunctionEval,ts,U,y,0);
440: if (rhsfunction) {
441: PetscStackPush("TS user right-hand-side function");
442: (*rhsfunction)(ts,t,U,y,ctx);
443: PetscStackPop;
444: } else {
445: VecZeroEntries(y);
446: }
448: PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);
449: return(0);
450: }
454: /*@
455: TSComputeSolutionFunction - Evaluates the solution function.
457: Collective on TS and Vec
459: Input Parameters:
460: + ts - the TS context
461: - t - current time
463: Output Parameter:
464: . U - the solution
466: Note:
467: Most users should not need to explicitly call this routine, as it
468: is used internally within the nonlinear solvers.
470: Level: developer
472: .keywords: TS, compute
474: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
475: @*/
476: PetscErrorCode TSComputeSolutionFunction(TS ts,PetscReal t,Vec U)
477: {
478: PetscErrorCode ierr;
479: TSSolutionFunction solutionfunction;
480: void *ctx;
481: DM dm;
486: TSGetDM(ts,&dm);
487: DMTSGetSolutionFunction(dm,&solutionfunction,&ctx);
489: if (solutionfunction) {
490: PetscStackPush("TS user solution function");
491: (*solutionfunction)(ts,t,U,ctx);
492: PetscStackPop;
493: }
494: return(0);
495: }
498: /*@
499: TSComputeForcingFunction - Evaluates the forcing function.
501: Collective on TS and Vec
503: Input Parameters:
504: + ts - the TS context
505: - t - current time
507: Output Parameter:
508: . U - the function value
510: Note:
511: Most users should not need to explicitly call this routine, as it
512: is used internally within the nonlinear solvers.
514: Level: developer
516: .keywords: TS, compute
518: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
519: @*/
520: PetscErrorCode TSComputeForcingFunction(TS ts,PetscReal t,Vec U)
521: {
522: PetscErrorCode ierr, (*forcing)(TS,PetscReal,Vec,void*);
523: void *ctx;
524: DM dm;
529: TSGetDM(ts,&dm);
530: DMTSGetForcingFunction(dm,&forcing,&ctx);
532: if (forcing) {
533: PetscStackPush("TS user forcing function");
534: (*forcing)(ts,t,U,ctx);
535: PetscStackPop;
536: }
537: return(0);
538: }
542: static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs)
543: {
544: Vec F;
548: *Frhs = NULL;
549: TSGetIFunction(ts,&F,NULL,NULL);
550: if (!ts->Frhs) {
551: VecDuplicate(F,&ts->Frhs);
552: }
553: *Frhs = ts->Frhs;
554: return(0);
555: }
559: static PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs)
560: {
561: Mat A,B;
565: TSGetIJacobian(ts,&A,&B,NULL,NULL);
566: if (Arhs) {
567: if (!ts->Arhs) {
568: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);
569: }
570: *Arhs = ts->Arhs;
571: }
572: if (Brhs) {
573: if (!ts->Brhs) {
574: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);
575: }
576: *Brhs = ts->Brhs;
577: }
578: return(0);
579: }
583: /*@
584: TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0
586: Collective on TS and Vec
588: Input Parameters:
589: + ts - the TS context
590: . t - current time
591: . U - state vector
592: . Udot - time derivative of state vector
593: - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate
595: Output Parameter:
596: . Y - right hand side
598: Note:
599: Most users should not need to explicitly call this routine, as it
600: is used internally within the nonlinear solvers.
602: If the user did did not write their equations in implicit form, this
603: function recasts them in implicit form.
605: Level: developer
607: .keywords: TS, compute
609: .seealso: TSSetIFunction(), TSComputeRHSFunction()
610: @*/
611: PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec Y,PetscBool imex)
612: {
614: TSIFunction ifunction;
615: TSRHSFunction rhsfunction;
616: void *ctx;
617: DM dm;
625: TSGetDM(ts,&dm);
626: DMTSGetIFunction(dm,&ifunction,&ctx);
627: DMTSGetRHSFunction(dm,&rhsfunction,NULL);
629: if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
631: PetscLogEventBegin(TS_FunctionEval,ts,U,Udot,Y);
632: if (ifunction) {
633: PetscStackPush("TS user implicit function");
634: (*ifunction)(ts,t,U,Udot,Y,ctx);
635: PetscStackPop;
636: }
637: if (imex) {
638: if (!ifunction) {
639: VecCopy(Udot,Y);
640: }
641: } else if (rhsfunction) {
642: if (ifunction) {
643: Vec Frhs;
644: TSGetRHSVec_Private(ts,&Frhs);
645: TSComputeRHSFunction(ts,t,U,Frhs);
646: VecAXPY(Y,-1,Frhs);
647: } else {
648: TSComputeRHSFunction(ts,t,U,Y);
649: VecAYPX(Y,-1,Udot);
650: }
651: }
652: PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);
653: return(0);
654: }
658: /*@
659: TSComputeIJacobian - Evaluates the Jacobian of the DAE
661: Collective on TS and Vec
663: Input
664: Input Parameters:
665: + ts - the TS context
666: . t - current timestep
667: . U - state vector
668: . Udot - time derivative of state vector
669: . shift - shift to apply, see note below
670: - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate
672: Output Parameters:
673: + A - Jacobian matrix
674: . B - optional preconditioning matrix
675: - flag - flag indicating matrix structure
677: Notes:
678: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
680: dF/dU + shift*dF/dUdot
682: Most users should not need to explicitly call this routine, as it
683: is used internally within the nonlinear solvers.
685: Level: developer
687: .keywords: TS, compute, Jacobian, matrix
689: .seealso: TSSetIJacobian()
690: @*/
691: PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,PetscBool imex)
692: {
694: TSIJacobian ijacobian;
695: TSRHSJacobian rhsjacobian;
696: DM dm;
697: void *ctx;
709: TSGetDM(ts,&dm);
710: DMTSGetIJacobian(dm,&ijacobian,&ctx);
711: DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);
713: if (!rhsjacobian && !ijacobian) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
715: *flg = SAME_NONZERO_PATTERN; /* In case we're solving a linear problem in which case it wouldn't get initialized below. */
716: PetscLogEventBegin(TS_JacobianEval,ts,U,*A,*B);
717: if (ijacobian) {
718: *flg = DIFFERENT_NONZERO_PATTERN;
719: PetscStackPush("TS user implicit Jacobian");
720: (*ijacobian)(ts,t,U,Udot,shift,A,B,flg,ctx);
721: PetscStackPop;
722: /* make sure user returned a correct Jacobian and preconditioner */
725: }
726: if (imex) {
727: if (!ijacobian) { /* system was written as Udot = G(t,U) */
728: MatZeroEntries(*A);
729: MatShift(*A,shift);
730: if (*A != *B) {
731: MatZeroEntries(*B);
732: MatShift(*B,shift);
733: }
734: *flg = SAME_PRECONDITIONER;
735: }
736: } else {
737: Mat Arhs = NULL,Brhs = NULL;
738: MatStructure flg2;
739: if (rhsjacobian) {
740: TSGetRHSMats_Private(ts,&Arhs,&Brhs);
741: TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);
742: }
743: if (Arhs == *A) { /* No IJacobian, so we only have the RHS matrix */
744: ts->rhsjacobian.scale = -1;
745: ts->rhsjacobian.shift = shift;
746: MatScale(*A,-1);
747: MatShift(*A,shift);
748: if (*A != *B) {
749: MatScale(*B,-1);
750: MatShift(*B,shift);
751: }
752: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
753: MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
754: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
755: MatZeroEntries(*A);
756: MatShift(*A,shift);
757: if (*A != *B) {
758: MatZeroEntries(*B);
759: MatShift(*B,shift);
760: }
761: }
762: MatAXPY(*A,-1,Arhs,axpy);
763: if (*A != *B) {
764: MatAXPY(*B,-1,Brhs,axpy);
765: }
766: *flg = PetscMin(*flg,flg2);
767: }
768: }
770: PetscLogEventEnd(TS_JacobianEval,ts,U,*A,*B);
771: return(0);
772: }
776: /*@C
777: TSSetRHSFunction - Sets the routine for evaluating the function,
778: where U_t = G(t,u).
780: Logically Collective on TS
782: Input Parameters:
783: + ts - the TS context obtained from TSCreate()
784: . r - vector to put the computed right hand side (or NULL to have it created)
785: . f - routine for evaluating the right-hand-side function
786: - ctx - [optional] user-defined context for private data for the
787: function evaluation routine (may be NULL)
789: Calling sequence of func:
790: $ func (TS ts,PetscReal t,Vec u,Vec F,void *ctx);
792: + t - current timestep
793: . u - input vector
794: . F - function vector
795: - ctx - [optional] user-defined function context
797: Level: beginner
799: .keywords: TS, timestep, set, right-hand-side, function
801: .seealso: TSSetRHSJacobian(), TSSetIJacobian()
802: @*/
803: PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx)
804: {
806: SNES snes;
807: Vec ralloc = NULL;
808: DM dm;
814: TSGetDM(ts,&dm);
815: DMTSSetRHSFunction(dm,f,ctx);
816: TSGetSNES(ts,&snes);
817: if (!r && !ts->dm && ts->vec_sol) {
818: VecDuplicate(ts->vec_sol,&ralloc);
819: r = ralloc;
820: }
821: SNESSetFunction(snes,r,SNESTSFormFunction,ts);
822: VecDestroy(&ralloc);
823: return(0);
824: }
828: /*@C
829: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
831: Logically Collective on TS
833: Input Parameters:
834: + ts - the TS context obtained from TSCreate()
835: . f - routine for evaluating the solution
836: - ctx - [optional] user-defined context for private data for the
837: function evaluation routine (may be NULL)
839: Calling sequence of func:
840: $ func (TS ts,PetscReal t,Vec u,void *ctx);
842: + t - current timestep
843: . u - output vector
844: - ctx - [optional] user-defined function context
846: Notes:
847: This routine is used for testing accuracy of time integration schemes when you already know the solution.
848: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
849: create closed-form solutions with non-physical forcing terms.
851: For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history.
853: Level: beginner
855: .keywords: TS, timestep, set, right-hand-side, function
857: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetForcingFunction()
858: @*/
859: PetscErrorCode TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx)
860: {
862: DM dm;
866: TSGetDM(ts,&dm);
867: DMTSSetSolutionFunction(dm,f,ctx);
868: return(0);
869: }
873: /*@C
874: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
876: Logically Collective on TS
878: Input Parameters:
879: + ts - the TS context obtained from TSCreate()
880: . f - routine for evaluating the forcing function
881: - ctx - [optional] user-defined context for private data for the
882: function evaluation routine (may be NULL)
884: Calling sequence of func:
885: $ func (TS ts,PetscReal t,Vec u,void *ctx);
887: + t - current timestep
888: . u - output vector
889: - ctx - [optional] user-defined function context
891: Notes:
892: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
893: create closed-form solutions with a non-physical forcing term.
895: For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history.
897: Level: beginner
899: .keywords: TS, timestep, set, right-hand-side, function
901: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetSolutionFunction()
902: @*/
903: PetscErrorCode TSSetForcingFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx)
904: {
906: DM dm;
910: TSGetDM(ts,&dm);
911: DMTSSetForcingFunction(dm,f,ctx);
912: return(0);
913: }
917: /*@C
918: TSSetRHSJacobian - Sets the function to compute the Jacobian of F,
919: where U_t = G(U,t), as well as the location to store the matrix.
921: Logically Collective on TS
923: Input Parameters:
924: + ts - the TS context obtained from TSCreate()
925: . Amat - (approximate) Jacobian matrix
926: . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
927: . f - the Jacobian evaluation routine
928: - ctx - [optional] user-defined context for private data for the
929: Jacobian evaluation routine (may be NULL)
931: Calling sequence of func:
932: $ func (TS ts,PetscReal t,Vec u,Mat *A,Mat *B,MatStructure *flag,void *ctx);
934: + t - current timestep
935: . u - input vector
936: . Amat - (approximate) Jacobian matrix
937: . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
938: . flag - flag indicating information about the preconditioner matrix
939: structure (same as flag in KSPSetOperators())
940: - ctx - [optional] user-defined context for matrix evaluation routine
942: Notes:
943: See KSPSetOperators() for important information about setting the flag
944: output parameter in the routine func(). Be sure to read this information!
946: The routine func() takes Mat * as the matrix arguments rather than Mat.
947: This allows the matrix evaluation routine to replace A and/or B with a
948: completely new matrix structure (not just different matrix elements)
949: when appropriate, for instance, if the nonzero structure is changing
950: throughout the global iterations.
952: Level: beginner
954: .keywords: TS, timestep, set, right-hand-side, Jacobian
956: .seealso: SNESComputeJacobianDefaultColor(), TSSetRHSFunction(), TSRHSJacobianSetReuse()
958: @*/
959: PetscErrorCode TSSetRHSJacobian(TS ts,Mat Amat,Mat Pmat,TSRHSJacobian f,void *ctx)
960: {
962: SNES snes;
963: DM dm;
964: TSIJacobian ijacobian;
973: TSGetDM(ts,&dm);
974: DMTSSetRHSJacobian(dm,f,ctx);
975: if (f == TSComputeRHSJacobianConstant) {
976: /* Handle this case automatically for the user; otherwise user should call themselves. */
977: TSRHSJacobianSetReuse(ts,PETSC_TRUE);
978: }
979: DMTSGetIJacobian(dm,&ijacobian,NULL);
980: TSGetSNES(ts,&snes);
981: if (!ijacobian) {
982: SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
983: }
984: if (Amat) {
985: PetscObjectReference((PetscObject)Amat);
986: MatDestroy(&ts->Arhs);
988: ts->Arhs = Amat;
989: }
990: if (Pmat) {
991: PetscObjectReference((PetscObject)Pmat);
992: MatDestroy(&ts->Brhs);
994: ts->Brhs = Pmat;
995: }
996: return(0);
997: }
1002: /*@C
1003: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1005: Logically Collective on TS
1007: Input Parameters:
1008: + ts - the TS context obtained from TSCreate()
1009: . r - vector to hold the residual (or NULL to have it created internally)
1010: . f - the function evaluation routine
1011: - ctx - user-defined context for private data for the function evaluation routine (may be NULL)
1013: Calling sequence of f:
1014: $ f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx);
1016: + t - time at step/stage being solved
1017: . u - state vector
1018: . u_t - time derivative of state vector
1019: . F - function vector
1020: - ctx - [optional] user-defined context for matrix evaluation routine
1022: Important:
1023: The user MUST call either this routine, TSSetRHSFunction(). This routine must be used when not solving an ODE, for example a DAE.
1025: Level: beginner
1027: .keywords: TS, timestep, set, DAE, Jacobian
1029: .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian()
1030: @*/
1031: PetscErrorCode TSSetIFunction(TS ts,Vec res,TSIFunction f,void *ctx)
1032: {
1034: SNES snes;
1035: Vec resalloc = NULL;
1036: DM dm;
1042: TSGetDM(ts,&dm);
1043: DMTSSetIFunction(dm,f,ctx);
1045: TSGetSNES(ts,&snes);
1046: if (!res && !ts->dm && ts->vec_sol) {
1047: VecDuplicate(ts->vec_sol,&resalloc);
1048: res = resalloc;
1049: }
1050: SNESSetFunction(snes,res,SNESTSFormFunction,ts);
1051: VecDestroy(&resalloc);
1052: return(0);
1053: }
1057: /*@C
1058: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it.
1060: Not Collective
1062: Input Parameter:
1063: . ts - the TS context
1065: Output Parameter:
1066: + r - vector to hold residual (or NULL)
1067: . func - the function to compute residual (or NULL)
1068: - ctx - the function context (or NULL)
1070: Level: advanced
1072: .keywords: TS, nonlinear, get, function
1074: .seealso: TSSetIFunction(), SNESGetFunction()
1075: @*/
1076: PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx)
1077: {
1079: SNES snes;
1080: DM dm;
1084: TSGetSNES(ts,&snes);
1085: SNESGetFunction(snes,r,NULL,NULL);
1086: TSGetDM(ts,&dm);
1087: DMTSGetIFunction(dm,func,ctx);
1088: return(0);
1089: }
1093: /*@C
1094: TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it.
1096: Not Collective
1098: Input Parameter:
1099: . ts - the TS context
1101: Output Parameter:
1102: + r - vector to hold computed right hand side (or NULL)
1103: . func - the function to compute right hand side (or NULL)
1104: - ctx - the function context (or NULL)
1106: Level: advanced
1108: .keywords: TS, nonlinear, get, function
1110: .seealso: TSSetRhsfunction(), SNESGetFunction()
1111: @*/
1112: PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx)
1113: {
1115: SNES snes;
1116: DM dm;
1120: TSGetSNES(ts,&snes);
1121: SNESGetFunction(snes,r,NULL,NULL);
1122: TSGetDM(ts,&dm);
1123: DMTSGetRHSFunction(dm,func,ctx);
1124: return(0);
1125: }
1129: /*@C
1130: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1131: you provided with TSSetIFunction().
1133: Logically Collective on TS
1135: Input Parameters:
1136: + ts - the TS context obtained from TSCreate()
1137: . Amat - (approximate) Jacobian matrix
1138: . Pmat - matrix used to compute preconditioner (usually the same as Amat)
1139: . f - the Jacobian evaluation routine
1140: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
1142: Calling sequence of f:
1143: $ f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat *Amat,Mat *Pmat,MatStructure *flag,void *ctx);
1145: + t - time at step/stage being solved
1146: . U - state vector
1147: . U_t - time derivative of state vector
1148: . a - shift
1149: . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1150: . Pmat - matrix used for constructing preconditioner, usually the same as Amat
1151: . flag - flag indicating information about the preconditioner matrix
1152: structure (same as flag in KSPSetOperators())
1153: - ctx - [optional] user-defined context for matrix evaluation routine
1155: Notes:
1156: The matrices Amat and Pmat are exactly the matrices that are used by SNES for the nonlinear solve.
1158: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1159: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1160: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1161: a and vector W depend on the integration method, step size, and past states. For example with
1162: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1163: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1165: Level: beginner
1167: .keywords: TS, timestep, DAE, Jacobian
1169: .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESComputeJacobianDefaultColor(), SNESComputeJacobianDefault()
1171: @*/
1172: PetscErrorCode TSSetIJacobian(TS ts,Mat Amat,Mat Pmat,TSIJacobian f,void *ctx)
1173: {
1175: SNES snes;
1176: DM dm;
1185: TSGetDM(ts,&dm);
1186: DMTSSetIJacobian(dm,f,ctx);
1188: TSGetSNES(ts,&snes);
1189: SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1190: return(0);
1191: }
1195: /*@
1196: TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating. Without this flag, TS will change the sign and
1197: shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute
1198: the entire Jacobian. The reuse flag must be set if the evaluation function will assume that the matrix entries have
1199: not been changed by the TS.
1201: Logically Collective
1203: Input Arguments:
1204: + ts - TS context obtained from TSCreate()
1205: - reuse - PETSC_TRUE if the RHS Jacobian
1207: Level: intermediate
1209: .seealso: TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
1210: @*/
1211: PetscErrorCode TSRHSJacobianSetReuse(TS ts,PetscBool reuse)
1212: {
1214: ts->rhsjacobian.reuse = reuse;
1215: return(0);
1216: }
1220: /*@C
1221: TSLoad - Loads a KSP that has been stored in binary with KSPView().
1223: Collective on PetscViewer
1225: Input Parameters:
1226: + newdm - the newly loaded TS, this needs to have been created with TSCreate() or
1227: some related function before a call to TSLoad().
1228: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1230: Level: intermediate
1232: Notes:
1233: The type is determined by the data in the file, any type set into the TS before this call is ignored.
1235: Notes for advanced users:
1236: Most users should not need to know the details of the binary storage
1237: format, since TSLoad() and TSView() completely hide these details.
1238: But for anyone who's interested, the standard binary matrix storage
1239: format is
1240: .vb
1241: has not yet been determined
1242: .ve
1244: .seealso: PetscViewerBinaryOpen(), TSView(), MatLoad(), VecLoad()
1245: @*/
1246: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1247: {
1249: PetscBool isbinary;
1250: PetscInt classid;
1251: char type[256];
1252: DMTS sdm;
1253: DM dm;
1258: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1259: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1261: PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
1262: if (classid != TS_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Not TS next in file");
1263: PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
1264: TSSetType(ts, type);
1265: if (ts->ops->load) {
1266: (*ts->ops->load)(ts,viewer);
1267: }
1268: DMCreate(PetscObjectComm((PetscObject)ts),&dm);
1269: DMLoad(dm,viewer);
1270: TSSetDM(ts,dm);
1271: DMCreateGlobalVector(ts->dm,&ts->vec_sol);
1272: VecLoad(ts->vec_sol,viewer);
1273: DMGetDMTS(ts->dm,&sdm);
1274: DMTSLoad(sdm,viewer);
1275: return(0);
1276: }
1278: #include <petscdraw.h>
1279: #if defined(PETSC_HAVE_AMS)
1280: #include <petscviewerams.h>
1281: #endif
1284: /*@C
1285: TSView - Prints the TS data structure.
1287: Collective on TS
1289: Input Parameters:
1290: + ts - the TS context obtained from TSCreate()
1291: - viewer - visualization context
1293: Options Database Key:
1294: . -ts_view - calls TSView() at end of TSStep()
1296: Notes:
1297: The available visualization contexts include
1298: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1299: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1300: output where only the first processor opens
1301: the file. All other processors send their
1302: data to the first processor to print.
1304: The user can open an alternative visualization context with
1305: PetscViewerASCIIOpen() - output to a specified file.
1307: Level: beginner
1309: .keywords: TS, timestep, view
1311: .seealso: PetscViewerASCIIOpen()
1312: @*/
1313: PetscErrorCode TSView(TS ts,PetscViewer viewer)
1314: {
1316: TSType type;
1317: PetscBool iascii,isstring,isundials,isbinary,isdraw;
1318: DMTS sdm;
1319: #if defined(PETSC_HAVE_AMS)
1320: PetscBool isams;
1321: #endif
1325: if (!viewer) {
1326: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts),&viewer);
1327: }
1331: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1332: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1333: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1334: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1335: #if defined(PETSC_HAVE_AMS)
1336: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);
1337: #endif
1338: if (iascii) {
1339: PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer,"TS Object");
1340: PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);
1341: PetscViewerASCIIPrintf(viewer," maximum time=%G\n",ts->max_time);
1342: if (ts->problem_type == TS_NONLINEAR) {
1343: PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->snes_its);
1344: PetscViewerASCIIPrintf(viewer," total number of nonlinear solve failures=%D\n",ts->num_snes_failures);
1345: }
1346: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->ksp_its);
1347: PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);
1348: DMGetDMTS(ts->dm,&sdm);
1349: DMTSView(sdm,viewer);
1350: if (ts->ops->view) {
1351: PetscViewerASCIIPushTab(viewer);
1352: (*ts->ops->view)(ts,viewer);
1353: PetscViewerASCIIPopTab(viewer);
1354: }
1355: } else if (isstring) {
1356: TSGetType(ts,&type);
1357: PetscViewerStringSPrintf(viewer," %-7.7s",type);
1358: } else if (isbinary) {
1359: PetscInt classid = TS_FILE_CLASSID;
1360: MPI_Comm comm;
1361: PetscMPIInt rank;
1362: char type[256];
1364: PetscObjectGetComm((PetscObject)ts,&comm);
1365: MPI_Comm_rank(comm,&rank);
1366: if (!rank) {
1367: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1368: PetscStrncpy(type,((PetscObject)ts)->type_name,256);
1369: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1370: }
1371: if (ts->ops->view) {
1372: (*ts->ops->view)(ts,viewer);
1373: }
1374: DMView(ts->dm,viewer);
1375: VecView(ts->vec_sol,viewer);
1376: DMGetDMTS(ts->dm,&sdm);
1377: DMTSView(sdm,viewer);
1378: } else if (isdraw) {
1379: PetscDraw draw;
1380: char str[36];
1381: PetscReal x,y,bottom,h;
1383: PetscViewerDrawGetDraw(viewer,0,&draw);
1384: PetscDrawGetCurrentPoint(draw,&x,&y);
1385: PetscStrcpy(str,"TS: ");
1386: PetscStrcat(str,((PetscObject)ts)->type_name);
1387: PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,NULL,&h);
1388: bottom = y - h;
1389: PetscDrawPushCurrentPoint(draw,x,bottom);
1390: if (ts->ops->view) {
1391: (*ts->ops->view)(ts,viewer);
1392: }
1393: PetscDrawPopCurrentPoint(draw);
1394: #if defined(PETSC_HAVE_AMS)
1395: } else if (isams) {
1396: if (((PetscObject)ts)->amsmem == -1) {
1397: PetscObjectViewAMS((PetscObject)ts,viewer);
1398: PetscStackCallAMS(AMS_Memory_take_access,(((PetscObject)ts)->amsmem));
1399: PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ts)->amsmem,"time step",&ts->steps,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
1400: PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ts)->amsmem,"time",&ts->ptime,1,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
1401: PetscStackCallAMS(AMS_Memory_grant_access,(((PetscObject)ts)->amsmem));
1402: }
1403: if (ts->ops->view) {
1404: (*ts->ops->view)(ts,viewer);
1405: }
1406: #endif
1407: }
1409: PetscViewerASCIIPushTab(viewer);
1410: PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);
1411: PetscViewerASCIIPopTab(viewer);
1412: return(0);
1413: }
1418: /*@
1419: TSSetApplicationContext - Sets an optional user-defined context for
1420: the timesteppers.
1422: Logically Collective on TS
1424: Input Parameters:
1425: + ts - the TS context obtained from TSCreate()
1426: - usrP - optional user context
1428: Level: intermediate
1430: .keywords: TS, timestep, set, application, context
1432: .seealso: TSGetApplicationContext()
1433: @*/
1434: PetscErrorCode TSSetApplicationContext(TS ts,void *usrP)
1435: {
1438: ts->user = usrP;
1439: return(0);
1440: }
1444: /*@
1445: TSGetApplicationContext - Gets the user-defined context for the
1446: timestepper.
1448: Not Collective
1450: Input Parameter:
1451: . ts - the TS context obtained from TSCreate()
1453: Output Parameter:
1454: . usrP - user context
1456: Level: intermediate
1458: .keywords: TS, timestep, get, application, context
1460: .seealso: TSSetApplicationContext()
1461: @*/
1462: PetscErrorCode TSGetApplicationContext(TS ts,void *usrP)
1463: {
1466: *(void**)usrP = ts->user;
1467: return(0);
1468: }
1472: /*@
1473: TSGetTimeStepNumber - Gets the number of time steps completed.
1475: Not Collective
1477: Input Parameter:
1478: . ts - the TS context obtained from TSCreate()
1480: Output Parameter:
1481: . iter - number of steps completed so far
1483: Level: intermediate
1485: .keywords: TS, timestep, get, iteration, number
1486: .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStep()
1487: @*/
1488: PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *iter)
1489: {
1493: *iter = ts->steps;
1494: return(0);
1495: }
1499: /*@
1500: TSSetInitialTimeStep - Sets the initial timestep to be used,
1501: as well as the initial time.
1503: Logically Collective on TS
1505: Input Parameters:
1506: + ts - the TS context obtained from TSCreate()
1507: . initial_time - the initial time
1508: - time_step - the size of the timestep
1510: Level: intermediate
1512: .seealso: TSSetTimeStep(), TSGetTimeStep()
1514: .keywords: TS, set, initial, timestep
1515: @*/
1516: PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step)
1517: {
1522: TSSetTimeStep(ts,time_step);
1523: TSSetTime(ts,initial_time);
1524: return(0);
1525: }
1529: /*@
1530: TSSetTimeStep - Allows one to reset the timestep at any time,
1531: useful for simple pseudo-timestepping codes.
1533: Logically Collective on TS
1535: Input Parameters:
1536: + ts - the TS context obtained from TSCreate()
1537: - time_step - the size of the timestep
1539: Level: intermediate
1541: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()
1543: .keywords: TS, set, timestep
1544: @*/
1545: PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step)
1546: {
1550: ts->time_step = time_step;
1551: ts->time_step_orig = time_step;
1552: return(0);
1553: }
1557: /*@
1558: TSSetExactFinalTime - Determines whether to adapt the final time step to
1559: match the exact final time, interpolate solution to the exact final time,
1560: or just return at the final time TS computed.
1562: Logically Collective on TS
1564: Input Parameter:
1565: + ts - the time-step context
1566: - eftopt - exact final time option
1568: Level: beginner
1570: .seealso: TSExactFinalTimeOption
1571: @*/
1572: PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt)
1573: {
1577: ts->exact_final_time = eftopt;
1578: return(0);
1579: }
1583: /*@
1584: TSGetTimeStep - Gets the current timestep size.
1586: Not Collective
1588: Input Parameter:
1589: . ts - the TS context obtained from TSCreate()
1591: Output Parameter:
1592: . dt - the current timestep size
1594: Level: intermediate
1596: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()
1598: .keywords: TS, get, timestep
1599: @*/
1600: PetscErrorCode TSGetTimeStep(TS ts,PetscReal *dt)
1601: {
1605: *dt = ts->time_step;
1606: return(0);
1607: }
1611: /*@
1612: TSGetSolution - Returns the solution at the present timestep. It
1613: is valid to call this routine inside the function that you are evaluating
1614: in order to move to the new timestep. This vector not changed until
1615: the solution at the next timestep has been calculated.
1617: Not Collective, but Vec returned is parallel if TS is parallel
1619: Input Parameter:
1620: . ts - the TS context obtained from TSCreate()
1622: Output Parameter:
1623: . v - the vector containing the solution
1625: Level: intermediate
1627: .seealso: TSGetTimeStep()
1629: .keywords: TS, timestep, get, solution
1630: @*/
1631: PetscErrorCode TSGetSolution(TS ts,Vec *v)
1632: {
1636: *v = ts->vec_sol;
1637: return(0);
1638: }
1640: /* ----- Routines to initialize and destroy a timestepper ---- */
1643: /*@
1644: TSSetProblemType - Sets the type of problem to be solved.
1646: Not collective
1648: Input Parameters:
1649: + ts - The TS
1650: - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
1651: .vb
1652: U_t - A U = 0 (linear)
1653: U_t - A(t) U = 0 (linear)
1654: F(t,U,U_t) = 0 (nonlinear)
1655: .ve
1657: Level: beginner
1659: .keywords: TS, problem type
1660: .seealso: TSSetUp(), TSProblemType, TS
1661: @*/
1662: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
1663: {
1668: ts->problem_type = type;
1669: if (type == TS_LINEAR) {
1670: SNES snes;
1671: TSGetSNES(ts,&snes);
1672: SNESSetType(snes,SNESKSPONLY);
1673: }
1674: return(0);
1675: }
1679: /*@C
1680: TSGetProblemType - Gets the type of problem to be solved.
1682: Not collective
1684: Input Parameter:
1685: . ts - The TS
1687: Output Parameter:
1688: . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
1689: .vb
1690: M U_t = A U
1691: M(t) U_t = A(t) U
1692: F(t,U,U_t)
1693: .ve
1695: Level: beginner
1697: .keywords: TS, problem type
1698: .seealso: TSSetUp(), TSProblemType, TS
1699: @*/
1700: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
1701: {
1705: *type = ts->problem_type;
1706: return(0);
1707: }
1711: /*@
1712: TSSetUp - Sets up the internal data structures for the later use
1713: of a timestepper.
1715: Collective on TS
1717: Input Parameter:
1718: . ts - the TS context obtained from TSCreate()
1720: Notes:
1721: For basic use of the TS solvers the user need not explicitly call
1722: TSSetUp(), since these actions will automatically occur during
1723: the call to TSStep(). However, if one wishes to control this
1724: phase separately, TSSetUp() should be called after TSCreate()
1725: and optional routines of the form TSSetXXX(), but before TSStep().
1727: Level: advanced
1729: .keywords: TS, timestep, setup
1731: .seealso: TSCreate(), TSStep(), TSDestroy()
1732: @*/
1733: PetscErrorCode TSSetUp(TS ts)
1734: {
1736: DM dm;
1737: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
1738: PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
1739: TSIJacobian ijac;
1740: TSRHSJacobian rhsjac;
1744: if (ts->setupcalled) return(0);
1746: if (!((PetscObject)ts)->type_name) {
1747: TSSetType(ts,TSEULER);
1748: }
1750: if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first");
1752: TSGetAdapt(ts,&ts->adapt);
1754: if (ts->rhsjacobian.reuse) {
1755: Mat Amat,Pmat;
1756: SNES snes;
1757: TSGetSNES(ts,&snes);
1758: SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);
1759: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
1760: * have displaced the RHS matrix */
1761: if (Amat == ts->Arhs) {
1762: MatDuplicate(ts->Arhs,MAT_DO_NOT_COPY_VALUES,&Amat);
1763: SNESSetJacobian(snes,Amat,NULL,NULL,NULL);
1764: MatDestroy(&Amat);
1765: }
1766: if (Pmat == ts->Brhs) {
1767: MatDuplicate(ts->Brhs,MAT_DO_NOT_COPY_VALUES,&Pmat);
1768: SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);
1769: MatDestroy(&Pmat);
1770: }
1771: }
1773: if (ts->ops->setup) {
1774: (*ts->ops->setup)(ts);
1775: }
1777: /* in the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
1778: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
1779: */
1780: TSGetDM(ts,&dm);
1781: DMSNESGetFunction(dm,&func,NULL);
1782: if (!func) {
1783: ierr =DMSNESSetFunction(dm,SNESTSFormFunction,ts);
1784: }
1785: /* if the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
1786: Otherwise, the SNES will use coloring internally to form the Jacobian.
1787: */
1788: DMSNESGetJacobian(dm,&jac,NULL);
1789: DMTSGetIJacobian(dm,&ijac,NULL);
1790: DMTSGetRHSJacobian(dm,&rhsjac,NULL);
1791: if (!jac && (ijac || rhsjac)) {
1792: DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);
1793: }
1794: ts->setupcalled = PETSC_TRUE;
1795: return(0);
1796: }
1800: /*@
1801: TSReset - Resets a TS context and removes any allocated Vecs and Mats.
1803: Collective on TS
1805: Input Parameter:
1806: . ts - the TS context obtained from TSCreate()
1808: Level: beginner
1810: .keywords: TS, timestep, reset
1812: .seealso: TSCreate(), TSSetup(), TSDestroy()
1813: @*/
1814: PetscErrorCode TSReset(TS ts)
1815: {
1820: if (ts->ops->reset) {
1821: (*ts->ops->reset)(ts);
1822: }
1823: if (ts->snes) {SNESReset(ts->snes);}
1825: MatDestroy(&ts->Arhs);
1826: MatDestroy(&ts->Brhs);
1827: VecDestroy(&ts->Frhs);
1828: VecDestroy(&ts->vec_sol);
1829: VecDestroy(&ts->vatol);
1830: VecDestroy(&ts->vrtol);
1831: VecDestroyVecs(ts->nwork,&ts->work);
1833: ts->setupcalled = PETSC_FALSE;
1834: return(0);
1835: }
1839: /*@
1840: TSDestroy - Destroys the timestepper context that was created
1841: with TSCreate().
1843: Collective on TS
1845: Input Parameter:
1846: . ts - the TS context obtained from TSCreate()
1848: Level: beginner
1850: .keywords: TS, timestepper, destroy
1852: .seealso: TSCreate(), TSSetUp(), TSSolve()
1853: @*/
1854: PetscErrorCode TSDestroy(TS *ts)
1855: {
1859: if (!*ts) return(0);
1861: if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; return(0);}
1863: TSReset((*ts));
1865: /* if memory was published with AMS then destroy it */
1866: PetscObjectAMSViewOff((PetscObject)*ts);
1867: if ((*ts)->ops->destroy) {(*(*ts)->ops->destroy)((*ts));}
1869: TSAdaptDestroy(&(*ts)->adapt);
1870: SNESDestroy(&(*ts)->snes);
1871: DMDestroy(&(*ts)->dm);
1872: TSMonitorCancel((*ts));
1874: PetscHeaderDestroy(ts);
1875: return(0);
1876: }
1880: /*@
1881: TSGetSNES - Returns the SNES (nonlinear solver) associated with
1882: a TS (timestepper) context. Valid only for nonlinear problems.
1884: Not Collective, but SNES is parallel if TS is parallel
1886: Input Parameter:
1887: . ts - the TS context obtained from TSCreate()
1889: Output Parameter:
1890: . snes - the nonlinear solver context
1892: Notes:
1893: The user can then directly manipulate the SNES context to set various
1894: options, etc. Likewise, the user can then extract and manipulate the
1895: KSP, KSP, and PC contexts as well.
1897: TSGetSNES() does not work for integrators that do not use SNES; in
1898: this case TSGetSNES() returns NULL in snes.
1900: Level: beginner
1902: .keywords: timestep, get, SNES
1903: @*/
1904: PetscErrorCode TSGetSNES(TS ts,SNES *snes)
1905: {
1911: if (!ts->snes) {
1912: SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);
1913: SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
1914: PetscLogObjectParent(ts,ts->snes);
1915: PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);
1916: if (ts->dm) {SNESSetDM(ts->snes,ts->dm);}
1917: if (ts->problem_type == TS_LINEAR) {
1918: SNESSetType(ts->snes,SNESKSPONLY);
1919: }
1920: }
1921: *snes = ts->snes;
1922: return(0);
1923: }
1927: /*@
1928: TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context
1930: Collective
1932: Input Parameter:
1933: + ts - the TS context obtained from TSCreate()
1934: - snes - the nonlinear solver context
1936: Notes:
1937: Most users should have the TS created by calling TSGetSNES()
1939: Level: developer
1941: .keywords: timestep, set, SNES
1942: @*/
1943: PetscErrorCode TSSetSNES(TS ts,SNES snes)
1944: {
1946: PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
1951: PetscObjectReference((PetscObject)snes);
1952: SNESDestroy(&ts->snes);
1954: ts->snes = snes;
1956: SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
1957: SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);
1958: if (func == SNESTSFormJacobian) {
1959: SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);
1960: }
1961: return(0);
1962: }
1966: /*@
1967: TSGetKSP - Returns the KSP (linear solver) associated with
1968: a TS (timestepper) context.
1970: Not Collective, but KSP is parallel if TS is parallel
1972: Input Parameter:
1973: . ts - the TS context obtained from TSCreate()
1975: Output Parameter:
1976: . ksp - the nonlinear solver context
1978: Notes:
1979: The user can then directly manipulate the KSP context to set various
1980: options, etc. Likewise, the user can then extract and manipulate the
1981: KSP and PC contexts as well.
1983: TSGetKSP() does not work for integrators that do not use KSP;
1984: in this case TSGetKSP() returns NULL in ksp.
1986: Level: beginner
1988: .keywords: timestep, get, KSP
1989: @*/
1990: PetscErrorCode TSGetKSP(TS ts,KSP *ksp)
1991: {
1993: SNES snes;
1998: if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first");
1999: if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()");
2000: TSGetSNES(ts,&snes);
2001: SNESGetKSP(snes,ksp);
2002: return(0);
2003: }
2005: /* ----------- Routines to set solver parameters ---------- */
2009: /*@
2010: TSGetDuration - Gets the maximum number of timesteps to use and
2011: maximum time for iteration.
2013: Not Collective
2015: Input Parameters:
2016: + ts - the TS context obtained from TSCreate()
2017: . maxsteps - maximum number of iterations to use, or NULL
2018: - maxtime - final time to iterate to, or NULL
2020: Level: intermediate
2022: .keywords: TS, timestep, get, maximum, iterations, time
2023: @*/
2024: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2025: {
2028: if (maxsteps) {
2030: *maxsteps = ts->max_steps;
2031: }
2032: if (maxtime) {
2034: *maxtime = ts->max_time;
2035: }
2036: return(0);
2037: }
2041: /*@
2042: TSSetDuration - Sets the maximum number of timesteps to use and
2043: maximum time for iteration.
2045: Logically Collective on TS
2047: Input Parameters:
2048: + ts - the TS context obtained from TSCreate()
2049: . maxsteps - maximum number of iterations to use
2050: - maxtime - final time to iterate to
2052: Options Database Keys:
2053: . -ts_max_steps <maxsteps> - Sets maxsteps
2054: . -ts_final_time <maxtime> - Sets maxtime
2056: Notes:
2057: The default maximum number of iterations is 5000. Default time is 5.0
2059: Level: intermediate
2061: .keywords: TS, timestep, set, maximum, iterations
2063: .seealso: TSSetExactFinalTime()
2064: @*/
2065: PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime)
2066: {
2071: if (maxsteps >= 0) ts->max_steps = maxsteps;
2072: if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
2073: return(0);
2074: }
2078: /*@
2079: TSSetSolution - Sets the initial solution vector
2080: for use by the TS routines.
2082: Logically Collective on TS and Vec
2084: Input Parameters:
2085: + ts - the TS context obtained from TSCreate()
2086: - u - the solution vector
2088: Level: beginner
2090: .keywords: TS, timestep, set, solution, initial conditions
2091: @*/
2092: PetscErrorCode TSSetSolution(TS ts,Vec u)
2093: {
2095: DM dm;
2100: PetscObjectReference((PetscObject)u);
2101: VecDestroy(&ts->vec_sol);
2103: ts->vec_sol = u;
2105: TSGetDM(ts,&dm);
2106: DMShellSetGlobalVector(dm,u);
2107: return(0);
2108: }
2112: /*@C
2113: TSSetPreStep - Sets the general-purpose function
2114: called once at the beginning of each time step.
2116: Logically Collective on TS
2118: Input Parameters:
2119: + ts - The TS context obtained from TSCreate()
2120: - func - The function
2122: Calling sequence of func:
2123: . func (TS ts);
2125: Level: intermediate
2127: Note:
2128: If a step is rejected, TSStep() will call this routine again before each attempt.
2129: The last completed time step number can be queried using TSGetTimeStepNumber(), the
2130: size of the step being attempted can be obtained using TSGetTimeStep().
2132: .keywords: TS, timestep
2133: .seealso: TSSetPreStage(), TSSetPostStep(), TSStep()
2134: @*/
2135: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
2136: {
2139: ts->prestep = func;
2140: return(0);
2141: }
2145: /*@
2146: TSPreStep - Runs the user-defined pre-step function.
2148: Collective on TS
2150: Input Parameters:
2151: . ts - The TS context obtained from TSCreate()
2153: Notes:
2154: TSPreStep() is typically used within time stepping implementations,
2155: so most users would not generally call this routine themselves.
2157: Level: developer
2159: .keywords: TS, timestep
2160: .seealso: TSSetPreStep(), TSPreStage(), TSPostStep()
2161: @*/
2162: PetscErrorCode TSPreStep(TS ts)
2163: {
2168: if (ts->prestep) {
2169: PetscStackCallStandard((*ts->prestep),(ts));
2170: }
2171: return(0);
2172: }
2176: /*@C
2177: TSSetPreStage - Sets the general-purpose function
2178: called once at the beginning of each stage.
2180: Logically Collective on TS
2182: Input Parameters:
2183: + ts - The TS context obtained from TSCreate()
2184: - func - The function
2186: Calling sequence of func:
2187: . PetscErrorCode func(TS ts, PetscReal stagetime);
2189: Level: intermediate
2191: Note:
2192: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
2193: The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being
2194: attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime().
2196: .keywords: TS, timestep
2197: .seealso: TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
2198: @*/
2199: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal))
2200: {
2203: ts->prestage = func;
2204: return(0);
2205: }
2209: /*@
2210: TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage()
2212: Collective on TS
2214: Input Parameters:
2215: . ts - The TS context obtained from TSCreate()
2217: Notes:
2218: TSPreStage() is typically used within time stepping implementations,
2219: most users would not generally call this routine themselves.
2221: Level: developer
2223: .keywords: TS, timestep
2224: .seealso: TSSetPreStep(), TSPreStep(), TSPostStep()
2225: @*/
2226: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
2227: {
2232: if (ts->prestage) {
2233: PetscStackCallStandard((*ts->prestage),(ts,stagetime));
2234: }
2235: return(0);
2236: }
2240: /*@C
2241: TSSetPostStep - Sets the general-purpose function
2242: called once at the end of each time step.
2244: Logically Collective on TS
2246: Input Parameters:
2247: + ts - The TS context obtained from TSCreate()
2248: - func - The function
2250: Calling sequence of func:
2251: $ func (TS ts);
2253: Level: intermediate
2255: .keywords: TS, timestep
2256: .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime()
2257: @*/
2258: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
2259: {
2262: ts->poststep = func;
2263: return(0);
2264: }
2268: /*@
2269: TSPostStep - Runs the user-defined post-step function.
2271: Collective on TS
2273: Input Parameters:
2274: . ts - The TS context obtained from TSCreate()
2276: Notes:
2277: TSPostStep() is typically used within time stepping implementations,
2278: so most users would not generally call this routine themselves.
2280: Level: developer
2282: .keywords: TS, timestep
2283: @*/
2284: PetscErrorCode TSPostStep(TS ts)
2285: {
2290: if (ts->poststep) {
2291: PetscStackCallStandard((*ts->poststep),(ts));
2292: }
2293: return(0);
2294: }
2296: /* ------------ Routines to set performance monitoring options ----------- */
2300: /*@C
2301: TSMonitorSet - Sets an ADDITIONAL function that is to be used at every
2302: timestep to display the iteration's progress.
2304: Logically Collective on TS
2306: Input Parameters:
2307: + ts - the TS context obtained from TSCreate()
2308: . monitor - monitoring routine
2309: . mctx - [optional] user-defined context for private data for the
2310: monitor routine (use NULL if no context is desired)
2311: - monitordestroy - [optional] routine that frees monitor context
2312: (may be NULL)
2314: Calling sequence of monitor:
2315: $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
2317: + ts - the TS context
2318: . steps - iteration number (after the final time step the monitor routine is called with a step of -1, this is at the final time which may have
2319: been interpolated to)
2320: . time - current time
2321: . u - current iterate
2322: - mctx - [optional] monitoring context
2324: Notes:
2325: This routine adds an additional monitor to the list of monitors that
2326: already has been loaded.
2328: Fortran notes: Only a single monitor function can be set for each TS object
2330: Level: intermediate
2332: .keywords: TS, timestep, set, monitor
2334: .seealso: TSMonitorDefault(), TSMonitorCancel()
2335: @*/
2336: PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**))
2337: {
2340: if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
2341: ts->monitor[ts->numbermonitors] = monitor;
2342: ts->monitordestroy[ts->numbermonitors] = mdestroy;
2343: ts->monitorcontext[ts->numbermonitors++] = (void*)mctx;
2344: return(0);
2345: }
2349: /*@C
2350: TSMonitorCancel - Clears all the monitors that have been set on a time-step object.
2352: Logically Collective on TS
2354: Input Parameters:
2355: . ts - the TS context obtained from TSCreate()
2357: Notes:
2358: There is no way to remove a single, specific monitor.
2360: Level: intermediate
2362: .keywords: TS, timestep, set, monitor
2364: .seealso: TSMonitorDefault(), TSMonitorSet()
2365: @*/
2366: PetscErrorCode TSMonitorCancel(TS ts)
2367: {
2369: PetscInt i;
2373: for (i=0; i<ts->numbermonitors; i++) {
2374: if (ts->monitordestroy[i]) {
2375: (*ts->monitordestroy[i])(&ts->monitorcontext[i]);
2376: }
2377: }
2378: ts->numbermonitors = 0;
2379: return(0);
2380: }
2384: /*@
2385: TSMonitorDefault - Sets the Default monitor
2387: Level: intermediate
2389: .keywords: TS, set, monitor
2391: .seealso: TSMonitorDefault(), TSMonitorSet()
2392: @*/
2393: PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
2394: {
2396: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts));
2399: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
2400: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);
2401: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
2402: return(0);
2403: }
2407: /*@
2408: TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available.
2410: Logically Collective on TS
2412: Input Argument:
2413: . ts - time stepping context
2415: Output Argument:
2416: . flg - PETSC_TRUE or PETSC_FALSE
2418: Level: intermediate
2420: .keywords: TS, set
2422: .seealso: TSInterpolate(), TSSetPostStep()
2423: @*/
2424: PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg)
2425: {
2428: ts->retain_stages = flg;
2429: return(0);
2430: }
2434: /*@
2435: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
2437: Collective on TS
2439: Input Argument:
2440: + ts - time stepping context
2441: - t - time to interpolate to
2443: Output Argument:
2444: . U - state at given time
2446: Notes:
2447: The user should call TSSetRetainStages() before taking a step in which interpolation will be requested.
2449: Level: intermediate
2451: Developer Notes:
2452: TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
2454: .keywords: TS, set
2456: .seealso: TSSetRetainStages(), TSSetPostStep()
2457: @*/
2458: PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U)
2459: {
2465: if (t < ts->ptime - ts->time_step_prev || t > ts->ptime) SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Requested time %G not in last time steps [%G,%G]",t,ts->ptime-ts->time_step_prev,ts->ptime);
2466: if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name);
2467: (*ts->ops->interpolate)(ts,t,U);
2468: return(0);
2469: }
2473: /*@
2474: TSStep - Steps one time step
2476: Collective on TS
2478: Input Parameter:
2479: . ts - the TS context obtained from TSCreate()
2481: Level: intermediate
2483: Notes:
2484: The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may
2485: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
2487: This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the
2488: time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep.
2490: .keywords: TS, timestep, solve
2492: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate()
2493: @*/
2494: PetscErrorCode TSStep(TS ts)
2495: {
2496: PetscReal ptime_prev;
2501: TSSetUp(ts);
2503: ts->reason = TS_CONVERGED_ITERATING;
2504: ptime_prev = ts->ptime;
2506: PetscLogEventBegin(TS_Step,ts,0,0,0);
2507: (*ts->ops->step)(ts);
2508: PetscLogEventEnd(TS_Step,ts,0,0,0);
2510: ts->time_step_prev = ts->ptime - ptime_prev;
2512: if (ts->reason < 0) {
2513: if (ts->errorifstepfailed) {
2514: if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) {
2515: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s, increase -ts_max_snes_failures or make negative to attempt recovery",TSConvergedReasons[ts->reason]);
2516: } else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
2517: }
2518: } else if (!ts->reason) {
2519: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
2520: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
2521: }
2522: return(0);
2523: }
2527: /*@
2528: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
2530: Collective on TS
2532: Input Arguments:
2533: + ts - time stepping context
2534: . order - desired order of accuracy
2535: - done - whether the step was evaluated at this order (pass NULL to generate an error if not available)
2537: Output Arguments:
2538: . U - state at the end of the current step
2540: Level: advanced
2542: Notes:
2543: This function cannot be called until all stages have been evaluated.
2544: It is normally called by adaptive controllers before a step has been accepted and may also be called by the user after TSStep() has returned.
2546: .seealso: TSStep(), TSAdapt
2547: @*/
2548: PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done)
2549: {
2556: if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name);
2557: (*ts->ops->evaluatestep)(ts,order,U,done);
2558: return(0);
2559: }
2563: /*@
2564: TSSolve - Steps the requested number of timesteps.
2566: Collective on TS
2568: Input Parameter:
2569: + ts - the TS context obtained from TSCreate()
2570: - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions)
2572: Level: beginner
2574: Notes:
2575: The final time returned by this function may be different from the time of the internally
2576: held state accessible by TSGetSolution() and TSGetTime() because the method may have
2577: stepped over the final time.
2579: .keywords: TS, timestep, solve
2581: .seealso: TSCreate(), TSSetSolution(), TSStep()
2582: @*/
2583: PetscErrorCode TSSolve(TS ts,Vec u)
2584: {
2585: PetscBool flg;
2586: PetscViewer viewer;
2587: Vec solution;
2588: PetscErrorCode ierr;
2589: PetscViewerFormat format;
2594: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE) { /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */
2596: if (!ts->vec_sol || u == ts->vec_sol) {
2597: VecDuplicate(u,&solution);
2598: TSSetSolution(ts,solution);
2599: VecDestroy(&solution); /* grant ownership */
2600: }
2601: VecCopy(u,ts->vec_sol);
2602: } else if (u) {
2603: TSSetSolution(ts,u);
2604: }
2605: TSSetUp(ts);
2606: /* reset time step and iteration counters */
2607: ts->steps = 0;
2608: ts->ksp_its = 0;
2609: ts->snes_its = 0;
2610: ts->num_snes_failures = 0;
2611: ts->reject = 0;
2612: ts->reason = TS_CONVERGED_ITERATING;
2614: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,"-ts_view_pre",&viewer,&format,&flg);
2615: if (flg && !PetscPreLoadingOn) {
2616: PetscViewerPushFormat(viewer,format);
2617: TSView(ts,viewer);
2618: PetscViewerPopFormat(viewer);
2619: PetscViewerDestroy(&viewer);
2620: }
2622: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
2623: (*ts->ops->solve)(ts);
2624: VecCopy(ts->vec_sol,u);
2625: ts->solvetime = ts->ptime;
2626: } else {
2627: /* steps the requested number of timesteps. */
2628: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
2629: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
2630: while (!ts->reason) {
2631: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
2632: TSStep(ts);
2633: TSPostStep(ts);
2634: }
2635: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
2636: TSInterpolate(ts,ts->max_time,u);
2637: ts->solvetime = ts->max_time;
2638: solution = u;
2639: } else {
2640: if (u) {VecCopy(ts->vec_sol,u);}
2641: ts->solvetime = ts->ptime;
2642: solution = ts->vec_sol;
2643: }
2644: TSMonitor(ts,ts->steps,ts->solvetime,solution);
2645: }
2646: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,"-ts_view",&viewer,&format,&flg);
2647: if (flg && !PetscPreLoadingOn) {
2648: PetscViewerPushFormat(viewer,format);
2649: TSView(ts,viewer);
2650: PetscViewerPopFormat(viewer);
2651: PetscViewerDestroy(&viewer);
2652: }
2653: return(0);
2654: }
2658: /*@
2659: TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet()
2661: Collective on TS
2663: Input Parameters:
2664: + ts - time stepping context obtained from TSCreate()
2665: . step - step number that has just completed
2666: . ptime - model time of the state
2667: - u - state at the current model time
2669: Notes:
2670: TSMonitor() is typically used within the time stepping implementations.
2671: Users might call this function when using the TSStep() interface instead of TSSolve().
2673: Level: advanced
2675: .keywords: TS, timestep
2676: @*/
2677: PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u)
2678: {
2680: PetscInt i,n = ts->numbermonitors;
2685: for (i=0; i<n; i++) {
2686: (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);
2687: }
2688: return(0);
2689: }
2691: /* ------------------------------------------------------------------------*/
2692: struct _n_TSMonitorLGCtx {
2693: PetscDrawLG lg;
2694: PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
2695: PetscInt ksp_its,snes_its;
2696: };
2701: /*@C
2702: TSMonitorLGCtxCreate - Creates a line graph context for use with
2703: TS to monitor the solution process graphically in various ways
2705: Collective on TS
2707: Input Parameters:
2708: + host - the X display to open, or null for the local machine
2709: . label - the title to put in the title bar
2710: . x, y - the screen coordinates of the upper left coordinate of the window
2711: . m, n - the screen width and height in pixels
2712: - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time
2714: Output Parameter:
2715: . ctx - the context
2717: Options Database Key:
2718: + -ts_monitor_lg_timestep - automatically sets line graph monitor
2719: . -ts_monitor_lg_solution -
2720: . -ts_monitor_lg_error -
2721: . -ts_monitor_lg_ksp_iterations -
2722: . -ts_monitor_lg_snes_iterations -
2723: - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true
2725: Notes:
2726: Use TSMonitorLGCtxDestroy() to destroy.
2728: Level: intermediate
2730: .keywords: TS, monitor, line graph, residual, seealso
2732: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()
2734: @*/
2735: PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx)
2736: {
2737: PetscDraw win;
2739: PetscBool flg = PETSC_TRUE;
2742: PetscNew(struct _n_TSMonitorLGCtx,ctx);
2743: PetscDrawCreate(comm,host,label,x,y,m,n,&win);
2744: PetscDrawSetFromOptions(win);
2745: PetscDrawLGCreate(win,1,&(*ctx)->lg);
2746: PetscOptionsGetBool(NULL,"-lg_indicate_data_points",&flg,NULL);
2747: if (flg) {
2748: PetscDrawLGIndicateDataPoints((*ctx)->lg);
2749: }
2750: PetscLogObjectParent((*ctx)->lg,win);
2751: (*ctx)->howoften = howoften;
2752: return(0);
2753: }
2757: PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
2758: {
2759: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
2760: PetscReal x = ptime,y;
2764: if (!step) {
2765: PetscDrawAxis axis;
2766: PetscDrawLGGetAxis(ctx->lg,&axis);
2767: PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");
2768: PetscDrawLGReset(ctx->lg);
2769: }
2770: TSGetTimeStep(ts,&y);
2771: PetscDrawLGAddPoint(ctx->lg,&x,&y);
2772: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
2773: PetscDrawLGDraw(ctx->lg);
2774: }
2775: return(0);
2776: }
2780: /*@C
2781: TSMonitorLGCtxDestroy - Destroys a line graph context that was created
2782: with TSMonitorLGCtxCreate().
2784: Collective on TSMonitorLGCtx
2786: Input Parameter:
2787: . ctx - the monitor context
2789: Level: intermediate
2791: .keywords: TS, monitor, line graph, destroy
2793: .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep();
2794: @*/
2795: PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx)
2796: {
2797: PetscDraw draw;
2801: PetscDrawLGGetDraw((*ctx)->lg,&draw);
2802: PetscDrawDestroy(&draw);
2803: PetscDrawLGDestroy(&(*ctx)->lg);
2804: PetscFree(*ctx);
2805: return(0);
2806: }
2810: /*@
2811: TSGetTime - Gets the time of the most recently completed step.
2813: Not Collective
2815: Input Parameter:
2816: . ts - the TS context obtained from TSCreate()
2818: Output Parameter:
2819: . t - the current time
2821: Level: beginner
2823: Note:
2824: When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(),
2825: TSSetPreStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated.
2827: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()
2829: .keywords: TS, get, time
2830: @*/
2831: PetscErrorCode TSGetTime(TS ts,PetscReal *t)
2832: {
2836: *t = ts->ptime;
2837: return(0);
2838: }
2842: /*@
2843: TSSetTime - Allows one to reset the time.
2845: Logically Collective on TS
2847: Input Parameters:
2848: + ts - the TS context obtained from TSCreate()
2849: - time - the time
2851: Level: intermediate
2853: .seealso: TSGetTime(), TSSetDuration()
2855: .keywords: TS, set, time
2856: @*/
2857: PetscErrorCode TSSetTime(TS ts, PetscReal t)
2858: {
2862: ts->ptime = t;
2863: return(0);
2864: }
2868: /*@C
2869: TSSetOptionsPrefix - Sets the prefix used for searching for all
2870: TS options in the database.
2872: Logically Collective on TS
2874: Input Parameter:
2875: + ts - The TS context
2876: - prefix - The prefix to prepend to all option names
2878: Notes:
2879: A hyphen (-) must NOT be given at the beginning of the prefix name.
2880: The first character of all runtime options is AUTOMATICALLY the
2881: hyphen.
2883: Level: advanced
2885: .keywords: TS, set, options, prefix, database
2887: .seealso: TSSetFromOptions()
2889: @*/
2890: PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[])
2891: {
2893: SNES snes;
2897: PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);
2898: TSGetSNES(ts,&snes);
2899: SNESSetOptionsPrefix(snes,prefix);
2900: return(0);
2901: }
2906: /*@C
2907: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
2908: TS options in the database.
2910: Logically Collective on TS
2912: Input Parameter:
2913: + ts - The TS context
2914: - prefix - The prefix to prepend to all option names
2916: Notes:
2917: A hyphen (-) must NOT be given at the beginning of the prefix name.
2918: The first character of all runtime options is AUTOMATICALLY the
2919: hyphen.
2921: Level: advanced
2923: .keywords: TS, append, options, prefix, database
2925: .seealso: TSGetOptionsPrefix()
2927: @*/
2928: PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[])
2929: {
2931: SNES snes;
2935: PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);
2936: TSGetSNES(ts,&snes);
2937: SNESAppendOptionsPrefix(snes,prefix);
2938: return(0);
2939: }
2943: /*@C
2944: TSGetOptionsPrefix - Sets the prefix used for searching for all
2945: TS options in the database.
2947: Not Collective
2949: Input Parameter:
2950: . ts - The TS context
2952: Output Parameter:
2953: . prefix - A pointer to the prefix string used
2955: Notes: On the fortran side, the user should pass in a string 'prifix' of
2956: sufficient length to hold the prefix.
2958: Level: intermediate
2960: .keywords: TS, get, options, prefix, database
2962: .seealso: TSAppendOptionsPrefix()
2963: @*/
2964: PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[])
2965: {
2971: PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);
2972: return(0);
2973: }
2977: /*@C
2978: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
2980: Not Collective, but parallel objects are returned if TS is parallel
2982: Input Parameter:
2983: . ts - The TS context obtained from TSCreate()
2985: Output Parameters:
2986: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL)
2987: . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL)
2988: . func - Function to compute the Jacobian of the RHS (or NULL)
2989: - ctx - User-defined context for Jacobian evaluation routine (or NULL)
2991: Notes: You can pass in NULL for any return argument you do not need.
2993: Level: intermediate
2995: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber()
2997: .keywords: TS, timestep, get, matrix, Jacobian
2998: @*/
2999: PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx)
3000: {
3002: SNES snes;
3003: DM dm;
3006: TSGetSNES(ts,&snes);
3007: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
3008: TSGetDM(ts,&dm);
3009: DMTSGetRHSJacobian(dm,func,ctx);
3010: return(0);
3011: }
3015: /*@C
3016: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
3018: Not Collective, but parallel objects are returned if TS is parallel
3020: Input Parameter:
3021: . ts - The TS context obtained from TSCreate()
3023: Output Parameters:
3024: + Amat - The (approximate) Jacobian of F(t,U,U_t)
3025: . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat
3026: . f - The function to compute the matrices
3027: - ctx - User-defined context for Jacobian evaluation routine
3029: Notes: You can pass in NULL for any return argument you do not need.
3031: Level: advanced
3033: .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber()
3035: .keywords: TS, timestep, get, matrix, Jacobian
3036: @*/
3037: PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx)
3038: {
3040: SNES snes;
3041: DM dm;
3044: TSGetSNES(ts,&snes);
3045: SNESSetUpMatrices(snes);
3046: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
3047: TSGetDM(ts,&dm);
3048: DMTSGetIJacobian(dm,f,ctx);
3049: return(0);
3050: }
3055: /*@C
3056: TSMonitorDrawSolution - Monitors progress of the TS solvers by calling
3057: VecView() for the solution at each timestep
3059: Collective on TS
3061: Input Parameters:
3062: + ts - the TS context
3063: . step - current time-step
3064: . ptime - current time
3065: - dummy - either a viewer or NULL
3067: Options Database:
3068: . -ts_monitor_draw_solution_initial - show initial solution as well as current solution
3070: Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial
3071: will look bad
3073: Level: intermediate
3075: .keywords: TS, vector, monitor, view
3077: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
3078: @*/
3079: PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
3080: {
3081: PetscErrorCode ierr;
3082: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
3083: PetscDraw draw;
3086: if (!step && ictx->showinitial) {
3087: if (!ictx->initialsolution) {
3088: VecDuplicate(u,&ictx->initialsolution);
3089: }
3090: VecCopy(u,ictx->initialsolution);
3091: }
3092: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);
3094: if (ictx->showinitial) {
3095: PetscReal pause;
3096: PetscViewerDrawGetPause(ictx->viewer,&pause);
3097: PetscViewerDrawSetPause(ictx->viewer,0.0);
3098: VecView(ictx->initialsolution,ictx->viewer);
3099: PetscViewerDrawSetPause(ictx->viewer,pause);
3100: PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);
3101: }
3102: VecView(u,ictx->viewer);
3103: if (ictx->showtimestepandtime) {
3104: PetscReal xl,yl,xr,yr,tw,w,h;
3105: char time[32];
3106: size_t len;
3108: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
3109: PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);
3110: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
3111: PetscStrlen(time,&len);
3112: PetscDrawStringGetSize(draw,&tw,NULL);
3113: w = xl + .5*(xr - xl) - .5*len*tw;
3114: h = yl + .95*(yr - yl);
3115: PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);
3116: PetscDrawFlush(draw);
3117: }
3119: if (ictx->showinitial) {
3120: PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);
3121: }
3122: return(0);
3123: }
3127: /*@C
3128: TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram
3130: Collective on TS
3132: Input Parameters:
3133: + ts - the TS context
3134: . step - current time-step
3135: . ptime - current time
3136: - dummy - either a viewer or NULL
3138: Level: intermediate
3140: .keywords: TS, vector, monitor, view
3142: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
3143: @*/
3144: PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
3145: {
3146: PetscErrorCode ierr;
3147: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
3148: PetscDraw draw;
3149: MPI_Comm comm;
3150: PetscInt n;
3151: PetscMPIInt size;
3152: PetscReal xl,yl,xr,yr,tw,w,h;
3153: char time[32];
3154: size_t len;
3155: const PetscScalar *U;
3158: PetscObjectGetComm((PetscObject)ts,&comm);
3159: MPI_Comm_size(comm,&size);
3160: if (size != 1) SETERRQ(comm,PETSC_ERR_SUP,"Only allowed for sequential runs");
3161: VecGetSize(u,&n);
3162: if (n != 2) SETERRQ(comm,PETSC_ERR_SUP,"Only for ODEs with two unknowns");
3164: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
3166: VecGetArrayRead(u,&U);
3167: PetscDrawAxisGetLimits(ictx->axis,&xl,&xr,&yl,&yr);
3168: if ((PetscRealPart(U[0]) < xl) || (PetscRealPart(U[1]) < yl) || (PetscRealPart(U[0]) > xr) || (PetscRealPart(U[1]) > yr)) {
3169: VecRestoreArrayRead(u,&U);
3170: return(0);
3171: }
3172: if (!step) ictx->color++;
3173: PetscDrawPoint(draw,PetscRealPart(U[0]),PetscRealPart(U[1]),ictx->color);
3174: VecRestoreArrayRead(u,&U);
3176: if (ictx->showtimestepandtime) {
3177: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
3178: PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);
3179: PetscStrlen(time,&len);
3180: PetscDrawStringGetSize(draw,&tw,NULL);
3181: w = xl + .5*(xr - xl) - .5*len*tw;
3182: h = yl + .95*(yr - yl);
3183: PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);
3184: }
3185: PetscDrawFlush(draw);
3186: return(0);
3187: }
3192: /*@C
3193: TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution()
3195: Collective on TS
3197: Input Parameters:
3198: . ctx - the monitor context
3200: Level: intermediate
3202: .keywords: TS, vector, monitor, view
3204: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError()
3205: @*/
3206: PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx)
3207: {
3211: PetscDrawAxisDestroy(&(*ictx)->axis);
3212: PetscViewerDestroy(&(*ictx)->viewer);
3213: VecDestroy(&(*ictx)->initialsolution);
3214: PetscFree(*ictx);
3215: return(0);
3216: }
3220: /*@C
3221: TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx
3223: Collective on TS
3225: Input Parameter:
3226: . ts - time-step context
3228: Output Patameter:
3229: . ctx - the monitor context
3231: Options Database:
3232: . -ts_monitor_draw_solution_initial - show initial solution as well as current solution
3234: Level: intermediate
3236: .keywords: TS, vector, monitor, view
3238: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx()
3239: @*/
3240: PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx)
3241: {
3242: PetscErrorCode ierr;
3245: PetscNew(struct _n_TSMonitorDrawCtx,ctx);
3246: PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);
3247: PetscViewerSetFromOptions((*ctx)->viewer);
3249: (*ctx)->howoften = howoften;
3250: (*ctx)->showinitial = PETSC_FALSE;
3251: PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);
3253: (*ctx)->showtimestepandtime = PETSC_FALSE;
3254: PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);
3255: (*ctx)->color = PETSC_DRAW_WHITE;
3256: return(0);
3257: }
3261: /*@C
3262: TSMonitorDrawError - Monitors progress of the TS solvers by calling
3263: VecView() for the error at each timestep
3265: Collective on TS
3267: Input Parameters:
3268: + ts - the TS context
3269: . step - current time-step
3270: . ptime - current time
3271: - dummy - either a viewer or NULL
3273: Level: intermediate
3275: .keywords: TS, vector, monitor, view
3277: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
3278: @*/
3279: PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
3280: {
3281: PetscErrorCode ierr;
3282: TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy;
3283: PetscViewer viewer = ctx->viewer;
3284: Vec work;
3287: if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return(0);
3288: VecDuplicate(u,&work);
3289: TSComputeSolutionFunction(ts,ptime,work);
3290: VecAXPY(work,-1.0,u);
3291: VecView(work,viewer);
3292: VecDestroy(&work);
3293: return(0);
3294: }
3296: #include <petsc-private/dmimpl.h>
3299: /*@
3300: TSSetDM - Sets the DM that may be used by some preconditioners
3302: Logically Collective on TS and DM
3304: Input Parameters:
3305: + ts - the preconditioner context
3306: - dm - the dm
3308: Level: intermediate
3311: .seealso: TSGetDM(), SNESSetDM(), SNESGetDM()
3312: @*/
3313: PetscErrorCode TSSetDM(TS ts,DM dm)
3314: {
3316: SNES snes;
3317: DMTS tsdm;
3321: PetscObjectReference((PetscObject)dm);
3322: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
3323: if (ts->dm->dmts && !dm->dmts) {
3324: DMCopyDMTS(ts->dm,dm);
3325: DMGetDMTS(ts->dm,&tsdm);
3326: if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */
3327: tsdm->originaldm = dm;
3328: }
3329: }
3330: DMDestroy(&ts->dm);
3331: }
3332: ts->dm = dm;
3334: TSGetSNES(ts,&snes);
3335: SNESSetDM(snes,dm);
3336: return(0);
3337: }
3341: /*@
3342: TSGetDM - Gets the DM that may be used by some preconditioners
3344: Not Collective
3346: Input Parameter:
3347: . ts - the preconditioner context
3349: Output Parameter:
3350: . dm - the dm
3352: Level: intermediate
3355: .seealso: TSSetDM(), SNESSetDM(), SNESGetDM()
3356: @*/
3357: PetscErrorCode TSGetDM(TS ts,DM *dm)
3358: {
3363: if (!ts->dm) {
3364: DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);
3365: if (ts->snes) {SNESSetDM(ts->snes,ts->dm);}
3366: }
3367: *dm = ts->dm;
3368: return(0);
3369: }
3373: /*@
3374: SNESTSFormFunction - Function to evaluate nonlinear residual
3376: Logically Collective on SNES
3378: Input Parameter:
3379: + snes - nonlinear solver
3380: . U - the current state at which to evaluate the residual
3381: - ctx - user context, must be a TS
3383: Output Parameter:
3384: . F - the nonlinear residual
3386: Notes:
3387: This function is not normally called by users and is automatically registered with the SNES used by TS.
3388: It is most frequently passed to MatFDColoringSetFunction().
3390: Level: advanced
3392: .seealso: SNESSetFunction(), MatFDColoringSetFunction()
3393: @*/
3394: PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx)
3395: {
3396: TS ts = (TS)ctx;
3404: (ts->ops->snesfunction)(snes,U,F,ts);
3405: return(0);
3406: }
3410: /*@
3411: SNESTSFormJacobian - Function to evaluate the Jacobian
3413: Collective on SNES
3415: Input Parameter:
3416: + snes - nonlinear solver
3417: . U - the current state at which to evaluate the residual
3418: - ctx - user context, must be a TS
3420: Output Parameter:
3421: + A - the Jacobian
3422: . B - the preconditioning matrix (may be the same as A)
3423: - flag - indicates any structure change in the matrix
3425: Notes:
3426: This function is not normally called by users and is automatically registered with the SNES used by TS.
3428: Level: developer
3430: .seealso: SNESSetJacobian()
3431: @*/
3432: PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *flag,void *ctx)
3433: {
3434: TS ts = (TS)ctx;
3446: (ts->ops->snesjacobian)(snes,U,A,B,flag,ts);
3447: return(0);
3448: }
3452: /*@C
3453: TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only
3455: Collective on TS
3457: Input Arguments:
3458: + ts - time stepping context
3459: . t - time at which to evaluate
3460: . U - state at which to evaluate
3461: - ctx - context
3463: Output Arguments:
3464: . F - right hand side
3466: Level: intermediate
3468: Notes:
3469: This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems.
3470: The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian().
3472: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
3473: @*/
3474: PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
3475: {
3477: Mat Arhs,Brhs;
3478: MatStructure flg2;
3481: TSGetRHSMats_Private(ts,&Arhs,&Brhs);
3482: TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);
3483: MatMult(Arhs,U,F);
3484: return(0);
3485: }
3489: /*@C
3490: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
3492: Collective on TS
3494: Input Arguments:
3495: + ts - time stepping context
3496: . t - time at which to evaluate
3497: . U - state at which to evaluate
3498: - ctx - context
3500: Output Arguments:
3501: + A - pointer to operator
3502: . B - pointer to preconditioning matrix
3503: - flg - matrix structure flag
3505: Level: intermediate
3507: Notes:
3508: This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems.
3510: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear()
3511: @*/
3512: PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat *A,Mat *B,MatStructure *flg,void *ctx)
3513: {
3515: *flg = SAME_PRECONDITIONER;
3516: return(0);
3517: }
3521: /*@C
3522: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
3524: Collective on TS
3526: Input Arguments:
3527: + ts - time stepping context
3528: . t - time at which to evaluate
3529: . U - state at which to evaluate
3530: . Udot - time derivative of state vector
3531: - ctx - context
3533: Output Arguments:
3534: . F - left hand side
3536: Level: intermediate
3538: Notes:
3539: The assumption here is that the left hand side is of the form A*Udot (and not A*Udot + B*U). For other cases, the
3540: user is required to write their own TSComputeIFunction.
3541: This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems.
3542: The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian().
3544: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant()
3545: @*/
3546: PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
3547: {
3549: Mat A,B;
3550: MatStructure flg2;
3553: TSGetIJacobian(ts,&A,&B,NULL,NULL);
3554: TSComputeIJacobian(ts,t,U,Udot,1.0,&A,&B,&flg2,PETSC_TRUE);
3555: MatMult(A,Udot,F);
3556: return(0);
3557: }
3561: /*@C
3562: TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE
3564: Collective on TS
3566: Input Arguments:
3567: + ts - time stepping context
3568: . t - time at which to evaluate
3569: . U - state at which to evaluate
3570: . Udot - time derivative of state vector
3571: . shift - shift to apply
3572: - ctx - context
3574: Output Arguments:
3575: + A - pointer to operator
3576: . B - pointer to preconditioning matrix
3577: - flg - matrix structure flag
3579: Level: advanced
3581: Notes:
3582: This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems.
3584: It is only appropriate for problems of the form
3586: $ M Udot = F(U,t)
3588: where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only
3589: works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing
3590: an implicit operator of the form
3592: $ shift*M + J
3594: where J is the Jacobian of -F(U). Support may be added in a future version of PETSc, but for now, the user must store
3595: a copy of M or reassemble it when requested.
3597: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear()
3598: @*/
3599: PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx)
3600: {
3604: MatScale(*A, shift / ts->ijacobian.shift);
3605: ts->ijacobian.shift = shift;
3606: *flg = SAME_PRECONDITIONER;
3607: return(0);
3608: }
3612: /*@
3613: TSGetEquationType - Gets the type of the equation that TS is solving.
3615: Not Collective
3617: Input Parameter:
3618: . ts - the TS context
3620: Output Parameter:
3621: . equation_type - see TSEquationType
3623: Level: beginner
3625: .keywords: TS, equation type
3627: .seealso: TSSetEquationType(), TSEquationType
3628: @*/
3629: PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type)
3630: {
3634: *equation_type = ts->equation_type;
3635: return(0);
3636: }
3640: /*@
3641: TSSetEquationType - Sets the type of the equation that TS is solving.
3643: Not Collective
3645: Input Parameter:
3646: + ts - the TS context
3647: . equation_type - see TSEquationType
3649: Level: advanced
3651: .keywords: TS, equation type
3653: .seealso: TSGetEquationType(), TSEquationType
3654: @*/
3655: PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type)
3656: {
3659: ts->equation_type = equation_type;
3660: return(0);
3661: }
3665: /*@
3666: TSGetConvergedReason - Gets the reason the TS iteration was stopped.
3668: Not Collective
3670: Input Parameter:
3671: . ts - the TS context
3673: Output Parameter:
3674: . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
3675: manual pages for the individual convergence tests for complete lists
3677: Level: beginner
3679: Notes:
3680: Can only be called after the call to TSSolve() is complete.
3682: .keywords: TS, nonlinear, set, convergence, test
3684: .seealso: TSSetConvergenceTest(), TSConvergedReason
3685: @*/
3686: PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason)
3687: {
3691: *reason = ts->reason;
3692: return(0);
3693: }
3697: /*@
3698: TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve.
3700: Not Collective
3702: Input Parameter:
3703: + ts - the TS context
3704: . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
3705: manual pages for the individual convergence tests for complete lists
3707: Level: advanced
3709: Notes:
3710: Can only be called during TSSolve() is active.
3712: .keywords: TS, nonlinear, set, convergence, test
3714: .seealso: TSConvergedReason
3715: @*/
3716: PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason)
3717: {
3720: ts->reason = reason;
3721: return(0);
3722: }
3726: /*@
3727: TSGetSolveTime - Gets the time after a call to TSSolve()
3729: Not Collective
3731: Input Parameter:
3732: . ts - the TS context
3734: Output Parameter:
3735: . ftime - the final time. This time should correspond to the final time set with TSSetDuration()
3737: Level: beginner
3739: Notes:
3740: Can only be called after the call to TSSolve() is complete.
3742: .keywords: TS, nonlinear, set, convergence, test
3744: .seealso: TSSetConvergenceTest(), TSConvergedReason
3745: @*/
3746: PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime)
3747: {
3751: *ftime = ts->solvetime;
3752: return(0);
3753: }
3757: /*@
3758: TSGetSNESIterations - Gets the total number of nonlinear iterations
3759: used by the time integrator.
3761: Not Collective
3763: Input Parameter:
3764: . ts - TS context
3766: Output Parameter:
3767: . nits - number of nonlinear iterations
3769: Notes:
3770: This counter is reset to zero for each successive call to TSSolve().
3772: Level: intermediate
3774: .keywords: TS, get, number, nonlinear, iterations
3776: .seealso: TSGetKSPIterations()
3777: @*/
3778: PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits)
3779: {
3783: *nits = ts->snes_its;
3784: return(0);
3785: }
3789: /*@
3790: TSGetKSPIterations - Gets the total number of linear iterations
3791: used by the time integrator.
3793: Not Collective
3795: Input Parameter:
3796: . ts - TS context
3798: Output Parameter:
3799: . lits - number of linear iterations
3801: Notes:
3802: This counter is reset to zero for each successive call to TSSolve().
3804: Level: intermediate
3806: .keywords: TS, get, number, linear, iterations
3808: .seealso: TSGetSNESIterations(), SNESGetKSPIterations()
3809: @*/
3810: PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits)
3811: {
3815: *lits = ts->ksp_its;
3816: return(0);
3817: }
3821: /*@
3822: TSGetStepRejections - Gets the total number of rejected steps.
3824: Not Collective
3826: Input Parameter:
3827: . ts - TS context
3829: Output Parameter:
3830: . rejects - number of steps rejected
3832: Notes:
3833: This counter is reset to zero for each successive call to TSSolve().
3835: Level: intermediate
3837: .keywords: TS, get, number
3839: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails()
3840: @*/
3841: PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects)
3842: {
3846: *rejects = ts->reject;
3847: return(0);
3848: }
3852: /*@
3853: TSGetSNESFailures - Gets the total number of failed SNES solves
3855: Not Collective
3857: Input Parameter:
3858: . ts - TS context
3860: Output Parameter:
3861: . fails - number of failed nonlinear solves
3863: Notes:
3864: This counter is reset to zero for each successive call to TSSolve().
3866: Level: intermediate
3868: .keywords: TS, get, number
3870: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures()
3871: @*/
3872: PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails)
3873: {
3877: *fails = ts->num_snes_failures;
3878: return(0);
3879: }
3883: /*@
3884: TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails
3886: Not Collective
3888: Input Parameter:
3889: + ts - TS context
3890: - rejects - maximum number of rejected steps, pass -1 for unlimited
3892: Notes:
3893: The counter is reset to zero for each step
3895: Options Database Key:
3896: . -ts_max_reject - Maximum number of step rejections before a step fails
3898: Level: intermediate
3900: .keywords: TS, set, maximum, number
3902: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
3903: @*/
3904: PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects)
3905: {
3908: ts->max_reject = rejects;
3909: return(0);
3910: }
3914: /*@
3915: TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves
3917: Not Collective
3919: Input Parameter:
3920: + ts - TS context
3921: - fails - maximum number of failed nonlinear solves, pass -1 for unlimited
3923: Notes:
3924: The counter is reset to zero for each successive call to TSSolve().
3926: Options Database Key:
3927: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
3929: Level: intermediate
3931: .keywords: TS, set, maximum, number
3933: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason()
3934: @*/
3935: PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails)
3936: {
3939: ts->max_snes_failures = fails;
3940: return(0);
3941: }
3945: /*@
3946: TSSetErrorIfStepFails - Error if no step succeeds
3948: Not Collective
3950: Input Parameter:
3951: + ts - TS context
3952: - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure
3954: Options Database Key:
3955: . -ts_error_if_step_fails - Error if no step succeeds
3957: Level: intermediate
3959: .keywords: TS, set, error
3961: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
3962: @*/
3963: PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err)
3964: {
3967: ts->errorifstepfailed = err;
3968: return(0);
3969: }
3973: /*@C
3974: TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file
3976: Collective on TS
3978: Input Parameters:
3979: + ts - the TS context
3980: . step - current time-step
3981: . ptime - current time
3982: . u - current state
3983: - viewer - binary viewer
3985: Level: intermediate
3987: .keywords: TS, vector, monitor, view
3989: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
3990: @*/
3991: PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer)
3992: {
3994: PetscViewer v = (PetscViewer)viewer;
3997: VecView(u,v);
3998: return(0);
3999: }
4003: /*@C
4004: TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep.
4006: Collective on TS
4008: Input Parameters:
4009: + ts - the TS context
4010: . step - current time-step
4011: . ptime - current time
4012: . u - current state
4013: - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
4015: Level: intermediate
4017: Notes:
4018: The VTK format does not allow writing multiple time steps in the same file, therefore a different file will be written for each time step.
4019: These are named according to the file name template.
4021: This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy().
4023: .keywords: TS, vector, monitor, view
4025: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4026: @*/
4027: PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate)
4028: {
4030: char filename[PETSC_MAX_PATH_LEN];
4031: PetscViewer viewer;
4034: PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);
4035: PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);
4036: VecView(u,viewer);
4037: PetscViewerDestroy(&viewer);
4038: return(0);
4039: }
4043: /*@C
4044: TSMonitorSolutionVTKDestroy - Destroy context for monitoring
4046: Collective on TS
4048: Input Parameters:
4049: . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
4051: Level: intermediate
4053: Note:
4054: This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK().
4056: .keywords: TS, vector, monitor, view
4058: .seealso: TSMonitorSet(), TSMonitorSolutionVTK()
4059: @*/
4060: PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate)
4061: {
4065: PetscFree(*(char**)filenametemplate);
4066: return(0);
4067: }
4071: /*@
4072: TSGetAdapt - Get the adaptive controller context for the current method
4074: Collective on TS if controller has not been created yet
4076: Input Arguments:
4077: . ts - time stepping context
4079: Output Arguments:
4080: . adapt - adaptive controller
4082: Level: intermediate
4084: .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose()
4085: @*/
4086: PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt)
4087: {
4093: if (!ts->adapt) {
4094: TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);
4095: PetscLogObjectParent(ts,ts->adapt);
4096: PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);
4097: }
4098: *adapt = ts->adapt;
4099: return(0);
4100: }
4104: /*@
4105: TSSetTolerances - Set tolerances for local truncation error when using adaptive controller
4107: Logically Collective
4109: Input Arguments:
4110: + ts - time integration context
4111: . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value
4112: . vatol - vector of absolute tolerances or NULL, used in preference to atol if present
4113: . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value
4114: - vrtol - vector of relative tolerances or NULL, used in preference to atol if present
4116: Level: beginner
4118: .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances()
4119: @*/
4120: PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol)
4121: {
4125: if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol;
4126: if (vatol) {
4127: PetscObjectReference((PetscObject)vatol);
4128: VecDestroy(&ts->vatol);
4130: ts->vatol = vatol;
4131: }
4132: if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol;
4133: if (vrtol) {
4134: PetscObjectReference((PetscObject)vrtol);
4135: VecDestroy(&ts->vrtol);
4137: ts->vrtol = vrtol;
4138: }
4139: return(0);
4140: }
4144: /*@
4145: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
4147: Logically Collective
4149: Input Arguments:
4150: . ts - time integration context
4152: Output Arguments:
4153: + atol - scalar absolute tolerances, NULL to ignore
4154: . vatol - vector of absolute tolerances, NULL to ignore
4155: . rtol - scalar relative tolerances, NULL to ignore
4156: - vrtol - vector of relative tolerances, NULL to ignore
4158: Level: beginner
4160: .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances()
4161: @*/
4162: PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol)
4163: {
4165: if (atol) *atol = ts->atol;
4166: if (vatol) *vatol = ts->vatol;
4167: if (rtol) *rtol = ts->rtol;
4168: if (vrtol) *vrtol = ts->vrtol;
4169: return(0);
4170: }
4174: /*@
4175: TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state
4177: Collective on TS
4179: Input Arguments:
4180: + ts - time stepping context
4181: - Y - state vector to be compared to ts->vec_sol
4183: Output Arguments:
4184: . norm - weighted norm, a value of 1.0 is considered small
4186: Level: developer
4188: .seealso: TSSetTolerances()
4189: @*/
4190: PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm)
4191: {
4192: PetscErrorCode ierr;
4193: PetscInt i,n,N;
4194: const PetscScalar *u,*y;
4195: Vec U;
4196: PetscReal sum,gsum;
4202: U = ts->vec_sol;
4204: if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector");
4206: VecGetSize(U,&N);
4207: VecGetLocalSize(U,&n);
4208: VecGetArrayRead(U,&u);
4209: VecGetArrayRead(Y,&y);
4210: sum = 0.;
4211: if (ts->vatol && ts->vrtol) {
4212: const PetscScalar *atol,*rtol;
4213: VecGetArrayRead(ts->vatol,&atol);
4214: VecGetArrayRead(ts->vrtol,&rtol);
4215: for (i=0; i<n; i++) {
4216: PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
4217: sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol);
4218: }
4219: VecRestoreArrayRead(ts->vatol,&atol);
4220: VecRestoreArrayRead(ts->vrtol,&rtol);
4221: } else if (ts->vatol) { /* vector atol, scalar rtol */
4222: const PetscScalar *atol;
4223: VecGetArrayRead(ts->vatol,&atol);
4224: for (i=0; i<n; i++) {
4225: PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
4226: sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol);
4227: }
4228: VecRestoreArrayRead(ts->vatol,&atol);
4229: } else if (ts->vrtol) { /* scalar atol, vector rtol */
4230: const PetscScalar *rtol;
4231: VecGetArrayRead(ts->vrtol,&rtol);
4232: for (i=0; i<n; i++) {
4233: PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
4234: sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol);
4235: }
4236: VecRestoreArrayRead(ts->vrtol,&rtol);
4237: } else { /* scalar atol, scalar rtol */
4238: for (i=0; i<n; i++) {
4239: PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
4240: sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol);
4241: }
4242: }
4243: VecRestoreArrayRead(U,&u);
4244: VecRestoreArrayRead(Y,&y);
4246: MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));
4247: *norm = PetscSqrtReal(gsum / N);
4248: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
4249: return(0);
4250: }
4254: /*@
4255: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
4257: Logically Collective on TS
4259: Input Arguments:
4260: + ts - time stepping context
4261: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
4263: Note:
4264: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
4266: Level: intermediate
4268: .seealso: TSGetCFLTime(), TSADAPTCFL
4269: @*/
4270: PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime)
4271: {
4274: ts->cfltime_local = cfltime;
4275: ts->cfltime = -1.;
4276: return(0);
4277: }
4281: /*@
4282: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
4284: Collective on TS
4286: Input Arguments:
4287: . ts - time stepping context
4289: Output Arguments:
4290: . cfltime - maximum stable time step for forward Euler
4292: Level: advanced
4294: .seealso: TSSetCFLTimeLocal()
4295: @*/
4296: PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime)
4297: {
4301: if (ts->cfltime < 0) {
4302: MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));
4303: }
4304: *cfltime = ts->cfltime;
4305: return(0);
4306: }
4310: /*@
4311: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
4313: Input Parameters:
4314: . ts - the TS context.
4315: . xl - lower bound.
4316: . xu - upper bound.
4318: Notes:
4319: If this routine is not called then the lower and upper bounds are set to
4320: SNES_VI_NINF and SNES_VI_INF respectively during SNESSetUp().
4322: Level: advanced
4324: @*/
4325: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
4326: {
4328: SNES snes;
4331: TSGetSNES(ts,&snes);
4332: SNESVISetVariableBounds(snes,xl,xu);
4333: return(0);
4334: }
4336: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4337: #include <mex.h>
4339: typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext;
4343: /*
4344: TSComputeFunction_Matlab - Calls the function that has been set with
4345: TSSetFunctionMatlab().
4347: Collective on TS
4349: Input Parameters:
4350: + snes - the TS context
4351: - u - input vector
4353: Output Parameter:
4354: . y - function vector, as set by TSSetFunction()
4356: Notes:
4357: TSComputeFunction() is typically used within nonlinear solvers
4358: implementations, so most users would not generally call this routine
4359: themselves.
4361: Level: developer
4363: .keywords: TS, nonlinear, compute, function
4365: .seealso: TSSetFunction(), TSGetFunction()
4366: */
4367: PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx)
4368: {
4369: PetscErrorCode ierr;
4370: TSMatlabContext *sctx = (TSMatlabContext*)ctx;
4371: int nlhs = 1,nrhs = 7;
4372: mxArray *plhs[1],*prhs[7];
4373: long long int lx = 0,lxdot = 0,ly = 0,ls = 0;
4383: PetscMemcpy(&ls,&snes,sizeof(snes));
4384: PetscMemcpy(&lx,&u,sizeof(u));
4385: PetscMemcpy(&lxdot,&udot,sizeof(udot));
4386: PetscMemcpy(&ly,&y,sizeof(u));
4388: prhs[0] = mxCreateDoubleScalar((double)ls);
4389: prhs[1] = mxCreateDoubleScalar(time);
4390: prhs[2] = mxCreateDoubleScalar((double)lx);
4391: prhs[3] = mxCreateDoubleScalar((double)lxdot);
4392: prhs[4] = mxCreateDoubleScalar((double)ly);
4393: prhs[5] = mxCreateString(sctx->funcname);
4394: prhs[6] = sctx->ctx;
4395: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");
4396: mxGetScalar(plhs[0]);
4397: mxDestroyArray(prhs[0]);
4398: mxDestroyArray(prhs[1]);
4399: mxDestroyArray(prhs[2]);
4400: mxDestroyArray(prhs[3]);
4401: mxDestroyArray(prhs[4]);
4402: mxDestroyArray(prhs[5]);
4403: mxDestroyArray(plhs[0]);
4404: return(0);
4405: }
4410: /*
4411: TSSetFunctionMatlab - Sets the function evaluation routine and function
4412: vector for use by the TS routines in solving ODEs
4413: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
4415: Logically Collective on TS
4417: Input Parameters:
4418: + ts - the TS context
4419: - func - function evaluation routine
4421: Calling sequence of func:
4422: $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx);
4424: Level: beginner
4426: .keywords: TS, nonlinear, set, function
4428: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
4429: */
4430: PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx)
4431: {
4432: PetscErrorCode ierr;
4433: TSMatlabContext *sctx;
4436: /* currently sctx is memory bleed */
4437: PetscMalloc(sizeof(TSMatlabContext),&sctx);
4438: PetscStrallocpy(func,&sctx->funcname);
4439: /*
4440: This should work, but it doesn't
4441: sctx->ctx = ctx;
4442: mexMakeArrayPersistent(sctx->ctx);
4443: */
4444: sctx->ctx = mxDuplicateArray(ctx);
4446: TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);
4447: return(0);
4448: }
4452: /*
4453: TSComputeJacobian_Matlab - Calls the function that has been set with
4454: TSSetJacobianMatlab().
4456: Collective on TS
4458: Input Parameters:
4459: + ts - the TS context
4460: . u - input vector
4461: . A, B - the matrices
4462: - ctx - user context
4464: Output Parameter:
4465: . flag - structure of the matrix
4467: Level: developer
4469: .keywords: TS, nonlinear, compute, function
4471: .seealso: TSSetFunction(), TSGetFunction()
4472: @*/
4473: PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx)
4474: {
4475: PetscErrorCode ierr;
4476: TSMatlabContext *sctx = (TSMatlabContext*)ctx;
4477: int nlhs = 2,nrhs = 9;
4478: mxArray *plhs[2],*prhs[9];
4479: long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0;
4485: /* call Matlab function in ctx with arguments u and y */
4487: PetscMemcpy(&ls,&ts,sizeof(ts));
4488: PetscMemcpy(&lx,&u,sizeof(u));
4489: PetscMemcpy(&lxdot,&udot,sizeof(u));
4490: PetscMemcpy(&lA,A,sizeof(u));
4491: PetscMemcpy(&lB,B,sizeof(u));
4493: prhs[0] = mxCreateDoubleScalar((double)ls);
4494: prhs[1] = mxCreateDoubleScalar((double)time);
4495: prhs[2] = mxCreateDoubleScalar((double)lx);
4496: prhs[3] = mxCreateDoubleScalar((double)lxdot);
4497: prhs[4] = mxCreateDoubleScalar((double)shift);
4498: prhs[5] = mxCreateDoubleScalar((double)lA);
4499: prhs[6] = mxCreateDoubleScalar((double)lB);
4500: prhs[7] = mxCreateString(sctx->funcname);
4501: prhs[8] = sctx->ctx;
4502: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");
4503: mxGetScalar(plhs[0]);
4504: *flag = (MatStructure) mxGetScalar(plhs[1]);
4505: mxDestroyArray(prhs[0]);
4506: mxDestroyArray(prhs[1]);
4507: mxDestroyArray(prhs[2]);
4508: mxDestroyArray(prhs[3]);
4509: mxDestroyArray(prhs[4]);
4510: mxDestroyArray(prhs[5]);
4511: mxDestroyArray(prhs[6]);
4512: mxDestroyArray(prhs[7]);
4513: mxDestroyArray(plhs[0]);
4514: mxDestroyArray(plhs[1]);
4515: return(0);
4516: }
4521: /*
4522: TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
4523: vector for use by the TS routines in solving ODEs from MATLAB. Here the function is a string containing the name of a MATLAB function
4525: Logically Collective on TS
4527: Input Parameters:
4528: + ts - the TS context
4529: . A,B - Jacobian matrices
4530: . func - function evaluation routine
4531: - ctx - user context
4533: Calling sequence of func:
4534: $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx);
4537: Level: developer
4539: .keywords: TS, nonlinear, set, function
4541: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
4542: */
4543: PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx)
4544: {
4545: PetscErrorCode ierr;
4546: TSMatlabContext *sctx;
4549: /* currently sctx is memory bleed */
4550: PetscMalloc(sizeof(TSMatlabContext),&sctx);
4551: PetscStrallocpy(func,&sctx->funcname);
4552: /*
4553: This should work, but it doesn't
4554: sctx->ctx = ctx;
4555: mexMakeArrayPersistent(sctx->ctx);
4556: */
4557: sctx->ctx = mxDuplicateArray(ctx);
4559: TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);
4560: return(0);
4561: }
4565: /*
4566: TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab().
4568: Collective on TS
4570: .seealso: TSSetFunction(), TSGetFunction()
4571: @*/
4572: PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx)
4573: {
4574: PetscErrorCode ierr;
4575: TSMatlabContext *sctx = (TSMatlabContext*)ctx;
4576: int nlhs = 1,nrhs = 6;
4577: mxArray *plhs[1],*prhs[6];
4578: long long int lx = 0,ls = 0;
4584: PetscMemcpy(&ls,&ts,sizeof(ts));
4585: PetscMemcpy(&lx,&u,sizeof(u));
4587: prhs[0] = mxCreateDoubleScalar((double)ls);
4588: prhs[1] = mxCreateDoubleScalar((double)it);
4589: prhs[2] = mxCreateDoubleScalar((double)time);
4590: prhs[3] = mxCreateDoubleScalar((double)lx);
4591: prhs[4] = mxCreateString(sctx->funcname);
4592: prhs[5] = sctx->ctx;
4593: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");
4594: mxGetScalar(plhs[0]);
4595: mxDestroyArray(prhs[0]);
4596: mxDestroyArray(prhs[1]);
4597: mxDestroyArray(prhs[2]);
4598: mxDestroyArray(prhs[3]);
4599: mxDestroyArray(prhs[4]);
4600: mxDestroyArray(plhs[0]);
4601: return(0);
4602: }
4607: /*
4608: TSMonitorSetMatlab - Sets the monitor function from Matlab
4610: Level: developer
4612: .keywords: TS, nonlinear, set, function
4614: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
4615: */
4616: PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx)
4617: {
4618: PetscErrorCode ierr;
4619: TSMatlabContext *sctx;
4622: /* currently sctx is memory bleed */
4623: PetscMalloc(sizeof(TSMatlabContext),&sctx);
4624: PetscStrallocpy(func,&sctx->funcname);
4625: /*
4626: This should work, but it doesn't
4627: sctx->ctx = ctx;
4628: mexMakeArrayPersistent(sctx->ctx);
4629: */
4630: sctx->ctx = mxDuplicateArray(ctx);
4632: TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);
4633: return(0);
4634: }
4635: #endif
4641: /*@C
4642: TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector
4643: in a time based line graph
4645: Collective on TS
4647: Input Parameters:
4648: + ts - the TS context
4649: . step - current time-step
4650: . ptime - current time
4651: - lg - a line graph object
4653: Level: intermediate
4655: Notes: each process in a parallel run displays its component solutions in a separate window
4657: .keywords: TS, vector, monitor, view
4659: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4660: @*/
4661: PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4662: {
4663: PetscErrorCode ierr;
4664: TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy;
4665: const PetscScalar *yy;
4666: PetscInt dim;
4669: if (!step) {
4670: PetscDrawAxis axis;
4671: PetscDrawLGGetAxis(ctx->lg,&axis);
4672: PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");
4673: VecGetLocalSize(u,&dim);
4674: PetscDrawLGSetDimension(ctx->lg,dim);
4675: PetscDrawLGReset(ctx->lg);
4676: }
4677: VecGetArrayRead(u,&yy);
4678: #if defined(PETSC_USE_COMPLEX)
4679: {
4680: PetscReal *yreal;
4681: PetscInt i,n;
4682: VecGetLocalSize(u,&n);
4683: PetscMalloc(n*sizeof(PetscReal),&yreal);
4684: for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
4685: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
4686: PetscFree(yreal);
4687: }
4688: #else
4689: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
4690: #endif
4691: VecRestoreArrayRead(u,&yy);
4692: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
4693: PetscDrawLGDraw(ctx->lg);
4694: }
4695: return(0);
4696: }
4700: /*@C
4701: TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector
4702: in a time based line graph
4704: Collective on TS
4706: Input Parameters:
4707: + ts - the TS context
4708: . step - current time-step
4709: . ptime - current time
4710: - lg - a line graph object
4712: Level: intermediate
4714: Notes:
4715: Only for sequential solves.
4717: The user must provide the solution using TSSetSolutionFunction() to use this monitor.
4719: Options Database Keys:
4720: . -ts_monitor_lg_error - create a graphical monitor of error history
4722: .keywords: TS, vector, monitor, view
4724: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
4725: @*/
4726: PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4727: {
4728: PetscErrorCode ierr;
4729: TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy;
4730: const PetscScalar *yy;
4731: Vec y;
4732: PetscInt dim;
4735: if (!step) {
4736: PetscDrawAxis axis;
4737: PetscDrawLGGetAxis(ctx->lg,&axis);
4738: PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");
4739: VecGetLocalSize(u,&dim);
4740: PetscDrawLGSetDimension(ctx->lg,dim);
4741: PetscDrawLGReset(ctx->lg);
4742: }
4743: VecDuplicate(u,&y);
4744: TSComputeSolutionFunction(ts,ptime,y);
4745: VecAXPY(y,-1.0,u);
4746: VecGetArrayRead(y,&yy);
4747: #if defined(PETSC_USE_COMPLEX)
4748: {
4749: PetscReal *yreal;
4750: PetscInt i,n;
4751: VecGetLocalSize(y,&n);
4752: PetscMalloc(n*sizeof(PetscReal),&yreal);
4753: for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
4754: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
4755: PetscFree(yreal);
4756: }
4757: #else
4758: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
4759: #endif
4760: VecRestoreArrayRead(y,&yy);
4761: VecDestroy(&y);
4762: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
4763: PetscDrawLGDraw(ctx->lg);
4764: }
4765: return(0);
4766: }
4770: PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
4771: {
4772: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
4773: PetscReal x = ptime,y;
4775: PetscInt its;
4778: if (!n) {
4779: PetscDrawAxis axis;
4781: PetscDrawLGGetAxis(ctx->lg,&axis);
4782: PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");
4783: PetscDrawLGReset(ctx->lg);
4785: ctx->snes_its = 0;
4786: }
4787: TSGetSNESIterations(ts,&its);
4788: y = its - ctx->snes_its;
4789: PetscDrawLGAddPoint(ctx->lg,&x,&y);
4790: if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
4791: PetscDrawLGDraw(ctx->lg);
4792: }
4793: ctx->snes_its = its;
4794: return(0);
4795: }
4799: PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
4800: {
4801: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
4802: PetscReal x = ptime,y;
4804: PetscInt its;
4807: if (!n) {
4808: PetscDrawAxis axis;
4810: PetscDrawLGGetAxis(ctx->lg,&axis);
4811: PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");
4812: PetscDrawLGReset(ctx->lg);
4814: ctx->ksp_its = 0;
4815: }
4816: TSGetKSPIterations(ts,&its);
4817: y = its - ctx->ksp_its;
4818: PetscDrawLGAddPoint(ctx->lg,&x,&y);
4819: if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
4820: PetscDrawLGDraw(ctx->lg);
4821: }
4822: ctx->ksp_its = its;
4823: return(0);
4824: }
4828: /*@
4829: TSComputeLinearStability - computes the linear stability function at a point
4831: Collective on TS and Vec
4833: Input Parameters:
4834: + ts - the TS context
4835: - xr,xi - real and imaginary part of input arguments
4837: Output Parameters:
4838: . yr,yi - real and imaginary part of function value
4840: Level: developer
4842: .keywords: TS, compute
4844: .seealso: TSSetRHSFunction(), TSComputeIFunction()
4845: @*/
4846: PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
4847: {
4852: if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method");
4853: (*ts->ops->linearstability)(ts,xr,xi,yr,yi);
4854: return(0);
4855: }