Actual source code: ts.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/tsimpl.h>
2: #include <petscdmshell.h>
3: #include <petscdmda.h>
4: #include <petscviewer.h>
5: #include <petscdraw.h>
7: /* Logging support */
8: PetscClassId TS_CLASSID, DMTS_CLASSID;
9: PetscLogEvent TS_AdjointStep, TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;
11: const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED","STEPOVER","INTERPOLATE","MATCHSTEP","TSExactFinalTimeOption","TS_EXACTFINALTIME_",0};
13: struct _n_TSMonitorDrawCtx {
14: PetscViewer viewer;
15: Vec initialsolution;
16: PetscBool showinitial;
17: PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
18: PetscBool showtimestepandtime;
19: };
21: /*@C
22: TSMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
24: Collective on TS
26: Input Parameters:
27: + ts - TS object you wish to monitor
28: . name - the monitor type one is seeking
29: . help - message indicating what monitoring is done
30: . manual - manual page for the monitor
31: . monitor - the monitor function
32: - monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the TS or PetscViewer objects
34: Level: developer
36: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
37: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
38: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
39: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
40: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
41: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
42: PetscOptionsFList(), PetscOptionsEList()
43: @*/
44: PetscErrorCode TSMonitorSetFromOptions(TS ts,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(TS,PetscViewerAndFormat*))
45: {
46: PetscErrorCode ierr;
47: PetscViewer viewer;
48: PetscViewerFormat format;
49: PetscBool flg;
52: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
53: if (flg) {
54: PetscViewerAndFormat *vf;
55: PetscViewerAndFormatCreate(viewer,format,&vf);
56: PetscObjectDereference((PetscObject)viewer);
57: if (monitorsetup) {
58: (*monitorsetup)(ts,vf);
59: }
60: TSMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
61: }
62: return(0);
63: }
65: /*@C
66: TSAdjointMonitorSensi - monitors the first lambda sensitivity
68: Level: intermediate
70: .keywords: TS, set, monitor
72: .seealso: TSAdjointMonitorSet()
73: @*/
74: PetscErrorCode TSAdjointMonitorSensi(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
75: {
77: PetscViewer viewer = vf->viewer;
81: PetscViewerPushFormat(viewer,vf->format);
82: VecView(lambda[0],viewer);
83: PetscViewerPopFormat(viewer);
84: return(0);
85: }
87: /*@C
88: TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
90: Collective on TS
92: Input Parameters:
93: + ts - TS object you wish to monitor
94: . name - the monitor type one is seeking
95: . help - message indicating what monitoring is done
96: . manual - manual page for the monitor
97: . monitor - the monitor function
98: - monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the TS or PetscViewer objects
100: Level: developer
102: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
103: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
104: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
105: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
106: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
107: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
108: PetscOptionsFList(), PetscOptionsEList()
109: @*/
110: PetscErrorCode TSAdjointMonitorSetFromOptions(TS ts,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(TS,PetscViewerAndFormat*))
111: {
112: PetscErrorCode ierr;
113: PetscViewer viewer;
114: PetscViewerFormat format;
115: PetscBool flg;
118: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
119: if (flg) {
120: PetscViewerAndFormat *vf;
121: PetscViewerAndFormatCreate(viewer,format,&vf);
122: PetscObjectDereference((PetscObject)viewer);
123: if (monitorsetup) {
124: (*monitorsetup)(ts,vf);
125: }
126: TSAdjointMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
127: }
128: return(0);
129: }
131: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt,TSAdaptType default_type)
132: {
138: if (!((PetscObject)adapt)->type_name) {
139: TSAdaptSetType(adapt,default_type);
140: }
141: return(0);
142: }
144: /*@
145: TSSetFromOptions - Sets various TS parameters from user options.
147: Collective on TS
149: Input Parameter:
150: . ts - the TS context obtained from TSCreate()
152: Options Database Keys:
153: + -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSALPHA, TSGLLE, TSSSP, TSGLEE
154: . -ts_save_trajectory - checkpoint the solution at each time-step
155: . -ts_max_time <time> - maximum time to compute to
156: . -ts_max_steps <steps> - maximum number of time-steps to take
157: . -ts_init_time <time> - initial time to start computation
158: . -ts_final_time <time> - final time to compute to
159: . -ts_dt <dt> - initial time step
160: . -ts_exact_final_time <stepover,interpolate,matchstep> whether to stop at the exact given final time and how to compute the solution at that ti,e
161: . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
162: . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
163: . -ts_error_if_step_fails <true,false> - Error if no step succeeds
164: . -ts_rtol <rtol> - relative tolerance for local truncation error
165: . -ts_atol <atol> Absolute tolerance for local truncation error
166: . -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
167: . -ts_fd_color - Use finite differences with coloring to compute IJacobian
168: . -ts_monitor - print information at each timestep
169: . -ts_monitor_lg_solution - Monitor solution graphically
170: . -ts_monitor_lg_error - Monitor error graphically
171: . -ts_monitor_lg_timestep - Monitor timestep size graphically
172: . -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
173: . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
174: . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
175: . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
176: . -ts_monitor_draw_solution - Monitor solution graphically
177: . -ts_monitor_draw_solution_phase <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
178: . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
179: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
180: . -ts_monitor_solution_vtk <filename.vts,filename.vtu> - Save each time step to a binary file, use filename-%%03D.vts (filename-%%03D.vtu)
181: . -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
182: . -ts_adjoint_monitor - print information at each adjoint time step
183: - -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically
185: Developer Note: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
187: Level: beginner
189: .keywords: TS, timestep, set, options, database
191: .seealso: TSGetType()
192: @*/
193: PetscErrorCode TSSetFromOptions(TS ts)
194: {
195: PetscBool opt,flg,tflg;
196: PetscErrorCode ierr;
197: char monfilename[PETSC_MAX_PATH_LEN];
198: PetscReal time_step;
199: TSExactFinalTimeOption eftopt;
200: char dir[16];
201: TSIFunction ifun;
202: const char *defaultType;
203: char typeName[256];
208: TSRegisterAll();
209: TSGetIFunction(ts,NULL,&ifun,NULL);
211: PetscObjectOptionsBegin((PetscObject)ts);
212: if (((PetscObject)ts)->type_name)
213: defaultType = ((PetscObject)ts)->type_name;
214: else
215: defaultType = ifun ? TSBEULER : TSEULER;
216: PetscOptionsFList("-ts_type","TS method","TSSetType",TSList,defaultType,typeName,256,&opt);
217: if (opt) {
218: TSSetType(ts,typeName);
219: } else {
220: TSSetType(ts,defaultType);
221: }
223: /* Handle generic TS options */
224: PetscOptionsReal("-ts_max_time","Maximum time to run to","TSSetMaxTime",ts->max_time,&ts->max_time,NULL);
225: PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetMaxSteps",ts->max_steps,&ts->max_steps,NULL);
226: PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,NULL);
227: PetscOptionsReal("-ts_final_time","Final time to run to","TSSetMaxTime",ts->max_time,&ts->max_time,NULL);
228: PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&time_step,&flg);
229: if (flg) {TSSetTimeStep(ts,time_step);}
230: PetscOptionsEnum("-ts_exact_final_time","Option for handling of final time step","TSSetExactFinalTime",TSExactFinalTimeOptions,(PetscEnum)ts->exact_final_time,(PetscEnum*)&eftopt,&flg);
231: if (flg) {TSSetExactFinalTime(ts,eftopt);}
232: PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,NULL);
233: PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,NULL);
234: PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",ts->errorifstepfailed,&ts->errorifstepfailed,NULL);
235: PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,NULL);
236: PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,NULL);
238: #if defined(PETSC_HAVE_SAWS)
239: {
240: PetscBool set;
241: flg = PETSC_FALSE;
242: PetscOptionsBool("-ts_saws_block","Block for SAWs memory snooper at end of TSSolve","PetscObjectSAWsBlock",((PetscObject)ts)->amspublishblock,&flg,&set);
243: if (set) {
244: PetscObjectSAWsSetBlock((PetscObject)ts,flg);
245: }
246: }
247: #endif
249: /* Monitor options */
250: TSMonitorSetFromOptions(ts,"-ts_monitor","Monitor time and timestep size","TSMonitorDefault",TSMonitorDefault,NULL);
251: TSMonitorSetFromOptions(ts,"-ts_monitor_solution","View the solution at each timestep","TSMonitorSolution",TSMonitorSolution,NULL);
252: TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor","Monitor adjoint timestep size","TSAdjointMonitorDefault",TSAdjointMonitorDefault,NULL);
253: TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor_sensi","Monitor sensitivity in the adjoint computation","TSAdjointMonitorSensi",TSAdjointMonitorSensi,NULL);
255: PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
256: if (flg) {PetscPythonMonitorSet((PetscObject)ts,monfilename);}
258: PetscOptionsName("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",&opt);
259: if (opt) {
260: TSMonitorLGCtx ctx;
261: PetscInt howoften = 1;
263: PetscOptionsInt("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",howoften,&howoften,NULL);
264: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
265: TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
266: }
268: PetscOptionsName("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",&opt);
269: if (opt) {
270: TSMonitorLGCtx ctx;
271: PetscInt howoften = 1;
273: PetscOptionsInt("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",howoften,&howoften,NULL);
274: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
275: TSMonitorSet(ts,TSMonitorLGError,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
276: }
278: PetscOptionsName("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",&opt);
279: if (opt) {
280: TSMonitorLGCtx ctx;
281: PetscInt howoften = 1;
283: PetscOptionsInt("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);
284: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
285: TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
286: }
287: PetscOptionsName("-ts_monitor_lg_timestep_log","Monitor log timestep size graphically","TSMonitorLGTimeStep",&opt);
288: if (opt) {
289: TSMonitorLGCtx ctx;
290: PetscInt howoften = 1;
292: PetscOptionsInt("-ts_monitor_lg_timestep_log","Monitor log timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);
293: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
294: TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
295: ctx->semilogy = PETSC_TRUE;
296: }
298: PetscOptionsName("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",&opt);
299: if (opt) {
300: TSMonitorLGCtx ctx;
301: PetscInt howoften = 1;
303: PetscOptionsInt("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",howoften,&howoften,NULL);
304: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
305: TSMonitorSet(ts,TSMonitorLGSNESIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
306: }
307: PetscOptionsName("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",&opt);
308: if (opt) {
309: TSMonitorLGCtx ctx;
310: PetscInt howoften = 1;
312: PetscOptionsInt("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",howoften,&howoften,NULL);
313: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
314: TSMonitorSet(ts,TSMonitorLGKSPIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
315: }
316: PetscOptionsName("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",&opt);
317: if (opt) {
318: TSMonitorSPEigCtx ctx;
319: PetscInt howoften = 1;
321: PetscOptionsInt("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",howoften,&howoften,NULL);
322: TSMonitorSPEigCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
323: TSMonitorSet(ts,TSMonitorSPEig,ctx,(PetscErrorCode (*)(void**))TSMonitorSPEigCtxDestroy);
324: }
325: opt = PETSC_FALSE;
326: PetscOptionsName("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",&opt);
327: if (opt) {
328: TSMonitorDrawCtx ctx;
329: PetscInt howoften = 1;
331: PetscOptionsInt("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",howoften,&howoften,NULL);
332: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
333: TSMonitorSet(ts,TSMonitorDrawSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
334: }
335: opt = PETSC_FALSE;
336: PetscOptionsName("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",&opt);
337: if (opt) {
338: TSMonitorDrawCtx ctx;
339: PetscInt howoften = 1;
341: PetscOptionsInt("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",howoften,&howoften,NULL);
342: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
343: TSAdjointMonitorSet(ts,TSAdjointMonitorDrawSensi,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
344: }
345: opt = PETSC_FALSE;
346: PetscOptionsName("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",&opt);
347: if (opt) {
348: TSMonitorDrawCtx ctx;
349: PetscReal bounds[4];
350: PetscInt n = 4;
351: PetscDraw draw;
352: PetscDrawAxis axis;
354: PetscOptionsRealArray("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",bounds,&n,NULL);
355: if (n != 4) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Must provide bounding box of phase field");
356: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,1,&ctx);
357: PetscViewerDrawGetDraw(ctx->viewer,0,&draw);
358: PetscViewerDrawGetDrawAxis(ctx->viewer,0,&axis);
359: PetscDrawAxisSetLimits(axis,bounds[0],bounds[2],bounds[1],bounds[3]);
360: PetscDrawAxisSetLabels(axis,"Phase Diagram","Variable 1","Variable 2");
361: TSMonitorSet(ts,TSMonitorDrawSolutionPhase,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
362: }
363: opt = PETSC_FALSE;
364: PetscOptionsName("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",&opt);
365: if (opt) {
366: TSMonitorDrawCtx ctx;
367: PetscInt howoften = 1;
369: PetscOptionsInt("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",howoften,&howoften,NULL);
370: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
371: TSMonitorSet(ts,TSMonitorDrawError,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
372: }
374: opt = PETSC_FALSE;
375: 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);
376: if (flg) {
377: const char *ptr,*ptr2;
378: char *filetemplate;
379: if (!monfilename[0]) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
380: /* Do some cursory validation of the input. */
381: PetscStrstr(monfilename,"%",(char**)&ptr);
382: if (!ptr) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
383: for (ptr++; ptr && *ptr; ptr++) {
384: PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
385: 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");
386: if (ptr2) break;
387: }
388: PetscStrallocpy(monfilename,&filetemplate);
389: TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);
390: }
392: PetscOptionsString("-ts_monitor_dmda_ray","Display a ray of the solution","None","y=0",dir,16,&flg);
393: if (flg) {
394: TSMonitorDMDARayCtx *rayctx;
395: int ray = 0;
396: DMDADirection ddir;
397: DM da;
398: PetscMPIInt rank;
400: if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
401: if (dir[0] == 'x') ddir = DMDA_X;
402: else if (dir[0] == 'y') ddir = DMDA_Y;
403: else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
404: sscanf(dir+2,"%d",&ray);
406: PetscInfo2(((PetscObject)ts),"Displaying DMDA ray %c = %D\n",dir[0],ray);
407: PetscNew(&rayctx);
408: TSGetDM(ts,&da);
409: DMDAGetRay(da,ddir,ray,&rayctx->ray,&rayctx->scatter);
410: MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);
411: if (!rank) {
412: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,600,300,&rayctx->viewer);
413: }
414: rayctx->lgctx = NULL;
415: TSMonitorSet(ts,TSMonitorDMDARay,rayctx,TSMonitorDMDARayDestroy);
416: }
417: PetscOptionsString("-ts_monitor_lg_dmda_ray","Display a ray of the solution","None","x=0",dir,16,&flg);
418: if (flg) {
419: TSMonitorDMDARayCtx *rayctx;
420: int ray = 0;
421: DMDADirection ddir;
422: DM da;
423: PetscInt howoften = 1;
425: if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
426: if (dir[0] == 'x') ddir = DMDA_X;
427: else if (dir[0] == 'y') ddir = DMDA_Y;
428: else SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
429: sscanf(dir+2, "%d", &ray);
431: PetscInfo2(((PetscObject) ts),"Displaying LG DMDA ray %c = %D\n", dir[0], ray);
432: PetscNew(&rayctx);
433: TSGetDM(ts, &da);
434: DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);
435: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&rayctx->lgctx);
436: TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);
437: }
439: PetscOptionsName("-ts_monitor_envelope","Monitor maximum and minimum value of each component of the solution","TSMonitorEnvelope",&opt);
440: if (opt) {
441: TSMonitorEnvelopeCtx ctx;
443: TSMonitorEnvelopeCtxCreate(ts,&ctx);
444: TSMonitorSet(ts,TSMonitorEnvelope,ctx,(PetscErrorCode (*)(void**))TSMonitorEnvelopeCtxDestroy);
445: }
447: flg = PETSC_FALSE;
448: PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeJacobianDefaultColor", flg, &flg, NULL);
449: if (flg) {
450: DM dm;
451: DMTS tdm;
453: TSGetDM(ts, &dm);
454: DMGetDMTS(dm, &tdm);
455: tdm->ijacobianctx = NULL;
456: TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, 0);
457: PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n");
458: }
460: /* Handle specific TS options */
461: if (ts->ops->setfromoptions) {
462: (*ts->ops->setfromoptions)(PetscOptionsObject,ts);
463: }
465: /* Handle TSAdapt options */
466: TSGetAdapt(ts,&ts->adapt);
467: TSAdaptSetDefaultType(ts->adapt,ts->default_adapt_type);
468: TSAdaptSetFromOptions(PetscOptionsObject,ts->adapt);
470: /* TS trajectory must be set after TS, since it may use some TS options above */
471: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
472: PetscOptionsBool("-ts_save_trajectory","Save the solution at each timestep","TSSetSaveTrajectory",tflg,&tflg,NULL);
473: if (tflg) {
474: TSSetSaveTrajectory(ts);
475: }
476: tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
477: PetscOptionsBool("-ts_adjoint_solve","Solve the adjoint problem immediately after solving the forward problem","",tflg,&tflg,&flg);
478: if (flg) {
479: TSSetSaveTrajectory(ts);
480: ts->adjoint_solve = tflg;
481: }
483: /* process any options handlers added with PetscObjectAddOptionsHandler() */
484: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)ts);
485: PetscOptionsEnd();
487: if (ts->trajectory) {
488: TSTrajectorySetFromOptions(ts->trajectory,ts);
489: }
491: TSGetSNES(ts,&ts->snes);
492: if (ts->problem_type == TS_LINEAR) {SNESSetType(ts->snes,SNESKSPONLY);}
493: SNESSetFromOptions(ts->snes);
494: return(0);
495: }
497: /*@
498: TSGetTrajectory - Gets the trajectory from a TS if it exists
500: Collective on TS
502: Input Parameters:
503: . ts - the TS context obtained from TSCreate()
505: Output Parameters;
506: . tr - the TSTrajectory object, if it exists
508: Note: This routine should be called after all TS options have been set
510: Level: advanced
512: .seealso: TSGetTrajectory(), TSAdjointSolve(), TSTrajectory, TSTrajectoryCreate()
514: .keywords: TS, set, checkpoint,
515: @*/
516: PetscErrorCode TSGetTrajectory(TS ts,TSTrajectory *tr)
517: {
520: *tr = ts->trajectory;
521: return(0);
522: }
524: /*@
525: TSSetSaveTrajectory - Causes the TS to save its solutions as it iterates forward in time in a TSTrajectory object
527: Collective on TS
529: Input Parameters:
530: . ts - the TS context obtained from TSCreate()
532: Options Database:
533: + -ts_save_trajectory - saves the trajectory to a file
534: - -ts_trajectory_type type
536: Note: This routine should be called after all TS options have been set
538: The TSTRAJECTORYVISUALIZATION files can be loaded into Python with $PETSC_DIR/bin/PetscBinaryIOTrajectory.py and
539: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
541: Level: intermediate
543: .seealso: TSGetTrajectory(), TSAdjointSolve(), TSTrajectoryType, TSSetTrajectoryType()
545: .keywords: TS, set, checkpoint,
546: @*/
547: PetscErrorCode TSSetSaveTrajectory(TS ts)
548: {
553: if (!ts->trajectory) {
554: TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);
555: TSTrajectorySetFromOptions(ts->trajectory,ts);
556: }
557: return(0);
558: }
560: /*@
561: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
562: set with TSSetRHSJacobian().
564: Collective on TS and Vec
566: Input Parameters:
567: + ts - the TS context
568: . t - current timestep
569: - U - input vector
571: Output Parameters:
572: + A - Jacobian matrix
573: . B - optional preconditioning matrix
574: - flag - flag indicating matrix structure
576: Notes:
577: Most users should not need to explicitly call this routine, as it
578: is used internally within the nonlinear solvers.
580: See KSPSetOperators() for important information about setting the
581: flag parameter.
583: Level: developer
585: .keywords: SNES, compute, Jacobian, matrix
587: .seealso: TSSetRHSJacobian(), KSPSetOperators()
588: @*/
589: PetscErrorCode TSComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B)
590: {
591: PetscErrorCode ierr;
592: PetscObjectState Ustate;
593: PetscObjectId Uid;
594: DM dm;
595: DMTS tsdm;
596: TSRHSJacobian rhsjacobianfunc;
597: void *ctx;
598: TSIJacobian ijacobianfunc;
599: TSRHSFunction rhsfunction;
605: TSGetDM(ts,&dm);
606: DMGetDMTS(dm,&tsdm);
607: DMTSGetRHSJacobian(dm,&rhsjacobianfunc,&ctx);
608: DMTSGetIJacobian(dm,&ijacobianfunc,NULL);
609: DMTSGetRHSFunction(dm,&rhsfunction,&ctx);
610: PetscObjectStateGet((PetscObject)U,&Ustate);
611: PetscObjectGetId((PetscObject)U,&Uid);
612: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) {
613: return(0);
614: }
616: if (!rhsjacobianfunc && !ijacobianfunc) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
618: if (ts->rhsjacobian.reuse) {
619: MatShift(A,-ts->rhsjacobian.shift);
620: MatScale(A,1./ts->rhsjacobian.scale);
621: if (B && A != B) {
622: MatShift(B,-ts->rhsjacobian.shift);
623: MatScale(B,1./ts->rhsjacobian.scale);
624: }
625: ts->rhsjacobian.shift = 0;
626: ts->rhsjacobian.scale = 1.;
627: }
629: if (rhsjacobianfunc) {
630: PetscBool missing;
631: PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);
632: PetscStackPush("TS user Jacobian function");
633: (*rhsjacobianfunc)(ts,t,U,A,B,ctx);
634: PetscStackPop;
635: PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);
636: if (A) {
637: MatMissingDiagonal(A,&missing,NULL);
638: if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Amat passed to TSSetRHSJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
639: }
640: if (B && B != A) {
641: MatMissingDiagonal(B,&missing,NULL);
642: if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Bmat passed to TSSetRHSJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
643: }
644: } else {
645: MatZeroEntries(A);
646: if (A != B) {MatZeroEntries(B);}
647: }
648: ts->rhsjacobian.time = t;
649: PetscObjectGetId((PetscObject)U,&ts->rhsjacobian.Xid);
650: PetscObjectStateGet((PetscObject)U,&ts->rhsjacobian.Xstate);
651: return(0);
652: }
654: /*@
655: TSComputeRHSFunction - Evaluates the right-hand-side function.
657: Collective on TS and Vec
659: Input Parameters:
660: + ts - the TS context
661: . t - current time
662: - U - state vector
664: Output Parameter:
665: . y - right hand side
667: Note:
668: Most users should not need to explicitly call this routine, as it
669: is used internally within the nonlinear solvers.
671: Level: developer
673: .keywords: TS, compute
675: .seealso: TSSetRHSFunction(), TSComputeIFunction()
676: @*/
677: PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y)
678: {
680: TSRHSFunction rhsfunction;
681: TSIFunction ifunction;
682: void *ctx;
683: DM dm;
689: TSGetDM(ts,&dm);
690: DMTSGetRHSFunction(dm,&rhsfunction,&ctx);
691: DMTSGetIFunction(dm,&ifunction,NULL);
693: if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
695: PetscLogEventBegin(TS_FunctionEval,ts,U,y,0);
696: if (rhsfunction) {
697: PetscStackPush("TS user right-hand-side function");
698: (*rhsfunction)(ts,t,U,y,ctx);
699: PetscStackPop;
700: } else {
701: VecZeroEntries(y);
702: }
704: PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);
705: return(0);
706: }
708: /*@
709: TSComputeSolutionFunction - Evaluates the solution function.
711: Collective on TS and Vec
713: Input Parameters:
714: + ts - the TS context
715: - t - current time
717: Output Parameter:
718: . U - the solution
720: Note:
721: Most users should not need to explicitly call this routine, as it
722: is used internally within the nonlinear solvers.
724: Level: developer
726: .keywords: TS, compute
728: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
729: @*/
730: PetscErrorCode TSComputeSolutionFunction(TS ts,PetscReal t,Vec U)
731: {
732: PetscErrorCode ierr;
733: TSSolutionFunction solutionfunction;
734: void *ctx;
735: DM dm;
740: TSGetDM(ts,&dm);
741: DMTSGetSolutionFunction(dm,&solutionfunction,&ctx);
743: if (solutionfunction) {
744: PetscStackPush("TS user solution function");
745: (*solutionfunction)(ts,t,U,ctx);
746: PetscStackPop;
747: }
748: return(0);
749: }
750: /*@
751: TSComputeForcingFunction - Evaluates the forcing function.
753: Collective on TS and Vec
755: Input Parameters:
756: + ts - the TS context
757: - t - current time
759: Output Parameter:
760: . U - the function value
762: Note:
763: Most users should not need to explicitly call this routine, as it
764: is used internally within the nonlinear solvers.
766: Level: developer
768: .keywords: TS, compute
770: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
771: @*/
772: PetscErrorCode TSComputeForcingFunction(TS ts,PetscReal t,Vec U)
773: {
774: PetscErrorCode ierr, (*forcing)(TS,PetscReal,Vec,void*);
775: void *ctx;
776: DM dm;
781: TSGetDM(ts,&dm);
782: DMTSGetForcingFunction(dm,&forcing,&ctx);
784: if (forcing) {
785: PetscStackPush("TS user forcing function");
786: (*forcing)(ts,t,U,ctx);
787: PetscStackPop;
788: }
789: return(0);
790: }
792: static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs)
793: {
794: Vec F;
798: *Frhs = NULL;
799: TSGetIFunction(ts,&F,NULL,NULL);
800: if (!ts->Frhs) {
801: VecDuplicate(F,&ts->Frhs);
802: }
803: *Frhs = ts->Frhs;
804: return(0);
805: }
807: static PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs)
808: {
809: Mat A,B;
811: TSIJacobian ijacobian;
814: if (Arhs) *Arhs = NULL;
815: if (Brhs) *Brhs = NULL;
816: TSGetIJacobian(ts,&A,&B,&ijacobian,NULL);
817: if (Arhs) {
818: if (!ts->Arhs) {
819: if (ijacobian) {
820: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);
821: } else {
822: ts->Arhs = A;
823: PetscObjectReference((PetscObject)A);
824: }
825: } else {
826: PetscBool flg;
827: SNESGetUseMatrixFree(ts->snes,NULL,&flg);
828: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
829: if (flg && !ijacobian && ts->Arhs == ts->Brhs){
830: PetscObjectDereference((PetscObject)ts->Arhs);
831: ts->Arhs = A;
832: PetscObjectReference((PetscObject)A);
833: }
834: }
835: *Arhs = ts->Arhs;
836: }
837: if (Brhs) {
838: if (!ts->Brhs) {
839: if (A != B) {
840: if (ijacobian) {
841: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);
842: } else {
843: ts->Brhs = B;
844: PetscObjectReference((PetscObject)B);
845: }
846: } else {
847: PetscObjectReference((PetscObject)ts->Arhs);
848: ts->Brhs = ts->Arhs;
849: }
850: }
851: *Brhs = ts->Brhs;
852: }
853: return(0);
854: }
856: /*@
857: TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0
859: Collective on TS and Vec
861: Input Parameters:
862: + ts - the TS context
863: . t - current time
864: . U - state vector
865: . Udot - time derivative of state vector
866: - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate
868: Output Parameter:
869: . Y - right hand side
871: Note:
872: Most users should not need to explicitly call this routine, as it
873: is used internally within the nonlinear solvers.
875: If the user did did not write their equations in implicit form, this
876: function recasts them in implicit form.
878: Level: developer
880: .keywords: TS, compute
882: .seealso: TSSetIFunction(), TSComputeRHSFunction()
883: @*/
884: PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec Y,PetscBool imex)
885: {
887: TSIFunction ifunction;
888: TSRHSFunction rhsfunction;
889: void *ctx;
890: DM dm;
898: TSGetDM(ts,&dm);
899: DMTSGetIFunction(dm,&ifunction,&ctx);
900: DMTSGetRHSFunction(dm,&rhsfunction,NULL);
902: if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
904: PetscLogEventBegin(TS_FunctionEval,ts,U,Udot,Y);
905: if (ifunction) {
906: PetscStackPush("TS user implicit function");
907: (*ifunction)(ts,t,U,Udot,Y,ctx);
908: PetscStackPop;
909: }
910: if (imex) {
911: if (!ifunction) {
912: VecCopy(Udot,Y);
913: }
914: } else if (rhsfunction) {
915: if (ifunction) {
916: Vec Frhs;
917: TSGetRHSVec_Private(ts,&Frhs);
918: TSComputeRHSFunction(ts,t,U,Frhs);
919: VecAXPY(Y,-1,Frhs);
920: } else {
921: TSComputeRHSFunction(ts,t,U,Y);
922: VecAYPX(Y,-1,Udot);
923: }
924: }
925: PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);
926: return(0);
927: }
929: /*@
930: TSComputeIJacobian - Evaluates the Jacobian of the DAE
932: Collective on TS and Vec
934: Input
935: Input Parameters:
936: + ts - the TS context
937: . t - current timestep
938: . U - state vector
939: . Udot - time derivative of state vector
940: . shift - shift to apply, see note below
941: - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate
943: Output Parameters:
944: + A - Jacobian matrix
945: - B - matrix from which the preconditioner is constructed; often the same as A
947: Notes:
948: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
950: dF/dU + shift*dF/dUdot
952: Most users should not need to explicitly call this routine, as it
953: is used internally within the nonlinear solvers.
955: Level: developer
957: .keywords: TS, compute, Jacobian, matrix
959: .seealso: TSSetIJacobian()
960: @*/
961: PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,PetscBool imex)
962: {
964: TSIJacobian ijacobian;
965: TSRHSJacobian rhsjacobian;
966: DM dm;
967: void *ctx;
978: TSGetDM(ts,&dm);
979: DMTSGetIJacobian(dm,&ijacobian,&ctx);
980: DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);
982: if (!rhsjacobian && !ijacobian) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
984: PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);
985: if (ijacobian) {
986: PetscBool missing;
987: PetscStackPush("TS user implicit Jacobian");
988: (*ijacobian)(ts,t,U,Udot,shift,A,B,ctx);
989: PetscStackPop;
990: MatMissingDiagonal(A,&missing,NULL);
991: if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Amat passed to TSSetIJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
992: if (B != A) {
993: MatMissingDiagonal(B,&missing,NULL);
994: if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Bmat passed to TSSetIJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
995: }
996: }
997: if (imex) {
998: if (!ijacobian) { /* system was written as Udot = G(t,U) */
999: PetscBool assembled;
1000: MatZeroEntries(A);
1001: MatAssembled(A,&assembled);
1002: if (!assembled) {
1003: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1004: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1005: }
1006: MatShift(A,shift);
1007: if (A != B) {
1008: MatZeroEntries(B);
1009: MatAssembled(B,&assembled);
1010: if (!assembled) {
1011: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1012: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1013: }
1014: MatShift(B,shift);
1015: }
1016: }
1017: } else {
1018: Mat Arhs = NULL,Brhs = NULL;
1019: if (rhsjacobian) {
1020: TSGetRHSMats_Private(ts,&Arhs,&Brhs);
1021: TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);
1022: }
1023: if (Arhs == A) { /* No IJacobian, so we only have the RHS matrix */
1024: PetscBool flg;
1025: ts->rhsjacobian.scale = -1;
1026: ts->rhsjacobian.shift = shift;
1027: SNESGetUseMatrixFree(ts->snes,NULL,&flg);
1028: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
1029: if (!flg) {
1030: MatScale(A,-1);
1031: MatShift(A,shift);
1032: }
1033: if (A != B) {
1034: MatScale(B,-1);
1035: MatShift(B,shift);
1036: }
1037: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
1038: MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
1039: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
1040: MatZeroEntries(A);
1041: MatShift(A,shift);
1042: if (A != B) {
1043: MatZeroEntries(B);
1044: MatShift(B,shift);
1045: }
1046: }
1047: MatAXPY(A,-1,Arhs,axpy);
1048: if (A != B) {
1049: MatAXPY(B,-1,Brhs,axpy);
1050: }
1051: }
1052: }
1053: PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);
1054: return(0);
1055: }
1057: /*@C
1058: TSSetRHSFunction - Sets the routine for evaluating the function,
1059: where U_t = G(t,u).
1061: Logically Collective on TS
1063: Input Parameters:
1064: + ts - the TS context obtained from TSCreate()
1065: . r - vector to put the computed right hand side (or NULL to have it created)
1066: . f - routine for evaluating the right-hand-side function
1067: - ctx - [optional] user-defined context for private data for the
1068: function evaluation routine (may be NULL)
1070: Calling sequence of func:
1071: $ func (TS ts,PetscReal t,Vec u,Vec F,void *ctx);
1073: + t - current timestep
1074: . u - input vector
1075: . F - function vector
1076: - ctx - [optional] user-defined function context
1078: Level: beginner
1080: Notes: You must call this function or TSSetIFunction() to define your ODE. You cannot use this function when solving a DAE.
1082: .keywords: TS, timestep, set, right-hand-side, function
1084: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSSetIFunction()
1085: @*/
1086: PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx)
1087: {
1089: SNES snes;
1090: Vec ralloc = NULL;
1091: DM dm;
1097: TSGetDM(ts,&dm);
1098: DMTSSetRHSFunction(dm,f,ctx);
1099: TSGetSNES(ts,&snes);
1100: if (!r && !ts->dm && ts->vec_sol) {
1101: VecDuplicate(ts->vec_sol,&ralloc);
1102: r = ralloc;
1103: }
1104: SNESSetFunction(snes,r,SNESTSFormFunction,ts);
1105: VecDestroy(&ralloc);
1106: return(0);
1107: }
1109: /*@C
1110: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1112: Logically Collective on TS
1114: Input Parameters:
1115: + ts - the TS context obtained from TSCreate()
1116: . f - routine for evaluating the solution
1117: - ctx - [optional] user-defined context for private data for the
1118: function evaluation routine (may be NULL)
1120: Calling sequence of func:
1121: $ func (TS ts,PetscReal t,Vec u,void *ctx);
1123: + t - current timestep
1124: . u - output vector
1125: - ctx - [optional] user-defined function context
1127: Notes:
1128: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1129: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1130: create closed-form solutions with non-physical forcing terms.
1132: For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history.
1134: Level: beginner
1136: .keywords: TS, timestep, set, right-hand-side, function
1138: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetForcingFunction()
1139: @*/
1140: PetscErrorCode TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx)
1141: {
1143: DM dm;
1147: TSGetDM(ts,&dm);
1148: DMTSSetSolutionFunction(dm,f,ctx);
1149: return(0);
1150: }
1152: /*@C
1153: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1155: Logically Collective on TS
1157: Input Parameters:
1158: + ts - the TS context obtained from TSCreate()
1159: . func - routine for evaluating the forcing function
1160: - ctx - [optional] user-defined context for private data for the
1161: function evaluation routine (may be NULL)
1163: Calling sequence of func:
1164: $ func (TS ts,PetscReal t,Vec f,void *ctx);
1166: + t - current timestep
1167: . f - output vector
1168: - ctx - [optional] user-defined function context
1170: Notes:
1171: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1172: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1173: definition of the problem you are solving and hence possibly introducing bugs.
1175: This replaces the ODE F(u,u_t,t) = 0 the TS is solving with F(u,u_t,t) - func(t) = 0
1177: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1178: parameters can be passed in the ctx variable.
1180: For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history.
1182: Level: beginner
1184: .keywords: TS, timestep, set, right-hand-side, function
1186: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetSolutionFunction()
1187: @*/
1188: PetscErrorCode TSSetForcingFunction(TS ts,TSForcingFunction func,void *ctx)
1189: {
1191: DM dm;
1195: TSGetDM(ts,&dm);
1196: DMTSSetForcingFunction(dm,func,ctx);
1197: return(0);
1198: }
1200: /*@C
1201: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1202: where U_t = G(U,t), as well as the location to store the matrix.
1204: Logically Collective on TS
1206: Input Parameters:
1207: + ts - the TS context obtained from TSCreate()
1208: . Amat - (approximate) Jacobian matrix
1209: . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1210: . f - the Jacobian evaluation routine
1211: - ctx - [optional] user-defined context for private data for the
1212: Jacobian evaluation routine (may be NULL)
1214: Calling sequence of f:
1215: $ func (TS ts,PetscReal t,Vec u,Mat A,Mat B,void *ctx);
1217: + t - current timestep
1218: . u - input vector
1219: . Amat - (approximate) Jacobian matrix
1220: . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1221: - ctx - [optional] user-defined context for matrix evaluation routine
1223: Notes:
1224: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1226: The TS solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1227: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1229: Level: beginner
1231: .keywords: TS, timestep, set, right-hand-side, Jacobian
1233: .seealso: SNESComputeJacobianDefaultColor(), TSSetRHSFunction(), TSRHSJacobianSetReuse(), TSSetIJacobian()
1235: @*/
1236: PetscErrorCode TSSetRHSJacobian(TS ts,Mat Amat,Mat Pmat,TSRHSJacobian f,void *ctx)
1237: {
1239: SNES snes;
1240: DM dm;
1241: TSIJacobian ijacobian;
1250: TSGetDM(ts,&dm);
1251: DMTSSetRHSJacobian(dm,f,ctx);
1252: if (f == TSComputeRHSJacobianConstant) {
1253: /* Handle this case automatically for the user; otherwise user should call themselves. */
1254: TSRHSJacobianSetReuse(ts,PETSC_TRUE);
1255: }
1256: DMTSGetIJacobian(dm,&ijacobian,NULL);
1257: TSGetSNES(ts,&snes);
1258: if (!ijacobian) {
1259: SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1260: }
1261: if (Amat) {
1262: PetscObjectReference((PetscObject)Amat);
1263: MatDestroy(&ts->Arhs);
1264: ts->Arhs = Amat;
1265: }
1266: if (Pmat) {
1267: PetscObjectReference((PetscObject)Pmat);
1268: MatDestroy(&ts->Brhs);
1269: ts->Brhs = Pmat;
1270: }
1271: return(0);
1272: }
1274: /*@C
1275: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1277: Logically Collective on TS
1279: Input Parameters:
1280: + ts - the TS context obtained from TSCreate()
1281: . r - vector to hold the residual (or NULL to have it created internally)
1282: . f - the function evaluation routine
1283: - ctx - user-defined context for private data for the function evaluation routine (may be NULL)
1285: Calling sequence of f:
1286: $ f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx);
1288: + t - time at step/stage being solved
1289: . u - state vector
1290: . u_t - time derivative of state vector
1291: . F - function vector
1292: - ctx - [optional] user-defined context for matrix evaluation routine
1294: Important:
1295: The user MUST call either this routine or TSSetRHSFunction() to define the ODE. When solving DAEs you must use this function.
1297: Level: beginner
1299: .keywords: TS, timestep, set, DAE, Jacobian
1301: .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian()
1302: @*/
1303: PetscErrorCode TSSetIFunction(TS ts,Vec r,TSIFunction f,void *ctx)
1304: {
1306: SNES snes;
1307: Vec ralloc = NULL;
1308: DM dm;
1314: TSGetDM(ts,&dm);
1315: DMTSSetIFunction(dm,f,ctx);
1317: TSGetSNES(ts,&snes);
1318: if (!r && !ts->dm && ts->vec_sol) {
1319: VecDuplicate(ts->vec_sol,&ralloc);
1320: r = ralloc;
1321: }
1322: SNESSetFunction(snes,r,SNESTSFormFunction,ts);
1323: VecDestroy(&ralloc);
1324: return(0);
1325: }
1327: /*@C
1328: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it.
1330: Not Collective
1332: Input Parameter:
1333: . ts - the TS context
1335: Output Parameter:
1336: + r - vector to hold residual (or NULL)
1337: . func - the function to compute residual (or NULL)
1338: - ctx - the function context (or NULL)
1340: Level: advanced
1342: .keywords: TS, nonlinear, get, function
1344: .seealso: TSSetIFunction(), SNESGetFunction()
1345: @*/
1346: PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx)
1347: {
1349: SNES snes;
1350: DM dm;
1354: TSGetSNES(ts,&snes);
1355: SNESGetFunction(snes,r,NULL,NULL);
1356: TSGetDM(ts,&dm);
1357: DMTSGetIFunction(dm,func,ctx);
1358: return(0);
1359: }
1361: /*@C
1362: TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it.
1364: Not Collective
1366: Input Parameter:
1367: . ts - the TS context
1369: Output Parameter:
1370: + r - vector to hold computed right hand side (or NULL)
1371: . func - the function to compute right hand side (or NULL)
1372: - ctx - the function context (or NULL)
1374: Level: advanced
1376: .keywords: TS, nonlinear, get, function
1378: .seealso: TSSetRHSFunction(), SNESGetFunction()
1379: @*/
1380: PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx)
1381: {
1383: SNES snes;
1384: DM dm;
1388: TSGetSNES(ts,&snes);
1389: SNESGetFunction(snes,r,NULL,NULL);
1390: TSGetDM(ts,&dm);
1391: DMTSGetRHSFunction(dm,func,ctx);
1392: return(0);
1393: }
1395: /*@C
1396: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1397: provided with TSSetIFunction().
1399: Logically Collective on TS
1401: Input Parameters:
1402: + ts - the TS context obtained from TSCreate()
1403: . Amat - (approximate) Jacobian matrix
1404: . Pmat - matrix used to compute preconditioner (usually the same as Amat)
1405: . f - the Jacobian evaluation routine
1406: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
1408: Calling sequence of f:
1409: $ f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat Amat,Mat Pmat,void *ctx);
1411: + t - time at step/stage being solved
1412: . U - state vector
1413: . U_t - time derivative of state vector
1414: . a - shift
1415: . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1416: . Pmat - matrix used for constructing preconditioner, usually the same as Amat
1417: - ctx - [optional] user-defined context for matrix evaluation routine
1419: Notes:
1420: The matrices Amat and Pmat are exactly the matrices that are used by SNES for the nonlinear solve.
1422: If you know the operator Amat has a null space you can use MatSetNullSpace() and MatSetTransposeNullSpace() to supply the null
1423: space to Amat and the KSP solvers will automatically use that null space as needed during the solution process.
1425: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1426: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1427: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1428: a and vector W depend on the integration method, step size, and past states. For example with
1429: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1430: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1432: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1434: The TS solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1435: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1437: Level: beginner
1439: .keywords: TS, timestep, DAE, Jacobian
1441: .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESComputeJacobianDefaultColor(), SNESComputeJacobianDefault(), TSSetRHSFunction()
1443: @*/
1444: PetscErrorCode TSSetIJacobian(TS ts,Mat Amat,Mat Pmat,TSIJacobian f,void *ctx)
1445: {
1447: SNES snes;
1448: DM dm;
1457: TSGetDM(ts,&dm);
1458: DMTSSetIJacobian(dm,f,ctx);
1460: TSGetSNES(ts,&snes);
1461: SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1462: return(0);
1463: }
1465: /*@
1466: TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating. Without this flag, TS will change the sign and
1467: shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute
1468: the entire Jacobian. The reuse flag must be set if the evaluation function will assume that the matrix entries have
1469: not been changed by the TS.
1471: Logically Collective
1473: Input Arguments:
1474: + ts - TS context obtained from TSCreate()
1475: - reuse - PETSC_TRUE if the RHS Jacobian
1477: Level: intermediate
1479: .seealso: TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
1480: @*/
1481: PetscErrorCode TSRHSJacobianSetReuse(TS ts,PetscBool reuse)
1482: {
1484: ts->rhsjacobian.reuse = reuse;
1485: return(0);
1486: }
1488: /*@C
1489: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1491: Logically Collective on TS
1493: Input Parameters:
1494: + ts - the TS context obtained from TSCreate()
1495: . F - vector to hold the residual (or NULL to have it created internally)
1496: . fun - the function evaluation routine
1497: - ctx - user-defined context for private data for the function evaluation routine (may be NULL)
1499: Calling sequence of fun:
1500: $ fun(TS ts,PetscReal t,Vec U,Vec U_t,Vec U_tt,Vec F,ctx);
1502: + t - time at step/stage being solved
1503: . U - state vector
1504: . U_t - time derivative of state vector
1505: . U_tt - second time derivative of state vector
1506: . F - function vector
1507: - ctx - [optional] user-defined context for matrix evaluation routine (may be NULL)
1509: Level: beginner
1511: .keywords: TS, timestep, set, ODE, DAE, Function
1513: .seealso: TSSetI2Jacobian()
1514: @*/
1515: PetscErrorCode TSSetI2Function(TS ts,Vec F,TSI2Function fun,void *ctx)
1516: {
1517: DM dm;
1523: TSSetIFunction(ts,F,NULL,NULL);
1524: TSGetDM(ts,&dm);
1525: DMTSSetI2Function(dm,fun,ctx);
1526: return(0);
1527: }
1529: /*@C
1530: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/contex to compute it.
1532: Not Collective
1534: Input Parameter:
1535: . ts - the TS context
1537: Output Parameter:
1538: + r - vector to hold residual (or NULL)
1539: . fun - the function to compute residual (or NULL)
1540: - ctx - the function context (or NULL)
1542: Level: advanced
1544: .keywords: TS, nonlinear, get, function
1546: .seealso: TSSetI2Function(), SNESGetFunction()
1547: @*/
1548: PetscErrorCode TSGetI2Function(TS ts,Vec *r,TSI2Function *fun,void **ctx)
1549: {
1551: SNES snes;
1552: DM dm;
1556: TSGetSNES(ts,&snes);
1557: SNESGetFunction(snes,r,NULL,NULL);
1558: TSGetDM(ts,&dm);
1559: DMTSGetI2Function(dm,fun,ctx);
1560: return(0);
1561: }
1563: /*@C
1564: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1565: where F(t,U,U_t,U_tt) is the function you provided with TSSetI2Function().
1567: Logically Collective on TS
1569: Input Parameters:
1570: + ts - the TS context obtained from TSCreate()
1571: . J - Jacobian matrix
1572: . P - preconditioning matrix for J (may be same as J)
1573: . jac - the Jacobian evaluation routine
1574: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
1576: Calling sequence of jac:
1577: $ jac(TS ts,PetscReal t,Vec U,Vec U_t,Vec U_tt,PetscReal v,PetscReal a,Mat J,Mat P,void *ctx);
1579: + t - time at step/stage being solved
1580: . U - state vector
1581: . U_t - time derivative of state vector
1582: . U_tt - second time derivative of state vector
1583: . v - shift for U_t
1584: . a - shift for U_tt
1585: . J - Jacobian of G(U) = F(t,U,W+v*U,W'+a*U), equivalent to dF/dU + v*dF/dU_t + a*dF/dU_tt
1586: . P - preconditioning matrix for J, may be same as J
1587: - ctx - [optional] user-defined context for matrix evaluation routine
1589: Notes:
1590: The matrices J and P are exactly the matrices that are used by SNES for the nonlinear solve.
1592: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1593: the Jacobian of G(U) = F(t,U,W+v*U,W'+a*U) where F(t,U,U_t,U_tt) = 0 is the DAE to be solved.
1594: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1595: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1597: Level: beginner
1599: .keywords: TS, timestep, set, ODE, DAE, Jacobian
1601: .seealso: TSSetI2Function()
1602: @*/
1603: PetscErrorCode TSSetI2Jacobian(TS ts,Mat J,Mat P,TSI2Jacobian jac,void *ctx)
1604: {
1605: DM dm;
1612: TSSetIJacobian(ts,J,P,NULL,NULL);
1613: TSGetDM(ts,&dm);
1614: DMTSSetI2Jacobian(dm,jac,ctx);
1615: return(0);
1616: }
1618: /*@C
1619: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1621: Not Collective, but parallel objects are returned if TS is parallel
1623: Input Parameter:
1624: . ts - The TS context obtained from TSCreate()
1626: Output Parameters:
1627: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1628: . P - The matrix from which the preconditioner is constructed, often the same as J
1629: . jac - The function to compute the Jacobian matrices
1630: - ctx - User-defined context for Jacobian evaluation routine
1632: Notes: You can pass in NULL for any return argument you do not need.
1634: Level: advanced
1636: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()
1638: .keywords: TS, timestep, get, matrix, Jacobian
1639: @*/
1640: PetscErrorCode TSGetI2Jacobian(TS ts,Mat *J,Mat *P,TSI2Jacobian *jac,void **ctx)
1641: {
1643: SNES snes;
1644: DM dm;
1647: TSGetSNES(ts,&snes);
1648: SNESSetUpMatrices(snes);
1649: SNESGetJacobian(snes,J,P,NULL,NULL);
1650: TSGetDM(ts,&dm);
1651: DMTSGetI2Jacobian(dm,jac,ctx);
1652: return(0);
1653: }
1655: /*@
1656: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1658: Collective on TS and Vec
1660: Input Parameters:
1661: + ts - the TS context
1662: . t - current time
1663: . U - state vector
1664: . V - time derivative of state vector (U_t)
1665: - A - second time derivative of state vector (U_tt)
1667: Output Parameter:
1668: . F - the residual vector
1670: Note:
1671: Most users should not need to explicitly call this routine, as it
1672: is used internally within the nonlinear solvers.
1674: Level: developer
1676: .keywords: TS, compute, function, vector
1678: .seealso: TSSetI2Function()
1679: @*/
1680: PetscErrorCode TSComputeI2Function(TS ts,PetscReal t,Vec U,Vec V,Vec A,Vec F)
1681: {
1682: DM dm;
1683: TSI2Function I2Function;
1684: void *ctx;
1685: TSRHSFunction rhsfunction;
1695: TSGetDM(ts,&dm);
1696: DMTSGetI2Function(dm,&I2Function,&ctx);
1697: DMTSGetRHSFunction(dm,&rhsfunction,NULL);
1699: if (!I2Function) {
1700: TSComputeIFunction(ts,t,U,A,F,PETSC_FALSE);
1701: return(0);
1702: }
1704: PetscLogEventBegin(TS_FunctionEval,ts,U,V,F);
1706: PetscStackPush("TS user implicit function");
1707: I2Function(ts,t,U,V,A,F,ctx);
1708: PetscStackPop;
1710: if (rhsfunction) {
1711: Vec Frhs;
1712: TSGetRHSVec_Private(ts,&Frhs);
1713: TSComputeRHSFunction(ts,t,U,Frhs);
1714: VecAXPY(F,-1,Frhs);
1715: }
1717: PetscLogEventEnd(TS_FunctionEval,ts,U,V,F);
1718: return(0);
1719: }
1721: /*@
1722: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1724: Collective on TS and Vec
1726: Input Parameters:
1727: + ts - the TS context
1728: . t - current timestep
1729: . U - state vector
1730: . V - time derivative of state vector
1731: . A - second time derivative of state vector
1732: . shiftV - shift to apply, see note below
1733: - shiftA - shift to apply, see note below
1735: Output Parameters:
1736: + J - Jacobian matrix
1737: - P - optional preconditioning matrix
1739: Notes:
1740: If F(t,U,V,A)=0 is the DAE, the required Jacobian is
1742: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1744: Most users should not need to explicitly call this routine, as it
1745: is used internally within the nonlinear solvers.
1747: Level: developer
1749: .keywords: TS, compute, Jacobian, matrix
1751: .seealso: TSSetI2Jacobian()
1752: @*/
1753: PetscErrorCode TSComputeI2Jacobian(TS ts,PetscReal t,Vec U,Vec V,Vec A,PetscReal shiftV,PetscReal shiftA,Mat J,Mat P)
1754: {
1755: DM dm;
1756: TSI2Jacobian I2Jacobian;
1757: void *ctx;
1758: TSRHSJacobian rhsjacobian;
1769: TSGetDM(ts,&dm);
1770: DMTSGetI2Jacobian(dm,&I2Jacobian,&ctx);
1771: DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);
1773: if (!I2Jacobian) {
1774: TSComputeIJacobian(ts,t,U,A,shiftA,J,P,PETSC_FALSE);
1775: return(0);
1776: }
1778: PetscLogEventBegin(TS_JacobianEval,ts,U,J,P);
1780: PetscStackPush("TS user implicit Jacobian");
1781: I2Jacobian(ts,t,U,V,A,shiftV,shiftA,J,P,ctx);
1782: PetscStackPop;
1784: if (rhsjacobian) {
1785: Mat Jrhs,Prhs; MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
1786: TSGetRHSMats_Private(ts,&Jrhs,&Prhs);
1787: TSComputeRHSJacobian(ts,t,U,Jrhs,Prhs);
1788: MatAXPY(J,-1,Jrhs,axpy);
1789: if (P != J) {MatAXPY(P,-1,Prhs,axpy);}
1790: }
1792: PetscLogEventEnd(TS_JacobianEval,ts,U,J,P);
1793: return(0);
1794: }
1796: /*@
1797: TS2SetSolution - Sets the initial solution and time derivative vectors
1798: for use by the TS routines handling second order equations.
1800: Logically Collective on TS and Vec
1802: Input Parameters:
1803: + ts - the TS context obtained from TSCreate()
1804: . u - the solution vector
1805: - v - the time derivative vector
1807: Level: beginner
1809: .keywords: TS, timestep, set, solution, initial conditions
1810: @*/
1811: PetscErrorCode TS2SetSolution(TS ts,Vec u,Vec v)
1812: {
1819: TSSetSolution(ts,u);
1820: PetscObjectReference((PetscObject)v);
1821: VecDestroy(&ts->vec_dot);
1822: ts->vec_dot = v;
1823: return(0);
1824: }
1826: /*@
1827: TS2GetSolution - Returns the solution and time derivative at the present timestep
1828: for second order equations. It is valid to call this routine inside the function
1829: that you are evaluating in order to move to the new timestep. This vector not
1830: changed until the solution at the next timestep has been calculated.
1832: Not Collective, but Vec returned is parallel if TS is parallel
1834: Input Parameter:
1835: . ts - the TS context obtained from TSCreate()
1837: Output Parameter:
1838: + u - the vector containing the solution
1839: - v - the vector containing the time derivative
1841: Level: intermediate
1843: .seealso: TS2SetSolution(), TSGetTimeStep(), TSGetTime()
1845: .keywords: TS, timestep, get, solution
1846: @*/
1847: PetscErrorCode TS2GetSolution(TS ts,Vec *u,Vec *v)
1848: {
1853: if (u) *u = ts->vec_sol;
1854: if (v) *v = ts->vec_dot;
1855: return(0);
1856: }
1858: /*@C
1859: TSLoad - Loads a KSP that has been stored in binary with KSPView().
1861: Collective on PetscViewer
1863: Input Parameters:
1864: + newdm - the newly loaded TS, this needs to have been created with TSCreate() or
1865: some related function before a call to TSLoad().
1866: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1868: Level: intermediate
1870: Notes:
1871: The type is determined by the data in the file, any type set into the TS before this call is ignored.
1873: Notes for advanced users:
1874: Most users should not need to know the details of the binary storage
1875: format, since TSLoad() and TSView() completely hide these details.
1876: But for anyone who's interested, the standard binary matrix storage
1877: format is
1878: .vb
1879: has not yet been determined
1880: .ve
1882: .seealso: PetscViewerBinaryOpen(), TSView(), MatLoad(), VecLoad()
1883: @*/
1884: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1885: {
1887: PetscBool isbinary;
1888: PetscInt classid;
1889: char type[256];
1890: DMTS sdm;
1891: DM dm;
1896: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1897: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1899: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1900: if (classid != TS_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Not TS next in file");
1901: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1902: TSSetType(ts, type);
1903: if (ts->ops->load) {
1904: (*ts->ops->load)(ts,viewer);
1905: }
1906: DMCreate(PetscObjectComm((PetscObject)ts),&dm);
1907: DMLoad(dm,viewer);
1908: TSSetDM(ts,dm);
1909: DMCreateGlobalVector(ts->dm,&ts->vec_sol);
1910: VecLoad(ts->vec_sol,viewer);
1911: DMGetDMTS(ts->dm,&sdm);
1912: DMTSLoad(sdm,viewer);
1913: return(0);
1914: }
1916: #include <petscdraw.h>
1917: #if defined(PETSC_HAVE_SAWS)
1918: #include <petscviewersaws.h>
1919: #endif
1920: /*@C
1921: TSView - Prints the TS data structure.
1923: Collective on TS
1925: Input Parameters:
1926: + ts - the TS context obtained from TSCreate()
1927: - viewer - visualization context
1929: Options Database Key:
1930: . -ts_view - calls TSView() at end of TSStep()
1932: Notes:
1933: The available visualization contexts include
1934: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1935: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1936: output where only the first processor opens
1937: the file. All other processors send their
1938: data to the first processor to print.
1940: The user can open an alternative visualization context with
1941: PetscViewerASCIIOpen() - output to a specified file.
1943: Level: beginner
1945: .keywords: TS, timestep, view
1947: .seealso: PetscViewerASCIIOpen()
1948: @*/
1949: PetscErrorCode TSView(TS ts,PetscViewer viewer)
1950: {
1952: TSType type;
1953: PetscBool iascii,isstring,isundials,isbinary,isdraw;
1954: DMTS sdm;
1955: #if defined(PETSC_HAVE_SAWS)
1956: PetscBool issaws;
1957: #endif
1961: if (!viewer) {
1962: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts),&viewer);
1963: }
1967: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1968: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1969: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1970: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1971: #if defined(PETSC_HAVE_SAWS)
1972: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1973: #endif
1974: if (iascii) {
1975: PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer);
1976: if (ts->ops->view) {
1977: PetscViewerASCIIPushTab(viewer);
1978: (*ts->ops->view)(ts,viewer);
1979: PetscViewerASCIIPopTab(viewer);
1980: }
1981: if (ts->max_steps < PETSC_MAX_INT) {
1982: PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);
1983: }
1984: if (ts->max_time < PETSC_MAX_REAL) {
1985: PetscViewerASCIIPrintf(viewer," maximum time=%g\n",(double)ts->max_time);
1986: }
1987: if (ts->usessnes) {
1988: PetscBool lin;
1989: if (ts->problem_type == TS_NONLINEAR) {
1990: PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->snes_its);
1991: }
1992: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->ksp_its);
1993: PetscObjectTypeCompare((PetscObject)ts->snes,SNESKSPONLY,&lin);
1994: PetscViewerASCIIPrintf(viewer," total number of %slinear solve failures=%D\n",lin ? "" : "non",ts->num_snes_failures);
1995: }
1996: PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);
1997: if (ts->vrtol) {
1998: PetscViewerASCIIPrintf(viewer," using vector of relative error tolerances, ");
1999: } else {
2000: PetscViewerASCIIPrintf(viewer," using relative error tolerance of %g, ",(double)ts->rtol);
2001: }
2002: if (ts->vatol) {
2003: PetscViewerASCIIPrintf(viewer," using vector of absolute error tolerances\n");
2004: } else {
2005: PetscViewerASCIIPrintf(viewer," using absolute error tolerance of %g\n",(double)ts->atol);
2006: }
2007: TSAdaptView(ts->adapt,viewer);
2008: if (ts->snes && ts->usessnes) {SNESView(ts->snes,viewer);}
2009: DMGetDMTS(ts->dm,&sdm);
2010: DMTSView(sdm,viewer);
2011: } else if (isstring) {
2012: TSGetType(ts,&type);
2013: PetscViewerStringSPrintf(viewer," %-7.7s",type);
2014: } else if (isbinary) {
2015: PetscInt classid = TS_FILE_CLASSID;
2016: MPI_Comm comm;
2017: PetscMPIInt rank;
2018: char type[256];
2020: PetscObjectGetComm((PetscObject)ts,&comm);
2021: MPI_Comm_rank(comm,&rank);
2022: if (!rank) {
2023: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
2024: PetscStrncpy(type,((PetscObject)ts)->type_name,256);
2025: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
2026: }
2027: if (ts->ops->view) {
2028: (*ts->ops->view)(ts,viewer);
2029: }
2030: if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
2031: DMView(ts->dm,viewer);
2032: VecView(ts->vec_sol,viewer);
2033: DMGetDMTS(ts->dm,&sdm);
2034: DMTSView(sdm,viewer);
2035: } else if (isdraw) {
2036: PetscDraw draw;
2037: char str[36];
2038: PetscReal x,y,bottom,h;
2040: PetscViewerDrawGetDraw(viewer,0,&draw);
2041: PetscDrawGetCurrentPoint(draw,&x,&y);
2042: PetscStrcpy(str,"TS: ");
2043: PetscStrcat(str,((PetscObject)ts)->type_name);
2044: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,NULL,&h);
2045: bottom = y - h;
2046: PetscDrawPushCurrentPoint(draw,x,bottom);
2047: if (ts->ops->view) {
2048: (*ts->ops->view)(ts,viewer);
2049: }
2050: if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
2051: if (ts->snes) {SNESView(ts->snes,viewer);}
2052: PetscDrawPopCurrentPoint(draw);
2053: #if defined(PETSC_HAVE_SAWS)
2054: } else if (issaws) {
2055: PetscMPIInt rank;
2056: const char *name;
2058: PetscObjectGetName((PetscObject)ts,&name);
2059: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
2060: if (!((PetscObject)ts)->amsmem && !rank) {
2061: char dir[1024];
2063: PetscObjectViewSAWs((PetscObject)ts,viewer);
2064: PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time_step",name);
2065: PetscStackCallSAWs(SAWs_Register,(dir,&ts->steps,1,SAWs_READ,SAWs_INT));
2066: PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time",name);
2067: PetscStackCallSAWs(SAWs_Register,(dir,&ts->ptime,1,SAWs_READ,SAWs_DOUBLE));
2068: }
2069: if (ts->ops->view) {
2070: (*ts->ops->view)(ts,viewer);
2071: }
2072: #endif
2073: }
2075: PetscViewerASCIIPushTab(viewer);
2076: PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);
2077: PetscViewerASCIIPopTab(viewer);
2078: return(0);
2079: }
2081: /*@
2082: TSSetApplicationContext - Sets an optional user-defined context for
2083: the timesteppers.
2085: Logically Collective on TS
2087: Input Parameters:
2088: + ts - the TS context obtained from TSCreate()
2089: - usrP - optional user context
2091: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
2092: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
2094: Level: intermediate
2096: .keywords: TS, timestep, set, application, context
2098: .seealso: TSGetApplicationContext()
2099: @*/
2100: PetscErrorCode TSSetApplicationContext(TS ts,void *usrP)
2101: {
2104: ts->user = usrP;
2105: return(0);
2106: }
2108: /*@
2109: TSGetApplicationContext - Gets the user-defined context for the
2110: timestepper.
2112: Not Collective
2114: Input Parameter:
2115: . ts - the TS context obtained from TSCreate()
2117: Output Parameter:
2118: . usrP - user context
2120: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
2121: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
2123: Level: intermediate
2125: .keywords: TS, timestep, get, application, context
2127: .seealso: TSSetApplicationContext()
2128: @*/
2129: PetscErrorCode TSGetApplicationContext(TS ts,void *usrP)
2130: {
2133: *(void**)usrP = ts->user;
2134: return(0);
2135: }
2137: /*@
2138: TSGetStepNumber - Gets the number of steps completed.
2140: Not Collective
2142: Input Parameter:
2143: . ts - the TS context obtained from TSCreate()
2145: Output Parameter:
2146: . steps - number of steps completed so far
2148: Level: intermediate
2150: .keywords: TS, timestep, get, iteration, number
2151: .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSSetPostStep()
2152: @*/
2153: PetscErrorCode TSGetStepNumber(TS ts,PetscInt *steps)
2154: {
2158: *steps = ts->steps;
2159: return(0);
2160: }
2162: /*@
2163: TSSetStepNumber - Sets the number of steps completed.
2165: Logically Collective on TS
2167: Input Parameters:
2168: + ts - the TS context
2169: - steps - number of steps completed so far
2171: Notes:
2172: For most uses of the TS solvers the user need not explicitly call
2173: TSSetStepNumber(), as the step counter is appropriately updated in
2174: TSSolve()/TSStep()/TSRollBack(). Power users may call this routine to
2175: reinitialize timestepping by setting the step counter to zero (and time
2176: to the initial time) to solve a similar problem with different initial
2177: conditions or parameters. Other possible use case is to continue
2178: timestepping from a previously interrupted run in such a way that TS
2179: monitors will be called with a initial nonzero step counter.
2181: Level: advanced
2183: .keywords: TS, timestep, set, iteration, number
2184: .seealso: TSGetStepNumber(), TSSetTime(), TSSetTimeStep(), TSSetSolution()
2185: @*/
2186: PetscErrorCode TSSetStepNumber(TS ts,PetscInt steps)
2187: {
2191: if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Step number must be non-negative");
2192: ts->steps = steps;
2193: return(0);
2194: }
2196: /*@
2197: TSSetTimeStep - Allows one to reset the timestep at any time,
2198: useful for simple pseudo-timestepping codes.
2200: Logically Collective on TS
2202: Input Parameters:
2203: + ts - the TS context obtained from TSCreate()
2204: - time_step - the size of the timestep
2206: Level: intermediate
2208: .seealso: TSGetTimeStep(), TSSetTime()
2210: .keywords: TS, set, timestep
2211: @*/
2212: PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step)
2213: {
2217: ts->time_step = time_step;
2218: return(0);
2219: }
2221: /*@
2222: TSSetExactFinalTime - Determines whether to adapt the final time step to
2223: match the exact final time, interpolate solution to the exact final time,
2224: or just return at the final time TS computed.
2226: Logically Collective on TS
2228: Input Parameter:
2229: + ts - the time-step context
2230: - eftopt - exact final time option
2232: $ TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded
2233: $ TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2234: $ TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2236: Options Database:
2237: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime
2239: Warning: If you use the option TS_EXACTFINALTIME_STEPOVER the solution may be at a very different time
2240: then the final time you selected.
2242: Level: beginner
2244: .seealso: TSExactFinalTimeOption, TSGetExactFinalTime()
2245: @*/
2246: PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt)
2247: {
2251: ts->exact_final_time = eftopt;
2252: return(0);
2253: }
2255: /*@
2256: TSGetExactFinalTime - Gets the exact final time option.
2258: Not Collective
2260: Input Parameter:
2261: . ts - the TS context
2263: Output Parameter:
2264: . eftopt - exact final time option
2266: Level: beginner
2268: .seealso: TSExactFinalTimeOption, TSSetExactFinalTime()
2269: @*/
2270: PetscErrorCode TSGetExactFinalTime(TS ts,TSExactFinalTimeOption *eftopt)
2271: {
2275: *eftopt = ts->exact_final_time;
2276: return(0);
2277: }
2279: /*@
2280: TSGetTimeStep - Gets the current timestep size.
2282: Not Collective
2284: Input Parameter:
2285: . ts - the TS context obtained from TSCreate()
2287: Output Parameter:
2288: . dt - the current timestep size
2290: Level: intermediate
2292: .seealso: TSSetTimeStep(), TSGetTime()
2294: .keywords: TS, get, timestep
2295: @*/
2296: PetscErrorCode TSGetTimeStep(TS ts,PetscReal *dt)
2297: {
2301: *dt = ts->time_step;
2302: return(0);
2303: }
2305: /*@
2306: TSGetSolution - Returns the solution at the present timestep. It
2307: is valid to call this routine inside the function that you are evaluating
2308: in order to move to the new timestep. This vector not changed until
2309: the solution at the next timestep has been calculated.
2311: Not Collective, but Vec returned is parallel if TS is parallel
2313: Input Parameter:
2314: . ts - the TS context obtained from TSCreate()
2316: Output Parameter:
2317: . v - the vector containing the solution
2319: Note: If you used TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); this does not return the solution at the requested
2320: final time. It returns the solution at the next timestep.
2322: Level: intermediate
2324: .seealso: TSGetTimeStep(), TSGetTime(), TSGetSolveTime(), TSGetSolutionComponents()
2326: .keywords: TS, timestep, get, solution
2327: @*/
2328: PetscErrorCode TSGetSolution(TS ts,Vec *v)
2329: {
2333: *v = ts->vec_sol;
2334: return(0);
2335: }
2337: /*@
2338: TSGetSolutionComponents - Returns any solution components at the present
2339: timestep, if available for the time integration method being used.
2340: Solution components are quantities that share the same size and
2341: structure as the solution vector.
2343: Not Collective, but Vec returned is parallel if TS is parallel
2345: Parameters :
2346: . ts - the TS context obtained from TSCreate() (input parameter).
2347: . n - If v is PETSC_NULL, then the number of solution components is
2348: returned through n, else the n-th solution component is
2349: returned in v.
2350: . v - the vector containing the n-th solution component
2351: (may be PETSC_NULL to use this function to find out
2352: the number of solutions components).
2354: Level: advanced
2356: .seealso: TSGetSolution()
2358: .keywords: TS, timestep, get, solution
2359: @*/
2360: PetscErrorCode TSGetSolutionComponents(TS ts,PetscInt *n,Vec *v)
2361: {
2366: if (!ts->ops->getsolutioncomponents) *n = 0;
2367: else {
2368: (*ts->ops->getsolutioncomponents)(ts,n,v);
2369: }
2370: return(0);
2371: }
2373: /*@
2374: TSGetAuxSolution - Returns an auxiliary solution at the present
2375: timestep, if available for the time integration method being used.
2377: Not Collective, but Vec returned is parallel if TS is parallel
2379: Parameters :
2380: . ts - the TS context obtained from TSCreate() (input parameter).
2381: . v - the vector containing the auxiliary solution
2383: Level: intermediate
2385: .seealso: TSGetSolution()
2387: .keywords: TS, timestep, get, solution
2388: @*/
2389: PetscErrorCode TSGetAuxSolution(TS ts,Vec *v)
2390: {
2395: if (ts->ops->getauxsolution) {
2396: (*ts->ops->getauxsolution)(ts,v);
2397: } else {
2398: VecZeroEntries(*v);
2399: }
2400: return(0);
2401: }
2403: /*@
2404: TSGetTimeError - Returns the estimated error vector, if the chosen
2405: TSType has an error estimation functionality.
2407: Not Collective, but Vec returned is parallel if TS is parallel
2409: Note: MUST call after TSSetUp()
2411: Parameters :
2412: . ts - the TS context obtained from TSCreate() (input parameter).
2413: . n - current estimate (n=0) or previous one (n=-1)
2414: . v - the vector containing the error (same size as the solution).
2416: Level: intermediate
2418: .seealso: TSGetSolution(), TSSetTimeError()
2420: .keywords: TS, timestep, get, error
2421: @*/
2422: PetscErrorCode TSGetTimeError(TS ts,PetscInt n,Vec *v)
2423: {
2428: if (ts->ops->gettimeerror) {
2429: (*ts->ops->gettimeerror)(ts,n,v);
2430: } else {
2431: VecZeroEntries(*v);
2432: }
2433: return(0);
2434: }
2436: /*@
2437: TSSetTimeError - Sets the estimated error vector, if the chosen
2438: TSType has an error estimation functionality. This can be used
2439: to restart such a time integrator with a given error vector.
2441: Not Collective, but Vec returned is parallel if TS is parallel
2443: Parameters :
2444: . ts - the TS context obtained from TSCreate() (input parameter).
2445: . v - the vector containing the error (same size as the solution).
2447: Level: intermediate
2449: .seealso: TSSetSolution(), TSGetTimeError)
2451: .keywords: TS, timestep, get, error
2452: @*/
2453: PetscErrorCode TSSetTimeError(TS ts,Vec v)
2454: {
2459: if (!ts->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetUp() first");
2460: if (ts->ops->settimeerror) {
2461: (*ts->ops->settimeerror)(ts,v);
2462: }
2463: return(0);
2464: }
2466: /*@
2467: TSGetCostGradients - Returns the gradients from the TSAdjointSolve()
2469: Not Collective, but Vec returned is parallel if TS is parallel
2471: Input Parameter:
2472: . ts - the TS context obtained from TSCreate()
2474: Output Parameter:
2475: + lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
2476: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
2478: Level: intermediate
2480: .seealso: TSGetTimeStep()
2482: .keywords: TS, timestep, get, sensitivity
2483: @*/
2484: PetscErrorCode TSGetCostGradients(TS ts,PetscInt *numcost,Vec **lambda,Vec **mu)
2485: {
2488: if (numcost) *numcost = ts->numcost;
2489: if (lambda) *lambda = ts->vecs_sensi;
2490: if (mu) *mu = ts->vecs_sensip;
2491: return(0);
2492: }
2494: /* ----- Routines to initialize and destroy a timestepper ---- */
2495: /*@
2496: TSSetProblemType - Sets the type of problem to be solved.
2498: Not collective
2500: Input Parameters:
2501: + ts - The TS
2502: - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2503: .vb
2504: U_t - A U = 0 (linear)
2505: U_t - A(t) U = 0 (linear)
2506: F(t,U,U_t) = 0 (nonlinear)
2507: .ve
2509: Level: beginner
2511: .keywords: TS, problem type
2512: .seealso: TSSetUp(), TSProblemType, TS
2513: @*/
2514: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2515: {
2520: ts->problem_type = type;
2521: if (type == TS_LINEAR) {
2522: SNES snes;
2523: TSGetSNES(ts,&snes);
2524: SNESSetType(snes,SNESKSPONLY);
2525: }
2526: return(0);
2527: }
2529: /*@C
2530: TSGetProblemType - Gets the type of problem to be solved.
2532: Not collective
2534: Input Parameter:
2535: . ts - The TS
2537: Output Parameter:
2538: . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2539: .vb
2540: M U_t = A U
2541: M(t) U_t = A(t) U
2542: F(t,U,U_t)
2543: .ve
2545: Level: beginner
2547: .keywords: TS, problem type
2548: .seealso: TSSetUp(), TSProblemType, TS
2549: @*/
2550: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2551: {
2555: *type = ts->problem_type;
2556: return(0);
2557: }
2559: /*@
2560: TSSetUp - Sets up the internal data structures for the later use
2561: of a timestepper.
2563: Collective on TS
2565: Input Parameter:
2566: . ts - the TS context obtained from TSCreate()
2568: Notes:
2569: For basic use of the TS solvers the user need not explicitly call
2570: TSSetUp(), since these actions will automatically occur during
2571: the call to TSStep() or TSSolve(). However, if one wishes to control this
2572: phase separately, TSSetUp() should be called after TSCreate()
2573: and optional routines of the form TSSetXXX(), but before TSStep() and TSSolve().
2575: Level: advanced
2577: .keywords: TS, timestep, setup
2579: .seealso: TSCreate(), TSStep(), TSDestroy(), TSSolve()
2580: @*/
2581: PetscErrorCode TSSetUp(TS ts)
2582: {
2584: DM dm;
2585: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2586: PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*);
2587: TSIFunction ifun;
2588: TSIJacobian ijac;
2589: TSI2Jacobian i2jac;
2590: TSRHSJacobian rhsjac;
2591: PetscBool isnone;
2595: if (ts->setupcalled) return(0);
2597: if (!((PetscObject)ts)->type_name) {
2598: TSGetIFunction(ts,NULL,&ifun,NULL);
2599: TSSetType(ts,ifun ? TSBEULER : TSEULER);
2600: }
2602: if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first");
2604: if (ts->rhsjacobian.reuse) {
2605: Mat Amat,Pmat;
2606: SNES snes;
2607: TSGetSNES(ts,&snes);
2608: SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);
2609: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2610: * have displaced the RHS matrix */
2611: if (Amat == ts->Arhs) {
2612: /* we need to copy the values of the matrix because for the constant Jacobian case the user will never set the numerical values in this new location */
2613: MatDuplicate(ts->Arhs,MAT_COPY_VALUES,&Amat);
2614: SNESSetJacobian(snes,Amat,NULL,NULL,NULL);
2615: MatDestroy(&Amat);
2616: }
2617: if (Pmat == ts->Brhs) {
2618: MatDuplicate(ts->Brhs,MAT_COPY_VALUES,&Pmat);
2619: SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);
2620: MatDestroy(&Pmat);
2621: }
2622: }
2624: TSGetAdapt(ts,&ts->adapt);
2625: TSAdaptSetDefaultType(ts->adapt,ts->default_adapt_type);
2627: if (ts->ops->setup) {
2628: (*ts->ops->setup)(ts);
2629: }
2631: /* Attempt to check/preset a default value for the exact final time option */
2632: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&isnone);
2633: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
2634: ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2636: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2637: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2638: */
2639: TSGetDM(ts,&dm);
2640: DMSNESGetFunction(dm,&func,NULL);
2641: if (!func) {
2642: DMSNESSetFunction(dm,SNESTSFormFunction,ts);
2643: }
2644: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2645: Otherwise, the SNES will use coloring internally to form the Jacobian.
2646: */
2647: DMSNESGetJacobian(dm,&jac,NULL);
2648: DMTSGetIJacobian(dm,&ijac,NULL);
2649: DMTSGetI2Jacobian(dm,&i2jac,NULL);
2650: DMTSGetRHSJacobian(dm,&rhsjac,NULL);
2651: if (!jac && (ijac || i2jac || rhsjac)) {
2652: DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);
2653: }
2655: /* if time integration scheme has a starting method, call it */
2656: if (ts->ops->startingmethod) {
2657: (*ts->ops->startingmethod)(ts);
2658: }
2660: ts->setupcalled = PETSC_TRUE;
2661: return(0);
2662: }
2664: /*@
2665: TSAdjointSetUp - Sets up the internal data structures for the later use
2666: of an adjoint solver
2668: Collective on TS
2670: Input Parameter:
2671: . ts - the TS context obtained from TSCreate()
2673: Level: advanced
2675: .keywords: TS, timestep, setup
2677: .seealso: TSCreate(), TSAdjointStep(), TSSetCostGradients()
2678: @*/
2679: PetscErrorCode TSAdjointSetUp(TS ts)
2680: {
2685: if (ts->adjointsetupcalled) return(0);
2686: if (!ts->vecs_sensi) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetCostGradients() first");
2687: if (ts->vecs_sensip && !ts->Jacp) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSAdjointSetRHSJacobian() first");
2689: if (ts->vec_costintegral) { /* if there is integral in the cost function */
2690: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&ts->vecs_drdy);
2691: if (ts->vecs_sensip){
2692: VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&ts->vecs_drdp);
2693: }
2694: }
2696: if (ts->ops->adjointsetup) {
2697: (*ts->ops->adjointsetup)(ts);
2698: }
2699: ts->adjointsetupcalled = PETSC_TRUE;
2700: return(0);
2701: }
2703: /*@
2704: TSReset - Resets a TS context and removes any allocated Vecs and Mats.
2706: Collective on TS
2708: Input Parameter:
2709: . ts - the TS context obtained from TSCreate()
2711: Level: beginner
2713: .keywords: TS, timestep, reset
2715: .seealso: TSCreate(), TSSetup(), TSDestroy()
2716: @*/
2717: PetscErrorCode TSReset(TS ts)
2718: {
2724: if (ts->ops->reset) {
2725: (*ts->ops->reset)(ts);
2726: }
2727: if (ts->snes) {SNESReset(ts->snes);}
2728: if (ts->adapt) {TSAdaptReset(ts->adapt);}
2730: MatDestroy(&ts->Arhs);
2731: MatDestroy(&ts->Brhs);
2732: VecDestroy(&ts->Frhs);
2733: VecDestroy(&ts->vec_sol);
2734: VecDestroy(&ts->vec_dot);
2735: VecDestroy(&ts->vatol);
2736: VecDestroy(&ts->vrtol);
2737: VecDestroyVecs(ts->nwork,&ts->work);
2739: VecDestroyVecs(ts->numcost,&ts->vecs_drdy);
2740: VecDestroyVecs(ts->numcost,&ts->vecs_drdp);
2742: MatDestroy(&ts->Jacp);
2743: VecDestroy(&ts->vec_costintegral);
2744: VecDestroy(&ts->vec_costintegrand);
2746: PetscFree(ts->vecs_fwdsensipacked);
2748: ts->setupcalled = PETSC_FALSE;
2749: return(0);
2750: }
2752: /*@
2753: TSDestroy - Destroys the timestepper context that was created
2754: with TSCreate().
2756: Collective on TS
2758: Input Parameter:
2759: . ts - the TS context obtained from TSCreate()
2761: Level: beginner
2763: .keywords: TS, timestepper, destroy
2765: .seealso: TSCreate(), TSSetUp(), TSSolve()
2766: @*/
2767: PetscErrorCode TSDestroy(TS *ts)
2768: {
2772: if (!*ts) return(0);
2774: if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; return(0);}
2776: TSReset((*ts));
2778: /* if memory was published with SAWs then destroy it */
2779: PetscObjectSAWsViewOff((PetscObject)*ts);
2780: if ((*ts)->ops->destroy) {(*(*ts)->ops->destroy)((*ts));}
2782: TSTrajectoryDestroy(&(*ts)->trajectory);
2784: TSAdaptDestroy(&(*ts)->adapt);
2785: TSEventDestroy(&(*ts)->event);
2787: SNESDestroy(&(*ts)->snes);
2788: DMDestroy(&(*ts)->dm);
2789: TSMonitorCancel((*ts));
2790: TSAdjointMonitorCancel((*ts));
2792: PetscHeaderDestroy(ts);
2793: return(0);
2794: }
2796: /*@
2797: TSGetSNES - Returns the SNES (nonlinear solver) associated with
2798: a TS (timestepper) context. Valid only for nonlinear problems.
2800: Not Collective, but SNES is parallel if TS is parallel
2802: Input Parameter:
2803: . ts - the TS context obtained from TSCreate()
2805: Output Parameter:
2806: . snes - the nonlinear solver context
2808: Notes:
2809: The user can then directly manipulate the SNES context to set various
2810: options, etc. Likewise, the user can then extract and manipulate the
2811: KSP, KSP, and PC contexts as well.
2813: TSGetSNES() does not work for integrators that do not use SNES; in
2814: this case TSGetSNES() returns NULL in snes.
2816: Level: beginner
2818: .keywords: timestep, get, SNES
2819: @*/
2820: PetscErrorCode TSGetSNES(TS ts,SNES *snes)
2821: {
2827: if (!ts->snes) {
2828: SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);
2829: SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
2830: PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);
2831: PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);
2832: if (ts->dm) {SNESSetDM(ts->snes,ts->dm);}
2833: if (ts->problem_type == TS_LINEAR) {
2834: SNESSetType(ts->snes,SNESKSPONLY);
2835: }
2836: }
2837: *snes = ts->snes;
2838: return(0);
2839: }
2841: /*@
2842: TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context
2844: Collective
2846: Input Parameter:
2847: + ts - the TS context obtained from TSCreate()
2848: - snes - the nonlinear solver context
2850: Notes:
2851: Most users should have the TS created by calling TSGetSNES()
2853: Level: developer
2855: .keywords: timestep, set, SNES
2856: @*/
2857: PetscErrorCode TSSetSNES(TS ts,SNES snes)
2858: {
2860: PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*);
2865: PetscObjectReference((PetscObject)snes);
2866: SNESDestroy(&ts->snes);
2868: ts->snes = snes;
2870: SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
2871: SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);
2872: if (func == SNESTSFormJacobian) {
2873: SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);
2874: }
2875: return(0);
2876: }
2878: /*@
2879: TSGetKSP - Returns the KSP (linear solver) associated with
2880: a TS (timestepper) context.
2882: Not Collective, but KSP is parallel if TS is parallel
2884: Input Parameter:
2885: . ts - the TS context obtained from TSCreate()
2887: Output Parameter:
2888: . ksp - the nonlinear solver context
2890: Notes:
2891: The user can then directly manipulate the KSP context to set various
2892: options, etc. Likewise, the user can then extract and manipulate the
2893: KSP and PC contexts as well.
2895: TSGetKSP() does not work for integrators that do not use KSP;
2896: in this case TSGetKSP() returns NULL in ksp.
2898: Level: beginner
2900: .keywords: timestep, get, KSP
2901: @*/
2902: PetscErrorCode TSGetKSP(TS ts,KSP *ksp)
2903: {
2905: SNES snes;
2910: if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first");
2911: if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()");
2912: TSGetSNES(ts,&snes);
2913: SNESGetKSP(snes,ksp);
2914: return(0);
2915: }
2917: /* ----------- Routines to set solver parameters ---------- */
2919: /*@
2920: TSSetMaxSteps - Sets the maximum number of steps to use.
2922: Logically Collective on TS
2924: Input Parameters:
2925: + ts - the TS context obtained from TSCreate()
2926: - maxsteps - maximum number of steps to use
2928: Options Database Keys:
2929: . -ts_max_steps <maxsteps> - Sets maxsteps
2931: Notes:
2932: The default maximum number of steps is 5000
2934: Level: intermediate
2936: .keywords: TS, timestep, set, maximum, steps
2938: .seealso: TSGetMaxSteps(), TSSetMaxTime(), TSSetExactFinalTime()
2939: @*/
2940: PetscErrorCode TSSetMaxSteps(TS ts,PetscInt maxsteps)
2941: {
2945: if (maxsteps < 0 ) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of steps must be non-negative");
2946: ts->max_steps = maxsteps;
2947: return(0);
2948: }
2950: /*@
2951: TSGetMaxSteps - Gets the maximum number of steps to use.
2953: Not Collective
2955: Input Parameters:
2956: . ts - the TS context obtained from TSCreate()
2958: Output Parameter:
2959: . maxsteps - maximum number of steps to use
2961: Level: advanced
2963: .keywords: TS, timestep, get, maximum, steps
2965: .seealso: TSSetMaxSteps(), TSGetMaxTime(), TSSetMaxTime()
2966: @*/
2967: PetscErrorCode TSGetMaxSteps(TS ts,PetscInt *maxsteps)
2968: {
2972: *maxsteps = ts->max_steps;
2973: return(0);
2974: }
2976: /*@
2977: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
2979: Logically Collective on TS
2981: Input Parameters:
2982: + ts - the TS context obtained from TSCreate()
2983: - maxtime - final time to step to
2985: Options Database Keys:
2986: . -ts_max_time <maxtime> - Sets maxtime
2988: Notes:
2989: The default maximum time is 5.0
2991: Level: intermediate
2993: .keywords: TS, timestep, set, maximum, time
2995: .seealso: TSGetMaxTime(), TSSetMaxSteps(), TSSetExactFinalTime()
2996: @*/
2997: PetscErrorCode TSSetMaxTime(TS ts,PetscReal maxtime)
2998: {
3002: ts->max_time = maxtime;
3003: return(0);
3004: }
3006: /*@
3007: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
3009: Not Collective
3011: Input Parameters:
3012: . ts - the TS context obtained from TSCreate()
3014: Output Parameter:
3015: . maxtime - final time to step to
3017: Level: advanced
3019: .keywords: TS, timestep, get, maximum, time
3021: .seealso: TSSetMaxTime(), TSGetMaxSteps(), TSSetMaxSteps()
3022: @*/
3023: PetscErrorCode TSGetMaxTime(TS ts,PetscReal *maxtime)
3024: {
3028: *maxtime = ts->max_time;
3029: return(0);
3030: }
3032: /*@
3033: TSSetInitialTimeStep - Deprecated, use TSSetTime() and TSSetTimeStep().
3035: Level: deprecated
3037: @*/
3038: PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step)
3039: {
3043: TSSetTime(ts,initial_time);
3044: TSSetTimeStep(ts,time_step);
3045: return(0);
3046: }
3048: /*@
3049: TSGetDuration - Deprecated, use TSGetMaxSteps() and TSGetMaxTime().
3051: Level: deprecated
3053: @*/
3054: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
3055: {
3058: if (maxsteps) {
3060: *maxsteps = ts->max_steps;
3061: }
3062: if (maxtime) {
3064: *maxtime = ts->max_time;
3065: }
3066: return(0);
3067: }
3069: /*@
3070: TSSetDuration - Deprecated, use TSSetMaxSteps() and TSSetMaxTime().
3072: Level: deprecated
3074: @*/
3075: PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime)
3076: {
3081: if (maxsteps >= 0) ts->max_steps = maxsteps;
3082: if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
3083: return(0);
3084: }
3086: /*@
3087: TSGetTimeStepNumber - Deprecated, use TSGetStepNumber().
3089: Level: deprecated
3091: @*/
3092: PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *steps) { return TSGetStepNumber(ts,steps); }
3094: /*@
3095: TSGetTotalSteps - Deprecated, use TSGetStepNumber().
3097: Level: deprecated
3099: @*/
3100: PetscErrorCode TSGetTotalSteps(TS ts,PetscInt *steps) { return TSGetStepNumber(ts,steps); }
3102: /*@
3103: TSSetSolution - Sets the initial solution vector
3104: for use by the TS routines.
3106: Logically Collective on TS and Vec
3108: Input Parameters:
3109: + ts - the TS context obtained from TSCreate()
3110: - u - the solution vector
3112: Level: beginner
3114: .keywords: TS, timestep, set, solution, initial values
3115: @*/
3116: PetscErrorCode TSSetSolution(TS ts,Vec u)
3117: {
3119: DM dm;
3124: PetscObjectReference((PetscObject)u);
3125: VecDestroy(&ts->vec_sol);
3126: ts->vec_sol = u;
3128: TSGetDM(ts,&dm);
3129: DMShellSetGlobalVector(dm,u);
3130: return(0);
3131: }
3133: /*@
3134: TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time
3136: Logically Collective on TS
3138: Input Parameters:
3139: + ts - the TS context obtained from TSCreate()
3140: . steps - number of steps to use
3142: Level: intermediate
3144: Notes: Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this
3145: so as to integrate back to less than the original timestep
3147: .keywords: TS, timestep, set, maximum, iterations
3149: .seealso: TSSetExactFinalTime()
3150: @*/
3151: PetscErrorCode TSAdjointSetSteps(TS ts,PetscInt steps)
3152: {
3156: if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back a negative number of steps");
3157: if (steps > ts->steps) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back more than the total number of forward steps");
3158: ts->adjoint_max_steps = steps;
3159: return(0);
3160: }
3162: /*@
3163: TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial values and w.r.t. the problem parameters
3164: for use by the TSAdjoint routines.
3166: Logically Collective on TS and Vec
3168: Input Parameters:
3169: + ts - the TS context obtained from TSCreate()
3170: . lambda - gradients with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
3171: - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
3173: Level: beginner
3175: Notes: the entries in these vectors must be correctly initialized with the values lamda_i = df/dy|finaltime mu_i = df/dp|finaltime
3177: .keywords: TS, timestep, set, sensitivity, initial values
3178: @*/
3179: PetscErrorCode TSSetCostGradients(TS ts,PetscInt numcost,Vec *lambda,Vec *mu)
3180: {
3184: ts->vecs_sensi = lambda;
3185: ts->vecs_sensip = mu;
3186: if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2rd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand");
3187: ts->numcost = numcost;
3188: return(0);
3189: }
3191: /*@C
3192: TSAdjointSetRHSJacobian - Sets the function that computes the Jacobian of G w.r.t. the parameters p where y_t = G(y,p,t), as well as the location to store the matrix.
3194: Logically Collective on TS
3196: Input Parameters:
3197: + ts - The TS context obtained from TSCreate()
3198: - func - The function
3200: Calling sequence of func:
3201: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
3202: + t - current timestep
3203: . y - input vector (current ODE solution)
3204: . A - output matrix
3205: - ctx - [optional] user-defined function context
3207: Level: intermediate
3209: Notes: Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p
3211: .keywords: TS, sensitivity
3212: .seealso:
3213: @*/
3214: PetscErrorCode TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
3215: {
3222: ts->rhsjacobianp = func;
3223: ts->rhsjacobianpctx = ctx;
3224: if(Amat) {
3225: PetscObjectReference((PetscObject)Amat);
3226: MatDestroy(&ts->Jacp);
3227: ts->Jacp = Amat;
3228: }
3229: return(0);
3230: }
3232: /*@C
3233: TSAdjointComputeRHSJacobian - Runs the user-defined Jacobian function.
3235: Collective on TS
3237: Input Parameters:
3238: . ts - The TS context obtained from TSCreate()
3240: Level: developer
3242: .keywords: TS, sensitivity
3243: .seealso: TSAdjointSetRHSJacobian()
3244: @*/
3245: PetscErrorCode TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat)
3246: {
3254: PetscStackPush("TS user JacobianP function for sensitivity analysis");
3255: (*ts->rhsjacobianp)(ts,t,X,Amat,ts->rhsjacobianpctx);
3256: PetscStackPop;
3257: return(0);
3258: }
3260: /*@C
3261: TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions
3263: Logically Collective on TS
3265: Input Parameters:
3266: + ts - the TS context obtained from TSCreate()
3267: . numcost - number of gradients to be computed, this is the number of cost functions
3268: . costintegral - vector that stores the integral values
3269: . rf - routine for evaluating the integrand function
3270: . drdyf - function that computes the gradients of the r's with respect to y,NULL if not a function y
3271: . drdpf - function that computes the gradients of the r's with respect to p, NULL if not a function of p
3272: . fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
3273: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
3275: Calling sequence of rf:
3276: $ PetscErrorCode rf(TS ts,PetscReal t,Vec y,Vec f,void *ctx);
3278: Calling sequence of drdyf:
3279: $ PetscErroCode drdyf(TS ts,PetscReal t,Vec y,Vec *drdy,void *ctx);
3281: Calling sequence of drdpf:
3282: $ PetscErroCode drdpf(TS ts,PetscReal t,Vec y,Vec *drdp,void *ctx);
3284: Level: intermediate
3286: Notes: For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions
3288: .keywords: TS, sensitivity analysis, timestep, set, quadrature, function
3290: .seealso: TSAdjointSetRHSJacobian(),TSGetCostGradients(), TSSetCostGradients()
3291: @*/
3292: PetscErrorCode TSSetCostIntegrand(TS ts,PetscInt numcost,Vec costintegral,PetscErrorCode (*rf)(TS,PetscReal,Vec,Vec,void*),
3293: PetscErrorCode (*drdyf)(TS,PetscReal,Vec,Vec*,void*),
3294: PetscErrorCode (*drdpf)(TS,PetscReal,Vec,Vec*,void*),
3295: PetscBool fwd,void *ctx)
3296: {
3302: if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2rd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostGradients() or TSForwardSetIntegralGradients()");
3303: if (!ts->numcost) ts->numcost=numcost;
3305: if (costintegral) {
3306: PetscObjectReference((PetscObject)costintegral);
3307: VecDestroy(&ts->vec_costintegral);
3308: ts->vec_costintegral = costintegral;
3309: } else {
3310: if (!ts->vec_costintegral) { /* Create a seq vec if user does not provide one */
3311: VecCreateSeq(PETSC_COMM_SELF,numcost,&ts->vec_costintegral);
3312: } else {
3313: VecSet(ts->vec_costintegral,0.0);
3314: }
3315: }
3316: if (!ts->vec_costintegrand) {
3317: VecDuplicate(ts->vec_costintegral,&ts->vec_costintegrand);
3318: } else {
3319: VecSet(ts->vec_costintegrand,0.0);
3320: }
3321: ts->costintegralfwd = fwd; /* Evaluate the cost integral in forward run if fwd is true */
3322: ts->costintegrand = rf;
3323: ts->costintegrandctx = ctx;
3324: ts->drdyfunction = drdyf;
3325: ts->drdpfunction = drdpf;
3326: return(0);
3327: }
3329: /*@
3330: TSGetCostIntegral - Returns the values of the integral term in the cost functions.
3331: It is valid to call the routine after a backward run.
3333: Not Collective
3335: Input Parameter:
3336: . ts - the TS context obtained from TSCreate()
3338: Output Parameter:
3339: . v - the vector containing the integrals for each cost function
3341: Level: intermediate
3343: .seealso: TSSetCostIntegrand()
3345: .keywords: TS, sensitivity analysis
3346: @*/
3347: PetscErrorCode TSGetCostIntegral(TS ts,Vec *v)
3348: {
3352: *v = ts->vec_costintegral;
3353: return(0);
3354: }
3356: /*@
3357: TSComputeCostIntegrand - Evaluates the integral function in the cost functions.
3359: Input Parameters:
3360: + ts - the TS context
3361: . t - current time
3362: - y - state vector, i.e. current solution
3364: Output Parameter:
3365: . q - vector of size numcost to hold the outputs
3367: Note:
3368: Most users should not need to explicitly call this routine, as it
3369: is used internally within the sensitivity analysis context.
3371: Level: developer
3373: .keywords: TS, compute
3375: .seealso: TSSetCostIntegrand()
3376: @*/
3377: PetscErrorCode TSComputeCostIntegrand(TS ts,PetscReal t,Vec y,Vec q)
3378: {
3386: PetscLogEventBegin(TS_FunctionEval,ts,y,q,0);
3387: if (ts->costintegrand) {
3388: PetscStackPush("TS user integrand in the cost function");
3389: (*ts->costintegrand)(ts,t,y,q,ts->costintegrandctx);
3390: PetscStackPop;
3391: } else {
3392: VecZeroEntries(q);
3393: }
3395: PetscLogEventEnd(TS_FunctionEval,ts,y,q,0);
3396: return(0);
3397: }
3399: /*@
3400: TSAdjointComputeDRDYFunction - Runs the user-defined DRDY function.
3402: Collective on TS
3404: Input Parameters:
3405: . ts - The TS context obtained from TSCreate()
3407: Notes:
3408: TSAdjointComputeDRDYFunction() is typically used for sensitivity implementation,
3409: so most users would not generally call this routine themselves.
3411: Level: developer
3413: .keywords: TS, sensitivity
3414: .seealso: TSAdjointComputeDRDYFunction()
3415: @*/
3416: PetscErrorCode TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec y,Vec *drdy)
3417: {
3424: PetscStackPush("TS user DRDY function for sensitivity analysis");
3425: (*ts->drdyfunction)(ts,t,y,drdy,ts->costintegrandctx);
3426: PetscStackPop;
3427: return(0);
3428: }
3430: /*@
3431: TSAdjointComputeDRDPFunction - Runs the user-defined DRDP function.
3433: Collective on TS
3435: Input Parameters:
3436: . ts - The TS context obtained from TSCreate()
3438: Notes:
3439: TSDRDPFunction() is typically used for sensitivity implementation,
3440: so most users would not generally call this routine themselves.
3442: Level: developer
3444: .keywords: TS, sensitivity
3445: .seealso: TSAdjointSetDRDPFunction()
3446: @*/
3447: PetscErrorCode TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec y,Vec *drdp)
3448: {
3455: PetscStackPush("TS user DRDP function for sensitivity analysis");
3456: (*ts->drdpfunction)(ts,t,y,drdp,ts->costintegrandctx);
3457: PetscStackPop;
3458: return(0);
3459: }
3461: /*@C
3462: TSSetPreStep - Sets the general-purpose function
3463: called once at the beginning of each time step.
3465: Logically Collective on TS
3467: Input Parameters:
3468: + ts - The TS context obtained from TSCreate()
3469: - func - The function
3471: Calling sequence of func:
3472: . func (TS ts);
3474: Level: intermediate
3476: .keywords: TS, timestep
3477: .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep(), TSRestartStep()
3478: @*/
3479: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
3480: {
3483: ts->prestep = func;
3484: return(0);
3485: }
3487: /*@
3488: TSPreStep - Runs the user-defined pre-step function.
3490: Collective on TS
3492: Input Parameters:
3493: . ts - The TS context obtained from TSCreate()
3495: Notes:
3496: TSPreStep() is typically used within time stepping implementations,
3497: so most users would not generally call this routine themselves.
3499: Level: developer
3501: .keywords: TS, timestep
3502: .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep()
3503: @*/
3504: PetscErrorCode TSPreStep(TS ts)
3505: {
3510: if (ts->prestep) {
3511: Vec U;
3512: PetscObjectState sprev,spost;
3514: TSGetSolution(ts,&U);
3515: PetscObjectStateGet((PetscObject)U,&sprev);
3516: PetscStackCallStandard((*ts->prestep),(ts));
3517: PetscObjectStateGet((PetscObject)U,&spost);
3518: if (sprev != spost) {TSRestartStep(ts);}
3519: }
3520: return(0);
3521: }
3523: /*@C
3524: TSSetPreStage - Sets the general-purpose function
3525: called once at the beginning of each stage.
3527: Logically Collective on TS
3529: Input Parameters:
3530: + ts - The TS context obtained from TSCreate()
3531: - func - The function
3533: Calling sequence of func:
3534: . PetscErrorCode func(TS ts, PetscReal stagetime);
3536: Level: intermediate
3538: Note:
3539: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3540: The time step number being computed can be queried using TSGetStepNumber() and the total size of the step being
3541: attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime().
3543: .keywords: TS, timestep
3544: .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3545: @*/
3546: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal))
3547: {
3550: ts->prestage = func;
3551: return(0);
3552: }
3554: /*@C
3555: TSSetPostStage - Sets the general-purpose function
3556: called once at the end of each stage.
3558: Logically Collective on TS
3560: Input Parameters:
3561: + ts - The TS context obtained from TSCreate()
3562: - func - The function
3564: Calling sequence of func:
3565: . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y);
3567: Level: intermediate
3569: Note:
3570: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3571: The time step number being computed can be queried using TSGetStepNumber() and the total size of the step being
3572: attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime().
3574: .keywords: TS, timestep
3575: .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3576: @*/
3577: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*))
3578: {
3581: ts->poststage = func;
3582: return(0);
3583: }
3585: /*@C
3586: TSSetPostEvaluate - Sets the general-purpose function
3587: called once at the end of each step evaluation.
3589: Logically Collective on TS
3591: Input Parameters:
3592: + ts - The TS context obtained from TSCreate()
3593: - func - The function
3595: Calling sequence of func:
3596: . PetscErrorCode func(TS ts);
3598: Level: intermediate
3600: Note:
3601: Semantically, TSSetPostEvaluate() differs from TSSetPostStep() since the function it sets is called before event-handling
3602: thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, TSPostStep()
3603: may be passed a different solution, possibly changed by the event handler. TSPostEvaluate() is called after the next step
3604: solution is evaluated allowing to modify it, if need be. The solution can be obtained with TSGetSolution(), the time step
3605: with TSGetTimeStep(), and the time at the start of the step is available via TSGetTime()
3607: .keywords: TS, timestep
3608: .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3609: @*/
3610: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS))
3611: {
3614: ts->postevaluate = func;
3615: return(0);
3616: }
3618: /*@
3619: TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage()
3621: Collective on TS
3623: Input Parameters:
3624: . ts - The TS context obtained from TSCreate()
3625: stagetime - The absolute time of the current stage
3627: Notes:
3628: TSPreStage() is typically used within time stepping implementations,
3629: most users would not generally call this routine themselves.
3631: Level: developer
3633: .keywords: TS, timestep
3634: .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep()
3635: @*/
3636: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3637: {
3642: if (ts->prestage) {
3643: PetscStackCallStandard((*ts->prestage),(ts,stagetime));
3644: }
3645: return(0);
3646: }
3648: /*@
3649: TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage()
3651: Collective on TS
3653: Input Parameters:
3654: . ts - The TS context obtained from TSCreate()
3655: stagetime - The absolute time of the current stage
3656: stageindex - Stage number
3657: Y - Array of vectors (of size = total number
3658: of stages) with the stage solutions
3660: Notes:
3661: TSPostStage() is typically used within time stepping implementations,
3662: most users would not generally call this routine themselves.
3664: Level: developer
3666: .keywords: TS, timestep
3667: .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep()
3668: @*/
3669: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3670: {
3675: if (ts->poststage) {
3676: PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y));
3677: }
3678: return(0);
3679: }
3681: /*@
3682: TSPostEvaluate - Runs the user-defined post-evaluate function set using TSSetPostEvaluate()
3684: Collective on TS
3686: Input Parameters:
3687: . ts - The TS context obtained from TSCreate()
3689: Notes:
3690: TSPostEvaluate() is typically used within time stepping implementations,
3691: most users would not generally call this routine themselves.
3693: Level: developer
3695: .keywords: TS, timestep
3696: .seealso: TSSetPostEvaluate(), TSSetPreStep(), TSPreStep(), TSPostStep()
3697: @*/
3698: PetscErrorCode TSPostEvaluate(TS ts)
3699: {
3704: if (ts->postevaluate) {
3705: Vec U;
3706: PetscObjectState sprev,spost;
3708: TSGetSolution(ts,&U);
3709: PetscObjectStateGet((PetscObject)U,&sprev);
3710: PetscStackCallStandard((*ts->postevaluate),(ts));
3711: PetscObjectStateGet((PetscObject)U,&spost);
3712: if (sprev != spost) {TSRestartStep(ts);}
3713: }
3714: return(0);
3715: }
3717: /*@C
3718: TSSetPostStep - Sets the general-purpose function
3719: called once at the end of each time step.
3721: Logically Collective on TS
3723: Input Parameters:
3724: + ts - The TS context obtained from TSCreate()
3725: - func - The function
3727: Calling sequence of func:
3728: $ func (TS ts);
3730: Notes:
3731: The function set by TSSetPostStep() is called after each successful step. The solution vector X
3732: obtained by TSGetSolution() may be different than that computed at the step end if the event handler
3733: locates an event and TSPostEvent() modifies it. Use TSSetPostEvaluate() if an unmodified solution is needed instead.
3735: Level: intermediate
3737: .keywords: TS, timestep
3738: .seealso: TSSetPreStep(), TSSetPreStage(), TSSetPostEvaluate(), TSGetTimeStep(), TSGetStepNumber(), TSGetTime(), TSRestartStep()
3739: @*/
3740: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3741: {
3744: ts->poststep = func;
3745: return(0);
3746: }
3748: /*@
3749: TSPostStep - Runs the user-defined post-step function.
3751: Collective on TS
3753: Input Parameters:
3754: . ts - The TS context obtained from TSCreate()
3756: Notes:
3757: TSPostStep() is typically used within time stepping implementations,
3758: so most users would not generally call this routine themselves.
3760: Level: developer
3762: .keywords: TS, timestep
3763: @*/
3764: PetscErrorCode TSPostStep(TS ts)
3765: {
3770: if (ts->poststep) {
3771: Vec U;
3772: PetscObjectState sprev,spost;
3774: TSGetSolution(ts,&U);
3775: PetscObjectStateGet((PetscObject)U,&sprev);
3776: PetscStackCallStandard((*ts->poststep),(ts));
3777: PetscObjectStateGet((PetscObject)U,&spost);
3778: if (sprev != spost) {TSRestartStep(ts);}
3779: }
3780: return(0);
3781: }
3783: /* ------------ Routines to set performance monitoring options ----------- */
3785: /*@C
3786: TSMonitorSet - Sets an ADDITIONAL function that is to be used at every
3787: timestep to display the iteration's progress.
3789: Logically Collective on TS
3791: Input Parameters:
3792: + ts - the TS context obtained from TSCreate()
3793: . monitor - monitoring routine
3794: . mctx - [optional] user-defined context for private data for the
3795: monitor routine (use NULL if no context is desired)
3796: - monitordestroy - [optional] routine that frees monitor context
3797: (may be NULL)
3799: Calling sequence of monitor:
3800: $ PetscErrorCode monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
3802: + ts - the TS context
3803: . steps - iteration number (after the final time step the monitor routine may be called with a step of -1, this indicates the solution has been interpolated to this time)
3804: . time - current time
3805: . u - current iterate
3806: - mctx - [optional] monitoring context
3808: Notes:
3809: This routine adds an additional monitor to the list of monitors that
3810: already has been loaded.
3812: Fortran notes: Only a single monitor function can be set for each TS object
3814: Level: intermediate
3816: .keywords: TS, timestep, set, monitor
3818: .seealso: TSMonitorDefault(), TSMonitorCancel()
3819: @*/
3820: PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**))
3821: {
3823: PetscInt i;
3824: PetscBool identical;
3828: for (i=0; i<ts->numbermonitors;i++) {
3829: PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,mdestroy,(PetscErrorCode (*)(void))ts->monitor[i],ts->monitorcontext[i],ts->monitordestroy[i],&identical);
3830: if (identical) return(0);
3831: }
3832: if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3833: ts->monitor[ts->numbermonitors] = monitor;
3834: ts->monitordestroy[ts->numbermonitors] = mdestroy;
3835: ts->monitorcontext[ts->numbermonitors++] = (void*)mctx;
3836: return(0);
3837: }
3839: /*@C
3840: TSMonitorCancel - Clears all the monitors that have been set on a time-step object.
3842: Logically Collective on TS
3844: Input Parameters:
3845: . ts - the TS context obtained from TSCreate()
3847: Notes:
3848: There is no way to remove a single, specific monitor.
3850: Level: intermediate
3852: .keywords: TS, timestep, set, monitor
3854: .seealso: TSMonitorDefault(), TSMonitorSet()
3855: @*/
3856: PetscErrorCode TSMonitorCancel(TS ts)
3857: {
3859: PetscInt i;
3863: for (i=0; i<ts->numbermonitors; i++) {
3864: if (ts->monitordestroy[i]) {
3865: (*ts->monitordestroy[i])(&ts->monitorcontext[i]);
3866: }
3867: }
3868: ts->numbermonitors = 0;
3869: return(0);
3870: }
3872: /*@C
3873: TSMonitorDefault - The Default monitor, prints the timestep and time for each step
3875: Level: intermediate
3877: .keywords: TS, set, monitor
3879: .seealso: TSMonitorSet()
3880: @*/
3881: PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscViewerAndFormat *vf)
3882: {
3884: PetscViewer viewer = vf->viewer;
3885: PetscBool iascii,ibinary;
3889: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
3890: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
3891: PetscViewerPushFormat(viewer,vf->format);
3892: if (iascii) {
3893: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
3894: if (step == -1){ /* this indicates it is an interpolated solution */
3895: PetscViewerASCIIPrintf(viewer,"Interpolated solution at time %g between steps %D and %D\n",(double)ptime,ts->steps-1,ts->steps);
3896: } else {
3897: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
3898: }
3899: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
3900: } else if (ibinary) {
3901: PetscMPIInt rank;
3902: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);
3903: if (!rank) {
3904: PetscBool skipHeader;
3905: PetscInt classid = REAL_FILE_CLASSID;
3907: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
3908: if (!skipHeader) {
3909: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
3910: }
3911: PetscRealView(1,&ptime,viewer);
3912: } else {
3913: PetscRealView(0,&ptime,viewer);
3914: }
3915: }
3916: PetscViewerPopFormat(viewer);
3917: return(0);
3918: }
3920: /*@C
3921: TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
3922: timestep to display the iteration's progress.
3924: Logically Collective on TS
3926: Input Parameters:
3927: + ts - the TS context obtained from TSCreate()
3928: . adjointmonitor - monitoring routine
3929: . adjointmctx - [optional] user-defined context for private data for the
3930: monitor routine (use NULL if no context is desired)
3931: - adjointmonitordestroy - [optional] routine that frees monitor context
3932: (may be NULL)
3934: Calling sequence of monitor:
3935: $ int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)
3937: + ts - the TS context
3938: . 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
3939: been interpolated to)
3940: . time - current time
3941: . u - current iterate
3942: . numcost - number of cost functionos
3943: . lambda - sensitivities to initial conditions
3944: . mu - sensitivities to parameters
3945: - adjointmctx - [optional] adjoint monitoring context
3947: Notes:
3948: This routine adds an additional monitor to the list of monitors that
3949: already has been loaded.
3951: Fortran notes: Only a single monitor function can be set for each TS object
3953: Level: intermediate
3955: .keywords: TS, timestep, set, adjoint, monitor
3957: .seealso: TSAdjointMonitorCancel()
3958: @*/
3959: PetscErrorCode TSAdjointMonitorSet(TS ts,PetscErrorCode (*adjointmonitor)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*),void *adjointmctx,PetscErrorCode (*adjointmdestroy)(void**))
3960: {
3962: PetscInt i;
3963: PetscBool identical;
3967: for (i=0; i<ts->numbermonitors;i++) {
3968: PetscMonitorCompare((PetscErrorCode (*)(void))adjointmonitor,adjointmctx,adjointmdestroy,(PetscErrorCode (*)(void))ts->adjointmonitor[i],ts->adjointmonitorcontext[i],ts->adjointmonitordestroy[i],&identical);
3969: if (identical) return(0);
3970: }
3971: if (ts->numberadjointmonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many adjoint monitors set");
3972: ts->adjointmonitor[ts->numberadjointmonitors] = adjointmonitor;
3973: ts->adjointmonitordestroy[ts->numberadjointmonitors] = adjointmdestroy;
3974: ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void*)adjointmctx;
3975: return(0);
3976: }
3978: /*@C
3979: TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.
3981: Logically Collective on TS
3983: Input Parameters:
3984: . ts - the TS context obtained from TSCreate()
3986: Notes:
3987: There is no way to remove a single, specific monitor.
3989: Level: intermediate
3991: .keywords: TS, timestep, set, adjoint, monitor
3993: .seealso: TSAdjointMonitorSet()
3994: @*/
3995: PetscErrorCode TSAdjointMonitorCancel(TS ts)
3996: {
3998: PetscInt i;
4002: for (i=0; i<ts->numberadjointmonitors; i++) {
4003: if (ts->adjointmonitordestroy[i]) {
4004: (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);
4005: }
4006: }
4007: ts->numberadjointmonitors = 0;
4008: return(0);
4009: }
4011: /*@C
4012: TSAdjointMonitorDefault - the default monitor of adjoint computations
4014: Level: intermediate
4016: .keywords: TS, set, monitor
4018: .seealso: TSAdjointMonitorSet()
4019: @*/
4020: PetscErrorCode TSAdjointMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
4021: {
4023: PetscViewer viewer = vf->viewer;
4027: PetscViewerPushFormat(viewer,vf->format);
4028: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
4029: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
4030: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
4031: PetscViewerPopFormat(viewer);
4032: return(0);
4033: }
4035: /*@
4036: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
4038: Collective on TS
4040: Input Argument:
4041: + ts - time stepping context
4042: - t - time to interpolate to
4044: Output Argument:
4045: . U - state at given time
4047: Level: intermediate
4049: Developer Notes:
4050: TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
4052: .keywords: TS, set
4054: .seealso: TSSetExactFinalTime(), TSSolve()
4055: @*/
4056: PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U)
4057: {
4063: if (t < ts->ptime_prev || t > ts->ptime) SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Requested time %g not in last time steps [%g,%g]",t,(double)ts->ptime_prev,(double)ts->ptime);
4064: if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name);
4065: (*ts->ops->interpolate)(ts,t,U);
4066: return(0);
4067: }
4069: /*@
4070: TSStep - Steps one time step
4072: Collective on TS
4074: Input Parameter:
4075: . ts - the TS context obtained from TSCreate()
4077: Level: developer
4079: Notes:
4080: The public interface for the ODE/DAE solvers is TSSolve(), you should almost for sure be using that routine and not this routine.
4082: The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may
4083: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
4085: This may over-step the final time provided in TSSetMaxTime() depending on the time-step used. TSSolve() interpolates to exactly the
4086: time provided in TSSetMaxTime(). One can use TSInterpolate() to determine an interpolated solution within the final timestep.
4088: .keywords: TS, timestep, solve
4090: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate()
4091: @*/
4092: PetscErrorCode TSStep(TS ts)
4093: {
4094: PetscErrorCode ierr;
4095: static PetscBool cite = PETSC_FALSE;
4096: PetscReal ptime;
4100: PetscCitationsRegister("@techreport{tspaper,\n"
4101: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
4102: " author = {Shrirang Abhyankar and Jed Brown and Emil Constantinescu and Debojyoti Ghosh and Barry F. Smith},\n"
4103: " type = {Preprint},\n"
4104: " number = {ANL/MCS-P5061-0114},\n"
4105: " institution = {Argonne National Laboratory},\n"
4106: " year = {2014}\n}\n",&cite);
4108: TSSetUp(ts);
4109: TSTrajectorySetUp(ts->trajectory,ts);
4111: if (ts->max_time >= PETSC_MAX_REAL && ts->max_steps == PETSC_MAX_INT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
4112: if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"You must call TSSetExactFinalTime() or use -ts_exact_final_time <stepover,interpolate,matchstep> before calling TSStep()");
4113: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP && !ts->adapt) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Since TS is not adaptive you cannot use TS_EXACTFINALTIME_MATCHSTEP, suggest TS_EXACTFINALTIME_INTERPOLATE");
4115: if (!ts->steps) ts->ptime_prev = ts->ptime;
4116: ptime = ts->ptime; ts->ptime_prev_rollback = ts->ptime_prev;
4117: ts->reason = TS_CONVERGED_ITERATING;
4118: if (!ts->ops->step) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSStep not implemented for type '%s'",((PetscObject)ts)->type_name);
4119: PetscLogEventBegin(TS_Step,ts,0,0,0);
4120: (*ts->ops->step)(ts);
4121: PetscLogEventEnd(TS_Step,ts,0,0,0);
4122: ts->ptime_prev = ptime;
4123: ts->steps++;
4124: ts->steprollback = PETSC_FALSE;
4125: ts->steprestart = PETSC_FALSE;
4127: if (ts->reason < 0) {
4128: if (ts->errorifstepfailed) {
4129: if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) 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]);
4130: else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
4131: }
4132: } else if (!ts->reason) {
4133: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4134: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4135: }
4136: return(0);
4137: }
4139: /*@
4140: TSAdjointStep - Steps one time step backward in the adjoint run
4142: Collective on TS
4144: Input Parameter:
4145: . ts - the TS context obtained from TSCreate()
4147: Level: intermediate
4149: .keywords: TS, adjoint, step
4151: .seealso: TSAdjointSetUp(), TSAdjointSolve()
4152: @*/
4153: PetscErrorCode TSAdjointStep(TS ts)
4154: {
4155: DM dm;
4156: PetscErrorCode ierr;
4160: TSGetDM(ts,&dm);
4161: TSAdjointSetUp(ts);
4163: VecViewFromOptions(ts->vec_sol,(PetscObject)ts,"-ts_view_solution");
4165: ts->reason = TS_CONVERGED_ITERATING;
4166: ts->ptime_prev = ts->ptime;
4167: if (!ts->ops->adjointstep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed because the adjoint of %s has not been implemented, try other time stepping methods for adjoint sensitivity analysis",((PetscObject)ts)->type_name);
4168: PetscLogEventBegin(TS_AdjointStep,ts,0,0,0);
4169: (*ts->ops->adjointstep)(ts);
4170: PetscLogEventEnd(TS_AdjointStep,ts,0,0,0);
4171: ts->adjoint_steps++; ts->steps--;
4173: if (ts->reason < 0) {
4174: if (ts->errorifstepfailed) {
4175: if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) 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]);
4176: else if (ts->reason == TS_DIVERGED_STEP_REJECTED) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s, increase -ts_max_reject or make negative to attempt recovery",TSConvergedReasons[ts->reason]);
4177: else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
4178: }
4179: } else if (!ts->reason) {
4180: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
4181: }
4182: return(0);
4183: }
4185: /*@
4186: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
4187: at the end of a time step with a given order of accuracy.
4189: Collective on TS
4191: Input Arguments:
4192: + ts - time stepping context
4193: . wnormtype - norm type, either NORM_2 or NORM_INFINITY
4194: - order - optional, desired order for the error evaluation or PETSC_DECIDE
4196: Output Arguments:
4197: + order - optional, the actual order of the error evaluation
4198: - wlte - the weighted local truncation error norm
4200: Level: advanced
4202: Notes:
4203: If the timestepper cannot evaluate the error in a particular step
4204: (eg. in the first step or restart steps after event handling),
4205: this routine returns wlte=-1.0 .
4207: .seealso: TSStep(), TSAdapt, TSErrorWeightedNorm()
4208: @*/
4209: PetscErrorCode TSEvaluateWLTE(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
4210: {
4220: if (wnormtype != NORM_2 && wnormtype != NORM_INFINITY) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
4221: if (!ts->ops->evaluatewlte) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateWLTE not implemented for type '%s'",((PetscObject)ts)->type_name);
4222: (*ts->ops->evaluatewlte)(ts,wnormtype,order,wlte);
4223: return(0);
4224: }
4226: /*@
4227: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
4229: Collective on TS
4231: Input Arguments:
4232: + ts - time stepping context
4233: . order - desired order of accuracy
4234: - done - whether the step was evaluated at this order (pass NULL to generate an error if not available)
4236: Output Arguments:
4237: . U - state at the end of the current step
4239: Level: advanced
4241: Notes:
4242: This function cannot be called until all stages have been evaluated.
4243: 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.
4245: .seealso: TSStep(), TSAdapt
4246: @*/
4247: PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done)
4248: {
4255: if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name);
4256: (*ts->ops->evaluatestep)(ts,order,U,done);
4257: return(0);
4258: }
4260: /*@
4261: TSForwardCostIntegral - Evaluate the cost integral in the forward run.
4263: Collective on TS
4265: Input Arguments:
4266: . ts - time stepping context
4268: Level: advanced
4270: Notes:
4271: This function cannot be called until TSStep() has been completed.
4273: .seealso: TSSolve(), TSAdjointCostIntegral()
4274: @*/
4275: PetscErrorCode TSForwardCostIntegral(TS ts)
4276: {
4279: if (!ts->ops->forwardintegral) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide integral evaluation in the forward run",((PetscObject)ts)->type_name);
4280: (*ts->ops->forwardintegral)(ts);
4281: return(0);
4282: }
4284: /*@
4285: TSSolve - Steps the requested number of timesteps.
4287: Collective on TS
4289: Input Parameter:
4290: + ts - the TS context obtained from TSCreate()
4291: - u - the solution vector (can be null if TSSetSolution() was used and TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP) was not used,
4292: otherwise must contain the initial conditions and will contain the solution at the final requested time
4294: Level: beginner
4296: Notes:
4297: The final time returned by this function may be different from the time of the internally
4298: held state accessible by TSGetSolution() and TSGetTime() because the method may have
4299: stepped over the final time.
4301: .keywords: TS, timestep, solve
4303: .seealso: TSCreate(), TSSetSolution(), TSStep(), TSGetTime(), TSGetSolveTime()
4304: @*/
4305: PetscErrorCode TSSolve(TS ts,Vec u)
4306: {
4307: Vec solution;
4308: PetscErrorCode ierr;
4314: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && u) { /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */
4315: if (!ts->vec_sol || u == ts->vec_sol) {
4316: VecDuplicate(u,&solution);
4317: TSSetSolution(ts,solution);
4318: VecDestroy(&solution); /* grant ownership */
4319: }
4320: VecCopy(u,ts->vec_sol);
4321: if (ts->forward_solve) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
4322: } else if (u) {
4323: TSSetSolution(ts,u);
4324: }
4325: TSSetUp(ts);
4326: TSTrajectorySetUp(ts->trajectory,ts);
4328: if (ts->max_time >= PETSC_MAX_REAL && ts->max_steps == PETSC_MAX_INT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
4329: if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"You must call TSSetExactFinalTime() or use -ts_exact_final_time <stepover,interpolate,matchstep> before calling TSSolve()");
4330: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP && !ts->adapt) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Since TS is not adaptive you cannot use TS_EXACTFINALTIME_MATCHSTEP, suggest TS_EXACTFINALTIME_INTERPOLATE");
4332: if (ts->forward_solve) {
4333: TSForwardSetUp(ts);
4334: }
4336: /* reset number of steps only when the step is not restarted. ARKIMEX
4337: restarts the step after an event. Resetting these counters in such case causes
4338: TSTrajectory to incorrectly save the output files
4339: */
4340: /* reset time step and iteration counters */
4341: if (!ts->steps) {
4342: ts->ksp_its = 0;
4343: ts->snes_its = 0;
4344: ts->num_snes_failures = 0;
4345: ts->reject = 0;
4346: ts->steprestart = PETSC_TRUE;
4347: ts->steprollback = PETSC_FALSE;
4348: }
4349: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP && ts->ptime + ts->time_step > ts->max_time) ts->time_step = ts->max_time - ts->ptime;
4350: ts->reason = TS_CONVERGED_ITERATING;
4352: TSViewFromOptions(ts,NULL,"-ts_view_pre");
4354: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4355: (*ts->ops->solve)(ts);
4356: if (u) {VecCopy(ts->vec_sol,u);}
4357: ts->solvetime = ts->ptime;
4358: solution = ts->vec_sol;
4359: } else { /* Step the requested number of timesteps. */
4360: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4361: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4363: if (!ts->steps) {
4364: TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);
4365: TSEventInitialize(ts->event,ts,ts->ptime,ts->vec_sol);
4366: }
4368: while (!ts->reason) {
4369: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
4370: if (!ts->steprollback) {
4371: TSPreStep(ts);
4372: }
4373: TSStep(ts);
4374: if (ts->vec_costintegral && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4375: TSForwardCostIntegral(ts);
4376: }
4377: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4378: TSForwardStep(ts);
4379: }
4380: TSPostEvaluate(ts);
4381: TSEventHandler(ts); /* The right-hand side may be changed due to event. Be careful with Any computation using the RHS information after this point. */
4382: if (ts->steprollback) {
4383: TSPostEvaluate(ts);
4384: }
4385: if (!ts->steprollback) {
4386: TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);
4387: TSPostStep(ts);
4388: }
4389: }
4390: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
4392: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4393: TSInterpolate(ts,ts->max_time,u);
4394: ts->solvetime = ts->max_time;
4395: solution = u;
4396: TSMonitor(ts,-1,ts->solvetime,solution);
4397: } else {
4398: if (u) {VecCopy(ts->vec_sol,u);}
4399: ts->solvetime = ts->ptime;
4400: solution = ts->vec_sol;
4401: }
4402: }
4404: TSViewFromOptions(ts,NULL,"-ts_view");
4405: VecViewFromOptions(solution,NULL,"-ts_view_solution");
4406: PetscObjectSAWsBlock((PetscObject)ts);
4407: if (ts->adjoint_solve) {
4408: TSAdjointSolve(ts);
4409: }
4410: return(0);
4411: }
4413: /*@
4414: TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
4416: Collective on TS
4418: Input Arguments:
4419: . ts - time stepping context
4421: Level: advanced
4423: Notes:
4424: This function cannot be called until TSAdjointStep() has been completed.
4426: .seealso: TSAdjointSolve(), TSAdjointStep
4427: @*/
4428: PetscErrorCode TSAdjointCostIntegral(TS ts)
4429: {
4432: if (!ts->ops->adjointintegral) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide integral evaluation in the adjoint run",((PetscObject)ts)->type_name);
4433: (*ts->ops->adjointintegral)(ts);
4434: return(0);
4435: }
4437: /*@
4438: TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE
4440: Collective on TS
4442: Input Parameter:
4443: . ts - the TS context obtained from TSCreate()
4445: Options Database:
4446: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values
4448: Level: intermediate
4450: Notes:
4451: This must be called after a call to TSSolve() that solves the forward problem
4453: By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time
4455: .keywords: TS, timestep, solve
4457: .seealso: TSCreate(), TSSetCostGradients(), TSSetSolution(), TSAdjointStep()
4458: @*/
4459: PetscErrorCode TSAdjointSolve(TS ts)
4460: {
4461: PetscErrorCode ierr;
4465: TSAdjointSetUp(ts);
4467: /* reset time step and iteration counters */
4468: ts->adjoint_steps = 0;
4469: ts->ksp_its = 0;
4470: ts->snes_its = 0;
4471: ts->num_snes_failures = 0;
4472: ts->reject = 0;
4473: ts->reason = TS_CONVERGED_ITERATING;
4475: if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
4476: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
4478: while (!ts->reason) {
4479: TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);
4480: TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
4481: TSAdjointEventHandler(ts);
4482: TSAdjointStep(ts);
4483: if (ts->vec_costintegral && !ts->costintegralfwd) {
4484: TSAdjointCostIntegral(ts);
4485: }
4486: }
4487: TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);
4488: TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
4489: ts->solvetime = ts->ptime;
4490: TSTrajectoryViewFromOptions(ts->trajectory,NULL,"-ts_trajectory_view");
4491: VecViewFromOptions(ts->vecs_sensi[0],(PetscObject) ts, "-ts_adjoint_view_solution");
4492: return(0);
4493: }
4495: /*@C
4496: TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet()
4498: Collective on TS
4500: Input Parameters:
4501: + ts - time stepping context obtained from TSCreate()
4502: . step - step number that has just completed
4503: . ptime - model time of the state
4504: - u - state at the current model time
4506: Notes:
4507: TSMonitor() is typically used automatically within the time stepping implementations.
4508: Users would almost never call this routine directly.
4510: A step of -1 indicates that the monitor is being called on a solution obtained by interpolating from computed solutions
4512: Level: developer
4514: .keywords: TS, timestep
4515: @*/
4516: PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u)
4517: {
4518: DM dm;
4519: PetscInt i,n = ts->numbermonitors;
4526: TSGetDM(ts,&dm);
4527: DMSetOutputSequenceNumber(dm,step,ptime);
4529: VecLockPush(u);
4530: for (i=0; i<n; i++) {
4531: (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);
4532: }
4533: VecLockPop(u);
4534: return(0);
4535: }
4537: /*@C
4538: TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using TSAdjointMonitorSet()
4540: Collective on TS
4542: Input Parameters:
4543: + ts - time stepping context obtained from TSCreate()
4544: . step - step number that has just completed
4545: . ptime - model time of the state
4546: . u - state at the current model time
4547: . numcost - number of cost functions (dimension of lambda or mu)
4548: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
4549: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
4551: Notes:
4552: TSAdjointMonitor() is typically used automatically within the time stepping implementations.
4553: Users would almost never call this routine directly.
4555: Level: developer
4557: .keywords: TS, timestep
4558: @*/
4559: PetscErrorCode TSAdjointMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda, Vec *mu)
4560: {
4562: PetscInt i,n = ts->numberadjointmonitors;
4567: VecLockPush(u);
4568: for (i=0; i<n; i++) {
4569: (*ts->adjointmonitor[i])(ts,step,ptime,u,numcost,lambda,mu,ts->adjointmonitorcontext[i]);
4570: }
4571: VecLockPop(u);
4572: return(0);
4573: }
4575: /* ------------------------------------------------------------------------*/
4576: /*@C
4577: TSMonitorLGCtxCreate - Creates a TSMonitorLGCtx context for use with
4578: TS to monitor the solution process graphically in various ways
4580: Collective on TS
4582: Input Parameters:
4583: + host - the X display to open, or null for the local machine
4584: . label - the title to put in the title bar
4585: . x, y - the screen coordinates of the upper left coordinate of the window
4586: . m, n - the screen width and height in pixels
4587: - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time
4589: Output Parameter:
4590: . ctx - the context
4592: Options Database Key:
4593: + -ts_monitor_lg_timestep - automatically sets line graph monitor
4594: + -ts_monitor_lg_timestep_log - automatically sets line graph monitor
4595: . -ts_monitor_lg_solution - monitor the solution (or certain values of the solution by calling TSMonitorLGSetDisplayVariables() or TSMonitorLGCtxSetDisplayVariables())
4596: . -ts_monitor_lg_error - monitor the error
4597: . -ts_monitor_lg_ksp_iterations - monitor the number of KSP iterations needed for each timestep
4598: . -ts_monitor_lg_snes_iterations - monitor the number of SNES iterations needed for each timestep
4599: - -lg_use_markers <true,false> - mark the data points (at each time step) on the plot; default is true
4601: Notes:
4602: Use TSMonitorLGCtxDestroy() to destroy.
4604: One can provide a function that transforms the solution before plotting it with TSMonitorLGCtxSetTransform() or TSMonitorLGSetTransform()
4606: Many of the functions that control the monitoring have two forms: TSMonitorLGSet/GetXXXX() and TSMonitorLGCtxSet/GetXXXX() the first take a TS object as the
4607: first argument (if that TS object does not have a TSMonitorLGCtx associated with it the function call is ignored) and the second takes a TSMonitorLGCtx object
4608: as the first argument.
4610: One can control the names displayed for each solution or error variable with TSMonitorLGCtxSetVariableNames() or TSMonitorLGSetVariableNames()
4612: Level: intermediate
4614: .keywords: TS, monitor, line graph, residual
4616: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError(), TSMonitorDefault(), VecView(),
4617: TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
4618: TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
4619: TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
4620: TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()
4622: @*/
4623: PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx)
4624: {
4625: PetscDraw draw;
4629: PetscNew(ctx);
4630: PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
4631: PetscDrawSetFromOptions(draw);
4632: PetscDrawLGCreate(draw,1,&(*ctx)->lg);
4633: PetscDrawLGSetFromOptions((*ctx)->lg);
4634: PetscDrawDestroy(&draw);
4635: (*ctx)->howoften = howoften;
4636: return(0);
4637: }
4639: PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
4640: {
4641: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
4642: PetscReal x = ptime,y;
4646: if (step < 0) return(0); /* -1 indicates an interpolated solution */
4647: if (!step) {
4648: PetscDrawAxis axis;
4649: const char *ylabel = ctx->semilogy ? "Log Time Step" : "Time Step";
4650: PetscDrawLGGetAxis(ctx->lg,&axis);
4651: PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time",ylabel);
4652: PetscDrawLGReset(ctx->lg);
4653: }
4654: TSGetTimeStep(ts,&y);
4655: if (ctx->semilogy) y = PetscLog10Real(y);
4656: PetscDrawLGAddPoint(ctx->lg,&x,&y);
4657: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
4658: PetscDrawLGDraw(ctx->lg);
4659: PetscDrawLGSave(ctx->lg);
4660: }
4661: return(0);
4662: }
4664: /*@C
4665: TSMonitorLGCtxDestroy - Destroys a line graph context that was created
4666: with TSMonitorLGCtxCreate().
4668: Collective on TSMonitorLGCtx
4670: Input Parameter:
4671: . ctx - the monitor context
4673: Level: intermediate
4675: .keywords: TS, monitor, line graph, destroy
4677: .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep();
4678: @*/
4679: PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx)
4680: {
4684: if ((*ctx)->transformdestroy) {
4685: ((*ctx)->transformdestroy)((*ctx)->transformctx);
4686: }
4687: PetscDrawLGDestroy(&(*ctx)->lg);
4688: PetscStrArrayDestroy(&(*ctx)->names);
4689: PetscStrArrayDestroy(&(*ctx)->displaynames);
4690: PetscFree((*ctx)->displayvariables);
4691: PetscFree((*ctx)->displayvalues);
4692: PetscFree(*ctx);
4693: return(0);
4694: }
4696: /*@
4697: TSGetTime - Gets the time of the most recently completed step.
4699: Not Collective
4701: Input Parameter:
4702: . ts - the TS context obtained from TSCreate()
4704: Output Parameter:
4705: . t - the current time. This time may not corresponds to the final time set with TSSetMaxTime(), use TSGetSolveTime().
4707: Level: beginner
4709: Note:
4710: When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(),
4711: TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated.
4713: .seealso: TSGetSolveTime(), TSSetTime(), TSGetTimeStep()
4715: .keywords: TS, get, time
4716: @*/
4717: PetscErrorCode TSGetTime(TS ts,PetscReal *t)
4718: {
4722: *t = ts->ptime;
4723: return(0);
4724: }
4726: /*@
4727: TSGetPrevTime - Gets the starting time of the previously completed step.
4729: Not Collective
4731: Input Parameter:
4732: . ts - the TS context obtained from TSCreate()
4734: Output Parameter:
4735: . t - the previous time
4737: Level: beginner
4739: .seealso: TSGetTime(), TSGetSolveTime(), TSGetTimeStep()
4741: .keywords: TS, get, time
4742: @*/
4743: PetscErrorCode TSGetPrevTime(TS ts,PetscReal *t)
4744: {
4748: *t = ts->ptime_prev;
4749: return(0);
4750: }
4752: /*@
4753: TSSetTime - Allows one to reset the time.
4755: Logically Collective on TS
4757: Input Parameters:
4758: + ts - the TS context obtained from TSCreate()
4759: - time - the time
4761: Level: intermediate
4763: .seealso: TSGetTime(), TSSetMaxSteps()
4765: .keywords: TS, set, time
4766: @*/
4767: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4768: {
4772: ts->ptime = t;
4773: return(0);
4774: }
4776: /*@C
4777: TSSetOptionsPrefix - Sets the prefix used for searching for all
4778: TS options in the database.
4780: Logically Collective on TS
4782: Input Parameter:
4783: + ts - The TS context
4784: - prefix - The prefix to prepend to all option names
4786: Notes:
4787: A hyphen (-) must NOT be given at the beginning of the prefix name.
4788: The first character of all runtime options is AUTOMATICALLY the
4789: hyphen.
4791: Level: advanced
4793: .keywords: TS, set, options, prefix, database
4795: .seealso: TSSetFromOptions()
4797: @*/
4798: PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[])
4799: {
4801: SNES snes;
4805: PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);
4806: TSGetSNES(ts,&snes);
4807: SNESSetOptionsPrefix(snes,prefix);
4808: return(0);
4809: }
4811: /*@C
4812: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4813: TS options in the database.
4815: Logically Collective on TS
4817: Input Parameter:
4818: + ts - The TS context
4819: - prefix - The prefix to prepend to all option names
4821: Notes:
4822: A hyphen (-) must NOT be given at the beginning of the prefix name.
4823: The first character of all runtime options is AUTOMATICALLY the
4824: hyphen.
4826: Level: advanced
4828: .keywords: TS, append, options, prefix, database
4830: .seealso: TSGetOptionsPrefix()
4832: @*/
4833: PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[])
4834: {
4836: SNES snes;
4840: PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);
4841: TSGetSNES(ts,&snes);
4842: SNESAppendOptionsPrefix(snes,prefix);
4843: return(0);
4844: }
4846: /*@C
4847: TSGetOptionsPrefix - Sets the prefix used for searching for all
4848: TS options in the database.
4850: Not Collective
4852: Input Parameter:
4853: . ts - The TS context
4855: Output Parameter:
4856: . prefix - A pointer to the prefix string used
4858: Notes: On the fortran side, the user should pass in a string 'prifix' of
4859: sufficient length to hold the prefix.
4861: Level: intermediate
4863: .keywords: TS, get, options, prefix, database
4865: .seealso: TSAppendOptionsPrefix()
4866: @*/
4867: PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[])
4868: {
4874: PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);
4875: return(0);
4876: }
4878: /*@C
4879: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4881: Not Collective, but parallel objects are returned if TS is parallel
4883: Input Parameter:
4884: . ts - The TS context obtained from TSCreate()
4886: Output Parameters:
4887: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL)
4888: . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL)
4889: . func - Function to compute the Jacobian of the RHS (or NULL)
4890: - ctx - User-defined context for Jacobian evaluation routine (or NULL)
4892: Notes: You can pass in NULL for any return argument you do not need.
4894: Level: intermediate
4896: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()
4898: .keywords: TS, timestep, get, matrix, Jacobian
4899: @*/
4900: PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx)
4901: {
4903: DM dm;
4906: if (Amat || Pmat) {
4907: SNES snes;
4908: TSGetSNES(ts,&snes);
4909: SNESSetUpMatrices(snes);
4910: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
4911: }
4912: TSGetDM(ts,&dm);
4913: DMTSGetRHSJacobian(dm,func,ctx);
4914: return(0);
4915: }
4917: /*@C
4918: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4920: Not Collective, but parallel objects are returned if TS is parallel
4922: Input Parameter:
4923: . ts - The TS context obtained from TSCreate()
4925: Output Parameters:
4926: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4927: . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat
4928: . f - The function to compute the matrices
4929: - ctx - User-defined context for Jacobian evaluation routine
4931: Notes: You can pass in NULL for any return argument you do not need.
4933: Level: advanced
4935: .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()
4937: .keywords: TS, timestep, get, matrix, Jacobian
4938: @*/
4939: PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx)
4940: {
4942: DM dm;
4945: if (Amat || Pmat) {
4946: SNES snes;
4947: TSGetSNES(ts,&snes);
4948: SNESSetUpMatrices(snes);
4949: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
4950: }
4951: TSGetDM(ts,&dm);
4952: DMTSGetIJacobian(dm,f,ctx);
4953: return(0);
4954: }
4956: /*@C
4957: TSMonitorDrawSolution - Monitors progress of the TS solvers by calling
4958: VecView() for the solution at each timestep
4960: Collective on TS
4962: Input Parameters:
4963: + ts - the TS context
4964: . step - current time-step
4965: . ptime - current time
4966: - dummy - either a viewer or NULL
4968: Options Database:
4969: . -ts_monitor_draw_solution_initial - show initial solution as well as current solution
4971: Notes: the initial solution and current solution are not display with a common axis scaling so generally the option -ts_monitor_draw_solution_initial
4972: will look bad
4974: Level: intermediate
4976: .keywords: TS, vector, monitor, view
4978: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4979: @*/
4980: PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4981: {
4982: PetscErrorCode ierr;
4983: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
4984: PetscDraw draw;
4987: if (!step && ictx->showinitial) {
4988: if (!ictx->initialsolution) {
4989: VecDuplicate(u,&ictx->initialsolution);
4990: }
4991: VecCopy(u,ictx->initialsolution);
4992: }
4993: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);
4995: if (ictx->showinitial) {
4996: PetscReal pause;
4997: PetscViewerDrawGetPause(ictx->viewer,&pause);
4998: PetscViewerDrawSetPause(ictx->viewer,0.0);
4999: VecView(ictx->initialsolution,ictx->viewer);
5000: PetscViewerDrawSetPause(ictx->viewer,pause);
5001: PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);
5002: }
5003: VecView(u,ictx->viewer);
5004: if (ictx->showtimestepandtime) {
5005: PetscReal xl,yl,xr,yr,h;
5006: char time[32];
5008: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
5009: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
5010: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
5011: h = yl + .95*(yr - yl);
5012: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
5013: PetscDrawFlush(draw);
5014: }
5016: if (ictx->showinitial) {
5017: PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);
5018: }
5019: return(0);
5020: }
5022: /*@C
5023: TSAdjointMonitorDrawSensi - Monitors progress of the adjoint TS solvers by calling
5024: VecView() for the sensitivities to initial states at each timestep
5026: Collective on TS
5028: Input Parameters:
5029: + ts - the TS context
5030: . step - current time-step
5031: . ptime - current time
5032: . u - current state
5033: . numcost - number of cost functions
5034: . lambda - sensitivities to initial conditions
5035: . mu - sensitivities to parameters
5036: - dummy - either a viewer or NULL
5038: Level: intermediate
5040: .keywords: TS, vector, adjoint, monitor, view
5042: .seealso: TSAdjointMonitorSet(), TSAdjointMonitorDefault(), VecView()
5043: @*/
5044: PetscErrorCode TSAdjointMonitorDrawSensi(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda,Vec *mu,void *dummy)
5045: {
5046: PetscErrorCode ierr;
5047: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
5048: PetscDraw draw;
5049: PetscReal xl,yl,xr,yr,h;
5050: char time[32];
5053: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);
5055: VecView(lambda[0],ictx->viewer);
5056: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
5057: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
5058: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
5059: h = yl + .95*(yr - yl);
5060: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
5061: PetscDrawFlush(draw);
5062: return(0);
5063: }
5065: /*@C
5066: TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram
5068: Collective on TS
5070: Input Parameters:
5071: + ts - the TS context
5072: . step - current time-step
5073: . ptime - current time
5074: - dummy - either a viewer or NULL
5076: Level: intermediate
5078: .keywords: TS, vector, monitor, view
5080: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5081: @*/
5082: PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
5083: {
5084: PetscErrorCode ierr;
5085: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
5086: PetscDraw draw;
5087: PetscDrawAxis axis;
5088: PetscInt n;
5089: PetscMPIInt size;
5090: PetscReal U0,U1,xl,yl,xr,yr,h;
5091: char time[32];
5092: const PetscScalar *U;
5095: MPI_Comm_size(PetscObjectComm((PetscObject)ts),&size);
5096: if (size != 1) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Only allowed for sequential runs");
5097: VecGetSize(u,&n);
5098: if (n != 2) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Only for ODEs with two unknowns");
5100: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
5101: PetscViewerDrawGetDrawAxis(ictx->viewer,0,&axis);
5102: PetscDrawAxisGetLimits(axis,&xl,&xr,&yl,&yr);
5103: if (!step) {
5104: PetscDrawClear(draw);
5105: PetscDrawAxisDraw(axis);
5106: }
5108: VecGetArrayRead(u,&U);
5109: U0 = PetscRealPart(U[0]);
5110: U1 = PetscRealPart(U[1]);
5111: VecRestoreArrayRead(u,&U);
5112: if ((U0 < xl) || (U1 < yl) || (U0 > xr) || (U1 > yr)) return(0);
5114: PetscDrawCollectiveBegin(draw);
5115: PetscDrawPoint(draw,U0,U1,PETSC_DRAW_BLACK);
5116: if (ictx->showtimestepandtime) {
5117: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
5118: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
5119: h = yl + .95*(yr - yl);
5120: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
5121: }
5122: PetscDrawCollectiveEnd(draw);
5123: PetscDrawFlush(draw);
5124: PetscDrawPause(draw);
5125: PetscDrawSave(draw);
5126: return(0);
5127: }
5129: /*@C
5130: TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution()
5132: Collective on TS
5134: Input Parameters:
5135: . ctx - the monitor context
5137: Level: intermediate
5139: .keywords: TS, vector, monitor, view
5141: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError()
5142: @*/
5143: PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx)
5144: {
5148: PetscViewerDestroy(&(*ictx)->viewer);
5149: VecDestroy(&(*ictx)->initialsolution);
5150: PetscFree(*ictx);
5151: return(0);
5152: }
5154: /*@C
5155: TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx
5157: Collective on TS
5159: Input Parameter:
5160: . ts - time-step context
5162: Output Patameter:
5163: . ctx - the monitor context
5165: Options Database:
5166: . -ts_monitor_draw_solution_initial - show initial solution as well as current solution
5168: Level: intermediate
5170: .keywords: TS, vector, monitor, view
5172: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx()
5173: @*/
5174: PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx)
5175: {
5176: PetscErrorCode ierr;
5179: PetscNew(ctx);
5180: PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);
5181: PetscViewerSetFromOptions((*ctx)->viewer);
5183: (*ctx)->howoften = howoften;
5184: (*ctx)->showinitial = PETSC_FALSE;
5185: PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);
5187: (*ctx)->showtimestepandtime = PETSC_FALSE;
5188: PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);
5189: return(0);
5190: }
5192: /*@C
5193: TSMonitorDrawError - Monitors progress of the TS solvers by calling
5194: VecView() for the error at each timestep
5196: Collective on TS
5198: Input Parameters:
5199: + ts - the TS context
5200: . step - current time-step
5201: . ptime - current time
5202: - dummy - either a viewer or NULL
5204: Level: intermediate
5206: .keywords: TS, vector, monitor, view
5208: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5209: @*/
5210: PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
5211: {
5212: PetscErrorCode ierr;
5213: TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy;
5214: PetscViewer viewer = ctx->viewer;
5215: Vec work;
5218: if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return(0);
5219: VecDuplicate(u,&work);
5220: TSComputeSolutionFunction(ts,ptime,work);
5221: VecAXPY(work,-1.0,u);
5222: VecView(work,viewer);
5223: VecDestroy(&work);
5224: return(0);
5225: }
5227: #include <petsc/private/dmimpl.h>
5228: /*@
5229: TSSetDM - Sets the DM that may be used by some nonlinear solvers or preconditioners under the TS
5231: Logically Collective on TS and DM
5233: Input Parameters:
5234: + ts - the ODE integrator object
5235: - dm - the dm, cannot be NULL
5237: Level: intermediate
5239: .seealso: TSGetDM(), SNESSetDM(), SNESGetDM()
5240: @*/
5241: PetscErrorCode TSSetDM(TS ts,DM dm)
5242: {
5244: SNES snes;
5245: DMTS tsdm;
5250: PetscObjectReference((PetscObject)dm);
5251: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
5252: if (ts->dm->dmts && !dm->dmts) {
5253: DMCopyDMTS(ts->dm,dm);
5254: DMGetDMTS(ts->dm,&tsdm);
5255: if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */
5256: tsdm->originaldm = dm;
5257: }
5258: }
5259: DMDestroy(&ts->dm);
5260: }
5261: ts->dm = dm;
5263: TSGetSNES(ts,&snes);
5264: SNESSetDM(snes,dm);
5265: return(0);
5266: }
5268: /*@
5269: TSGetDM - Gets the DM that may be used by some preconditioners
5271: Not Collective
5273: Input Parameter:
5274: . ts - the preconditioner context
5276: Output Parameter:
5277: . dm - the dm
5279: Level: intermediate
5281: .seealso: TSSetDM(), SNESSetDM(), SNESGetDM()
5282: @*/
5283: PetscErrorCode TSGetDM(TS ts,DM *dm)
5284: {
5289: if (!ts->dm) {
5290: DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);
5291: if (ts->snes) {SNESSetDM(ts->snes,ts->dm);}
5292: }
5293: *dm = ts->dm;
5294: return(0);
5295: }
5297: /*@
5298: SNESTSFormFunction - Function to evaluate nonlinear residual
5300: Logically Collective on SNES
5302: Input Parameter:
5303: + snes - nonlinear solver
5304: . U - the current state at which to evaluate the residual
5305: - ctx - user context, must be a TS
5307: Output Parameter:
5308: . F - the nonlinear residual
5310: Notes:
5311: This function is not normally called by users and is automatically registered with the SNES used by TS.
5312: It is most frequently passed to MatFDColoringSetFunction().
5314: Level: advanced
5316: .seealso: SNESSetFunction(), MatFDColoringSetFunction()
5317: @*/
5318: PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx)
5319: {
5320: TS ts = (TS)ctx;
5328: (ts->ops->snesfunction)(snes,U,F,ts);
5329: return(0);
5330: }
5332: /*@
5333: SNESTSFormJacobian - Function to evaluate the Jacobian
5335: Collective on SNES
5337: Input Parameter:
5338: + snes - nonlinear solver
5339: . U - the current state at which to evaluate the residual
5340: - ctx - user context, must be a TS
5342: Output Parameter:
5343: + A - the Jacobian
5344: . B - the preconditioning matrix (may be the same as A)
5345: - flag - indicates any structure change in the matrix
5347: Notes:
5348: This function is not normally called by users and is automatically registered with the SNES used by TS.
5350: Level: developer
5352: .seealso: SNESSetJacobian()
5353: @*/
5354: PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx)
5355: {
5356: TS ts = (TS)ctx;
5367: (ts->ops->snesjacobian)(snes,U,A,B,ts);
5368: return(0);
5369: }
5371: /*@C
5372: TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems Udot = A U only
5374: Collective on TS
5376: Input Arguments:
5377: + ts - time stepping context
5378: . t - time at which to evaluate
5379: . U - state at which to evaluate
5380: - ctx - context
5382: Output Arguments:
5383: . F - right hand side
5385: Level: intermediate
5387: Notes:
5388: This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems.
5389: The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian().
5391: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
5392: @*/
5393: PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
5394: {
5396: Mat Arhs,Brhs;
5399: TSGetRHSMats_Private(ts,&Arhs,&Brhs);
5400: TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);
5401: MatMult(Arhs,U,F);
5402: return(0);
5403: }
5405: /*@C
5406: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
5408: Collective on TS
5410: Input Arguments:
5411: + ts - time stepping context
5412: . t - time at which to evaluate
5413: . U - state at which to evaluate
5414: - ctx - context
5416: Output Arguments:
5417: + A - pointer to operator
5418: . B - pointer to preconditioning matrix
5419: - flg - matrix structure flag
5421: Level: intermediate
5423: Notes:
5424: This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems.
5426: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear()
5427: @*/
5428: PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
5429: {
5431: return(0);
5432: }
5434: /*@C
5435: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
5437: Collective on TS
5439: Input Arguments:
5440: + ts - time stepping context
5441: . t - time at which to evaluate
5442: . U - state at which to evaluate
5443: . Udot - time derivative of state vector
5444: - ctx - context
5446: Output Arguments:
5447: . F - left hand side
5449: Level: intermediate
5451: Notes:
5452: 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
5453: user is required to write their own TSComputeIFunction.
5454: This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems.
5455: The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian().
5457: Note that using this function is NOT equivalent to using TSComputeRHSFunctionLinear() since that solves Udot = A U
5459: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant(), TSComputeRHSFunctionLinear()
5460: @*/
5461: PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
5462: {
5464: Mat A,B;
5467: TSGetIJacobian(ts,&A,&B,NULL,NULL);
5468: TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);
5469: MatMult(A,Udot,F);
5470: return(0);
5471: }
5473: /*@C
5474: TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE
5476: Collective on TS
5478: Input Arguments:
5479: + ts - time stepping context
5480: . t - time at which to evaluate
5481: . U - state at which to evaluate
5482: . Udot - time derivative of state vector
5483: . shift - shift to apply
5484: - ctx - context
5486: Output Arguments:
5487: + A - pointer to operator
5488: . B - pointer to preconditioning matrix
5489: - flg - matrix structure flag
5491: Level: advanced
5493: Notes:
5494: This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems.
5496: It is only appropriate for problems of the form
5498: $ M Udot = F(U,t)
5500: where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only
5501: works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing
5502: an implicit operator of the form
5504: $ shift*M + J
5506: 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
5507: a copy of M or reassemble it when requested.
5509: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear()
5510: @*/
5511: PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx)
5512: {
5516: MatScale(A, shift / ts->ijacobian.shift);
5517: ts->ijacobian.shift = shift;
5518: return(0);
5519: }
5521: /*@
5522: TSGetEquationType - Gets the type of the equation that TS is solving.
5524: Not Collective
5526: Input Parameter:
5527: . ts - the TS context
5529: Output Parameter:
5530: . equation_type - see TSEquationType
5532: Level: beginner
5534: .keywords: TS, equation type
5536: .seealso: TSSetEquationType(), TSEquationType
5537: @*/
5538: PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type)
5539: {
5543: *equation_type = ts->equation_type;
5544: return(0);
5545: }
5547: /*@
5548: TSSetEquationType - Sets the type of the equation that TS is solving.
5550: Not Collective
5552: Input Parameter:
5553: + ts - the TS context
5554: - equation_type - see TSEquationType
5556: Level: advanced
5558: .keywords: TS, equation type
5560: .seealso: TSGetEquationType(), TSEquationType
5561: @*/
5562: PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type)
5563: {
5566: ts->equation_type = equation_type;
5567: return(0);
5568: }
5570: /*@
5571: TSGetConvergedReason - Gets the reason the TS iteration was stopped.
5573: Not Collective
5575: Input Parameter:
5576: . ts - the TS context
5578: Output Parameter:
5579: . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
5580: manual pages for the individual convergence tests for complete lists
5582: Level: beginner
5584: Notes:
5585: Can only be called after the call to TSSolve() is complete.
5587: .keywords: TS, nonlinear, set, convergence, test
5589: .seealso: TSSetConvergenceTest(), TSConvergedReason
5590: @*/
5591: PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason)
5592: {
5596: *reason = ts->reason;
5597: return(0);
5598: }
5600: /*@
5601: TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve.
5603: Not Collective
5605: Input Parameter:
5606: + ts - the TS context
5607: . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
5608: manual pages for the individual convergence tests for complete lists
5610: Level: advanced
5612: Notes:
5613: Can only be called during TSSolve() is active.
5615: .keywords: TS, nonlinear, set, convergence, test
5617: .seealso: TSConvergedReason
5618: @*/
5619: PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason)
5620: {
5623: ts->reason = reason;
5624: return(0);
5625: }
5627: /*@
5628: TSGetSolveTime - Gets the time after a call to TSSolve()
5630: Not Collective
5632: Input Parameter:
5633: . ts - the TS context
5635: Output Parameter:
5636: . ftime - the final time. This time corresponds to the final time set with TSSetMaxTime()
5638: Level: beginner
5640: Notes:
5641: Can only be called after the call to TSSolve() is complete.
5643: .keywords: TS, nonlinear, set, convergence, test
5645: .seealso: TSSetConvergenceTest(), TSConvergedReason
5646: @*/
5647: PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime)
5648: {
5652: *ftime = ts->solvetime;
5653: return(0);
5654: }
5656: /*@
5657: TSGetSNESIterations - Gets the total number of nonlinear iterations
5658: used by the time integrator.
5660: Not Collective
5662: Input Parameter:
5663: . ts - TS context
5665: Output Parameter:
5666: . nits - number of nonlinear iterations
5668: Notes:
5669: This counter is reset to zero for each successive call to TSSolve().
5671: Level: intermediate
5673: .keywords: TS, get, number, nonlinear, iterations
5675: .seealso: TSGetKSPIterations()
5676: @*/
5677: PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits)
5678: {
5682: *nits = ts->snes_its;
5683: return(0);
5684: }
5686: /*@
5687: TSGetKSPIterations - Gets the total number of linear iterations
5688: used by the time integrator.
5690: Not Collective
5692: Input Parameter:
5693: . ts - TS context
5695: Output Parameter:
5696: . lits - number of linear iterations
5698: Notes:
5699: This counter is reset to zero for each successive call to TSSolve().
5701: Level: intermediate
5703: .keywords: TS, get, number, linear, iterations
5705: .seealso: TSGetSNESIterations(), SNESGetKSPIterations()
5706: @*/
5707: PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits)
5708: {
5712: *lits = ts->ksp_its;
5713: return(0);
5714: }
5716: /*@
5717: TSGetStepRejections - Gets the total number of rejected steps.
5719: Not Collective
5721: Input Parameter:
5722: . ts - TS context
5724: Output Parameter:
5725: . rejects - number of steps rejected
5727: Notes:
5728: This counter is reset to zero for each successive call to TSSolve().
5730: Level: intermediate
5732: .keywords: TS, get, number
5734: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails()
5735: @*/
5736: PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects)
5737: {
5741: *rejects = ts->reject;
5742: return(0);
5743: }
5745: /*@
5746: TSGetSNESFailures - Gets the total number of failed SNES solves
5748: Not Collective
5750: Input Parameter:
5751: . ts - TS context
5753: Output Parameter:
5754: . fails - number of failed nonlinear solves
5756: Notes:
5757: This counter is reset to zero for each successive call to TSSolve().
5759: Level: intermediate
5761: .keywords: TS, get, number
5763: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures()
5764: @*/
5765: PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails)
5766: {
5770: *fails = ts->num_snes_failures;
5771: return(0);
5772: }
5774: /*@
5775: TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails
5777: Not Collective
5779: Input Parameter:
5780: + ts - TS context
5781: - rejects - maximum number of rejected steps, pass -1 for unlimited
5783: Notes:
5784: The counter is reset to zero for each step
5786: Options Database Key:
5787: . -ts_max_reject - Maximum number of step rejections before a step fails
5789: Level: intermediate
5791: .keywords: TS, set, maximum, number
5793: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
5794: @*/
5795: PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects)
5796: {
5799: ts->max_reject = rejects;
5800: return(0);
5801: }
5803: /*@
5804: TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves
5806: Not Collective
5808: Input Parameter:
5809: + ts - TS context
5810: - fails - maximum number of failed nonlinear solves, pass -1 for unlimited
5812: Notes:
5813: The counter is reset to zero for each successive call to TSSolve().
5815: Options Database Key:
5816: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
5818: Level: intermediate
5820: .keywords: TS, set, maximum, number
5822: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason()
5823: @*/
5824: PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails)
5825: {
5828: ts->max_snes_failures = fails;
5829: return(0);
5830: }
5832: /*@
5833: TSSetErrorIfStepFails - Error if no step succeeds
5835: Not Collective
5837: Input Parameter:
5838: + ts - TS context
5839: - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure
5841: Options Database Key:
5842: . -ts_error_if_step_fails - Error if no step succeeds
5844: Level: intermediate
5846: .keywords: TS, set, error
5848: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
5849: @*/
5850: PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err)
5851: {
5854: ts->errorifstepfailed = err;
5855: return(0);
5856: }
5858: /*@C
5859: TSMonitorSolution - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file or a PetscDraw object
5861: Collective on TS
5863: Input Parameters:
5864: + ts - the TS context
5865: . step - current time-step
5866: . ptime - current time
5867: . u - current state
5868: - vf - viewer and its format
5870: Level: intermediate
5872: .keywords: TS, vector, monitor, view
5874: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5875: @*/
5876: PetscErrorCode TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscViewerAndFormat *vf)
5877: {
5881: PetscViewerPushFormat(vf->viewer,vf->format);
5882: VecView(u,vf->viewer);
5883: PetscViewerPopFormat(vf->viewer);
5884: return(0);
5885: }
5887: /*@C
5888: TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep.
5890: Collective on TS
5892: Input Parameters:
5893: + ts - the TS context
5894: . step - current time-step
5895: . ptime - current time
5896: . u - current state
5897: - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
5899: Level: intermediate
5901: Notes:
5902: 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.
5903: These are named according to the file name template.
5905: This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy().
5907: .keywords: TS, vector, monitor, view
5909: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5910: @*/
5911: PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate)
5912: {
5914: char filename[PETSC_MAX_PATH_LEN];
5915: PetscViewer viewer;
5918: if (step < 0) return(0); /* -1 indicates interpolated solution */
5919: PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);
5920: PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);
5921: VecView(u,viewer);
5922: PetscViewerDestroy(&viewer);
5923: return(0);
5924: }
5926: /*@C
5927: TSMonitorSolutionVTKDestroy - Destroy context for monitoring
5929: Collective on TS
5931: Input Parameters:
5932: . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
5934: Level: intermediate
5936: Note:
5937: This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK().
5939: .keywords: TS, vector, monitor, view
5941: .seealso: TSMonitorSet(), TSMonitorSolutionVTK()
5942: @*/
5943: PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate)
5944: {
5948: PetscFree(*(char**)filenametemplate);
5949: return(0);
5950: }
5952: /*@
5953: TSGetAdapt - Get the adaptive controller context for the current method
5955: Collective on TS if controller has not been created yet
5957: Input Arguments:
5958: . ts - time stepping context
5960: Output Arguments:
5961: . adapt - adaptive controller
5963: Level: intermediate
5965: .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose()
5966: @*/
5967: PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt)
5968: {
5974: if (!ts->adapt) {
5975: TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);
5976: PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);
5977: PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);
5978: }
5979: *adapt = ts->adapt;
5980: return(0);
5981: }
5983: /*@
5984: TSSetTolerances - Set tolerances for local truncation error when using adaptive controller
5986: Logically Collective
5988: Input Arguments:
5989: + ts - time integration context
5990: . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value
5991: . vatol - vector of absolute tolerances or NULL, used in preference to atol if present
5992: . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value
5993: - vrtol - vector of relative tolerances or NULL, used in preference to atol if present
5995: Options Database keys:
5996: + -ts_rtol <rtol> - relative tolerance for local truncation error
5997: - -ts_atol <atol> Absolute tolerance for local truncation error
5999: Notes:
6000: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
6001: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
6002: computed only for the differential or the algebraic part then this can be done using the vector of
6003: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
6004: differential part and infinity for the algebraic part, the LTE calculation will include only the
6005: differential variables.
6007: Level: beginner
6009: .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances()
6010: @*/
6011: PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol)
6012: {
6016: if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol;
6017: if (vatol) {
6018: PetscObjectReference((PetscObject)vatol);
6019: VecDestroy(&ts->vatol);
6020: ts->vatol = vatol;
6021: }
6022: if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol;
6023: if (vrtol) {
6024: PetscObjectReference((PetscObject)vrtol);
6025: VecDestroy(&ts->vrtol);
6026: ts->vrtol = vrtol;
6027: }
6028: return(0);
6029: }
6031: /*@
6032: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
6034: Logically Collective
6036: Input Arguments:
6037: . ts - time integration context
6039: Output Arguments:
6040: + atol - scalar absolute tolerances, NULL to ignore
6041: . vatol - vector of absolute tolerances, NULL to ignore
6042: . rtol - scalar relative tolerances, NULL to ignore
6043: - vrtol - vector of relative tolerances, NULL to ignore
6045: Level: beginner
6047: .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances()
6048: @*/
6049: PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol)
6050: {
6052: if (atol) *atol = ts->atol;
6053: if (vatol) *vatol = ts->vatol;
6054: if (rtol) *rtol = ts->rtol;
6055: if (vrtol) *vrtol = ts->vrtol;
6056: return(0);
6057: }
6059: /*@
6060: TSErrorWeightedNorm2 - compute a weighted 2-norm of the difference between two state vectors
6062: Collective on TS
6064: Input Arguments:
6065: + ts - time stepping context
6066: . U - state vector, usually ts->vec_sol
6067: - Y - state vector to be compared to U
6069: Output Arguments:
6070: . norm - weighted norm, a value of 1.0 means that the error matches the tolerances
6071: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
6072: . normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
6074: Level: developer
6076: .seealso: TSErrorWeightedNorm(), TSErrorWeightedNormInfinity()
6077: @*/
6078: PetscErrorCode TSErrorWeightedNorm2(TS ts,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6079: {
6080: PetscErrorCode ierr;
6081: PetscInt i,n,N,rstart;
6082: PetscInt n_loc,na_loc,nr_loc;
6083: PetscReal n_glb,na_glb,nr_glb;
6084: const PetscScalar *u,*y;
6085: PetscReal sum,suma,sumr,gsum,gsuma,gsumr,diff;
6086: PetscReal tol,tola,tolr;
6087: PetscReal err_loc[6],err_glb[6];
6099: if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector");
6101: VecGetSize(U,&N);
6102: VecGetLocalSize(U,&n);
6103: VecGetOwnershipRange(U,&rstart,NULL);
6104: VecGetArrayRead(U,&u);
6105: VecGetArrayRead(Y,&y);
6106: sum = 0.; n_loc = 0;
6107: suma = 0.; na_loc = 0;
6108: sumr = 0.; nr_loc = 0;
6109: if (ts->vatol && ts->vrtol) {
6110: const PetscScalar *atol,*rtol;
6111: VecGetArrayRead(ts->vatol,&atol);
6112: VecGetArrayRead(ts->vrtol,&rtol);
6113: for (i=0; i<n; i++) {
6114: diff = PetscAbsScalar(y[i] - u[i]);
6115: tola = PetscRealPart(atol[i]);
6116: if(tola>0.){
6117: suma += PetscSqr(diff/tola);
6118: na_loc++;
6119: }
6120: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6121: if(tolr>0.){
6122: sumr += PetscSqr(diff/tolr);
6123: nr_loc++;
6124: }
6125: tol=tola+tolr;
6126: if(tol>0.){
6127: sum += PetscSqr(diff/tol);
6128: n_loc++;
6129: }
6130: }
6131: VecRestoreArrayRead(ts->vatol,&atol);
6132: VecRestoreArrayRead(ts->vrtol,&rtol);
6133: } else if (ts->vatol) { /* vector atol, scalar rtol */
6134: const PetscScalar *atol;
6135: VecGetArrayRead(ts->vatol,&atol);
6136: for (i=0; i<n; i++) {
6137: diff = PetscAbsScalar(y[i] - u[i]);
6138: tola = PetscRealPart(atol[i]);
6139: if(tola>0.){
6140: suma += PetscSqr(diff/tola);
6141: na_loc++;
6142: }
6143: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6144: if(tolr>0.){
6145: sumr += PetscSqr(diff/tolr);
6146: nr_loc++;
6147: }
6148: tol=tola+tolr;
6149: if(tol>0.){
6150: sum += PetscSqr(diff/tol);
6151: n_loc++;
6152: }
6153: }
6154: VecRestoreArrayRead(ts->vatol,&atol);
6155: } else if (ts->vrtol) { /* scalar atol, vector rtol */
6156: const PetscScalar *rtol;
6157: VecGetArrayRead(ts->vrtol,&rtol);
6158: for (i=0; i<n; i++) {
6159: diff = PetscAbsScalar(y[i] - u[i]);
6160: tola = ts->atol;
6161: if(tola>0.){
6162: suma += PetscSqr(diff/tola);
6163: na_loc++;
6164: }
6165: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6166: if(tolr>0.){
6167: sumr += PetscSqr(diff/tolr);
6168: nr_loc++;
6169: }
6170: tol=tola+tolr;
6171: if(tol>0.){
6172: sum += PetscSqr(diff/tol);
6173: n_loc++;
6174: }
6175: }
6176: VecRestoreArrayRead(ts->vrtol,&rtol);
6177: } else { /* scalar atol, scalar rtol */
6178: for (i=0; i<n; i++) {
6179: diff = PetscAbsScalar(y[i] - u[i]);
6180: tola = ts->atol;
6181: if(tola>0.){
6182: suma += PetscSqr(diff/tola);
6183: na_loc++;
6184: }
6185: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6186: if(tolr>0.){
6187: sumr += PetscSqr(diff/tolr);
6188: nr_loc++;
6189: }
6190: tol=tola+tolr;
6191: if(tol>0.){
6192: sum += PetscSqr(diff/tol);
6193: n_loc++;
6194: }
6195: }
6196: }
6197: VecRestoreArrayRead(U,&u);
6198: VecRestoreArrayRead(Y,&y);
6200: err_loc[0] = sum;
6201: err_loc[1] = suma;
6202: err_loc[2] = sumr;
6203: err_loc[3] = (PetscReal)n_loc;
6204: err_loc[4] = (PetscReal)na_loc;
6205: err_loc[5] = (PetscReal)nr_loc;
6207: MPIU_Allreduce(err_loc,err_glb,6,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));
6209: gsum = err_glb[0];
6210: gsuma = err_glb[1];
6211: gsumr = err_glb[2];
6212: n_glb = err_glb[3];
6213: na_glb = err_glb[4];
6214: nr_glb = err_glb[5];
6216: *norm = 0.;
6217: if(n_glb>0. ){*norm = PetscSqrtReal(gsum / n_glb );}
6218: *norma = 0.;
6219: if(na_glb>0.){*norma = PetscSqrtReal(gsuma / na_glb);}
6220: *normr = 0.;
6221: if(nr_glb>0.){*normr = PetscSqrtReal(gsumr / nr_glb);}
6223: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
6224: if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
6225: if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
6226: return(0);
6227: }
6229: /*@
6230: TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between two state vectors
6232: Collective on TS
6234: Input Arguments:
6235: + ts - time stepping context
6236: . U - state vector, usually ts->vec_sol
6237: - Y - state vector to be compared to U
6239: Output Arguments:
6240: . norm - weighted norm, a value of 1.0 means that the error matches the tolerances
6241: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
6242: . normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
6244: Level: developer
6246: .seealso: TSErrorWeightedNorm(), TSErrorWeightedNorm2()
6247: @*/
6248: PetscErrorCode TSErrorWeightedNormInfinity(TS ts,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6249: {
6250: PetscErrorCode ierr;
6251: PetscInt i,n,N,rstart;
6252: const PetscScalar *u,*y;
6253: PetscReal max,gmax,maxa,gmaxa,maxr,gmaxr;
6254: PetscReal tol,tola,tolr,diff;
6255: PetscReal err_loc[3],err_glb[3];
6267: if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector");
6269: VecGetSize(U,&N);
6270: VecGetLocalSize(U,&n);
6271: VecGetOwnershipRange(U,&rstart,NULL);
6272: VecGetArrayRead(U,&u);
6273: VecGetArrayRead(Y,&y);
6275: max=0.;
6276: maxa=0.;
6277: maxr=0.;
6279: if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
6280: const PetscScalar *atol,*rtol;
6281: VecGetArrayRead(ts->vatol,&atol);
6282: VecGetArrayRead(ts->vrtol,&rtol);
6284: for (i=0; i<n; i++) {
6285: diff = PetscAbsScalar(y[i] - u[i]);
6286: tola = PetscRealPart(atol[i]);
6287: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6288: tol = tola+tolr;
6289: if(tola>0.){
6290: maxa = PetscMax(maxa,diff / tola);
6291: }
6292: if(tolr>0.){
6293: maxr = PetscMax(maxr,diff / tolr);
6294: }
6295: if(tol>0.){
6296: max = PetscMax(max,diff / tol);
6297: }
6298: }
6299: VecRestoreArrayRead(ts->vatol,&atol);
6300: VecRestoreArrayRead(ts->vrtol,&rtol);
6301: } else if (ts->vatol) { /* vector atol, scalar rtol */
6302: const PetscScalar *atol;
6303: VecGetArrayRead(ts->vatol,&atol);
6304: for (i=0; i<n; i++) {
6305: diff = PetscAbsScalar(y[i] - u[i]);
6306: tola = PetscRealPart(atol[i]);
6307: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6308: tol = tola+tolr;
6309: if(tola>0.){
6310: maxa = PetscMax(maxa,diff / tola);
6311: }
6312: if(tolr>0.){
6313: maxr = PetscMax(maxr,diff / tolr);
6314: }
6315: if(tol>0.){
6316: max = PetscMax(max,diff / tol);
6317: }
6318: }
6319: VecRestoreArrayRead(ts->vatol,&atol);
6320: } else if (ts->vrtol) { /* scalar atol, vector rtol */
6321: const PetscScalar *rtol;
6322: VecGetArrayRead(ts->vrtol,&rtol);
6324: for (i=0; i<n; i++) {
6325: diff = PetscAbsScalar(y[i] - u[i]);
6326: tola = ts->atol;
6327: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6328: tol = tola+tolr;
6329: if(tola>0.){
6330: maxa = PetscMax(maxa,diff / tola);
6331: }
6332: if(tolr>0.){
6333: maxr = PetscMax(maxr,diff / tolr);
6334: }
6335: if(tol>0.){
6336: max = PetscMax(max,diff / tol);
6337: }
6338: }
6339: VecRestoreArrayRead(ts->vrtol,&rtol);
6340: } else { /* scalar atol, scalar rtol */
6342: for (i=0; i<n; i++) {
6343: diff = PetscAbsScalar(y[i] - u[i]);
6344: tola = ts->atol;
6345: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6346: tol = tola+tolr;
6347: if(tola>0.){
6348: maxa = PetscMax(maxa,diff / tola);
6349: }
6350: if(tolr>0.){
6351: maxr = PetscMax(maxr,diff / tolr);
6352: }
6353: if(tol>0.){
6354: max = PetscMax(max,diff / tol);
6355: }
6356: }
6357: }
6358: VecRestoreArrayRead(U,&u);
6359: VecRestoreArrayRead(Y,&y);
6360: err_loc[0] = max;
6361: err_loc[1] = maxa;
6362: err_loc[2] = maxr;
6363: MPIU_Allreduce(err_loc,err_glb,3,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));
6364: gmax = err_glb[0];
6365: gmaxa = err_glb[1];
6366: gmaxr = err_glb[2];
6368: *norm = gmax;
6369: *norma = gmaxa;
6370: *normr = gmaxr;
6371: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
6372: if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
6373: if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
6374: return(0);
6375: }
6377: /*@
6378: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
6380: Collective on TS
6382: Input Arguments:
6383: + ts - time stepping context
6384: . U - state vector, usually ts->vec_sol
6385: . Y - state vector to be compared to U
6386: - wnormtype - norm type, either NORM_2 or NORM_INFINITY
6388: Output Arguments:
6389: . norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
6390: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
6391: . normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
6393: Options Database Keys:
6394: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
6396: Level: developer
6398: .seealso: TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2(), TSErrorWeightedENorm
6399: @*/
6400: PetscErrorCode TSErrorWeightedNorm(TS ts,Vec U,Vec Y,NormType wnormtype,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6401: {
6405: if (wnormtype == NORM_2) {
6406: TSErrorWeightedNorm2(ts,U,Y,norm,norma,normr);
6407: } else if(wnormtype == NORM_INFINITY) {
6408: TSErrorWeightedNormInfinity(ts,U,Y,norm,norma,normr);
6409: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
6410: return(0);
6411: }
6414: /*@
6415: TSErrorWeightedENorm2 - compute a weighted 2 error norm based on supplied absolute and relative tolerances
6417: Collective on TS
6419: Input Arguments:
6420: + ts - time stepping context
6421: . E - error vector
6422: . U - state vector, usually ts->vec_sol
6423: - Y - state vector, previous time step
6425: Output Arguments:
6426: . norm - weighted norm, a value of 1.0 means that the error matches the tolerances
6427: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
6428: . normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
6430: Level: developer
6432: .seealso: TSErrorWeightedENorm(), TSErrorWeightedENormInfinity()
6433: @*/
6434: PetscErrorCode TSErrorWeightedENorm2(TS ts,Vec E,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6435: {
6436: PetscErrorCode ierr;
6437: PetscInt i,n,N,rstart;
6438: PetscInt n_loc,na_loc,nr_loc;
6439: PetscReal n_glb,na_glb,nr_glb;
6440: const PetscScalar *e,*u,*y;
6441: PetscReal err,sum,suma,sumr,gsum,gsuma,gsumr;
6442: PetscReal tol,tola,tolr;
6443: PetscReal err_loc[6],err_glb[6];
6459: VecGetSize(E,&N);
6460: VecGetLocalSize(E,&n);
6461: VecGetOwnershipRange(E,&rstart,NULL);
6462: VecGetArrayRead(E,&e);
6463: VecGetArrayRead(U,&u);
6464: VecGetArrayRead(Y,&y);
6465: sum = 0.; n_loc = 0;
6466: suma = 0.; na_loc = 0;
6467: sumr = 0.; nr_loc = 0;
6468: if (ts->vatol && ts->vrtol) {
6469: const PetscScalar *atol,*rtol;
6470: VecGetArrayRead(ts->vatol,&atol);
6471: VecGetArrayRead(ts->vrtol,&rtol);
6472: for (i=0; i<n; i++) {
6473: err = PetscAbsScalar(e[i]);
6474: tola = PetscRealPart(atol[i]);
6475: if(tola>0.){
6476: suma += PetscSqr(err/tola);
6477: na_loc++;
6478: }
6479: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6480: if(tolr>0.){
6481: sumr += PetscSqr(err/tolr);
6482: nr_loc++;
6483: }
6484: tol=tola+tolr;
6485: if(tol>0.){
6486: sum += PetscSqr(err/tol);
6487: n_loc++;
6488: }
6489: }
6490: VecRestoreArrayRead(ts->vatol,&atol);
6491: VecRestoreArrayRead(ts->vrtol,&rtol);
6492: } else if (ts->vatol) { /* vector atol, scalar rtol */
6493: const PetscScalar *atol;
6494: VecGetArrayRead(ts->vatol,&atol);
6495: for (i=0; i<n; i++) {
6496: err = PetscAbsScalar(e[i]);
6497: tola = PetscRealPart(atol[i]);
6498: if(tola>0.){
6499: suma += PetscSqr(err/tola);
6500: na_loc++;
6501: }
6502: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6503: if(tolr>0.){
6504: sumr += PetscSqr(err/tolr);
6505: nr_loc++;
6506: }
6507: tol=tola+tolr;
6508: if(tol>0.){
6509: sum += PetscSqr(err/tol);
6510: n_loc++;
6511: }
6512: }
6513: VecRestoreArrayRead(ts->vatol,&atol);
6514: } else if (ts->vrtol) { /* scalar atol, vector rtol */
6515: const PetscScalar *rtol;
6516: VecGetArrayRead(ts->vrtol,&rtol);
6517: for (i=0; i<n; i++) {
6518: err = PetscAbsScalar(e[i]);
6519: tola = ts->atol;
6520: if(tola>0.){
6521: suma += PetscSqr(err/tola);
6522: na_loc++;
6523: }
6524: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6525: if(tolr>0.){
6526: sumr += PetscSqr(err/tolr);
6527: nr_loc++;
6528: }
6529: tol=tola+tolr;
6530: if(tol>0.){
6531: sum += PetscSqr(err/tol);
6532: n_loc++;
6533: }
6534: }
6535: VecRestoreArrayRead(ts->vrtol,&rtol);
6536: } else { /* scalar atol, scalar rtol */
6537: for (i=0; i<n; i++) {
6538: err = PetscAbsScalar(e[i]);
6539: tola = ts->atol;
6540: if(tola>0.){
6541: suma += PetscSqr(err/tola);
6542: na_loc++;
6543: }
6544: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6545: if(tolr>0.){
6546: sumr += PetscSqr(err/tolr);
6547: nr_loc++;
6548: }
6549: tol=tola+tolr;
6550: if(tol>0.){
6551: sum += PetscSqr(err/tol);
6552: n_loc++;
6553: }
6554: }
6555: }
6556: VecRestoreArrayRead(E,&e);
6557: VecRestoreArrayRead(U,&u);
6558: VecRestoreArrayRead(Y,&y);
6560: err_loc[0] = sum;
6561: err_loc[1] = suma;
6562: err_loc[2] = sumr;
6563: err_loc[3] = (PetscReal)n_loc;
6564: err_loc[4] = (PetscReal)na_loc;
6565: err_loc[5] = (PetscReal)nr_loc;
6567: MPIU_Allreduce(err_loc,err_glb,6,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));
6569: gsum = err_glb[0];
6570: gsuma = err_glb[1];
6571: gsumr = err_glb[2];
6572: n_glb = err_glb[3];
6573: na_glb = err_glb[4];
6574: nr_glb = err_glb[5];
6576: *norm = 0.;
6577: if(n_glb>0. ){*norm = PetscSqrtReal(gsum / n_glb );}
6578: *norma = 0.;
6579: if(na_glb>0.){*norma = PetscSqrtReal(gsuma / na_glb);}
6580: *normr = 0.;
6581: if(nr_glb>0.){*normr = PetscSqrtReal(gsumr / nr_glb);}
6583: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
6584: if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
6585: if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
6586: return(0);
6587: }
6589: /*@
6590: TSErrorWeightedENormInfinity - compute a weighted infinity error norm based on supplied absolute and relative tolerances
6591: Collective on TS
6593: Input Arguments:
6594: + ts - time stepping context
6595: . E - error vector
6596: . U - state vector, usually ts->vec_sol
6597: - Y - state vector, previous time step
6599: Output Arguments:
6600: . norm - weighted norm, a value of 1.0 means that the error matches the tolerances
6601: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
6602: . normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
6604: Level: developer
6606: .seealso: TSErrorWeightedENorm(), TSErrorWeightedENorm2()
6607: @*/
6608: PetscErrorCode TSErrorWeightedENormInfinity(TS ts,Vec E,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6609: {
6610: PetscErrorCode ierr;
6611: PetscInt i,n,N,rstart;
6612: const PetscScalar *e,*u,*y;
6613: PetscReal err,max,gmax,maxa,gmaxa,maxr,gmaxr;
6614: PetscReal tol,tola,tolr;
6615: PetscReal err_loc[3],err_glb[3];
6631: VecGetSize(E,&N);
6632: VecGetLocalSize(E,&n);
6633: VecGetOwnershipRange(E,&rstart,NULL);
6634: VecGetArrayRead(E,&e);
6635: VecGetArrayRead(U,&u);
6636: VecGetArrayRead(Y,&y);
6638: max=0.;
6639: maxa=0.;
6640: maxr=0.;
6642: if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
6643: const PetscScalar *atol,*rtol;
6644: VecGetArrayRead(ts->vatol,&atol);
6645: VecGetArrayRead(ts->vrtol,&rtol);
6647: for (i=0; i<n; i++) {
6648: err = PetscAbsScalar(e[i]);
6649: tola = PetscRealPart(atol[i]);
6650: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6651: tol = tola+tolr;
6652: if(tola>0.){
6653: maxa = PetscMax(maxa,err / tola);
6654: }
6655: if(tolr>0.){
6656: maxr = PetscMax(maxr,err / tolr);
6657: }
6658: if(tol>0.){
6659: max = PetscMax(max,err / tol);
6660: }
6661: }
6662: VecRestoreArrayRead(ts->vatol,&atol);
6663: VecRestoreArrayRead(ts->vrtol,&rtol);
6664: } else if (ts->vatol) { /* vector atol, scalar rtol */
6665: const PetscScalar *atol;
6666: VecGetArrayRead(ts->vatol,&atol);
6667: for (i=0; i<n; i++) {
6668: err = PetscAbsScalar(e[i]);
6669: tola = PetscRealPart(atol[i]);
6670: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6671: tol = tola+tolr;
6672: if(tola>0.){
6673: maxa = PetscMax(maxa,err / tola);
6674: }
6675: if(tolr>0.){
6676: maxr = PetscMax(maxr,err / tolr);
6677: }
6678: if(tol>0.){
6679: max = PetscMax(max,err / tol);
6680: }
6681: }
6682: VecRestoreArrayRead(ts->vatol,&atol);
6683: } else if (ts->vrtol) { /* scalar atol, vector rtol */
6684: const PetscScalar *rtol;
6685: VecGetArrayRead(ts->vrtol,&rtol);
6687: for (i=0; i<n; i++) {
6688: err = PetscAbsScalar(e[i]);
6689: tola = ts->atol;
6690: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6691: tol = tola+tolr;
6692: if(tola>0.){
6693: maxa = PetscMax(maxa,err / tola);
6694: }
6695: if(tolr>0.){
6696: maxr = PetscMax(maxr,err / tolr);
6697: }
6698: if(tol>0.){
6699: max = PetscMax(max,err / tol);
6700: }
6701: }
6702: VecRestoreArrayRead(ts->vrtol,&rtol);
6703: } else { /* scalar atol, scalar rtol */
6705: for (i=0; i<n; i++) {
6706: err = PetscAbsScalar(e[i]);
6707: tola = ts->atol;
6708: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6709: tol = tola+tolr;
6710: if(tola>0.){
6711: maxa = PetscMax(maxa,err / tola);
6712: }
6713: if(tolr>0.){
6714: maxr = PetscMax(maxr,err / tolr);
6715: }
6716: if(tol>0.){
6717: max = PetscMax(max,err / tol);
6718: }
6719: }
6720: }
6721: VecRestoreArrayRead(E,&e);
6722: VecRestoreArrayRead(U,&u);
6723: VecRestoreArrayRead(Y,&y);
6724: err_loc[0] = max;
6725: err_loc[1] = maxa;
6726: err_loc[2] = maxr;
6727: MPIU_Allreduce(err_loc,err_glb,3,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));
6728: gmax = err_glb[0];
6729: gmaxa = err_glb[1];
6730: gmaxr = err_glb[2];
6732: *norm = gmax;
6733: *norma = gmaxa;
6734: *normr = gmaxr;
6735: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
6736: if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
6737: if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
6738: return(0);
6739: }
6741: /*@
6742: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
6744: Collective on TS
6746: Input Arguments:
6747: + ts - time stepping context
6748: . E - error vector
6749: . U - state vector, usually ts->vec_sol
6750: . Y - state vector, previous time step
6751: - wnormtype - norm type, either NORM_2 or NORM_INFINITY
6753: Output Arguments:
6754: . norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
6755: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
6756: . normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
6758: Options Database Keys:
6759: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
6761: Level: developer
6763: .seealso: TSErrorWeightedENormInfinity(), TSErrorWeightedENorm2(), TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2()
6764: @*/
6765: PetscErrorCode TSErrorWeightedENorm(TS ts,Vec E,Vec U,Vec Y,NormType wnormtype,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6766: {
6770: if (wnormtype == NORM_2) {
6771: TSErrorWeightedENorm2(ts,E,U,Y,norm,norma,normr);
6772: } else if(wnormtype == NORM_INFINITY) {
6773: TSErrorWeightedENormInfinity(ts,E,U,Y,norm,norma,normr);
6774: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
6775: return(0);
6776: }
6779: /*@
6780: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
6782: Logically Collective on TS
6784: Input Arguments:
6785: + ts - time stepping context
6786: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
6788: Note:
6789: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
6791: Level: intermediate
6793: .seealso: TSGetCFLTime(), TSADAPTCFL
6794: @*/
6795: PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime)
6796: {
6799: ts->cfltime_local = cfltime;
6800: ts->cfltime = -1.;
6801: return(0);
6802: }
6804: /*@
6805: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
6807: Collective on TS
6809: Input Arguments:
6810: . ts - time stepping context
6812: Output Arguments:
6813: . cfltime - maximum stable time step for forward Euler
6815: Level: advanced
6817: .seealso: TSSetCFLTimeLocal()
6818: @*/
6819: PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime)
6820: {
6824: if (ts->cfltime < 0) {
6825: MPIU_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));
6826: }
6827: *cfltime = ts->cfltime;
6828: return(0);
6829: }
6831: /*@
6832: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
6834: Input Parameters:
6835: . ts - the TS context.
6836: . xl - lower bound.
6837: . xu - upper bound.
6839: Notes:
6840: If this routine is not called then the lower and upper bounds are set to
6841: PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().
6843: Level: advanced
6845: @*/
6846: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
6847: {
6849: SNES snes;
6852: TSGetSNES(ts,&snes);
6853: SNESVISetVariableBounds(snes,xl,xu);
6854: return(0);
6855: }
6857: #if defined(PETSC_HAVE_MATLAB_ENGINE)
6858: #include <mex.h>
6860: typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext;
6862: /*
6863: TSComputeFunction_Matlab - Calls the function that has been set with
6864: TSSetFunctionMatlab().
6866: Collective on TS
6868: Input Parameters:
6869: + snes - the TS context
6870: - u - input vector
6872: Output Parameter:
6873: . y - function vector, as set by TSSetFunction()
6875: Notes:
6876: TSComputeFunction() is typically used within nonlinear solvers
6877: implementations, so most users would not generally call this routine
6878: themselves.
6880: Level: developer
6882: .keywords: TS, nonlinear, compute, function
6884: .seealso: TSSetFunction(), TSGetFunction()
6885: */
6886: PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx)
6887: {
6888: PetscErrorCode ierr;
6889: TSMatlabContext *sctx = (TSMatlabContext*)ctx;
6890: int nlhs = 1,nrhs = 7;
6891: mxArray *plhs[1],*prhs[7];
6892: long long int lx = 0,lxdot = 0,ly = 0,ls = 0;
6902: PetscMemcpy(&ls,&snes,sizeof(snes));
6903: PetscMemcpy(&lx,&u,sizeof(u));
6904: PetscMemcpy(&lxdot,&udot,sizeof(udot));
6905: PetscMemcpy(&ly,&y,sizeof(u));
6907: prhs[0] = mxCreateDoubleScalar((double)ls);
6908: prhs[1] = mxCreateDoubleScalar(time);
6909: prhs[2] = mxCreateDoubleScalar((double)lx);
6910: prhs[3] = mxCreateDoubleScalar((double)lxdot);
6911: prhs[4] = mxCreateDoubleScalar((double)ly);
6912: prhs[5] = mxCreateString(sctx->funcname);
6913: prhs[6] = sctx->ctx;
6914: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");
6915: mxGetScalar(plhs[0]);
6916: mxDestroyArray(prhs[0]);
6917: mxDestroyArray(prhs[1]);
6918: mxDestroyArray(prhs[2]);
6919: mxDestroyArray(prhs[3]);
6920: mxDestroyArray(prhs[4]);
6921: mxDestroyArray(prhs[5]);
6922: mxDestroyArray(plhs[0]);
6923: return(0);
6924: }
6926: /*
6927: TSSetFunctionMatlab - Sets the function evaluation routine and function
6928: vector for use by the TS routines in solving ODEs
6929: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
6931: Logically Collective on TS
6933: Input Parameters:
6934: + ts - the TS context
6935: - func - function evaluation routine
6937: Calling sequence of func:
6938: $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx);
6940: Level: beginner
6942: .keywords: TS, nonlinear, set, function
6944: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
6945: */
6946: PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx)
6947: {
6948: PetscErrorCode ierr;
6949: TSMatlabContext *sctx;
6952: /* currently sctx is memory bleed */
6953: PetscNew(&sctx);
6954: PetscStrallocpy(func,&sctx->funcname);
6955: /*
6956: This should work, but it doesn't
6957: sctx->ctx = ctx;
6958: mexMakeArrayPersistent(sctx->ctx);
6959: */
6960: sctx->ctx = mxDuplicateArray(ctx);
6962: TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);
6963: return(0);
6964: }
6966: /*
6967: TSComputeJacobian_Matlab - Calls the function that has been set with
6968: TSSetJacobianMatlab().
6970: Collective on TS
6972: Input Parameters:
6973: + ts - the TS context
6974: . u - input vector
6975: . A, B - the matrices
6976: - ctx - user context
6978: Level: developer
6980: .keywords: TS, nonlinear, compute, function
6982: .seealso: TSSetFunction(), TSGetFunction()
6983: @*/
6984: PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat A,Mat B,void *ctx)
6985: {
6986: PetscErrorCode ierr;
6987: TSMatlabContext *sctx = (TSMatlabContext*)ctx;
6988: int nlhs = 2,nrhs = 9;
6989: mxArray *plhs[2],*prhs[9];
6990: long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0;
6996: /* call Matlab function in ctx with arguments u and y */
6998: PetscMemcpy(&ls,&ts,sizeof(ts));
6999: PetscMemcpy(&lx,&u,sizeof(u));
7000: PetscMemcpy(&lxdot,&udot,sizeof(u));
7001: PetscMemcpy(&lA,A,sizeof(u));
7002: PetscMemcpy(&lB,B,sizeof(u));
7004: prhs[0] = mxCreateDoubleScalar((double)ls);
7005: prhs[1] = mxCreateDoubleScalar((double)time);
7006: prhs[2] = mxCreateDoubleScalar((double)lx);
7007: prhs[3] = mxCreateDoubleScalar((double)lxdot);
7008: prhs[4] = mxCreateDoubleScalar((double)shift);
7009: prhs[5] = mxCreateDoubleScalar((double)lA);
7010: prhs[6] = mxCreateDoubleScalar((double)lB);
7011: prhs[7] = mxCreateString(sctx->funcname);
7012: prhs[8] = sctx->ctx;
7013: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");
7014: mxGetScalar(plhs[0]);
7015: mxDestroyArray(prhs[0]);
7016: mxDestroyArray(prhs[1]);
7017: mxDestroyArray(prhs[2]);
7018: mxDestroyArray(prhs[3]);
7019: mxDestroyArray(prhs[4]);
7020: mxDestroyArray(prhs[5]);
7021: mxDestroyArray(prhs[6]);
7022: mxDestroyArray(prhs[7]);
7023: mxDestroyArray(plhs[0]);
7024: mxDestroyArray(plhs[1]);
7025: return(0);
7026: }
7028: /*
7029: TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
7030: 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
7032: Logically Collective on TS
7034: Input Parameters:
7035: + ts - the TS context
7036: . A,B - Jacobian matrices
7037: . func - function evaluation routine
7038: - ctx - user context
7040: Calling sequence of func:
7041: $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx);
7043: Level: developer
7045: .keywords: TS, nonlinear, set, function
7047: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
7048: */
7049: PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx)
7050: {
7051: PetscErrorCode ierr;
7052: TSMatlabContext *sctx;
7055: /* currently sctx is memory bleed */
7056: PetscNew(&sctx);
7057: PetscStrallocpy(func,&sctx->funcname);
7058: /*
7059: This should work, but it doesn't
7060: sctx->ctx = ctx;
7061: mexMakeArrayPersistent(sctx->ctx);
7062: */
7063: sctx->ctx = mxDuplicateArray(ctx);
7065: TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);
7066: return(0);
7067: }
7069: /*
7070: TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab().
7072: Collective on TS
7074: .seealso: TSSetFunction(), TSGetFunction()
7075: @*/
7076: PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx)
7077: {
7078: PetscErrorCode ierr;
7079: TSMatlabContext *sctx = (TSMatlabContext*)ctx;
7080: int nlhs = 1,nrhs = 6;
7081: mxArray *plhs[1],*prhs[6];
7082: long long int lx = 0,ls = 0;
7088: PetscMemcpy(&ls,&ts,sizeof(ts));
7089: PetscMemcpy(&lx,&u,sizeof(u));
7091: prhs[0] = mxCreateDoubleScalar((double)ls);
7092: prhs[1] = mxCreateDoubleScalar((double)it);
7093: prhs[2] = mxCreateDoubleScalar((double)time);
7094: prhs[3] = mxCreateDoubleScalar((double)lx);
7095: prhs[4] = mxCreateString(sctx->funcname);
7096: prhs[5] = sctx->ctx;
7097: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");
7098: mxGetScalar(plhs[0]);
7099: mxDestroyArray(prhs[0]);
7100: mxDestroyArray(prhs[1]);
7101: mxDestroyArray(prhs[2]);
7102: mxDestroyArray(prhs[3]);
7103: mxDestroyArray(prhs[4]);
7104: mxDestroyArray(plhs[0]);
7105: return(0);
7106: }
7108: /*
7109: TSMonitorSetMatlab - Sets the monitor function from Matlab
7111: Level: developer
7113: .keywords: TS, nonlinear, set, function
7115: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
7116: */
7117: PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx)
7118: {
7119: PetscErrorCode ierr;
7120: TSMatlabContext *sctx;
7123: /* currently sctx is memory bleed */
7124: PetscNew(&sctx);
7125: PetscStrallocpy(func,&sctx->funcname);
7126: /*
7127: This should work, but it doesn't
7128: sctx->ctx = ctx;
7129: mexMakeArrayPersistent(sctx->ctx);
7130: */
7131: sctx->ctx = mxDuplicateArray(ctx);
7133: TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);
7134: return(0);
7135: }
7136: #endif
7138: /*@C
7139: TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector
7140: in a time based line graph
7142: Collective on TS
7144: Input Parameters:
7145: + ts - the TS context
7146: . step - current time-step
7147: . ptime - current time
7148: . u - current solution
7149: - dctx - the TSMonitorLGCtx object that contains all the options for the monitoring, this is created with TSMonitorLGCtxCreate()
7151: Options Database:
7152: . -ts_monitor_lg_solution_variables
7154: Level: intermediate
7156: Notes: Each process in a parallel run displays its component solutions in a separate window
7158: .keywords: TS, vector, monitor, view
7160: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
7161: TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
7162: TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
7163: TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()
7164: @*/
7165: PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
7166: {
7167: PetscErrorCode ierr;
7168: TSMonitorLGCtx ctx = (TSMonitorLGCtx)dctx;
7169: const PetscScalar *yy;
7170: Vec v;
7173: if (step < 0) return(0); /* -1 indicates interpolated solution */
7174: if (!step) {
7175: PetscDrawAxis axis;
7176: PetscInt dim;
7177: PetscDrawLGGetAxis(ctx->lg,&axis);
7178: PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");
7179: if (!ctx->names) {
7180: PetscBool flg;
7181: /* user provides names of variables to plot but no names has been set so assume names are integer values */
7182: PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",&flg);
7183: if (flg) {
7184: PetscInt i,n;
7185: char **names;
7186: VecGetSize(u,&n);
7187: PetscMalloc1(n+1,&names);
7188: for (i=0; i<n; i++) {
7189: PetscMalloc1(5,&names[i]);
7190: PetscSNPrintf(names[i],5,"%D",i);
7191: }
7192: names[n] = NULL;
7193: ctx->names = names;
7194: }
7195: }
7196: if (ctx->names && !ctx->displaynames) {
7197: char **displaynames;
7198: PetscBool flg;
7199: VecGetLocalSize(u,&dim);
7200: PetscMalloc1(dim+1,&displaynames);
7201: PetscMemzero(displaynames,(dim+1)*sizeof(char*));
7202: PetscOptionsGetStringArray(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",displaynames,&dim,&flg);
7203: if (flg) {
7204: TSMonitorLGCtxSetDisplayVariables(ctx,(const char *const *)displaynames);
7205: }
7206: PetscStrArrayDestroy(&displaynames);
7207: }
7208: if (ctx->displaynames) {
7209: PetscDrawLGSetDimension(ctx->lg,ctx->ndisplayvariables);
7210: PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->displaynames);
7211: } else if (ctx->names) {
7212: VecGetLocalSize(u,&dim);
7213: PetscDrawLGSetDimension(ctx->lg,dim);
7214: PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->names);
7215: } else {
7216: VecGetLocalSize(u,&dim);
7217: PetscDrawLGSetDimension(ctx->lg,dim);
7218: }
7219: PetscDrawLGReset(ctx->lg);
7220: }
7222: if (!ctx->transform) v = u;
7223: else {(*ctx->transform)(ctx->transformctx,u,&v);}
7224: VecGetArrayRead(v,&yy);
7225: if (ctx->displaynames) {
7226: PetscInt i;
7227: for (i=0; i<ctx->ndisplayvariables; i++)
7228: ctx->displayvalues[i] = PetscRealPart(yy[ctx->displayvariables[i]]);
7229: PetscDrawLGAddCommonPoint(ctx->lg,ptime,ctx->displayvalues);
7230: } else {
7231: #if defined(PETSC_USE_COMPLEX)
7232: PetscInt i,n;
7233: PetscReal *yreal;
7234: VecGetLocalSize(v,&n);
7235: PetscMalloc1(n,&yreal);
7236: for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
7237: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
7238: PetscFree(yreal);
7239: #else
7240: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
7241: #endif
7242: }
7243: VecRestoreArrayRead(v,&yy);
7244: if (ctx->transform) {VecDestroy(&v);}
7246: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
7247: PetscDrawLGDraw(ctx->lg);
7248: PetscDrawLGSave(ctx->lg);
7249: }
7250: return(0);
7251: }
7253: /*@C
7254: TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot
7256: Collective on TS
7258: Input Parameters:
7259: + ts - the TS context
7260: - names - the names of the components, final string must be NULL
7262: Level: intermediate
7264: Notes: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
7266: .keywords: TS, vector, monitor, view
7268: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetVariableNames()
7269: @*/
7270: PetscErrorCode TSMonitorLGSetVariableNames(TS ts,const char * const *names)
7271: {
7272: PetscErrorCode ierr;
7273: PetscInt i;
7276: for (i=0; i<ts->numbermonitors; i++) {
7277: if (ts->monitor[i] == TSMonitorLGSolution) {
7278: TSMonitorLGCtxSetVariableNames((TSMonitorLGCtx)ts->monitorcontext[i],names);
7279: break;
7280: }
7281: }
7282: return(0);
7283: }
7285: /*@C
7286: TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot
7288: Collective on TS
7290: Input Parameters:
7291: + ts - the TS context
7292: - names - the names of the components, final string must be NULL
7294: Level: intermediate
7296: .keywords: TS, vector, monitor, view
7298: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGSetVariableNames()
7299: @*/
7300: PetscErrorCode TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx,const char * const *names)
7301: {
7302: PetscErrorCode ierr;
7305: PetscStrArrayDestroy(&ctx->names);
7306: PetscStrArrayallocpy(names,&ctx->names);
7307: return(0);
7308: }
7310: /*@C
7311: TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot
7313: Collective on TS
7315: Input Parameter:
7316: . ts - the TS context
7318: Output Parameter:
7319: . names - the names of the components, final string must be NULL
7321: Level: intermediate
7323: Notes: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
7325: .keywords: TS, vector, monitor, view
7327: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
7328: @*/
7329: PetscErrorCode TSMonitorLGGetVariableNames(TS ts,const char *const **names)
7330: {
7331: PetscInt i;
7334: *names = NULL;
7335: for (i=0; i<ts->numbermonitors; i++) {
7336: if (ts->monitor[i] == TSMonitorLGSolution) {
7337: TSMonitorLGCtx ctx = (TSMonitorLGCtx) ts->monitorcontext[i];
7338: *names = (const char *const *)ctx->names;
7339: break;
7340: }
7341: }
7342: return(0);
7343: }
7345: /*@C
7346: TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor
7348: Collective on TS
7350: Input Parameters:
7351: + ctx - the TSMonitorLG context
7352: . displaynames - the names of the components, final string must be NULL
7354: Level: intermediate
7356: .keywords: TS, vector, monitor, view
7358: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
7359: @*/
7360: PetscErrorCode TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx,const char * const *displaynames)
7361: {
7362: PetscInt j = 0,k;
7363: PetscErrorCode ierr;
7366: if (!ctx->names) return(0);
7367: PetscStrArrayDestroy(&ctx->displaynames);
7368: PetscStrArrayallocpy(displaynames,&ctx->displaynames);
7369: while (displaynames[j]) j++;
7370: ctx->ndisplayvariables = j;
7371: PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvariables);
7372: PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvalues);
7373: j = 0;
7374: while (displaynames[j]) {
7375: k = 0;
7376: while (ctx->names[k]) {
7377: PetscBool flg;
7378: PetscStrcmp(displaynames[j],ctx->names[k],&flg);
7379: if (flg) {
7380: ctx->displayvariables[j] = k;
7381: break;
7382: }
7383: k++;
7384: }
7385: j++;
7386: }
7387: return(0);
7388: }
7390: /*@C
7391: TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor
7393: Collective on TS
7395: Input Parameters:
7396: + ts - the TS context
7397: . displaynames - the names of the components, final string must be NULL
7399: Notes: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
7401: Level: intermediate
7403: .keywords: TS, vector, monitor, view
7405: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
7406: @*/
7407: PetscErrorCode TSMonitorLGSetDisplayVariables(TS ts,const char * const *displaynames)
7408: {
7409: PetscInt i;
7410: PetscErrorCode ierr;
7413: for (i=0; i<ts->numbermonitors; i++) {
7414: if (ts->monitor[i] == TSMonitorLGSolution) {
7415: TSMonitorLGCtxSetDisplayVariables((TSMonitorLGCtx)ts->monitorcontext[i],displaynames);
7416: break;
7417: }
7418: }
7419: return(0);
7420: }
7422: /*@C
7423: TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed
7425: Collective on TS
7427: Input Parameters:
7428: + ts - the TS context
7429: . transform - the transform function
7430: . destroy - function to destroy the optional context
7431: - ctx - optional context used by transform function
7433: Notes: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
7435: Level: intermediate
7437: .keywords: TS, vector, monitor, view
7439: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGCtxSetTransform()
7440: @*/
7441: PetscErrorCode TSMonitorLGSetTransform(TS ts,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
7442: {
7443: PetscInt i;
7444: PetscErrorCode ierr;
7447: for (i=0; i<ts->numbermonitors; i++) {
7448: if (ts->monitor[i] == TSMonitorLGSolution) {
7449: TSMonitorLGCtxSetTransform((TSMonitorLGCtx)ts->monitorcontext[i],transform,destroy,tctx);
7450: }
7451: }
7452: return(0);
7453: }
7455: /*@C
7456: TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed
7458: Collective on TSLGCtx
7460: Input Parameters:
7461: + ts - the TS context
7462: . transform - the transform function
7463: . destroy - function to destroy the optional context
7464: - ctx - optional context used by transform function
7466: Level: intermediate
7468: .keywords: TS, vector, monitor, view
7470: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGSetTransform()
7471: @*/
7472: PetscErrorCode TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
7473: {
7475: ctx->transform = transform;
7476: ctx->transformdestroy = destroy;
7477: ctx->transformctx = tctx;
7478: return(0);
7479: }
7481: /*@C
7482: TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the error
7483: in a time based line graph
7485: Collective on TS
7487: Input Parameters:
7488: + ts - the TS context
7489: . step - current time-step
7490: . ptime - current time
7491: . u - current solution
7492: - dctx - TSMonitorLGCtx object created with TSMonitorLGCtxCreate()
7494: Level: intermediate
7496: Notes: Each process in a parallel run displays its component errors in a separate window
7498: The user must provide the solution using TSSetSolutionFunction() to use this monitor.
7500: Options Database Keys:
7501: . -ts_monitor_lg_error - create a graphical monitor of error history
7503: .keywords: TS, vector, monitor, view
7505: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
7506: @*/
7507: PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
7508: {
7509: PetscErrorCode ierr;
7510: TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy;
7511: const PetscScalar *yy;
7512: Vec y;
7515: if (!step) {
7516: PetscDrawAxis axis;
7517: PetscInt dim;
7518: PetscDrawLGGetAxis(ctx->lg,&axis);
7519: PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Error");
7520: VecGetLocalSize(u,&dim);
7521: PetscDrawLGSetDimension(ctx->lg,dim);
7522: PetscDrawLGReset(ctx->lg);
7523: }
7524: VecDuplicate(u,&y);
7525: TSComputeSolutionFunction(ts,ptime,y);
7526: VecAXPY(y,-1.0,u);
7527: VecGetArrayRead(y,&yy);
7528: #if defined(PETSC_USE_COMPLEX)
7529: {
7530: PetscReal *yreal;
7531: PetscInt i,n;
7532: VecGetLocalSize(y,&n);
7533: PetscMalloc1(n,&yreal);
7534: for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
7535: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
7536: PetscFree(yreal);
7537: }
7538: #else
7539: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
7540: #endif
7541: VecRestoreArrayRead(y,&yy);
7542: VecDestroy(&y);
7543: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
7544: PetscDrawLGDraw(ctx->lg);
7545: PetscDrawLGSave(ctx->lg);
7546: }
7547: return(0);
7548: }
7550: PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
7551: {
7552: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
7553: PetscReal x = ptime,y;
7555: PetscInt its;
7558: if (n < 0) return(0); /* -1 indicates interpolated solution */
7559: if (!n) {
7560: PetscDrawAxis axis;
7561: PetscDrawLGGetAxis(ctx->lg,&axis);
7562: PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");
7563: PetscDrawLGReset(ctx->lg);
7564: ctx->snes_its = 0;
7565: }
7566: TSGetSNESIterations(ts,&its);
7567: y = its - ctx->snes_its;
7568: PetscDrawLGAddPoint(ctx->lg,&x,&y);
7569: if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
7570: PetscDrawLGDraw(ctx->lg);
7571: PetscDrawLGSave(ctx->lg);
7572: }
7573: ctx->snes_its = its;
7574: return(0);
7575: }
7577: PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
7578: {
7579: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
7580: PetscReal x = ptime,y;
7582: PetscInt its;
7585: if (n < 0) return(0); /* -1 indicates interpolated solution */
7586: if (!n) {
7587: PetscDrawAxis axis;
7588: PetscDrawLGGetAxis(ctx->lg,&axis);
7589: PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");
7590: PetscDrawLGReset(ctx->lg);
7591: ctx->ksp_its = 0;
7592: }
7593: TSGetKSPIterations(ts,&its);
7594: y = its - ctx->ksp_its;
7595: PetscDrawLGAddPoint(ctx->lg,&x,&y);
7596: if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
7597: PetscDrawLGDraw(ctx->lg);
7598: PetscDrawLGSave(ctx->lg);
7599: }
7600: ctx->ksp_its = its;
7601: return(0);
7602: }
7604: /*@
7605: TSComputeLinearStability - computes the linear stability function at a point
7607: Collective on TS and Vec
7609: Input Parameters:
7610: + ts - the TS context
7611: - xr,xi - real and imaginary part of input arguments
7613: Output Parameters:
7614: . yr,yi - real and imaginary part of function value
7616: Level: developer
7618: .keywords: TS, compute
7620: .seealso: TSSetRHSFunction(), TSComputeIFunction()
7621: @*/
7622: PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
7623: {
7628: if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method");
7629: (*ts->ops->linearstability)(ts,xr,xi,yr,yi);
7630: return(0);
7631: }
7633: /* ------------------------------------------------------------------------*/
7634: /*@C
7635: TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope()
7637: Collective on TS
7639: Input Parameters:
7640: . ts - the ODE solver object
7642: Output Parameter:
7643: . ctx - the context
7645: Level: intermediate
7647: .keywords: TS, monitor, line graph, residual, seealso
7649: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()
7651: @*/
7652: PetscErrorCode TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx)
7653: {
7657: PetscNew(ctx);
7658: return(0);
7659: }
7661: /*@C
7662: TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution
7664: Collective on TS
7666: Input Parameters:
7667: + ts - the TS context
7668: . step - current time-step
7669: . ptime - current time
7670: . u - current solution
7671: - dctx - the envelope context
7673: Options Database:
7674: . -ts_monitor_envelope
7676: Level: intermediate
7678: Notes: after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope
7680: .keywords: TS, vector, monitor, view
7682: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxCreate()
7683: @*/
7684: PetscErrorCode TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
7685: {
7686: PetscErrorCode ierr;
7687: TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dctx;
7690: if (!ctx->max) {
7691: VecDuplicate(u,&ctx->max);
7692: VecDuplicate(u,&ctx->min);
7693: VecCopy(u,ctx->max);
7694: VecCopy(u,ctx->min);
7695: } else {
7696: VecPointwiseMax(ctx->max,u,ctx->max);
7697: VecPointwiseMin(ctx->min,u,ctx->min);
7698: }
7699: return(0);
7700: }
7702: /*@C
7703: TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution
7705: Collective on TS
7707: Input Parameter:
7708: . ts - the TS context
7710: Output Parameter:
7711: + max - the maximum values
7712: - min - the minimum values
7714: Notes: If the TS does not have a TSMonitorEnvelopeCtx associated with it then this function is ignored
7716: Level: intermediate
7718: .keywords: TS, vector, monitor, view
7720: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
7721: @*/
7722: PetscErrorCode TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min)
7723: {
7724: PetscInt i;
7727: if (max) *max = NULL;
7728: if (min) *min = NULL;
7729: for (i=0; i<ts->numbermonitors; i++) {
7730: if (ts->monitor[i] == TSMonitorEnvelope) {
7731: TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx) ts->monitorcontext[i];
7732: if (max) *max = ctx->max;
7733: if (min) *min = ctx->min;
7734: break;
7735: }
7736: }
7737: return(0);
7738: }
7740: /*@C
7741: TSMonitorEnvelopeCtxDestroy - Destroys a context that was created with TSMonitorEnvelopeCtxCreate().
7743: Collective on TSMonitorEnvelopeCtx
7745: Input Parameter:
7746: . ctx - the monitor context
7748: Level: intermediate
7750: .keywords: TS, monitor, line graph, destroy
7752: .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep()
7753: @*/
7754: PetscErrorCode TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx)
7755: {
7759: VecDestroy(&(*ctx)->min);
7760: VecDestroy(&(*ctx)->max);
7761: PetscFree(*ctx);
7762: return(0);
7763: }
7765: /*@
7766: TSRestartStep - Flags the solver to restart the next step
7768: Collective on TS
7770: Input Parameter:
7771: . ts - the TS context obtained from TSCreate()
7773: Level: advanced
7775: Notes:
7776: Multistep methods like BDF or Runge-Kutta methods with FSAL property require restarting the solver in the event of
7777: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
7778: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
7779: the sake of correctness and maximum safety, users are expected to call TSRestart() whenever they introduce
7780: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
7781: discontinuous source terms).
7783: .keywords: TS, timestep, restart
7785: .seealso: TSSolve(), TSSetPreStep(), TSSetPostStep()
7786: @*/
7787: PetscErrorCode TSRestartStep(TS ts)
7788: {
7791: ts->steprestart = PETSC_TRUE;
7792: return(0);
7793: }
7795: /*@
7796: TSRollBack - Rolls back one time step
7798: Collective on TS
7800: Input Parameter:
7801: . ts - the TS context obtained from TSCreate()
7803: Level: advanced
7805: .keywords: TS, timestep, rollback
7807: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate()
7808: @*/
7809: PetscErrorCode TSRollBack(TS ts)
7810: {
7815: if (ts->steprollback) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TSRollBack already called");
7816: if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name);
7817: (*ts->ops->rollback)(ts);
7818: ts->time_step = ts->ptime - ts->ptime_prev;
7819: ts->ptime = ts->ptime_prev;
7820: ts->ptime_prev = ts->ptime_prev_rollback;
7821: ts->steps--;
7822: ts->steprollback = PETSC_TRUE;
7823: return(0);
7824: }
7826: /*@
7827: TSGetStages - Get the number of stages and stage values
7829: Input Parameter:
7830: . ts - the TS context obtained from TSCreate()
7832: Level: advanced
7834: .keywords: TS, getstages
7836: .seealso: TSCreate()
7837: @*/
7838: PetscErrorCode TSGetStages(TS ts,PetscInt *ns,Vec **Y)
7839: {
7846: if (!ts->ops->getstages) *ns=0;
7847: else {
7848: (*ts->ops->getstages)(ts,ns,Y);
7849: }
7850: return(0);
7851: }
7853: /*@C
7854: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
7856: Collective on SNES
7858: Input Parameters:
7859: + ts - the TS context
7860: . t - current timestep
7861: . U - state vector
7862: . Udot - time derivative of state vector
7863: . shift - shift to apply, see note below
7864: - ctx - an optional user context
7866: Output Parameters:
7867: + J - Jacobian matrix (not altered in this routine)
7868: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)
7870: Level: intermediate
7872: Notes:
7873: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
7875: dF/dU + shift*dF/dUdot
7877: Most users should not need to explicitly call this routine, as it
7878: is used internally within the nonlinear solvers.
7880: This will first try to get the coloring from the DM. If the DM type has no coloring
7881: routine, then it will try to get the coloring from the matrix. This requires that the
7882: matrix have nonzero entries precomputed.
7884: .keywords: TS, finite differences, Jacobian, coloring, sparse
7885: .seealso: TSSetIJacobian(), MatFDColoringCreate(), MatFDColoringSetFunction()
7886: @*/
7887: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat J,Mat B,void *ctx)
7888: {
7889: SNES snes;
7890: MatFDColoring color;
7891: PetscBool hascolor, matcolor = PETSC_FALSE;
7895: PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject) ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL);
7896: PetscObjectQuery((PetscObject) B, "TSMatFDColoring", (PetscObject *) &color);
7897: if (!color) {
7898: DM dm;
7899: ISColoring iscoloring;
7901: TSGetDM(ts, &dm);
7902: DMHasColoring(dm, &hascolor);
7903: if (hascolor && !matcolor) {
7904: DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);
7905: MatFDColoringCreate(B, iscoloring, &color);
7906: MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);
7907: MatFDColoringSetFromOptions(color);
7908: MatFDColoringSetUp(B, iscoloring, color);
7909: ISColoringDestroy(&iscoloring);
7910: } else {
7911: MatColoring mc;
7913: MatColoringCreate(B, &mc);
7914: MatColoringSetDistance(mc, 2);
7915: MatColoringSetType(mc, MATCOLORINGSL);
7916: MatColoringSetFromOptions(mc);
7917: MatColoringApply(mc, &iscoloring);
7918: MatColoringDestroy(&mc);
7919: MatFDColoringCreate(B, iscoloring, &color);
7920: MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);
7921: MatFDColoringSetFromOptions(color);
7922: MatFDColoringSetUp(B, iscoloring, color);
7923: ISColoringDestroy(&iscoloring);
7924: }
7925: PetscObjectCompose((PetscObject) B, "TSMatFDColoring", (PetscObject) color);
7926: PetscObjectDereference((PetscObject) color);
7927: }
7928: TSGetSNES(ts, &snes);
7929: MatFDColoringApply(B, color, U, snes);
7930: if (J != B) {
7931: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
7932: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
7933: }
7934: return(0);
7935: }
7937: /*@
7938: TSSetFunctionDomainError - Set the function testing if the current state vector is valid
7940: Input Parameters:
7941: ts - the TS context
7942: func - function called within TSFunctionDomainError
7944: Level: intermediate
7946: .keywords: TS, state, domain
7947: .seealso: TSAdaptCheckStage(), TSFunctionDomainError()
7948: @*/
7950: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS,PetscReal,Vec,PetscBool*))
7951: {
7954: ts->functiondomainerror = func;
7955: return(0);
7956: }
7958: /*@
7959: TSFunctionDomainError - Check if the current state is valid
7961: Input Parameters:
7962: ts - the TS context
7963: stagetime - time of the simulation
7964: Y - state vector to check.
7966: Output Parameter:
7967: accept - Set to PETSC_FALSE if the current state vector is valid.
7969: Note:
7970: This function should be used to ensure the state is in a valid part of the space.
7971: For example, one can ensure here all values are positive.
7973: Level: advanced
7974: @*/
7975: PetscErrorCode TSFunctionDomainError(TS ts,PetscReal stagetime,Vec Y,PetscBool* accept)
7976: {
7982: *accept = PETSC_TRUE;
7983: if (ts->functiondomainerror) {
7984: PetscStackCallStandard((*ts->functiondomainerror),(ts,stagetime,Y,accept));
7985: }
7986: return(0);
7987: }
7989: /*@C
7990: TSClone - This function clones a time step object.
7992: Collective on MPI_Comm
7994: Input Parameter:
7995: . tsin - The input TS
7997: Output Parameter:
7998: . tsout - The output TS (cloned)
8000: Notes:
8001: This function is used to create a clone of a TS object. It is used in ARKIMEX for initializing the slope for first stage explicit methods. It will likely be replaced in the future with a mechanism of switching methods on the fly.
8003: When using TSDestroy() on a clone the user has to first reset the correct TS reference in the embedded SNES object: e.g.: by running SNES snes_dup=NULL; TSGetSNES(ts,&snes_dup); TSSetSNES(ts,snes_dup);
8005: Level: developer
8007: .keywords: TS, clone
8008: .seealso: TSCreate(), TSSetType(), TSSetUp(), TSDestroy(), TSSetProblemType()
8009: @*/
8010: PetscErrorCode TSClone(TS tsin, TS *tsout)
8011: {
8012: TS t;
8014: SNES snes_start;
8015: DM dm;
8016: TSType type;
8020: *tsout = NULL;
8022: PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView);
8024: /* General TS description */
8025: t->numbermonitors = 0;
8026: t->setupcalled = 0;
8027: t->ksp_its = 0;
8028: t->snes_its = 0;
8029: t->nwork = 0;
8030: t->rhsjacobian.time = -1e20;
8031: t->rhsjacobian.scale = 1.;
8032: t->ijacobian.shift = 1.;
8034: TSGetSNES(tsin,&snes_start);
8035: TSSetSNES(t,snes_start);
8037: TSGetDM(tsin,&dm);
8038: TSSetDM(t,dm);
8040: t->adapt = tsin->adapt;
8041: PetscObjectReference((PetscObject)t->adapt);
8043: t->trajectory = tsin->trajectory;
8044: PetscObjectReference((PetscObject)t->trajectory);
8046: t->event = tsin->event;
8047: if (t->event) t->event->refct++;
8049: t->problem_type = tsin->problem_type;
8050: t->ptime = tsin->ptime;
8051: t->ptime_prev = tsin->ptime_prev;
8052: t->time_step = tsin->time_step;
8053: t->max_time = tsin->max_time;
8054: t->steps = tsin->steps;
8055: t->max_steps = tsin->max_steps;
8056: t->equation_type = tsin->equation_type;
8057: t->atol = tsin->atol;
8058: t->rtol = tsin->rtol;
8059: t->max_snes_failures = tsin->max_snes_failures;
8060: t->max_reject = tsin->max_reject;
8061: t->errorifstepfailed = tsin->errorifstepfailed;
8063: TSGetType(tsin,&type);
8064: TSSetType(t,type);
8066: t->vec_sol = NULL;
8068: t->cfltime = tsin->cfltime;
8069: t->cfltime_local = tsin->cfltime_local;
8070: t->exact_final_time = tsin->exact_final_time;
8072: PetscMemcpy(t->ops,tsin->ops,sizeof(struct _TSOps));
8074: if (((PetscObject)tsin)->fortran_func_pointers) {
8075: PetscInt i;
8076: PetscMalloc((10)*sizeof(void(*)(void)),&((PetscObject)t)->fortran_func_pointers);
8077: for (i=0; i<10; i++) {
8078: ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
8079: }
8080: }
8081: *tsout = t;
8082: return(0);
8083: }