Actual source code: ts.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/tsimpl.h>
2: #include <petscdmshell.h>
3: #include <petscdmda.h>
4: #include <petscviewer.h>
5: #include <petscdraw.h>
6: #include <petscconvest.h>
8: #define SkipSmallValue(a,b,tol) if(PetscAbsScalar(a)< tol || PetscAbsScalar(b)< tol) continue;
10: /* Logging support */
11: PetscClassId TS_CLASSID, DMTS_CLASSID;
12: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;
14: const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED","STEPOVER","INTERPOLATE","MATCHSTEP","TSExactFinalTimeOption","TS_EXACTFINALTIME_",0};
17: /*@C
18: TSMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
20: Collective on TS
22: Input Parameters:
23: + ts - TS object you wish to monitor
24: . name - the monitor type one is seeking
25: . help - message indicating what monitoring is done
26: . manual - manual page for the monitor
27: . monitor - the monitor function
28: - 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
30: Level: developer
32: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
33: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
34: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
35: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
36: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
37: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
38: PetscOptionsFList(), PetscOptionsEList()
39: @*/
40: PetscErrorCode TSMonitorSetFromOptions(TS ts,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(TS,PetscViewerAndFormat*))
41: {
42: PetscErrorCode ierr;
43: PetscViewer viewer;
44: PetscViewerFormat format;
45: PetscBool flg;
48: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject) ts)->options,((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
49: if (flg) {
50: PetscViewerAndFormat *vf;
51: PetscViewerAndFormatCreate(viewer,format,&vf);
52: PetscObjectDereference((PetscObject)viewer);
53: if (monitorsetup) {
54: (*monitorsetup)(ts,vf);
55: }
56: TSMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
57: }
58: return(0);
59: }
61: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt,TSAdaptType default_type)
62: {
68: if (!((PetscObject)adapt)->type_name) {
69: TSAdaptSetType(adapt,default_type);
70: }
71: return(0);
72: }
74: /*@
75: TSSetFromOptions - Sets various TS parameters from user options.
77: Collective on TS
79: Input Parameter:
80: . ts - the TS context obtained from TSCreate()
82: Options Database Keys:
83: + -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSALPHA, TSGLLE, TSSSP, TSGLEE, TSBSYMP
84: . -ts_save_trajectory - checkpoint the solution at each time-step
85: . -ts_max_time <time> - maximum time to compute to
86: . -ts_max_steps <steps> - maximum number of time-steps to take
87: . -ts_init_time <time> - initial time to start computation
88: . -ts_final_time <time> - final time to compute to (deprecated: use -ts_max_time)
89: . -ts_dt <dt> - initial time step
90: . -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
91: . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
92: . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
93: . -ts_error_if_step_fails <true,false> - Error if no step succeeds
94: . -ts_rtol <rtol> - relative tolerance for local truncation error
95: . -ts_atol <atol> Absolute tolerance for local truncation error
96: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
97: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - test the Jacobian at each iteration against finite difference with RHS function
98: . -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
99: . -ts_fd_color - Use finite differences with coloring to compute IJacobian
100: . -ts_monitor - print information at each timestep
101: . -ts_monitor_lg_solution - Monitor solution graphically
102: . -ts_monitor_lg_error - Monitor error graphically
103: . -ts_monitor_error - Monitors norm of error
104: . -ts_monitor_lg_timestep - Monitor timestep size graphically
105: . -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
106: . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
107: . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
108: . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
109: . -ts_monitor_draw_solution - Monitor solution graphically
110: . -ts_monitor_draw_solution_phase <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
111: . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
112: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
113: . -ts_monitor_solution_vtk <filename.vts,filename.vtu> - Save each time step to a binary file, use filename-%%03D.vts (filename-%%03D.vtu)
114: - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
116: Developer Note:
117: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
119: Level: beginner
121: .seealso: TSGetType()
122: @*/
123: PetscErrorCode TSSetFromOptions(TS ts)
124: {
125: PetscBool opt,flg,tflg;
126: PetscErrorCode ierr;
127: char monfilename[PETSC_MAX_PATH_LEN];
128: PetscReal time_step;
129: TSExactFinalTimeOption eftopt;
130: char dir[16];
131: TSIFunction ifun;
132: const char *defaultType;
133: char typeName[256];
138: TSRegisterAll();
139: TSGetIFunction(ts,NULL,&ifun,NULL);
141: PetscObjectOptionsBegin((PetscObject)ts);
142: if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
143: else defaultType = ifun ? TSBEULER : TSEULER;
144: PetscOptionsFList("-ts_type","TS method","TSSetType",TSList,defaultType,typeName,256,&opt);
145: if (opt) {
146: TSSetType(ts,typeName);
147: } else {
148: TSSetType(ts,defaultType);
149: }
151: /* Handle generic TS options */
152: PetscOptionsDeprecated("-ts_final_time","-ts_max_time","3.10",NULL);
153: PetscOptionsReal("-ts_max_time","Maximum time to run to","TSSetMaxTime",ts->max_time,&ts->max_time,NULL);
154: PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetMaxSteps",ts->max_steps,&ts->max_steps,NULL);
155: PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,NULL);
156: PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&time_step,&flg);
157: if (flg) {TSSetTimeStep(ts,time_step);}
158: PetscOptionsEnum("-ts_exact_final_time","Option for handling of final time step","TSSetExactFinalTime",TSExactFinalTimeOptions,(PetscEnum)ts->exact_final_time,(PetscEnum*)&eftopt,&flg);
159: if (flg) {TSSetExactFinalTime(ts,eftopt);}
160: PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,NULL);
161: PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,NULL);
162: PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",ts->errorifstepfailed,&ts->errorifstepfailed,NULL);
163: PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,NULL);
164: PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,NULL);
166: PetscOptionsBool("-ts_rhs_jacobian_test_mult","Test the RHS Jacobian for consistency with RHS at each solve ","None",ts->testjacobian,&ts->testjacobian,NULL);
167: PetscOptionsBool("-ts_rhs_jacobian_test_mult_transpose","Test the RHS Jacobian transpose for consistency with RHS at each solve ","None",ts->testjacobiantranspose,&ts->testjacobiantranspose,NULL);
168: PetscOptionsBool("-ts_use_splitrhsfunction","Use the split RHS function for multirate solvers ","TSSetUseSplitRHSFunction",ts->use_splitrhsfunction,&ts->use_splitrhsfunction,NULL);
169: #if defined(PETSC_HAVE_SAWS)
170: {
171: PetscBool set;
172: flg = PETSC_FALSE;
173: PetscOptionsBool("-ts_saws_block","Block for SAWs memory snooper at end of TSSolve","PetscObjectSAWsBlock",((PetscObject)ts)->amspublishblock,&flg,&set);
174: if (set) {
175: PetscObjectSAWsSetBlock((PetscObject)ts,flg);
176: }
177: }
178: #endif
180: /* Monitor options */
181: TSMonitorSetFromOptions(ts,"-ts_monitor","Monitor time and timestep size","TSMonitorDefault",TSMonitorDefault,NULL);
182: TSMonitorSetFromOptions(ts,"-ts_monitor_extreme","Monitor extreme values of the solution","TSMonitorExtreme",TSMonitorExtreme,NULL);
183: TSMonitorSetFromOptions(ts,"-ts_monitor_solution","View the solution at each timestep","TSMonitorSolution",TSMonitorSolution,NULL);
185: PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
186: if (flg) {PetscPythonMonitorSet((PetscObject)ts,monfilename);}
188: PetscOptionsName("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",&opt);
189: if (opt) {
190: TSMonitorLGCtx ctx;
191: PetscInt howoften = 1;
193: PetscOptionsInt("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",howoften,&howoften,NULL);
194: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
195: TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
196: }
198: PetscOptionsName("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",&opt);
199: if (opt) {
200: TSMonitorLGCtx ctx;
201: PetscInt howoften = 1;
203: PetscOptionsInt("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",howoften,&howoften,NULL);
204: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
205: TSMonitorSet(ts,TSMonitorLGError,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
206: }
207: TSMonitorSetFromOptions(ts,"-ts_monitor_error","View the error at each timestep","TSMonitorError",TSMonitorError,NULL);
209: PetscOptionsName("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",&opt);
210: if (opt) {
211: TSMonitorLGCtx ctx;
212: PetscInt howoften = 1;
214: PetscOptionsInt("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);
215: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
216: TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
217: }
218: PetscOptionsName("-ts_monitor_lg_timestep_log","Monitor log timestep size graphically","TSMonitorLGTimeStep",&opt);
219: if (opt) {
220: TSMonitorLGCtx ctx;
221: PetscInt howoften = 1;
223: PetscOptionsInt("-ts_monitor_lg_timestep_log","Monitor log timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);
224: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
225: TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
226: ctx->semilogy = PETSC_TRUE;
227: }
229: PetscOptionsName("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",&opt);
230: if (opt) {
231: TSMonitorLGCtx ctx;
232: PetscInt howoften = 1;
234: PetscOptionsInt("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",howoften,&howoften,NULL);
235: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
236: TSMonitorSet(ts,TSMonitorLGSNESIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
237: }
238: PetscOptionsName("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",&opt);
239: if (opt) {
240: TSMonitorLGCtx ctx;
241: PetscInt howoften = 1;
243: PetscOptionsInt("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",howoften,&howoften,NULL);
244: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
245: TSMonitorSet(ts,TSMonitorLGKSPIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
246: }
247: PetscOptionsName("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",&opt);
248: if (opt) {
249: TSMonitorSPEigCtx ctx;
250: PetscInt howoften = 1;
252: PetscOptionsInt("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",howoften,&howoften,NULL);
253: TSMonitorSPEigCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
254: TSMonitorSet(ts,TSMonitorSPEig,ctx,(PetscErrorCode (*)(void**))TSMonitorSPEigCtxDestroy);
255: }
256: PetscOptionsName("-ts_monitor_sp_swarm","Display particle phase from the DMSwarm","TSMonitorSPSwarm",&opt);
257: if (opt) {
258: TSMonitorSPCtx ctx;
259: PetscInt howoften = 1;
260: PetscOptionsInt("-ts_monitor_sp_swarm","Display particles phase from the DMSwarm","TSMonitorSPSwarm",howoften,&howoften,NULL);
261: TSMonitorSPCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
262: TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode (*)(void**))TSMonitorSPCtxDestroy);
263: }
264: opt = PETSC_FALSE;
265: PetscOptionsName("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",&opt);
266: if (opt) {
267: TSMonitorDrawCtx ctx;
268: PetscInt howoften = 1;
270: PetscOptionsInt("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",howoften,&howoften,NULL);
271: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,"Computed Solution",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
272: TSMonitorSet(ts,TSMonitorDrawSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
273: }
274: opt = PETSC_FALSE;
275: PetscOptionsName("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",&opt);
276: if (opt) {
277: TSMonitorDrawCtx ctx;
278: PetscReal bounds[4];
279: PetscInt n = 4;
280: PetscDraw draw;
281: PetscDrawAxis axis;
283: PetscOptionsRealArray("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",bounds,&n,NULL);
284: if (n != 4) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Must provide bounding box of phase field");
285: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,1,&ctx);
286: PetscViewerDrawGetDraw(ctx->viewer,0,&draw);
287: PetscViewerDrawGetDrawAxis(ctx->viewer,0,&axis);
288: PetscDrawAxisSetLimits(axis,bounds[0],bounds[2],bounds[1],bounds[3]);
289: PetscDrawAxisSetLabels(axis,"Phase Diagram","Variable 1","Variable 2");
290: TSMonitorSet(ts,TSMonitorDrawSolutionPhase,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
291: }
292: opt = PETSC_FALSE;
293: PetscOptionsName("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",&opt);
294: if (opt) {
295: TSMonitorDrawCtx ctx;
296: PetscInt howoften = 1;
298: PetscOptionsInt("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",howoften,&howoften,NULL);
299: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,"Error",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
300: TSMonitorSet(ts,TSMonitorDrawError,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
301: }
302: opt = PETSC_FALSE;
303: PetscOptionsName("-ts_monitor_draw_solution_function","Monitor solution provided by TSMonitorSetSolutionFunction() graphically","TSMonitorDrawSolutionFunction",&opt);
304: if (opt) {
305: TSMonitorDrawCtx ctx;
306: PetscInt howoften = 1;
308: PetscOptionsInt("-ts_monitor_draw_solution_function","Monitor solution provided by TSMonitorSetSolutionFunction() graphically","TSMonitorDrawSolutionFunction",howoften,&howoften,NULL);
309: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,"Solution provided by user function",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
310: TSMonitorSet(ts,TSMonitorDrawSolutionFunction,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
311: }
313: opt = PETSC_FALSE;
314: 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);
315: if (flg) {
316: const char *ptr,*ptr2;
317: char *filetemplate;
318: if (!monfilename[0]) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
319: /* Do some cursory validation of the input. */
320: PetscStrstr(monfilename,"%",(char**)&ptr);
321: if (!ptr) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
322: for (ptr++; ptr && *ptr; ptr++) {
323: PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
324: 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");
325: if (ptr2) break;
326: }
327: PetscStrallocpy(monfilename,&filetemplate);
328: TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);
329: }
331: PetscOptionsString("-ts_monitor_dmda_ray","Display a ray of the solution","None","y=0",dir,16,&flg);
332: if (flg) {
333: TSMonitorDMDARayCtx *rayctx;
334: int ray = 0;
335: DMDirection ddir;
336: DM da;
337: PetscMPIInt rank;
339: if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
340: if (dir[0] == 'x') ddir = DM_X;
341: else if (dir[0] == 'y') ddir = DM_Y;
342: else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
343: sscanf(dir+2,"%d",&ray);
345: PetscInfo2(((PetscObject)ts),"Displaying DMDA ray %c = %d\n",dir[0],ray);
346: PetscNew(&rayctx);
347: TSGetDM(ts,&da);
348: DMDAGetRay(da,ddir,ray,&rayctx->ray,&rayctx->scatter);
349: MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);
350: if (!rank) {
351: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,600,300,&rayctx->viewer);
352: }
353: rayctx->lgctx = NULL;
354: TSMonitorSet(ts,TSMonitorDMDARay,rayctx,TSMonitorDMDARayDestroy);
355: }
356: PetscOptionsString("-ts_monitor_lg_dmda_ray","Display a ray of the solution","None","x=0",dir,16,&flg);
357: if (flg) {
358: TSMonitorDMDARayCtx *rayctx;
359: int ray = 0;
360: DMDirection ddir;
361: DM da;
362: PetscInt howoften = 1;
364: if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
365: if (dir[0] == 'x') ddir = DM_X;
366: else if (dir[0] == 'y') ddir = DM_Y;
367: else SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
368: sscanf(dir+2, "%d", &ray);
370: PetscInfo2(((PetscObject) ts),"Displaying LG DMDA ray %c = %d\n", dir[0], ray);
371: PetscNew(&rayctx);
372: TSGetDM(ts, &da);
373: DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);
374: TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&rayctx->lgctx);
375: TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);
376: }
378: PetscOptionsName("-ts_monitor_envelope","Monitor maximum and minimum value of each component of the solution","TSMonitorEnvelope",&opt);
379: if (opt) {
380: TSMonitorEnvelopeCtx ctx;
382: TSMonitorEnvelopeCtxCreate(ts,&ctx);
383: TSMonitorSet(ts,TSMonitorEnvelope,ctx,(PetscErrorCode (*)(void**))TSMonitorEnvelopeCtxDestroy);
384: }
386: flg = PETSC_FALSE;
387: PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeJacobianDefaultColor", flg, &flg, NULL);
388: if (flg) {
389: DM dm;
390: DMTS tdm;
392: TSGetDM(ts, &dm);
393: DMGetDMTS(dm, &tdm);
394: tdm->ijacobianctx = NULL;
395: TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, 0);
396: PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n");
397: }
399: /* Handle specific TS options */
400: if (ts->ops->setfromoptions) {
401: (*ts->ops->setfromoptions)(PetscOptionsObject,ts);
402: }
404: /* Handle TSAdapt options */
405: TSGetAdapt(ts,&ts->adapt);
406: TSAdaptSetDefaultType(ts->adapt,ts->default_adapt_type);
407: TSAdaptSetFromOptions(PetscOptionsObject,ts->adapt);
409: /* TS trajectory must be set after TS, since it may use some TS options above */
410: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
411: PetscOptionsBool("-ts_save_trajectory","Save the solution at each timestep","TSSetSaveTrajectory",tflg,&tflg,NULL);
412: if (tflg) {
413: TSSetSaveTrajectory(ts);
414: }
416: TSAdjointSetFromOptions(PetscOptionsObject,ts);
418: /* process any options handlers added with PetscObjectAddOptionsHandler() */
419: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)ts);
420: PetscOptionsEnd();
422: if (ts->trajectory) {
423: TSTrajectorySetFromOptions(ts->trajectory,ts);
424: }
426: /* why do we have to do this here and not during TSSetUp? */
427: TSGetSNES(ts,&ts->snes);
428: if (ts->problem_type == TS_LINEAR) {
429: PetscObjectTypeCompareAny((PetscObject)ts->snes,&flg,SNESKSPONLY,SNESKSPTRANSPOSEONLY,"");
430: if (!flg) { SNESSetType(ts->snes,SNESKSPONLY); }
431: }
432: SNESSetFromOptions(ts->snes);
433: return(0);
434: }
436: /*@
437: TSGetTrajectory - Gets the trajectory from a TS if it exists
439: Collective on TS
441: Input Parameters:
442: . ts - the TS context obtained from TSCreate()
444: Output Parameters:
445: . tr - the TSTrajectory object, if it exists
447: Note: This routine should be called after all TS options have been set
449: Level: advanced
451: .seealso: TSGetTrajectory(), TSAdjointSolve(), TSTrajectory, TSTrajectoryCreate()
453: @*/
454: PetscErrorCode TSGetTrajectory(TS ts,TSTrajectory *tr)
455: {
458: *tr = ts->trajectory;
459: return(0);
460: }
462: /*@
463: TSSetSaveTrajectory - Causes the TS to save its solutions as it iterates forward in time in a TSTrajectory object
465: Collective on TS
467: Input Parameters:
468: . ts - the TS context obtained from TSCreate()
470: Options Database:
471: + -ts_save_trajectory - saves the trajectory to a file
472: - -ts_trajectory_type type
474: Note: This routine should be called after all TS options have been set
476: The TSTRAJECTORYVISUALIZATION files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
477: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
479: Level: intermediate
481: .seealso: TSGetTrajectory(), TSAdjointSolve()
483: @*/
484: PetscErrorCode TSSetSaveTrajectory(TS ts)
485: {
490: if (!ts->trajectory) {
491: TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);
492: }
493: return(0);
494: }
496: /*@
497: TSResetTrajectory - Destroys and recreates the internal TSTrajectory object
499: Collective on TS
501: Input Parameters:
502: . ts - the TS context obtained from TSCreate()
504: Level: intermediate
506: .seealso: TSGetTrajectory(), TSAdjointSolve()
508: @*/
509: PetscErrorCode TSResetTrajectory(TS ts)
510: {
515: if (ts->trajectory) {
516: TSTrajectoryDestroy(&ts->trajectory);
517: TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);
518: }
519: return(0);
520: }
522: /*@
523: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
524: set with TSSetRHSJacobian().
526: Collective on TS
528: Input Parameters:
529: + ts - the TS context
530: . t - current timestep
531: - U - input vector
533: Output Parameters:
534: + A - Jacobian matrix
535: . B - optional preconditioning matrix
536: - flag - flag indicating matrix structure
538: Notes:
539: Most users should not need to explicitly call this routine, as it
540: is used internally within the nonlinear solvers.
542: See KSPSetOperators() for important information about setting the
543: flag parameter.
545: Level: developer
547: .seealso: TSSetRHSJacobian(), KSPSetOperators()
548: @*/
549: PetscErrorCode TSComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B)
550: {
551: PetscErrorCode ierr;
552: PetscObjectState Ustate;
553: PetscObjectId Uid;
554: DM dm;
555: DMTS tsdm;
556: TSRHSJacobian rhsjacobianfunc;
557: void *ctx;
558: TSIJacobian ijacobianfunc;
559: TSRHSFunction rhsfunction;
565: TSGetDM(ts,&dm);
566: DMGetDMTS(dm,&tsdm);
567: DMTSGetRHSJacobian(dm,&rhsjacobianfunc,&ctx);
568: DMTSGetIJacobian(dm,&ijacobianfunc,NULL);
569: DMTSGetRHSFunction(dm,&rhsfunction,NULL);
570: PetscObjectStateGet((PetscObject)U,&Ustate);
571: PetscObjectGetId((PetscObject)U,&Uid);
573: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) {
574: /* restore back RHS Jacobian matrices if they have been shifted/scaled */
575: if (A == ts->Arhs) {
576: if (ts->rhsjacobian.shift != 0) {
577: MatShift(A,-ts->rhsjacobian.shift);
578: }
579: if (ts->rhsjacobian.scale != 1.) {
580: MatScale(A,1./ts->rhsjacobian.scale);
581: }
582: }
583: if (B && B == ts->Brhs && A != B) {
584: if (ts->rhsjacobian.shift != 0) {
585: MatShift(B,-ts->rhsjacobian.shift);
586: }
587: if (ts->rhsjacobian.scale != 1.) {
588: MatScale(B,1./ts->rhsjacobian.scale);
589: }
590: }
591: ts->rhsjacobian.shift = 0;
592: ts->rhsjacobian.scale = 1.;
593: return(0);
594: }
596: if (!rhsjacobianfunc && !ijacobianfunc) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
598: if (ts->rhsjacobian.reuse) {
599: if (A == ts->Arhs) {
600: /* MatScale has a short path for this case.
601: However, this code path is taken the first time TSComputeRHSJacobian is called
602: and the matrices have not assembled yet */
603: if (ts->rhsjacobian.shift != 0) {
604: MatShift(A,-ts->rhsjacobian.shift);
605: }
606: if (ts->rhsjacobian.scale != 1.) {
607: MatScale(A,1./ts->rhsjacobian.scale);
608: }
609: }
610: if (B && B == ts->Brhs && A != B) {
611: if (ts->rhsjacobian.shift != 0) {
612: MatShift(B,-ts->rhsjacobian.shift);
613: }
614: if (ts->rhsjacobian.scale != 1.) {
615: MatScale(B,1./ts->rhsjacobian.scale);
616: }
617: }
618: }
620: if (rhsjacobianfunc) {
621: PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);
622: PetscStackPush("TS user Jacobian function");
623: (*rhsjacobianfunc)(ts,t,U,A,B,ctx);
624: PetscStackPop;
625: PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);
626: } else {
627: MatZeroEntries(A);
628: if (B && A != B) {MatZeroEntries(B);}
629: }
630: ts->rhsjacobian.time = t;
631: ts->rhsjacobian.shift = 0;
632: ts->rhsjacobian.scale = 1.;
633: PetscObjectGetId((PetscObject)U,&ts->rhsjacobian.Xid);
634: PetscObjectStateGet((PetscObject)U,&ts->rhsjacobian.Xstate);
635: return(0);
636: }
638: /*@
639: TSComputeRHSFunction - Evaluates the right-hand-side function.
641: Collective on TS
643: Input Parameters:
644: + ts - the TS context
645: . t - current time
646: - U - state vector
648: Output Parameter:
649: . y - right hand side
651: Note:
652: Most users should not need to explicitly call this routine, as it
653: is used internally within the nonlinear solvers.
655: Level: developer
657: .seealso: TSSetRHSFunction(), TSComputeIFunction()
658: @*/
659: PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y)
660: {
662: TSRHSFunction rhsfunction;
663: TSIFunction ifunction;
664: void *ctx;
665: DM dm;
671: TSGetDM(ts,&dm);
672: DMTSGetRHSFunction(dm,&rhsfunction,&ctx);
673: DMTSGetIFunction(dm,&ifunction,NULL);
675: if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
677: PetscLogEventBegin(TS_FunctionEval,ts,U,y,0);
678: if (rhsfunction) {
679: VecLockReadPush(U);
680: PetscStackPush("TS user right-hand-side function");
681: (*rhsfunction)(ts,t,U,y,ctx);
682: PetscStackPop;
683: VecLockReadPop(U);
684: } else {
685: VecZeroEntries(y);
686: }
688: PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);
689: return(0);
690: }
692: /*@
693: TSComputeSolutionFunction - Evaluates the solution function.
695: Collective on TS
697: Input Parameters:
698: + ts - the TS context
699: - t - current time
701: Output Parameter:
702: . U - the solution
704: Note:
705: Most users should not need to explicitly call this routine, as it
706: is used internally within the nonlinear solvers.
708: Level: developer
710: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
711: @*/
712: PetscErrorCode TSComputeSolutionFunction(TS ts,PetscReal t,Vec U)
713: {
714: PetscErrorCode ierr;
715: TSSolutionFunction solutionfunction;
716: void *ctx;
717: DM dm;
722: TSGetDM(ts,&dm);
723: DMTSGetSolutionFunction(dm,&solutionfunction,&ctx);
725: if (solutionfunction) {
726: PetscStackPush("TS user solution function");
727: (*solutionfunction)(ts,t,U,ctx);
728: PetscStackPop;
729: }
730: return(0);
731: }
732: /*@
733: TSComputeForcingFunction - Evaluates the forcing function.
735: Collective on TS
737: Input Parameters:
738: + ts - the TS context
739: - t - current time
741: Output Parameter:
742: . U - the function value
744: Note:
745: Most users should not need to explicitly call this routine, as it
746: is used internally within the nonlinear solvers.
748: Level: developer
750: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
751: @*/
752: PetscErrorCode TSComputeForcingFunction(TS ts,PetscReal t,Vec U)
753: {
754: PetscErrorCode ierr, (*forcing)(TS,PetscReal,Vec,void*);
755: void *ctx;
756: DM dm;
761: TSGetDM(ts,&dm);
762: DMTSGetForcingFunction(dm,&forcing,&ctx);
764: if (forcing) {
765: PetscStackPush("TS user forcing function");
766: (*forcing)(ts,t,U,ctx);
767: PetscStackPop;
768: }
769: return(0);
770: }
772: static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs)
773: {
774: Vec F;
778: *Frhs = NULL;
779: TSGetIFunction(ts,&F,NULL,NULL);
780: if (!ts->Frhs) {
781: VecDuplicate(F,&ts->Frhs);
782: }
783: *Frhs = ts->Frhs;
784: return(0);
785: }
787: PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs)
788: {
789: Mat A,B;
791: TSIJacobian ijacobian;
794: if (Arhs) *Arhs = NULL;
795: if (Brhs) *Brhs = NULL;
796: TSGetIJacobian(ts,&A,&B,&ijacobian,NULL);
797: if (Arhs) {
798: if (!ts->Arhs) {
799: if (ijacobian) {
800: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);
801: } else {
802: ts->Arhs = A;
803: PetscObjectReference((PetscObject)A);
804: }
805: } else {
806: PetscBool flg;
807: SNESGetUseMatrixFree(ts->snes,NULL,&flg);
808: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
809: if (flg && !ijacobian && ts->Arhs == ts->Brhs){
810: PetscObjectDereference((PetscObject)ts->Arhs);
811: ts->Arhs = A;
812: PetscObjectReference((PetscObject)A);
813: }
814: }
815: *Arhs = ts->Arhs;
816: }
817: if (Brhs) {
818: if (!ts->Brhs) {
819: if (A != B) {
820: if (ijacobian) {
821: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);
822: } else {
823: ts->Brhs = B;
824: PetscObjectReference((PetscObject)B);
825: }
826: } else {
827: PetscObjectReference((PetscObject)ts->Arhs);
828: ts->Brhs = ts->Arhs;
829: }
830: }
831: *Brhs = ts->Brhs;
832: }
833: return(0);
834: }
836: /*@
837: TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0
839: Collective on TS
841: Input Parameters:
842: + ts - the TS context
843: . t - current time
844: . U - state vector
845: . Udot - time derivative of state vector
846: - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate
848: Output Parameter:
849: . Y - right hand side
851: Note:
852: Most users should not need to explicitly call this routine, as it
853: is used internally within the nonlinear solvers.
855: If the user did did not write their equations in implicit form, this
856: function recasts them in implicit form.
858: Level: developer
860: .seealso: TSSetIFunction(), TSComputeRHSFunction()
861: @*/
862: PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec Y,PetscBool imex)
863: {
865: TSIFunction ifunction;
866: TSRHSFunction rhsfunction;
867: void *ctx;
868: DM dm;
876: TSGetDM(ts,&dm);
877: DMTSGetIFunction(dm,&ifunction,&ctx);
878: DMTSGetRHSFunction(dm,&rhsfunction,NULL);
880: if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
882: PetscLogEventBegin(TS_FunctionEval,ts,U,Udot,Y);
883: if (ifunction) {
884: PetscStackPush("TS user implicit function");
885: (*ifunction)(ts,t,U,Udot,Y,ctx);
886: PetscStackPop;
887: }
888: if (imex) {
889: if (!ifunction) {
890: VecCopy(Udot,Y);
891: }
892: } else if (rhsfunction) {
893: if (ifunction) {
894: Vec Frhs;
895: TSGetRHSVec_Private(ts,&Frhs);
896: TSComputeRHSFunction(ts,t,U,Frhs);
897: VecAXPY(Y,-1,Frhs);
898: } else {
899: TSComputeRHSFunction(ts,t,U,Y);
900: VecAYPX(Y,-1,Udot);
901: }
902: }
903: PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);
904: return(0);
905: }
907: /*@
908: TSComputeIJacobian - Evaluates the Jacobian of the DAE
910: Collective on TS
912: Input
913: Input Parameters:
914: + ts - the TS context
915: . t - current timestep
916: . U - state vector
917: . Udot - time derivative of state vector
918: . shift - shift to apply, see note below
919: - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate
921: Output Parameters:
922: + A - Jacobian matrix
923: - B - matrix from which the preconditioner is constructed; often the same as A
925: Notes:
926: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
928: dF/dU + shift*dF/dUdot
930: Most users should not need to explicitly call this routine, as it
931: is used internally within the nonlinear solvers.
933: Level: developer
935: .seealso: TSSetIJacobian()
936: @*/
937: PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,PetscBool imex)
938: {
940: TSIJacobian ijacobian;
941: TSRHSJacobian rhsjacobian;
942: DM dm;
943: void *ctx;
954: TSGetDM(ts,&dm);
955: DMTSGetIJacobian(dm,&ijacobian,&ctx);
956: DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);
958: if (!rhsjacobian && !ijacobian) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
960: PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);
961: if (ijacobian) {
962: PetscStackPush("TS user implicit Jacobian");
963: (*ijacobian)(ts,t,U,Udot,shift,A,B,ctx);
964: PetscStackPop;
965: }
966: if (imex) {
967: if (!ijacobian) { /* system was written as Udot = G(t,U) */
968: PetscBool assembled;
969: if (rhsjacobian) {
970: Mat Arhs = NULL;
971: TSGetRHSMats_Private(ts,&Arhs,NULL);
972: if (A == Arhs) {
973: if (rhsjacobian == TSComputeRHSJacobianConstant) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Unsupported operation! cannot use TSComputeRHSJacobianConstant");
974: ts->rhsjacobian.time = PETSC_MIN_REAL;
975: }
976: }
977: MatZeroEntries(A);
978: MatAssembled(A,&assembled);
979: if (!assembled) {
980: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
981: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
982: }
983: MatShift(A,shift);
984: if (A != B) {
985: MatZeroEntries(B);
986: MatAssembled(B,&assembled);
987: if (!assembled) {
988: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
989: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
990: }
991: MatShift(B,shift);
992: }
993: }
994: } else {
995: Mat Arhs = NULL,Brhs = NULL;
996: if (rhsjacobian) {
997: TSGetRHSMats_Private(ts,&Arhs,&Brhs);
998: TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);
999: }
1000: if (Arhs == A) { /* No IJacobian, so we only have the RHS matrix */
1001: PetscBool flg;
1002: ts->rhsjacobian.scale = -1;
1003: ts->rhsjacobian.shift = shift;
1004: SNESGetUseMatrixFree(ts->snes,NULL,&flg);
1005: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
1006: if (!flg) {
1007: MatScale(A,-1);
1008: MatShift(A,shift);
1009: }
1010: if (A != B) {
1011: MatScale(B,-1);
1012: MatShift(B,shift);
1013: }
1014: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
1015: MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
1016: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
1017: MatZeroEntries(A);
1018: MatShift(A,shift);
1019: if (A != B) {
1020: MatZeroEntries(B);
1021: MatShift(B,shift);
1022: }
1023: }
1024: MatAXPY(A,-1,Arhs,axpy);
1025: if (A != B) {
1026: MatAXPY(B,-1,Brhs,axpy);
1027: }
1028: }
1029: }
1030: PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);
1031: return(0);
1032: }
1034: /*@C
1035: TSSetRHSFunction - Sets the routine for evaluating the function,
1036: where U_t = G(t,u).
1038: Logically Collective on TS
1040: Input Parameters:
1041: + ts - the TS context obtained from TSCreate()
1042: . r - vector to put the computed right hand side (or NULL to have it created)
1043: . f - routine for evaluating the right-hand-side function
1044: - ctx - [optional] user-defined context for private data for the
1045: function evaluation routine (may be NULL)
1047: Calling sequence of func:
1048: $ PetscErrorCode func (TS ts,PetscReal t,Vec u,Vec F,void *ctx);
1050: + t - current timestep
1051: . u - input vector
1052: . F - function vector
1053: - ctx - [optional] user-defined function context
1055: Level: beginner
1057: Notes:
1058: You must call this function or TSSetIFunction() to define your ODE. You cannot use this function when solving a DAE.
1060: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSSetIFunction()
1061: @*/
1062: PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx)
1063: {
1065: SNES snes;
1066: Vec ralloc = NULL;
1067: DM dm;
1073: TSGetDM(ts,&dm);
1074: DMTSSetRHSFunction(dm,f,ctx);
1075: TSGetSNES(ts,&snes);
1076: if (!r && !ts->dm && ts->vec_sol) {
1077: VecDuplicate(ts->vec_sol,&ralloc);
1078: r = ralloc;
1079: }
1080: SNESSetFunction(snes,r,SNESTSFormFunction,ts);
1081: VecDestroy(&ralloc);
1082: return(0);
1083: }
1085: /*@C
1086: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1088: Logically Collective on TS
1090: Input Parameters:
1091: + ts - the TS context obtained from TSCreate()
1092: . f - routine for evaluating the solution
1093: - ctx - [optional] user-defined context for private data for the
1094: function evaluation routine (may be NULL)
1096: Calling sequence of func:
1097: $ PetscErrorCode func (TS ts,PetscReal t,Vec u,void *ctx);
1099: + t - current timestep
1100: . u - output vector
1101: - ctx - [optional] user-defined function context
1103: Options Database:
1104: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided TSSetSolutionFunction()
1105: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided TSSetSolutionFunction()
1107: Notes:
1108: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1109: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1110: create closed-form solutions with non-physical forcing terms.
1112: For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history.
1114: Level: beginner
1116: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetForcingFunction(), TSSetSolution(), TSGetSolution(), TSMonitorLGError(), TSMonitorDrawError()
1117: @*/
1118: PetscErrorCode TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx)
1119: {
1121: DM dm;
1125: TSGetDM(ts,&dm);
1126: DMTSSetSolutionFunction(dm,f,ctx);
1127: return(0);
1128: }
1130: /*@C
1131: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1133: Logically Collective on TS
1135: Input Parameters:
1136: + ts - the TS context obtained from TSCreate()
1137: . func - routine for evaluating the forcing function
1138: - ctx - [optional] user-defined context for private data for the
1139: function evaluation routine (may be NULL)
1141: Calling sequence of func:
1142: $ PetscErrorCode func (TS ts,PetscReal t,Vec f,void *ctx);
1144: + t - current timestep
1145: . f - output vector
1146: - ctx - [optional] user-defined function context
1148: Notes:
1149: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1150: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1151: definition of the problem you are solving and hence possibly introducing bugs.
1153: This replaces the ODE F(u,u_t,t) = 0 the TS is solving with F(u,u_t,t) - func(t) = 0
1155: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1156: parameters can be passed in the ctx variable.
1158: For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history.
1160: Level: beginner
1162: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetSolutionFunction()
1163: @*/
1164: PetscErrorCode TSSetForcingFunction(TS ts,TSForcingFunction func,void *ctx)
1165: {
1167: DM dm;
1171: TSGetDM(ts,&dm);
1172: DMTSSetForcingFunction(dm,func,ctx);
1173: return(0);
1174: }
1176: /*@C
1177: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1178: where U_t = G(U,t), as well as the location to store the matrix.
1180: Logically Collective on TS
1182: Input Parameters:
1183: + ts - the TS context obtained from TSCreate()
1184: . Amat - (approximate) Jacobian matrix
1185: . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1186: . f - the Jacobian evaluation routine
1187: - ctx - [optional] user-defined context for private data for the
1188: Jacobian evaluation routine (may be NULL)
1190: Calling sequence of f:
1191: $ PetscErrorCode func (TS ts,PetscReal t,Vec u,Mat A,Mat B,void *ctx);
1193: + t - current timestep
1194: . u - input vector
1195: . Amat - (approximate) Jacobian matrix
1196: . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1197: - ctx - [optional] user-defined context for matrix evaluation routine
1199: Notes:
1200: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1202: The TS solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1203: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1205: Level: beginner
1207: .seealso: SNESComputeJacobianDefaultColor(), TSSetRHSFunction(), TSRHSJacobianSetReuse(), TSSetIJacobian()
1209: @*/
1210: PetscErrorCode TSSetRHSJacobian(TS ts,Mat Amat,Mat Pmat,TSRHSJacobian f,void *ctx)
1211: {
1213: SNES snes;
1214: DM dm;
1215: TSIJacobian ijacobian;
1224: TSGetDM(ts,&dm);
1225: DMTSSetRHSJacobian(dm,f,ctx);
1226: if (f == TSComputeRHSJacobianConstant) {
1227: /* Handle this case automatically for the user; otherwise user should call themselves. */
1228: TSRHSJacobianSetReuse(ts,PETSC_TRUE);
1229: }
1230: DMTSGetIJacobian(dm,&ijacobian,NULL);
1231: TSGetSNES(ts,&snes);
1232: if (!ijacobian) {
1233: SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1234: }
1235: if (Amat) {
1236: PetscObjectReference((PetscObject)Amat);
1237: MatDestroy(&ts->Arhs);
1238: ts->Arhs = Amat;
1239: }
1240: if (Pmat) {
1241: PetscObjectReference((PetscObject)Pmat);
1242: MatDestroy(&ts->Brhs);
1243: ts->Brhs = Pmat;
1244: }
1245: return(0);
1246: }
1248: /*@C
1249: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1251: Logically Collective on TS
1253: Input Parameters:
1254: + ts - the TS context obtained from TSCreate()
1255: . r - vector to hold the residual (or NULL to have it created internally)
1256: . f - the function evaluation routine
1257: - ctx - user-defined context for private data for the function evaluation routine (may be NULL)
1259: Calling sequence of f:
1260: $ PetscErrorCode f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx);
1262: + t - time at step/stage being solved
1263: . u - state vector
1264: . u_t - time derivative of state vector
1265: . F - function vector
1266: - ctx - [optional] user-defined context for matrix evaluation routine
1268: Important:
1269: The user MUST call either this routine or TSSetRHSFunction() to define the ODE. When solving DAEs you must use this function.
1271: Level: beginner
1273: .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian()
1274: @*/
1275: PetscErrorCode TSSetIFunction(TS ts,Vec r,TSIFunction f,void *ctx)
1276: {
1278: SNES snes;
1279: Vec ralloc = NULL;
1280: DM dm;
1286: TSGetDM(ts,&dm);
1287: DMTSSetIFunction(dm,f,ctx);
1289: TSGetSNES(ts,&snes);
1290: if (!r && !ts->dm && ts->vec_sol) {
1291: VecDuplicate(ts->vec_sol,&ralloc);
1292: r = ralloc;
1293: }
1294: SNESSetFunction(snes,r,SNESTSFormFunction,ts);
1295: VecDestroy(&ralloc);
1296: return(0);
1297: }
1299: /*@C
1300: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it.
1302: Not Collective
1304: Input Parameter:
1305: . ts - the TS context
1307: Output Parameter:
1308: + r - vector to hold residual (or NULL)
1309: . func - the function to compute residual (or NULL)
1310: - ctx - the function context (or NULL)
1312: Level: advanced
1314: .seealso: TSSetIFunction(), SNESGetFunction()
1315: @*/
1316: PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx)
1317: {
1319: SNES snes;
1320: DM dm;
1324: TSGetSNES(ts,&snes);
1325: SNESGetFunction(snes,r,NULL,NULL);
1326: TSGetDM(ts,&dm);
1327: DMTSGetIFunction(dm,func,ctx);
1328: return(0);
1329: }
1331: /*@C
1332: TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it.
1334: Not Collective
1336: Input Parameter:
1337: . ts - the TS context
1339: Output Parameter:
1340: + r - vector to hold computed right hand side (or NULL)
1341: . func - the function to compute right hand side (or NULL)
1342: - ctx - the function context (or NULL)
1344: Level: advanced
1346: .seealso: TSSetRHSFunction(), SNESGetFunction()
1347: @*/
1348: PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx)
1349: {
1351: SNES snes;
1352: DM dm;
1356: TSGetSNES(ts,&snes);
1357: SNESGetFunction(snes,r,NULL,NULL);
1358: TSGetDM(ts,&dm);
1359: DMTSGetRHSFunction(dm,func,ctx);
1360: return(0);
1361: }
1363: /*@C
1364: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1365: provided with TSSetIFunction().
1367: Logically Collective on TS
1369: Input Parameters:
1370: + ts - the TS context obtained from TSCreate()
1371: . Amat - (approximate) Jacobian matrix
1372: . Pmat - matrix used to compute preconditioner (usually the same as Amat)
1373: . f - the Jacobian evaluation routine
1374: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
1376: Calling sequence of f:
1377: $ PetscErrorCode f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat Amat,Mat Pmat,void *ctx);
1379: + t - time at step/stage being solved
1380: . U - state vector
1381: . U_t - time derivative of state vector
1382: . a - shift
1383: . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1384: . Pmat - matrix used for constructing preconditioner, usually the same as Amat
1385: - ctx - [optional] user-defined context for matrix evaluation routine
1387: Notes:
1388: The matrices Amat and Pmat are exactly the matrices that are used by SNES for the nonlinear solve.
1390: If you know the operator Amat has a null space you can use MatSetNullSpace() and MatSetTransposeNullSpace() to supply the null
1391: space to Amat and the KSP solvers will automatically use that null space as needed during the solution process.
1393: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1394: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1395: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1396: a and vector W depend on the integration method, step size, and past states. For example with
1397: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1398: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1400: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1402: The TS solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1403: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1405: Level: beginner
1407: .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESComputeJacobianDefaultColor(), SNESComputeJacobianDefault(), TSSetRHSFunction()
1409: @*/
1410: PetscErrorCode TSSetIJacobian(TS ts,Mat Amat,Mat Pmat,TSIJacobian f,void *ctx)
1411: {
1413: SNES snes;
1414: DM dm;
1423: TSGetDM(ts,&dm);
1424: DMTSSetIJacobian(dm,f,ctx);
1426: TSGetSNES(ts,&snes);
1427: SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1428: return(0);
1429: }
1431: /*@
1432: TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating. Without this flag, TS will change the sign and
1433: shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute
1434: the entire Jacobian. The reuse flag must be set if the evaluation function will assume that the matrix entries have
1435: not been changed by the TS.
1437: Logically Collective
1439: Input Arguments:
1440: + ts - TS context obtained from TSCreate()
1441: - reuse - PETSC_TRUE if the RHS Jacobian
1443: Level: intermediate
1445: .seealso: TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
1446: @*/
1447: PetscErrorCode TSRHSJacobianSetReuse(TS ts,PetscBool reuse)
1448: {
1450: ts->rhsjacobian.reuse = reuse;
1451: return(0);
1452: }
1454: /*@C
1455: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1457: Logically Collective on TS
1459: Input Parameters:
1460: + ts - the TS context obtained from TSCreate()
1461: . F - vector to hold the residual (or NULL to have it created internally)
1462: . fun - the function evaluation routine
1463: - ctx - user-defined context for private data for the function evaluation routine (may be NULL)
1465: Calling sequence of fun:
1466: $ PetscErrorCode fun(TS ts,PetscReal t,Vec U,Vec U_t,Vec U_tt,Vec F,ctx);
1468: + t - time at step/stage being solved
1469: . U - state vector
1470: . U_t - time derivative of state vector
1471: . U_tt - second time derivative of state vector
1472: . F - function vector
1473: - ctx - [optional] user-defined context for matrix evaluation routine (may be NULL)
1475: Level: beginner
1477: .seealso: TSSetI2Jacobian()
1478: @*/
1479: PetscErrorCode TSSetI2Function(TS ts,Vec F,TSI2Function fun,void *ctx)
1480: {
1481: DM dm;
1487: TSSetIFunction(ts,F,NULL,NULL);
1488: TSGetDM(ts,&dm);
1489: DMTSSetI2Function(dm,fun,ctx);
1490: return(0);
1491: }
1493: /*@C
1494: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/contex to compute it.
1496: Not Collective
1498: Input Parameter:
1499: . ts - the TS context
1501: Output Parameter:
1502: + r - vector to hold residual (or NULL)
1503: . fun - the function to compute residual (or NULL)
1504: - ctx - the function context (or NULL)
1506: Level: advanced
1508: .seealso: TSSetI2Function(), SNESGetFunction()
1509: @*/
1510: PetscErrorCode TSGetI2Function(TS ts,Vec *r,TSI2Function *fun,void **ctx)
1511: {
1513: SNES snes;
1514: DM dm;
1518: TSGetSNES(ts,&snes);
1519: SNESGetFunction(snes,r,NULL,NULL);
1520: TSGetDM(ts,&dm);
1521: DMTSGetI2Function(dm,fun,ctx);
1522: return(0);
1523: }
1525: /*@C
1526: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1527: where F(t,U,U_t,U_tt) is the function you provided with TSSetI2Function().
1529: Logically Collective on TS
1531: Input Parameters:
1532: + ts - the TS context obtained from TSCreate()
1533: . J - Jacobian matrix
1534: . P - preconditioning matrix for J (may be same as J)
1535: . jac - the Jacobian evaluation routine
1536: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
1538: Calling sequence of jac:
1539: $ PetscErrorCode jac(TS ts,PetscReal t,Vec U,Vec U_t,Vec U_tt,PetscReal v,PetscReal a,Mat J,Mat P,void *ctx);
1541: + t - time at step/stage being solved
1542: . U - state vector
1543: . U_t - time derivative of state vector
1544: . U_tt - second time derivative of state vector
1545: . v - shift for U_t
1546: . a - shift for U_tt
1547: . 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
1548: . P - preconditioning matrix for J, may be same as J
1549: - ctx - [optional] user-defined context for matrix evaluation routine
1551: Notes:
1552: The matrices J and P are exactly the matrices that are used by SNES for the nonlinear solve.
1554: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1555: 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.
1556: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1557: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1559: Level: beginner
1561: .seealso: TSSetI2Function()
1562: @*/
1563: PetscErrorCode TSSetI2Jacobian(TS ts,Mat J,Mat P,TSI2Jacobian jac,void *ctx)
1564: {
1565: DM dm;
1572: TSSetIJacobian(ts,J,P,NULL,NULL);
1573: TSGetDM(ts,&dm);
1574: DMTSSetI2Jacobian(dm,jac,ctx);
1575: return(0);
1576: }
1578: /*@C
1579: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1581: Not Collective, but parallel objects are returned if TS is parallel
1583: Input Parameter:
1584: . ts - The TS context obtained from TSCreate()
1586: Output Parameters:
1587: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1588: . P - The matrix from which the preconditioner is constructed, often the same as J
1589: . jac - The function to compute the Jacobian matrices
1590: - ctx - User-defined context for Jacobian evaluation routine
1592: Notes:
1593: You can pass in NULL for any return argument you do not need.
1595: Level: advanced
1597: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()
1599: @*/
1600: PetscErrorCode TSGetI2Jacobian(TS ts,Mat *J,Mat *P,TSI2Jacobian *jac,void **ctx)
1601: {
1603: SNES snes;
1604: DM dm;
1607: TSGetSNES(ts,&snes);
1608: SNESSetUpMatrices(snes);
1609: SNESGetJacobian(snes,J,P,NULL,NULL);
1610: TSGetDM(ts,&dm);
1611: DMTSGetI2Jacobian(dm,jac,ctx);
1612: return(0);
1613: }
1615: /*@
1616: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1618: Collective on TS
1620: Input Parameters:
1621: + ts - the TS context
1622: . t - current time
1623: . U - state vector
1624: . V - time derivative of state vector (U_t)
1625: - A - second time derivative of state vector (U_tt)
1627: Output Parameter:
1628: . F - the residual vector
1630: Note:
1631: Most users should not need to explicitly call this routine, as it
1632: is used internally within the nonlinear solvers.
1634: Level: developer
1636: .seealso: TSSetI2Function()
1637: @*/
1638: PetscErrorCode TSComputeI2Function(TS ts,PetscReal t,Vec U,Vec V,Vec A,Vec F)
1639: {
1640: DM dm;
1641: TSI2Function I2Function;
1642: void *ctx;
1643: TSRHSFunction rhsfunction;
1653: TSGetDM(ts,&dm);
1654: DMTSGetI2Function(dm,&I2Function,&ctx);
1655: DMTSGetRHSFunction(dm,&rhsfunction,NULL);
1657: if (!I2Function) {
1658: TSComputeIFunction(ts,t,U,A,F,PETSC_FALSE);
1659: return(0);
1660: }
1662: PetscLogEventBegin(TS_FunctionEval,ts,U,V,F);
1664: PetscStackPush("TS user implicit function");
1665: I2Function(ts,t,U,V,A,F,ctx);
1666: PetscStackPop;
1668: if (rhsfunction) {
1669: Vec Frhs;
1670: TSGetRHSVec_Private(ts,&Frhs);
1671: TSComputeRHSFunction(ts,t,U,Frhs);
1672: VecAXPY(F,-1,Frhs);
1673: }
1675: PetscLogEventEnd(TS_FunctionEval,ts,U,V,F);
1676: return(0);
1677: }
1679: /*@
1680: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1682: Collective on TS
1684: Input Parameters:
1685: + ts - the TS context
1686: . t - current timestep
1687: . U - state vector
1688: . V - time derivative of state vector
1689: . A - second time derivative of state vector
1690: . shiftV - shift to apply, see note below
1691: - shiftA - shift to apply, see note below
1693: Output Parameters:
1694: + J - Jacobian matrix
1695: - P - optional preconditioning matrix
1697: Notes:
1698: If F(t,U,V,A)=0 is the DAE, the required Jacobian is
1700: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1702: Most users should not need to explicitly call this routine, as it
1703: is used internally within the nonlinear solvers.
1705: Level: developer
1707: .seealso: TSSetI2Jacobian()
1708: @*/
1709: PetscErrorCode TSComputeI2Jacobian(TS ts,PetscReal t,Vec U,Vec V,Vec A,PetscReal shiftV,PetscReal shiftA,Mat J,Mat P)
1710: {
1711: DM dm;
1712: TSI2Jacobian I2Jacobian;
1713: void *ctx;
1714: TSRHSJacobian rhsjacobian;
1725: TSGetDM(ts,&dm);
1726: DMTSGetI2Jacobian(dm,&I2Jacobian,&ctx);
1727: DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);
1729: if (!I2Jacobian) {
1730: TSComputeIJacobian(ts,t,U,A,shiftA,J,P,PETSC_FALSE);
1731: return(0);
1732: }
1734: PetscLogEventBegin(TS_JacobianEval,ts,U,J,P);
1736: PetscStackPush("TS user implicit Jacobian");
1737: I2Jacobian(ts,t,U,V,A,shiftV,shiftA,J,P,ctx);
1738: PetscStackPop;
1740: if (rhsjacobian) {
1741: Mat Jrhs,Prhs; MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
1742: TSGetRHSMats_Private(ts,&Jrhs,&Prhs);
1743: TSComputeRHSJacobian(ts,t,U,Jrhs,Prhs);
1744: MatAXPY(J,-1,Jrhs,axpy);
1745: if (P != J) {MatAXPY(P,-1,Prhs,axpy);}
1746: }
1748: PetscLogEventEnd(TS_JacobianEval,ts,U,J,P);
1749: return(0);
1750: }
1752: /*@C
1753: TSSetTransientVariable - sets function to transform from state to transient variables
1755: Logically Collective
1757: Input Arguments:
1758: + ts - time stepping context on which to change the transient variable
1759: . tvar - a function that transforms in-place to transient variables
1760: - ctx - a context for tvar
1762: Level: advanced
1764: Notes:
1765: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., TSBDF)
1766: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1767: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1768: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1769: evaluated via the chain rule, as in
1771: dF/dP + shift * dF/dCdot dC/dP.
1773: .seealso: DMTSSetTransientVariable(), DMTSGetTransientVariable(), TSSetIFunction(), TSSetIJacobian()
1774: @*/
1775: PetscErrorCode TSSetTransientVariable(TS ts,TSTransientVariable tvar,void *ctx)
1776: {
1778: DM dm;
1782: TSGetDM(ts,&dm);
1783: DMTSSetTransientVariable(dm,tvar,ctx);
1784: return(0);
1785: }
1787: /*@
1788: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1790: Logically Collective
1792: Input Parameters:
1793: + ts - TS on which to compute
1794: - U - state vector to be transformed to transient variables
1796: Output Parameters:
1797: . C - transient (conservative) variable
1799: Developer Notes:
1800: If DMTSSetTransientVariable() has not been called, then C is not modified in this routine and C=NULL is allowed.
1801: This makes it safe to call without a guard. One can use TSHasTransientVariable() to check if transient variables are
1802: being used.
1804: Level: developer
1806: .seealso: DMTSSetTransientVariable(), TSComputeIFunction(), TSComputeIJacobian()
1807: @*/
1808: PetscErrorCode TSComputeTransientVariable(TS ts,Vec U,Vec C)
1809: {
1811: DM dm;
1812: DMTS dmts;
1817: TSGetDM(ts,&dm);
1818: DMGetDMTS(dm,&dmts);
1819: if (dmts->ops->transientvar) {
1821: (*dmts->ops->transientvar)(ts,U,C,dmts->transientvarctx);
1822: }
1823: return(0);
1824: }
1826: /*@
1827: TSHasTransientVariable - determine whether transient variables have been set
1829: Logically Collective
1831: Input Parameters:
1832: . ts - TS on which to compute
1834: Output Parameters:
1835: . has - PETSC_TRUE if transient variables have been set
1837: Level: developer
1839: .seealso: DMTSSetTransientVariable(), TSComputeTransientVariable()
1840: @*/
1841: PetscErrorCode TSHasTransientVariable(TS ts,PetscBool *has)
1842: {
1844: DM dm;
1845: DMTS dmts;
1849: TSGetDM(ts,&dm);
1850: DMGetDMTS(dm,&dmts);
1851: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1852: return(0);
1853: }
1855: /*@
1856: TS2SetSolution - Sets the initial solution and time derivative vectors
1857: for use by the TS routines handling second order equations.
1859: Logically Collective on TS
1861: Input Parameters:
1862: + ts - the TS context obtained from TSCreate()
1863: . u - the solution vector
1864: - v - the time derivative vector
1866: Level: beginner
1868: @*/
1869: PetscErrorCode TS2SetSolution(TS ts,Vec u,Vec v)
1870: {
1877: TSSetSolution(ts,u);
1878: PetscObjectReference((PetscObject)v);
1879: VecDestroy(&ts->vec_dot);
1880: ts->vec_dot = v;
1881: return(0);
1882: }
1884: /*@
1885: TS2GetSolution - Returns the solution and time derivative at the present timestep
1886: for second order equations. It is valid to call this routine inside the function
1887: that you are evaluating in order to move to the new timestep. This vector not
1888: changed until the solution at the next timestep has been calculated.
1890: Not Collective, but Vec returned is parallel if TS is parallel
1892: Input Parameter:
1893: . ts - the TS context obtained from TSCreate()
1895: Output Parameter:
1896: + u - the vector containing the solution
1897: - v - the vector containing the time derivative
1899: Level: intermediate
1901: .seealso: TS2SetSolution(), TSGetTimeStep(), TSGetTime()
1903: @*/
1904: PetscErrorCode TS2GetSolution(TS ts,Vec *u,Vec *v)
1905: {
1910: if (u) *u = ts->vec_sol;
1911: if (v) *v = ts->vec_dot;
1912: return(0);
1913: }
1915: /*@C
1916: TSLoad - Loads a KSP that has been stored in binary with KSPView().
1918: Collective on PetscViewer
1920: Input Parameters:
1921: + newdm - the newly loaded TS, this needs to have been created with TSCreate() or
1922: some related function before a call to TSLoad().
1923: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1925: Level: intermediate
1927: Notes:
1928: The type is determined by the data in the file, any type set into the TS before this call is ignored.
1930: Notes for advanced users:
1931: Most users should not need to know the details of the binary storage
1932: format, since TSLoad() and TSView() completely hide these details.
1933: But for anyone who's interested, the standard binary matrix storage
1934: format is
1935: .vb
1936: has not yet been determined
1937: .ve
1939: .seealso: PetscViewerBinaryOpen(), TSView(), MatLoad(), VecLoad()
1940: @*/
1941: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1942: {
1944: PetscBool isbinary;
1945: PetscInt classid;
1946: char type[256];
1947: DMTS sdm;
1948: DM dm;
1953: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1954: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1956: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1957: if (classid != TS_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Not TS next in file");
1958: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1959: TSSetType(ts, type);
1960: if (ts->ops->load) {
1961: (*ts->ops->load)(ts,viewer);
1962: }
1963: DMCreate(PetscObjectComm((PetscObject)ts),&dm);
1964: DMLoad(dm,viewer);
1965: TSSetDM(ts,dm);
1966: DMCreateGlobalVector(ts->dm,&ts->vec_sol);
1967: VecLoad(ts->vec_sol,viewer);
1968: DMGetDMTS(ts->dm,&sdm);
1969: DMTSLoad(sdm,viewer);
1970: return(0);
1971: }
1973: #include <petscdraw.h>
1974: #if defined(PETSC_HAVE_SAWS)
1975: #include <petscviewersaws.h>
1976: #endif
1978: /*@C
1979: TSViewFromOptions - View from Options
1981: Collective on TS
1983: Input Parameters:
1984: + A - the Section 1.5 Writing Application Codes with PETSc ordering context
1985: . obj - Optional object
1986: - name - command line option
1988: Level: intermediate
1989: .seealso: TS, TSView, PetscObjectViewFromOptions(), TSCreate()
1990: @*/
1991: PetscErrorCode TSViewFromOptions(TS A,PetscObject obj,const char name[])
1992: {
1997: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1998: return(0);
1999: }
2001: /*@C
2002: TSView - Prints the TS data structure.
2004: Collective on TS
2006: Input Parameters:
2007: + ts - the TS context obtained from TSCreate()
2008: - viewer - visualization context
2010: Options Database Key:
2011: . -ts_view - calls TSView() at end of TSStep()
2013: Notes:
2014: The available visualization contexts include
2015: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
2016: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
2017: output where only the first processor opens
2018: the file. All other processors send their
2019: data to the first processor to print.
2021: The user can open an alternative visualization context with
2022: PetscViewerASCIIOpen() - output to a specified file.
2024: Level: beginner
2026: .seealso: PetscViewerASCIIOpen()
2027: @*/
2028: PetscErrorCode TSView(TS ts,PetscViewer viewer)
2029: {
2031: TSType type;
2032: PetscBool iascii,isstring,isundials,isbinary,isdraw;
2033: DMTS sdm;
2034: #if defined(PETSC_HAVE_SAWS)
2035: PetscBool issaws;
2036: #endif
2040: if (!viewer) {
2041: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts),&viewer);
2042: }
2046: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2047: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
2048: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2049: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
2050: #if defined(PETSC_HAVE_SAWS)
2051: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
2052: #endif
2053: if (iascii) {
2054: PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer);
2055: if (ts->ops->view) {
2056: PetscViewerASCIIPushTab(viewer);
2057: (*ts->ops->view)(ts,viewer);
2058: PetscViewerASCIIPopTab(viewer);
2059: }
2060: if (ts->max_steps < PETSC_MAX_INT) {
2061: PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);
2062: }
2063: if (ts->max_time < PETSC_MAX_REAL) {
2064: PetscViewerASCIIPrintf(viewer," maximum time=%g\n",(double)ts->max_time);
2065: }
2066: if (ts->usessnes) {
2067: PetscBool lin;
2068: if (ts->problem_type == TS_NONLINEAR) {
2069: PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->snes_its);
2070: }
2071: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->ksp_its);
2072: PetscObjectTypeCompareAny((PetscObject)ts->snes,&lin,SNESKSPONLY,SNESKSPTRANSPOSEONLY,"");
2073: PetscViewerASCIIPrintf(viewer," total number of %slinear solve failures=%D\n",lin ? "" : "non",ts->num_snes_failures);
2074: }
2075: PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);
2076: if (ts->vrtol) {
2077: PetscViewerASCIIPrintf(viewer," using vector of relative error tolerances, ");
2078: } else {
2079: PetscViewerASCIIPrintf(viewer," using relative error tolerance of %g, ",(double)ts->rtol);
2080: }
2081: if (ts->vatol) {
2082: PetscViewerASCIIPrintf(viewer," using vector of absolute error tolerances\n");
2083: } else {
2084: PetscViewerASCIIPrintf(viewer," using absolute error tolerance of %g\n",(double)ts->atol);
2085: }
2086: PetscViewerASCIIPushTab(viewer);
2087: TSAdaptView(ts->adapt,viewer);
2088: PetscViewerASCIIPopTab(viewer);
2089: } else if (isstring) {
2090: TSGetType(ts,&type);
2091: PetscViewerStringSPrintf(viewer," TSType: %-7.7s",type);
2092: if (ts->ops->view) {(*ts->ops->view)(ts,viewer);}
2093: } else if (isbinary) {
2094: PetscInt classid = TS_FILE_CLASSID;
2095: MPI_Comm comm;
2096: PetscMPIInt rank;
2097: char type[256];
2099: PetscObjectGetComm((PetscObject)ts,&comm);
2100: MPI_Comm_rank(comm,&rank);
2101: if (!rank) {
2102: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);
2103: PetscStrncpy(type,((PetscObject)ts)->type_name,256);
2104: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);
2105: }
2106: if (ts->ops->view) {
2107: (*ts->ops->view)(ts,viewer);
2108: }
2109: if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
2110: DMView(ts->dm,viewer);
2111: VecView(ts->vec_sol,viewer);
2112: DMGetDMTS(ts->dm,&sdm);
2113: DMTSView(sdm,viewer);
2114: } else if (isdraw) {
2115: PetscDraw draw;
2116: char str[36];
2117: PetscReal x,y,bottom,h;
2119: PetscViewerDrawGetDraw(viewer,0,&draw);
2120: PetscDrawGetCurrentPoint(draw,&x,&y);
2121: PetscStrcpy(str,"TS: ");
2122: PetscStrcat(str,((PetscObject)ts)->type_name);
2123: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,NULL,&h);
2124: bottom = y - h;
2125: PetscDrawPushCurrentPoint(draw,x,bottom);
2126: if (ts->ops->view) {
2127: (*ts->ops->view)(ts,viewer);
2128: }
2129: if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
2130: if (ts->snes) {SNESView(ts->snes,viewer);}
2131: PetscDrawPopCurrentPoint(draw);
2132: #if defined(PETSC_HAVE_SAWS)
2133: } else if (issaws) {
2134: PetscMPIInt rank;
2135: const char *name;
2137: PetscObjectGetName((PetscObject)ts,&name);
2138: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
2139: if (!((PetscObject)ts)->amsmem && !rank) {
2140: char dir[1024];
2142: PetscObjectViewSAWs((PetscObject)ts,viewer);
2143: PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time_step",name);
2144: PetscStackCallSAWs(SAWs_Register,(dir,&ts->steps,1,SAWs_READ,SAWs_INT));
2145: PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time",name);
2146: PetscStackCallSAWs(SAWs_Register,(dir,&ts->ptime,1,SAWs_READ,SAWs_DOUBLE));
2147: }
2148: if (ts->ops->view) {
2149: (*ts->ops->view)(ts,viewer);
2150: }
2151: #endif
2152: }
2153: if (ts->snes && ts->usessnes) {
2154: PetscViewerASCIIPushTab(viewer);
2155: SNESView(ts->snes,viewer);
2156: PetscViewerASCIIPopTab(viewer);
2157: }
2158: DMGetDMTS(ts->dm,&sdm);
2159: DMTSView(sdm,viewer);
2161: PetscViewerASCIIPushTab(viewer);
2162: PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);
2163: PetscViewerASCIIPopTab(viewer);
2164: return(0);
2165: }
2167: /*@
2168: TSSetApplicationContext - Sets an optional user-defined context for
2169: the timesteppers.
2171: Logically Collective on TS
2173: Input Parameters:
2174: + ts - the TS context obtained from TSCreate()
2175: - usrP - optional user context
2177: Fortran Notes:
2178: To use this from Fortran you must write a Fortran interface definition for this
2179: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
2181: Level: intermediate
2183: .seealso: TSGetApplicationContext()
2184: @*/
2185: PetscErrorCode TSSetApplicationContext(TS ts,void *usrP)
2186: {
2189: ts->user = usrP;
2190: return(0);
2191: }
2193: /*@
2194: TSGetApplicationContext - Gets the user-defined context for the
2195: timestepper.
2197: Not Collective
2199: Input Parameter:
2200: . ts - the TS context obtained from TSCreate()
2202: Output Parameter:
2203: . usrP - user context
2205: Fortran Notes:
2206: To use this from Fortran you must write a Fortran interface definition for this
2207: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
2209: Level: intermediate
2211: .seealso: TSSetApplicationContext()
2212: @*/
2213: PetscErrorCode TSGetApplicationContext(TS ts,void *usrP)
2214: {
2217: *(void**)usrP = ts->user;
2218: return(0);
2219: }
2221: /*@
2222: TSGetStepNumber - Gets the number of steps completed.
2224: Not Collective
2226: Input Parameter:
2227: . ts - the TS context obtained from TSCreate()
2229: Output Parameter:
2230: . steps - number of steps completed so far
2232: Level: intermediate
2234: .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSSetPostStep()
2235: @*/
2236: PetscErrorCode TSGetStepNumber(TS ts,PetscInt *steps)
2237: {
2241: *steps = ts->steps;
2242: return(0);
2243: }
2245: /*@
2246: TSSetStepNumber - Sets the number of steps completed.
2248: Logically Collective on TS
2250: Input Parameters:
2251: + ts - the TS context
2252: - steps - number of steps completed so far
2254: Notes:
2255: For most uses of the TS solvers the user need not explicitly call
2256: TSSetStepNumber(), as the step counter is appropriately updated in
2257: TSSolve()/TSStep()/TSRollBack(). Power users may call this routine to
2258: reinitialize timestepping by setting the step counter to zero (and time
2259: to the initial time) to solve a similar problem with different initial
2260: conditions or parameters. Other possible use case is to continue
2261: timestepping from a previously interrupted run in such a way that TS
2262: monitors will be called with a initial nonzero step counter.
2264: Level: advanced
2266: .seealso: TSGetStepNumber(), TSSetTime(), TSSetTimeStep(), TSSetSolution()
2267: @*/
2268: PetscErrorCode TSSetStepNumber(TS ts,PetscInt steps)
2269: {
2273: if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Step number must be non-negative");
2274: ts->steps = steps;
2275: return(0);
2276: }
2278: /*@
2279: TSSetTimeStep - Allows one to reset the timestep at any time,
2280: useful for simple pseudo-timestepping codes.
2282: Logically Collective on TS
2284: Input Parameters:
2285: + ts - the TS context obtained from TSCreate()
2286: - time_step - the size of the timestep
2288: Level: intermediate
2290: .seealso: TSGetTimeStep(), TSSetTime()
2292: @*/
2293: PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step)
2294: {
2298: ts->time_step = time_step;
2299: return(0);
2300: }
2302: /*@
2303: TSSetExactFinalTime - Determines whether to adapt the final time step to
2304: match the exact final time, interpolate solution to the exact final time,
2305: or just return at the final time TS computed.
2307: Logically Collective on TS
2309: Input Parameter:
2310: + ts - the time-step context
2311: - eftopt - exact final time option
2313: $ TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded
2314: $ TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2315: $ TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2317: Options Database:
2318: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime
2320: Warning: If you use the option TS_EXACTFINALTIME_STEPOVER the solution may be at a very different time
2321: then the final time you selected.
2323: Level: beginner
2325: .seealso: TSExactFinalTimeOption, TSGetExactFinalTime()
2326: @*/
2327: PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt)
2328: {
2332: ts->exact_final_time = eftopt;
2333: return(0);
2334: }
2336: /*@
2337: TSGetExactFinalTime - Gets the exact final time option.
2339: Not Collective
2341: Input Parameter:
2342: . ts - the TS context
2344: Output Parameter:
2345: . eftopt - exact final time option
2347: Level: beginner
2349: .seealso: TSExactFinalTimeOption, TSSetExactFinalTime()
2350: @*/
2351: PetscErrorCode TSGetExactFinalTime(TS ts,TSExactFinalTimeOption *eftopt)
2352: {
2356: *eftopt = ts->exact_final_time;
2357: return(0);
2358: }
2360: /*@
2361: TSGetTimeStep - Gets the current timestep size.
2363: Not Collective
2365: Input Parameter:
2366: . ts - the TS context obtained from TSCreate()
2368: Output Parameter:
2369: . dt - the current timestep size
2371: Level: intermediate
2373: .seealso: TSSetTimeStep(), TSGetTime()
2375: @*/
2376: PetscErrorCode TSGetTimeStep(TS ts,PetscReal *dt)
2377: {
2381: *dt = ts->time_step;
2382: return(0);
2383: }
2385: /*@
2386: TSGetSolution - Returns the solution at the present timestep. It
2387: is valid to call this routine inside the function that you are evaluating
2388: in order to move to the new timestep. This vector not changed until
2389: the solution at the next timestep has been calculated.
2391: Not Collective, but Vec returned is parallel if TS is parallel
2393: Input Parameter:
2394: . ts - the TS context obtained from TSCreate()
2396: Output Parameter:
2397: . v - the vector containing the solution
2399: Note: If you used TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); this does not return the solution at the requested
2400: final time. It returns the solution at the next timestep.
2402: Level: intermediate
2404: .seealso: TSGetTimeStep(), TSGetTime(), TSGetSolveTime(), TSGetSolutionComponents(), TSSetSolutionFunction()
2406: @*/
2407: PetscErrorCode TSGetSolution(TS ts,Vec *v)
2408: {
2412: *v = ts->vec_sol;
2413: return(0);
2414: }
2416: /*@
2417: TSGetSolutionComponents - Returns any solution components at the present
2418: timestep, if available for the time integration method being used.
2419: Solution components are quantities that share the same size and
2420: structure as the solution vector.
2422: Not Collective, but Vec returned is parallel if TS is parallel
2424: Parameters :
2425: + ts - the TS context obtained from TSCreate() (input parameter).
2426: . n - If v is PETSC_NULL, then the number of solution components is
2427: returned through n, else the n-th solution component is
2428: returned in v.
2429: - v - the vector containing the n-th solution component
2430: (may be PETSC_NULL to use this function to find out
2431: the number of solutions components).
2433: Level: advanced
2435: .seealso: TSGetSolution()
2437: @*/
2438: PetscErrorCode TSGetSolutionComponents(TS ts,PetscInt *n,Vec *v)
2439: {
2444: if (!ts->ops->getsolutioncomponents) *n = 0;
2445: else {
2446: (*ts->ops->getsolutioncomponents)(ts,n,v);
2447: }
2448: return(0);
2449: }
2451: /*@
2452: TSGetAuxSolution - Returns an auxiliary solution at the present
2453: timestep, if available for the time integration method being used.
2455: Not Collective, but Vec returned is parallel if TS is parallel
2457: Parameters :
2458: + ts - the TS context obtained from TSCreate() (input parameter).
2459: - v - the vector containing the auxiliary solution
2461: Level: intermediate
2463: .seealso: TSGetSolution()
2465: @*/
2466: PetscErrorCode TSGetAuxSolution(TS ts,Vec *v)
2467: {
2472: if (ts->ops->getauxsolution) {
2473: (*ts->ops->getauxsolution)(ts,v);
2474: } else {
2475: VecZeroEntries(*v);
2476: }
2477: return(0);
2478: }
2480: /*@
2481: TSGetTimeError - Returns the estimated error vector, if the chosen
2482: TSType has an error estimation functionality.
2484: Not Collective, but Vec returned is parallel if TS is parallel
2486: Note: MUST call after TSSetUp()
2488: Parameters :
2489: + ts - the TS context obtained from TSCreate() (input parameter).
2490: . n - current estimate (n=0) or previous one (n=-1)
2491: - v - the vector containing the error (same size as the solution).
2493: Level: intermediate
2495: .seealso: TSGetSolution(), TSSetTimeError()
2497: @*/
2498: PetscErrorCode TSGetTimeError(TS ts,PetscInt n,Vec *v)
2499: {
2504: if (ts->ops->gettimeerror) {
2505: (*ts->ops->gettimeerror)(ts,n,v);
2506: } else {
2507: VecZeroEntries(*v);
2508: }
2509: return(0);
2510: }
2512: /*@
2513: TSSetTimeError - Sets the estimated error vector, if the chosen
2514: TSType has an error estimation functionality. This can be used
2515: to restart such a time integrator with a given error vector.
2517: Not Collective, but Vec returned is parallel if TS is parallel
2519: Parameters :
2520: + ts - the TS context obtained from TSCreate() (input parameter).
2521: - v - the vector containing the error (same size as the solution).
2523: Level: intermediate
2525: .seealso: TSSetSolution(), TSGetTimeError)
2527: @*/
2528: PetscErrorCode TSSetTimeError(TS ts,Vec v)
2529: {
2534: if (!ts->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetUp() first");
2535: if (ts->ops->settimeerror) {
2536: (*ts->ops->settimeerror)(ts,v);
2537: }
2538: return(0);
2539: }
2541: /* ----- Routines to initialize and destroy a timestepper ---- */
2542: /*@
2543: TSSetProblemType - Sets the type of problem to be solved.
2545: Not collective
2547: Input Parameters:
2548: + ts - The TS
2549: - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2550: .vb
2551: U_t - A U = 0 (linear)
2552: U_t - A(t) U = 0 (linear)
2553: F(t,U,U_t) = 0 (nonlinear)
2554: .ve
2556: Level: beginner
2558: .seealso: TSSetUp(), TSProblemType, TS
2559: @*/
2560: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2561: {
2566: ts->problem_type = type;
2567: if (type == TS_LINEAR) {
2568: SNES snes;
2569: TSGetSNES(ts,&snes);
2570: SNESSetType(snes,SNESKSPONLY);
2571: }
2572: return(0);
2573: }
2575: /*@C
2576: TSGetProblemType - Gets the type of problem to be solved.
2578: Not collective
2580: Input Parameter:
2581: . ts - The TS
2583: Output Parameter:
2584: . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2585: .vb
2586: M U_t = A U
2587: M(t) U_t = A(t) U
2588: F(t,U,U_t)
2589: .ve
2591: Level: beginner
2593: .seealso: TSSetUp(), TSProblemType, TS
2594: @*/
2595: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2596: {
2600: *type = ts->problem_type;
2601: return(0);
2602: }
2604: /*@
2605: TSSetUp - Sets up the internal data structures for the later use
2606: of a timestepper.
2608: Collective on TS
2610: Input Parameter:
2611: . ts - the TS context obtained from TSCreate()
2613: Notes:
2614: For basic use of the TS solvers the user need not explicitly call
2615: TSSetUp(), since these actions will automatically occur during
2616: the call to TSStep() or TSSolve(). However, if one wishes to control this
2617: phase separately, TSSetUp() should be called after TSCreate()
2618: and optional routines of the form TSSetXXX(), but before TSStep() and TSSolve().
2620: Level: advanced
2622: .seealso: TSCreate(), TSStep(), TSDestroy(), TSSolve()
2623: @*/
2624: PetscErrorCode TSSetUp(TS ts)
2625: {
2627: DM dm;
2628: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2629: PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*);
2630: TSIFunction ifun;
2631: TSIJacobian ijac;
2632: TSI2Jacobian i2jac;
2633: TSRHSJacobian rhsjac;
2634: PetscBool isnone;
2638: if (ts->setupcalled) return(0);
2640: if (!((PetscObject)ts)->type_name) {
2641: TSGetIFunction(ts,NULL,&ifun,NULL);
2642: TSSetType(ts,ifun ? TSBEULER : TSEULER);
2643: }
2645: if (!ts->vec_sol) {
2646: if (ts->dm) {
2647: DMCreateGlobalVector(ts->dm,&ts->vec_sol);
2648: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first");
2649: }
2651: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2652: PetscObjectReference((PetscObject)ts->Jacprhs);
2653: ts->Jacp = ts->Jacprhs;
2654: }
2656: if (ts->quadraturets) {
2657: TSSetUp(ts->quadraturets);
2658: VecDestroy(&ts->vec_costintegrand);
2659: VecDuplicate(ts->quadraturets->vec_sol,&ts->vec_costintegrand);
2660: }
2662: TSGetRHSJacobian(ts,NULL,NULL,&rhsjac,NULL);
2663: if (ts->rhsjacobian.reuse && rhsjac == TSComputeRHSJacobianConstant) {
2664: Mat Amat,Pmat;
2665: SNES snes;
2666: TSGetSNES(ts,&snes);
2667: SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);
2668: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2669: * have displaced the RHS matrix */
2670: if (Amat && Amat == ts->Arhs) {
2671: /* 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 */
2672: MatDuplicate(ts->Arhs,MAT_COPY_VALUES,&Amat);
2673: SNESSetJacobian(snes,Amat,NULL,NULL,NULL);
2674: MatDestroy(&Amat);
2675: }
2676: if (Pmat && Pmat == ts->Brhs) {
2677: MatDuplicate(ts->Brhs,MAT_COPY_VALUES,&Pmat);
2678: SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);
2679: MatDestroy(&Pmat);
2680: }
2681: }
2683: TSGetAdapt(ts,&ts->adapt);
2684: TSAdaptSetDefaultType(ts->adapt,ts->default_adapt_type);
2686: if (ts->ops->setup) {
2687: (*ts->ops->setup)(ts);
2688: }
2690: /* Attempt to check/preset a default value for the exact final time option */
2691: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&isnone);
2692: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
2693: ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2695: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2696: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2697: */
2698: TSGetDM(ts,&dm);
2699: DMSNESGetFunction(dm,&func,NULL);
2700: if (!func) {
2701: DMSNESSetFunction(dm,SNESTSFormFunction,ts);
2702: }
2703: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2704: Otherwise, the SNES will use coloring internally to form the Jacobian.
2705: */
2706: DMSNESGetJacobian(dm,&jac,NULL);
2707: DMTSGetIJacobian(dm,&ijac,NULL);
2708: DMTSGetI2Jacobian(dm,&i2jac,NULL);
2709: DMTSGetRHSJacobian(dm,&rhsjac,NULL);
2710: if (!jac && (ijac || i2jac || rhsjac)) {
2711: DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);
2712: }
2714: /* if time integration scheme has a starting method, call it */
2715: if (ts->ops->startingmethod) {
2716: (*ts->ops->startingmethod)(ts);
2717: }
2719: ts->setupcalled = PETSC_TRUE;
2720: return(0);
2721: }
2723: /*@
2724: TSReset - Resets a TS context and removes any allocated Vecs and Mats.
2726: Collective on TS
2728: Input Parameter:
2729: . ts - the TS context obtained from TSCreate()
2731: Level: beginner
2733: .seealso: TSCreate(), TSSetup(), TSDestroy()
2734: @*/
2735: PetscErrorCode TSReset(TS ts)
2736: {
2737: TS_RHSSplitLink ilink = ts->tsrhssplit,next;
2738: PetscErrorCode ierr;
2743: if (ts->ops->reset) {
2744: (*ts->ops->reset)(ts);
2745: }
2746: if (ts->snes) {SNESReset(ts->snes);}
2747: if (ts->adapt) {TSAdaptReset(ts->adapt);}
2749: MatDestroy(&ts->Arhs);
2750: MatDestroy(&ts->Brhs);
2751: VecDestroy(&ts->Frhs);
2752: VecDestroy(&ts->vec_sol);
2753: VecDestroy(&ts->vec_dot);
2754: VecDestroy(&ts->vatol);
2755: VecDestroy(&ts->vrtol);
2756: VecDestroyVecs(ts->nwork,&ts->work);
2758: MatDestroy(&ts->Jacprhs);
2759: MatDestroy(&ts->Jacp);
2760: if (ts->forward_solve) {
2761: TSForwardReset(ts);
2762: }
2763: if (ts->quadraturets) {
2764: TSReset(ts->quadraturets);
2765: VecDestroy(&ts->vec_costintegrand);
2766: }
2767: while (ilink) {
2768: next = ilink->next;
2769: TSDestroy(&ilink->ts);
2770: PetscFree(ilink->splitname);
2771: ISDestroy(&ilink->is);
2772: PetscFree(ilink);
2773: ilink = next;
2774: }
2775: ts->num_rhs_splits = 0;
2776: ts->setupcalled = PETSC_FALSE;
2777: return(0);
2778: }
2780: /*@
2781: TSDestroy - Destroys the timestepper context that was created
2782: with TSCreate().
2784: Collective on TS
2786: Input Parameter:
2787: . ts - the TS context obtained from TSCreate()
2789: Level: beginner
2791: .seealso: TSCreate(), TSSetUp(), TSSolve()
2792: @*/
2793: PetscErrorCode TSDestroy(TS *ts)
2794: {
2798: if (!*ts) return(0);
2800: if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; return(0);}
2802: TSReset(*ts);
2803: TSAdjointReset(*ts);
2804: if ((*ts)->forward_solve) {
2805: TSForwardReset(*ts);
2806: }
2807: /* if memory was published with SAWs then destroy it */
2808: PetscObjectSAWsViewOff((PetscObject)*ts);
2809: if ((*ts)->ops->destroy) {(*(*ts)->ops->destroy)((*ts));}
2811: TSTrajectoryDestroy(&(*ts)->trajectory);
2813: TSAdaptDestroy(&(*ts)->adapt);
2814: TSEventDestroy(&(*ts)->event);
2816: SNESDestroy(&(*ts)->snes);
2817: DMDestroy(&(*ts)->dm);
2818: TSMonitorCancel((*ts));
2819: TSAdjointMonitorCancel((*ts));
2821: TSDestroy(&(*ts)->quadraturets);
2822: PetscHeaderDestroy(ts);
2823: return(0);
2824: }
2826: /*@
2827: TSGetSNES - Returns the SNES (nonlinear solver) associated with
2828: a TS (timestepper) context. Valid only for nonlinear problems.
2830: Not Collective, but SNES is parallel if TS is parallel
2832: Input Parameter:
2833: . ts - the TS context obtained from TSCreate()
2835: Output Parameter:
2836: . snes - the nonlinear solver context
2838: Notes:
2839: The user can then directly manipulate the SNES context to set various
2840: options, etc. Likewise, the user can then extract and manipulate the
2841: KSP, KSP, and PC contexts as well.
2843: TSGetSNES() does not work for integrators that do not use SNES; in
2844: this case TSGetSNES() returns NULL in snes.
2846: Level: beginner
2848: @*/
2849: PetscErrorCode TSGetSNES(TS ts,SNES *snes)
2850: {
2856: if (!ts->snes) {
2857: SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);
2858: PetscObjectSetOptions((PetscObject)ts->snes,((PetscObject)ts)->options);
2859: SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
2860: PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);
2861: PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);
2862: if (ts->dm) {SNESSetDM(ts->snes,ts->dm);}
2863: if (ts->problem_type == TS_LINEAR) {
2864: SNESSetType(ts->snes,SNESKSPONLY);
2865: }
2866: }
2867: *snes = ts->snes;
2868: return(0);
2869: }
2871: /*@
2872: TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context
2874: Collective
2876: Input Parameter:
2877: + ts - the TS context obtained from TSCreate()
2878: - snes - the nonlinear solver context
2880: Notes:
2881: Most users should have the TS created by calling TSGetSNES()
2883: Level: developer
2885: @*/
2886: PetscErrorCode TSSetSNES(TS ts,SNES snes)
2887: {
2889: PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*);
2894: PetscObjectReference((PetscObject)snes);
2895: SNESDestroy(&ts->snes);
2897: ts->snes = snes;
2899: SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
2900: SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);
2901: if (func == SNESTSFormJacobian) {
2902: SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);
2903: }
2904: return(0);
2905: }
2907: /*@
2908: TSGetKSP - Returns the KSP (linear solver) associated with
2909: a TS (timestepper) context.
2911: Not Collective, but KSP is parallel if TS is parallel
2913: Input Parameter:
2914: . ts - the TS context obtained from TSCreate()
2916: Output Parameter:
2917: . ksp - the nonlinear solver context
2919: Notes:
2920: The user can then directly manipulate the KSP context to set various
2921: options, etc. Likewise, the user can then extract and manipulate the
2922: KSP and PC contexts as well.
2924: TSGetKSP() does not work for integrators that do not use KSP;
2925: in this case TSGetKSP() returns NULL in ksp.
2927: Level: beginner
2929: @*/
2930: PetscErrorCode TSGetKSP(TS ts,KSP *ksp)
2931: {
2933: SNES snes;
2938: if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first");
2939: if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()");
2940: TSGetSNES(ts,&snes);
2941: SNESGetKSP(snes,ksp);
2942: return(0);
2943: }
2945: /* ----------- Routines to set solver parameters ---------- */
2947: /*@
2948: TSSetMaxSteps - Sets the maximum number of steps to use.
2950: Logically Collective on TS
2952: Input Parameters:
2953: + ts - the TS context obtained from TSCreate()
2954: - maxsteps - maximum number of steps to use
2956: Options Database Keys:
2957: . -ts_max_steps <maxsteps> - Sets maxsteps
2959: Notes:
2960: The default maximum number of steps is 5000
2962: Level: intermediate
2964: .seealso: TSGetMaxSteps(), TSSetMaxTime(), TSSetExactFinalTime()
2965: @*/
2966: PetscErrorCode TSSetMaxSteps(TS ts,PetscInt maxsteps)
2967: {
2971: if (maxsteps < 0 ) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of steps must be non-negative");
2972: ts->max_steps = maxsteps;
2973: return(0);
2974: }
2976: /*@
2977: TSGetMaxSteps - Gets the maximum number of steps to use.
2979: Not Collective
2981: Input Parameters:
2982: . ts - the TS context obtained from TSCreate()
2984: Output Parameter:
2985: . maxsteps - maximum number of steps to use
2987: Level: advanced
2989: .seealso: TSSetMaxSteps(), TSGetMaxTime(), TSSetMaxTime()
2990: @*/
2991: PetscErrorCode TSGetMaxSteps(TS ts,PetscInt *maxsteps)
2992: {
2996: *maxsteps = ts->max_steps;
2997: return(0);
2998: }
3000: /*@
3001: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
3003: Logically Collective on TS
3005: Input Parameters:
3006: + ts - the TS context obtained from TSCreate()
3007: - maxtime - final time to step to
3009: Options Database Keys:
3010: . -ts_max_time <maxtime> - Sets maxtime
3012: Notes:
3013: The default maximum time is 5.0
3015: Level: intermediate
3017: .seealso: TSGetMaxTime(), TSSetMaxSteps(), TSSetExactFinalTime()
3018: @*/
3019: PetscErrorCode TSSetMaxTime(TS ts,PetscReal maxtime)
3020: {
3024: ts->max_time = maxtime;
3025: return(0);
3026: }
3028: /*@
3029: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
3031: Not Collective
3033: Input Parameters:
3034: . ts - the TS context obtained from TSCreate()
3036: Output Parameter:
3037: . maxtime - final time to step to
3039: Level: advanced
3041: .seealso: TSSetMaxTime(), TSGetMaxSteps(), TSSetMaxSteps()
3042: @*/
3043: PetscErrorCode TSGetMaxTime(TS ts,PetscReal *maxtime)
3044: {
3048: *maxtime = ts->max_time;
3049: return(0);
3050: }
3052: /*@
3053: TSSetInitialTimeStep - Deprecated, use TSSetTime() and TSSetTimeStep().
3055: Level: deprecated
3057: @*/
3058: PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step)
3059: {
3063: TSSetTime(ts,initial_time);
3064: TSSetTimeStep(ts,time_step);
3065: return(0);
3066: }
3068: /*@
3069: TSGetDuration - Deprecated, use TSGetMaxSteps() and TSGetMaxTime().
3071: Level: deprecated
3073: @*/
3074: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
3075: {
3078: if (maxsteps) {
3080: *maxsteps = ts->max_steps;
3081: }
3082: if (maxtime) {
3084: *maxtime = ts->max_time;
3085: }
3086: return(0);
3087: }
3089: /*@
3090: TSSetDuration - Deprecated, use TSSetMaxSteps() and TSSetMaxTime().
3092: Level: deprecated
3094: @*/
3095: PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime)
3096: {
3101: if (maxsteps >= 0) ts->max_steps = maxsteps;
3102: if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
3103: return(0);
3104: }
3106: /*@
3107: TSGetTimeStepNumber - Deprecated, use TSGetStepNumber().
3109: Level: deprecated
3111: @*/
3112: PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *steps) { return TSGetStepNumber(ts,steps); }
3114: /*@
3115: TSGetTotalSteps - Deprecated, use TSGetStepNumber().
3117: Level: deprecated
3119: @*/
3120: PetscErrorCode TSGetTotalSteps(TS ts,PetscInt *steps) { return TSGetStepNumber(ts,steps); }
3122: /*@
3123: TSSetSolution - Sets the initial solution vector
3124: for use by the TS routines.
3126: Logically Collective on TS
3128: Input Parameters:
3129: + ts - the TS context obtained from TSCreate()
3130: - u - the solution vector
3132: Level: beginner
3134: .seealso: TSSetSolutionFunction(), TSGetSolution(), TSCreate()
3135: @*/
3136: PetscErrorCode TSSetSolution(TS ts,Vec u)
3137: {
3139: DM dm;
3144: PetscObjectReference((PetscObject)u);
3145: VecDestroy(&ts->vec_sol);
3146: ts->vec_sol = u;
3148: TSGetDM(ts,&dm);
3149: DMShellSetGlobalVector(dm,u);
3150: return(0);
3151: }
3153: /*@C
3154: TSSetPreStep - Sets the general-purpose function
3155: called once at the beginning of each time step.
3157: Logically Collective on TS
3159: Input Parameters:
3160: + ts - The TS context obtained from TSCreate()
3161: - func - The function
3163: Calling sequence of func:
3164: . PetscErrorCode func (TS ts);
3166: Level: intermediate
3168: .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep(), TSRestartStep()
3169: @*/
3170: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
3171: {
3174: ts->prestep = func;
3175: return(0);
3176: }
3178: /*@
3179: TSPreStep - Runs the user-defined pre-step function.
3181: Collective on TS
3183: Input Parameters:
3184: . ts - The TS context obtained from TSCreate()
3186: Notes:
3187: TSPreStep() is typically used within time stepping implementations,
3188: so most users would not generally call this routine themselves.
3190: Level: developer
3192: .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep()
3193: @*/
3194: PetscErrorCode TSPreStep(TS ts)
3195: {
3200: if (ts->prestep) {
3201: Vec U;
3202: PetscObjectState sprev,spost;
3204: TSGetSolution(ts,&U);
3205: PetscObjectStateGet((PetscObject)U,&sprev);
3206: PetscStackCallStandard((*ts->prestep),(ts));
3207: PetscObjectStateGet((PetscObject)U,&spost);
3208: if (sprev != spost) {TSRestartStep(ts);}
3209: }
3210: return(0);
3211: }
3213: /*@C
3214: TSSetPreStage - Sets the general-purpose function
3215: called once at the beginning of each stage.
3217: Logically Collective on TS
3219: Input Parameters:
3220: + ts - The TS context obtained from TSCreate()
3221: - func - The function
3223: Calling sequence of func:
3224: . PetscErrorCode func(TS ts, PetscReal stagetime);
3226: Level: intermediate
3228: Note:
3229: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3230: The time step number being computed can be queried using TSGetStepNumber() and the total size of the step being
3231: attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime().
3233: .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3234: @*/
3235: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal))
3236: {
3239: ts->prestage = func;
3240: return(0);
3241: }
3243: /*@C
3244: TSSetPostStage - Sets the general-purpose function
3245: called once at the end of each stage.
3247: Logically Collective on TS
3249: Input Parameters:
3250: + ts - The TS context obtained from TSCreate()
3251: - func - The function
3253: Calling sequence of func:
3254: . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y);
3256: Level: intermediate
3258: Note:
3259: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3260: The time step number being computed can be queried using TSGetStepNumber() and the total size of the step being
3261: attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime().
3263: .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3264: @*/
3265: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*))
3266: {
3269: ts->poststage = func;
3270: return(0);
3271: }
3273: /*@C
3274: TSSetPostEvaluate - Sets the general-purpose function
3275: called once at the end of each step evaluation.
3277: Logically Collective on TS
3279: Input Parameters:
3280: + ts - The TS context obtained from TSCreate()
3281: - func - The function
3283: Calling sequence of func:
3284: . PetscErrorCode func(TS ts);
3286: Level: intermediate
3288: Note:
3289: Semantically, TSSetPostEvaluate() differs from TSSetPostStep() since the function it sets is called before event-handling
3290: thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, TSPostStep()
3291: may be passed a different solution, possibly changed by the event handler. TSPostEvaluate() is called after the next step
3292: solution is evaluated allowing to modify it, if need be. The solution can be obtained with TSGetSolution(), the time step
3293: with TSGetTimeStep(), and the time at the start of the step is available via TSGetTime()
3295: .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3296: @*/
3297: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS))
3298: {
3301: ts->postevaluate = func;
3302: return(0);
3303: }
3305: /*@
3306: TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage()
3308: Collective on TS
3310: Input Parameters:
3311: . ts - The TS context obtained from TSCreate()
3312: stagetime - The absolute time of the current stage
3314: Notes:
3315: TSPreStage() is typically used within time stepping implementations,
3316: most users would not generally call this routine themselves.
3318: Level: developer
3320: .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep()
3321: @*/
3322: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3323: {
3326: if (ts->prestage) {
3327: PetscStackCallStandard((*ts->prestage),(ts,stagetime));
3328: }
3329: return(0);
3330: }
3332: /*@
3333: TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage()
3335: Collective on TS
3337: Input Parameters:
3338: . ts - The TS context obtained from TSCreate()
3339: stagetime - The absolute time of the current stage
3340: stageindex - Stage number
3341: Y - Array of vectors (of size = total number
3342: of stages) with the stage solutions
3344: Notes:
3345: TSPostStage() is typically used within time stepping implementations,
3346: most users would not generally call this routine themselves.
3348: Level: developer
3350: .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep()
3351: @*/
3352: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3353: {
3356: if (ts->poststage) {
3357: PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y));
3358: }
3359: return(0);
3360: }
3362: /*@
3363: TSPostEvaluate - Runs the user-defined post-evaluate function set using TSSetPostEvaluate()
3365: Collective on TS
3367: Input Parameters:
3368: . ts - The TS context obtained from TSCreate()
3370: Notes:
3371: TSPostEvaluate() is typically used within time stepping implementations,
3372: most users would not generally call this routine themselves.
3374: Level: developer
3376: .seealso: TSSetPostEvaluate(), TSSetPreStep(), TSPreStep(), TSPostStep()
3377: @*/
3378: PetscErrorCode TSPostEvaluate(TS ts)
3379: {
3384: if (ts->postevaluate) {
3385: Vec U;
3386: PetscObjectState sprev,spost;
3388: TSGetSolution(ts,&U);
3389: PetscObjectStateGet((PetscObject)U,&sprev);
3390: PetscStackCallStandard((*ts->postevaluate),(ts));
3391: PetscObjectStateGet((PetscObject)U,&spost);
3392: if (sprev != spost) {TSRestartStep(ts);}
3393: }
3394: return(0);
3395: }
3397: /*@C
3398: TSSetPostStep - Sets the general-purpose function
3399: called once at the end of each time step.
3401: Logically Collective on TS
3403: Input Parameters:
3404: + ts - The TS context obtained from TSCreate()
3405: - func - The function
3407: Calling sequence of func:
3408: $ func (TS ts);
3410: Notes:
3411: The function set by TSSetPostStep() is called after each successful step. The solution vector X
3412: obtained by TSGetSolution() may be different than that computed at the step end if the event handler
3413: locates an event and TSPostEvent() modifies it. Use TSSetPostEvaluate() if an unmodified solution is needed instead.
3415: Level: intermediate
3417: .seealso: TSSetPreStep(), TSSetPreStage(), TSSetPostEvaluate(), TSGetTimeStep(), TSGetStepNumber(), TSGetTime(), TSRestartStep()
3418: @*/
3419: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3420: {
3423: ts->poststep = func;
3424: return(0);
3425: }
3427: /*@
3428: TSPostStep - Runs the user-defined post-step function.
3430: Collective on TS
3432: Input Parameters:
3433: . ts - The TS context obtained from TSCreate()
3435: Notes:
3436: TSPostStep() is typically used within time stepping implementations,
3437: so most users would not generally call this routine themselves.
3439: Level: developer
3441: @*/
3442: PetscErrorCode TSPostStep(TS ts)
3443: {
3448: if (ts->poststep) {
3449: Vec U;
3450: PetscObjectState sprev,spost;
3452: TSGetSolution(ts,&U);
3453: PetscObjectStateGet((PetscObject)U,&sprev);
3454: PetscStackCallStandard((*ts->poststep),(ts));
3455: PetscObjectStateGet((PetscObject)U,&spost);
3456: if (sprev != spost) {TSRestartStep(ts);}
3457: }
3458: return(0);
3459: }
3461: /* ------------ Routines to set performance monitoring options ----------- */
3463: /*@C
3464: TSMonitorSet - Sets an ADDITIONAL function that is to be used at every
3465: timestep to display the iteration's progress.
3467: Logically Collective on TS
3469: Input Parameters:
3470: + ts - the TS context obtained from TSCreate()
3471: . monitor - monitoring routine
3472: . mctx - [optional] user-defined context for private data for the
3473: monitor routine (use NULL if no context is desired)
3474: - monitordestroy - [optional] routine that frees monitor context
3475: (may be NULL)
3477: Calling sequence of monitor:
3478: $ PetscErrorCode monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
3480: + ts - the TS context
3481: . 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)
3482: . time - current time
3483: . u - current iterate
3484: - mctx - [optional] monitoring context
3486: Notes:
3487: This routine adds an additional monitor to the list of monitors that
3488: already has been loaded.
3490: Fortran Notes:
3491: Only a single monitor function can be set for each TS object
3493: Level: intermediate
3495: .seealso: TSMonitorDefault(), TSMonitorCancel()
3496: @*/
3497: PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**))
3498: {
3500: PetscInt i;
3501: PetscBool identical;
3505: for (i=0; i<ts->numbermonitors;i++) {
3506: PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,mdestroy,(PetscErrorCode (*)(void))ts->monitor[i],ts->monitorcontext[i],ts->monitordestroy[i],&identical);
3507: if (identical) return(0);
3508: }
3509: if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3510: ts->monitor[ts->numbermonitors] = monitor;
3511: ts->monitordestroy[ts->numbermonitors] = mdestroy;
3512: ts->monitorcontext[ts->numbermonitors++] = (void*)mctx;
3513: return(0);
3514: }
3516: /*@C
3517: TSMonitorCancel - Clears all the monitors that have been set on a time-step object.
3519: Logically Collective on TS
3521: Input Parameters:
3522: . ts - the TS context obtained from TSCreate()
3524: Notes:
3525: There is no way to remove a single, specific monitor.
3527: Level: intermediate
3529: .seealso: TSMonitorDefault(), TSMonitorSet()
3530: @*/
3531: PetscErrorCode TSMonitorCancel(TS ts)
3532: {
3534: PetscInt i;
3538: for (i=0; i<ts->numbermonitors; i++) {
3539: if (ts->monitordestroy[i]) {
3540: (*ts->monitordestroy[i])(&ts->monitorcontext[i]);
3541: }
3542: }
3543: ts->numbermonitors = 0;
3544: return(0);
3545: }
3547: /*@C
3548: TSMonitorDefault - The Default monitor, prints the timestep and time for each step
3550: Level: intermediate
3552: .seealso: TSMonitorSet()
3553: @*/
3554: PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscViewerAndFormat *vf)
3555: {
3557: PetscViewer viewer = vf->viewer;
3558: PetscBool iascii,ibinary;
3562: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
3563: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
3564: PetscViewerPushFormat(viewer,vf->format);
3565: if (iascii) {
3566: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
3567: if (step == -1){ /* this indicates it is an interpolated solution */
3568: PetscViewerASCIIPrintf(viewer,"Interpolated solution at time %g between steps %D and %D\n",(double)ptime,ts->steps-1,ts->steps);
3569: } else {
3570: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
3571: }
3572: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
3573: } else if (ibinary) {
3574: PetscMPIInt rank;
3575: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);
3576: if (!rank) {
3577: PetscBool skipHeader;
3578: PetscInt classid = REAL_FILE_CLASSID;
3580: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
3581: if (!skipHeader) {
3582: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);
3583: }
3584: PetscRealView(1,&ptime,viewer);
3585: } else {
3586: PetscRealView(0,&ptime,viewer);
3587: }
3588: }
3589: PetscViewerPopFormat(viewer);
3590: return(0);
3591: }
3593: /*@C
3594: TSMonitorExtreme - Prints the extreme values of the solution at each timestep
3596: Level: intermediate
3598: .seealso: TSMonitorSet()
3599: @*/
3600: PetscErrorCode TSMonitorExtreme(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscViewerAndFormat *vf)
3601: {
3603: PetscViewer viewer = vf->viewer;
3604: PetscBool iascii;
3605: PetscReal max,min;
3610: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
3611: PetscViewerPushFormat(viewer,vf->format);
3612: if (iascii) {
3613: VecMax(v,NULL,&max);
3614: VecMin(v,NULL,&min);
3615: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
3616: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s max %g min %g\n",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)" : "",(double)max,(double)min);
3617: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
3618: }
3619: PetscViewerPopFormat(viewer);
3620: return(0);
3621: }
3623: /*@
3624: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3626: Collective on TS
3628: Input Argument:
3629: + ts - time stepping context
3630: - t - time to interpolate to
3632: Output Argument:
3633: . U - state at given time
3635: Level: intermediate
3637: Developer Notes:
3638: TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3640: .seealso: TSSetExactFinalTime(), TSSolve()
3641: @*/
3642: PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U)
3643: {
3649: 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);
3650: if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name);
3651: (*ts->ops->interpolate)(ts,t,U);
3652: return(0);
3653: }
3655: /*@
3656: TSStep - Steps one time step
3658: Collective on TS
3660: Input Parameter:
3661: . ts - the TS context obtained from TSCreate()
3663: Level: developer
3665: Notes:
3666: The public interface for the ODE/DAE solvers is TSSolve(), you should almost for sure be using that routine and not this routine.
3668: The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may
3669: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3671: This may over-step the final time provided in TSSetMaxTime() depending on the time-step used. TSSolve() interpolates to exactly the
3672: time provided in TSSetMaxTime(). One can use TSInterpolate() to determine an interpolated solution within the final timestep.
3674: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate()
3675: @*/
3676: PetscErrorCode TSStep(TS ts)
3677: {
3678: PetscErrorCode ierr;
3679: static PetscBool cite = PETSC_FALSE;
3680: PetscReal ptime;
3684: PetscCitationsRegister("@techreport{tspaper,\n"
3685: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3686: " author = {Shrirang Abhyankar and Jed Brown and Emil Constantinescu and Debojyoti Ghosh and Barry F. Smith},\n"
3687: " type = {Preprint},\n"
3688: " number = {ANL/MCS-P5061-0114},\n"
3689: " institution = {Argonne National Laboratory},\n"
3690: " year = {2014}\n}\n",&cite);
3692: TSSetUp(ts);
3693: TSTrajectorySetUp(ts->trajectory,ts);
3695: if (!ts->ops->step) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSStep not implemented for type '%s'",((PetscObject)ts)->type_name);
3696: 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>");
3697: 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()");
3698: 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");
3700: if (!ts->steps) ts->ptime_prev = ts->ptime;
3701: ptime = ts->ptime; ts->ptime_prev_rollback = ts->ptime_prev;
3702: ts->reason = TS_CONVERGED_ITERATING;
3704: PetscLogEventBegin(TS_Step,ts,0,0,0);
3705: (*ts->ops->step)(ts);
3706: PetscLogEventEnd(TS_Step,ts,0,0,0);
3708: if (ts->reason >= 0) {
3709: ts->ptime_prev = ptime;
3710: ts->steps++;
3711: ts->steprollback = PETSC_FALSE;
3712: ts->steprestart = PETSC_FALSE;
3713: }
3715: if (!ts->reason) {
3716: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3717: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3718: }
3720: if (ts->reason < 0 && ts->errorifstepfailed && 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]);
3721: if (ts->reason < 0 && ts->errorifstepfailed) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
3722: return(0);
3723: }
3725: /*@
3726: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3727: at the end of a time step with a given order of accuracy.
3729: Collective on TS
3731: Input Arguments:
3732: + ts - time stepping context
3733: . wnormtype - norm type, either NORM_2 or NORM_INFINITY
3734: - order - optional, desired order for the error evaluation or PETSC_DECIDE
3736: Output Arguments:
3737: + order - optional, the actual order of the error evaluation
3738: - wlte - the weighted local truncation error norm
3740: Level: advanced
3742: Notes:
3743: If the timestepper cannot evaluate the error in a particular step
3744: (eg. in the first step or restart steps after event handling),
3745: this routine returns wlte=-1.0 .
3747: .seealso: TSStep(), TSAdapt, TSErrorWeightedNorm()
3748: @*/
3749: PetscErrorCode TSEvaluateWLTE(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
3750: {
3760: if (wnormtype != NORM_2 && wnormtype != NORM_INFINITY) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
3761: if (!ts->ops->evaluatewlte) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateWLTE not implemented for type '%s'",((PetscObject)ts)->type_name);
3762: (*ts->ops->evaluatewlte)(ts,wnormtype,order,wlte);
3763: return(0);
3764: }
3766: /*@
3767: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3769: Collective on TS
3771: Input Arguments:
3772: + ts - time stepping context
3773: . order - desired order of accuracy
3774: - done - whether the step was evaluated at this order (pass NULL to generate an error if not available)
3776: Output Arguments:
3777: . U - state at the end of the current step
3779: Level: advanced
3781: Notes:
3782: This function cannot be called until all stages have been evaluated.
3783: 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.
3785: .seealso: TSStep(), TSAdapt
3786: @*/
3787: PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done)
3788: {
3795: if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name);
3796: (*ts->ops->evaluatestep)(ts,order,U,done);
3797: return(0);
3798: }
3800: /*@C
3801: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3803: Not collective
3805: Input Argument:
3806: . ts - time stepping context
3808: Output Argument:
3809: . initConditions - The function which computes an initial condition
3811: Level: advanced
3813: Notes:
3814: The calling sequence for the function is
3815: $ initCondition(TS ts, Vec u)
3816: $ ts - The timestepping context
3817: $ u - The input vector in which the initial condition is stored
3819: .seealso: TSSetComputeInitialCondition(), TSComputeInitialCondition()
3820: @*/
3821: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS, Vec))
3822: {
3826: *initCondition = ts->ops->initcondition;
3827: return(0);
3828: }
3830: /*@C
3831: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3833: Logically collective on ts
3835: Input Arguments:
3836: + ts - time stepping context
3837: - initCondition - The function which computes an initial condition
3839: Level: advanced
3841: Notes:
3842: The calling sequence for the function is
3843: $ initCondition(TS ts, Vec u)
3844: $ ts - The timestepping context
3845: $ u - The input vector in which the initial condition is stored
3847: .seealso: TSGetComputeInitialCondition(), TSComputeInitialCondition()
3848: @*/
3849: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS, Vec))
3850: {
3854: ts->ops->initcondition = initCondition;
3855: return(0);
3856: }
3858: /*@
3859: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set.
3861: Collective on ts
3863: Input Arguments:
3864: + ts - time stepping context
3865: - u - The Vec to store the condition in which will be used in TSSolve()
3867: Level: advanced
3869: Notes:
3870: The calling sequence for the function is
3871: $ initCondition(TS ts, Vec u)
3872: $ ts - The timestepping context
3873: $ u - The input vector in which the initial condition is stored
3875: .seealso: TSGetComputeInitialCondition(), TSSetComputeInitialCondition(), TSSolve()
3876: @*/
3877: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3878: {
3884: if (ts->ops->initcondition) {(*ts->ops->initcondition)(ts, u);}
3885: return(0);
3886: }
3888: /*@C
3889: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3891: Not collective
3893: Input Argument:
3894: . ts - time stepping context
3896: Output Argument:
3897: . exactError - The function which computes the solution error
3899: Level: advanced
3901: Notes:
3902: The calling sequence for the function is
3903: $ exactError(TS ts, Vec u)
3904: $ ts - The timestepping context
3905: $ u - The approximate solution vector
3906: $ e - The input vector in which the error is stored
3908: .seealso: TSGetComputeExactError(), TSComputeExactError()
3909: @*/
3910: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS, Vec, Vec))
3911: {
3915: *exactError = ts->ops->exacterror;
3916: return(0);
3917: }
3919: /*@C
3920: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3922: Logically collective on ts
3924: Input Arguments:
3925: + ts - time stepping context
3926: - exactError - The function which computes the solution error
3928: Level: advanced
3930: Notes:
3931: The calling sequence for the function is
3932: $ exactError(TS ts, Vec u)
3933: $ ts - The timestepping context
3934: $ u - The approximate solution vector
3935: $ e - The input vector in which the error is stored
3937: .seealso: TSGetComputeExactError(), TSComputeExactError()
3938: @*/
3939: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS, Vec, Vec))
3940: {
3944: ts->ops->exacterror = exactError;
3945: return(0);
3946: }
3948: /*@
3949: TSComputeExactError - Compute the solution error for the timestepping using the function previously set.
3951: Collective on ts
3953: Input Arguments:
3954: + ts - time stepping context
3955: . u - The approximate solution
3956: - e - The Vec used to store the error
3958: Level: advanced
3960: Notes:
3961: The calling sequence for the function is
3962: $ exactError(TS ts, Vec u)
3963: $ ts - The timestepping context
3964: $ u - The approximate solution vector
3965: $ e - The input vector in which the error is stored
3967: .seealso: TSGetComputeInitialCondition(), TSSetComputeInitialCondition(), TSSolve()
3968: @*/
3969: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3970: {
3977: if (ts->ops->exacterror) {(*ts->ops->exacterror)(ts, u, e);}
3978: return(0);
3979: }
3981: /*@
3982: TSSolve - Steps the requested number of timesteps.
3984: Collective on TS
3986: Input Parameter:
3987: + ts - the TS context obtained from TSCreate()
3988: - u - the solution vector (can be null if TSSetSolution() was used and TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP) was not used,
3989: otherwise must contain the initial conditions and will contain the solution at the final requested time
3991: Level: beginner
3993: Notes:
3994: The final time returned by this function may be different from the time of the internally
3995: held state accessible by TSGetSolution() and TSGetTime() because the method may have
3996: stepped over the final time.
3998: .seealso: TSCreate(), TSSetSolution(), TSStep(), TSGetTime(), TSGetSolveTime()
3999: @*/
4000: PetscErrorCode TSSolve(TS ts,Vec u)
4001: {
4002: Vec solution;
4003: PetscErrorCode ierr;
4008: 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 */
4009: if (!ts->vec_sol || u == ts->vec_sol) {
4010: VecDuplicate(u,&solution);
4011: TSSetSolution(ts,solution);
4012: VecDestroy(&solution); /* grant ownership */
4013: }
4014: VecCopy(u,ts->vec_sol);
4015: if (ts->forward_solve) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
4016: } else if (u) {
4017: TSSetSolution(ts,u);
4018: }
4019: TSSetUp(ts);
4020: TSTrajectorySetUp(ts->trajectory,ts);
4022: 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>");
4023: 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()");
4024: 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");
4026: if (ts->forward_solve) {
4027: TSForwardSetUp(ts);
4028: }
4030: /* reset number of steps only when the step is not restarted. ARKIMEX
4031: restarts the step after an event. Resetting these counters in such case causes
4032: TSTrajectory to incorrectly save the output files
4033: */
4034: /* reset time step and iteration counters */
4035: if (!ts->steps) {
4036: ts->ksp_its = 0;
4037: ts->snes_its = 0;
4038: ts->num_snes_failures = 0;
4039: ts->reject = 0;
4040: ts->steprestart = PETSC_TRUE;
4041: ts->steprollback = PETSC_FALSE;
4042: ts->rhsjacobian.time = PETSC_MIN_REAL;
4043: }
4045: /* make sure initial time step does not overshoot final time */
4046: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
4047: PetscReal maxdt = ts->max_time-ts->ptime;
4048: PetscReal dt = ts->time_step;
4050: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt,maxdt,10*PETSC_MACHINE_EPSILON,0) ? maxdt : dt);
4051: }
4052: ts->reason = TS_CONVERGED_ITERATING;
4054: {
4055: PetscViewer viewer;
4056: PetscViewerFormat format;
4057: PetscBool flg;
4058: static PetscBool incall = PETSC_FALSE;
4060: if (!incall) {
4061: /* Estimate the convergence rate of the time discretization */
4062: PetscOptionsGetViewer(PetscObjectComm((PetscObject) ts),((PetscObject)ts)->options, ((PetscObject) ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg);
4063: if (flg) {
4064: PetscConvEst conv;
4065: DM dm;
4066: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4067: PetscInt Nf;
4069: incall = PETSC_TRUE;
4070: TSGetDM(ts, &dm);
4071: DMGetNumFields(dm, &Nf);
4072: PetscCalloc1(PetscMax(Nf, 1), &alpha);
4073: PetscConvEstCreate(PetscObjectComm((PetscObject) ts), &conv);
4074: PetscConvEstUseTS(conv);
4075: PetscConvEstSetSolver(conv, (PetscObject) ts);
4076: PetscConvEstSetFromOptions(conv);
4077: PetscConvEstSetUp(conv);
4078: PetscConvEstGetConvRate(conv, alpha);
4079: PetscViewerPushFormat(viewer, format);
4080: PetscConvEstRateView(conv, alpha, viewer);
4081: PetscViewerPopFormat(viewer);
4082: PetscViewerDestroy(&viewer);
4083: PetscConvEstDestroy(&conv);
4084: PetscFree(alpha);
4085: incall = PETSC_FALSE;
4086: }
4087: }
4088: }
4090: TSViewFromOptions(ts,NULL,"-ts_view_pre");
4092: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4093: (*ts->ops->solve)(ts);
4094: if (u) {VecCopy(ts->vec_sol,u);}
4095: ts->solvetime = ts->ptime;
4096: solution = ts->vec_sol;
4097: } else { /* Step the requested number of timesteps. */
4098: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4099: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4101: if (!ts->steps) {
4102: TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);
4103: TSEventInitialize(ts->event,ts,ts->ptime,ts->vec_sol);
4104: }
4106: while (!ts->reason) {
4107: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
4108: if (!ts->steprollback) {
4109: TSPreStep(ts);
4110: }
4111: TSStep(ts);
4112: if (ts->testjacobian) {
4113: TSRHSJacobianTest(ts,NULL);
4114: }
4115: if (ts->testjacobiantranspose) {
4116: TSRHSJacobianTestTranspose(ts,NULL);
4117: }
4118: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4119: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4120: TSForwardCostIntegral(ts);
4121: if (ts->reason >= 0) ts->steps++;
4122: }
4123: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4124: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4125: TSForwardStep(ts);
4126: if (ts->reason >= 0) ts->steps++;
4127: }
4128: TSPostEvaluate(ts);
4129: TSEventHandler(ts); /* The right-hand side may be changed due to event. Be careful with Any computation using the RHS information after this point. */
4130: if (ts->steprollback) {
4131: TSPostEvaluate(ts);
4132: }
4133: if (!ts->steprollback) {
4134: TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);
4135: TSPostStep(ts);
4136: }
4137: }
4138: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
4140: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4141: TSInterpolate(ts,ts->max_time,u);
4142: ts->solvetime = ts->max_time;
4143: solution = u;
4144: TSMonitor(ts,-1,ts->solvetime,solution);
4145: } else {
4146: if (u) {VecCopy(ts->vec_sol,u);}
4147: ts->solvetime = ts->ptime;
4148: solution = ts->vec_sol;
4149: }
4150: }
4152: TSViewFromOptions(ts,NULL,"-ts_view");
4153: VecViewFromOptions(solution,NULL,"-ts_view_solution");
4154: PetscObjectSAWsBlock((PetscObject)ts);
4155: if (ts->adjoint_solve) {
4156: TSAdjointSolve(ts);
4157: }
4158: return(0);
4159: }
4161: /*@C
4162: TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet()
4164: Collective on TS
4166: Input Parameters:
4167: + ts - time stepping context obtained from TSCreate()
4168: . step - step number that has just completed
4169: . ptime - model time of the state
4170: - u - state at the current model time
4172: Notes:
4173: TSMonitor() is typically used automatically within the time stepping implementations.
4174: Users would almost never call this routine directly.
4176: A step of -1 indicates that the monitor is being called on a solution obtained by interpolating from computed solutions
4178: Level: developer
4180: @*/
4181: PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u)
4182: {
4183: DM dm;
4184: PetscInt i,n = ts->numbermonitors;
4191: TSGetDM(ts,&dm);
4192: DMSetOutputSequenceNumber(dm,step,ptime);
4194: VecLockReadPush(u);
4195: for (i=0; i<n; i++) {
4196: (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);
4197: }
4198: VecLockReadPop(u);
4199: return(0);
4200: }
4202: /* ------------------------------------------------------------------------*/
4203: /*@C
4204: TSMonitorLGCtxCreate - Creates a TSMonitorLGCtx context for use with
4205: TS to monitor the solution process graphically in various ways
4207: Collective on TS
4209: Input Parameters:
4210: + host - the X display to open, or null for the local machine
4211: . label - the title to put in the title bar
4212: . x, y - the screen coordinates of the upper left coordinate of the window
4213: . m, n - the screen width and height in pixels
4214: - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time
4216: Output Parameter:
4217: . ctx - the context
4219: Options Database Key:
4220: + -ts_monitor_lg_timestep - automatically sets line graph monitor
4221: + -ts_monitor_lg_timestep_log - automatically sets line graph monitor
4222: . -ts_monitor_lg_solution - monitor the solution (or certain values of the solution by calling TSMonitorLGSetDisplayVariables() or TSMonitorLGCtxSetDisplayVariables())
4223: . -ts_monitor_lg_error - monitor the error
4224: . -ts_monitor_lg_ksp_iterations - monitor the number of KSP iterations needed for each timestep
4225: . -ts_monitor_lg_snes_iterations - monitor the number of SNES iterations needed for each timestep
4226: - -lg_use_markers <true,false> - mark the data points (at each time step) on the plot; default is true
4228: Notes:
4229: Use TSMonitorLGCtxDestroy() to destroy.
4231: One can provide a function that transforms the solution before plotting it with TSMonitorLGCtxSetTransform() or TSMonitorLGSetTransform()
4233: Many of the functions that control the monitoring have two forms: TSMonitorLGSet/GetXXXX() and TSMonitorLGCtxSet/GetXXXX() the first take a TS object as the
4234: 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
4235: as the first argument.
4237: One can control the names displayed for each solution or error variable with TSMonitorLGCtxSetVariableNames() or TSMonitorLGSetVariableNames()
4239: Level: intermediate
4241: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError(), TSMonitorDefault(), VecView(),
4242: TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
4243: TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
4244: TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
4245: TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()
4247: @*/
4248: PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx)
4249: {
4250: PetscDraw draw;
4254: PetscNew(ctx);
4255: PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
4256: PetscDrawSetFromOptions(draw);
4257: PetscDrawLGCreate(draw,1,&(*ctx)->lg);
4258: PetscDrawLGSetFromOptions((*ctx)->lg);
4259: PetscDrawDestroy(&draw);
4260: (*ctx)->howoften = howoften;
4261: return(0);
4262: }
4264: PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
4265: {
4266: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
4267: PetscReal x = ptime,y;
4271: if (step < 0) return(0); /* -1 indicates an interpolated solution */
4272: if (!step) {
4273: PetscDrawAxis axis;
4274: const char *ylabel = ctx->semilogy ? "Log Time Step" : "Time Step";
4275: PetscDrawLGGetAxis(ctx->lg,&axis);
4276: PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time",ylabel);
4277: PetscDrawLGReset(ctx->lg);
4278: }
4279: TSGetTimeStep(ts,&y);
4280: if (ctx->semilogy) y = PetscLog10Real(y);
4281: PetscDrawLGAddPoint(ctx->lg,&x,&y);
4282: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
4283: PetscDrawLGDraw(ctx->lg);
4284: PetscDrawLGSave(ctx->lg);
4285: }
4286: return(0);
4287: }
4289: /*@C
4290: TSMonitorLGCtxDestroy - Destroys a line graph context that was created
4291: with TSMonitorLGCtxCreate().
4293: Collective on TSMonitorLGCtx
4295: Input Parameter:
4296: . ctx - the monitor context
4298: Level: intermediate
4300: .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep();
4301: @*/
4302: PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx)
4303: {
4307: if ((*ctx)->transformdestroy) {
4308: ((*ctx)->transformdestroy)((*ctx)->transformctx);
4309: }
4310: PetscDrawLGDestroy(&(*ctx)->lg);
4311: PetscStrArrayDestroy(&(*ctx)->names);
4312: PetscStrArrayDestroy(&(*ctx)->displaynames);
4313: PetscFree((*ctx)->displayvariables);
4314: PetscFree((*ctx)->displayvalues);
4315: PetscFree(*ctx);
4316: return(0);
4317: }
4319: /*
4321: Creates a TS Monitor SPCtx for use with DM Swarm particle visualizations
4323: */
4324: PetscErrorCode TSMonitorSPCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorSPCtx *ctx)
4325: {
4326: PetscDraw draw;
4330: PetscNew(ctx);
4331: PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
4332: PetscDrawSetFromOptions(draw);
4333: PetscDrawSPCreate(draw,1,&(*ctx)->sp);
4334: PetscDrawDestroy(&draw);
4335: (*ctx)->howoften = howoften;
4336: return(0);
4338: }
4340: /*
4341: Destroys a TSMonitorSPCtx that was created with TSMonitorSPCtxCreate
4342: */
4343: PetscErrorCode TSMonitorSPCtxDestroy(TSMonitorSPCtx *ctx)
4344: {
4349: PetscDrawSPDestroy(&(*ctx)->sp);
4350: PetscFree(*ctx);
4352: return(0);
4354: }
4356: /*@
4357: TSGetTime - Gets the time of the most recently completed step.
4359: Not Collective
4361: Input Parameter:
4362: . ts - the TS context obtained from TSCreate()
4364: Output Parameter:
4365: . t - the current time. This time may not corresponds to the final time set with TSSetMaxTime(), use TSGetSolveTime().
4367: Level: beginner
4369: Note:
4370: When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(),
4371: TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated.
4373: .seealso: TSGetSolveTime(), TSSetTime(), TSGetTimeStep(), TSGetStepNumber()
4375: @*/
4376: PetscErrorCode TSGetTime(TS ts,PetscReal *t)
4377: {
4381: *t = ts->ptime;
4382: return(0);
4383: }
4385: /*@
4386: TSGetPrevTime - Gets the starting time of the previously completed step.
4388: Not Collective
4390: Input Parameter:
4391: . ts - the TS context obtained from TSCreate()
4393: Output Parameter:
4394: . t - the previous time
4396: Level: beginner
4398: .seealso: TSGetTime(), TSGetSolveTime(), TSGetTimeStep()
4400: @*/
4401: PetscErrorCode TSGetPrevTime(TS ts,PetscReal *t)
4402: {
4406: *t = ts->ptime_prev;
4407: return(0);
4408: }
4410: /*@
4411: TSSetTime - Allows one to reset the time.
4413: Logically Collective on TS
4415: Input Parameters:
4416: + ts - the TS context obtained from TSCreate()
4417: - time - the time
4419: Level: intermediate
4421: .seealso: TSGetTime(), TSSetMaxSteps()
4423: @*/
4424: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4425: {
4429: ts->ptime = t;
4430: return(0);
4431: }
4433: /*@C
4434: TSSetOptionsPrefix - Sets the prefix used for searching for all
4435: TS options in the database.
4437: Logically Collective on TS
4439: Input Parameter:
4440: + ts - The TS context
4441: - prefix - The prefix to prepend to all option names
4443: Notes:
4444: A hyphen (-) must NOT be given at the beginning of the prefix name.
4445: The first character of all runtime options is AUTOMATICALLY the
4446: hyphen.
4448: Level: advanced
4450: .seealso: TSSetFromOptions()
4452: @*/
4453: PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[])
4454: {
4456: SNES snes;
4460: PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);
4461: TSGetSNES(ts,&snes);
4462: SNESSetOptionsPrefix(snes,prefix);
4463: return(0);
4464: }
4466: /*@C
4467: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4468: TS options in the database.
4470: Logically Collective on TS
4472: Input Parameter:
4473: + ts - The TS context
4474: - prefix - The prefix to prepend to all option names
4476: Notes:
4477: A hyphen (-) must NOT be given at the beginning of the prefix name.
4478: The first character of all runtime options is AUTOMATICALLY the
4479: hyphen.
4481: Level: advanced
4483: .seealso: TSGetOptionsPrefix()
4485: @*/
4486: PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[])
4487: {
4489: SNES snes;
4493: PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);
4494: TSGetSNES(ts,&snes);
4495: SNESAppendOptionsPrefix(snes,prefix);
4496: return(0);
4497: }
4499: /*@C
4500: TSGetOptionsPrefix - Sets the prefix used for searching for all
4501: TS options in the database.
4503: Not Collective
4505: Input Parameter:
4506: . ts - The TS context
4508: Output Parameter:
4509: . prefix - A pointer to the prefix string used
4511: Notes:
4512: On the fortran side, the user should pass in a string 'prifix' of
4513: sufficient length to hold the prefix.
4515: Level: intermediate
4517: .seealso: TSAppendOptionsPrefix()
4518: @*/
4519: PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[])
4520: {
4526: PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);
4527: return(0);
4528: }
4530: /*@C
4531: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4533: Not Collective, but parallel objects are returned if TS is parallel
4535: Input Parameter:
4536: . ts - The TS context obtained from TSCreate()
4538: Output Parameters:
4539: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL)
4540: . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL)
4541: . func - Function to compute the Jacobian of the RHS (or NULL)
4542: - ctx - User-defined context for Jacobian evaluation routine (or NULL)
4544: Notes:
4545: You can pass in NULL for any return argument you do not need.
4547: Level: intermediate
4549: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()
4551: @*/
4552: PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx)
4553: {
4555: DM dm;
4558: if (Amat || Pmat) {
4559: SNES snes;
4560: TSGetSNES(ts,&snes);
4561: SNESSetUpMatrices(snes);
4562: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
4563: }
4564: TSGetDM(ts,&dm);
4565: DMTSGetRHSJacobian(dm,func,ctx);
4566: return(0);
4567: }
4569: /*@C
4570: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4572: Not Collective, but parallel objects are returned if TS is parallel
4574: Input Parameter:
4575: . ts - The TS context obtained from TSCreate()
4577: Output Parameters:
4578: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4579: . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat
4580: . f - The function to compute the matrices
4581: - ctx - User-defined context for Jacobian evaluation routine
4583: Notes:
4584: You can pass in NULL for any return argument you do not need.
4586: Level: advanced
4588: .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()
4590: @*/
4591: PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx)
4592: {
4594: DM dm;
4597: if (Amat || Pmat) {
4598: SNES snes;
4599: TSGetSNES(ts,&snes);
4600: SNESSetUpMatrices(snes);
4601: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
4602: }
4603: TSGetDM(ts,&dm);
4604: DMTSGetIJacobian(dm,f,ctx);
4605: return(0);
4606: }
4608: /*@C
4609: TSMonitorDrawSolution - Monitors progress of the TS solvers by calling
4610: VecView() for the solution at each timestep
4612: Collective on TS
4614: Input Parameters:
4615: + ts - the TS context
4616: . step - current time-step
4617: . ptime - current time
4618: - dummy - either a viewer or NULL
4620: Options Database:
4621: . -ts_monitor_draw_solution_initial - show initial solution as well as current solution
4623: Notes:
4624: the initial solution and current solution are not display with a common axis scaling so generally the option -ts_monitor_draw_solution_initial
4625: will look bad
4627: Level: intermediate
4629: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4630: @*/
4631: PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4632: {
4633: PetscErrorCode ierr;
4634: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
4635: PetscDraw draw;
4638: if (!step && ictx->showinitial) {
4639: if (!ictx->initialsolution) {
4640: VecDuplicate(u,&ictx->initialsolution);
4641: }
4642: VecCopy(u,ictx->initialsolution);
4643: }
4644: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);
4646: if (ictx->showinitial) {
4647: PetscReal pause;
4648: PetscViewerDrawGetPause(ictx->viewer,&pause);
4649: PetscViewerDrawSetPause(ictx->viewer,0.0);
4650: VecView(ictx->initialsolution,ictx->viewer);
4651: PetscViewerDrawSetPause(ictx->viewer,pause);
4652: PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);
4653: }
4654: VecView(u,ictx->viewer);
4655: if (ictx->showtimestepandtime) {
4656: PetscReal xl,yl,xr,yr,h;
4657: char time[32];
4659: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
4660: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
4661: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
4662: h = yl + .95*(yr - yl);
4663: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
4664: PetscDrawFlush(draw);
4665: }
4667: if (ictx->showinitial) {
4668: PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);
4669: }
4670: return(0);
4671: }
4673: /*@C
4674: TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram
4676: Collective on TS
4678: Input Parameters:
4679: + ts - the TS context
4680: . step - current time-step
4681: . ptime - current time
4682: - dummy - either a viewer or NULL
4684: Level: intermediate
4686: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4687: @*/
4688: PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4689: {
4690: PetscErrorCode ierr;
4691: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
4692: PetscDraw draw;
4693: PetscDrawAxis axis;
4694: PetscInt n;
4695: PetscMPIInt size;
4696: PetscReal U0,U1,xl,yl,xr,yr,h;
4697: char time[32];
4698: const PetscScalar *U;
4701: MPI_Comm_size(PetscObjectComm((PetscObject)ts),&size);
4702: if (size != 1) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Only allowed for sequential runs");
4703: VecGetSize(u,&n);
4704: if (n != 2) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Only for ODEs with two unknowns");
4706: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
4707: PetscViewerDrawGetDrawAxis(ictx->viewer,0,&axis);
4708: PetscDrawAxisGetLimits(axis,&xl,&xr,&yl,&yr);
4709: if (!step) {
4710: PetscDrawClear(draw);
4711: PetscDrawAxisDraw(axis);
4712: }
4714: VecGetArrayRead(u,&U);
4715: U0 = PetscRealPart(U[0]);
4716: U1 = PetscRealPart(U[1]);
4717: VecRestoreArrayRead(u,&U);
4718: if ((U0 < xl) || (U1 < yl) || (U0 > xr) || (U1 > yr)) return(0);
4720: PetscDrawCollectiveBegin(draw);
4721: PetscDrawPoint(draw,U0,U1,PETSC_DRAW_BLACK);
4722: if (ictx->showtimestepandtime) {
4723: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
4724: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
4725: h = yl + .95*(yr - yl);
4726: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
4727: }
4728: PetscDrawCollectiveEnd(draw);
4729: PetscDrawFlush(draw);
4730: PetscDrawPause(draw);
4731: PetscDrawSave(draw);
4732: return(0);
4733: }
4735: /*@C
4736: TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution()
4738: Collective on TS
4740: Input Parameters:
4741: . ctx - the monitor context
4743: Level: intermediate
4745: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError()
4746: @*/
4747: PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx)
4748: {
4752: PetscViewerDestroy(&(*ictx)->viewer);
4753: VecDestroy(&(*ictx)->initialsolution);
4754: PetscFree(*ictx);
4755: return(0);
4756: }
4758: /*@C
4759: TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx
4761: Collective on TS
4763: Input Parameter:
4764: . ts - time-step context
4766: Output Patameter:
4767: . ctx - the monitor context
4769: Options Database:
4770: . -ts_monitor_draw_solution_initial - show initial solution as well as current solution
4772: Level: intermediate
4774: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx()
4775: @*/
4776: PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx)
4777: {
4778: PetscErrorCode ierr;
4781: PetscNew(ctx);
4782: PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);
4783: PetscViewerSetFromOptions((*ctx)->viewer);
4785: (*ctx)->howoften = howoften;
4786: (*ctx)->showinitial = PETSC_FALSE;
4787: PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);
4789: (*ctx)->showtimestepandtime = PETSC_FALSE;
4790: PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);
4791: return(0);
4792: }
4794: /*@C
4795: TSMonitorDrawSolutionFunction - Monitors progress of the TS solvers by calling
4796: VecView() for the solution provided by TSSetSolutionFunction() at each timestep
4798: Collective on TS
4800: Input Parameters:
4801: + ts - the TS context
4802: . step - current time-step
4803: . ptime - current time
4804: - dummy - either a viewer or NULL
4806: Options Database:
4807: . -ts_monitor_draw_solution_function - Monitor error graphically, requires user to have provided TSSetSolutionFunction()
4809: Level: intermediate
4811: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
4812: @*/
4813: PetscErrorCode TSMonitorDrawSolutionFunction(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4814: {
4815: PetscErrorCode ierr;
4816: TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy;
4817: PetscViewer viewer = ctx->viewer;
4818: Vec work;
4821: if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return(0);
4822: VecDuplicate(u,&work);
4823: TSComputeSolutionFunction(ts,ptime,work);
4824: VecView(work,viewer);
4825: VecDestroy(&work);
4826: return(0);
4827: }
4829: /*@C
4830: TSMonitorDrawError - Monitors progress of the TS solvers by calling
4831: VecView() for the error at each timestep
4833: Collective on TS
4835: Input Parameters:
4836: + ts - the TS context
4837: . step - current time-step
4838: . ptime - current time
4839: - dummy - either a viewer or NULL
4841: Options Database:
4842: . -ts_monitor_draw_error - Monitor error graphically, requires user to have provided TSSetSolutionFunction()
4844: Level: intermediate
4846: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
4847: @*/
4848: PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4849: {
4850: PetscErrorCode ierr;
4851: TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy;
4852: PetscViewer viewer = ctx->viewer;
4853: Vec work;
4856: if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return(0);
4857: VecDuplicate(u,&work);
4858: TSComputeSolutionFunction(ts,ptime,work);
4859: VecAXPY(work,-1.0,u);
4860: VecView(work,viewer);
4861: VecDestroy(&work);
4862: return(0);
4863: }
4865: #include <petsc/private/dmimpl.h>
4866: /*@
4867: TSSetDM - Sets the DM that may be used by some nonlinear solvers or preconditioners under the TS
4869: Logically Collective on ts
4871: Input Parameters:
4872: + ts - the ODE integrator object
4873: - dm - the dm, cannot be NULL
4875: Notes:
4876: A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
4877: even when not using interfaces like DMTSSetIFunction(). Use DMClone() to get a distinct DM when solving
4878: different problems using the same function space.
4880: Level: intermediate
4882: .seealso: TSGetDM(), SNESSetDM(), SNESGetDM()
4883: @*/
4884: PetscErrorCode TSSetDM(TS ts,DM dm)
4885: {
4887: SNES snes;
4888: DMTS tsdm;
4893: PetscObjectReference((PetscObject)dm);
4894: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4895: if (ts->dm->dmts && !dm->dmts) {
4896: DMCopyDMTS(ts->dm,dm);
4897: DMGetDMTS(ts->dm,&tsdm);
4898: if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */
4899: tsdm->originaldm = dm;
4900: }
4901: }
4902: DMDestroy(&ts->dm);
4903: }
4904: ts->dm = dm;
4906: TSGetSNES(ts,&snes);
4907: SNESSetDM(snes,dm);
4908: return(0);
4909: }
4911: /*@
4912: TSGetDM - Gets the DM that may be used by some preconditioners
4914: Not Collective
4916: Input Parameter:
4917: . ts - the preconditioner context
4919: Output Parameter:
4920: . dm - the dm
4922: Level: intermediate
4924: .seealso: TSSetDM(), SNESSetDM(), SNESGetDM()
4925: @*/
4926: PetscErrorCode TSGetDM(TS ts,DM *dm)
4927: {
4932: if (!ts->dm) {
4933: DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);
4934: if (ts->snes) {SNESSetDM(ts->snes,ts->dm);}
4935: }
4936: *dm = ts->dm;
4937: return(0);
4938: }
4940: /*@
4941: SNESTSFormFunction - Function to evaluate nonlinear residual
4943: Logically Collective on SNES
4945: Input Parameter:
4946: + snes - nonlinear solver
4947: . U - the current state at which to evaluate the residual
4948: - ctx - user context, must be a TS
4950: Output Parameter:
4951: . F - the nonlinear residual
4953: Notes:
4954: This function is not normally called by users and is automatically registered with the SNES used by TS.
4955: It is most frequently passed to MatFDColoringSetFunction().
4957: Level: advanced
4959: .seealso: SNESSetFunction(), MatFDColoringSetFunction()
4960: @*/
4961: PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx)
4962: {
4963: TS ts = (TS)ctx;
4971: (ts->ops->snesfunction)(snes,U,F,ts);
4972: return(0);
4973: }
4975: /*@
4976: SNESTSFormJacobian - Function to evaluate the Jacobian
4978: Collective on SNES
4980: Input Parameter:
4981: + snes - nonlinear solver
4982: . U - the current state at which to evaluate the residual
4983: - ctx - user context, must be a TS
4985: Output Parameter:
4986: + A - the Jacobian
4987: . B - the preconditioning matrix (may be the same as A)
4988: - flag - indicates any structure change in the matrix
4990: Notes:
4991: This function is not normally called by users and is automatically registered with the SNES used by TS.
4993: Level: developer
4995: .seealso: SNESSetJacobian()
4996: @*/
4997: PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx)
4998: {
4999: TS ts = (TS)ctx;
5010: (ts->ops->snesjacobian)(snes,U,A,B,ts);
5011: return(0);
5012: }
5014: /*@C
5015: TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems Udot = A U only
5017: Collective on TS
5019: Input Arguments:
5020: + ts - time stepping context
5021: . t - time at which to evaluate
5022: . U - state at which to evaluate
5023: - ctx - context
5025: Output Arguments:
5026: . F - right hand side
5028: Level: intermediate
5030: Notes:
5031: This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems.
5032: The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian().
5034: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
5035: @*/
5036: PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
5037: {
5039: Mat Arhs,Brhs;
5042: TSGetRHSMats_Private(ts,&Arhs,&Brhs);
5043: TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);
5044: MatMult(Arhs,U,F);
5045: return(0);
5046: }
5048: /*@C
5049: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
5051: Collective on TS
5053: Input Arguments:
5054: + ts - time stepping context
5055: . t - time at which to evaluate
5056: . U - state at which to evaluate
5057: - ctx - context
5059: Output Arguments:
5060: + A - pointer to operator
5061: . B - pointer to preconditioning matrix
5062: - flg - matrix structure flag
5064: Level: intermediate
5066: Notes:
5067: This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems.
5069: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear()
5070: @*/
5071: PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
5072: {
5074: return(0);
5075: }
5077: /*@C
5078: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
5080: Collective on TS
5082: Input Arguments:
5083: + ts - time stepping context
5084: . t - time at which to evaluate
5085: . U - state at which to evaluate
5086: . Udot - time derivative of state vector
5087: - ctx - context
5089: Output Arguments:
5090: . F - left hand side
5092: Level: intermediate
5094: Notes:
5095: 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
5096: user is required to write their own TSComputeIFunction.
5097: This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems.
5098: The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian().
5100: Note that using this function is NOT equivalent to using TSComputeRHSFunctionLinear() since that solves Udot = A U
5102: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant(), TSComputeRHSFunctionLinear()
5103: @*/
5104: PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
5105: {
5107: Mat A,B;
5110: TSGetIJacobian(ts,&A,&B,NULL,NULL);
5111: TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);
5112: MatMult(A,Udot,F);
5113: return(0);
5114: }
5116: /*@C
5117: TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE
5119: Collective on TS
5121: Input Arguments:
5122: + ts - time stepping context
5123: . t - time at which to evaluate
5124: . U - state at which to evaluate
5125: . Udot - time derivative of state vector
5126: . shift - shift to apply
5127: - ctx - context
5129: Output Arguments:
5130: + A - pointer to operator
5131: . B - pointer to preconditioning matrix
5132: - flg - matrix structure flag
5134: Level: advanced
5136: Notes:
5137: This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems.
5139: It is only appropriate for problems of the form
5141: $ M Udot = F(U,t)
5143: where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only
5144: works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing
5145: an implicit operator of the form
5147: $ shift*M + J
5149: 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
5150: a copy of M or reassemble it when requested.
5152: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear()
5153: @*/
5154: PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx)
5155: {
5159: MatScale(A, shift / ts->ijacobian.shift);
5160: ts->ijacobian.shift = shift;
5161: return(0);
5162: }
5164: /*@
5165: TSGetEquationType - Gets the type of the equation that TS is solving.
5167: Not Collective
5169: Input Parameter:
5170: . ts - the TS context
5172: Output Parameter:
5173: . equation_type - see TSEquationType
5175: Level: beginner
5177: .seealso: TSSetEquationType(), TSEquationType
5178: @*/
5179: PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type)
5180: {
5184: *equation_type = ts->equation_type;
5185: return(0);
5186: }
5188: /*@
5189: TSSetEquationType - Sets the type of the equation that TS is solving.
5191: Not Collective
5193: Input Parameter:
5194: + ts - the TS context
5195: - equation_type - see TSEquationType
5197: Level: advanced
5199: .seealso: TSGetEquationType(), TSEquationType
5200: @*/
5201: PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type)
5202: {
5205: ts->equation_type = equation_type;
5206: return(0);
5207: }
5209: /*@
5210: TSGetConvergedReason - Gets the reason the TS iteration was stopped.
5212: Not Collective
5214: Input Parameter:
5215: . ts - the TS context
5217: Output Parameter:
5218: . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
5219: manual pages for the individual convergence tests for complete lists
5221: Level: beginner
5223: Notes:
5224: Can only be called after the call to TSSolve() is complete.
5226: .seealso: TSSetConvergenceTest(), TSConvergedReason
5227: @*/
5228: PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason)
5229: {
5233: *reason = ts->reason;
5234: return(0);
5235: }
5237: /*@
5238: TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve.
5240: Logically Collective; reason must contain common value
5242: Input Parameters:
5243: + ts - the TS context
5244: - reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
5245: manual pages for the individual convergence tests for complete lists
5247: Level: advanced
5249: Notes:
5250: Can only be called while TSSolve() is active.
5252: .seealso: TSConvergedReason
5253: @*/
5254: PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason)
5255: {
5258: ts->reason = reason;
5259: return(0);
5260: }
5262: /*@
5263: TSGetSolveTime - Gets the time after a call to TSSolve()
5265: Not Collective
5267: Input Parameter:
5268: . ts - the TS context
5270: Output Parameter:
5271: . ftime - the final time. This time corresponds to the final time set with TSSetMaxTime()
5273: Level: beginner
5275: Notes:
5276: Can only be called after the call to TSSolve() is complete.
5278: .seealso: TSSetConvergenceTest(), TSConvergedReason
5279: @*/
5280: PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime)
5281: {
5285: *ftime = ts->solvetime;
5286: return(0);
5287: }
5289: /*@
5290: TSGetSNESIterations - Gets the total number of nonlinear iterations
5291: used by the time integrator.
5293: Not Collective
5295: Input Parameter:
5296: . ts - TS context
5298: Output Parameter:
5299: . nits - number of nonlinear iterations
5301: Notes:
5302: This counter is reset to zero for each successive call to TSSolve().
5304: Level: intermediate
5306: .seealso: TSGetKSPIterations()
5307: @*/
5308: PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits)
5309: {
5313: *nits = ts->snes_its;
5314: return(0);
5315: }
5317: /*@
5318: TSGetKSPIterations - Gets the total number of linear iterations
5319: used by the time integrator.
5321: Not Collective
5323: Input Parameter:
5324: . ts - TS context
5326: Output Parameter:
5327: . lits - number of linear iterations
5329: Notes:
5330: This counter is reset to zero for each successive call to TSSolve().
5332: Level: intermediate
5334: .seealso: TSGetSNESIterations(), SNESGetKSPIterations()
5335: @*/
5336: PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits)
5337: {
5341: *lits = ts->ksp_its;
5342: return(0);
5343: }
5345: /*@
5346: TSGetStepRejections - Gets the total number of rejected steps.
5348: Not Collective
5350: Input Parameter:
5351: . ts - TS context
5353: Output Parameter:
5354: . rejects - number of steps rejected
5356: Notes:
5357: This counter is reset to zero for each successive call to TSSolve().
5359: Level: intermediate
5361: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails()
5362: @*/
5363: PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects)
5364: {
5368: *rejects = ts->reject;
5369: return(0);
5370: }
5372: /*@
5373: TSGetSNESFailures - Gets the total number of failed SNES solves
5375: Not Collective
5377: Input Parameter:
5378: . ts - TS context
5380: Output Parameter:
5381: . fails - number of failed nonlinear solves
5383: Notes:
5384: This counter is reset to zero for each successive call to TSSolve().
5386: Level: intermediate
5388: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures()
5389: @*/
5390: PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails)
5391: {
5395: *fails = ts->num_snes_failures;
5396: return(0);
5397: }
5399: /*@
5400: TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails
5402: Not Collective
5404: Input Parameter:
5405: + ts - TS context
5406: - rejects - maximum number of rejected steps, pass -1 for unlimited
5408: Notes:
5409: The counter is reset to zero for each step
5411: Options Database Key:
5412: . -ts_max_reject - Maximum number of step rejections before a step fails
5414: Level: intermediate
5416: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
5417: @*/
5418: PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects)
5419: {
5422: ts->max_reject = rejects;
5423: return(0);
5424: }
5426: /*@
5427: TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves
5429: Not Collective
5431: Input Parameter:
5432: + ts - TS context
5433: - fails - maximum number of failed nonlinear solves, pass -1 for unlimited
5435: Notes:
5436: The counter is reset to zero for each successive call to TSSolve().
5438: Options Database Key:
5439: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
5441: Level: intermediate
5443: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason()
5444: @*/
5445: PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails)
5446: {
5449: ts->max_snes_failures = fails;
5450: return(0);
5451: }
5453: /*@
5454: TSSetErrorIfStepFails - Error if no step succeeds
5456: Not Collective
5458: Input Parameter:
5459: + ts - TS context
5460: - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure
5462: Options Database Key:
5463: . -ts_error_if_step_fails - Error if no step succeeds
5465: Level: intermediate
5467: .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
5468: @*/
5469: PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err)
5470: {
5473: ts->errorifstepfailed = err;
5474: return(0);
5475: }
5477: /*@C
5478: 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
5480: Collective on TS
5482: Input Parameters:
5483: + ts - the TS context
5484: . step - current time-step
5485: . ptime - current time
5486: . u - current state
5487: - vf - viewer and its format
5489: Level: intermediate
5491: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5492: @*/
5493: PetscErrorCode TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscViewerAndFormat *vf)
5494: {
5498: PetscViewerPushFormat(vf->viewer,vf->format);
5499: VecView(u,vf->viewer);
5500: PetscViewerPopFormat(vf->viewer);
5501: return(0);
5502: }
5504: /*@C
5505: TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep.
5507: Collective on TS
5509: Input Parameters:
5510: + ts - the TS context
5511: . step - current time-step
5512: . ptime - current time
5513: . u - current state
5514: - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
5516: Level: intermediate
5518: Notes:
5519: 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.
5520: These are named according to the file name template.
5522: This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy().
5524: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5525: @*/
5526: PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate)
5527: {
5529: char filename[PETSC_MAX_PATH_LEN];
5530: PetscViewer viewer;
5533: if (step < 0) return(0); /* -1 indicates interpolated solution */
5534: PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);
5535: PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);
5536: VecView(u,viewer);
5537: PetscViewerDestroy(&viewer);
5538: return(0);
5539: }
5541: /*@C
5542: TSMonitorSolutionVTKDestroy - Destroy context for monitoring
5544: Collective on TS
5546: Input Parameters:
5547: . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
5549: Level: intermediate
5551: Note:
5552: This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK().
5554: .seealso: TSMonitorSet(), TSMonitorSolutionVTK()
5555: @*/
5556: PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate)
5557: {
5561: PetscFree(*(char**)filenametemplate);
5562: return(0);
5563: }
5565: /*@
5566: TSGetAdapt - Get the adaptive controller context for the current method
5568: Collective on TS if controller has not been created yet
5570: Input Arguments:
5571: . ts - time stepping context
5573: Output Arguments:
5574: . adapt - adaptive controller
5576: Level: intermediate
5578: .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose()
5579: @*/
5580: PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt)
5581: {
5587: if (!ts->adapt) {
5588: TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);
5589: PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);
5590: PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);
5591: }
5592: *adapt = ts->adapt;
5593: return(0);
5594: }
5596: /*@
5597: TSSetTolerances - Set tolerances for local truncation error when using adaptive controller
5599: Logically Collective
5601: Input Arguments:
5602: + ts - time integration context
5603: . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value
5604: . vatol - vector of absolute tolerances or NULL, used in preference to atol if present
5605: . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value
5606: - vrtol - vector of relative tolerances or NULL, used in preference to atol if present
5608: Options Database keys:
5609: + -ts_rtol <rtol> - relative tolerance for local truncation error
5610: - -ts_atol <atol> Absolute tolerance for local truncation error
5612: Notes:
5613: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5614: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5615: computed only for the differential or the algebraic part then this can be done using the vector of
5616: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5617: differential part and infinity for the algebraic part, the LTE calculation will include only the
5618: differential variables.
5620: Level: beginner
5622: .seealso: TS, TSAdapt, TSErrorWeightedNorm(), TSGetTolerances()
5623: @*/
5624: PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol)
5625: {
5629: if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol;
5630: if (vatol) {
5631: PetscObjectReference((PetscObject)vatol);
5632: VecDestroy(&ts->vatol);
5633: ts->vatol = vatol;
5634: }
5635: if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol;
5636: if (vrtol) {
5637: PetscObjectReference((PetscObject)vrtol);
5638: VecDestroy(&ts->vrtol);
5639: ts->vrtol = vrtol;
5640: }
5641: return(0);
5642: }
5644: /*@
5645: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5647: Logically Collective
5649: Input Arguments:
5650: . ts - time integration context
5652: Output Arguments:
5653: + atol - scalar absolute tolerances, NULL to ignore
5654: . vatol - vector of absolute tolerances, NULL to ignore
5655: . rtol - scalar relative tolerances, NULL to ignore
5656: - vrtol - vector of relative tolerances, NULL to ignore
5658: Level: beginner
5660: .seealso: TS, TSAdapt, TSErrorWeightedNorm(), TSSetTolerances()
5661: @*/
5662: PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol)
5663: {
5665: if (atol) *atol = ts->atol;
5666: if (vatol) *vatol = ts->vatol;
5667: if (rtol) *rtol = ts->rtol;
5668: if (vrtol) *vrtol = ts->vrtol;
5669: return(0);
5670: }
5672: /*@
5673: TSErrorWeightedNorm2 - compute a weighted 2-norm of the difference between two state vectors
5675: Collective on TS
5677: Input Arguments:
5678: + ts - time stepping context
5679: . U - state vector, usually ts->vec_sol
5680: - Y - state vector to be compared to U
5682: Output Arguments:
5683: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5684: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5685: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
5687: Level: developer
5689: .seealso: TSErrorWeightedNorm(), TSErrorWeightedNormInfinity()
5690: @*/
5691: PetscErrorCode TSErrorWeightedNorm2(TS ts,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
5692: {
5693: PetscErrorCode ierr;
5694: PetscInt i,n,N,rstart;
5695: PetscInt n_loc,na_loc,nr_loc;
5696: PetscReal n_glb,na_glb,nr_glb;
5697: const PetscScalar *u,*y;
5698: PetscReal sum,suma,sumr,gsum,gsuma,gsumr,diff;
5699: PetscReal tol,tola,tolr;
5700: PetscReal err_loc[6],err_glb[6];
5712: if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector");
5714: VecGetSize(U,&N);
5715: VecGetLocalSize(U,&n);
5716: VecGetOwnershipRange(U,&rstart,NULL);
5717: VecGetArrayRead(U,&u);
5718: VecGetArrayRead(Y,&y);
5719: sum = 0.; n_loc = 0;
5720: suma = 0.; na_loc = 0;
5721: sumr = 0.; nr_loc = 0;
5722: if (ts->vatol && ts->vrtol) {
5723: const PetscScalar *atol,*rtol;
5724: VecGetArrayRead(ts->vatol,&atol);
5725: VecGetArrayRead(ts->vrtol,&rtol);
5726: for (i=0; i<n; i++) {
5727: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
5728: diff = PetscAbsScalar(y[i] - u[i]);
5729: tola = PetscRealPart(atol[i]);
5730: if(tola>0.){
5731: suma += PetscSqr(diff/tola);
5732: na_loc++;
5733: }
5734: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5735: if(tolr>0.){
5736: sumr += PetscSqr(diff/tolr);
5737: nr_loc++;
5738: }
5739: tol=tola+tolr;
5740: if(tol>0.){
5741: sum += PetscSqr(diff/tol);
5742: n_loc++;
5743: }
5744: }
5745: VecRestoreArrayRead(ts->vatol,&atol);
5746: VecRestoreArrayRead(ts->vrtol,&rtol);
5747: } else if (ts->vatol) { /* vector atol, scalar rtol */
5748: const PetscScalar *atol;
5749: VecGetArrayRead(ts->vatol,&atol);
5750: for (i=0; i<n; i++) {
5751: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
5752: diff = PetscAbsScalar(y[i] - u[i]);
5753: tola = PetscRealPart(atol[i]);
5754: if(tola>0.){
5755: suma += PetscSqr(diff/tola);
5756: na_loc++;
5757: }
5758: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5759: if(tolr>0.){
5760: sumr += PetscSqr(diff/tolr);
5761: nr_loc++;
5762: }
5763: tol=tola+tolr;
5764: if(tol>0.){
5765: sum += PetscSqr(diff/tol);
5766: n_loc++;
5767: }
5768: }
5769: VecRestoreArrayRead(ts->vatol,&atol);
5770: } else if (ts->vrtol) { /* scalar atol, vector rtol */
5771: const PetscScalar *rtol;
5772: VecGetArrayRead(ts->vrtol,&rtol);
5773: for (i=0; i<n; i++) {
5774: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
5775: diff = PetscAbsScalar(y[i] - u[i]);
5776: tola = ts->atol;
5777: if(tola>0.){
5778: suma += PetscSqr(diff/tola);
5779: na_loc++;
5780: }
5781: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5782: if(tolr>0.){
5783: sumr += PetscSqr(diff/tolr);
5784: nr_loc++;
5785: }
5786: tol=tola+tolr;
5787: if(tol>0.){
5788: sum += PetscSqr(diff/tol);
5789: n_loc++;
5790: }
5791: }
5792: VecRestoreArrayRead(ts->vrtol,&rtol);
5793: } else { /* scalar atol, scalar rtol */
5794: for (i=0; i<n; i++) {
5795: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
5796: diff = PetscAbsScalar(y[i] - u[i]);
5797: tola = ts->atol;
5798: if(tola>0.){
5799: suma += PetscSqr(diff/tola);
5800: na_loc++;
5801: }
5802: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5803: if(tolr>0.){
5804: sumr += PetscSqr(diff/tolr);
5805: nr_loc++;
5806: }
5807: tol=tola+tolr;
5808: if(tol>0.){
5809: sum += PetscSqr(diff/tol);
5810: n_loc++;
5811: }
5812: }
5813: }
5814: VecRestoreArrayRead(U,&u);
5815: VecRestoreArrayRead(Y,&y);
5817: err_loc[0] = sum;
5818: err_loc[1] = suma;
5819: err_loc[2] = sumr;
5820: err_loc[3] = (PetscReal)n_loc;
5821: err_loc[4] = (PetscReal)na_loc;
5822: err_loc[5] = (PetscReal)nr_loc;
5824: MPIU_Allreduce(err_loc,err_glb,6,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));
5826: gsum = err_glb[0];
5827: gsuma = err_glb[1];
5828: gsumr = err_glb[2];
5829: n_glb = err_glb[3];
5830: na_glb = err_glb[4];
5831: nr_glb = err_glb[5];
5833: *norm = 0.;
5834: if(n_glb>0. ){*norm = PetscSqrtReal(gsum / n_glb );}
5835: *norma = 0.;
5836: if(na_glb>0.){*norma = PetscSqrtReal(gsuma / na_glb);}
5837: *normr = 0.;
5838: if(nr_glb>0.){*normr = PetscSqrtReal(gsumr / nr_glb);}
5840: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
5841: if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
5842: if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
5843: return(0);
5844: }
5846: /*@
5847: TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between two state vectors
5849: Collective on TS
5851: Input Arguments:
5852: + ts - time stepping context
5853: . U - state vector, usually ts->vec_sol
5854: - Y - state vector to be compared to U
5856: Output Arguments:
5857: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5858: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5859: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
5861: Level: developer
5863: .seealso: TSErrorWeightedNorm(), TSErrorWeightedNorm2()
5864: @*/
5865: PetscErrorCode TSErrorWeightedNormInfinity(TS ts,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
5866: {
5867: PetscErrorCode ierr;
5868: PetscInt i,n,N,rstart;
5869: const PetscScalar *u,*y;
5870: PetscReal max,gmax,maxa,gmaxa,maxr,gmaxr;
5871: PetscReal tol,tola,tolr,diff;
5872: PetscReal err_loc[3],err_glb[3];
5884: if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector");
5886: VecGetSize(U,&N);
5887: VecGetLocalSize(U,&n);
5888: VecGetOwnershipRange(U,&rstart,NULL);
5889: VecGetArrayRead(U,&u);
5890: VecGetArrayRead(Y,&y);
5892: max=0.;
5893: maxa=0.;
5894: maxr=0.;
5896: if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
5897: const PetscScalar *atol,*rtol;
5898: VecGetArrayRead(ts->vatol,&atol);
5899: VecGetArrayRead(ts->vrtol,&rtol);
5901: for (i=0; i<n; i++) {
5902: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
5903: diff = PetscAbsScalar(y[i] - u[i]);
5904: tola = PetscRealPart(atol[i]);
5905: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5906: tol = tola+tolr;
5907: if(tola>0.){
5908: maxa = PetscMax(maxa,diff / tola);
5909: }
5910: if(tolr>0.){
5911: maxr = PetscMax(maxr,diff / tolr);
5912: }
5913: if(tol>0.){
5914: max = PetscMax(max,diff / tol);
5915: }
5916: }
5917: VecRestoreArrayRead(ts->vatol,&atol);
5918: VecRestoreArrayRead(ts->vrtol,&rtol);
5919: } else if (ts->vatol) { /* vector atol, scalar rtol */
5920: const PetscScalar *atol;
5921: VecGetArrayRead(ts->vatol,&atol);
5922: for (i=0; i<n; i++) {
5923: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
5924: diff = PetscAbsScalar(y[i] - u[i]);
5925: tola = PetscRealPart(atol[i]);
5926: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5927: tol = tola+tolr;
5928: if(tola>0.){
5929: maxa = PetscMax(maxa,diff / tola);
5930: }
5931: if(tolr>0.){
5932: maxr = PetscMax(maxr,diff / tolr);
5933: }
5934: if(tol>0.){
5935: max = PetscMax(max,diff / tol);
5936: }
5937: }
5938: VecRestoreArrayRead(ts->vatol,&atol);
5939: } else if (ts->vrtol) { /* scalar atol, vector rtol */
5940: const PetscScalar *rtol;
5941: VecGetArrayRead(ts->vrtol,&rtol);
5943: for (i=0; i<n; i++) {
5944: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
5945: diff = PetscAbsScalar(y[i] - u[i]);
5946: tola = ts->atol;
5947: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5948: tol = tola+tolr;
5949: if(tola>0.){
5950: maxa = PetscMax(maxa,diff / tola);
5951: }
5952: if(tolr>0.){
5953: maxr = PetscMax(maxr,diff / tolr);
5954: }
5955: if(tol>0.){
5956: max = PetscMax(max,diff / tol);
5957: }
5958: }
5959: VecRestoreArrayRead(ts->vrtol,&rtol);
5960: } else { /* scalar atol, scalar rtol */
5962: for (i=0; i<n; i++) {
5963: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
5964: diff = PetscAbsScalar(y[i] - u[i]);
5965: tola = ts->atol;
5966: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5967: tol = tola+tolr;
5968: if(tola>0.){
5969: maxa = PetscMax(maxa,diff / tola);
5970: }
5971: if(tolr>0.){
5972: maxr = PetscMax(maxr,diff / tolr);
5973: }
5974: if(tol>0.){
5975: max = PetscMax(max,diff / tol);
5976: }
5977: }
5978: }
5979: VecRestoreArrayRead(U,&u);
5980: VecRestoreArrayRead(Y,&y);
5981: err_loc[0] = max;
5982: err_loc[1] = maxa;
5983: err_loc[2] = maxr;
5984: MPIU_Allreduce(err_loc,err_glb,3,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));
5985: gmax = err_glb[0];
5986: gmaxa = err_glb[1];
5987: gmaxr = err_glb[2];
5989: *norm = gmax;
5990: *norma = gmaxa;
5991: *normr = gmaxr;
5992: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
5993: if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
5994: if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
5995: return(0);
5996: }
5998: /*@
5999: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
6001: Collective on TS
6003: Input Arguments:
6004: + ts - time stepping context
6005: . U - state vector, usually ts->vec_sol
6006: . Y - state vector to be compared to U
6007: - wnormtype - norm type, either NORM_2 or NORM_INFINITY
6009: Output Arguments:
6010: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
6011: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
6012: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
6014: Options Database Keys:
6015: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
6017: Level: developer
6019: .seealso: TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2(), TSErrorWeightedENorm
6020: @*/
6021: PetscErrorCode TSErrorWeightedNorm(TS ts,Vec U,Vec Y,NormType wnormtype,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6022: {
6026: if (wnormtype == NORM_2) {
6027: TSErrorWeightedNorm2(ts,U,Y,norm,norma,normr);
6028: } else if(wnormtype == NORM_INFINITY) {
6029: TSErrorWeightedNormInfinity(ts,U,Y,norm,norma,normr);
6030: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
6031: return(0);
6032: }
6035: /*@
6036: TSErrorWeightedENorm2 - compute a weighted 2 error norm based on supplied absolute and relative tolerances
6038: Collective on TS
6040: Input Arguments:
6041: + ts - time stepping context
6042: . E - error vector
6043: . U - state vector, usually ts->vec_sol
6044: - Y - state vector, previous time step
6046: Output Arguments:
6047: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
6048: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
6049: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
6051: Level: developer
6053: .seealso: TSErrorWeightedENorm(), TSErrorWeightedENormInfinity()
6054: @*/
6055: PetscErrorCode TSErrorWeightedENorm2(TS ts,Vec E,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6056: {
6057: PetscErrorCode ierr;
6058: PetscInt i,n,N,rstart;
6059: PetscInt n_loc,na_loc,nr_loc;
6060: PetscReal n_glb,na_glb,nr_glb;
6061: const PetscScalar *e,*u,*y;
6062: PetscReal err,sum,suma,sumr,gsum,gsuma,gsumr;
6063: PetscReal tol,tola,tolr;
6064: PetscReal err_loc[6],err_glb[6];
6080: VecGetSize(E,&N);
6081: VecGetLocalSize(E,&n);
6082: VecGetOwnershipRange(E,&rstart,NULL);
6083: VecGetArrayRead(E,&e);
6084: VecGetArrayRead(U,&u);
6085: VecGetArrayRead(Y,&y);
6086: sum = 0.; n_loc = 0;
6087: suma = 0.; na_loc = 0;
6088: sumr = 0.; nr_loc = 0;
6089: if (ts->vatol && ts->vrtol) {
6090: const PetscScalar *atol,*rtol;
6091: VecGetArrayRead(ts->vatol,&atol);
6092: VecGetArrayRead(ts->vrtol,&rtol);
6093: for (i=0; i<n; i++) {
6094: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
6095: err = PetscAbsScalar(e[i]);
6096: tola = PetscRealPart(atol[i]);
6097: if(tola>0.){
6098: suma += PetscSqr(err/tola);
6099: na_loc++;
6100: }
6101: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6102: if(tolr>0.){
6103: sumr += PetscSqr(err/tolr);
6104: nr_loc++;
6105: }
6106: tol=tola+tolr;
6107: if(tol>0.){
6108: sum += PetscSqr(err/tol);
6109: n_loc++;
6110: }
6111: }
6112: VecRestoreArrayRead(ts->vatol,&atol);
6113: VecRestoreArrayRead(ts->vrtol,&rtol);
6114: } else if (ts->vatol) { /* vector atol, scalar rtol */
6115: const PetscScalar *atol;
6116: VecGetArrayRead(ts->vatol,&atol);
6117: for (i=0; i<n; i++) {
6118: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
6119: err = PetscAbsScalar(e[i]);
6120: tola = PetscRealPart(atol[i]);
6121: if(tola>0.){
6122: suma += PetscSqr(err/tola);
6123: na_loc++;
6124: }
6125: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6126: if(tolr>0.){
6127: sumr += PetscSqr(err/tolr);
6128: nr_loc++;
6129: }
6130: tol=tola+tolr;
6131: if(tol>0.){
6132: sum += PetscSqr(err/tol);
6133: n_loc++;
6134: }
6135: }
6136: VecRestoreArrayRead(ts->vatol,&atol);
6137: } else if (ts->vrtol) { /* scalar atol, vector rtol */
6138: const PetscScalar *rtol;
6139: VecGetArrayRead(ts->vrtol,&rtol);
6140: for (i=0; i<n; i++) {
6141: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
6142: err = PetscAbsScalar(e[i]);
6143: tola = ts->atol;
6144: if(tola>0.){
6145: suma += PetscSqr(err/tola);
6146: na_loc++;
6147: }
6148: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6149: if(tolr>0.){
6150: sumr += PetscSqr(err/tolr);
6151: nr_loc++;
6152: }
6153: tol=tola+tolr;
6154: if(tol>0.){
6155: sum += PetscSqr(err/tol);
6156: n_loc++;
6157: }
6158: }
6159: VecRestoreArrayRead(ts->vrtol,&rtol);
6160: } else { /* scalar atol, scalar rtol */
6161: for (i=0; i<n; i++) {
6162: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
6163: err = PetscAbsScalar(e[i]);
6164: tola = ts->atol;
6165: if(tola>0.){
6166: suma += PetscSqr(err/tola);
6167: na_loc++;
6168: }
6169: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6170: if(tolr>0.){
6171: sumr += PetscSqr(err/tolr);
6172: nr_loc++;
6173: }
6174: tol=tola+tolr;
6175: if(tol>0.){
6176: sum += PetscSqr(err/tol);
6177: n_loc++;
6178: }
6179: }
6180: }
6181: VecRestoreArrayRead(E,&e);
6182: VecRestoreArrayRead(U,&u);
6183: VecRestoreArrayRead(Y,&y);
6185: err_loc[0] = sum;
6186: err_loc[1] = suma;
6187: err_loc[2] = sumr;
6188: err_loc[3] = (PetscReal)n_loc;
6189: err_loc[4] = (PetscReal)na_loc;
6190: err_loc[5] = (PetscReal)nr_loc;
6192: MPIU_Allreduce(err_loc,err_glb,6,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));
6194: gsum = err_glb[0];
6195: gsuma = err_glb[1];
6196: gsumr = err_glb[2];
6197: n_glb = err_glb[3];
6198: na_glb = err_glb[4];
6199: nr_glb = err_glb[5];
6201: *norm = 0.;
6202: if(n_glb>0. ){*norm = PetscSqrtReal(gsum / n_glb );}
6203: *norma = 0.;
6204: if(na_glb>0.){*norma = PetscSqrtReal(gsuma / na_glb);}
6205: *normr = 0.;
6206: if(nr_glb>0.){*normr = PetscSqrtReal(gsumr / nr_glb);}
6208: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
6209: if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
6210: if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
6211: return(0);
6212: }
6214: /*@
6215: TSErrorWeightedENormInfinity - compute a weighted infinity error norm based on supplied absolute and relative tolerances
6216: Collective on TS
6218: Input Arguments:
6219: + ts - time stepping context
6220: . E - error vector
6221: . U - state vector, usually ts->vec_sol
6222: - Y - state vector, previous time step
6224: Output Arguments:
6225: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
6226: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
6227: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
6229: Level: developer
6231: .seealso: TSErrorWeightedENorm(), TSErrorWeightedENorm2()
6232: @*/
6233: PetscErrorCode TSErrorWeightedENormInfinity(TS ts,Vec E,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6234: {
6235: PetscErrorCode ierr;
6236: PetscInt i,n,N,rstart;
6237: const PetscScalar *e,*u,*y;
6238: PetscReal err,max,gmax,maxa,gmaxa,maxr,gmaxr;
6239: PetscReal tol,tola,tolr;
6240: PetscReal err_loc[3],err_glb[3];
6256: VecGetSize(E,&N);
6257: VecGetLocalSize(E,&n);
6258: VecGetOwnershipRange(E,&rstart,NULL);
6259: VecGetArrayRead(E,&e);
6260: VecGetArrayRead(U,&u);
6261: VecGetArrayRead(Y,&y);
6263: max=0.;
6264: maxa=0.;
6265: maxr=0.;
6267: if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
6268: const PetscScalar *atol,*rtol;
6269: VecGetArrayRead(ts->vatol,&atol);
6270: VecGetArrayRead(ts->vrtol,&rtol);
6272: for (i=0; i<n; i++) {
6273: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
6274: err = PetscAbsScalar(e[i]);
6275: tola = PetscRealPart(atol[i]);
6276: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6277: tol = tola+tolr;
6278: if(tola>0.){
6279: maxa = PetscMax(maxa,err / tola);
6280: }
6281: if(tolr>0.){
6282: maxr = PetscMax(maxr,err / tolr);
6283: }
6284: if(tol>0.){
6285: max = PetscMax(max,err / tol);
6286: }
6287: }
6288: VecRestoreArrayRead(ts->vatol,&atol);
6289: VecRestoreArrayRead(ts->vrtol,&rtol);
6290: } else if (ts->vatol) { /* vector atol, scalar rtol */
6291: const PetscScalar *atol;
6292: VecGetArrayRead(ts->vatol,&atol);
6293: for (i=0; i<n; i++) {
6294: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
6295: err = PetscAbsScalar(e[i]);
6296: tola = PetscRealPart(atol[i]);
6297: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6298: tol = tola+tolr;
6299: if(tola>0.){
6300: maxa = PetscMax(maxa,err / tola);
6301: }
6302: if(tolr>0.){
6303: maxr = PetscMax(maxr,err / tolr);
6304: }
6305: if(tol>0.){
6306: max = PetscMax(max,err / tol);
6307: }
6308: }
6309: VecRestoreArrayRead(ts->vatol,&atol);
6310: } else if (ts->vrtol) { /* scalar atol, vector rtol */
6311: const PetscScalar *rtol;
6312: VecGetArrayRead(ts->vrtol,&rtol);
6314: for (i=0; i<n; i++) {
6315: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
6316: err = PetscAbsScalar(e[i]);
6317: tola = ts->atol;
6318: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6319: tol = tola+tolr;
6320: if(tola>0.){
6321: maxa = PetscMax(maxa,err / tola);
6322: }
6323: if(tolr>0.){
6324: maxr = PetscMax(maxr,err / tolr);
6325: }
6326: if(tol>0.){
6327: max = PetscMax(max,err / tol);
6328: }
6329: }
6330: VecRestoreArrayRead(ts->vrtol,&rtol);
6331: } else { /* scalar atol, scalar rtol */
6333: for (i=0; i<n; i++) {
6334: SkipSmallValue(y[i],u[i],ts->adapt->ignore_max);
6335: err = PetscAbsScalar(e[i]);
6336: tola = ts->atol;
6337: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6338: tol = tola+tolr;
6339: if(tola>0.){
6340: maxa = PetscMax(maxa,err / tola);
6341: }
6342: if(tolr>0.){
6343: maxr = PetscMax(maxr,err / tolr);
6344: }
6345: if(tol>0.){
6346: max = PetscMax(max,err / tol);
6347: }
6348: }
6349: }
6350: VecRestoreArrayRead(E,&e);
6351: VecRestoreArrayRead(U,&u);
6352: VecRestoreArrayRead(Y,&y);
6353: err_loc[0] = max;
6354: err_loc[1] = maxa;
6355: err_loc[2] = maxr;
6356: MPIU_Allreduce(err_loc,err_glb,3,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));
6357: gmax = err_glb[0];
6358: gmaxa = err_glb[1];
6359: gmaxr = err_glb[2];
6361: *norm = gmax;
6362: *norma = gmaxa;
6363: *normr = gmaxr;
6364: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
6365: if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
6366: if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
6367: return(0);
6368: }
6370: /*@
6371: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
6373: Collective on TS
6375: Input Arguments:
6376: + ts - time stepping context
6377: . E - error vector
6378: . U - state vector, usually ts->vec_sol
6379: . Y - state vector, previous time step
6380: - wnormtype - norm type, either NORM_2 or NORM_INFINITY
6382: Output Arguments:
6383: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
6384: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
6385: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
6387: Options Database Keys:
6388: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
6390: Level: developer
6392: .seealso: TSErrorWeightedENormInfinity(), TSErrorWeightedENorm2(), TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2()
6393: @*/
6394: PetscErrorCode TSErrorWeightedENorm(TS ts,Vec E,Vec U,Vec Y,NormType wnormtype,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6395: {
6399: if (wnormtype == NORM_2) {
6400: TSErrorWeightedENorm2(ts,E,U,Y,norm,norma,normr);
6401: } else if(wnormtype == NORM_INFINITY) {
6402: TSErrorWeightedENormInfinity(ts,E,U,Y,norm,norma,normr);
6403: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
6404: return(0);
6405: }
6408: /*@
6409: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
6411: Logically Collective on TS
6413: Input Arguments:
6414: + ts - time stepping context
6415: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
6417: Note:
6418: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
6420: Level: intermediate
6422: .seealso: TSGetCFLTime(), TSADAPTCFL
6423: @*/
6424: PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime)
6425: {
6428: ts->cfltime_local = cfltime;
6429: ts->cfltime = -1.;
6430: return(0);
6431: }
6433: /*@
6434: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
6436: Collective on TS
6438: Input Arguments:
6439: . ts - time stepping context
6441: Output Arguments:
6442: . cfltime - maximum stable time step for forward Euler
6444: Level: advanced
6446: .seealso: TSSetCFLTimeLocal()
6447: @*/
6448: PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime)
6449: {
6453: if (ts->cfltime < 0) {
6454: MPIU_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));
6455: }
6456: *cfltime = ts->cfltime;
6457: return(0);
6458: }
6460: /*@
6461: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
6463: Input Parameters:
6464: + ts - the TS context.
6465: . xl - lower bound.
6466: - xu - upper bound.
6468: Notes:
6469: If this routine is not called then the lower and upper bounds are set to
6470: PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().
6472: Level: advanced
6474: @*/
6475: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
6476: {
6478: SNES snes;
6481: TSGetSNES(ts,&snes);
6482: SNESVISetVariableBounds(snes,xl,xu);
6483: return(0);
6484: }
6486: /*@C
6487: TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector
6488: in a time based line graph
6490: Collective on TS
6492: Input Parameters:
6493: + ts - the TS context
6494: . step - current time-step
6495: . ptime - current time
6496: . u - current solution
6497: - dctx - the TSMonitorLGCtx object that contains all the options for the monitoring, this is created with TSMonitorLGCtxCreate()
6499: Options Database:
6500: . -ts_monitor_lg_solution_variables
6502: Level: intermediate
6504: Notes:
6505: Each process in a parallel run displays its component solutions in a separate window
6507: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
6508: TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
6509: TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
6510: TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()
6511: @*/
6512: PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
6513: {
6514: PetscErrorCode ierr;
6515: TSMonitorLGCtx ctx = (TSMonitorLGCtx)dctx;
6516: const PetscScalar *yy;
6517: Vec v;
6520: if (step < 0) return(0); /* -1 indicates interpolated solution */
6521: if (!step) {
6522: PetscDrawAxis axis;
6523: PetscInt dim;
6524: PetscDrawLGGetAxis(ctx->lg,&axis);
6525: PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");
6526: if (!ctx->names) {
6527: PetscBool flg;
6528: /* user provides names of variables to plot but no names has been set so assume names are integer values */
6529: PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",&flg);
6530: if (flg) {
6531: PetscInt i,n;
6532: char **names;
6533: VecGetSize(u,&n);
6534: PetscMalloc1(n+1,&names);
6535: for (i=0; i<n; i++) {
6536: PetscMalloc1(5,&names[i]);
6537: PetscSNPrintf(names[i],5,"%D",i);
6538: }
6539: names[n] = NULL;
6540: ctx->names = names;
6541: }
6542: }
6543: if (ctx->names && !ctx->displaynames) {
6544: char **displaynames;
6545: PetscBool flg;
6546: VecGetLocalSize(u,&dim);
6547: PetscCalloc1(dim+1,&displaynames);
6548: PetscOptionsGetStringArray(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",displaynames,&dim,&flg);
6549: if (flg) {
6550: TSMonitorLGCtxSetDisplayVariables(ctx,(const char *const *)displaynames);
6551: }
6552: PetscStrArrayDestroy(&displaynames);
6553: }
6554: if (ctx->displaynames) {
6555: PetscDrawLGSetDimension(ctx->lg,ctx->ndisplayvariables);
6556: PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->displaynames);
6557: } else if (ctx->names) {
6558: VecGetLocalSize(u,&dim);
6559: PetscDrawLGSetDimension(ctx->lg,dim);
6560: PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->names);
6561: } else {
6562: VecGetLocalSize(u,&dim);
6563: PetscDrawLGSetDimension(ctx->lg,dim);
6564: }
6565: PetscDrawLGReset(ctx->lg);
6566: }
6568: if (!ctx->transform) v = u;
6569: else {(*ctx->transform)(ctx->transformctx,u,&v);}
6570: VecGetArrayRead(v,&yy);
6571: if (ctx->displaynames) {
6572: PetscInt i;
6573: for (i=0; i<ctx->ndisplayvariables; i++)
6574: ctx->displayvalues[i] = PetscRealPart(yy[ctx->displayvariables[i]]);
6575: PetscDrawLGAddCommonPoint(ctx->lg,ptime,ctx->displayvalues);
6576: } else {
6577: #if defined(PETSC_USE_COMPLEX)
6578: PetscInt i,n;
6579: PetscReal *yreal;
6580: VecGetLocalSize(v,&n);
6581: PetscMalloc1(n,&yreal);
6582: for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
6583: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
6584: PetscFree(yreal);
6585: #else
6586: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
6587: #endif
6588: }
6589: VecRestoreArrayRead(v,&yy);
6590: if (ctx->transform) {VecDestroy(&v);}
6592: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
6593: PetscDrawLGDraw(ctx->lg);
6594: PetscDrawLGSave(ctx->lg);
6595: }
6596: return(0);
6597: }
6599: /*@C
6600: TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot
6602: Collective on TS
6604: Input Parameters:
6605: + ts - the TS context
6606: - names - the names of the components, final string must be NULL
6608: Level: intermediate
6610: Notes:
6611: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
6613: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetVariableNames()
6614: @*/
6615: PetscErrorCode TSMonitorLGSetVariableNames(TS ts,const char * const *names)
6616: {
6617: PetscErrorCode ierr;
6618: PetscInt i;
6621: for (i=0; i<ts->numbermonitors; i++) {
6622: if (ts->monitor[i] == TSMonitorLGSolution) {
6623: TSMonitorLGCtxSetVariableNames((TSMonitorLGCtx)ts->monitorcontext[i],names);
6624: break;
6625: }
6626: }
6627: return(0);
6628: }
6630: /*@C
6631: TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot
6633: Collective on TS
6635: Input Parameters:
6636: + ts - the TS context
6637: - names - the names of the components, final string must be NULL
6639: Level: intermediate
6641: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGSetVariableNames()
6642: @*/
6643: PetscErrorCode TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx,const char * const *names)
6644: {
6645: PetscErrorCode ierr;
6648: PetscStrArrayDestroy(&ctx->names);
6649: PetscStrArrayallocpy(names,&ctx->names);
6650: return(0);
6651: }
6653: /*@C
6654: TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot
6656: Collective on TS
6658: Input Parameter:
6659: . ts - the TS context
6661: Output Parameter:
6662: . names - the names of the components, final string must be NULL
6664: Level: intermediate
6666: Notes:
6667: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
6669: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
6670: @*/
6671: PetscErrorCode TSMonitorLGGetVariableNames(TS ts,const char *const **names)
6672: {
6673: PetscInt i;
6676: *names = NULL;
6677: for (i=0; i<ts->numbermonitors; i++) {
6678: if (ts->monitor[i] == TSMonitorLGSolution) {
6679: TSMonitorLGCtx ctx = (TSMonitorLGCtx) ts->monitorcontext[i];
6680: *names = (const char *const *)ctx->names;
6681: break;
6682: }
6683: }
6684: return(0);
6685: }
6687: /*@C
6688: TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor
6690: Collective on TS
6692: Input Parameters:
6693: + ctx - the TSMonitorLG context
6694: - displaynames - the names of the components, final string must be NULL
6696: Level: intermediate
6698: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
6699: @*/
6700: PetscErrorCode TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx,const char * const *displaynames)
6701: {
6702: PetscInt j = 0,k;
6703: PetscErrorCode ierr;
6706: if (!ctx->names) return(0);
6707: PetscStrArrayDestroy(&ctx->displaynames);
6708: PetscStrArrayallocpy(displaynames,&ctx->displaynames);
6709: while (displaynames[j]) j++;
6710: ctx->ndisplayvariables = j;
6711: PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvariables);
6712: PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvalues);
6713: j = 0;
6714: while (displaynames[j]) {
6715: k = 0;
6716: while (ctx->names[k]) {
6717: PetscBool flg;
6718: PetscStrcmp(displaynames[j],ctx->names[k],&flg);
6719: if (flg) {
6720: ctx->displayvariables[j] = k;
6721: break;
6722: }
6723: k++;
6724: }
6725: j++;
6726: }
6727: return(0);
6728: }
6730: /*@C
6731: TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor
6733: Collective on TS
6735: Input Parameters:
6736: + ts - the TS context
6737: - displaynames - the names of the components, final string must be NULL
6739: Notes:
6740: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
6742: Level: intermediate
6744: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
6745: @*/
6746: PetscErrorCode TSMonitorLGSetDisplayVariables(TS ts,const char * const *displaynames)
6747: {
6748: PetscInt i;
6749: PetscErrorCode ierr;
6752: for (i=0; i<ts->numbermonitors; i++) {
6753: if (ts->monitor[i] == TSMonitorLGSolution) {
6754: TSMonitorLGCtxSetDisplayVariables((TSMonitorLGCtx)ts->monitorcontext[i],displaynames);
6755: break;
6756: }
6757: }
6758: return(0);
6759: }
6761: /*@C
6762: TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed
6764: Collective on TS
6766: Input Parameters:
6767: + ts - the TS context
6768: . transform - the transform function
6769: . destroy - function to destroy the optional context
6770: - ctx - optional context used by transform function
6772: Notes:
6773: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored
6775: Level: intermediate
6777: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGCtxSetTransform()
6778: @*/
6779: PetscErrorCode TSMonitorLGSetTransform(TS ts,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
6780: {
6781: PetscInt i;
6782: PetscErrorCode ierr;
6785: for (i=0; i<ts->numbermonitors; i++) {
6786: if (ts->monitor[i] == TSMonitorLGSolution) {
6787: TSMonitorLGCtxSetTransform((TSMonitorLGCtx)ts->monitorcontext[i],transform,destroy,tctx);
6788: }
6789: }
6790: return(0);
6791: }
6793: /*@C
6794: TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed
6796: Collective on TSLGCtx
6798: Input Parameters:
6799: + ts - the TS context
6800: . transform - the transform function
6801: . destroy - function to destroy the optional context
6802: - ctx - optional context used by transform function
6804: Level: intermediate
6806: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGSetTransform()
6807: @*/
6808: PetscErrorCode TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
6809: {
6811: ctx->transform = transform;
6812: ctx->transformdestroy = destroy;
6813: ctx->transformctx = tctx;
6814: return(0);
6815: }
6817: /*@C
6818: TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the error
6819: in a time based line graph
6821: Collective on TS
6823: Input Parameters:
6824: + ts - the TS context
6825: . step - current time-step
6826: . ptime - current time
6827: . u - current solution
6828: - dctx - TSMonitorLGCtx object created with TSMonitorLGCtxCreate()
6830: Level: intermediate
6832: Notes:
6833: Each process in a parallel run displays its component errors in a separate window
6835: The user must provide the solution using TSSetSolutionFunction() to use this monitor.
6837: Options Database Keys:
6838: . -ts_monitor_lg_error - create a graphical monitor of error history
6840: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
6841: @*/
6842: PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
6843: {
6844: PetscErrorCode ierr;
6845: TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy;
6846: const PetscScalar *yy;
6847: Vec y;
6850: if (!step) {
6851: PetscDrawAxis axis;
6852: PetscInt dim;
6853: PetscDrawLGGetAxis(ctx->lg,&axis);
6854: PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Error");
6855: VecGetLocalSize(u,&dim);
6856: PetscDrawLGSetDimension(ctx->lg,dim);
6857: PetscDrawLGReset(ctx->lg);
6858: }
6859: VecDuplicate(u,&y);
6860: TSComputeSolutionFunction(ts,ptime,y);
6861: VecAXPY(y,-1.0,u);
6862: VecGetArrayRead(y,&yy);
6863: #if defined(PETSC_USE_COMPLEX)
6864: {
6865: PetscReal *yreal;
6866: PetscInt i,n;
6867: VecGetLocalSize(y,&n);
6868: PetscMalloc1(n,&yreal);
6869: for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
6870: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
6871: PetscFree(yreal);
6872: }
6873: #else
6874: PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
6875: #endif
6876: VecRestoreArrayRead(y,&yy);
6877: VecDestroy(&y);
6878: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
6879: PetscDrawLGDraw(ctx->lg);
6880: PetscDrawLGSave(ctx->lg);
6881: }
6882: return(0);
6883: }
6885: /*@C
6886: TSMonitorSPSwarmSolution - Graphically displays phase plots of DMSwarm particles on a scatter plot
6888: Input Parameters:
6889: + ts - the TS context
6890: . step - current time-step
6891: . ptime - current time
6892: . u - current solution
6893: - dctx - the TSMonitorSPCtx object that contains all the options for the monitoring, this is created with TSMonitorSPCtxCreate()
6895: Options Database:
6896: . -ts_monitor_sp_swarm
6898: Level: intermediate
6900: @*/
6901: PetscErrorCode TSMonitorSPSwarmSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
6902: {
6903: PetscErrorCode ierr;
6904: TSMonitorSPCtx ctx = (TSMonitorSPCtx)dctx;
6905: const PetscScalar *yy;
6906: PetscReal *y,*x;
6907: PetscInt Np, p, dim=2;
6908: DM dm;
6912: if (step < 0) return(0); /* -1 indicates interpolated solution */
6913: if (!step) {
6914: PetscDrawAxis axis;
6915: PetscDrawSPGetAxis(ctx->sp,&axis);
6916: PetscDrawAxisSetLabels(axis,"Particles","X","Y");
6917: PetscDrawAxisSetLimits(axis, -5, 5, -5, 5);
6918: PetscDrawAxisSetHoldLimits(axis, PETSC_TRUE);
6919: TSGetDM(ts, &dm);
6920: DMGetDimension(dm, &dim);
6921: if(dim!=2) SETERRQ(PETSC_COMM_SELF, ierr, "Dimensions improper for monitor arguments! Current support: two dimensions.");
6922: VecGetLocalSize(u, &Np);
6923: Np /= 2*dim;
6924: PetscDrawSPSetDimension(ctx->sp, Np);
6925: PetscDrawSPReset(ctx->sp);
6926: }
6928: VecGetLocalSize(u, &Np);
6929: Np /= 2*dim;
6930: VecGetArrayRead(u,&yy);
6931: PetscMalloc2(Np, &x, Np, &y);
6932: /* get points from solution vector */
6933: for (p=0; p<Np; ++p){
6934: x[p] = PetscRealPart(yy[2*dim*p]);
6935: y[p] = PetscRealPart(yy[2*dim*p+1]);
6936: }
6937: VecRestoreArrayRead(u,&yy);
6939: if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
6940: PetscDrawSPAddPoint(ctx->sp,x,y);
6941: PetscDrawSPDraw(ctx->sp,PETSC_FALSE);
6942: PetscDrawSPSave(ctx->sp);
6943: }
6945: PetscFree2(x, y);
6947: return(0);
6948: }
6952: /*@C
6953: TSMonitorError - Monitors progress of the TS solvers by printing the 2 norm of the error at each timestep
6955: Collective on TS
6957: Input Parameters:
6958: + ts - the TS context
6959: . step - current time-step
6960: . ptime - current time
6961: . u - current solution
6962: - dctx - unused context
6964: Level: intermediate
6966: The user must provide the solution using TSSetSolutionFunction() to use this monitor.
6968: Options Database Keys:
6969: . -ts_monitor_error - create a graphical monitor of error history
6971: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
6972: @*/
6973: PetscErrorCode TSMonitorError(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscViewerAndFormat *vf)
6974: {
6975: PetscErrorCode ierr;
6976: Vec y;
6977: PetscReal nrm;
6978: PetscBool flg;
6981: VecDuplicate(u,&y);
6982: TSComputeSolutionFunction(ts,ptime,y);
6983: VecAXPY(y,-1.0,u);
6984: PetscObjectTypeCompare((PetscObject)vf->viewer,PETSCVIEWERASCII,&flg);
6985: if (flg) {
6986: VecNorm(y,NORM_2,&nrm);
6987: PetscViewerASCIIPrintf(vf->viewer,"2-norm of error %g\n",(double)nrm);
6988: }
6989: PetscObjectTypeCompare((PetscObject)vf->viewer,PETSCVIEWERDRAW,&flg);
6990: if (flg) {
6991: VecView(y,vf->viewer);
6992: }
6993: VecDestroy(&y);
6994: return(0);
6995: }
6997: PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
6998: {
6999: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
7000: PetscReal x = ptime,y;
7002: PetscInt its;
7005: if (n < 0) return(0); /* -1 indicates interpolated solution */
7006: if (!n) {
7007: PetscDrawAxis axis;
7008: PetscDrawLGGetAxis(ctx->lg,&axis);
7009: PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");
7010: PetscDrawLGReset(ctx->lg);
7011: ctx->snes_its = 0;
7012: }
7013: TSGetSNESIterations(ts,&its);
7014: y = its - ctx->snes_its;
7015: PetscDrawLGAddPoint(ctx->lg,&x,&y);
7016: if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
7017: PetscDrawLGDraw(ctx->lg);
7018: PetscDrawLGSave(ctx->lg);
7019: }
7020: ctx->snes_its = its;
7021: return(0);
7022: }
7024: PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
7025: {
7026: TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
7027: PetscReal x = ptime,y;
7029: PetscInt its;
7032: if (n < 0) return(0); /* -1 indicates interpolated solution */
7033: if (!n) {
7034: PetscDrawAxis axis;
7035: PetscDrawLGGetAxis(ctx->lg,&axis);
7036: PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");
7037: PetscDrawLGReset(ctx->lg);
7038: ctx->ksp_its = 0;
7039: }
7040: TSGetKSPIterations(ts,&its);
7041: y = its - ctx->ksp_its;
7042: PetscDrawLGAddPoint(ctx->lg,&x,&y);
7043: if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
7044: PetscDrawLGDraw(ctx->lg);
7045: PetscDrawLGSave(ctx->lg);
7046: }
7047: ctx->ksp_its = its;
7048: return(0);
7049: }
7051: /*@
7052: TSComputeLinearStability - computes the linear stability function at a point
7054: Collective on TS
7056: Input Parameters:
7057: + ts - the TS context
7058: - xr,xi - real and imaginary part of input arguments
7060: Output Parameters:
7061: . yr,yi - real and imaginary part of function value
7063: Level: developer
7065: .seealso: TSSetRHSFunction(), TSComputeIFunction()
7066: @*/
7067: PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
7068: {
7073: if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method");
7074: (*ts->ops->linearstability)(ts,xr,xi,yr,yi);
7075: return(0);
7076: }
7078: /* ------------------------------------------------------------------------*/
7079: /*@C
7080: TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope()
7082: Collective on TS
7084: Input Parameters:
7085: . ts - the ODE solver object
7087: Output Parameter:
7088: . ctx - the context
7090: Level: intermediate
7092: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()
7094: @*/
7095: PetscErrorCode TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx)
7096: {
7100: PetscNew(ctx);
7101: return(0);
7102: }
7104: /*@C
7105: TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution
7107: Collective on TS
7109: Input Parameters:
7110: + ts - the TS context
7111: . step - current time-step
7112: . ptime - current time
7113: . u - current solution
7114: - dctx - the envelope context
7116: Options Database:
7117: . -ts_monitor_envelope
7119: Level: intermediate
7121: Notes:
7122: after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope
7124: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxCreate()
7125: @*/
7126: PetscErrorCode TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
7127: {
7128: PetscErrorCode ierr;
7129: TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dctx;
7132: if (!ctx->max) {
7133: VecDuplicate(u,&ctx->max);
7134: VecDuplicate(u,&ctx->min);
7135: VecCopy(u,ctx->max);
7136: VecCopy(u,ctx->min);
7137: } else {
7138: VecPointwiseMax(ctx->max,u,ctx->max);
7139: VecPointwiseMin(ctx->min,u,ctx->min);
7140: }
7141: return(0);
7142: }
7144: /*@C
7145: TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution
7147: Collective on TS
7149: Input Parameter:
7150: . ts - the TS context
7152: Output Parameter:
7153: + max - the maximum values
7154: - min - the minimum values
7156: Notes:
7157: If the TS does not have a TSMonitorEnvelopeCtx associated with it then this function is ignored
7159: Level: intermediate
7161: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
7162: @*/
7163: PetscErrorCode TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min)
7164: {
7165: PetscInt i;
7168: if (max) *max = NULL;
7169: if (min) *min = NULL;
7170: for (i=0; i<ts->numbermonitors; i++) {
7171: if (ts->monitor[i] == TSMonitorEnvelope) {
7172: TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx) ts->monitorcontext[i];
7173: if (max) *max = ctx->max;
7174: if (min) *min = ctx->min;
7175: break;
7176: }
7177: }
7178: return(0);
7179: }
7181: /*@C
7182: TSMonitorEnvelopeCtxDestroy - Destroys a context that was created with TSMonitorEnvelopeCtxCreate().
7184: Collective on TSMonitorEnvelopeCtx
7186: Input Parameter:
7187: . ctx - the monitor context
7189: Level: intermediate
7191: .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep()
7192: @*/
7193: PetscErrorCode TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx)
7194: {
7198: VecDestroy(&(*ctx)->min);
7199: VecDestroy(&(*ctx)->max);
7200: PetscFree(*ctx);
7201: return(0);
7202: }
7204: /*@
7205: TSRestartStep - Flags the solver to restart the next step
7207: Collective on TS
7209: Input Parameter:
7210: . ts - the TS context obtained from TSCreate()
7212: Level: advanced
7214: Notes:
7215: Multistep methods like BDF or Runge-Kutta methods with FSAL property require restarting the solver in the event of
7216: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
7217: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
7218: the sake of correctness and maximum safety, users are expected to call TSRestart() whenever they introduce
7219: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
7220: discontinuous source terms).
7222: .seealso: TSSolve(), TSSetPreStep(), TSSetPostStep()
7223: @*/
7224: PetscErrorCode TSRestartStep(TS ts)
7225: {
7228: ts->steprestart = PETSC_TRUE;
7229: return(0);
7230: }
7232: /*@
7233: TSRollBack - Rolls back one time step
7235: Collective on TS
7237: Input Parameter:
7238: . ts - the TS context obtained from TSCreate()
7240: Level: advanced
7242: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate()
7243: @*/
7244: PetscErrorCode TSRollBack(TS ts)
7245: {
7250: if (ts->steprollback) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TSRollBack already called");
7251: if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name);
7252: (*ts->ops->rollback)(ts);
7253: ts->time_step = ts->ptime - ts->ptime_prev;
7254: ts->ptime = ts->ptime_prev;
7255: ts->ptime_prev = ts->ptime_prev_rollback;
7256: ts->steps--;
7257: ts->steprollback = PETSC_TRUE;
7258: return(0);
7259: }
7261: /*@
7262: TSGetStages - Get the number of stages and stage values
7264: Input Parameter:
7265: . ts - the TS context obtained from TSCreate()
7267: Output Parameters:
7268: + ns - the number of stages
7269: - Y - the current stage vectors
7271: Level: advanced
7273: Notes: Both ns and Y can be NULL.
7275: .seealso: TSCreate()
7276: @*/
7277: PetscErrorCode TSGetStages(TS ts,PetscInt *ns,Vec **Y)
7278: {
7285: if (!ts->ops->getstages) {
7286: if (ns) *ns = 0;
7287: if (Y) *Y = NULL;
7288: } else {
7289: (*ts->ops->getstages)(ts,ns,Y);
7290: }
7291: return(0);
7292: }
7294: /*@C
7295: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
7297: Collective on SNES
7299: Input Parameters:
7300: + ts - the TS context
7301: . t - current timestep
7302: . U - state vector
7303: . Udot - time derivative of state vector
7304: . shift - shift to apply, see note below
7305: - ctx - an optional user context
7307: Output Parameters:
7308: + J - Jacobian matrix (not altered in this routine)
7309: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)
7311: Level: intermediate
7313: Notes:
7314: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
7316: dF/dU + shift*dF/dUdot
7318: Most users should not need to explicitly call this routine, as it
7319: is used internally within the nonlinear solvers.
7321: This will first try to get the coloring from the DM. If the DM type has no coloring
7322: routine, then it will try to get the coloring from the matrix. This requires that the
7323: matrix have nonzero entries precomputed.
7325: .seealso: TSSetIJacobian(), MatFDColoringCreate(), MatFDColoringSetFunction()
7326: @*/
7327: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat J,Mat B,void *ctx)
7328: {
7329: SNES snes;
7330: MatFDColoring color;
7331: PetscBool hascolor, matcolor = PETSC_FALSE;
7335: PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject) ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL);
7336: PetscObjectQuery((PetscObject) B, "TSMatFDColoring", (PetscObject *) &color);
7337: if (!color) {
7338: DM dm;
7339: ISColoring iscoloring;
7341: TSGetDM(ts, &dm);
7342: DMHasColoring(dm, &hascolor);
7343: if (hascolor && !matcolor) {
7344: DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);
7345: MatFDColoringCreate(B, iscoloring, &color);
7346: MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);
7347: MatFDColoringSetFromOptions(color);
7348: MatFDColoringSetUp(B, iscoloring, color);
7349: ISColoringDestroy(&iscoloring);
7350: } else {
7351: MatColoring mc;
7353: MatColoringCreate(B, &mc);
7354: MatColoringSetDistance(mc, 2);
7355: MatColoringSetType(mc, MATCOLORINGSL);
7356: MatColoringSetFromOptions(mc);
7357: MatColoringApply(mc, &iscoloring);
7358: MatColoringDestroy(&mc);
7359: MatFDColoringCreate(B, iscoloring, &color);
7360: MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);
7361: MatFDColoringSetFromOptions(color);
7362: MatFDColoringSetUp(B, iscoloring, color);
7363: ISColoringDestroy(&iscoloring);
7364: }
7365: PetscObjectCompose((PetscObject) B, "TSMatFDColoring", (PetscObject) color);
7366: PetscObjectDereference((PetscObject) color);
7367: }
7368: TSGetSNES(ts, &snes);
7369: MatFDColoringApply(B, color, U, snes);
7370: if (J != B) {
7371: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
7372: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
7373: }
7374: return(0);
7375: }
7377: /*@
7378: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
7380: Input Parameters:
7381: + ts - the TS context
7382: - func - function called within TSFunctionDomainError
7384: Calling sequence of func:
7385: $ PetscErrorCode func(TS ts,PetscReal time,Vec state,PetscBool reject)
7387: + ts - the TS context
7388: . time - the current time (of the stage)
7389: . state - the state to check if it is valid
7390: - reject - (output parameter) PETSC_FALSE if the state is acceptable, PETSC_TRUE if not acceptable
7392: Level: intermediate
7394: Notes:
7395: If an implicit ODE solver is being used then, in addition to providing this routine, the
7396: user's code should call SNESSetFunctionDomainError() when domain errors occur during
7397: function evaluations where the functions are provided by TSSetIFunction() or TSSetRHSFunction().
7398: Use TSGetSNES() to obtain the SNES object
7400: Developer Notes:
7401: The naming of this function is inconsistent with the SNESSetFunctionDomainError()
7402: since one takes a function pointer and the other does not.
7404: .seealso: TSAdaptCheckStage(), TSFunctionDomainError(), SNESSetFunctionDomainError(), TSGetSNES()
7405: @*/
7407: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS,PetscReal,Vec,PetscBool*))
7408: {
7411: ts->functiondomainerror = func;
7412: return(0);
7413: }
7415: /*@
7416: TSFunctionDomainError - Checks if the current state is valid
7418: Input Parameters:
7419: + ts - the TS context
7420: . stagetime - time of the simulation
7421: - Y - state vector to check.
7423: Output Parameter:
7424: . accept - Set to PETSC_FALSE if the current state vector is valid.
7426: Note:
7427: This function is called by the TS integration routines and calls the user provided function (set with TSSetFunctionDomainError())
7428: to check if the current state is valid.
7430: Level: developer
7432: .seealso: TSSetFunctionDomainError()
7433: @*/
7434: PetscErrorCode TSFunctionDomainError(TS ts,PetscReal stagetime,Vec Y,PetscBool* accept)
7435: {
7438: *accept = PETSC_TRUE;
7439: if (ts->functiondomainerror) {
7440: PetscStackCallStandard((*ts->functiondomainerror),(ts,stagetime,Y,accept));
7441: }
7442: return(0);
7443: }
7445: /*@C
7446: TSClone - This function clones a time step object.
7448: Collective
7450: Input Parameter:
7451: . tsin - The input TS
7453: Output Parameter:
7454: . tsout - The output TS (cloned)
7456: Notes:
7457: 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.
7459: 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);
7461: Level: developer
7463: .seealso: TSCreate(), TSSetType(), TSSetUp(), TSDestroy(), TSSetProblemType()
7464: @*/
7465: PetscErrorCode TSClone(TS tsin, TS *tsout)
7466: {
7467: TS t;
7469: SNES snes_start;
7470: DM dm;
7471: TSType type;
7475: *tsout = NULL;
7477: PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView);
7479: /* General TS description */
7480: t->numbermonitors = 0;
7481: t->setupcalled = 0;
7482: t->ksp_its = 0;
7483: t->snes_its = 0;
7484: t->nwork = 0;
7485: t->rhsjacobian.time = PETSC_MIN_REAL;
7486: t->rhsjacobian.scale = 1.;
7487: t->ijacobian.shift = 1.;
7489: TSGetSNES(tsin,&snes_start);
7490: TSSetSNES(t,snes_start);
7492: TSGetDM(tsin,&dm);
7493: TSSetDM(t,dm);
7495: t->adapt = tsin->adapt;
7496: PetscObjectReference((PetscObject)t->adapt);
7498: t->trajectory = tsin->trajectory;
7499: PetscObjectReference((PetscObject)t->trajectory);
7501: t->event = tsin->event;
7502: if (t->event) t->event->refct++;
7504: t->problem_type = tsin->problem_type;
7505: t->ptime = tsin->ptime;
7506: t->ptime_prev = tsin->ptime_prev;
7507: t->time_step = tsin->time_step;
7508: t->max_time = tsin->max_time;
7509: t->steps = tsin->steps;
7510: t->max_steps = tsin->max_steps;
7511: t->equation_type = tsin->equation_type;
7512: t->atol = tsin->atol;
7513: t->rtol = tsin->rtol;
7514: t->max_snes_failures = tsin->max_snes_failures;
7515: t->max_reject = tsin->max_reject;
7516: t->errorifstepfailed = tsin->errorifstepfailed;
7518: TSGetType(tsin,&type);
7519: TSSetType(t,type);
7521: t->vec_sol = NULL;
7523: t->cfltime = tsin->cfltime;
7524: t->cfltime_local = tsin->cfltime_local;
7525: t->exact_final_time = tsin->exact_final_time;
7527: PetscMemcpy(t->ops,tsin->ops,sizeof(struct _TSOps));
7529: if (((PetscObject)tsin)->fortran_func_pointers) {
7530: PetscInt i;
7531: PetscMalloc((10)*sizeof(void(*)(void)),&((PetscObject)t)->fortran_func_pointers);
7532: for (i=0; i<10; i++) {
7533: ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
7534: }
7535: }
7536: *tsout = t;
7537: return(0);
7538: }
7540: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void* ctx,Vec x,Vec y)
7541: {
7543: TS ts = (TS) ctx;
7546: TSComputeRHSFunction(ts,0,x,y);
7547: return(0);
7548: }
7550: /*@
7551: TSRHSJacobianTest - Compares the multiply routine provided to the MATSHELL with differencing on the TS given RHS function.
7553: Logically Collective on TS
7555: Input Parameters:
7556: TS - the time stepping routine
7558: Output Parameter:
7559: . flg - PETSC_TRUE if the multiply is likely correct
7561: Options Database:
7562: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
7564: Level: advanced
7566: Notes:
7567: This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian
7569: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose(), TSRHSJacobianTestTranspose()
7570: @*/
7571: PetscErrorCode TSRHSJacobianTest(TS ts,PetscBool *flg)
7572: {
7573: Mat J,B;
7575: TSRHSJacobian func;
7576: void* ctx;
7579: TSGetRHSJacobian(ts,&J,&B,&func,&ctx);
7580: (*func)(ts,0.0,ts->vec_sol,J,B,ctx);
7581: MatShellTestMult(J,RHSWrapperFunction_TSRHSJacobianTest,ts->vec_sol,ts,flg);
7582: return(0);
7583: }
7585: /*@C
7586: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the MATSHELL with differencing on the TS given RHS function.
7588: Logically Collective on TS
7590: Input Parameters:
7591: TS - the time stepping routine
7593: Output Parameter:
7594: . flg - PETSC_TRUE if the multiply is likely correct
7596: Options Database:
7597: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
7599: Notes:
7600: This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian
7602: Level: advanced
7604: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose(), TSRHSJacobianTest()
7605: @*/
7606: PetscErrorCode TSRHSJacobianTestTranspose(TS ts,PetscBool *flg)
7607: {
7608: Mat J,B;
7610: void *ctx;
7611: TSRHSJacobian func;
7614: TSGetRHSJacobian(ts,&J,&B,&func,&ctx);
7615: (*func)(ts,0.0,ts->vec_sol,J,B,ctx);
7616: MatShellTestMultTranspose(J,RHSWrapperFunction_TSRHSJacobianTest,ts->vec_sol,ts,flg);
7617: return(0);
7618: }
7620: /*@
7621: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
7623: Logically collective
7625: Input Parameter:
7626: + ts - timestepping context
7627: - use_splitrhsfunction - PETSC_TRUE indicates that the split RHSFunction will be used
7629: Options Database:
7630: . -ts_use_splitrhsfunction - <true,false>
7632: Notes:
7633: This is only useful for multirate methods
7635: Level: intermediate
7637: .seealso: TSGetUseSplitRHSFunction()
7638: @*/
7639: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
7640: {
7643: ts->use_splitrhsfunction = use_splitrhsfunction;
7644: return(0);
7645: }
7647: /*@
7648: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
7650: Not collective
7652: Input Parameter:
7653: . ts - timestepping context
7655: Output Parameter:
7656: . use_splitrhsfunction - PETSC_TRUE indicates that the split RHSFunction will be used
7658: Level: intermediate
7660: .seealso: TSSetUseSplitRHSFunction()
7661: @*/
7662: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
7663: {
7666: *use_splitrhsfunction = ts->use_splitrhsfunction;
7667: return(0);
7668: }