Actual source code: ts.c
2: #include <private/tsimpl.h> /*I "petscts.h" I*/
4: /* Logging support */
5: PetscClassId TS_CLASSID;
6: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;
10: /*
11: TSSetTypeFromOptions - Sets the type of ts from user options.
13: Collective on TS
15: Input Parameter:
16: . ts - The ts
18: Level: intermediate
20: .keywords: TS, set, options, database, type
21: .seealso: TSSetFromOptions(), TSSetType()
22: */
23: static PetscErrorCode TSSetTypeFromOptions(TS ts)
24: {
25: PetscBool opt;
26: const char *defaultType;
27: char typeName[256];
31: if (((PetscObject)ts)->type_name) {
32: defaultType = ((PetscObject)ts)->type_name;
33: } else {
34: defaultType = TSEULER;
35: }
37: if (!TSRegisterAllCalled) {TSRegisterAll(PETSC_NULL);}
38: PetscOptionsList("-ts_type", "TS method"," TSSetType", TSList, defaultType, typeName, 256, &opt);
39: if (opt) {
40: TSSetType(ts, typeName);
41: } else {
42: TSSetType(ts, defaultType);
43: }
44: return(0);
45: }
49: /*@
50: TSSetFromOptions - Sets various TS parameters from user options.
52: Collective on TS
54: Input Parameter:
55: . ts - the TS context obtained from TSCreate()
57: Options Database Keys:
58: + -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSGL, TSSSP
59: . -ts_max_steps maxsteps - maximum number of time-steps to take
60: . -ts_final_time time - maximum time to compute to
61: . -ts_dt dt - initial time step
62: . -ts_monitor - print information at each timestep
63: - -ts_monitor_draw - plot information at each timestep
65: Level: beginner
67: .keywords: TS, timestep, set, options, database
69: .seealso: TSGetType()
70: @*/
71: PetscErrorCode TSSetFromOptions(TS ts)
72: {
73: PetscBool opt,flg;
75: PetscViewer monviewer;
76: char monfilename[PETSC_MAX_PATH_LEN];
77: SNES snes;
78: TSAdapt adapt;
82: PetscObjectOptionsBegin((PetscObject)ts);
83: /* Handle TS type options */
84: TSSetTypeFromOptions(ts);
86: /* Handle generic TS options */
87: PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetDuration",ts->max_steps,&ts->max_steps,PETSC_NULL);
88: PetscOptionsReal("-ts_final_time","Time to run to","TSSetDuration",ts->max_time,&ts->max_time,PETSC_NULL);
89: PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,PETSC_NULL);
90: PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&ts->time_step,PETSC_NULL);
91: opt = ts->exact_final_time == PETSC_DECIDE ? PETSC_FALSE : (PetscBool)ts->exact_final_time;
92: PetscOptionsBool("-ts_exact_final_time","Interpolate output to stop exactly at the final time","TSSetExactFinalTime",opt,&opt,&flg);
93: if (flg) {TSSetExactFinalTime(ts,opt);}
94: PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","",ts->max_snes_failures,&ts->max_snes_failures,PETSC_NULL);
95: PetscOptionsInt("-ts_max_reject","Maximum number of step rejections","",ts->max_reject,&ts->max_reject,PETSC_NULL);
96: PetscOptionsBool("-ts_error_if_step_failed","Error if no step succeeds","",ts->errorifstepfailed,&ts->errorifstepfailed,PETSC_NULL);
97: PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,PETSC_NULL);
98: PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,PETSC_NULL);
100: /* Monitor options */
101: PetscOptionsString("-ts_monitor","Monitor timestep size","TSMonitorDefault","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
102: if (flg) {
103: PetscViewerASCIIOpen(((PetscObject)ts)->comm,monfilename,&monviewer);
104: TSMonitorSet(ts,TSMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
105: }
106: PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
107: if (flg) {PetscPythonMonitorSet((PetscObject)ts,monfilename);}
109: opt = PETSC_FALSE;
110: PetscOptionsBool("-ts_monitor_draw","Monitor timestep size graphically","TSMonitorLG",opt,&opt,PETSC_NULL);
111: if (opt) {
112: TSMonitorSet(ts,TSMonitorLG,PETSC_NULL,PETSC_NULL);
113: }
114: opt = PETSC_FALSE;
115: PetscOptionsBool("-ts_monitor_solution","Monitor solution graphically","TSMonitorSolution",opt,&opt,PETSC_NULL);
116: if (opt) {
117: void *ctx;
118: TSMonitorSolutionCreate(ts,PETSC_NULL,&ctx);
119: TSMonitorSet(ts,TSMonitorSolution,ctx,TSMonitorSolutionDestroy);
120: }
121: opt = PETSC_FALSE;
122: PetscOptionsString("-ts_monitor_solution_binary","Save each solution to a binary file","TSMonitorSolutionBinary",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
123: if (flg) {
124: PetscViewer ctx;
125: if (monfilename[0]) {
126: PetscViewerBinaryOpen(((PetscObject)ts)->comm,monfilename,FILE_MODE_WRITE,&ctx);
127: } else {
128: ctx = PETSC_VIEWER_BINARY_(((PetscObject)ts)->comm);
129: }
130: TSMonitorSet(ts,TSMonitorSolutionBinary,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);
131: }
132: opt = PETSC_FALSE;
133: 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);
134: if (flg) {
135: const char *ptr,*ptr2;
136: char *filetemplate;
137: if (!monfilename[0]) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
138: /* Do some cursory validation of the input. */
139: PetscStrstr(monfilename,"%",(char**)&ptr);
140: if (!ptr) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
141: for (ptr++ ; ptr && *ptr; ptr++) {
142: PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
143: if (!ptr2 && (*ptr < '0' || '9' < *ptr)) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03D.vts");
144: if (ptr2) break;
145: }
146: PetscStrallocpy(monfilename,&filetemplate);
147: TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);
148: }
150: TSGetAdapt(ts,&adapt);
151: TSAdaptSetFromOptions(adapt);
153: TSGetSNES(ts,&snes);
154: if (ts->problem_type == TS_LINEAR) {SNESSetType(snes,SNESKSPONLY);}
156: /* Handle specific TS options */
157: if (ts->ops->setfromoptions) {
158: (*ts->ops->setfromoptions)(ts);
159: }
161: /* process any options handlers added with PetscObjectAddOptionsHandler() */
162: PetscObjectProcessOptionsHandlers((PetscObject)ts);
163: PetscOptionsEnd();
164: return(0);
165: }
170: /*@
171: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
172: set with TSSetRHSJacobian().
174: Collective on TS and Vec
176: Input Parameters:
177: + ts - the TS context
178: . t - current timestep
179: - x - input vector
181: Output Parameters:
182: + A - Jacobian matrix
183: . B - optional preconditioning matrix
184: - flag - flag indicating matrix structure
186: Notes:
187: Most users should not need to explicitly call this routine, as it
188: is used internally within the nonlinear solvers.
190: See KSPSetOperators() for important information about setting the
191: flag parameter.
193: Level: developer
195: .keywords: SNES, compute, Jacobian, matrix
197: .seealso: TSSetRHSJacobian(), KSPSetOperators()
198: @*/
199: PetscErrorCode TSComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat *A,Mat *B,MatStructure *flg)
200: {
202: PetscInt Xstate;
208: PetscObjectStateQuery((PetscObject)X,&Xstate);
209: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.X == X && ts->rhsjacobian.Xstate == Xstate))) {
210: *flg = ts->rhsjacobian.mstructure;
211: return(0);
212: }
214: if (!ts->userops->rhsjacobian && !ts->userops->ijacobian) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
216: if (ts->userops->rhsjacobian) {
217: PetscLogEventBegin(TS_JacobianEval,ts,X,*A,*B);
218: *flg = DIFFERENT_NONZERO_PATTERN;
219: PetscStackPush("TS user Jacobian function");
220: (*ts->userops->rhsjacobian)(ts,t,X,A,B,flg,ts->jacP);
221: PetscStackPop;
222: PetscLogEventEnd(TS_JacobianEval,ts,X,*A,*B);
223: /* make sure user returned a correct Jacobian and preconditioner */
226: } else {
227: MatZeroEntries(*A);
228: if (*A != *B) {MatZeroEntries(*B);}
229: *flg = SAME_NONZERO_PATTERN;
230: }
231: ts->rhsjacobian.time = t;
232: ts->rhsjacobian.X = X;
233: PetscObjectStateQuery((PetscObject)X,&ts->rhsjacobian.Xstate);
234: ts->rhsjacobian.mstructure = *flg;
235: return(0);
236: }
240: /*@
241: TSComputeRHSFunction - Evaluates the right-hand-side function.
243: Collective on TS and Vec
245: Input Parameters:
246: + ts - the TS context
247: . t - current time
248: - x - state vector
250: Output Parameter:
251: . y - right hand side
253: Note:
254: Most users should not need to explicitly call this routine, as it
255: is used internally within the nonlinear solvers.
257: Level: developer
259: .keywords: TS, compute
261: .seealso: TSSetRHSFunction(), TSComputeIFunction()
262: @*/
263: PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec x,Vec y)
264: {
272: if (!ts->userops->rhsfunction && !ts->userops->ifunction) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
274: PetscLogEventBegin(TS_FunctionEval,ts,x,y,0);
275: if (ts->userops->rhsfunction) {
276: PetscStackPush("TS user right-hand-side function");
277: (*ts->userops->rhsfunction)(ts,t,x,y,ts->funP);
278: PetscStackPop;
279: } else {
280: VecZeroEntries(y);
281: }
283: PetscLogEventEnd(TS_FunctionEval,ts,x,y,0);
284: return(0);
285: }
289: static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs)
290: {
291: Vec F;
295: *Frhs = PETSC_NULL;
296: TSGetIFunction(ts,&F,PETSC_NULL,PETSC_NULL);
297: if (!ts->Frhs) {
298: VecDuplicate(F,&ts->Frhs);
299: }
300: *Frhs = ts->Frhs;
301: return(0);
302: }
306: static PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs)
307: {
308: Mat A,B;
312: TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);
313: if (Arhs) {
314: if (!ts->Arhs) {
315: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);
316: }
317: *Arhs = ts->Arhs;
318: }
319: if (Brhs) {
320: if (!ts->Brhs) {
321: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);
322: }
323: *Brhs = ts->Brhs;
324: }
325: return(0);
326: }
330: /*@
331: TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,X,Xdot)=0
333: Collective on TS and Vec
335: Input Parameters:
336: + ts - the TS context
337: . t - current time
338: . X - state vector
339: . Xdot - time derivative of state vector
340: - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate
342: Output Parameter:
343: . Y - right hand side
345: Note:
346: Most users should not need to explicitly call this routine, as it
347: is used internally within the nonlinear solvers.
349: If the user did did not write their equations in implicit form, this
350: function recasts them in implicit form.
352: Level: developer
354: .keywords: TS, compute
356: .seealso: TSSetIFunction(), TSComputeRHSFunction()
357: @*/
358: PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec Y,PetscBool imex)
359: {
368: if (!ts->userops->rhsfunction && !ts->userops->ifunction) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
370: PetscLogEventBegin(TS_FunctionEval,ts,X,Xdot,Y);
371: if (ts->userops->ifunction) {
372: PetscStackPush("TS user implicit function");
373: (*ts->userops->ifunction)(ts,t,X,Xdot,Y,ts->funP);
374: PetscStackPop;
375: }
376: if (imex) {
377: if (!ts->userops->ifunction) {
378: VecCopy(Xdot,Y);
379: }
380: } else if (ts->userops->rhsfunction) {
381: if (ts->userops->ifunction) {
382: Vec Frhs;
383: TSGetRHSVec_Private(ts,&Frhs);
384: TSComputeRHSFunction(ts,t,X,Frhs);
385: VecAXPY(Y,-1,Frhs);
386: } else {
387: TSComputeRHSFunction(ts,t,X,Y);
388: VecAYPX(Y,-1,Xdot);
389: }
390: }
391: PetscLogEventEnd(TS_FunctionEval,ts,X,Xdot,Y);
392: return(0);
393: }
397: /*@
398: TSComputeIJacobian - Evaluates the Jacobian of the DAE
400: Collective on TS and Vec
402: Input
403: Input Parameters:
404: + ts - the TS context
405: . t - current timestep
406: . X - state vector
407: . Xdot - time derivative of state vector
408: . shift - shift to apply, see note below
409: - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate
411: Output Parameters:
412: + A - Jacobian matrix
413: . B - optional preconditioning matrix
414: - flag - flag indicating matrix structure
416: Notes:
417: If F(t,X,Xdot)=0 is the DAE, the required Jacobian is
419: dF/dX + shift*dF/dXdot
421: Most users should not need to explicitly call this routine, as it
422: is used internally within the nonlinear solvers.
424: Level: developer
426: .keywords: TS, compute, Jacobian, matrix
428: .seealso: TSSetIJacobian()
429: @*/
430: PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,PetscBool imex)
431: {
432: PetscInt Xstate, Xdotstate;
444: PetscObjectStateQuery((PetscObject)X,&Xstate);
445: PetscObjectStateQuery((PetscObject)Xdot,&Xdotstate);
446: if (ts->ijacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->ijacobian.X == X && ts->ijacobian.Xstate == Xstate && ts->ijacobian.Xdot == Xdot && ts->ijacobian.Xdotstate == Xdotstate && ts->ijacobian.imex == imex))) {
447: *flg = ts->ijacobian.mstructure;
448: MatScale(*A, shift / ts->ijacobian.shift);
449: return(0);
450: }
452: if (!ts->userops->rhsjacobian && !ts->userops->ijacobian) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
454: *flg = SAME_NONZERO_PATTERN; /* In case we're solving a linear problem in which case it wouldn't get initialized below. */
455: PetscLogEventBegin(TS_JacobianEval,ts,X,*A,*B);
456: if (ts->userops->ijacobian) {
457: *flg = DIFFERENT_NONZERO_PATTERN;
458: PetscStackPush("TS user implicit Jacobian");
459: (*ts->userops->ijacobian)(ts,t,X,Xdot,shift,A,B,flg,ts->jacP);
460: PetscStackPop;
461: /* make sure user returned a correct Jacobian and preconditioner */
464: }
465: if (imex) {
466: if (!ts->userops->ijacobian) { /* system was written as Xdot = F(t,X) */
467: MatZeroEntries(*A);
468: MatShift(*A,shift);
469: if (*A != *B) {
470: MatZeroEntries(*B);
471: MatShift(*B,shift);
472: }
473: *flg = SAME_PRECONDITIONER;
474: }
475: } else {
476: if (!ts->userops->ijacobian) {
477: TSComputeRHSJacobian(ts,t,X,A,B,flg);
478: MatScale(*A,-1);
479: MatShift(*A,shift);
480: if (*A != *B) {
481: MatScale(*B,-1);
482: MatShift(*B,shift);
483: }
484: } else if (ts->userops->rhsjacobian) {
485: Mat Arhs,Brhs;
486: MatStructure axpy,flg2 = DIFFERENT_NONZERO_PATTERN;
487: TSGetRHSMats_Private(ts,&Arhs,&Brhs);
488: TSComputeRHSJacobian(ts,t,X,&Arhs,&Brhs,&flg2);
489: axpy = (*flg == flg2) ? SAME_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;
490: MatAXPY(*A,-1,Arhs,axpy);
491: if (*A != *B) {
492: MatAXPY(*B,-1,Brhs,axpy);
493: }
494: *flg = PetscMin(*flg,flg2);
495: }
496: }
498: ts->ijacobian.time = t;
499: ts->ijacobian.X = X;
500: ts->ijacobian.Xdot = Xdot;
501: PetscObjectStateQuery((PetscObject)X,&ts->ijacobian.Xstate);
502: PetscObjectStateQuery((PetscObject)Xdot,&ts->ijacobian.Xdotstate);
503: ts->ijacobian.shift = shift;
504: ts->ijacobian.imex = imex;
505: ts->ijacobian.mstructure = *flg;
506: PetscLogEventEnd(TS_JacobianEval,ts,X,*A,*B);
507: return(0);
508: }
512: /*@C
513: TSSetRHSFunction - Sets the routine for evaluating the function,
514: F(t,u), where U_t = F(t,u).
516: Logically Collective on TS
518: Input Parameters:
519: + ts - the TS context obtained from TSCreate()
520: . r - vector to put the computed right hand side (or PETSC_NULL to have it created)
521: . f - routine for evaluating the right-hand-side function
522: - ctx - [optional] user-defined context for private data for the
523: function evaluation routine (may be PETSC_NULL)
525: Calling sequence of func:
526: $ func (TS ts,PetscReal t,Vec u,Vec F,void *ctx);
528: + t - current timestep
529: . u - input vector
530: . F - function vector
531: - ctx - [optional] user-defined function context
533: Level: beginner
535: .keywords: TS, timestep, set, right-hand-side, function
537: .seealso: TSSetRHSJacobian(), TSSetIJacobian()
538: @*/
539: PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx)
540: {
542: SNES snes;
543: Vec ralloc = PETSC_NULL;
548: if (f) ts->userops->rhsfunction = f;
549: if (ctx) ts->funP = ctx;
550: TSGetSNES(ts,&snes);
551: if (!r && !ts->dm && ts->vec_sol) {
552: VecDuplicate(ts->vec_sol,&ralloc);
553: r = ralloc;
554: }
555: SNESSetFunction(snes,r,SNESTSFormFunction,ts);
556: VecDestroy(&ralloc);
557: return(0);
558: }
562: /*@C
563: TSSetRHSJacobian - Sets the function to compute the Jacobian of F,
564: where U_t = F(U,t), as well as the location to store the matrix.
566: Logically Collective on TS
568: Input Parameters:
569: + ts - the TS context obtained from TSCreate()
570: . A - Jacobian matrix
571: . B - preconditioner matrix (usually same as A)
572: . f - the Jacobian evaluation routine
573: - ctx - [optional] user-defined context for private data for the
574: Jacobian evaluation routine (may be PETSC_NULL)
576: Calling sequence of func:
577: $ func (TS ts,PetscReal t,Vec u,Mat *A,Mat *B,MatStructure *flag,void *ctx);
579: + t - current timestep
580: . u - input vector
581: . A - matrix A, where U_t = A(t)u
582: . B - preconditioner matrix, usually the same as A
583: . flag - flag indicating information about the preconditioner matrix
584: structure (same as flag in KSPSetOperators())
585: - ctx - [optional] user-defined context for matrix evaluation routine
587: Notes:
588: See KSPSetOperators() for important information about setting the flag
589: output parameter in the routine func(). Be sure to read this information!
591: The routine func() takes Mat * as the matrix arguments rather than Mat.
592: This allows the matrix evaluation routine to replace A and/or B with a
593: completely new matrix structure (not just different matrix elements)
594: when appropriate, for instance, if the nonzero structure is changing
595: throughout the global iterations.
597: Level: beginner
598:
599: .keywords: TS, timestep, set, right-hand-side, Jacobian
601: .seealso: SNESDefaultComputeJacobianColor(), TSSetRHSFunction()
603: @*/
604: PetscErrorCode TSSetRHSJacobian(TS ts,Mat A,Mat B,TSRHSJacobian f,void *ctx)
605: {
607: SNES snes;
616: if (f) ts->userops->rhsjacobian = f;
617: if (ctx) ts->jacP = ctx;
618: TSGetSNES(ts,&snes);
619: if (!ts->userops->ijacobian) {
620: SNESSetJacobian(snes,A,B,SNESTSFormJacobian,ts);
621: }
622: if (A) {
623: PetscObjectReference((PetscObject)A);
624: MatDestroy(&ts->Arhs);
625: ts->Arhs = A;
626: }
627: if (B) {
628: PetscObjectReference((PetscObject)B);
629: MatDestroy(&ts->Brhs);
630: ts->Brhs = B;
631: }
632: return(0);
633: }
638: /*@C
639: TSSetIFunction - Set the function to compute F(t,U,U_t) where F = 0 is the DAE to be solved.
641: Logically Collective on TS
643: Input Parameters:
644: + ts - the TS context obtained from TSCreate()
645: . r - vector to hold the residual (or PETSC_NULL to have it created internally)
646: . f - the function evaluation routine
647: - ctx - user-defined context for private data for the function evaluation routine (may be PETSC_NULL)
649: Calling sequence of f:
650: $ f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx);
652: + t - time at step/stage being solved
653: . u - state vector
654: . u_t - time derivative of state vector
655: . F - function vector
656: - ctx - [optional] user-defined context for matrix evaluation routine
658: Important:
659: The user MUST call either this routine, TSSetRHSFunction(). This routine must be used when not solving an ODE, for example a DAE.
661: Level: beginner
663: .keywords: TS, timestep, set, DAE, Jacobian
665: .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian()
666: @*/
667: PetscErrorCode TSSetIFunction(TS ts,Vec res,TSIFunction f,void *ctx)
668: {
670: SNES snes;
671: Vec resalloc = PETSC_NULL;
676: if (f) ts->userops->ifunction = f;
677: if (ctx) ts->funP = ctx;
678: TSGetSNES(ts,&snes);
679: if (!res && !ts->dm && ts->vec_sol) {
680: VecDuplicate(ts->vec_sol,&resalloc);
681: res = resalloc;
682: }
683: SNESSetFunction(snes,res,SNESTSFormFunction,ts);
684: VecDestroy(&resalloc);
685: return(0);
686: }
690: /*@C
691: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it.
693: Not Collective
695: Input Parameter:
696: . ts - the TS context
698: Output Parameter:
699: + r - vector to hold residual (or PETSC_NULL)
700: . func - the function to compute residual (or PETSC_NULL)
701: - ctx - the function context (or PETSC_NULL)
703: Level: advanced
705: .keywords: TS, nonlinear, get, function
707: .seealso: TSSetIFunction(), SNESGetFunction()
708: @*/
709: PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx)
710: {
712: SNES snes;
716: TSGetSNES(ts,&snes);
717: SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);
718: if (func) *func = ts->userops->ifunction;
719: if (ctx) *ctx = ts->funP;
720: return(0);
721: }
725: /*@C
726: TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it.
728: Not Collective
730: Input Parameter:
731: . ts - the TS context
733: Output Parameter:
734: + r - vector to hold computed right hand side (or PETSC_NULL)
735: . func - the function to compute right hand side (or PETSC_NULL)
736: - ctx - the function context (or PETSC_NULL)
738: Level: advanced
740: .keywords: TS, nonlinear, get, function
742: .seealso: TSSetRhsfunction(), SNESGetFunction()
743: @*/
744: PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx)
745: {
747: SNES snes;
751: TSGetSNES(ts,&snes);
752: SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);
753: if (func) *func = ts->userops->rhsfunction;
754: if (ctx) *ctx = ts->funP;
755: return(0);
756: }
760: /*@C
761: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
762: you provided with TSSetIFunction().
764: Logically Collective on TS
766: Input Parameters:
767: + ts - the TS context obtained from TSCreate()
768: . A - Jacobian matrix
769: . B - preconditioning matrix for A (may be same as A)
770: . f - the Jacobian evaluation routine
771: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be PETSC_NULL)
773: Calling sequence of f:
774: $ f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx);
776: + t - time at step/stage being solved
777: . U - state vector
778: . U_t - time derivative of state vector
779: . a - shift
780: . A - Jacobian of G(U) = F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
781: . B - preconditioning matrix for A, may be same as A
782: . flag - flag indicating information about the preconditioner matrix
783: structure (same as flag in KSPSetOperators())
784: - ctx - [optional] user-defined context for matrix evaluation routine
786: Notes:
787: The matrices A and B are exactly the matrices that are used by SNES for the nonlinear solve.
789: The matrix dF/dU + a*dF/dU_t you provide turns out to be
790: the Jacobian of G(U) = F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
791: The time integrator internally approximates U_t by W+a*U where the positive "shift"
792: a and vector W depend on the integration method, step size, and past states. For example with
793: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
794: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
796: Level: beginner
798: .keywords: TS, timestep, DAE, Jacobian
800: .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESDefaultComputeJacobianColor(), SNESDefaultComputeJacobian()
802: @*/
803: PetscErrorCode TSSetIJacobian(TS ts,Mat A,Mat B,TSIJacobian f,void *ctx)
804: {
806: SNES snes;
814: if (f) ts->userops->ijacobian = f;
815: if (ctx) ts->jacP = ctx;
816: TSGetSNES(ts,&snes);
817: SNESSetJacobian(snes,A,B,SNESTSFormJacobian,ts);
818: return(0);
819: }
823: /*@C
824: TSView - Prints the TS data structure.
826: Collective on TS
828: Input Parameters:
829: + ts - the TS context obtained from TSCreate()
830: - viewer - visualization context
832: Options Database Key:
833: . -ts_view - calls TSView() at end of TSStep()
835: Notes:
836: The available visualization contexts include
837: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
838: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
839: output where only the first processor opens
840: the file. All other processors send their
841: data to the first processor to print.
843: The user can open an alternative visualization context with
844: PetscViewerASCIIOpen() - output to a specified file.
846: Level: beginner
848: .keywords: TS, timestep, view
850: .seealso: PetscViewerASCIIOpen()
851: @*/
852: PetscErrorCode TSView(TS ts,PetscViewer viewer)
853: {
855: const TSType type;
856: PetscBool iascii,isstring,isundials;
860: if (!viewer) {
861: PetscViewerASCIIGetStdout(((PetscObject)ts)->comm,&viewer);
862: }
866: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
867: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
868: if (iascii) {
869: PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer,"TS Object");
870: PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);
871: PetscViewerASCIIPrintf(viewer," maximum time=%G\n",ts->max_time);
872: if (ts->problem_type == TS_NONLINEAR) {
873: PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->nonlinear_its);
874: PetscViewerASCIIPrintf(viewer," total number of nonlinear solve failures=%D\n",ts->num_snes_failures);
875: }
876: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->linear_its);
877: PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);
878: if (ts->ops->view) {
879: PetscViewerASCIIPushTab(viewer);
880: (*ts->ops->view)(ts,viewer);
881: PetscViewerASCIIPopTab(viewer);
882: }
883: } else if (isstring) {
884: TSGetType(ts,&type);
885: PetscViewerStringSPrintf(viewer," %-7.7s",type);
886: }
887: PetscViewerASCIIPushTab(viewer);
888: PetscTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);
889: PetscViewerASCIIPopTab(viewer);
890: return(0);
891: }
896: /*@
897: TSSetApplicationContext - Sets an optional user-defined context for
898: the timesteppers.
900: Logically Collective on TS
902: Input Parameters:
903: + ts - the TS context obtained from TSCreate()
904: - usrP - optional user context
906: Level: intermediate
908: .keywords: TS, timestep, set, application, context
910: .seealso: TSGetApplicationContext()
911: @*/
912: PetscErrorCode TSSetApplicationContext(TS ts,void *usrP)
913: {
916: ts->user = usrP;
917: return(0);
918: }
922: /*@
923: TSGetApplicationContext - Gets the user-defined context for the
924: timestepper.
926: Not Collective
928: Input Parameter:
929: . ts - the TS context obtained from TSCreate()
931: Output Parameter:
932: . usrP - user context
934: Level: intermediate
936: .keywords: TS, timestep, get, application, context
938: .seealso: TSSetApplicationContext()
939: @*/
940: PetscErrorCode TSGetApplicationContext(TS ts,void *usrP)
941: {
944: *(void**)usrP = ts->user;
945: return(0);
946: }
950: /*@
951: TSGetTimeStepNumber - Gets the current number of timesteps.
953: Not Collective
955: Input Parameter:
956: . ts - the TS context obtained from TSCreate()
958: Output Parameter:
959: . iter - number steps so far
961: Level: intermediate
963: .keywords: TS, timestep, get, iteration, number
964: @*/
965: PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt* iter)
966: {
970: *iter = ts->steps;
971: return(0);
972: }
976: /*@
977: TSSetInitialTimeStep - Sets the initial timestep to be used,
978: as well as the initial time.
980: Logically Collective on TS
982: Input Parameters:
983: + ts - the TS context obtained from TSCreate()
984: . initial_time - the initial time
985: - time_step - the size of the timestep
987: Level: intermediate
989: .seealso: TSSetTimeStep(), TSGetTimeStep()
991: .keywords: TS, set, initial, timestep
992: @*/
993: PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step)
994: {
999: TSSetTimeStep(ts,time_step);
1000: TSSetTime(ts,initial_time);
1001: return(0);
1002: }
1006: /*@
1007: TSSetTimeStep - Allows one to reset the timestep at any time,
1008: useful for simple pseudo-timestepping codes.
1010: Logically Collective on TS
1012: Input Parameters:
1013: + ts - the TS context obtained from TSCreate()
1014: - time_step - the size of the timestep
1016: Level: intermediate
1018: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()
1020: .keywords: TS, set, timestep
1021: @*/
1022: PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step)
1023: {
1027: ts->time_step = time_step;
1028: return(0);
1029: }
1033: /*@
1034: TSSetExactFinalTime - Determines whether to interpolate solution to the
1035: exact final time requested by the user or just returns it at the final time
1036: it computed.
1038: Logically Collective on TS
1040: Input Parameter:
1041: + ts - the time-step context
1042: - ft - PETSC_TRUE if interpolates, else PETSC_FALSE
1044: Level: beginner
1046: .seealso: TSSetDuration()
1047: @*/
1048: PetscErrorCode TSSetExactFinalTime(TS ts,PetscBool flg)
1049: {
1054: ts->exact_final_time = flg;
1055: return(0);
1056: }
1060: /*@
1061: TSGetTimeStep - Gets the current timestep size.
1063: Not Collective
1065: Input Parameter:
1066: . ts - the TS context obtained from TSCreate()
1068: Output Parameter:
1069: . dt - the current timestep size
1071: Level: intermediate
1073: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()
1075: .keywords: TS, get, timestep
1076: @*/
1077: PetscErrorCode TSGetTimeStep(TS ts,PetscReal* dt)
1078: {
1082: *dt = ts->time_step;
1083: return(0);
1084: }
1088: /*@
1089: TSGetSolution - Returns the solution at the present timestep. It
1090: is valid to call this routine inside the function that you are evaluating
1091: in order to move to the new timestep. This vector not changed until
1092: the solution at the next timestep has been calculated.
1094: Not Collective, but Vec returned is parallel if TS is parallel
1096: Input Parameter:
1097: . ts - the TS context obtained from TSCreate()
1099: Output Parameter:
1100: . v - the vector containing the solution
1102: Level: intermediate
1104: .seealso: TSGetTimeStep()
1106: .keywords: TS, timestep, get, solution
1107: @*/
1108: PetscErrorCode TSGetSolution(TS ts,Vec *v)
1109: {
1113: *v = ts->vec_sol;
1114: return(0);
1115: }
1117: /* ----- Routines to initialize and destroy a timestepper ---- */
1120: /*@
1121: TSSetProblemType - Sets the type of problem to be solved.
1123: Not collective
1125: Input Parameters:
1126: + ts - The TS
1127: - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
1128: .vb
1129: U_t = A U
1130: U_t = A(t) U
1131: U_t = F(t,U)
1132: .ve
1134: Level: beginner
1136: .keywords: TS, problem type
1137: .seealso: TSSetUp(), TSProblemType, TS
1138: @*/
1139: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
1140: {
1145: ts->problem_type = type;
1146: if (type == TS_LINEAR) {
1147: SNES snes;
1148: TSGetSNES(ts,&snes);
1149: SNESSetType(snes,SNESKSPONLY);
1150: }
1151: return(0);
1152: }
1156: /*@C
1157: TSGetProblemType - Gets the type of problem to be solved.
1159: Not collective
1161: Input Parameter:
1162: . ts - The TS
1164: Output Parameter:
1165: . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
1166: .vb
1167: M U_t = A U
1168: M(t) U_t = A(t) U
1169: U_t = F(t,U)
1170: .ve
1172: Level: beginner
1174: .keywords: TS, problem type
1175: .seealso: TSSetUp(), TSProblemType, TS
1176: @*/
1177: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
1178: {
1182: *type = ts->problem_type;
1183: return(0);
1184: }
1188: /*@
1189: TSSetUp - Sets up the internal data structures for the later use
1190: of a timestepper.
1192: Collective on TS
1194: Input Parameter:
1195: . ts - the TS context obtained from TSCreate()
1197: Notes:
1198: For basic use of the TS solvers the user need not explicitly call
1199: TSSetUp(), since these actions will automatically occur during
1200: the call to TSStep(). However, if one wishes to control this
1201: phase separately, TSSetUp() should be called after TSCreate()
1202: and optional routines of the form TSSetXXX(), but before TSStep().
1204: Level: advanced
1206: .keywords: TS, timestep, setup
1208: .seealso: TSCreate(), TSStep(), TSDestroy()
1209: @*/
1210: PetscErrorCode TSSetUp(TS ts)
1211: {
1216: if (ts->setupcalled) return(0);
1218: if (!((PetscObject)ts)->type_name) {
1219: TSSetType(ts,TSEULER);
1220: }
1221: if (ts->exact_final_time == PETSC_DECIDE) ts->exact_final_time = PETSC_FALSE;
1223: if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first");
1225: TSGetAdapt(ts,&ts->adapt);
1227: if (ts->ops->setup) {
1228: (*ts->ops->setup)(ts);
1229: }
1231: ts->setupcalled = PETSC_TRUE;
1232: return(0);
1233: }
1237: /*@
1238: TSReset - Resets a TS context and removes any allocated Vecs and Mats.
1240: Collective on TS
1242: Input Parameter:
1243: . ts - the TS context obtained from TSCreate()
1245: Level: beginner
1247: .keywords: TS, timestep, reset
1249: .seealso: TSCreate(), TSSetup(), TSDestroy()
1250: @*/
1251: PetscErrorCode TSReset(TS ts)
1252: {
1257: if (ts->ops->reset) {
1258: (*ts->ops->reset)(ts);
1259: }
1260: if (ts->snes) {SNESReset(ts->snes);}
1261: MatDestroy(&ts->Arhs);
1262: MatDestroy(&ts->Brhs);
1263: VecDestroy(&ts->Frhs);
1264: VecDestroy(&ts->vec_sol);
1265: VecDestroyVecs(ts->nwork,&ts->work);
1266: ts->setupcalled = PETSC_FALSE;
1267: return(0);
1268: }
1272: /*@
1273: TSDestroy - Destroys the timestepper context that was created
1274: with TSCreate().
1276: Collective on TS
1278: Input Parameter:
1279: . ts - the TS context obtained from TSCreate()
1281: Level: beginner
1283: .keywords: TS, timestepper, destroy
1285: .seealso: TSCreate(), TSSetUp(), TSSolve()
1286: @*/
1287: PetscErrorCode TSDestroy(TS *ts)
1288: {
1292: if (!*ts) return(0);
1294: if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; return(0);}
1296: TSReset((*ts));
1298: /* if memory was published with AMS then destroy it */
1299: PetscObjectDepublish((*ts));
1300: if ((*ts)->ops->destroy) {(*(*ts)->ops->destroy)((*ts));}
1302: TSAdaptDestroy(&(*ts)->adapt);
1303: SNESDestroy(&(*ts)->snes);
1304: DMDestroy(&(*ts)->dm);
1305: TSMonitorCancel((*ts));
1307: PetscFree((*ts)->userops);
1309: PetscHeaderDestroy(ts);
1310: return(0);
1311: }
1315: /*@
1316: TSGetSNES - Returns the SNES (nonlinear solver) associated with
1317: a TS (timestepper) context. Valid only for nonlinear problems.
1319: Not Collective, but SNES is parallel if TS is parallel
1321: Input Parameter:
1322: . ts - the TS context obtained from TSCreate()
1324: Output Parameter:
1325: . snes - the nonlinear solver context
1327: Notes:
1328: The user can then directly manipulate the SNES context to set various
1329: options, etc. Likewise, the user can then extract and manipulate the
1330: KSP, KSP, and PC contexts as well.
1332: TSGetSNES() does not work for integrators that do not use SNES; in
1333: this case TSGetSNES() returns PETSC_NULL in snes.
1335: Level: beginner
1337: .keywords: timestep, get, SNES
1338: @*/
1339: PetscErrorCode TSGetSNES(TS ts,SNES *snes)
1340: {
1346: if (!ts->snes) {
1347: SNESCreate(((PetscObject)ts)->comm,&ts->snes);
1348: PetscLogObjectParent(ts,ts->snes);
1349: PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);
1350: if (ts->problem_type == TS_LINEAR) {
1351: SNESSetType(ts->snes,SNESKSPONLY);
1352: }
1353: }
1354: *snes = ts->snes;
1355: return(0);
1356: }
1360: /*@
1361: TSGetKSP - Returns the KSP (linear solver) associated with
1362: a TS (timestepper) context.
1364: Not Collective, but KSP is parallel if TS is parallel
1366: Input Parameter:
1367: . ts - the TS context obtained from TSCreate()
1369: Output Parameter:
1370: . ksp - the nonlinear solver context
1372: Notes:
1373: The user can then directly manipulate the KSP context to set various
1374: options, etc. Likewise, the user can then extract and manipulate the
1375: KSP and PC contexts as well.
1377: TSGetKSP() does not work for integrators that do not use KSP;
1378: in this case TSGetKSP() returns PETSC_NULL in ksp.
1380: Level: beginner
1382: .keywords: timestep, get, KSP
1383: @*/
1384: PetscErrorCode TSGetKSP(TS ts,KSP *ksp)
1385: {
1387: SNES snes;
1392: if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first");
1393: if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()");
1394: TSGetSNES(ts,&snes);
1395: SNESGetKSP(snes,ksp);
1396: return(0);
1397: }
1399: /* ----------- Routines to set solver parameters ---------- */
1403: /*@
1404: TSGetDuration - Gets the maximum number of timesteps to use and
1405: maximum time for iteration.
1407: Not Collective
1409: Input Parameters:
1410: + ts - the TS context obtained from TSCreate()
1411: . maxsteps - maximum number of iterations to use, or PETSC_NULL
1412: - maxtime - final time to iterate to, or PETSC_NULL
1414: Level: intermediate
1416: .keywords: TS, timestep, get, maximum, iterations, time
1417: @*/
1418: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
1419: {
1422: if (maxsteps) {
1424: *maxsteps = ts->max_steps;
1425: }
1426: if (maxtime) {
1428: *maxtime = ts->max_time;
1429: }
1430: return(0);
1431: }
1435: /*@
1436: TSSetDuration - Sets the maximum number of timesteps to use and
1437: maximum time for iteration.
1439: Logically Collective on TS
1441: Input Parameters:
1442: + ts - the TS context obtained from TSCreate()
1443: . maxsteps - maximum number of iterations to use
1444: - maxtime - final time to iterate to
1446: Options Database Keys:
1447: . -ts_max_steps <maxsteps> - Sets maxsteps
1448: . -ts_final_time <maxtime> - Sets maxtime
1450: Notes:
1451: The default maximum number of iterations is 5000. Default time is 5.0
1453: Level: intermediate
1455: .keywords: TS, timestep, set, maximum, iterations
1457: .seealso: TSSetExactFinalTime()
1458: @*/
1459: PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime)
1460: {
1465: if (maxsteps >= 0) ts->max_steps = maxsteps;
1466: if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
1467: return(0);
1468: }
1472: /*@
1473: TSSetSolution - Sets the initial solution vector
1474: for use by the TS routines.
1476: Logically Collective on TS and Vec
1478: Input Parameters:
1479: + ts - the TS context obtained from TSCreate()
1480: - x - the solution vector
1482: Level: beginner
1484: .keywords: TS, timestep, set, solution, initial conditions
1485: @*/
1486: PetscErrorCode TSSetSolution(TS ts,Vec x)
1487: {
1493: PetscObjectReference((PetscObject)x);
1494: VecDestroy(&ts->vec_sol);
1495: ts->vec_sol = x;
1496: return(0);
1497: }
1501: /*@C
1502: TSSetPreStep - Sets the general-purpose function
1503: called once at the beginning of each time step.
1505: Logically Collective on TS
1507: Input Parameters:
1508: + ts - The TS context obtained from TSCreate()
1509: - func - The function
1511: Calling sequence of func:
1512: . func (TS ts);
1514: Level: intermediate
1516: .keywords: TS, timestep
1517: @*/
1518: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
1519: {
1522: ts->ops->prestep = func;
1523: return(0);
1524: }
1528: /*@
1529: TSPreStep - Runs the user-defined pre-step function.
1531: Collective on TS
1533: Input Parameters:
1534: . ts - The TS context obtained from TSCreate()
1536: Notes:
1537: TSPreStep() is typically used within time stepping implementations,
1538: so most users would not generally call this routine themselves.
1540: Level: developer
1542: .keywords: TS, timestep
1543: @*/
1544: PetscErrorCode TSPreStep(TS ts)
1545: {
1550: if (ts->ops->prestep) {
1551: PetscStackPush("TS PreStep function");
1552: (*ts->ops->prestep)(ts);
1553: PetscStackPop;
1554: }
1555: return(0);
1556: }
1560: /*@C
1561: TSSetPostStep - Sets the general-purpose function
1562: called once at the end of each time step.
1564: Logically Collective on TS
1566: Input Parameters:
1567: + ts - The TS context obtained from TSCreate()
1568: - func - The function
1570: Calling sequence of func:
1571: . func (TS ts);
1573: Level: intermediate
1575: .keywords: TS, timestep
1576: @*/
1577: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
1578: {
1581: ts->ops->poststep = func;
1582: return(0);
1583: }
1587: /*@
1588: TSPostStep - Runs the user-defined post-step function.
1590: Collective on TS
1592: Input Parameters:
1593: . ts - The TS context obtained from TSCreate()
1595: Notes:
1596: TSPostStep() is typically used within time stepping implementations,
1597: so most users would not generally call this routine themselves.
1599: Level: developer
1601: .keywords: TS, timestep
1602: @*/
1603: PetscErrorCode TSPostStep(TS ts)
1604: {
1609: if (ts->ops->poststep) {
1610: PetscStackPush("TS PostStep function");
1611: (*ts->ops->poststep)(ts);
1612: PetscStackPop;
1613: }
1614: return(0);
1615: }
1617: /* ------------ Routines to set performance monitoring options ----------- */
1621: /*@C
1622: TSMonitorSet - Sets an ADDITIONAL function that is to be used at every
1623: timestep to display the iteration's progress.
1625: Logically Collective on TS
1627: Input Parameters:
1628: + ts - the TS context obtained from TSCreate()
1629: . monitor - monitoring routine
1630: . mctx - [optional] user-defined context for private data for the
1631: monitor routine (use PETSC_NULL if no context is desired)
1632: - monitordestroy - [optional] routine that frees monitor context
1633: (may be PETSC_NULL)
1635: Calling sequence of monitor:
1636: $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec x,void *mctx)
1638: + ts - the TS context
1639: . steps - iteration number
1640: . time - current time
1641: . x - current iterate
1642: - mctx - [optional] monitoring context
1644: Notes:
1645: This routine adds an additional monitor to the list of monitors that
1646: already has been loaded.
1648: Fortran notes: Only a single monitor function can be set for each TS object
1650: Level: intermediate
1652: .keywords: TS, timestep, set, monitor
1654: .seealso: TSMonitorDefault(), TSMonitorCancel()
1655: @*/
1656: PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**))
1657: {
1660: if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
1661: ts->monitor[ts->numbermonitors] = monitor;
1662: ts->mdestroy[ts->numbermonitors] = mdestroy;
1663: ts->monitorcontext[ts->numbermonitors++] = (void*)mctx;
1664: return(0);
1665: }
1669: /*@C
1670: TSMonitorCancel - Clears all the monitors that have been set on a time-step object.
1672: Logically Collective on TS
1674: Input Parameters:
1675: . ts - the TS context obtained from TSCreate()
1677: Notes:
1678: There is no way to remove a single, specific monitor.
1680: Level: intermediate
1682: .keywords: TS, timestep, set, monitor
1684: .seealso: TSMonitorDefault(), TSMonitorSet()
1685: @*/
1686: PetscErrorCode TSMonitorCancel(TS ts)
1687: {
1689: PetscInt i;
1693: for (i=0; i<ts->numbermonitors; i++) {
1694: if (ts->mdestroy[i]) {
1695: (*ts->mdestroy[i])(&ts->monitorcontext[i]);
1696: }
1697: }
1698: ts->numbermonitors = 0;
1699: return(0);
1700: }
1704: /*@
1705: TSMonitorDefault - Sets the Default monitor
1707: Level: intermediate
1709: .keywords: TS, set, monitor
1711: .seealso: TSMonitorDefault(), TSMonitorSet()
1712: @*/
1713: PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
1714: {
1716: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ts)->comm);
1719: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
1720: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);
1721: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
1722: return(0);
1723: }
1727: /*@
1728: TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available.
1730: Logically Collective on TS
1732: Input Argument:
1733: . ts - time stepping context
1735: Output Argument:
1736: . flg - PETSC_TRUE or PETSC_FALSE
1738: Level: intermediate
1740: .keywords: TS, set
1742: .seealso: TSInterpolate(), TSSetPostStep()
1743: @*/
1744: PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg)
1745: {
1749: ts->retain_stages = flg;
1750: return(0);
1751: }
1755: /*@
1756: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
1758: Collective on TS
1760: Input Argument:
1761: + ts - time stepping context
1762: - t - time to interpolate to
1764: Output Argument:
1765: . X - state at given time
1767: Notes:
1768: The user should call TSSetRetainStages() before taking a step in which interpolation will be requested.
1770: Level: intermediate
1772: Developer Notes:
1773: TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
1775: .keywords: TS, set
1777: .seealso: TSSetRetainStages(), TSSetPostStep()
1778: @*/
1779: PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec X)
1780: {
1785: if (t < ts->ptime - ts->time_step_prev || t > ts->ptime) SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Requested time %G not in last time steps [%G,%G]",t,ts->ptime-ts->time_step_prev,ts->ptime);
1786: if (!ts->ops->interpolate) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name);
1787: (*ts->ops->interpolate)(ts,t,X);
1788: return(0);
1789: }
1793: /*@
1794: TSStep - Steps one time step
1796: Collective on TS
1798: Input Parameter:
1799: . ts - the TS context obtained from TSCreate()
1801: Level: intermediate
1803: .keywords: TS, timestep, solve
1805: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve()
1806: @*/
1807: PetscErrorCode TSStep(TS ts)
1808: {
1809: PetscReal ptime_prev;
1814: TSSetUp(ts);
1816: ts->reason = TS_CONVERGED_ITERATING;
1818: ptime_prev = ts->ptime;
1819: PetscLogEventBegin(TS_Step,ts,0,0,0);
1820: (*ts->ops->step)(ts);
1821: PetscLogEventEnd(TS_Step,ts,0,0,0);
1822: ts->time_step_prev = ts->ptime - ptime_prev;
1824: if (ts->reason < 0) {
1825: if (ts->errorifstepfailed) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
1826: } else if (!ts->reason) {
1827: if (ts->steps >= ts->max_steps)
1828: ts->reason = TS_CONVERGED_ITS;
1829: else if (ts->ptime >= ts->max_time)
1830: ts->reason = TS_CONVERGED_TIME;
1831: }
1833: return(0);
1834: }
1838: /*@
1839: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
1841: Collective on TS
1843: Input Arguments:
1844: + ts - time stepping context
1845: . order - desired order of accuracy
1846: - done - whether the step was evaluated at this order (pass PETSC_NULL to generate an error if not available)
1848: Output Arguments:
1849: . X - state at the end of the current step
1851: Level: advanced
1853: Notes:
1854: This function cannot be called until all stages have been evaluated.
1855: 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.
1857: .seealso: TSStep(), TSAdapt
1858: @*/
1859: PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec X,PetscBool *done)
1860: {
1867: if (!ts->ops->evaluatestep) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name);
1868: (*ts->ops->evaluatestep)(ts,order,X,done);
1869: return(0);
1870: }
1874: /*@
1875: TSSolve - Steps the requested number of timesteps.
1877: Collective on TS
1879: Input Parameter:
1880: + ts - the TS context obtained from TSCreate()
1881: - x - the solution vector
1883: Output Parameter:
1884: . ftime - time of the state vector x upon completion
1886: Level: beginner
1888: Notes:
1889: The final time returned by this function may be different from the time of the internally
1890: held state accessible by TSGetSolution() and TSGetTime() because the method may have
1891: stepped over the final time.
1893: .keywords: TS, timestep, solve
1895: .seealso: TSCreate(), TSSetSolution(), TSStep()
1896: @*/
1897: PetscErrorCode TSSolve(TS ts,Vec x,PetscReal *ftime)
1898: {
1899: PetscBool flg;
1900: char filename[PETSC_MAX_PATH_LEN];
1901: PetscViewer viewer;
1907: if (ts->exact_final_time) { /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */
1908: if (!ts->vec_sol || x == ts->vec_sol) {
1909: Vec y;
1910: VecDuplicate(x,&y);
1911: TSSetSolution(ts,y);
1912: VecDestroy(&y); /* grant ownership */
1913: }
1914: VecCopy(x,ts->vec_sol);
1915: } else {
1916: TSSetSolution(ts,x);
1917: }
1918: TSSetUp(ts);
1919: /* reset time step and iteration counters */
1920: ts->steps = 0;
1921: ts->linear_its = 0;
1922: ts->nonlinear_its = 0;
1923: ts->num_snes_failures = 0;
1924: ts->reject = 0;
1925: ts->reason = TS_CONVERGED_ITERATING;
1927: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
1928: (*ts->ops->solve)(ts);
1929: VecCopy(ts->vec_sol,x);
1930: if (ftime) *ftime = ts->ptime;
1931: } else {
1932: /* steps the requested number of timesteps. */
1933: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
1934: if (ts->steps >= ts->max_steps)
1935: ts->reason = TS_CONVERGED_ITS;
1936: else if (ts->ptime >= ts->max_time)
1937: ts->reason = TS_CONVERGED_TIME;
1938: while (!ts->reason) {
1939: TSPreStep(ts);
1940: TSStep(ts);
1941: TSPostStep(ts);
1942: TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
1943: }
1944: if (ts->exact_final_time && ts->ptime > ts->max_time) {
1945: TSInterpolate(ts,ts->max_time,x);
1946: if (ftime) *ftime = ts->max_time;
1947: } else {
1948: VecCopy(ts->vec_sol,x);
1949: if (ftime) *ftime = ts->ptime;
1950: }
1951: }
1952: PetscOptionsGetString(((PetscObject)ts)->prefix,"-ts_view",filename,PETSC_MAX_PATH_LEN,&flg);
1953: if (flg && !PetscPreLoadingOn) {
1954: PetscViewerASCIIOpen(((PetscObject)ts)->comm,filename,&viewer);
1955: TSView(ts,viewer);
1956: PetscViewerDestroy(&viewer);
1957: }
1958: return(0);
1959: }
1963: /*@
1964: TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet()
1966: Collective on TS
1968: Input Parameters:
1969: + ts - time stepping context obtained from TSCreate()
1970: . step - step number that has just completed
1971: . ptime - model time of the state
1972: - x - state at the current model time
1974: Notes:
1975: TSMonitor() is typically used within the time stepping implementations.
1976: Users might call this function when using the TSStep() interface instead of TSSolve().
1978: Level: advanced
1980: .keywords: TS, timestep
1981: @*/
1982: PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec x)
1983: {
1985: PetscInt i,n = ts->numbermonitors;
1988: for (i=0; i<n; i++) {
1989: (*ts->monitor[i])(ts,step,ptime,x,ts->monitorcontext[i]);
1990: }
1991: return(0);
1992: }
1994: /* ------------------------------------------------------------------------*/
1998: /*@C
1999: TSMonitorLGCreate - Creates a line graph context for use with
2000: TS to monitor convergence of preconditioned residual norms.
2002: Collective on TS
2004: Input Parameters:
2005: + host - the X display to open, or null for the local machine
2006: . label - the title to put in the title bar
2007: . x, y - the screen coordinates of the upper left coordinate of the window
2008: - m, n - the screen width and height in pixels
2010: Output Parameter:
2011: . draw - the drawing context
2013: Options Database Key:
2014: . -ts_monitor_draw - automatically sets line graph monitor
2016: Notes:
2017: Use TSMonitorLGDestroy() to destroy this line graph, not PetscDrawLGDestroy().
2019: Level: intermediate
2021: .keywords: TS, monitor, line graph, residual, seealso
2023: .seealso: TSMonitorLGDestroy(), TSMonitorSet()
2025: @*/
2026: PetscErrorCode TSMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw)
2027: {
2028: PetscDraw win;
2032: PetscDrawCreate(PETSC_COMM_SELF,host,label,x,y,m,n,&win);
2033: PetscDrawSetType(win,PETSC_DRAW_X);
2034: PetscDrawLGCreate(win,1,draw);
2035: PetscDrawLGIndicateDataPoints(*draw);
2037: PetscLogObjectParent(*draw,win);
2038: return(0);
2039: }
2043: PetscErrorCode TSMonitorLG(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
2044: {
2045: PetscDrawLG lg = (PetscDrawLG) monctx;
2046: PetscReal x,y = ptime;
2050: if (!monctx) {
2051: MPI_Comm comm;
2052: PetscViewer viewer;
2054: PetscObjectGetComm((PetscObject)ts,&comm);
2055: viewer = PETSC_VIEWER_DRAW_(comm);
2056: PetscViewerDrawGetDrawLG(viewer,0,&lg);
2057: }
2059: if (!n) {PetscDrawLGReset(lg);}
2060: x = (PetscReal)n;
2061: PetscDrawLGAddPoint(lg,&x,&y);
2062: if (n < 20 || (n % 5)) {
2063: PetscDrawLGDraw(lg);
2064: }
2065: return(0);
2066: }
2070: /*@C
2071: TSMonitorLGDestroy - Destroys a line graph context that was created
2072: with TSMonitorLGCreate().
2074: Collective on PetscDrawLG
2076: Input Parameter:
2077: . draw - the drawing context
2079: Level: intermediate
2081: .keywords: TS, monitor, line graph, destroy
2083: .seealso: TSMonitorLGCreate(), TSMonitorSet(), TSMonitorLG();
2084: @*/
2085: PetscErrorCode TSMonitorLGDestroy(PetscDrawLG *drawlg)
2086: {
2087: PetscDraw draw;
2091: PetscDrawLGGetDraw(*drawlg,&draw);
2092: PetscDrawDestroy(&draw);
2093: PetscDrawLGDestroy(drawlg);
2094: return(0);
2095: }
2099: /*@
2100: TSGetTime - Gets the current time.
2102: Not Collective
2104: Input Parameter:
2105: . ts - the TS context obtained from TSCreate()
2107: Output Parameter:
2108: . t - the current time
2110: Level: beginner
2112: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()
2114: .keywords: TS, get, time
2115: @*/
2116: PetscErrorCode TSGetTime(TS ts,PetscReal* t)
2117: {
2121: *t = ts->ptime;
2122: return(0);
2123: }
2127: /*@
2128: TSSetTime - Allows one to reset the time.
2130: Logically Collective on TS
2132: Input Parameters:
2133: + ts - the TS context obtained from TSCreate()
2134: - time - the time
2136: Level: intermediate
2138: .seealso: TSGetTime(), TSSetDuration()
2140: .keywords: TS, set, time
2141: @*/
2142: PetscErrorCode TSSetTime(TS ts, PetscReal t)
2143: {
2147: ts->ptime = t;
2148: return(0);
2149: }
2153: /*@C
2154: TSSetOptionsPrefix - Sets the prefix used for searching for all
2155: TS options in the database.
2157: Logically Collective on TS
2159: Input Parameter:
2160: + ts - The TS context
2161: - prefix - The prefix to prepend to all option names
2163: Notes:
2164: A hyphen (-) must NOT be given at the beginning of the prefix name.
2165: The first character of all runtime options is AUTOMATICALLY the
2166: hyphen.
2168: Level: advanced
2170: .keywords: TS, set, options, prefix, database
2172: .seealso: TSSetFromOptions()
2174: @*/
2175: PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[])
2176: {
2178: SNES snes;
2182: PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);
2183: TSGetSNES(ts,&snes);
2184: SNESSetOptionsPrefix(snes,prefix);
2185: return(0);
2186: }
2191: /*@C
2192: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
2193: TS options in the database.
2195: Logically Collective on TS
2197: Input Parameter:
2198: + ts - The TS context
2199: - prefix - The prefix to prepend to all option names
2201: Notes:
2202: A hyphen (-) must NOT be given at the beginning of the prefix name.
2203: The first character of all runtime options is AUTOMATICALLY the
2204: hyphen.
2206: Level: advanced
2208: .keywords: TS, append, options, prefix, database
2210: .seealso: TSGetOptionsPrefix()
2212: @*/
2213: PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[])
2214: {
2216: SNES snes;
2220: PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);
2221: TSGetSNES(ts,&snes);
2222: SNESAppendOptionsPrefix(snes,prefix);
2223: return(0);
2224: }
2228: /*@C
2229: TSGetOptionsPrefix - Sets the prefix used for searching for all
2230: TS options in the database.
2232: Not Collective
2234: Input Parameter:
2235: . ts - The TS context
2237: Output Parameter:
2238: . prefix - A pointer to the prefix string used
2240: Notes: On the fortran side, the user should pass in a string 'prifix' of
2241: sufficient length to hold the prefix.
2243: Level: intermediate
2245: .keywords: TS, get, options, prefix, database
2247: .seealso: TSAppendOptionsPrefix()
2248: @*/
2249: PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[])
2250: {
2256: PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);
2257: return(0);
2258: }
2262: /*@C
2263: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
2265: Not Collective, but parallel objects are returned if TS is parallel
2267: Input Parameter:
2268: . ts - The TS context obtained from TSCreate()
2270: Output Parameters:
2271: + J - The Jacobian J of F, where U_t = F(U,t)
2272: . M - The preconditioner matrix, usually the same as J
2273: . func - Function to compute the Jacobian of the RHS
2274: - ctx - User-defined context for Jacobian evaluation routine
2276: Notes: You can pass in PETSC_NULL for any return argument you do not need.
2278: Level: intermediate
2280: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber()
2282: .keywords: TS, timestep, get, matrix, Jacobian
2283: @*/
2284: PetscErrorCode TSGetRHSJacobian(TS ts,Mat *J,Mat *M,TSRHSJacobian *func,void **ctx)
2285: {
2287: SNES snes;
2290: TSGetSNES(ts,&snes);
2291: SNESGetJacobian(snes,J,M,PETSC_NULL,PETSC_NULL);
2292: if (func) *func = ts->userops->rhsjacobian;
2293: if (ctx) *ctx = ts->jacP;
2294: return(0);
2295: }
2299: /*@C
2300: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
2302: Not Collective, but parallel objects are returned if TS is parallel
2304: Input Parameter:
2305: . ts - The TS context obtained from TSCreate()
2307: Output Parameters:
2308: + A - The Jacobian of F(t,U,U_t)
2309: . B - The preconditioner matrix, often the same as A
2310: . f - The function to compute the matrices
2311: - ctx - User-defined context for Jacobian evaluation routine
2313: Notes: You can pass in PETSC_NULL for any return argument you do not need.
2315: Level: advanced
2317: .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber()
2319: .keywords: TS, timestep, get, matrix, Jacobian
2320: @*/
2321: PetscErrorCode TSGetIJacobian(TS ts,Mat *A,Mat *B,TSIJacobian *f,void **ctx)
2322: {
2324: SNES snes;
2327: TSGetSNES(ts,&snes);
2328: SNESGetJacobian(snes,A,B,PETSC_NULL,PETSC_NULL);
2329: if (f) *f = ts->userops->ijacobian;
2330: if (ctx) *ctx = ts->jacP;
2331: return(0);
2332: }
2334: typedef struct {
2335: PetscViewer viewer;
2336: Vec initialsolution;
2337: PetscBool showinitial;
2338: } TSMonitorSolutionCtx;
2342: /*@C
2343: TSMonitorSolution - Monitors progress of the TS solvers by calling
2344: VecView() for the solution at each timestep
2346: Collective on TS
2348: Input Parameters:
2349: + ts - the TS context
2350: . step - current time-step
2351: . ptime - current time
2352: - dummy - either a viewer or PETSC_NULL
2354: Level: intermediate
2356: .keywords: TS, vector, monitor, view
2358: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
2359: @*/
2360: PetscErrorCode TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec x,void *dummy)
2361: {
2362: PetscErrorCode ierr;
2363: TSMonitorSolutionCtx *ictx = (TSMonitorSolutionCtx*)dummy;
2366: if (!step && ictx->showinitial) {
2367: if (!ictx->initialsolution) {
2368: VecDuplicate(x,&ictx->initialsolution);
2369: }
2370: VecCopy(x,ictx->initialsolution);
2371: }
2372: if (ictx->showinitial) {
2373: PetscReal pause;
2374: PetscViewerDrawGetPause(ictx->viewer,&pause);
2375: PetscViewerDrawSetPause(ictx->viewer,0.0);
2376: VecView(ictx->initialsolution,ictx->viewer);
2377: PetscViewerDrawSetPause(ictx->viewer,pause);
2378: PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);
2379: }
2380: VecView(x,ictx->viewer);
2381: if (ictx->showinitial) {
2382: PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);
2383: }
2384: return(0);
2385: }
2390: /*@C
2391: TSMonitorSolutionDestroy - Destroys the monitor context for TSMonitorSolution
2393: Collective on TS
2395: Input Parameters:
2396: . ctx - the monitor context
2398: Level: intermediate
2400: .keywords: TS, vector, monitor, view
2402: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorSolution()
2403: @*/
2404: PetscErrorCode TSMonitorSolutionDestroy(void **ctx)
2405: {
2406: PetscErrorCode ierr;
2407: TSMonitorSolutionCtx *ictx = *(TSMonitorSolutionCtx**)ctx;
2408:
2410: PetscViewerDestroy(&ictx->viewer);
2411: VecDestroy(&ictx->initialsolution);
2412: PetscFree(ictx);
2413: return(0);
2414: }
2418: /*@C
2419: TSMonitorSolutionCreate - Creates the monitor context for TSMonitorSolution
2421: Collective on TS
2423: Input Parameter:
2424: . ts - time-step context
2426: Output Patameter:
2427: . ctx - the monitor context
2429: Level: intermediate
2431: .keywords: TS, vector, monitor, view
2433: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorSolution()
2434: @*/
2435: PetscErrorCode TSMonitorSolutionCreate(TS ts,PetscViewer viewer,void **ctx)
2436: {
2437: PetscErrorCode ierr;
2438: TSMonitorSolutionCtx *ictx;
2439:
2441: PetscNew(TSMonitorSolutionCtx,&ictx);
2442: *ctx = (void*)ictx;
2443: if (!viewer) {
2444: viewer = PETSC_VIEWER_DRAW_(((PetscObject)ts)->comm);
2445: }
2446: PetscObjectReference((PetscObject)viewer);
2447: ictx->viewer = viewer;
2448: ictx->showinitial = PETSC_FALSE;
2449: PetscOptionsGetBool(((PetscObject)ts)->prefix,"-ts_monitor_solution_initial",&ictx->showinitial,PETSC_NULL);
2450: return(0);
2451: }
2455: /*@
2456: TSSetDM - Sets the DM that may be used by some preconditioners
2458: Logically Collective on TS and DM
2460: Input Parameters:
2461: + ts - the preconditioner context
2462: - dm - the dm
2464: Level: intermediate
2467: .seealso: TSGetDM(), SNESSetDM(), SNESGetDM()
2468: @*/
2469: PetscErrorCode TSSetDM(TS ts,DM dm)
2470: {
2472: SNES snes;
2476: PetscObjectReference((PetscObject)dm);
2477: DMDestroy(&ts->dm);
2478: ts->dm = dm;
2479: TSGetSNES(ts,&snes);
2480: SNESSetDM(snes,dm);
2481: return(0);
2482: }
2486: /*@
2487: TSGetDM - Gets the DM that may be used by some preconditioners
2489: Not Collective
2491: Input Parameter:
2492: . ts - the preconditioner context
2494: Output Parameter:
2495: . dm - the dm
2497: Level: intermediate
2500: .seealso: TSSetDM(), SNESSetDM(), SNESGetDM()
2501: @*/
2502: PetscErrorCode TSGetDM(TS ts,DM *dm)
2503: {
2506: *dm = ts->dm;
2507: return(0);
2508: }
2512: /*@
2513: SNESTSFormFunction - Function to evaluate nonlinear residual
2515: Logically Collective on SNES
2517: Input Parameter:
2518: + snes - nonlinear solver
2519: . X - the current state at which to evaluate the residual
2520: - ctx - user context, must be a TS
2522: Output Parameter:
2523: . F - the nonlinear residual
2525: Notes:
2526: This function is not normally called by users and is automatically registered with the SNES used by TS.
2527: It is most frequently passed to MatFDColoringSetFunction().
2529: Level: advanced
2531: .seealso: SNESSetFunction(), MatFDColoringSetFunction()
2532: @*/
2533: PetscErrorCode SNESTSFormFunction(SNES snes,Vec X,Vec F,void *ctx)
2534: {
2535: TS ts = (TS)ctx;
2543: (ts->ops->snesfunction)(snes,X,F,ts);
2544: return(0);
2545: }
2549: /*@
2550: SNESTSFormJacobian - Function to evaluate the Jacobian
2552: Collective on SNES
2554: Input Parameter:
2555: + snes - nonlinear solver
2556: . X - the current state at which to evaluate the residual
2557: - ctx - user context, must be a TS
2559: Output Parameter:
2560: + A - the Jacobian
2561: . B - the preconditioning matrix (may be the same as A)
2562: - flag - indicates any structure change in the matrix
2564: Notes:
2565: This function is not normally called by users and is automatically registered with the SNES used by TS.
2567: Level: developer
2569: .seealso: SNESSetJacobian()
2570: @*/
2571: PetscErrorCode SNESTSFormJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flag,void *ctx)
2572: {
2573: TS ts = (TS)ctx;
2585: (ts->ops->snesjacobian)(snes,X,A,B,flag,ts);
2586: return(0);
2587: }
2591: /*@C
2592: TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only
2594: Collective on TS
2596: Input Arguments:
2597: + ts - time stepping context
2598: . t - time at which to evaluate
2599: . X - state at which to evaluate
2600: - ctx - context
2602: Output Arguments:
2603: . F - right hand side
2605: Level: intermediate
2607: Notes:
2608: This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems.
2609: The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian().
2611: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
2612: @*/
2613: PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec X,Vec F,void *ctx)
2614: {
2616: Mat Arhs,Brhs;
2617: MatStructure flg2;
2620: TSGetRHSMats_Private(ts,&Arhs,&Brhs);
2621: TSComputeRHSJacobian(ts,t,X,&Arhs,&Brhs,&flg2);
2622: MatMult(Arhs,X,F);
2623: return(0);
2624: }
2628: /*@C
2629: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
2631: Collective on TS
2633: Input Arguments:
2634: + ts - time stepping context
2635: . t - time at which to evaluate
2636: . X - state at which to evaluate
2637: - ctx - context
2639: Output Arguments:
2640: + A - pointer to operator
2641: . B - pointer to preconditioning matrix
2642: - flg - matrix structure flag
2644: Level: intermediate
2646: Notes:
2647: This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems.
2649: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear()
2650: @*/
2651: PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec X,Mat *A,Mat *B,MatStructure *flg,void *ctx)
2652: {
2655: *flg = SAME_PRECONDITIONER;
2656: return(0);
2657: }
2661: /*@C
2662: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
2664: Collective on TS
2666: Input Arguments:
2667: + ts - time stepping context
2668: . t - time at which to evaluate
2669: . X - state at which to evaluate
2670: . Xdot - time derivative of state vector
2671: - ctx - context
2673: Output Arguments:
2674: . F - left hand side
2676: Level: intermediate
2678: Notes:
2679: The assumption here is that the left hand side is of the form A*Xdot (and not A*Xdot + B*X). For other cases, the
2680: user is required to write their own TSComputeIFunction.
2681: This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems.
2682: The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian().
2684: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant()
2685: @*/
2686: PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
2687: {
2689: Mat A,B;
2690: MatStructure flg2;
2693: TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);
2694: TSComputeIJacobian(ts,t,X,Xdot,1.0,&A,&B,&flg2,PETSC_TRUE);
2695: MatMult(A,Xdot,F);
2696: return(0);
2697: }
2701: /*@C
2702: TSComputeIJacobianConstant - Reuses a Jacobian that is time-independent.
2704: Collective on TS
2706: Input Arguments:
2707: + ts - time stepping context
2708: . t - time at which to evaluate
2709: . X - state at which to evaluate
2710: . Xdot - time derivative of state vector
2711: . shift - shift to apply
2712: - ctx - context
2714: Output Arguments:
2715: + A - pointer to operator
2716: . B - pointer to preconditioning matrix
2717: - flg - matrix structure flag
2719: Level: intermediate
2721: Notes:
2722: This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems.
2724: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear()
2725: @*/
2726: PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx)
2727: {
2730: *flg = SAME_PRECONDITIONER;
2731: return(0);
2732: }
2737: /*@
2738: TSGetConvergedReason - Gets the reason the TS iteration was stopped.
2740: Not Collective
2742: Input Parameter:
2743: . ts - the TS context
2745: Output Parameter:
2746: . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
2747: manual pages for the individual convergence tests for complete lists
2749: Level: intermediate
2751: Notes:
2752: Can only be called after the call to TSSolve() is complete.
2754: .keywords: TS, nonlinear, set, convergence, test
2756: .seealso: TSSetConvergenceTest(), TSConvergedReason
2757: @*/
2758: PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason)
2759: {
2763: *reason = ts->reason;
2764: return(0);
2765: }
2769: /*@
2770: TSGetNonlinearSolveIterations - Gets the total number of linear iterations
2771: used by the time integrator.
2773: Not Collective
2775: Input Parameter:
2776: . ts - TS context
2778: Output Parameter:
2779: . nits - number of nonlinear iterations
2781: Notes:
2782: This counter is reset to zero for each successive call to TSSolve().
2784: Level: intermediate
2786: .keywords: TS, get, number, nonlinear, iterations
2788: .seealso: TSGetLinearSolveIterations()
2789: @*/
2790: PetscErrorCode TSGetNonlinearSolveIterations(TS ts,PetscInt *nits)
2791: {
2795: *nits = ts->nonlinear_its;
2796: return(0);
2797: }
2801: /*@
2802: TSGetLinearSolveIterations - Gets the total number of linear iterations
2803: used by the time integrator.
2805: Not Collective
2807: Input Parameter:
2808: . ts - TS context
2810: Output Parameter:
2811: . lits - number of linear iterations
2813: Notes:
2814: This counter is reset to zero for each successive call to TSSolve().
2816: Level: intermediate
2818: .keywords: TS, get, number, linear, iterations
2820: .seealso: TSGetNonlinearSolveIterations(), SNESGetLinearSolveIterations()
2821: @*/
2822: PetscErrorCode TSGetLinearSolveIterations(TS ts,PetscInt *lits)
2823: {
2827: *lits = ts->linear_its;
2828: return(0);
2829: }
2833: /*@C
2834: TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file
2836: Collective on TS
2838: Input Parameters:
2839: + ts - the TS context
2840: . step - current time-step
2841: . ptime - current time
2842: . x - current state
2843: - viewer - binary viewer
2845: Level: intermediate
2847: .keywords: TS, vector, monitor, view
2849: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
2850: @*/
2851: PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec x,void *viewer)
2852: {
2853: PetscErrorCode ierr;
2854: PetscViewer v = (PetscViewer)viewer;
2857: VecView(x,v);
2858: return(0);
2859: }
2863: /*@C
2864: TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep.
2866: Collective on TS
2868: Input Parameters:
2869: + ts - the TS context
2870: . step - current time-step
2871: . ptime - current time
2872: . x - current state
2873: - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
2875: Level: intermediate
2877: Notes:
2878: 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.
2879: These are named according to the file name template.
2881: This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy().
2883: .keywords: TS, vector, monitor, view
2885: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
2886: @*/
2887: PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec x,void *filenametemplate)
2888: {
2890: char filename[PETSC_MAX_PATH_LEN];
2891: PetscViewer viewer;
2894: PetscSNPrintf(filename,sizeof filename,(const char*)filenametemplate,step);
2895: PetscViewerVTKOpen(((PetscObject)ts)->comm,filename,FILE_MODE_WRITE,&viewer);
2896: VecView(x,viewer);
2897: PetscViewerDestroy(&viewer);
2898: return(0);
2899: }
2903: /*@C
2904: TSMonitorSolutionVTKDestroy - Destroy context for monitoring
2906: Collective on TS
2908: Input Parameters:
2909: . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)
2911: Level: intermediate
2913: Note:
2914: This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK().
2916: .keywords: TS, vector, monitor, view
2918: .seealso: TSMonitorSet(), TSMonitorSolutionVTK()
2919: @*/
2920: PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate)
2921: {
2925: PetscFree(*(char**)filenametemplate);
2926: return(0);
2927: }
2931: /*@
2932: TSGetAdapt - Get the adaptive controller context for the current method
2934: Collective on TS if controller has not been created yet
2936: Input Arguments:
2937: . ts - time stepping context
2939: Output Arguments:
2940: . adapt - adaptive controller
2942: Level: intermediate
2944: .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose()
2945: @*/
2946: PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt)
2947: {
2953: if (!ts->adapt) {
2954: TSAdaptCreate(((PetscObject)ts)->comm,&ts->adapt);
2955: PetscLogObjectParent(ts,ts->adapt);
2956: PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);
2957: }
2958: *adapt = ts->adapt;
2959: return(0);
2960: }
2964: /*@
2965: TSSetTolerances - Set tolerances for local truncation error when using adaptive controller
2967: Logically Collective
2969: Input Arguments:
2970: + ts - time integration context
2971: . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value
2972: . vatol - vector of absolute tolerances or PETSC_NULL, used in preference to atol if present
2973: . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value
2974: - vrtol - vector of relative tolerances or PETSC_NULL, used in preference to atol if present
2976: Level: beginner
2978: .seealso: TS, TSAdapt, TSVecNormWRMS()
2979: @*/
2980: PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol)
2981: {
2985: if (atol != PETSC_DECIDE) ts->atol = atol;
2986: if (vatol) {
2987: PetscObjectReference((PetscObject)vatol);
2988: VecDestroy(&ts->vatol);
2989: ts->vatol = vatol;
2990: }
2991: if (rtol != PETSC_DECIDE) ts->rtol = rtol;
2992: if (vrtol) {
2993: PetscObjectReference((PetscObject)vrtol);
2994: VecDestroy(&ts->vrtol);
2995: ts->vrtol = vrtol;
2996: }
2997: return(0);
2998: }
3002: /*@
3003: TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state
3005: Collective on TS
3007: Input Arguments:
3008: + ts - time stepping context
3009: - Y - state vector to be compared to ts->vec_sol
3011: Output Arguments:
3012: . norm - weighted norm, a value of 1.0 is considered small
3014: Level: developer
3016: .seealso: TSSetTolerances()
3017: @*/
3018: PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm)
3019: {
3021: PetscInt i,n,N;
3022: const PetscScalar *x,*y;
3023: Vec X;
3024: PetscReal sum,gsum;
3030: X = ts->vec_sol;
3032: if (X == Y) SETERRQ(((PetscObject)X)->comm,PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector");
3034: /* This is simple to implement, just not done yet */
3035: if (ts->vatol || ts->vrtol) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_SUP,"No support for vector scaling yet");
3037: VecGetSize(X,&N);
3038: VecGetLocalSize(X,&n);
3039: VecGetArrayRead(X,&x);
3040: VecGetArrayRead(Y,&y);
3041: sum = 0.;
3042: for (i=0; i<n; i++) {
3043: PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(x[i]),PetscAbsScalar(y[i]));
3044: sum += PetscSqr(PetscAbsScalar(y[i] - x[i]) / tol);
3045: }
3046: VecRestoreArrayRead(X,&x);
3047: VecRestoreArrayRead(Y,&y);
3049: MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,((PetscObject)ts)->comm);
3050: *norm = PetscSqrtReal(gsum / N);
3051: if (PetscIsInfOrNanScalar(*norm)) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
3052: return(0);
3053: }
3057: /*@
3058: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
3060: Logically Collective on TS
3062: Input Arguments:
3063: + ts - time stepping context
3064: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
3066: Note:
3067: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
3069: Level: intermediate
3071: .seealso: TSGetCFLTime(), TSADAPTCFL
3072: @*/
3073: PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime)
3074: {
3078: ts->cfltime_local = cfltime;
3079: ts->cfltime = -1.;
3080: return(0);
3081: }
3085: /*@
3086: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
3088: Collective on TS
3090: Input Arguments:
3091: . ts - time stepping context
3093: Output Arguments:
3094: . cfltime - maximum stable time step for forward Euler
3096: Level: advanced
3098: .seealso: TSSetCFLTimeLocal()
3099: @*/
3100: PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime)
3101: {
3105: if (ts->cfltime < 0) {
3106: MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,((PetscObject)ts)->comm);
3107: }
3108: *cfltime = ts->cfltime;
3109: return(0);
3110: }
3114: /*@
3115: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
3117: Input Parameters:
3118: . ts - the TS context.
3119: . xl - lower bound.
3120: . xu - upper bound.
3122: Notes:
3123: If this routine is not called then the lower and upper bounds are set to
3124: SNES_VI_NINF and SNES_VI_INF respectively during SNESSetUp().
3126: Level: advanced
3128: @*/
3129: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
3130: {
3132: SNES snes;
3135: TSGetSNES(ts,&snes);
3136: SNESVISetVariableBounds(snes,xl,xu);
3137: return(0);
3138: }
3140: #if defined(PETSC_HAVE_MATLAB_ENGINE)
3141: #include <mex.h>
3143: typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext;
3147: /*
3148: TSComputeFunction_Matlab - Calls the function that has been set with
3149: TSSetFunctionMatlab().
3151: Collective on TS
3153: Input Parameters:
3154: + snes - the TS context
3155: - x - input vector
3157: Output Parameter:
3158: . y - function vector, as set by TSSetFunction()
3160: Notes:
3161: TSComputeFunction() is typically used within nonlinear solvers
3162: implementations, so most users would not generally call this routine
3163: themselves.
3165: Level: developer
3167: .keywords: TS, nonlinear, compute, function
3169: .seealso: TSSetFunction(), TSGetFunction()
3170: */
3171: PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec x,Vec xdot,Vec y, void *ctx)
3172: {
3173: PetscErrorCode ierr;
3174: TSMatlabContext *sctx = (TSMatlabContext *)ctx;
3175: int nlhs = 1,nrhs = 7;
3176: mxArray *plhs[1],*prhs[7];
3177: long long int lx = 0,lxdot = 0,ly = 0,ls = 0;
3187: PetscMemcpy(&ls,&snes,sizeof(snes));
3188: PetscMemcpy(&lx,&x,sizeof(x));
3189: PetscMemcpy(&lxdot,&xdot,sizeof(xdot));
3190: PetscMemcpy(&ly,&y,sizeof(x));
3191: prhs[0] = mxCreateDoubleScalar((double)ls);
3192: prhs[1] = mxCreateDoubleScalar(time);
3193: prhs[2] = mxCreateDoubleScalar((double)lx);
3194: prhs[3] = mxCreateDoubleScalar((double)lxdot);
3195: prhs[4] = mxCreateDoubleScalar((double)ly);
3196: prhs[5] = mxCreateString(sctx->funcname);
3197: prhs[6] = sctx->ctx;
3198: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");
3199: mxGetScalar(plhs[0]);
3200: mxDestroyArray(prhs[0]);
3201: mxDestroyArray(prhs[1]);
3202: mxDestroyArray(prhs[2]);
3203: mxDestroyArray(prhs[3]);
3204: mxDestroyArray(prhs[4]);
3205: mxDestroyArray(prhs[5]);
3206: mxDestroyArray(plhs[0]);
3207: return(0);
3208: }
3213: /*
3214: TSSetFunctionMatlab - Sets the function evaluation routine and function
3215: vector for use by the TS routines in solving ODEs
3216: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
3218: Logically Collective on TS
3220: Input Parameters:
3221: + ts - the TS context
3222: - func - function evaluation routine
3224: Calling sequence of func:
3225: $ func (TS ts,PetscReal time,Vec x,Vec xdot,Vec f,void *ctx);
3227: Level: beginner
3229: .keywords: TS, nonlinear, set, function
3231: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
3232: */
3233: PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx)
3234: {
3235: PetscErrorCode ierr;
3236: TSMatlabContext *sctx;
3239: /* currently sctx is memory bleed */
3240: PetscMalloc(sizeof(TSMatlabContext),&sctx);
3241: PetscStrallocpy(func,&sctx->funcname);
3242: /*
3243: This should work, but it doesn't
3244: sctx->ctx = ctx;
3245: mexMakeArrayPersistent(sctx->ctx);
3246: */
3247: sctx->ctx = mxDuplicateArray(ctx);
3248: TSSetIFunction(ts,PETSC_NULL,TSComputeFunction_Matlab,sctx);
3249: return(0);
3250: }
3254: /*
3255: TSComputeJacobian_Matlab - Calls the function that has been set with
3256: TSSetJacobianMatlab().
3258: Collective on TS
3260: Input Parameters:
3261: + ts - the TS context
3262: . x - input vector
3263: . A, B - the matrices
3264: - ctx - user context
3266: Output Parameter:
3267: . flag - structure of the matrix
3269: Level: developer
3271: .keywords: TS, nonlinear, compute, function
3273: .seealso: TSSetFunction(), TSGetFunction()
3274: @*/
3275: PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec x,Vec xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx)
3276: {
3277: PetscErrorCode ierr;
3278: TSMatlabContext *sctx = (TSMatlabContext *)ctx;
3279: int nlhs = 2,nrhs = 9;
3280: mxArray *plhs[2],*prhs[9];
3281: long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0;
3287: /* call Matlab function in ctx with arguments x and y */
3289: PetscMemcpy(&ls,&ts,sizeof(ts));
3290: PetscMemcpy(&lx,&x,sizeof(x));
3291: PetscMemcpy(&lxdot,&xdot,sizeof(x));
3292: PetscMemcpy(&lA,A,sizeof(x));
3293: PetscMemcpy(&lB,B,sizeof(x));
3294: prhs[0] = mxCreateDoubleScalar((double)ls);
3295: prhs[1] = mxCreateDoubleScalar((double)time);
3296: prhs[2] = mxCreateDoubleScalar((double)lx);
3297: prhs[3] = mxCreateDoubleScalar((double)lxdot);
3298: prhs[4] = mxCreateDoubleScalar((double)shift);
3299: prhs[5] = mxCreateDoubleScalar((double)lA);
3300: prhs[6] = mxCreateDoubleScalar((double)lB);
3301: prhs[7] = mxCreateString(sctx->funcname);
3302: prhs[8] = sctx->ctx;
3303: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");
3304: mxGetScalar(plhs[0]);
3305: *flag = (MatStructure) mxGetScalar(plhs[1]);
3306: mxDestroyArray(prhs[0]);
3307: mxDestroyArray(prhs[1]);
3308: mxDestroyArray(prhs[2]);
3309: mxDestroyArray(prhs[3]);
3310: mxDestroyArray(prhs[4]);
3311: mxDestroyArray(prhs[5]);
3312: mxDestroyArray(prhs[6]);
3313: mxDestroyArray(prhs[7]);
3314: mxDestroyArray(plhs[0]);
3315: mxDestroyArray(plhs[1]);
3316: return(0);
3317: }
3322: /*
3323: TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
3324: vector for use by the TS routines in solving ODEs from MATLAB. Here the function is a string containing the name of a MATLAB function
3326: Logically Collective on TS
3328: Input Parameters:
3329: + ts - the TS context
3330: . A,B - Jacobian matrices
3331: . func - function evaluation routine
3332: - ctx - user context
3334: Calling sequence of func:
3335: $ flag = func (TS ts,PetscReal time,Vec x,Vec xdot,Mat A,Mat B,void *ctx);
3338: Level: developer
3340: .keywords: TS, nonlinear, set, function
3342: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
3343: */
3344: PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx)
3345: {
3346: PetscErrorCode ierr;
3347: TSMatlabContext *sctx;
3350: /* currently sctx is memory bleed */
3351: PetscMalloc(sizeof(TSMatlabContext),&sctx);
3352: PetscStrallocpy(func,&sctx->funcname);
3353: /*
3354: This should work, but it doesn't
3355: sctx->ctx = ctx;
3356: mexMakeArrayPersistent(sctx->ctx);
3357: */
3358: sctx->ctx = mxDuplicateArray(ctx);
3359: TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);
3360: return(0);
3361: }
3365: /*
3366: TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab().
3368: Collective on TS
3370: .seealso: TSSetFunction(), TSGetFunction()
3371: @*/
3372: PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec x, void *ctx)
3373: {
3374: PetscErrorCode ierr;
3375: TSMatlabContext *sctx = (TSMatlabContext *)ctx;
3376: int nlhs = 1,nrhs = 6;
3377: mxArray *plhs[1],*prhs[6];
3378: long long int lx = 0,ls = 0;
3384: PetscMemcpy(&ls,&ts,sizeof(ts));
3385: PetscMemcpy(&lx,&x,sizeof(x));
3386: prhs[0] = mxCreateDoubleScalar((double)ls);
3387: prhs[1] = mxCreateDoubleScalar((double)it);
3388: prhs[2] = mxCreateDoubleScalar((double)time);
3389: prhs[3] = mxCreateDoubleScalar((double)lx);
3390: prhs[4] = mxCreateString(sctx->funcname);
3391: prhs[5] = sctx->ctx;
3392: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");
3393: mxGetScalar(plhs[0]);
3394: mxDestroyArray(prhs[0]);
3395: mxDestroyArray(prhs[1]);
3396: mxDestroyArray(prhs[2]);
3397: mxDestroyArray(prhs[3]);
3398: mxDestroyArray(prhs[4]);
3399: mxDestroyArray(plhs[0]);
3400: return(0);
3401: }
3406: /*
3407: TSMonitorSetMatlab - Sets the monitor function from Matlab
3409: Level: developer
3411: .keywords: TS, nonlinear, set, function
3413: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
3414: */
3415: PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx)
3416: {
3417: PetscErrorCode ierr;
3418: TSMatlabContext *sctx;
3421: /* currently sctx is memory bleed */
3422: PetscMalloc(sizeof(TSMatlabContext),&sctx);
3423: PetscStrallocpy(func,&sctx->funcname);
3424: /*
3425: This should work, but it doesn't
3426: sctx->ctx = ctx;
3427: mexMakeArrayPersistent(sctx->ctx);
3428: */
3429: sctx->ctx = mxDuplicateArray(ctx);
3430: TSMonitorSet(ts,TSMonitor_Matlab,sctx,PETSC_NULL);
3431: return(0);
3432: }
3433: #endif